Supersymmetric Analysis of Stochastic Micro-Bending in Optical Waveguides
SSupersymmetric Analysis of Stochastic Micro-Bending in Optical Waveguides ∗ Stuart Ward , , Rouzbeh Allahverdi , and Arash Mafi , † Department of Physics & Astronomy, University of New Mexico, Albuquerque, NM 87106, USA. Center for High Technology Materials, University of New Mexico, Albuquerque, New Mexico 87106, USA. (Dated: October 2, 2020)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 equivalentto a change in the refractive index profile, resulting in a larger or smaller attenuation. For a broadclass of monomial index profiles, it is always possible to obtain an index profile with a largermicro-bending attenuation using a supersymmetric transformation. However, obtaining a smallerattenuation is not always possible and is restricted to a subset of index profiles.
Introduction:
Supersymmetric Quantum Mechanics(SUSY-QM) was first introduced as a toy model to studySUSY breaking [1], and has since evolved into a powerfultheoretical framework with many applications in theoret-ical physics [2–7]. Such application include the classifi-cation of potentials with identical spectra [8], analyticalderivation of spectra in certain classes of quantum Hamil-tonians [9], extensions to the WKB approximation [10],analysis of the Fokker-Planck equation (FPE) in statis-tical mechanics [7], and random matrix theory [11], dis-order, and chaos [12]. Using the close analogy betweenthe Schr¨odinger equation for a particle in a potential andlight propagation in a dielectric medium, SUSY-QM hasalso been employed extensively in various optical waveg-uide designs and optical scattering problems over the pastfew years [13–26].Another application of SUSY-QM is in the analysis ofclassical stochastic dynamics with one Cartesian degreeof freedom [7, 27, 28]. It uses the fact that the FPE canbe mapped into an imaginary-time Schr¨odinger equation − ∂∂t Ψ ± ( t, x ) = H ± Ψ ± ( t, x ) , (1)where H ± = − ∂ x + W ( x ) ± W (cid:48) ( x ) are the partner Hamil-tonians in the Witten’s model, and W ( x ) is the super-potential. The eigenvalues of H ± are the decay ratesassociated with the eigenfunctions of the imaginary-timeSchr¨odinger equation. In the presence of an unbrokenSUSY, H + and H − have identical positive eigenvalues(positive decay rates), while H − has an additional van-ishing eigenvalue, corresponding to its ground-state. Inthis Letter, we explore how this formalism can be used tosolve practical problems in classical stochastic optics. Inparticular, we consider the problem of micro-bending inoptical waveguides, which is an important issue affectingthe optical communications networks [29]. Undesirablefiber-fiber and fiber-cable interactions, caused by man-ufacturing of fiber cables or thermal-induced variations,result in random microscopic bends leading to signal loss. ∗ This research is supported by grant number W911NF-19-1-0352from the United States Army Research Office. † mafi@unm.edu This Letter is framed around using SUSY-QM to designan optical waveguide with a reduced micro-bending loss;however, the overall formalism and ideas can be readilyapplied to a wide range of problems involving stochasticdynamics.
Optical waveguide with micro-bending:
An optical waveg-uide is a transversely inhomogeneous structure, usuallycontaining a region of increased refractive index, com-pared with the surrounding medium. Here, to illustratethe main ideas, we consider the micro-bending of a planarwaveguide structure. Generalization to an optical fiber ispossible but is beyond the scope of the present work. Forthe problem of interest, we assume a planar waveguidestructure 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 inho-mogeneous refractive index, while propagating freely inthe z direction. The trajectories of light rays are deter-mined by Fermat’s principle [30], δ (cid:82) ba n ( 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 variationresults in the following differential equation for the raytrajectory y ( z ): dds (cid:18) n dyds (cid:19) = ∂n∂y , ds = (cid:112) dy + dz . (2)To model the micro-bending, we assume that the op-tical waveguide is bent in the y - z plane with a local cur-vature of κ ( z ). A bent waveguide can be conformallymapped to a straight waveguide, where the effect of thelocal bending can be modeled as a variation in the refrac-tive index [31, 32] given by n ( y, z ) ≈ n s ( y ) [1 − κ ( z ) y ]. n s is the effective refractive index for the straight waveg-uide (without micro-bending). In the paraxial approxi-mation, the ray trajectory is almost parallel to the z axis,so we can assume ds ≈ dz . Moreover, the maximum ofthe refractive index is commonly at or near the waveg-uide center at y = 0, and the refractive index graduallydecreases away from the center. Taking these points intoconsideration and in the absence of large transverse vari-ations, the ray propagation equation can be simplifiedinto the dynamical equation of a point particle subjectto a net force of F ( y ) − κ ( z ), where the conservative force F ( y ) is exclusively due to the inhomogeneous refractive a r X i v : . [ phy s i c s . op ti c s ] O c t index n s ( y ) [33, 34]: d ydz ≈ F ( y ) − κ ( z ) . (3)Note that the longitudinal coordinate z plays the role oftime for the point particle. The conservative force F ( y )is related to the effective potential U ( y ) by F ( y ) = − dU/dy, U ( y ) = 1 − n s ( y ) /n , (4)where n = n s ( y = 0) and U (0) = 0. Because the refrac-tive index decreases away from the center of the waveg-uide, the potential increases at large | y | and U ( y ) is aconfining potential.We assume that the waveguide curvature κ ( z ) is a zeromean ( (cid:104) κ ( z ) (cid:105) = 0) white stochastic process with the au-tocorrelation given by (cid:104) κ ( z ) κ ( z (cid:48) ) (cid:105) = κ δ ( z − z (cid:48) ) . (5)The ray propagation equation (3) can be expressed as astochastic Langevin equation in the following form: dydz = θ ( z ) , dθdz = F ( y ) − κ ( z ) , (6)where θ ( z ) is the local angle of the ray relative to the z axis. Instead of solving Eq. (6) for each random realiza-tion and then averaging the results, it is more convenientto use the corresponding FPE, from which we can deter-mine the probability P of an optical ray to be at lateralposition y and angle θ at point z along the fiber [33, 34]: (cid:20) ∂∂z + θ ∂∂y + F ∂∂θ − κ ∂ ∂θ (cid:21) P ( y, θ, z ) = 0 . (7)The FPE (7) describes many stochastic problems in thephase-space; therefore, the following results can be ap-plied to a wide range of physical systems beyond themicro-bending problem that is analyzed here.This can be further simplified by noting that theHamiltonian of the conservative system (in the absenceof stochastic noise) is given by E = 12 θ + U ( y ) , (8)where θ plays the role of the velocity of the point particlewith the mass m = 1. The phase-space action, which isa conserved quantity, is the area enclosed in the phase-space by the ray trajectory and is given as I ( E ) = (cid:90) (cid:90) dθdy = (cid:73) T (cid:112) E − U ( y )] dy = T θ , (9)where T is the ray period and overbar denotes the averageover one ray period. The derivative of I with respect to E also takes the simple form of I (cid:48) ( E ) = (cid:73) T (cid:112) E − U ) dy = T. (10) Using this information, the FPE in (7) can be expressedin a much simpler form, where y and θ are swapped witha single variable E to obtain2 κ ∂P ( E, z ) ∂z = − ∂P ( E, z ) ∂E + ∂ [ Z ( E ) P ( E, z )] ∂E , (11)where P ( E, z ) is the probability of the ray to have en-ergy E at point z along the waveguide, and Z ( E ) = I ( E ) /I (cid:48) ( E ) is twice the average kinetic energy over oneray period.We now make a change of variable (both in E and P ) and transform the FPE in (11) to an imaginary-time Schr¨odinger equation, similar to Eq. (1). We write E = E ( E ) and P ( E, z ) = B ( E ) Q ( E, z ) and choose thefunction E ( E ) and the form of B ( E ) to satisfy Z E (cid:48) = 1 , Z = ˙ E , Z B (cid:48) + (cid:0) Z (cid:48) − (cid:1) B = 0 , (12)where (cid:48) and · denote differentiation with respect to E and E , respectively. The solutions to these equations areformally given by E = (cid:90) dE √ Z , E = (cid:90) √ Z d E , (13) B = Z − / exp (cid:16) (cid:90) dE Z (cid:17) . Using these changes of variables, we arrive at theimaginary-time Schr¨odinger equation − ∂Q∂Z = − ∂ Q∂ E + V ( E ) Q, (14)where Z = k z/ V ( E ) = F − ∂ ∂ E F , F = (cid:2) I ( E ) (cid:3) (cid:48) = 2 (cid:2) ˙ I ( E ) (cid:3) . (15)We now have all the required ingredients to make aSUSY transformation on Eq. (14). Once we find theSUSY partner potential to V ( E ), we can work our waybackward 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 I (cid:48) ( E ) = (cid:90) E (cid:32) √ ∂ y U (cid:33) dU √ E − U , (16)2 √ ∂ y U = 1 π ∂∂U (cid:90) U I (cid:48) ( E ) dE √ U − E . (17)Note that Eq. (16) is essentially the same as Eq. (10),albeit with a change of variable. Equation (17) gives usthe required information about the new U ( y ). Example of a monomial potenial:
As a concrete exam-ple we explore a monomial potential of the form U ( y ) =∆ | y/y c | α in some detail. This model arises in graded-index optical waveguides [35, 36], where y c is the waveg-uide half-width, α controls the shape of the index profile,and ∆ determines the refractive index contrast of thecenter relative to the core edges at y = ± y c . Note thatit is common to take α ≈ y → − y , weonly consider the y ≥ I ( E ) using I ( E ) = (cid:90) E (cid:32) √ ∂ y U (cid:33) √ E − U dU (18)= √ π Γ(1 + 1 /α )Γ(3 / /α ) y c ∆ − /α E /ζ , where ζ = 2 αα + 2 . In this special case, the ratio Z ( E ) takes the simple formof Z ( E ) = ζE , and we obtain 4 E = ζ E and B α ∝ E (1 − α ) . Using these results and Eq. (15), we obtain: V ( E ) = 1 − αα E = ν − / E , where ν = 2 − α α . (19)Before we present the solution to Eq. (14), we notethat lossy rays are those which can breach the poten-tial barrier into the free space, i.e. with an energy largerthan a given E . This results in a boundary conditionof P ( E , z ) = 0, or equivalently Q ( (cid:15) , Z ) = 0, where4 E = ζ (cid:15) . Attempting an eigensolution of the form Q ( E ) = Q m ( E ) e − λ m Z , we obtain Q m ( E ) = (cid:112) E / (cid:15) J ν ( (cid:112) λ m E ) , λ m = u νm / (cid:15) , (20)where J ν is the Bessel function of order ν and u νm is the m th zero of J ν , where m = 0 , , · · · .To perform a SUSY transformation, we need to calcu-late the superpotential W = − ˙ Q /Q , where Q is theeigensolution with the smallest decay coefficient, i.e. λ .We obtain W = − ν + 1 / E + (cid:112) λ J ν +1 ( √ λ E ) J ν ( √ λ E ) . (21)The partner potentials are derived from V ± = W ± ˙ W .We obtain V − = − λ + ν − / E , (22) V + = + λ + GE , (23) G = 34 + ν ( ν + 2) − (4 ν + 2) √ λ E J ν +1 J ν + 2 λ E J ν +1 J ν , where we have suppressed the argument of the Besselfunctions, all being √ λ E . Note that V − is the same asthe potential in Eq. (19) except it is shifted by λ asexpected [37]. We also have λ = u ν / (cid:15) ; therefore, weare only interested in V ± in the range of 0 ≤ E < (cid:15) , where V + blows up at the upper boundary. To establish a concrete numerical example, we considerthe case of a quadratic graded index with α = 2 and ν = 0, where u ≈ . y ≥ y c are lost to free space, resultingin E = ∆ or equivalently (cid:15) = 2 √ ∆. We also assume∆ = 0 .
01, which is a reasonable value in graded indexoptics, resulting in (cid:15) = 0 . λ ≈ .
58. In this case,we plot V ± from Eqs. (22) and (23), except we redefinethese functions by adding λ to get the original V − back.These shifted potentials are plotted in FIG. 1. V + V - - - λ ℰ V - / λ a nd V + / λ FIG. 1. V ± from Eqs. (22) and (23) are plotted as a functionof √ λ E for ν = 0. The potentials are shifted up by λ andare expressed in units of λ . 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 − λ m Z and theDirichlet boundary condition: Q + (0) = Q + ( (cid:15) ) = 0 (su-perscript + signifies its correspondence to V + ). Theeigenvalues must match those of the excited states ob-tained from V − . We plot, in FIG. 2(a), the lowest 4eigenfunctions Q − m ( m = 0 , , ,
3) corresponding to V − from Eq. (20), where m corresponds to how many timesthe eigenfunction crosses the Q = 0 horizontal axis. Thecorresponding eigenvalues are λ m = ( u m /u ) λ for m = 0 , , ,
3, respectively. Similarly, in FIG. 2(b), weplot the lowest 3 eigenfunctions Q + m ( m = 1 , ,
3) cor-responding to V + that are evaluated numerically, where m − 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 − , Q − , and Q − , respectively, as ex-pected from an unbroken SUSY.We now outline the procedure to evaluate U ( y ) corre-sponding to V + . We first insert V + from Eq. (23) intoEq. (15) and obtain a numerical solution for F . Theboundary conditions for F are set at E →
0: we make aTaylor expansion of V + near this point and analyticallysolve Eq. (15) near E = 0 to find that F ( E ) ∝ E / and˙ F ( E ) ∝ (3 / E / . These expressions guide us to set theboundary conditions for F and ˙ F for E →
0. Note thatEq. (15) can determine F only up to an overall factor,and this is rooted in the fact that I ( E ) in Eq. (9) canbe multiplied by a constant factor without changing theFPE. The next step is to use √ I = F , assuming that I (0) = 0, to numerically evaluate I ( E ). It can be shown ( a ) λ ℰ ( b ) λ ℰ FIG. 2. (a) The lowest 4 eigenfunctions Q − m ( m = 0 , , , V − , and (b) the lowest 3 eigenfunctions Q + m ( m = 1 , ,
3) corresponding to V + are plotted. that Z ( E ) = 2 I / F , so we use this expression to evalu-ate Z ( E ). Next, Z ( E ) along with the transformations inEq. (13) are used to relate E to E . Using this relation-ship and I ( E ), we can find I ( E ) and consequently I (cid:48) ( E ),to be used in Eq. (17) to find U ( y ). Considering that U ( y = 0) = 0, we can rewrite Eq. (17) in the simplerform of y = 12 √ π (cid:90) U I (cid:48) ( E ) dE √ U − E . (24) U - U + ( a ) / y c U ( y ) U - U + ( b ) - - - / y c FIG. 3. (a) U − = ∆( y/y c ) (corresponding to V − ) is plottedin comparison to U + (corresponding to V + ), where U + is ob-tained according to the outlined numerical procedure. Plotsin (b) are the same as (a) but the vertical axis is logarithmic. Because the redundant multiplicative factor mentionedearlier propagates to I (cid:48) ( E ) in Eq. (24), the right-handside of Eq. (24) can only be determined up to an overallfactor. This ambiguity is there because nowhere in ourformalism did we use the actual value of y c , so we cannormalize the result obtained in Eq. (24) to y c . Thefinal result is shown in Fig. 3, where U − corresponds to V − and is given by ∆( y/y c ) , and U + corresponds to V + and is obtained according to the numerical procedureoutlined above. We remind that because these potentials are symmetric under y → − y , only the y ≥ U − is quadratic and concave;however, U + is convex and its derivative at y = 0 blowsup. In practical situations, a slightly rounded potentialat y = 0 does not alter the behavior of the system in anotable way.Let us recap the behavior of the two systems as E → V − , F ∝ E / , so ˙ I ∝ E and I ∝ E . Therefore, Z ∝E , which results in E ∝ E . As such, I ∝ E and I (cid:48) = const. , which results in ∂ y U = 0 for y → V + , F ∝ E / , so ˙ I ∝ E and I ∝ E . Therefore, Z ∝ E , which results in E ∝ E .As such, I ∝ E and I (cid:48) ∝ E , which results in | ∂ y U | → ∞ for | y | → U + is λ = ( u /u ) λ , which is largerthan that of U − ( λ ) and this agrees with the form of U + and the fact that it is shallower than U − . Lowering the loss:
We would now like to carry out thereverse operation and find out if there exists a SUSY-constructed potential, which we refer to as U − , − , thatcan result in a lower minimum attenuation than λ . Toachieve this, one must find a potential V − , − that is thesuperpartner potential to (a slightly shifted) V − but witha lower (zero energy) ground state. This implies that weneed to find a W that satisfies W + ˙ W = V − , wherethe positive sign behind ˙ W implies that V − is the higherenergy superpartner potential, and W − ˙ W = V − , − .It is seen from Eq. (19) that regardless of a finite shiftin V − , the behavior of V − for E → ν − / / E . Previously, setting W − ˙ W = V − re-sulted in W ∼ − ( ν + 1 / / E , consistent with Eq. (21) for E → Q ∼ E / ν with a Dirichlet boundary con-dition for Q at E = 0. For α > ν > − /
2, and hence Q is always well-behaved for E →
0, which is reassuring.Now, setting W + ˙ W = V − results in W ∼ (1 / − ν ) / E .As such, the ground state of V − , − , which would have tosatisfy W = − ˙ Q /Q gives Q ∼ E − / ν . This form of Q is only non-diverging for E →
0, and hence normal-izable, if ν > /
2. This happens if α <
1, which implies V − , − can be found only if ν > / α < y/y c ) , with ν = 0 and α = 2, there exists no V − , − .In other words, one cannot find a lower energy super-partner potential that results in a lower micro-bendingattenuation. However, for ∆( y/y c ) α with α < V − , − can be found. It can be usedto construct a potential U − , − , and hence a refractive in-dex that has a lower micro-bending attenuation. Theevidence of this already exists in our calculations. Previ-ously, we used ν = 0 for V − in Eq. (22), corresponding to U − = ∆( y/y c ) , to construct V + in Eq. (23) correspond-ing to U + . As such, one expects to be able to construct V − from V + . To reverse the logic, we note that the be-havior of V + for E → V + ∼ (3 / / E for E →
0. Because ν − / / ν = 1 > /
2, soit is no wonder that V − exists as the lower energy partnerof V + .Using these results, we can make general argumentson how the shape of U changes under SUSY transforma-tions, both going up ( V − to V + ) and possibly going down( V − to V − , − ). Recall that W − ˙ W = ( ν − / / E resultsin W ∼ − ( ν +1 / / E , while W + ˙ W = ( ν (cid:48) − / / E re-sults in W ∼ (1 / − ν (cid:48) ) / E . Thus, for an “up-ward” SUSYtransformation ν (cid:48) = ν +1 or, equivalently, α (cid:48) = α/ (1+ α ).We have already seen this for ν = 0 and α = 2 in V − , which resulted in ν (cid:48) = 1 and α (cid:48) = 2 / V + . Infact, it can be shown that U + scales as ∼ ( y/y c ) / near y = 0. The relation α (cid:48) = α/ (1 + α ) shows that in an “up-ward” SUSY transformation the power in the potentialof the form U ∼ ( y/y c ) α decreases and the potential be-comes (more) concave. Another important point is thatthe “up-ward” SUSY transformation is always permit-ted [37]. On the other hand, for a “down-ward” SUSYtransformation, i.e. for V − going down to V − , − , we have α (cid:48) = α/ (1 − α ). First, α < α = 1 / V − gives α (cid:48) = 1 for V − , − . Conclusion:
We have employed SUSY-QM in the analy-sis of a classical stochastic optics problem: the randommicro-bending loss of a waveguide. For a general classof monomial-shaped potentials (refractive index profile),we showed that a SUSY transformation can always beused to make the waveguide more lossy ( V − to V + ), butonly certain refractive index profiles are amenable to a(reverse) SUSY transformation to make the waveguideless lossy ( V − to V − , − ). We emphasize that this pro-cedure may not be the most practical way to addressthe micro-bending problem in real-world scenarios and aSUSY transformation is a subset of all possible deforma-tions one can apply to a waveguide. However, the resultscan potentially be used as a seed for the optimization orinverse design of optical waveguides. Our work may alsoserve as a framework to investigate similar phenomena inoptics or other areas of physics described by a FPE thatcan be investigated by using the SUSY-QM formalism. [1] E. Witten, “Dynamical breaking of supersymmetry,”Nucl. Phys. B , 513–554 (1981).[2] E. Witten, “Supersymmetry and Morse theory,” J. Diff.Geom. , 661–692 (1982).[3] F. Cooper, A. Khare and U. Sukhatme, “Supersymme-try and Quantum Mechanics,” Phys. Rept. , 267–385(1995).[4] F. Cooper and B. Freedman, “Aspects of supersymmetricquantum mechanics,” Ann. Phys. , 262–288 (1983).[5] F. Cooper, A. Khare, and U. Sukhatme, Supersymme-try in quantum mechanics (World Scientific, Singapore,2001).[6] S. M. Chumakov and K. B. Wolf, “Supersymmetry inHelmholtz optics,” Phys. Lett. A , 51–53 (1994).[7] G. Junker,
Supersymmetric methods in quantum and sta-tistical physics, (Springer-Verlag, Berlin, 2012).[8] M. A. Jafarizadeh and H. Fakhri, “Supersymmetry andshape invariance in differential equations of mathematicalphysics,” Phys. Lett. A , 164–170 (1997).[9] R. Dutt, A. Khare, and U. P. Sukhatme, “Supersymme-try, shape invariance, and exactly solvable potentials,”Am. J. Phys. , 163–168 (1988).[10] R. Dutt, A. Khare, and U. P. Sukhatme,“Supersymmetry-inspired WKB approximation inquantum mechanics,” Am. J. Phys. , 723–727 (1991).[11] M. R. Zirnbauer, “The supersymmetry method of ran-dom matrix theory,” arXiv preprint math-ph/0404057(2004).[12] K. Efetov, Supersymmetry in disorder and chaos, (Cam-bridge University Press, New York, 1999).[13] R. El-Ganainy, K. G. Makris, and D. N. Christodoulides,“Local PT invariance and supersymmetric parametric os-cillators,” Phys. Rev. A , 033813 (2012)[14] M. A. Miri, M. Heinrich, R. El-Ganainy, and D. N. Christodoulides, “Supersymmetric optical struc-tures,” Phys. Rev. Lett. , 233902 (2013).[15] M. A. Miri, M. Heinrich, and D. N. Christodoulides,“Supersymmetry-generated complex optical potentialswith real spectra,” Phys. Rev. A , 043819 (2013).[16] M. Heinrich, M. A. Miri, S. St¨utzer, R. El-Ganainy,S. Nolte, A. Szameit, and D. N. Christodoulides, “Su-persymmetric mode converters,” Nat. Commun. , 3698(2014).[17] M. Heinrich, M. A. Miri, S. St¨utzer, R. El-Ganainy,S. Nolte, A. Szameit, and D. N. Christodoulides, “Su-persymmetric mode converters,” Nat. Commun. , 3698(2014).[18] M. Heinrich, M. A. Miri, S. St¨utzer, S. Nolte,D. N. Christodoulides, and A. Szameit, “Observationof supersymmetric scattering in photonic lattices,” Opt.Lett. , 6130–6133 (2014).[19] M. A. Miri, M. Heinrich, and D. N. Christodoulides,“SUSY-inspired one-dimensional transformation optics,”Optica , 89–95 (2014).[20] R. El-Ganainy, L. Ge, M. Khajavikhan, andD. N. Christodoulides, “Supersymmetric laser arrays,”Phys. Rev. A , 033818 (2015).[21] W. Walasik, B. Midya, L. Feng, and N. M. Litchinitser,“Supersymmetry-guided method for mode selection andoptimization in coupled systems,” Opt. Lett. , 3758–3761 (2018).[22] B. Midya, W. Walasik, N. M. Litchinitser, and L. Feng,“Supercharge optical arrays,” Opt. Lett. , 4927–4930(2018).[23] B. Midya, H. Zhao, X. Qiao, P. Miao, W. Walasik,Z. Zhang, N. M. Litchinitser, and L. Feng, “Supersym-metric microring laser arrays,” Photon. Res. , 363–367(2019) [24] W. Walasik, N. Chandra, B. Midya, L. Feng, andN. M. Litchinitser, “Mode-sorter design using continu-ous supersymmetric transformation,” Opt. Express ,22429–22438 (2019)[25] Q. Zhong, S. Nelson, M. Khajavikhan,D. N. Christodoulides, and R. El-Ganainy, “Bosonic dis-crete supersymmetry for quasi-two-dimensional opticalarrays,” Photon. Res. , 1240–1243 (2019)[26] M. P. Hokmabadi, N. S. Nye, R. El-Ganainy,D. N. Christodoulides, and M. Khajavikhan, “Supersym-metric laser arrays,” Science (6427), 623–626 (2019).[27] G. Parisi, and N. Sourlas, “Random magnetic fields,supersymmetry, and negative dimensions,” Phys. Rev.Lett. , 744 (1979).[28] G. Parisi, and N. Sourlas, “Supersymmetric field theoriesand stochastic differential equations,” Nucl. Phys. B ,321–332 (1982).[29] W. B. Gardner, “Microbending loss in optical fibers,”Bell Labs Tech. J. , 457–465 (1975).[30] B. E. A. Saleh and M. C. Teich, Fundamentals of pho-tonics, (John Wiley & sons, New Jersey, 2019). [31] M. Heiblum and J. Harris, “Analysis of curved opti-cal waveguides by conformal transformation,” IEEE J.Quantum Electron. , 75–83 (1975).[32] R. A. Herrera, C. T. Law, and A. Mafi, “Calculation ofthe acousto-optic coupling coefficients in optical fibers,”Opt. Commun. , 217–220 (2013).[33] M. Rousseau and J. Arnaud, “Ray theory of microbend-ing,” Opt. Commun. , 333–336 (1978).[34] J. Arnaud and M. Rousseau, “Ray theory of randomlybent multimode optical fibers,” Opt. Lett. , 63–65(1978).[35] D. Gloge and E. A. J. Marcatili, “Multimode theory ofgraded-core fibers,” Bell. Syst. Tech. J. , 1563–1578(1973).[36] A. Mafi, “Pulse propagation in a short nonlinear graded-index multimode optical fiber,” J. Lightw. Technol.18