Determination of modes of elliptical waveguides with ellipse transformation perturbation theory

All optical fibers, including single-mode, multimode, and multicore fibers, exhibit some degree of birefringence, either purposely, such as in polarization-maintaining fibers, or inadvertently due to material or fabrication imperfections. Finding a low-complexity method to accurately calculate the modal characteristics of elliptical fibers has been a long-standing problem. We present a novel accurate perturbative method that avoids the difficulties associated with the traditional Mathieu function treatment. The method is also applicable to a broader class of oscillating systems with elliptical geometry. © 2017 Optical Society of America

The accurate determination of transverse modes of elliptical waveguides is important for both deliberately elliptical waveguides, such as polarization-maintaining fibers, but also in the context of circular fibers due to core shape fluctuations caused by manufacturing imperfections [1].Recently, multicore fibers with elliptical cores have been shown to be promising candidates for next-generation space-division multiplexing applications in optical communication networks [2].Traditionally, elliptical waveguides have been treated using the elliptical coordinate system that invariably leads to Mathieu differential equations [3].The normal modes are then expressed as linear combinations of products of angular and radial Mathieu functions, both of which are notoriously difficult to compute numerically [4,5].The angular Mathieu functions are expressed as Fourier series and the radial functions as series of products of Bessel functions, and the coefficients for the terms are obtained by calculating the eigenvectors of an infinite matrix.
The determination of transverse modes of step-index waveguides requires the matching of tangential fields at the corecladding interface.For circular waveguides, this is quite trivial and boils down to matching an oscillating Bessel function J m in the core with a decaying modified Bessel function K m in the cladding to ensure continuity of the tangential fields.For elliptical waveguides, this procedure is significantly more cumbersome as it requires the matching of an infinite series of products of Mathieu functions in the core with an infinite series of products of different Mathieu functions in the cladding.For this reason, numerous alternative approaches have been suggested.Schneider and Marquardt [6] used a power series to compute Mathieu functions, but the series becomes slowly converging for near-circular geometry.Nor does this approach avoid the tedious matching of two infinite series.Other methods [7][8][9][10][11][12][13][14][15] rely on heavy approximations and assumptions, such as weak guidance, and are therefore incapable of providing arbitrarily accurate results.
In this Letter we outline a method, dubbed the ellipse transformation perturbation method (ETPM) [16], to determine the modes of elliptical waveguides.ETPM is based on treating ellipticity as a perturbation from the circular case and therefore works especially well for near-circular geometry.The method is significantly simpler than the Mathieu function treatment, more versatile than previous approximate methods, and even in its coarsest form equally or more accurate than any past approximate method.Unlike previous approximate methods, ETPM is not restricted to small relative refractive index difference between the core and the cladding.Furthermore, ETPM is generally applicable to any oscillating physical system with elliptical geometry [17][18][19][20].
The refractive index profile of a waveguide with an elliptical cross section oriented along the z axis can be written as Now, instead of switching to the conventional elliptical cylindrical coordinates [3], we apply ETPM and re-scale the x axis by defining w b∕ax.In the new coordinates w; y; z, the elliptical core-cladding interface becomes circular-like as it takes the form w 2 y 2 b 2 .
To determine the waveguide modes, we solve for the longitudinal electric and magnetic field components E z F w; ye iβz−ωt and H z Gw; ye iβz−ωt .In the new coordinates w; y; z, The Helmholtz equation satisfied by E z becomes where we have defined δ b∕a 2 − 1.The parameter δ serves as a measure of the ellipticity of the waveguide core, vanishing for a circular waveguide.Note that δ can be negative or positive depending on whether the major axis of the ellipse is parallel to the x or the y axis, and either orientation is fine as long as jδj < 1.
By treating δ as a perturbation from the circular case, F (and G) can be expanded as without loss of generality.The condition jδj < 1 is needed for the series to converge.For Eq. ( 3) to satisfy Eq. ( 2), the equations have to be satisfied individually for each order of δ.This is analogous to perturbation theory commonly used in quantum mechanics with the key difference being that the perturbation affects the coordinate system used.Substituting Eq. (3) into Eq.( 2) yields an infinite set of equations: and where Eq. ( 5) was divided by δ n and holds for all n ≥ 1. Equation ( 4) is mathematically identical to the Helmholtz equation for a circular fiber, with the general solution given by where Z n is a Bessel function (J n in the core region and K n in the cladding region), γ is a constant related to k and β (with different values in the core and cladding regions), and a n and ϕ n are constants found by matching the boundary conditions.We used w r cos ϕ and y r sin ϕ, but r and ϕ are not the circular polar coordinates (a constant r corresponds to an ellipse instead of a circle).We refer to them as quasi-polar coordinates (QPCs).
In physical terms, the QPCs are polar coordinates after scaling the x axis by the axis ratio of the ellipse.The deviations of the field profile F from the circular case are given by the higher order terms F n obtained by solving Eq. ( 5), which is a driven Helmholtz equation.However, note that for δ ≠ 0 the zeroth-order term F 0 is already different from the circular waveguide modes as it has been scaled in the x dimension.The homogeneous part of Eq. ( 5) is identical to Eq. ( 4), and its general solution has the exact same form as the solution of Eq. ( 4) given in Eq. ( 6).The general form of F in Eq. ( 3) can thus be written as where F p n r; ϕ are particular solutions of Eq. ( 5).Even though an elliptical waveguide is not rotationally symmetric, it exhibits reflection symmetry with respect to both the x and y axes.As a result, F must be an even or odd function with respect to both x and y.Therefore, either ϕ n 0 or ϕ n π∕2 for all n.Furthermore, c n 0 for even n or c n 0 for odd n.Since sinx π∕2 cosx, the ϕ n 0 case will be referred to as the sine case and ϕ n π∕2 as the cosine case.In what follows, the treatment will be restricted to the sine case only.The treatment of the cosine case is identical with sines replaced by cosines.Equation (7) provides the general solution for the waveguide modes as an infinite series in powers of the perturbation parameter δ.To solve Eq. ( 5) for the first-order correction (n 1), we note that F 0 needs to converge toward the circular-waveguide solution in the limit δ → 0, say F 0 A sinnϕZ n γr.This means that for each mode the sum in Eq. ( 6) consists of only one term.Substituting this simple form of F 0 into Eq.( 5) then gives the equation for F 1 , the particular solution of which, after some algebra, is found to be where p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k 0 n 1 2 − β 2 p in the core, and where q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi β 2 − k 0 n 2 2 p in the cladding.The obtained result makes intuitive sense.When a circular-waveguide mode of angular dependence sinnϕ is perturbed by making the core slightly elliptical, the first-order correction comes from the sinn 2ϕ terms as they are the nearest-order non-zero perturbative terms.In addition to these particular solution terms, the sum in Eq. ( 7) also has terms of the form Z n2 γr sinn 2ϕ {or only the Z n2 γr sinn 2ϕ term if n < 2 since the sum in Eq. ( 7) starts from n 0} that originate from the homogeneous parts of the equations for F n .Since F 1 already has terms of the form sinn 2ϕ from the particular solution, we will interpret the Z n2 γr sinn 2ϕ terms that exist in Eq. ( 7) as belonging to F 1 , and other terms of the form Z n2 γr sinn 2mϕ that are present in the sum as being a part of F m .Even though the function F k affects the particular solution of F k1 , it does not matter which homogeneous solution terms are or are not included in F k , as the coefficients c n in Eq. ( 7) are arbitrary (for the time being) and the particular solutions give only the magnitudes of the corrective terms with respect to these coefficients.Without specifying which homogeneous solution terms belong to which order F n , only the zeroth-order term F 0 would be uniquely defined, aside from its normalization factor, as it has to coincide with the circular-waveguide mode for δ 0.
Since circular-waveguide modes have an azimuthal mode order n associated with them, and since each elliptical mode converges toward a unique circular mode in the limit δ → 0, it is meaningful to assign a mode order to the elliptical modes as well.Roughly speaking, the mode order of an elliptical mode is the index n of the dominant term in the sum in Eq. ( 7) (except for the TE-like n 0 mode).Now consider an elliptical mode of some order n.In general, all terms in Eq. ( 7) are needed, but for nearcircular waveguides we have jδj ≪ 1, and hence the first-order approximation F ≈ F 0 δF 1 can be expected to yield accurate results, as each successive higher order term is smaller by a factor of jδj compared to the previous one.To the first order in the perturbation expansion, the transverse profile F of the electric field's z component is found to be (10) in the core, and (11) in the cladding, where we have defined the short-hand notations j l J l pr, k l K l qr, and s l sinl ϕ, and where x 1 to x 6 are the amplitudes of various terms.As mentioned earlier, we can follow the same procedure for the magnetic field's z component and will obtain identical terms for G, the only difference being that all sines are replaced by cosines.Let us denote the magnetic field term amplitudes by y i instead of x i , so y i are to G what x i are to F .Essentially, Eqs. ( 10) and ( 11) are approximations to the transverse profile of E z , and the error in the approximation is of the order of δ 2 .For highly elliptical cores, δ 2 will be comparable to δ, in which case second-and possibly higher order corrective terms δ n F n should be added for improved accuracy.
In the general case, 12 variables are thus needed to characterize both E z and H z .However, depending on the order n, some of the terms in the expressions above might be redundant either because they vanish or because they are of the same functional form as some other term.Again, considering the sine case only, we have x 1 x 3 x 4 x 6 y 3 y 6 0 when n 0. For n 1 we have x 3 x 6 y 3 y 6 0, as these terms are of the same form as the ones with x 1 ; x 4 ; y 1 , and y 4 .For n 2 the sinn − 2ϕ terms become identically zero, which means we must set x 3 x 6 0. For modes with n > 2, all 12 variables are needed.The equations to remove redundant variables for the cosine case can be obtained by interchanging every x i with y i .
To solve the modal problem, we need to satisfy the boundary conditions along an elliptical cross section.For this purpose, we need the transverse field components, E T and H T , which can be written in terms of the longitudinal fields E z Fe iβz and H z Ge iβz .However, unlike the circular case, both fields now depend on all four derivatives ∂ r E z , ∂ ϕ E z , ∂ r H z , and ∂ ϕ H z : where k 2 c p 2 in the core and k 2 c −q 2 in the cladding.Even though these expressions appear lengthy, it should be kept in mind that they involve only trigonometric and Bessel functions and are thus considerably simpler than the Mathieu function treatment [3].It is important to note that the transverse tangential field components are not parallel to the unit vector φ in the QPC system.However, the transverse component tangential to the interface simplifies considerably and can be expressed as and where N ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 2 sin 2 ϕ b 2 cos 2 ϕ p .Note that Eqs. ( 12)-( 15) are the same for the cosine case, as the trigonometric functions in them originate from the coordinate system used and not the fields.
We need to match the tangential components of the electric and magnetic fields (E z ; H z ; E t ; H t ) at the core-cladding interface.
Matching of E z at the interface yields three equations in the general case because of the three different ϕ terms in Eqs.(10) and (11).Similarly, matching of H z , E t , and H t yields nine more equations.The equations can be written in matrix form as Mv 0, where v x 1 ; …; x 6 ; y 1 ; …; y 6 T .The modal propagation constants β are found by requiring detM 0, and the eigenvector of M corresponding to the zero eigenvalue then gives the 12 coefficients x i and y i up to a multiplicative factor.For n < 3, the number of linearly independent equations is smaller.The same procedure can be used to include second-and higher order perturbation terms, but the number of equations grows for higher order terms because of the increasing number of variables x i and y i .
Depending on the effective area of the fundamental mode compared to the core area, even slight ellipticity in the core of an optical fiber can have a significant effect on the birefringence between the slow and fast fundamental modes [3].Manufacturing imperfections thus cause fibers to be randomly birefringent due to core ellipticity and stress-induced anisotropic changes in the refractive index.This is generally detrimental in optical communication systems and a limiting factor for networks based on optical fibers [1].Varying core ellipticity may also lead to linear modecoupling in multimode fibers, which would be problematic for next-generation fiber networks utilizing space-division multiplexing [21].Accurate modeling of random changes in the effective indices and linear coupling coefficients would require the computation of propagation constants and modal field profiles for all possible core ellipticities and orientations, which would be very time consuming without simplifications or approximations.The ETPM is no more complex than previous approximate methods [7][8][9][10][11][12][13][14][15] and has the advantage of accurately determining the modal field profiles as well as the propagation constants of waveguide modes of any order.To compare first-order ETPM with the other methods, Fig. 1 shows the normalized birefringence of the fundamental mode defined as where d and D are the lengths of the semi-minor and -major axes, respectively, and Δβ is the difference between the propagation constants of the slow and fast fundamental modes.The refractive indices of the core and cladding are n 1 1.5 and n 2 1.49, respectively.The fiber core has an axis ratio of d ∕D 0.9, meaning δ −0.19 if we align the major axis along the x direction (a D), or δ ≈ 0.23 if the major axis is along the y axis (b D).The average of the predictions using the two possible orientations for ETPM has been used for Fig. 1.The normalized frequency V for the elliptical fiber is defined as , where k 0 is the vacuum wave number.It should be noted that Eq. ( 8) in the paper by Adams et al. [15] is erroneous, as described by Love et al. [22].We further point out that the equation provided by Love and coauthors also contains a minor error: Eq. (1) in Ref. [22] should have 2Δ 3∕2 instead of 2Δ 3 .The corrected equation has been used for Fig. 1. Figure 1 provides a clear example of the capabilities of the technique proposed in this Letter.The first-order ETPM involves only eight variables and yet is very accurate over the whole range, as Fig. 1 shows.Furthermore, more perturbation orders could be included beyond the first-order term F 1 , which means that ETPM is the only method, aside from the tedious "exact" Mathieu function treatment, capable of yielding arbitrarily accurate results.
In conclusion, we introduced a method, called the ellipse transformation perturbation method (ETPM), to determine the normal modes of oscillating systems with elliptical geometry.The method is mathematically considerably simpler than the conventional Mathieu function treatment and gives quite accurate results for nearly circular geometries, even with just the first-order corrective term, as was demonstrated in the context of birefringence in optical fibers with elliptical cores.Even for a relatively high core ellipticity (axis ratio 0.90), the inclusion of first-order corrections yielded results for fiber birefringence that were comparable to or more accurate than several previous approximate methods.ETPM avoids assumptions, such as weak guidance, and has the additional advantage of being able to provide both the modal eigenvalues and the eigenfunctions (corresponding to propagation constants and modal field distributions in the context of optical waveguides) for any mode to any desired accuracy.