Fourier decomposition algorithm for leaky modes of fibres with arbitrary geometry

(b)The right to post and update his or her Work on any internet site (other than the Author(s’) personal web home page) provided that the following conditions are met: (i) access to the server does not depend on payment for access, subscription or membership fees; and (ii) any such posting made or updated after acceptance of the Work for publication includes and prominently displays the correct bibliographic data and an OSA copyright notice (e.g. "© 2009 The Optical Society").


Introduction
Microstructured or holey optical fibres have attracted much interest recently because they can exhibit properties (dispersion, mode size) very different from conventional fibres.These structures can consist exclusively of air holes in a uniform host material.The central region can confine light if it is surrounded by air holes, or alternatively, it can be effectively air-suspended by fine supporting filaments (Fig. 1).Although the majority of microstructured fibres have been fabricated by stacking capillaries, and so the first approach has been predominantly used to date, the second approach has recently been demonstrated to be practical using alternative fabrication techniques such as extrusion [1].Regardless of the fabrication technique, single-material fibres do not have a raised index core, and so cannot support true bound modes.Instead, leaky modes describe the propagation and account for the gradual leakage of light between and across the holes.
Most existing algorithms rely on basis function expansions [2,3] and accurately yield many fibre properties (spot-size, dispersion) for low leakage.However, they cannot determine confinement loss since they assume either periodic or zero boundary conditions, thus ignoring the outward radiating leakage field.One formulation [4] matches the solution to the external field while only using an azimuthal expansion but restricts attention to bound modes.Scalar beam propagation is numerically intensive but has nevertheless been successfully used [5] to model confinement loss in microstructures where the wavelength was smaller than the hole spacing.Attempts to use vector beam propagation [6] in structures with elliptical holes and longer wavelengths has been less successful in finding confinement losses.Recently, a multipole method has been used [7] which correctly models the outward radiating leakage field but is limited to circular holes, and so is of restricted usefulness for the structures like those shown in Fig. 1.The algorithm presented here has the simplicity of a basis function expansion technique but can also find confinement loss by correctly determining the field outside the computational domain satisfying the outward radiating boundary condition.The confinement loss in dB per unit length is related to the imaginary part of the complex effective mode index by −20 log 10 (e)kIm[n eff ].This paper explores the new algorithm in the context of the scalar wave equation, since it was easier to investigate technical issues such as the best choice of basis functions, convergence and accuracy with the simpler system first.Generalisation to the vector wave equation will be presented in a subsequent paper.

Structure of algorithm
A circular computational domain D of radius R is used to encapsulate the waveguide structure, and the index is uniform and equal to n cl in the infinite region outside the domain.The scalar wave equation in dimensionless coordinates x = r/R is where The field outside the computational domain D can be expressed exactly in the form where K m (W x) are modified Bessel functions.For bound modes W is real and these fields are evanescent.For radiating and leaky modes W is complex and the fields express an outward propagating field.In either case, the expansion outside D has the correct physical behaviour at infinity.The field inside D is expanded in some complete set of basis functions A set of weighting functions φ mn (x, θ) is chosen [8], and the following generalised eigenvalue problem for the unknown coefficients A mn and the eigenvalue W 2 is obtained: where the matrix elements are given by The choice of inner product and basis functions is discussed in Appendix A. The complex eigenvalues of leaky modes arise from the non-self-adjoint boundary conditions applied to the system, which forces the matrix operators of Eq. ( 4) to be non-Hermitian.Thus, techniques which exploit Hermitian operator properties and variational theorems are unavailable to us, and as discussed below, we use inverse iteration to find the eigenvalues of interest.The essential properties of the expansion functions ψ mn are that they: have the correct behaviour at the origin (such that the coordinate singularity at the origin is not problematic); contain adjustable parameters that allow continuity with the expansion outside the domain (allowing adjustable boundary conditions); and allow computationally efficient evaluation of the matrix elements.Our choice is a polar-coordinate harmonic Fourier decomposition( Eq. ( 6) ).The basis functions can be chosen to exploit rotational symmetries of the structure for computational advantage.
The adjustable boundary condition Fourier decomposition method (ABC-FDM) proceeds as follows.To start the procedure, an initial guess for the value of n eff is used to determine the parameter W appearing in the external field expansion.The adjustable parameters in the basis functions (Eq.( 7)) are chosen to match the external field.The matrix problem is solved for the eigenvalue and eigenvector of interest.Inverse iteration methods can be used to narrow in on the eigenvalue closest to some chosen value, and additional eigenvalues and eigenvectors can subsequently be found by working in the subspace orthogonal to the modes already found.The eigenvalue then gives an improved estimate of n eff for the mode of interest.If required, this new estimate can be used to re-adjust the external field and the adjustable boundary conditions again to obtain an improved estimate and continued iteration can be used to increase the accuracy of the answers.Such re-iteration will converge quickly for any mode with the majority of its guided power within the computational domain.Since the imaginary part is often orders of magnitude smaller than the real part our experience was that, in most cases of interest, a single reiteration was sufficient to accurately determine the imaginary part of the index unless the loss is extremely large.
Although we apply our algorithm to structures which only have leaky modes, it can be used directly with structures that have both bound and leaky modes.The algorithm proceeds exactly as above except that the adjustable parameter W remains real for all bound modes.However, if one is only interested in bound modes, formulations that exploit self-adjointness of bound mode problems will in general be more efficient.

Numerical results
As our first example we calculate the first two modes of the three fold symmetric structure shown in Fig. 1(a).The modes are calculated at λ = 1.55µm and the inner and outer radii of the annular sector shaped holes are 1 and 2 µm respectively with an angular width of 108 degrees.Throughout we assume the material to be undoped silica, although this is not a restriction of the method.The contours are in decibels, and the fields are also displayed out beyond the domain radius.The phase of the fields is denoted by the colour.These modes were found using 1050 basis functions (50 radial terms, 21 azimuthal terms) which took of the order of 6 mins on a desktop Pentium III.
Note that in Fig. 2(a) the central confined region of the field of the first mode has essentially a constant phase (i.e. a planar phase front similar to a bound mode), but that beyond the confining holes the phase front becomes an outward travelling wave.The effective index is 1.374+0.000048iyielding a confinement loss of 1.7 dB/mm.The higher order mode (b) shows the expected angular variation in phase of the form exp(iφ). Its effective index and loss are 1.255 + 0.00075i and 27 dB/mm respectively.
The second structure we analyse is the hollow core structure shown in Fig. 1(d).The inner air hole has radius 1.05 µm.The first ring of elliptical holes have radii 0.3 and 0.9 µm at distance 1.65 µm from the origin.The second ring of circular holes have radii 0.3 µm and are 2.55 µm from the centre.These modes were calculated using 100 radial and 41 azimuthal terms taking a few hours.Two different modes are shown: one which is guided in the material around the inner hole Fig. 2(c), and one which is guided within the hole Fig. 2(d).The mode in (c) has effective index and loss of 1.243+0.0075iand 270 dB/mm respectively.The mode in (d) has effective index and loss of 0.848+0.00095iand 35 dB/mm respectively.All of these modes are extremely lossy and would be unlikely to be observed in practice without the addition of many more holes to improve the confinement.Despite this, the algorithm is still able to find the leaky modes, and can describe modes that are guided either by average index or band-gap effects.We next show how loss varies with wavelength, the number of holes and hole dimensions.loss for the structures shown in Fig. 1(a) and (c).Notice that increasing the outer radius to 3 µm, as in structure (c) substantially reduces the confinement loss as expected.Note these values are comparable (accounting for material and scattering loss and the scalar approximation) with the measurement [1] of 4 dB/m for a fibre with the same geometry as Fig. 1(c) with f = 0.9 but a higher glass index of 1.8 at 1.55 µm.Fig. 3(c) shows the real part of the effective index of the fundamental mode.Note that it depends weakly on f suggesting such properties can be approximated without including the bridges.

Inverse Error Analysis
If we substitute the calculated solution back into the original wave equation we can calculate a residual and use this as an error estimate.Moreover, if we rearrange the wave equation and solve for the index profile then we can use the calculated field to reconstruct the index profile.The departure of this reconstructed profile from the original profile reveals the error in the calculation (within the scalar approximation).Fig. 4 shows the reconstructed waveguides for the simplest and most complicated examples from Fig. 1.A contour plot is also shown for the second reconstruction to aid in visualisation.All the images reveal the tell-tale signs of a Fourier decomposition.The presence of small amplitude radial rings is related to the number of radial Fourier components used in the calculation.The effect of truncating the Fourier expansions can also be seen in the ripples in the angular direction as well as Gibbs phenomena associated with discontinuous profiles.As the number of basis functions is increased the reconstructed waveguide will converge to the original waveguide.The small size of these features give a good indication of the convergence of the solution and don't affect the mode solution appreciably.In particular, the calculated field is the exact mode of this reconstructed guide.Treating the reconstructed guide as a perturbation should allow the use of perturbation theorems [9] to estimate the error in the effective index.

Conclusion
The ABC-FDM can be used to find leaky modes of a large variety of structures.We intend to explore more efficient algorithms for finding eigenvectors of non-self-adjoint systems in future work.Detailed studies of convergence, efficiency and the generalisation to the vector case will also be forthcoming.

Fig. 1 .
Fig. 1.Four different fibre structures investigated with the ABC-FDM.The structures are all to scale and the core radius in structure (a) is 1 micron.

Fig. 2 .
Fig. 2. (a)Fundamental and (b) first higher order mode of 3 hole structure shown in Fig 1(a).(c)Ring mode and (d) hollow core mode of the structure shown in Fig 1(d).The contours are equally spaced at 3 dB intervals.The color indicates phase.The radial phase variation reveals outward propagating waves associated with confinement loss.

Fig. 3 .
Fig. 3. (a) Comparison of leakage loss for 3 and 6 hole structures shown in Fig. 1(a) and (b) as a function of wavelength and air fraction f .(b) Similar comparison of leakage loss for 3 holes with outer radii r = 2.0 and 3.0 µm.(c) Effective index.

Fig. 4 .
Fig. 4. (a) Reconstruction of waveguide from computed solution of Fig. 2(a).(b) Reconstruction of hollow core waveguide from solution shown in Fig. 2(d).(c) Contour map of the same reconstruction as in (b).Dimensions are in microns.