Supersymmetric Analysis of Stochastic Micro-Bending in Optical Waveguides

Micro-bending attenuation in an optical waveguide can be modeled by a Fokker-Planck equation. It is shown that a supersymmetric transformation applied to the Fokker-Planck equation is equivalent to a change in the refractive index profile, resulting in a larger or smaller attenuation. For a broad class of monomial index profiles, it is always possible to obtain an index profile with a larger micro-bending attenuation using a supersymmetric transformation. However, obtaining a smaller attenuation is not always possible and is restricted to a subset of index profiles.


Micro-bending attenuation in an optical waveguide can be modeled by a Fokker-Planck equation.
It is shown that a supersymmetric transformation applied to the Fokker-Planck equation is equivalent to a change in the refractive index profile, resulting in a larger or smaller attenuation. For a broad class of monomial index profiles, it is always possible to obtain an index profile with a larger micro-bending attenuation using a supersymmetric transformation. However, obtaining a smaller attenuation is not always possible and is restricted to a subset of index profiles.
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 where H ± = −∂ 2 x +W 2 (x)±W (x) are the partner Hamiltonians in the Witten's model, and W (x) is the superpotential. The eigenvalues of H ± are the decay rates associated with the eigenfunctions of the imaginary-time Schrödinger equation. In the presence of an unbroken SUSY, H + and H − have identical positive eigenvalues (positive decay rates), while H − has an additional vanishing eigenvalue, corresponding to its ground-state. In this paper, we explore how this formalism can be used to solve practical problems in classical stochastic optics. In particular, we consider the problem of micro-bending in optical waveguides, which is an important issue affecting the optical communications networks [28,29]. Undesirable * This research is supported by grant number W911NF-19-1-0352 from the United States Army Research Office. † mafi@unm.edu fiber-fiber and fiber-cable interactions, caused by manufacturing of fiber cables or thermal-induced variations, result in random microscopic bends leading to signal loss. The estimation of microbending loss has been performed in various publications. In particular, Jin and Payne [29] have presented a detailed numerical investigation, based on an analytical model, and have validated the model by comparing with experimental measurements. Therefore, there is a good coverage on the estimation of microbending loss in the existing literature. Our intention here is not to present a new way of estimating microbending attenuation; rather, this paper is framed around using SUSY-QM to design an optical waveguide with a reduced micro-bending loss. Moreover, the overall formalism and ideas can be readily applied to a wide range of problems involving stochastic dynamics.

II. 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 microbending 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], δ b a n(y)ds = 0, arXiv:2009.11847v3 [physics.optics] 6 Aug 2021 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): d ds n dy ds = ∂n ∂y , ds = dy 2 + dz 2 . (2) To model the micro-bending, we assume that the optical waveguide is bent in the y-z plane with a local curvature of κ(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) ≈ n s (y) [1 − κ(z) y]. 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 ≈ 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)−κ(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 where n 0 = n s (y = 0) and U (0) = 0. Because the refractive index decreases away from the center of the waveguide, the potential increases at large |y| and U (y) is a confining potential. We assume that the waveguide curvature κ(z) is a zero mean ( κ(z) = 0) white stochastic process with the autocorrelation given by The ray propagation equation (3) can be expressed as a stochastic Langevin equation in the following form: where θ(z) is the local angle of the ray relative to the z axis. Instead of solving Eq. (6) for each random realization and then averaging the results, it is more convenient to use the corresponding FPE, from which we can determine the probability P of an optical ray to be at lateral position y and angle θ at point z along the fiber [33,34]: 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 where θ plays the role of the velocity of the point particle with the mass m = 1. In a typical graded-index optical waveguide, where the refractive index peaks in the center of the waveguide and gradually decreases away from the center, light rays follow sinusoidal-like paths along the waveguide. In the center where U is zero, the angle θ and the kinetic term assume their maximum values. As a ray reaches its maximum transverse separation from the center and bends back towards the center, θ and the kinetic term become zero and U reaches its maximum value equal to E. In other words, E determines the maximum ray angle at the center and also the maximum distance a ray can reach away from the center. Of course, if the energy E of a light ray is too large, it can breach the potential barrier into free space. In the modal language of wave optics, the energy E can be related to the mode number, where a larger E corresponds to a higher order mode and the energy at which a light ray breaches the potential barrier corresponds to the highest order mode supported by the waveguide. In Eq. (7), the microbending attenuation appears as the exponential decay of the probability P (y, θ, z) with respect to the propagation distance z as will be discussed in detail later in the paper. 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 where T is the ray period and overbar denotes the average over one ray period. The derivative of I with respect to E also takes the simple form of Using this information, the FPE in (7) can be expressed in a much simpler form, where y and θ are swapped with a single variable E to obtain where P (E, z) is the probability of the ray to have energy E at point z along the waveguide, and Z(E) = I(E)/I (E) is twice the average kinetic energy over one ray period.
We now make a change of variable (both in E and P ) and transform the FPE in (11) to an imaginarytime Schrödinger equation, similar to Eq. (1). We write E = E(E) and P (E, z) = B(E)Q(E, z) and choose the function E(E) and the form of B(E) to satisfy (12) where and · denote differentiation with respect to E and E, respectively. The solutions to these equations are formally given by Using these changes of variables, we arrive at the imaginary-time Schrödinger equation where Z = κ 2 0 z/2 and We now have all the required ingredients to make a SUSY transformation on Eq. (14). Once we find the SUSY partner potential to V (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) = ∆|y/y c | α in some detail. This model arises in gradedindex optical waveguides [35,36], where y c is the waveguide half-width, α controls the shape of the index profile, and ∆ determines the refractive index contrast of the center relative to the core edges at y = ±y c . Note that it is common to take α ≈ 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 α ≈ 2, all rays of different energy E (and maximum value of θ) reach the end of the waveguide at the same time.
In the following discussion, because the potential is assumed to be symmetric under y → −y, we only consider the y ≥ 0 region without loss of generality. For this problem, we can evaluate I(E) using In this special case, the ratio Z(E) takes the simple form of Z(E) = ζE, and we obtain 4E = ζE 2 and B 2α ∝ E (1−α) . 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 E. This results in a boundary condition of P (E, z) = 0, or equivalently Q( , Z) = 0, where 4E = ζ 2 . Attempting an eigensolution of the form where J ν is the Bessel function of order ν and u νm is the mth zero of J ν , where m = 0, 1, · · · . Recall that Z = κ 2 0 z/2, where κ 0 is defined through the autocorrelation of the curvature which is defined in Eq. (5). The relationship between attenuation and curvature can now be seen from the term e −λmZ of the above eigensolution. Furthermore, recall that the goal of performing the SUSY transformation is to find a new potential which possesses an additional eigensolution with a decay constant that is smaller than the smallest decay constant (λ 0 ) of the above eigensolution.
To perform a SUSY transformation, we need to calculate the superpotential W = −Q 0 /Q 0 , where Q 0 is the eigensolution with the smallest decay coefficient, i.e. λ 0 . We obtain The partner potentials are derived from V ± = W 2 ±Ẇ . We obtain where we have suppressed the argument of the Bessel functions, all being √ λ 0 E. Note that V − is the same as the potential in Eq. (19) except it is shifted by λ 0 as expected [37]. We also have λ 0 = u 2 ν0 / 2 ; therefore, we are only interested in V ± in the range of 0 ≤ E < , where V + blows up at the upper boundary.
To establish a concrete numerical example, we consider the case of a quadratic graded index with α = 2 and ν = 0, where u 00 ≈ 2.405. We assume that those rays that can make it to y ≥ y c are lost to free space, resulting in E = ∆ or equivalently = 2 √ ∆. We also assume ∆ = 0.01, which is a reasonable value in graded index optics, resulting in = 0.2 and λ 0 ≈ 144.58. In this case, we plot V ± from Eqs. (22) and (23) ficient corresponds to λ 0 κ 2 0 /2. Considering a waveguide with y c ≈ 1µm and κ 2 0 ≈ 4.8 × 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.
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 −λmZ and the Dirichlet boundary condition: Q + (0) = Q + ( ) = 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) match those of Q − 1 , Q − 2 , and Q − 3 , respectively, as expected from an unbroken SUSY.
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 F. The boundary conditions for F are set at E → 0: we make a Taylor expansion of V + near this point and analytically solve Eq. (15) near E = 0 to find that F(E) ∝ E 3/2 anḋ F(E) ∝ (3/2)E 1/2 . These expressions guide us to set the boundary conditions for F andḞ for E → 0. Note that Eq. (15) can determine 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 √ 2İ = F 2 , assuming that I(0) = 0, to numerically evaluate I(E). It can be shown that Z(E) = 2I 2 /F 4 , so we use this expression to evaluate Z(E). Next, Z(E) along with the transformations in Eq. (13) are used to relate E to E. Using this relationship and I(E), we can find I(E) and consequently I (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 (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 ∆(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 → −y, only the y ≥ 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.
Let us recap the behavior of the two systems as E → 0. For V − , F ∝ E 1/2 , soİ ∝ E and I ∝ E 2 . Therefore, Z ∝ E 2 , which results in E ∝ E 2 . As such, I ∝ E and I = const., which results in ∂ y U = 0 for y → 0 according to Eq. (17). On the other hand, for V + , F ∝ E 3/2 , soİ ∝ E 3 and I ∝ E 4 . Therefore, Z ∝ E 2 , which results in E ∝ E 2 . As such, I ∝ E 2 and I ∝ E, which results in |∂ y U | → ∞ for |y| → 0 according to Eq. (17). The lowest order decay coefficient for U + is λ 1 = (u 01 /u 00 ) 2 λ 0 , which is larger than that of U − (λ 0 ) and this agrees with the form of U + and the fact that it is shallower than U − .

III. 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 λ 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 +Ẇ = V − , where the positive sign behinḋ W implies that V − is the higher energy superpartner potential, and It is seen from Eq. (19) that regardless of a finite shift in V − , the behavior of V − for E → 0 is of the form (ν 2 − 1/4)/E 2 . Previously, setting W 2 −Ẇ = V − resulted in W ∼ −(ν + 1/2)/E, consistent with Eq. (21) for E → 0 and Q 0 ∼ E 1/2+ν with a Dirichlet boundary condition for Q 0 at E = 0. For α > 0, ν > −1/2, and hence Q 0 is always well-behaved for E → 0, which is reassuring. Now, setting W 2 +Ẇ = V − results in W ∼ (1/2 − ν)/E. As such, the ground state of V −,− , which would have to satisfy W = −Q 0 /Q 0 gives Q 0 ∼ E −1/2+ν . This form of Q 0 is only non-diverging for E → 0, and hence normalizable, if ν > 1/2. This happens if α < 1, which implies V −,− can be found only if ν > 1/2 (equivalently α < 1).
This clearly shows that for the quadratic potential ∆(y/y c ) 2 , with ν = 0 and α = 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 ∆(y/y c ) α with α < 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 ν = 0 for V − in Eq. (22), corresponding to U − = ∆(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 E → 0 is of the form V + ∼ (3/4)/E 2 for E → 0. Because ν 2 −1/4 = 3/4 results in ν = 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 Thus, for an "up-ward" SUSY transformation ν = ν +1 or, equivalently, α = α/(1+α). We have already seen this for ν = 0 and α = 2 in V − , which resulted in ν = 1 and α = 2/3 in V + . In fact, it can be shown that U + scales as ∼ (y/y c ) 2/3 near y = 0. The relation α = α/(1 + α) shows that in an "upward" SUSY transformation the power in the potential of the form U ∼ (y/y c ) α 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 α = α/(1 − α). First, α < 1 is required for this to make sense, as we showed above, and the potential becomes (more) convex under this operation. For example, 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 α < 1. In particular, for the interesting case of α ≈ 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 α > 1 by leveraging other non-SUSY transformations. In practical optical waveguide designs (especially for optical fibers), some common values of α are α 1 for a step-index waveguide, α ≈ 2 for a graded-index waveguide to reduce the modal dispersion, and α ≈ 1 for a dispersionshifted waveguide [38]. One can envision α < 1 profiles for dispersion shifting as well. However, considering that α < 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.

IV. CONCLUSIONS
We have employed SUSY-QM in the analysis of a classical stochastic optics problem: the random microbending 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 (re-verse) 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 κ 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. Funding. Grant number W911NF-19-1-0352 from the United States Army Research Office.