Open Access
Add to BibSonomyAbstract
A new algorithm for calculating the confinement loss of leaky modes in arbitrary fibre structures is presented within the scalar wave approximation. The algorithm uses a polar-coordinate Fourier decomposition method with adjustable boundary conditions (ABC-FDM) to model the outward radiating fields. Leaky modes are calculated for different examples of microstructured fibres with various shaped holes.
©2002 Optical Society of America
1 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 - 20log10(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.
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.
2 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 ncl in the infinite region outside the domain. The scalar wave equation in dimensionless coordinates x = r/R is
where V 2(x,θ) = k 2 R 2[n 2(x,θ) - ] depends on the index profile n(x,θ) of the fibre and W 2 = k 2 R 2[ -] is a constant that depends on the effective index of the mode n eff. Note that V(x,θ) ≡ 0 outside D.
The field outside the computational domain D can be expressed exactly in the form
where Km (Wx) 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 Amn 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.
3 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.000048i yielding 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.0075i and 270 dB/mm respectively. The mode in (d) has effective index and loss of 0.848+0.00095i and 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.
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.
The solid and dashed curves in Fig. 3(a) compare the leakage loss for the lowest order mode as a function of wavelength for the structures shown in Fig. 1(a) and (b) respectively. The loss is shown for three different air fractions: f=0.8, 0.9 and 1.0. Notice that for the complete annulus (f=1.0) the 3 and 6 bridge structures become identical. Observe that the leakage loss produced by the bridges becomes relatively more important for shorter wavelengths, and that 6 thin bridges produces less loss than 3 thick bridges for the same air fraction. Likewise, the solid and dashed curves in Fig. 3(b) compare the 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.
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.
4 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.
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.
5 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.
Acknowledgments
We acknowledge the support of the Australian Research Council, Sydney VisLab and our colleagues Ian Bassett, Steven Manos and Whayne Padden.
A Choice of Basis Functions
We us a generalised Galerkin technique [8] wish expansion functions ψmn and weight functions ϕmn . The presence of leakage destroys the self-adjointness of M regardless of the choice of basis functions or inner product. Thus, we have the freedom to choose functions purely on the basis of computational efficiency. We investigated a large variety of different basis function sets and finally chose
where
The adjustable parameters α and β provide for the correct behaviour at the origin and also satisfy the boundary condition ψ′(x, 0)/ψ(x, 0) = WK′m(Wx)/Km (W) at x = 1. Our treatment differs from that in [4] because we expand both the azimuthal and radial dependence in complete Fourier series whereas in [4] the radial dependence is treated using Runge-Kutta techniques. Both treatments allow matching to arbitrary boundary conditions.
The azimuthal indices m vary over all integers and the radial indices n vary over all positive integers. Appropriate subsets of the azimuthal indices can be used if the modes have rotational symmetries. The second set of functions is complete over the space of bounded functions within the computational domain using the standard inner product
and the first set can be shown to be complete by using standard tests.
The integrals for 〈ϕj |∇2 ψi 〉 and 〈ϕj |ψi 〉 do not depend on the specific structure and involve simple trigonometric functions and so can be pre-tabulated for computational efficiency. The integrals for 〈ϕj |V 2(x,θ)ψi 〉 depend on the structure and can also be done analytically if the holes are in the shape of annular sectors, otherwise they can be done numerically. Most importantly though, because the trig functions appearing are all harmonically related, the number of distinct sum and difference frequencies that arise when products of trig functions appear in the inner products is only of order 2n max where n max is the largest radial index used in the basis set. This number is much less than the number of elements in the matrices and makes the scheme computationally viable even if the integrals need to be calculated numerically.
References and links
1. T.M. Monro, K.M. Kiang, J.H. Lee, K. Frampton, Z. Yusoff, R. Moore, J. Tucknott, D.W. Hewak, H.N. Rutt, and D.J. Richardson, “High nonlinearity extruded single-mode holey optical fibers” in Optical Fiber Communications Conference 2002 PD-FA1.
2. A. Ferrando, E. Silvestre, J. J. Miret, P. Andrés, and M. V. Andrés, “Full-vector analysis of a realistic photonic crystal fiber,” Opt. Lett. 24, 276–278 (1999). [CrossRef]
3. T. M. Monro, D. J. Richardson, N. G. R. Broderick, and P.J. Bennet, “Holey optical fibers: An efficient modal method,” J. Lightwave Technol. 17, 1093–1102 (1999). [CrossRef]
4. Y-H. Wang and C. Vassallo, “Circular fourier analysis of arbitrarily shaped optical fibers,” Opt. Lett. 14, 1377–1379 (1989). [CrossRef] [PubMed]
5. B. J. Eggleton, P. S. Westbrook, C. A. White, C. Kerbage, R. S. Windeler, and G. L. Burdge, “Cladding-mode-resonances in air-silica microstructure optical fibers,” J. Lightwave Technol. 18, 1084–1100 (2000). [CrossRef]
6. M. J. Steel and R. M. Osgood Jr., “Polarization and dispersive properties of elliptical-hole photonic crystal fibers,” J. Lightwave Technol. 19, 495–503 (2001). [CrossRef]
7. T. P. White, R. C. McPhedran, C. M. de Sterke, L. C. Botten, and M. J. Steel,“Confinement losses in microstructured optical fibres,” Opt. Lett. 26, 1660–1662 (2001). [CrossRef]
8. C. A. J. Fletcher, Computational Galerkin Methods, (Springer Verlag,1984). [CrossRef]
9. A. W. Snyder and J. D. Love, Optical Waveguide Theory, (Chapman and Hall,1983) p. 376.
