1. Introduction

Supersymmetric Quantum Mechanics (SUSY-QM) was first introduced as a toy model to study SUSY breaking [1], and has since evolved into a powerful theoretical framework with many applications in theoretical physics [2–7]. Such application include the classification of potentials with identical spectra [8], analytical derivation of spectra in certain classes of quantum Hamiltonians [9], extensions to the WKB approximation [10], analysis of the Fokker-Planck equation (FPE) in statistical mechanics [7], and random matrix theory [11], disorder, and chaos [12]. Using the close analogy between the Schrödinger equation for a particle in a potential and light propagation in a dielectric medium, SUSY-QM has also been employed extensively in various optical waveguide designs and optical scattering problems over the past few years [13–25].

Another application of SUSY-QM is in the analysis of classical stochastic dynamics with one Cartesian degree of freedom [7,26,27]. It uses the fact that the FPE can be mapped into an imaginary-time Schrödinger equation

2. Optical waveguide with micro-bending

An optical waveguide is a transversely inhomogeneous structure, usually containing a region of increased refractive index, compared with the surrounding medium. Here, to illustrate the main ideas, we consider the micro-bending of a planar waveguide structure. Generalization to an optical fiber is possible but is beyond the scope of the present work. We note that intrinsic fluctuations in the refractive index, including those stemming from random glass density fluctuations or the presence of impurities, lead to unwanted back-scattering and mode mixing. For the case of micro-bending, the refractive index fluctuation is not intrinsic but can be modeled as a modification of the refractive index by the addition of a linear term to the first order in the radial coordinate over the radius of the curvature. Such a tilt in the refractive index profile results in leaking radiation.

For the problem of interest, we assume a planar waveguide structure defined by an inhomogeneous refractive index $n(y)$, which is independent of the $x$ and $z$ coordinates. The light rays are confined in the $y$ direction by the inhomogeneous refractive index, while propagating freely in the $z$ direction. The trajectories of light rays are determined by Fermat’s principle [30], $\delta \int ^{b}_an(y)ds=0$, where $ds$ is the differential length on a trajectory between points $a$ and $b$ in the $y-z$ plane. Using the calculus of variation results in the following differential equation for the ray trajectory $y(z)$:

To model the micro-bending, we assume that the optical waveguide is bent in the $y$-$z$ plane with a local curvature of $\kappa (z)$. A bent waveguide can be conformally mapped to a straight waveguide, where the effect of the local bending can be modeled as a variation in the refractive index [31,32] given by $n(y,z)\approx n_{s}(y)\left [1-\,\kappa (z)\,y\right ]$. $n_{s}$ is the effective refractive index for the straight waveguide (without micro-bending). In the paraxial approximation, the ray trajectory is almost parallel to the $z$ axis, so we can assume $ds\approx dz$. Moreover, the maximum of the refractive index is commonly at or near the waveguide center at $y=0$, and the refractive index gradually decreases away from the center. Taking these points into consideration and in the absence of large transverse variations, the ray propagation equation can be simplified into the dynamical equation of a point particle subject to a net force of $F(y)-\kappa (z)$, where the conservative force $F(y)$ is exclusively due to the inhomogeneous refractive index $n_s(y)$ [33,34]:

Note that the longitudinal coordinate $z$ plays the role of time for the point particle. The conservative force $F(y)$ is related to the effective potential $U(y)$ by

We assume that the waveguide curvature $\kappa (z)$ is a zero mean ($\langle {\kappa (z)}\rangle =0$) white stochastic process with the autocorrelation given by

The ray propagation Eq. (3) can be expressed as a stochastic Langevin equation in the following form:

The FPE (7) describes many stochastic problems in the phase-space; therefore, the following results can be applied to a wide range of physical systems beyond the micro-bending problem that is analyzed here.

This can be further simplified by noting that the Hamiltonian of the conservative system (in the absence of stochastic noise) is given by

The phase-space action, which is a conserved quantity, is the area enclosed in the phase-space by the ray trajectory and is given as

Using this information, the FPE in (7) can be expressed in a much simpler form, where $y$ and $\theta$ are swapped with a single variable $E$ to obtain

We now make a change of variable (both in $E$ and $P$) and transform the FPE in (11) to an imaginary-time Schrödinger equation, similar to Eq. (1). We write $E=E(\mathcal {E})$ and $P(E,z)=B(E)Q(E,z)$ and choose the function $E(\mathcal {E})$ and the form of $B(E)$ to satisfy

Using these changes of variables, we arrive at the imaginary-time Schrödinger equation

We now have all the required ingredients to make a SUSY transformation on Eq. (14). Once we find the SUSY partner potential to $V(\mathcal {E})$, we can work our way backward to find the corresponding $I(E)$ from Eq. (15) and consequently $Z(E)$, from which we can find $U(y)$. This will be achieved by using the Abel transform-pairs, which can be expressed as

Note that Eq. (16) is essentially the same as Eq. (10), albeit with a change of variable. Equation (17) gives us the required information about the new $U(y)$. In other words, the Abel transform-pairs allow us to relate the phase-space action to the refractive index profile of the optical waveguide and vice versa.

Example of a monomial potenial: As a concrete example we explore a monomial potential of the form $U(y)=\Delta |y/y_c|^{\alpha }$ in some detail. This model arises in graded-index optical waveguides [35,36], where $y_c$ is the waveguide half-width, $\alpha$ controls the shape of the index profile, and $\Delta$ determines the refractive index contrast of the center relative to the core edges at $y=\pm y_c$. Note that it is common to take $\alpha \approx 2$ in practical optical waveguides [35,36], because it reduces the modal disperion in optical waveguides by equalizing the group velocities of the modes. In the ray picture, for $\alpha \approx 2$, all rays of different energy $E$ (and maximum value of $\theta$) reach the end of the waveguide at the same time.

In the following discussion, because the potential is assumed to be symmetric under $y\to -y$, we only consider the $y\ge 0$ region without loss of generality. For this problem, we can evaluate $I(E)$ using

In this special case, the ratio $\mathbb {Z}(E)$ takes the simple form of $\mathbb {Z}(E)=\zeta E$, and we obtain $4E=\zeta \mathcal {E}^{2}$ and $B^{2\alpha }\propto E^{(1-\alpha )}$. Using these results and Eq. (15), we obtain:

Before we present the solution to Eq. (14), we note that lossy rays are those which can breach the potential barrier into the free space, i.e. with an energy larger than a given $\mathbb {E}$. This results in a boundary condition of $P(\mathbb {E},z)=0$, or equivalently $Q(\boldsymbol {\epsilon },Z)=0$, where $4\mathbb {E}=\zeta \boldsymbol {\epsilon }^{2}$. Attempting an eigensolution of the form $Q(\mathcal {E})=Q_m(\mathcal {E})e^{-\lambda _m Z}$, we obtain

To perform a SUSY transformation, we need to calculate the superpotential $W=-\dot {Q}_0/Q_0$, where $Q_0$ is the eigensolution with the smallest decay coefficient, i.e. $\lambda _0$. We obtain

The partner potentials are derived from $V_{\pm } = W^{2} \pm \dot {W}$. We obtain

To establish a concrete numerical example, we consider the case of a quadratic graded index with $\alpha =2$ and $\nu =0$, where $u_{00}\approx 2.405$. We assume that those rays that can make it to $y\ge y_c$ are lost to free space, resulting in $\mathbb {E}=\Delta$ or equivalently $\boldsymbol {\epsilon }=2\sqrt {\Delta }$. We also assume $\Delta =0.01$, which is a reasonable value in graded index optics, resulting in $\boldsymbol {\epsilon }=0.2$ and $\lambda _0\approx 144.58$. In this case, we plot $V_{\pm }$ from Eqs. (22) and (23), except we redefine these functions by adding $\lambda _0$ to get the original $V_{-}$ back. These shifted potentials are plotted in FIG. 1. In physical units, the lowest attenuation coefficient corresponds to $\lambda _0\kappa _0^{2}/2$. Considering a waveguide with $y_c\approx 1$µm and $\kappa ^{2}_0\approx 4.8\times 10^{-6}$µm$^{-1}$, which corresponds to a fluctuation in radius of curvature of 457µm over a longitudinal distance comparable to $y_c$, the lowest attenuation is obtained to be around 3dB/mm. Less fluctuation corresponds to larger curvature radii fluctuating at a slower rate with the propagation distance.

 

Fig. 1. $V_{\pm }$ from Eqs. (22) and (23) are plotted as a function of $\sqrt {\lambda _{0}}\mathcal {E}$ for $\nu =0$. The potentials are shifted up by $\lambda _0$ and are expressed in units of $\lambda _0$.

Download Full Size | PDF

To verify that $V_+$ is a true superpartner potential to $V_-$, we solve the eigenvalue problem with $V_+$ in Eq. (14) assuming the $Z$-dependence of the form $e^{-\lambda _m Z}$ and the Dirichlet boundary condition: $Q^{+}(0)=Q^{+}(\boldsymbol {\epsilon })=0$ (superscript $+$ signifies its correspondence to $V_+$). The eigenvalues must match those of the excited states obtained from $V_-$. We plot, in FIG. 2(a), the lowest 4 eigenfunctions $Q^{-}_m$ ($m=0,1,2,3$) corresponding to $V_-$ from Eq. (20), where $m$ corresponds to how many times the eigenfunction crosses the $Q=0$ horizontal axis. The corresponding eigenvalues are $\lambda _m=(u_{0m}/u_{00})^{2}\lambda _0$ for $m=0,1,2,3$, respectively. Similarly, in FIG. 2(b), we plot the lowest 3 eigenfunctions $Q^{+}_m$ ($m=1,2,3$) corresponding to $V_+$ that are evaluated numerically, where $m-1$ corresponds to the number of crossings with the $Q=0$ horizontal axis. It can be confirmed, numerically, that the eigenvalues corresponding to $Q^{+}_1$, $Q^{+}_2$, and $Q^{+}_3$ match those of $Q^{-}_1$, $Q^{-}_2$, and $Q^{-}_3$, respectively, as expected from an unbroken SUSY.

 

Fig. 2. (a) The lowest 4 eigenfunctions $Q^{-}_m$ ($m=0,1,2,3$) corresponding to $V_-$, and (b) the lowest 3 eigenfunctions $Q^{+}_m$ ($m=1,2,3$) corresponding to $V_+$ are plotted.

Download Full Size | PDF

We now outline the procedure to evaluate $U(y)$ corresponding to $V_+$. We first insert $V_+$ from Eq. (23) into Eq. (15) and obtain a numerical solution for $\mathbb {F}$. The boundary conditions for $\mathbb {F}$ are set at $\mathcal {E}\to 0$: we make a Taylor expansion of $V_+$ near this point and analytically solve Eq. (15) near $\mathcal {E}=0$ to find that $\mathbb {F}(\mathcal {E})\propto \mathcal {E}^{3/2}$ and $\dot {\mathbb {F}}(\mathcal {E})\propto (3/2)\mathcal {E}^{1/2}$. These expressions guide us to set the boundary conditions for $\mathbb {F}$ and $\dot {\mathbb {F}}$ for $\mathcal {E}\to 0$. Note that Eq. (15) can determine $\mathbb {F}$ only up to an overall factor, and this is rooted in the fact that $I(E)$ in Eq. (9) can be multiplied by a constant factor without changing the FPE. The next step is to use $\sqrt {2}\dot {I}=\mathbb {F}^{2}$, assuming that $I(0)=0$, to numerically evaluate $I(\mathcal {E})$. It can be shown that $\mathbb {Z}(\mathcal {E})=2I^{2}/\mathbb {F}^{4}$, so we use this expression to evaluate $\mathbb {Z}(\mathcal {E})$. Next, $\mathbb {Z}(\mathcal {E})$ along with the transformations in Eq. (13) are used to relate $E$ to $\mathcal {E}$. Using this relationship and $I(\mathcal {E})$, we can find $I(E)$ and consequently $I^{\prime }(E)$, to be used in Eq. (17) to find $U(y)$. Considering that $U(y=0)=0$, we can rewrite Eq. (17) in the simpler form of

Because the redundant multiplicative factor mentioned earlier propagates to $I^{\prime }(E)$ in Eq. (24), the right-hand side of Eq. (24) can only be determined up to an overall factor. This ambiguity is there because nowhere in our formalism did we use the actual value of $y_c$, so we can normalize the result obtained in Eq. (24) to $y_c$. The final result is shown in Fig. 3, where $U_-$ corresponds to $V_-$ and is given by $\Delta (y/y_c)^{2}$, and $U_+$ corresponds to $V_+$ and is obtained according to the numerical procedure outlined above. We remind that because these potentials are symmetric under $y\to -y$, only the $y\ge 0$ region is shown in Fig. 3. Note that $U_-$ is quadratic and concave; however, $U_+$ is convex and its derivative at $y=0$ blows up. In practical situations, a slightly rounded potential at $y=0$ does not alter the behavior of the system in a notable way.

 

Fig. 3. (a) $U_-=\Delta (y/y_c)^{2}$ (corresponding to $V_-$) is plotted in comparison to $U_+$ (corresponding to $V_+$), where $U_+$ is obtained according to the outlined numerical procedure. Plots in (b) are the same as (a) but the vertical axis is logarithmic. Note that $V_{-,-}$ and $U_{-,-}$ do not exist for the case of $\alpha =2$ treated here.

Download Full Size | PDF

Let us recap the behavior of the two systems as $\mathcal {E}\to 0$. For $V_-$, $\mathbb {F}\propto \mathcal {E}^{1/2}$, so $\dot {I}\propto \mathcal {E}$ and $I\propto \mathcal {E}^{2}$. Therefore, $\mathbb {Z}\propto \mathcal {E}^{2}$, which results in $E\propto \mathcal {E}^{2}$. As such, $I\propto E$ and $I^{\prime }={\,const.}$, which results in $\partial _yU=0$ for $y\to 0$ according to Eq. (17). On the other hand, for $V_+$, $\mathbb {F}\propto \mathcal {E}^{3/2}$, so $\dot {I}\propto \mathcal {E}^{3}$ and $I\propto \mathcal {E}^{4}$. Therefore, $\mathbb {Z}\propto \mathcal {E}^{2}$, which results in $E\propto \mathcal {E}^{2}$. As such, $I\propto E^{2}$ and $I^{\prime }\propto E$, which results in $|\partial _yU|\to \infty$ for $|y|\to 0$ according to Eq. (17). The lowest order decay coefficient for $U_+$ is $\lambda _1=(u_{01}/u_{00})^{2}\lambda _0$, which is larger than that of $U_-$ ($\lambda _0$) and this agrees with the form of $U_+$ and the fact that it is shallower than $U_-$.

3. Lowering the loss

We would now like to carry out the reverse operation and find out if there exists a SUSY-constructed potential, which we refer to as $U_{-,-}$, that can result in a lower minimum attenuation than $\lambda _0$. To achieve this, one must find a potential $V_{-,-}$ that is the superpartner potential to (a slightly shifted) $V_{-}$ but with a lower (zero energy) ground state. This implies that we need to find a $W$ that satisfies $W^{2}+\dot {W}=V_-$, where the positive sign behind $\dot {W}$ implies that $V_-$ is the higher energy superpartner potential, and $W^{2}-\dot {W}=V_{-,-}$.

It is seen from Eq. (19) that regardless of a finite shift in $V_-$, the behavior of $V_-$ for $\mathcal {E}\to 0$ is of the form $(\nu ^{2}-1/4)/\mathcal {E}^{2}$. Previously, setting $W^{2}-\dot {W}=V_-$ resulted in $W\sim -(\nu +1/2)/\mathcal {E}$, consistent with Eq. (21) for $\mathcal {E}\to 0$ and $Q_0\sim \mathcal {E}^{1/2+\nu }$ with a Dirichlet boundary condition for $Q_0$ at $\mathcal {E}=0$. For $\alpha > 0$, $\nu >-1/2$, and hence $Q_0$ is always well-behaved for $\mathcal {E}\to 0$, which is reassuring. Now, setting $W^{2}+\dot {W}=V_-$ results in $W\sim (1/2-\nu )/\mathcal {E}$. As such, the ground state of $V_{-,-}$, which would have to satisfy $W=-\dot {Q}_0/Q_0$ gives $Q_0\sim \mathcal {E}^{-1/2+\nu }$. This form of $Q_0$ is only non-diverging for $\mathcal {E}\to 0$, and hence normalizable, if $\nu >1/2$. This happens if $\alpha <1$, which implies $V_{-,-}$ can be found only if $\nu >1/2$ (equivalently $\alpha <1$).

This clearly shows that for the quadratic potential $\Delta (y/y_c)^{2}$, with $\nu =0$ and $\alpha =2$, there exists no $V_{-,-}$. In other words, one cannot find a lower energy superpartner potential that results in a lower micro-bending attenuation. However, for $\Delta (y/y_c)^{\alpha }$ with $\alpha <1$ which is of the concave form, $V_{-,-}$ can be found. It can be used to construct a potential $U_{-,-}$, and hence a refractive index that has a lower micro-bending attenuation. The evidence of this already exists in our calculations. Previously, we used $\nu =0$ for $V_-$ in Eq. (22), corresponding to $U_-=\Delta (y/y_c)^{2}$, to construct $V_+$ in Eq. (23) corresponding to $U_+$. As such, one expects to be able to construct $V_-$ from $V_+$. To reverse the logic, we note that the behavior of $V_+$ for $\mathcal {E}\to 0$ is of the form $V_+\sim (3/4)/\mathcal {E}^{2}$ for $\mathcal {E}\to 0$. Because $\nu ^{2}-1/4=3/4$ results in $\nu =1>1/2$, so it is no wonder that $V_-$ exists as the lower energy partner of $V_+$.

Using these results, we can make general arguments on how the shape of $U$ changes under SUSY transformations, both going up ($V_-$ to $V_+$) and possibly going down ($V_-$ to $V_{-,-}$). Recall that $W^{2}-\dot {W}=(\nu ^{2}-1/4)/\mathcal {E}^{2}$ results in $W\sim -(\nu +1/2)/\mathcal {E}$, while $W^{2}+\dot {W}=(\nu ^{\prime 2}-1/4)/\mathcal {E}^{2}$ results in $W\sim (1/2-\nu ^{\prime })/\mathcal {E}$. Thus, for an “up-ward” SUSY transformation $\nu ^{\prime } = \nu + 1$ or, equivalently, $\alpha ^{\prime }=\alpha /(1+\alpha )$. We have already seen this for $\nu =0$ and $\alpha =2$ in $V_-$, which resulted in $\nu ^{\prime }=1$ and $\alpha ^{\prime }=2/3$ in $V_+$. In fact, it can be shown that $U_+$ scales as $\sim (y/y_c)^{2/3}$ near $y=0$. The relation $\alpha ^{\prime } = \alpha /(1 + \alpha )$ shows that in an “up-ward” SUSY transformation the power in the potential of the form $U \sim (y/y_c)^{\alpha }$ decreases and the potential becomes (more) concave. Another important point is that the “up-ward” SUSY transformation is always permitted [37]. On the other hand, for a “down-ward” SUSY transformation, i.e. for $V_-$ going down to $V_{-,-}$, we have $\alpha ^{\prime }=\alpha /(1-\alpha )$. First, $\alpha <1$ is required for this to make sense, as we showed above, and the potential becomes (more) convex under this operation. For example, $\alpha =1/2$ in $V_-$ gives $\alpha ^{\prime }=1$ for $V_{-,-}$.

We would like to emphasize two important points regarding the results obtained in this paper. First, we make a general observation that the attenuation can be lowered only in profiles with $\alpha <1$. In particular, for the interesting case of $\alpha \approx 2$, lowering the attenuation is not possible. However, our results are strictly true only for monomial index profiles and we have not explored other forms. Therefore, we are not ruling out the possibility that a SUSY transformation can be used to lower the attenuation for a non-monomial near-quadratic index profile. Second, our conclusions are based only on SUSY transformations and we are not ruling out the possibility of lowering the attenuation for the monomial index profile with $\alpha >1$ by leveraging other non-SUSY transformations. In practical optical waveguide designs (especially for optical fibers), some common values of $\alpha$ are $\alpha \gg 1$ for a step-index waveguide, $\alpha \approx 2$ for a graded-index waveguide to reduce the modal dispersion, and $\alpha \approx 1$ for a dispersion-shifted waveguide [38]. One can envision $\alpha <1$ profiles for dispersion shifting as well. However, considering that $\alpha <1$ is not a commonly used profile, our finding can be considered as a “negative result”, showing that the attenuation of commonly used optical waveguides cannot be reduced by a SUSY transformation, subject to the limitations stated above.

4. Conclusions

We have employed SUSY-QM in the analysis of a classical stochastic optics problem: the random micro-bending loss of a waveguide. For a general class of monomial-shaped potentials (refractive index profile), we showed that a SUSY transformation can always be used to make the waveguide more lossy ($V_-$ to $V_+$), but only certain refractive index profiles are amenable to a (reverse) SUSY transformation to make the waveguide less lossy ($V_-$ to $V_{-,-}$). The results of this paper can potentially be used as a seed for the optimization or inverse design of optical waveguides. Our work may also serve as a framework to investigate similar phenomena in optics or other areas of physics described by a FPE that can be investigated by using the SUSY-QM formalism.

We would like to emphasize that our work is performed using the geometrical ray picture in what effectively constitutes a highly multimode waveguide. We would like to acknowledge the interesting result obtained in Ref. [39] on the loss due to sidewall roughness for single mode propagation. Their attenuation is proportional to the standard deviation of the autocorrelation function, just as our calculated attenuation is proportional to $\kappa ^{2}_0$. However, they find that the attenuation is also proportional to the wavelength, which is not the case in our highly multimode ray-based picture. Another point that differentiates our work is our emphasis on refractive index-based optimization strategies to reduce the attenuation.