Balanced Gain-and-Loss Optical Waveguides: Exact Solutions for Guided Modes in Susy-QM

Abstract

The construction of exactly solvable refractive indices allowing guided TE modes in optical waveguides is investigated within the formalism of Darboux–Crum transformations. We apply the finite-difference algorithm for higher-order supersymmetric quantum mechanics to obtain complex-valued refractive indices admitting all-real eigenvalues in their point spectrum. The new refractive indices are such that their imaginary part gives zero if it is integrated over the entire domain of definition. This property, called condition of zero total area, ensures the conservation of optical power so the refractive index shows balanced gain and loss. Consequently, the complex-valued refractive indices reported in this work include but are not limited to the parity-time invariant case.

1. Introduction

Optical waveguides with arbitrarily shaped index profiles have been the subject of intensive research over recent decades. The accelerated development of fiber optics technologies and optical communications has prompted new ways of visualizing the interaction of light with matter. In addition, advances in material and optical technologies open wider possibilities to design index profiles on demand, so more sophisticated waveguiding-based optical devices may be achieved in a short time.

The behavior of light beams in spatially inhomogeneous media of refractive index n is formally studied through the wave equation derived from the Maxwell equations [1,2,3,4]. As the pulse spread at the receiving end of a waveguide can be reduced for refractive indices of quadratic (parabolic) profile, it is expected that shaping the profile more freely we might be able to reduce further this undesired phenomenon [1]. Such expectation motivates the investigation of arbitrarily shaped index profiles addressed to optical waveguiding. Nevertheless, analytical solutions of the equation that rules the propagation of light in this kind of media exist only for limited refractive index profiles. The WKB, Rayleigh-Ritz, power-series expansion, finite element or perturbation methods are then applied to obtain some approximated solutions [1,2,3,4], even numerical integration is usable. In this respect, the transformation introduced by Darboux in 1882 to study Sturm-Liouville problems deserves special attention.

The Darboux method permits the intertwining of two differential equations in such a form that their solutions are connected through a simple differential relationship; further improvements are due to Crum [5,6]. In the 1970s, such a method received a lot of interest in soliton theory since it is linked to the Bäcklund transformation, used to study the soliton-like solutions of diverse nonlinear differential equations [6]. The Darboux approach attracted renewed attention in the 1980s, after the inauguration of supersymmetric quantum mechanics by Witten [7], because it is also in the core of the factorization method introduced by Dirac ‘as a little stratagem to solve the spectral problem for the one-dimensional quantum oscillator’ [5]. Diverse supersymmetric formulations were immediately addressed to the spectral design of exactly solvable models in quantum mechanics [8,9,10,11,12].

Supersymmetric quantum mechanics (Susy-QM) is presently a very robust formulation embracing a very wide set of applications [5,13,14], including optics [15,16,17,18,19,20,21,22,23,24,25], where supersymmetry may be interpreted as describing two light beams of different colors that form a standing periodic interference pattern along the waveguide axis [15]. In the Helmholtz regime, the supersymmetric partners can be constructed to display parabolic [16] or sech-like index profile [17,18,19,20]. They also show identical coefficients of transmission and reflection for any angle of incidence [5,18,19,21,22], which may render them perfectly indistinguishable to an external observer [23]. The impressing breakthrough of supersymmetry in optics is that superpartner configurations are experimentally realized in coupled optical networks, where the light dynamics is directly observed [24] and the conditions for degeneracy breaking can be studied [25]. The theoretical modeling of diverse refractive indices using the powerful Darboux–Crum (supersymmetric) formulation [5,8,9,10,11,12,19,26] is therefore interesting by itself.

The Darboux-deformation of Hermitian systems is just one of the options permitted by the method. Indeed, non-Hermitian structures may be also developed [27,28,29,30,31], where some initial insights connecting supersymmetry with parity-time (PT) symmetry [32] shown additional applications of the supersymmetric approach. In this context, important theoretical achievements within the PT-symmetric formulation [33] envisioned the observation of parity-time symmetry in optical laboratories [34], so the practical implementation of Susy-QM could also be achieved in optics. Moreover, as the reality of the spectrum is not granted a priori for non-Hermitian models, the demonstration that PT symmetry implies real spectrum [35] stimulated the systematic search of PT-symmetric systems in quantum mechanics [36]. Nevertheless, further improvements shown that $\mathrm{PT}$-symmetry is not a necessary condition for the spectrum reality [37,38], a fact confirmed in diverse supersymmetric models of non-Hermiticity [39,40,41,42,43,44], where the imaginary part of a wide class of complex-valued potentials lead to balanced gain and loss probability without the necessity of PT symmetry. The latter property opens new possibilities in optical design, where the gain-loss manipulation includes but is not limited to PT symmetry [45].

In the present work we generate graded refractive indices with both real- and complex-valued profiles. Our interest is addressed to the study of guided TE modes propagating in media with balanced gain-and-loss refractive index. Taking into account the limited set of analytical solutions in optical waveguiding, the method presented here represents a very efficient resource of theoretical models on the matter. The approach is based on the finite-difference algorithm for higher-order supersymmetry introduced in [19], which summarizes the fundamental ingredients of supersymmetric quantum mechanics of any order [5]. It is a generalization of the supersymmetric formulation introduced by Mielnik in 1984 [8], which coincides with the Darboux transformation connecting two exactly solvable potentials in quantum mechanics.

The organization of the paper is as follows. In Section 2 we briefly revisit the mathematical structure underlying the propagation of light in graded-index waveguides. The well-known identification of the Helmholtz equation with the Schrödinger eigenvalue problem is the point of departure to apply the supersymmetric finite-difference algorithm. In Section 2.1 we provide general formulae to add guided TE modes one at a time in the dynamics of a given waveguide. Then, in Section 2.2, a new mechanism to produce complex-valued refractive indices with all-real eigenvalues is discussed. The achievement of this method is that the imaginary part of the new refractive index is such that its integration over the entire domain of definition is equal to zero, preserving in this form the total optical power of the system. Such property permits the acquisition of PT symmetry as a particular case. We also discuss briefly the bi-orthogonal structure that is necessary to properly normalize the guided TE modes of the new system. In Section 3 we present immediate applications. We show that the method produces a wide variety of exactly solvable refractive indices in both the real- and complex-valued configurations. The former case permits the manipulation of the index profile to better-fit it with (possible) experimental observations. The second case leads to PT-symmetric as well as to non-PT-symmetric refractive indices, both configurations admitting all-real propagation constants in the point spectrum. Some conclusions are given in Section 4. For the sake of completeness, we briefly revisit the generalities of the supersymmetric finite-difference algorithm in Appendix A.

2. Mathematical Physics of Graded-Index Waveguides

The Helmholtz equation

$$\frac{{\partial }^2}{\partial x^2}{\Phi }_y(x,z)+\frac{{\partial }^2}{\partial z^2}{\Phi }_y(x,z)+k_0^2{n}^2\left(x\right){\Phi }_y(x,z)=0$$

(1)

is useful to study TE modes propagating through inhomogeneous dielectric materials with negligible magnetic polarizability and dispersion. Here $k_0(=w/c)=2\pi {\lambda }_0^{-1}$ stands for the wavenumber in vacuum, and the refractive index n is assumed to depend only on the x-coordinate. Hence, ${\Phi }_y(x,z)$ represents an electric wave polarized in the y-direction, which propagates along the positive z-direction.

On the other hand, the paraxial Helmholtz equation

$$\frac{1}{2k_0^2{n}_{*}}\frac{{\partial }^2}{\partial x^2}{E}_y(x,z)+\frac{i}{k_0}\frac{\partial }{\partial z}{E}_y(x,z)+[n\left(x\right)-n_{*}]E_y(x,z)=0$$

(2)

results within the paraxial regime for weakly guiding media ($|n-n_{*}|\ll 1$), after using the ansatz ${\Phi }_y(x,z)=E_y(x,z)exp\left(i{k}_0{n}_{*}z\right)$. The number $n_{*}>0$ defines a reference refractive index. As $n\left(x\right)$ is independent of the z-coordinate, let us assume $E_y(x,z)=E\left(x\right)exp(-i{k}_0\epsilon z)$, where $\epsilon$ defines the spectrum propagation constants [46], then (2) is reduced to the eigenvalue problem

$$\left[-\frac{1}{2k_0^2{n}_{*}}\frac{d^2}{d{x}^2}-n\left(x\right)+n_{*}\right]E\left(x\right)=\epsilon E\left(x\right),{\left.E\left(x\right):=E_y(x,z)\right|}_{z=0}.$$

(3)

Comparing Equation (3) with the eigenvalue equation for one-dimensional stationary systems in quantum mechanics

$$\left[-\frac{{\hslash }^2}{2m}\frac{d^2}{d{x}^2}+V\left(x\right)\right]\psi \left(x\right)=\mathcal{E}\psi \left(x\right),$$

(4)

we obtain the identification [47]

$$k_0^2⟷\frac{mw}{\hslash },\left[-n\left(x\right)+n_{*}\right]n_{*}⟷\frac{V\left(x\right)}{\hslash w},\epsilon n_{*}⟷\frac{\mathcal{E}}{\hslash w}.$$

(5)

Then, the TE modes $E_y(x,z)=E\left(x\right)exp(-i{k}_0\epsilon z)$ may be associated with the solutions $\psi (x,t)=\psi \left(x\right)exp(-i\mathcal{E}t/\hslash )$ of the related Schrödinger equation, where $z⟷t\sqrt{w\hslash /m}$.

The link between (3) and (4) is even more clear after introducing the changes $x\to \chi /\left(k_0\sqrt{2}\right)$ and $x\to \chi \sqrt{\hslash /\left(2mw\right)}$, which gives

$$\left[-\frac{d^2}{d{\chi }^2}+(-n\left(\chi \right)+n_{*})n_{*}\right]E\left(\chi \right)=\epsilon n_{*}E\left(\chi \right),\left[-\frac{d^2}{d{\chi }^2}+\mathcal{V}\left(\chi \right)\right]\psi \left(\chi \right)=\tilde{\mathcal{E}}\psi \left(\chi \right),$$

(6)

with $\hslash w\mathcal{V}=V$ and $\hslash w\tilde{\mathcal{E}}=\mathcal{E}$.

In the sequel we take full advantage of the above mathematical relationship, providing solutions to the Helmholtz Equation (3) from the space state of the quantum mechanical problem (4). Our aim is to design refractive indices $n\left(x\right)$ on demand, locating propagation constants at concrete positions of the point spectrum, and producing balanced gain-and-loss. Keeping this in mind we will apply the Darboux method [6], as it has been developed in supersymmetric quantum mechanics [5], to deform the eigenfunctions of a given refractive index into the eigenfunctions of another one. The approach is useful to add/eliminate a concrete number of eigenvalues to/from the point spectrum of $n\left(x\right)$, at the price of transforming $n\left(x\right)$ and $E\left(x\right)$ into new functions. Changing the point spectrum of $n\left(x\right)$ by only one eigenvalue $\epsilon$ corresponds to the first-order Susy (Darboux) transformation. If the modification involves $k\ge 2$ eigenvalues, then one works with the kth-order Susy transformation, which may be performed iterating k-times the Darboux transformation or deforming $n\left(x\right)$ in a single step (Darboux–Crum) [5]. The Darboux method has been elegantly summarized in the finite-difference algorithm for higher-order supersymmetry introduced in [19], which is briefly revisited in Appendix A for completeness.

2.1. Adding Propagation Constants under Prescription

We are interested in generating refractive indices $n\left(x\right)$ with exact solutions to the paraxial Helmholtz Equation (3). An elegant way to achieve this is offered by the finite-difference algorithm for higher-order supersymmetry [19] revisited in Appendix A. The keystone is to have at hand at least one exactly solvable refractive index $n_0\left(x\right)$, the point spectrum of which is inherited to another refractive index $n_1(x;ϵ)$, defined by the Darboux transformation

$$(-n_1(x;ϵ)+n_{1,*})\frac{n_{1,*}}{n_{0,*}}=-n_0\left(x\right)+n_{0,*}+\frac{\sqrt{2}}{k_0{n}_{0,*}}{\beta }_1^{\prime }(x;ϵ),$$

(7)

where $f^{\prime }\left(x\right)=\frac{d}{dx}f\left(x\right)$, the superpotential ${\beta }_1(x;ϵ)$ is solution of the nonlinear Riccati equation

$$-\frac{1}{k_0\sqrt{2}}{\beta }_1^{\prime }(x;ϵ)+{\beta }_1^2(x;ϵ)=n_{0,*}\left[-n_0\left(x\right)+n_{0,*}-ϵ\right],$$

(8)

and the factorization constant $ϵ$ is to be determined (see Appendix A for details).

To simplify notation, without loss of generality, hereafter we make $n_{k,*}=n_{0,*}:=n_{*}$, with $k=1,2,\dots$

The refractive index (7) defines a new paraxial Helmholtz equation

$$\left[-\frac{1}{2k_0^2{n}_{*}}\frac{d^2}{d{x}^2}-n_1(x;ϵ)+n_{*}\right]E_{\left(1\right)}(x;ϵ)=\epsilon E_{\left(1\right)}(x;ϵ),$$

(9)

the solutions of which are constructed from the eigenfunctions $E_{\left(0\right)}\left(x\right)$ of $n_0\left(x\right)$ as follows

$${\mathcal{N}}_{\left(1\right)}^{-1}{E}_{\left(1\right)}(x;ϵ)=E_{\left(0\right)}^{\prime }\left(x\right)+{\beta }_1(x,ϵ)E_{\left(0\right)}\left(x\right),$$

(10)

where ${\mathcal{N}}_{\left(1\right)}$ obeys normalization.

The point spectrum of $n_1\left(x\right)$ will be exactly the same as the one of $n_0\left(x\right)$, or it may include the factorization constant $ϵ$ as an additional eigenvalue. In the latter case, the corresponding eigenfunction is written as follows

$$E_{\left(1\right)}^M(x;ϵ)={\mathcal{N}}_{\left(1\right)}^{M}exp\left[\int {\beta }_1(x;ϵ)dx\right].$$

(11)

The procedure described above can be iterated at will, see Appendix A. After k steps one obtains a refractive index $n_k(x;ϵ)$ that admits k additional eigenvalues in its point spectrum with respect to the point spectrum of $n_0\left(x\right)$.

Clearly, the profile of $n_k(x;ϵ)$ may be manipulated to obtain a concrete number of guided TE modes under prescription. For instance, suppose that $n_0\left(x\right)$ in (7) admits only N guided modes, the finite-difference algorithm provides $n_1(x;ϵ)$ with exactly the same propagation constants as $n_0\left(x\right)$ plus an additional one, located at $ϵ$ in the point spectrum of $n_1(x;ϵ)$ but absent in the spectrum of $n_0\left(x\right)$.

Depending on the prescription, the propagation constants of the guided modes may be added one at a time, iterating the finite-difference algorithm as necessary, or using a single transformation if just one additional propagation constant is required. Conventional supersymmetric approaches include the new eigenvalue below the lowest eigenvalue of the previous spectrum. The latter obeys the fact that the oscillation theorem prohibits constructing regular superpotentials ${\beta }_k$ if the factorization constant $ϵ$ is above the lowest propagation constant of $n_{k-1}$. Thus, viable factorization constants $ϵ$ are at most equal to the lowest propagation constant of $n_{k-1}$. In contraposition, a mechanism producing complex-valued refractive indices as Darboux-deformed versions of real-valued ones can be managed by following the method introduced in [39]. The novelty is that the eigenvalues of the new refractive indices will be all-real. Moreover, the additional eigenvalue $ϵ$ can be incorporated at any position of the spectrum [41]. The formulae (7)–(11) still apply for complex-valued refractive indices, and are particularly important to generate balanced gain-and-loss.

2.2. Balanced Gain-and-Loss Waveguides

Conventional supersymmetric approaches assume that the solution of the Riccati Equation (8) is real-valued. However, complex-valued solutions are feasible even for real-valued refractive indices $n_0\left(x\right)$ and real factorization energies $ϵ$. Indeed, the real and imaginary parts of Equation (8) give rise to the nonlinear differential equation

$$\left[-\frac{1}{2k_0^2{n}_{*}}\frac{d^2}{d{x}^2}-n_0\left(x\right)+n_{*}\right]\alpha (x;ϵ)=ϵ\alpha (x;ϵ)-\frac{1}{2k_0^2{n}_{*}}\frac{{\lambda }^2}{{\alpha }^3(x;ϵ)},\lambda \in \mathbb{R},$$

(12)

which is named after Ermakov [48]. Please note that (12) coincides with the paraxial Helmholtz Equation (3) when $\lambda =0$. Then, the solutions may be expressed as the nonlinear superposition [39]

$$\alpha (x;ϵ)={\left[a{u}_{\left(1\right);1}^2(x;ϵ)+b{u}_{\left(1\right);1}(x;ϵ)u_{\left(1\right);2}(x;ϵ)+c{u}_{\left(1\right);2}^2(x;ϵ)\right]}^{1/2},$$

(13)

where $u_{\left(1\right);1}$ and $u_{\left(1\right);2}$ are two-linearly independent solutions of (3) for $\epsilon =ϵ$. The $\alpha$-function is real-valued and free of zeros in $\mathrm{Dom}{n}_0$ if the set $\{a,b,c\}$ is composited by positive numbers fulfilling

$$b^2-4ac=-4{\lambda }^2/W_0^2,$$

(14)

where $W_0=W(u_{\left(1\right);1},u_{\left(1\right);2})=const$ is the Wronskian of $u_{\left(1\right);1}$ and $u_{\left(1\right);2}$.

Using (13), the complex-valued superpotential acquires a simple form

$${\beta }_1(x;ϵ)=-\frac{1}{2}\frac{d}{dx}lnv(x;ϵ)+i\frac{\lambda }{v(x;ϵ)}=-\left[\frac{v^{\prime }(x;ϵ)-i2\lambda }{2v(x;ϵ)}\right],$$

(15)

where

$$v(x;ϵ)=a{u}_{\left(1\right);1}^2(x;ϵ)+b{u}_{\left(1\right);1}(x;ϵ)u_{\left(1\right);2}(x;ϵ)+c{u}_{\left(1\right);2}^2(x;ϵ).$$

(16)

Then, the Darboux transformation (7) gives the complex-valued refractive index

$$n_1(x;ϵ)=n_0\left(x\right)+\frac{\sqrt{2}}{k_0{n}_{*}}\frac{d}{dx}\left[\frac{v^{\prime }(x;ϵ)-i2\lambda }{2v(x;ϵ)}\right].$$

(17)

On the other hand, it may be shown that the imaginary part of $n_1$ satisfies the condition of zero total area [49]:

$${\int }_{\mathbb{R}}\mathrm{Im}{n}_1(x;ϵ)dx={\left.-\frac{\sqrt{2}\lambda }{k_0{n}_{*}}\frac{1}{v(x;ϵ)}\right|}_{-\infty }^{+\infty }=0,$$

(18)

so the total optical power is conserved. Equation (18) implies a balanced interplay between gain and loss that does not depend on any symmetry of either $Im{n}_1(x;ϵ)$ or $Re{n}_1(x;ϵ)$.

In contraposition to conventional supersymmetric approaches, the factorization energy $ϵ$ can be positioned at any place in the point spectrum of $n_1$ [41]. Moreover, the missing state (11), now written as

$$E_{\left(1\right)}^M(x;ϵ)=\frac{{\mathcal{N}}_{\left(1\right)}^M}{\sqrt{v(x;ϵ)}}exp\left[i\lambda \int v^{-1}(x;ϵ)dx\right],$$

(19)

is such that its real and imaginary parts are even and odd functions of the position-variable x, respectively.

We would like to emphasize that the nonlinear superposition (16) marks a distance with conventional supersymmetric approaches. Indeed, we have already shown [43] that the superpotential (15) can be also written in the conventional logarithmic form $\beta (x;ϵ)=-\frac{d}{dx}lnw(x;ϵ)$, where the coefficients of the linear superposition

$$w(x;ϵ)=a{u}_{\left(1\right),1}(x,ϵ)+\left(\frac{b}{2}-i\frac{\lambda }{W_0}\right)u_{\left(1\right),2}(x,ϵ)$$

(20)

are ruled by the constraint (14), with a and b complex numbers in general. Clearly, such a concrete combination of coefficients permits $n_1$ to satisfy the condition of zero total area (18), which defines it as a balanced gain-and-loss refractive index.

2.2.1. Bi-Orthogonality

The solutions of the paraxial Helmholtz Equation (9), with $n_1(x;ϵ)$ given in (17), are obtainable from (10), (15) and (19). However, while $E_{\left(1\right)}^M$ and all the TE modes $E_{\left(1\right)}$ are normalizable, they form a peculiar set since $E_{\left(1\right)}^M$ is orthogonal to all the $E_{\left(1\right)}$ but these last are not mutually orthogonal [39] (such property is not a problem in the Hermitian case since all the new functions satisfy the conventional oscillation theorems). Nevertheless, the eigenfunctions of $n_1$ satisfy some properties of interlacing of zeros that permit the study of the related systems as if they were Hermitian [49]. In this context, the bi-orthogonal set formed by the eigenstates $E_{\left(1\right)}$ of $n_1$, together with the eigenstates ${\tilde{E}}_{\left(1\right)}$ of the complex-conjugated refractive index, written $n_1^C$, provide an extended space of states where all the basis elements are bi-orthonormal [39,40,42]. Indeed, the bi-product

$$\left({\tilde{E}}_{\left(1\right);m},E_{\left(1\right);n}\right)={\int }_{\mathbb{R}}{\tilde{E}}_{\left(1\right);m}^C(x;ϵ)E_{\left(1\right);n}(x;ϵ)dx$$

(21)

is equal to zero if $n\ne m$, and serves to define the bi-norm $\left|\right|E_{\left(1\right);n}{\left|\right|}_B=\left|\right|{\tilde{E}}_{\left(1\right);n}{\left|\right|}_B$ if $n=m$ [39]. With two possible normalizations at our disposal, $E_{\left(1\right)}/\left|\right|E_{\left(1\right)}\left|\right|$ and $E_{\left(1\right)}/\left|\right|E_{\left(1\right)}{\left|\right|}_B$, it is important to emphasize that the real and imaginary parts of the modes $E_{\left(1\right)}$ behave qualitatively equal in both normalizations [42], although their bi-normalized values are usually larger than those obtained with the conventional normalization. Nevertheless, the differences become negligible as the excitation of the TE mode increases, see [42] for details. Note also that the notions of bi-product and bi-norm introduced above coincide with the conventional definitions if $\lambda =0$.

2.2.2. PT-Symmetric Case

The expression (17) represents a wide family of balanced gain-and-loss refractive indices. A very interesting subset of such family is integrated by the so-called parity-time (PT) symmetric refractive indices. Recalling that invariance under parity and time-reversal transformations requires $n\left(x\right)=n^C(-x)$ in quantum mechanics [35], we realize that the initial refractive index $n_0\left(x\right)$ should be parity-invariant $n_0\left(x\right)=n_0(-x)$ to facilitate the construction of PT-symmetric refractive indices $n_1(x;ϵ)$. On the other hand, assuming real-valued transformation functions $u_{\left(1\right);1}$ and $u_{\left(1\right);2}$, we may take $b=0$ in (16) to obtain the quadratic form

$$v_{PT}(x;ϵ)=a{u}_{\left(1\right);1}^2(x;ϵ)+c{u}_{\left(1\right);2}^2(x;ϵ).$$

(22)

The straightforward calculation shows that using this function in (17) one obtains a complex-valued graded index that is PT-symmetric.

Please note that $b=0$ implies $ac={\lambda }^2/W_0^2$ in (14). To simplify notation, without loss of generality, let us make $a=c=\lambda /|W_0|$ in (22), where $|W_0|$ stands for the modulus of $W_0$. Then (17) yields

$$\begin{array}{c}{n}_{PT}(x;ϵ)=n_0\left(x\right)+\frac{1}{\sqrt{2}{k}_0{n}_{*}}\frac{d^2}{d{x}^2}ln\left[u_{\left(1\right);1}^2(x;ϵ)+u_{\left(1\right);2}^2(x;ϵ)\right]\hfill \\ \hfill -i\frac{\sqrt{2}\left|W_0\right|}{k_0{n}_{*}}\frac{d}{dx}{\left[u_{\left(1\right);1}^2(x;ϵ)+u_{\left(1\right);2}^2(x;ϵ)\right]}^{-1}.\end{array}$$

(23)

Observe that the w-configuration (20) leads to the same result provided $a=-i\lambda /W_0^C$.

2.2.3. Recovering the Real-Valued Case

As indicated above, if $\lambda =0$ the superpotential (15) is reduced to its Hermitian configuration, which produces real-valued indices only. Revisiting the constraint (14) we see that $\lambda =0$ implies $b^2=4ac$, and thus $b=\pm 2\sqrt{ac}$. We obtain the linear superpositions ${\alpha }_{\pm }=\sqrt{a}{u}_{\left(1\right);1}\pm \sqrt{c}{u}_{\left(1\right);2}$, so we arrive at the conventional superpotentials

$${\beta }_R(x;ϵ;\pm )=-\frac{d}{dx}ln\left[\sqrt{a}{u}_{\left(1\right);1}(x;ϵ)\pm \sqrt{c}{u}_{\left(1\right);2}(x;ϵ)\right],$$

(24)

where a and c are such that $\beta$ is free of singularities in $Dom{n}_0$. Therefore, from (7) one has the two-parametric family of real-valued graded indices

$$n_R(x;ϵ;\pm )=n_0\left(x\right)+\frac{\sqrt{2}}{k_0{n}_{*}}\frac{d^2}{d{x}^2}ln\left[\sqrt{a}{u}_{\left(1\right);1}(x;ϵ)\pm \sqrt{c}{u}_{\left(1\right);2}(x;ϵ)\right].$$

(25)

3. Applications

The method developed in previous sections may be applied to practically any exactly solvable refractive index $n_0\left(x\right)$. The expression $n_k(x;ϵ)$, obtained at the kth step, actually represents a very wide family of new exactly solvable refractive indices with k additional propagation constants in their point spectrum with respect to $n_0\left(x\right)$. Even more important is the fact that $n_k(x;ϵ)$ can be constructed to be a real- or complex-valued function. In any case, the propagation constants belonging to the point spectrum of $n_k(x;ϵ)$ are all-real. Next, we provide a very useful example of the applicability of our approach. We use the mathematical solutions associated with $n_0=0$ to produce diverse families of cosh-like refractive indices admitting the presence of a given number of guided TE modes. The results include complex-valued refractive indices that are not limited to the parity-invariant case.

From now on, for the sake of simplicity, the expressions of the refractive index profiles $n_0\left(x\right)$ and $n_k(x;ϵ)$ are considered up to the additive constant $n_{*}$.

3.1. Adding Guided Modes One at a Time

The fundamental solutions of the paraxial Helmholtz equation for $n_0=0$ are well known. For positive factorization constants $ϵ=k^2>0$ we write $u_1=e^{ik(x-x_0)}$ and $u_2=e^{-ik(x-x_0)}$, with $W_0=-i2k$. However, the above expressions yield sinusoidal refractive indices $n_1$ [39], which are out of the scope of the present work. Here we make $k=i\kappa$ to obtain negative factorization constants $ϵ=-{\kappa }^2$, therefore

$$v(x;\kappa )=a{e}^{2\kappa (x-x_0)}+c{e}^{-2\kappa (x-x_0)}+b,b^2-4ac=-{\lambda }^2/{\kappa }^2.$$

(26)

To simplify notation let us make $a=c$. Then $v(x;\kappa )=2acosh\left[2\kappa (x-x_0)\right]+b$, with $b^2=4a^2-{\lambda }^2/{\kappa }^2$. The superpotential (15) acquires the form

$${\beta }_1(x;\kappa )=-\left[\frac{\kappa sinh\left[2\kappa (x-x_0)\right]-i\frac{\lambda }{2a}}{cosh\left[2\kappa (x-x_0)\right]+\frac{b}{2a}}\right],b^2=4a^2-\frac{{\lambda }^2}{{\kappa }^2},$$

(27)

so the refractive index (17) is in this case

$$n_1(x;\kappa )=\frac{{\left(2\kappa \right)}^2\left(1+\frac{b}{2a}cosh\left[2\kappa (x-x_0)\right]\right)+i\frac{\lambda }{a}\left(2\kappa \right)sinh\left[2\kappa (x-x_0)\right]}{\sqrt{2}{k}_0{n}_{*}{\left(cosh\left[2\kappa (x-x_0)\right]+\frac{b}{2a}\right)}^2}.$$

(28)

The complex-valued graded refractive index (28) allows the presence of only one guided TE mode, obtained from (11) in the form

$$E_{\left(1\right)}^M(x;\kappa )=\frac{{\mathcal{N}}_{\left(1\right)}^M}{\sqrt{cosh\left[2\kappa (x-x_0)\right]+\frac{b}{2a}}}exp\left\{-\frac{i}{4a}arctan\left(\left(\frac{b}{2a}-1\right)tanh\left[\kappa (x-x_0)\right]\right)\right\}.$$

(29)

Following the indications of the previous section, let us make $\lambda =0$ and $b=2a$ in (27) and (28) to obtain

$${\beta }_R(x;\kappa )=-\kappa tanh\left[\kappa (x-x_0)\right],n_R(x;\kappa )=\frac{2{\kappa }^2}{\sqrt{2}{k}_0{n}_{*}{cosh}^2\left[\kappa (x-x_0)\right]},$$

(30)

which are the well-known expressions for the cosh-like refractive index. The function $n_R(x;\kappa )$ is depicted in Figure 1a for a representative propagation constant $\epsilon$, which may be located at any position ${\epsilon }_1=-{\kappa }_1^2$ since it is the only one eigenvalue in the point spectrum.

For $b=0$ and $a=\frac{\lambda }{2\kappa }$ the formulae (27) and (28) give rise to the PT-symmetric expressions

$${\beta }_{PT}(x;\kappa )=-\left(\frac{\kappa sinh\left[2\kappa (x-x_0)\right]-i\kappa }{cosh\left[2\kappa (x-x_0)\right]}\right),n_{PT}(x;\kappa )=\frac{{\left(2\kappa \right)}^2\left(1+isinh\left[2\kappa (x-x_0)\right]\right)}{\sqrt{2}{k}_0{n}_{*}{cosh}^2\left[2\kappa (x-x_0)\right]}.$$

(31)

The PT-symmetric refractive index (31) is shown in Figure 1b for a representative propagation constant $\epsilon$. As in the previous case, this eigenvalue can be positioned at will in the negative part of the real axis.

If we repeat the procedure, assuming now that $n_1(x;{\kappa }_1)$ has been already fixed, the finite-difference algorithm provides an immediate superpotential (see details in Appendix A):

$${\beta }_2(x;{\kappa }_1,\kappa )=-{\beta }_1(x;{\kappa }_1)+\frac{\sqrt{2}{k}_0{n}_{*}({\kappa }_1^2-{\kappa }^2)}{{\beta }_1(x;{\kappa }_1)-{\beta }_1(x;\kappa )},$$

(32)

where ${\kappa }_1$ and ${\beta }_1(x;{\kappa }_1)$ have been fixed in the previous step. Deciding the concrete value of $\kappa$, as well as the analytical form of ${\beta }_1(x;\kappa )$ in (27), the above equation provides the new refractive index

$$n_2(x;{\kappa }_1,\kappa )=n_1(x;{\kappa }_1)-\frac{\sqrt{2}}{k_0{n}_{*}}{\beta }_2^{\prime }(x;{\kappa }_1,\kappa )=-\frac{d}{dx}\left[\frac{2({\kappa }_1^2-{\kappa }^2)}{{\beta }_1(x;{\kappa }_1)-{\beta }_1(x;\kappa )}\right],$$

(33)

where we have used (7) with $n_0\left(x\right)=0$.

At the present stage, we have incorporated two propagation constants, so the point spectrum of $n_2(x;{\kappa }_1,\kappa )$ is composited by the eigenvalues ${\epsilon }_1=-{\kappa }_1^2$ and $\epsilon =-{\kappa }^2$. However, some caution is necessary if the first step was addressed to produce $n_R(x;{\kappa }_1)$ and we are looking for a second real-valued refractive index $n_R(x;{\kappa }_1,\kappa )$. In such a case the inequality $\epsilon <{\epsilon }_1$ must be satisfied to obtain regular functions $n_R(x;{\kappa }_1,\kappa )$. Moreover, in such case, it may be shown [19] that it is better to combine the two different real-valued superpotentials ${\beta }_R(x;\kappa ;\pm )$. The case “+” is reported in Equation (30), the case “−” corresponds to the complementary expression ${\beta }_R(x;\kappa ;-)=-\kappa coth\left[\kappa (x-x_0)\right]$, see Section 2.2.3. We therefore arrive at the real-valued graded index

$$n_R(x;{\kappa }_1,\kappa )=\frac{2({\kappa }_1^2-{\kappa }^2)\left({\kappa }_1^2{csch}^2\left[{\kappa }_1(x-x_1)\right]+{\kappa }^2{sech}^2\left[\kappa (x-x_0)\right]\right)}{{\left(-{\kappa }_{1}coth\left[{\kappa }_1(x-x_1)\right]+\kappa tanh\left[\kappa (x-x_0)\right]\right)}^2}.$$

(34)

The behavior of $n_R(x;{\kappa }_1,\kappa )$ is shown in Figure 2 for different spectra $\{\epsilon ,{\epsilon }_1\}$ and constants $x_0$ and $x_1$.

Remarkably, $\epsilon$ and ${\epsilon }_1$ characterize the global profile of the function (34). Indeed, for ${\kappa }_1\gg \kappa$ and $x_0=x_1=0$, the refractive index $n_R(x;{\kappa }_1,\kappa )$ acquires a bell-shaped form. However, a valley arises at the top of such curve if $\kappa ={\kappa }_1-\varrho$, with $0\le \varrho \ll 1$. The dent is more pronounced as $\varrho \to 0$, separating the initial bell-like curve into a pair of bell-shaped ones. At the very limit, the new curves have moved in opposite directions toward the domain edges $\pm \infty$. Quite interestingly, actual waveguides are manufactured by including such dent, “sometimes for reducing the internal mechanical stress due to the gradient of dopant concentration, and sometimes for reducing the multimode dispersion” [1], p. 83. With this in mind, Figure 3 shows the exploration of the parameters that characterize $n_R(x;{\kappa }_1,\kappa )$, addressed to produce different dent configurations in the refractive index. These may be completely symmetrical as in Figure 3a or asymmetrical, as shown in Figure 3b. For ${\kappa }_1>\kappa$, local deformations may be produced by tuning the displacement parameters $x_0$ and $x_1$, see Figure 3c.

The ordering problem suffered by the propagation constants in the construction of $n_R(x;{\kappa }_1,\kappa )$ is easily avoided by considering any superpotential (27) with $\lambda \ne 0$ in either of the two steps. For instance, as in the previous example, assume that $n_R(x;{\kappa }_1)$ has been fixed in the first step. To include the second eigenvalue $\epsilon$ this time we use the complex-valued superpotential ${\beta }_{PT}(x;\kappa )$ introduced in (31). The new refractive index (33) is now complex-valued, given by

$$n_2(x;{\kappa }_1,\kappa )=\frac{2({\kappa }_1^2-{\kappa }^2)f(x;{\kappa }_1,\kappa )}{g^2(x;{\kappa }_1,\kappa )},$$

(35)

where

$$f(x;{\kappa }_1,\kappa )=-{\kappa }_1^2{sech}^2\left({\kappa }_{1}x\right)+2{\kappa }^2{sech}^2\left(2\kappa x\right)+i2{\kappa }^{2}tanh\left(2\kappa x\right)sech\left(2\kappa x\right)$$

(36)

and

$$g(x;{\kappa }_1,\kappa )=-{\kappa }_{1}tanh\left({\kappa }_{1}x\right)+\kappa tanh\left(2\kappa x\right)-i\kappa sech\left(2\kappa x\right).$$

(37)

In (36) and (37) we have omitted the displacement constants $x_0$ and $x_1$ for the sake of simplicity.

As $n_R(x;{\kappa }_1)$ is parity-invariant $n_R(x;{\kappa }_1)=n_R(-x;{\kappa }_1)$, the parameters of $n_2(x;{\kappa }_1,\kappa )$ in (35) can be managed to obtain a PT-symmetric refractive index $n_{PT}(x;{\kappa }_1,\kappa )$. The result is shown in Figure 4a for the process in which we add first ${\epsilon }_1$ and then $\epsilon$, with ${\epsilon }_1>\epsilon$. The reversed process is shown in Figure 4b. Please note that although the profile of $n_{PT}(x;{\kappa }_1,\kappa )$ changes, the PT symmetry is preserved under the change ${\epsilon }_1\leftrightarrow \epsilon$. The same expression (35) gives rise to refractive indices that are not invariant under parity-time transformations, as exhibited in Figure 5.

We have already mentioned that the procedure may be repeated at will. At the kth step, the method provides a set of superpotentials ${\beta }_k$ that are available for the finite-difference algorithm, addressed to elaborate the step $k+1$. The case considered in this section takes the null function $n_0\left(x\right)=0$ as the initial refractive index. The propagation constants are added one at a time to arrive at the point spectrum $\{{\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_{k-1},\epsilon \}$, which may be decided under prescription. The refractive indices constructed in this form admit k guided TE modes, generated from the initial missing state $E_{\left(1\right)}^M(x;\epsilon )$, via the rule (10). These modes obey the bi-product introduced in Section 2.2.1, which also defines a proper bi-norm that coincides with the conventional norm if $\lambda =0$.

3.2. Manipulating a Set of Guided Modes at Once

To complete the revision of immediate applications, consider the cosh-like refractive index

$$n_0(x,m)=\frac{m(m+1){\kappa }^2}{{cosh}^2\left(\kappa x\right)},m=1,2,\dots$$

(38)

Potentials $V\left(x\right)=-n_0(x,m)$ form the subset of transparent Pöschl–Teller potentials in quantum mechanics. The solutions of the Schrödinger equation for the entire family are well known [18,50,51], including resonances and anti-bound states [52,53]. Our interest in $n_0(x,m)$ obeys the fact that this refractive index admits exactly m guided TE modes, defined by the quadratic rule [18,50,51]

$${\epsilon }_{m,\ell }=-{\kappa }^2{(m-\ell )}^2,\ell =0,1,\dots ,m-1,m=fixed.$$

(39)

The finite-difference algorithm will provide k additional eigenvalues at the kth iteration, so the spectrum of $n_k(x,m;\kappa )$ is composited by two finite subsets $\left\{{\epsilon }_i\right\}\cup \left\{{\epsilon }_{m,\ell }\right\}$, with $i=1,2,\dots ,k$. As we have shown in the previous section, depending on the 1-step superpotentials ${\beta }_1(x,m;\epsilon )$ and the factorization constants $\epsilon$, the new eigenvalues ${\epsilon }_i$ may be positioned at arbitrary places of the initial spectrum $\left\{{\epsilon }_{m,\ell }\right\}$.

The fundamental basis of solutions is in this case provided by the functions [18,50]

$$u_1(x;\kappa )={(cosh\kappa x)}^{m+1}{}_2{F}_1(a,b,\frac{1}{2};-{sinh}^2\kappa x),$$

(40)

and

$$u_2(x;\kappa )=(sinh\kappa x){(cosh\kappa x)}^{m+1}{}_2{F}_1(a+\frac{1}{2},b+\frac{1}{2},\frac{3}{2};-{sinh}^2\kappa x),$$

(41)

where

$$a=\frac{m+1}{2}-\frac{\sqrt{\left|\epsilon \right|}}{2\kappa },b=\frac{m+1}{2}+\frac{\sqrt{\left|\epsilon \right|}}{2\kappa }.$$

(42)

To construct the complex-valued superpotential (15), a first function $v(x;\kappa )$ is easily achieved by noticing that the hypergeometric function ${}_2{F}_1$ is reduced to the identity if $a=0$. From (42) we immediately realize that $\left|\epsilon \right|={\kappa }^2{(m+1)}^2$ produces such a result. Remarkably, the latter value is in correspondence with the spectral rule (39) if $\ell =-1$. Thus, we are in position of adding the eigenvalue $\epsilon ={\epsilon }_{m,-1}$ to the initial spectrum $\left\{{\epsilon }_{m,\ell }\right\}$. The resulting refractive index $n_1(x,m;\kappa )$, obtained from Equation (17), may be chosen to be either real, complex-valued or PT-symmetric.

In Figure 6 we have depicted the case in which $n_1(x,m;\kappa )$ exhibits PT symmetry. In Figure 6a we started with $n_0(x,1)$, which admits only one guided TE mode, the one associated with ${\epsilon }_{1,0}$. The spectrum of the resulting refractive index $n_1(x,1;\kappa )$ is therefore integrated by ${\epsilon }_M=-{\left(2\kappa \right)}^2$ and ${\epsilon }_{1,0}=-{\kappa }^2$. Figure 6b considers the initial spectrum ${\epsilon }_{2,0}=-{\left(2\kappa \right)}^2$, ${\epsilon }_{2,1}=-{\kappa }^2$, and includes the missing state ${\epsilon }_M=-{\left(3\kappa \right)}^2$. Similarly for Figure 6c. The configuration where the new refractive index is not PT-symmetric is shown in Figure 7.

4. Discussion of Results and Conclusions

We have provided new exactly solvable models for optical waveguiding. Applying the supersymmetric finite-difference algorithm [19], we have generated a wide family of refractive indices whose point spectrum can be designed under prescription. The family includes refractive indices in both the real- and complex-valued configurations, the latter admitting all-real eigenvalues (propagation constants) in their point spectrum. We have shown that the spectral distribution may be organized in arbitrary form if it is constructed adding one at a time the eigenvalues such that complex-valued superpotentials are included. The result is relevant since such a property seems to be unnoticed in optical supersymmetry until the present work, although we have already reported this possibility in quantum mechanics [41].

One of the main results presented in this work shows that the index profile strongly depends on the factorization constants that are incorporated. In particular, adding two of them, either in a single step or in a twice iterated movement, one may produce a dent over the top of the profile that is in complete agreement with actual manufacture of refractive indices [1], see Figure 3. The phenomenon is not exclusive of the real-valued indices so produced since it is also admissible in the complex-valued case for the real part of the PT-symmetric refractive indices, see Figure 2. Considering that “in the manufacture and evaluation of optical fibers, the measurement of the index profile is one of the most important steps” [1], our results may be useful to analyze the data obtained from such measurements.

Another of our results shows that refractive indices admitting a given number of guided TE modes, such as the sech-like ones, can be deformed to admit an additional guided mode, the propagation of which can be positioned anywhere in the point spectrum of the initial refractive index. In addition, the new indices are not required to be PT-symmetric to allow all-real eigenvalues in their point spectrum. The transition from these results to the time-dependent case is straightforward [54], where PT-symmetric structures find interesting applications [55,56].

We have addressed the investigation to obtain guided TE modes in the new waveguiding-structures. This is the reason for which we started from initial refractive indices admitting no leaky modes. In previous works we have studied such possibility by analyzing the resonances of the initial structure [16,47,57]. It is viable to construct the supersymmetric partners using resonances of the initial system [31], a technique implemented also in the cosh-like case [52,53] and for soliton-like models [58]. However, the transformation of resonances and/or using resonances is elaborated, so it will be analyzed elsewhere. An important point to notice is that although generated from transparent refractive indices, the new structures presented here lack this property as a consequence of the non-Hermiticity (the clear exception is the real-valued case, since it is well known that supersymmetry leaves transparency invariant for such systems). Insights on the matter have been presented by the PT symmetry community and will be considered for future progress of our model.