A finite-element-based vectorial optical mode solver is used to analyze microstructured optical waveguides. By employing 1st-order Bayliss-Gunzburger-Turkel-like transparent boundary conditions, both the real and imaginary part of the modal indices can be calculated in a relatively small computational domain. Results for waveguides with either circular or non-circular microstructured holes, solid- or air-core will be presented, including the silica-air Bragg fiber recently demonstrated by Vienne et al. (Post-deadline Paper PDP25, OFC 2004). The results of solid-core structures are in good agreement with the results of other methods while the results of air-core structure agree to the experimental results.
© 2004 Optical Society of America
Since the introduction of the holey fiber , various waveguiding structures that utilize the arrangement of microstructured holes  or thin layers  have been realized. Due to the finite number of holes or layers in the cladding, the realistic structures are usually leaky. In this class of structures, both the real and imaginary parts of the modal indices are essential parameters. The real part of the modal indices will determine among others, the dispersion properties [4,5], while the imaginary part will determine the confinement loss . The confinement loss discriminates one mode from the others; hence, one can have structures with wide effectively-single-mode operation wavelength range  or non-polarization-degenerate effectively-single-mode operation [8–10]. Therefore, the ability to calculate both the real and imaginary part of the modal indices is important.
The large variety of possible hole shapes and arrangements demand the use of numerical methods that can handle arbitrary cross-section shape to analyze this kind of structures. Besides, the existence of interfaces with high index-contrast between the host material and air holes requires a vectorial method to accurately model the structure. Among others; plane-wave expansion method , supercell lattice method , multipole method , Fourier decomposition method , beam propagation method (BPM) [15,16], and also mode solver based on finite difference method (FDM) , and finite element method (FEM) [18–20] have been used to model the microstructured waveguide. Some of these methods can calculate only the real part of the mode effective indices. With the multipole method one can calculate both the real and imaginary part of the modal indices, but it is restricted to circular shape of holes. The Fourier decomposition method with adjustable boundary conditions (FDM-ABC) is suitable for calculation of both the real and imaginary part of the modal indices of structure with either circular or non-circular holes arranged in a circularly oriented setting and homogeneous exterior domain. The BPM with perfectly matched layers (PML) is able to calculate both the real and imaginary part of the modal indices of structures with arbitrary cross-section, but might be rather intricate for multimode calculation. FDM and FEM mode solvers employing PML [21,22] can be used to calculate both the real and imaginary part of the modal indices of structure with arbitrary cross-section shape, but the PML itself occupies extra memory space and requires some skills to determine the optimal parameters to work effectively.
In this work, we use a vectorial optical mode solver based on Galerkin finite element method, which is furnished with a 1st-order Bayliss-Gunzburger-Turkel-like (BGT-like) transparent boundary conditions (TBC)  to model various kinds of microstructured optical waveguides. Thanks to the boundary conditions, the structure can be analyzed in a relatively small computational domain for its complex-valued modal indices and field profiles. The structures being considered include those with either solid material or air as the core; circular or non-circular microstructured holes arranged around the core. Although not being demonstrated in this paper; in principle, structures with arbitrary cross-section shapes can be handled, including those with non-circularly oriented hole arrangements and quasi-homogeneous exterior domain. This feature might be useful in the future, if the advancement of the fabrication technology would enable the application of the same principle of waveguide engineering in integrated optics.
2. Description of the method
The details of the scheme, used for mode solving of ordinary waveguides supporting guided and leaky modes has been given elsewhere . For convenience, it will be briefly explained here. Starting from the H-field-based vectorial wave-equation, for longitudinally-invariant structures composed of non-magnetic anisotropic materials with diagonal permittivity tensors and exp(jωt) time dependent of the field; it is possible to get the vectorial wave-equation expressed only in terms of the transverse components of the magnetic field as follows:
Here the x and y denote the transverse Cartesian coordinates associated with the structure cross-section, k 0 the vacuum wavenumber, n eff the complex modal index, Hx and Hy the x and y components of the magnetic field, while , , and are the non-zero entries of the permittivity tensor associated with the x, y, and z components of the electric field. Using the Galerkin procedure and discretizing the computational domain into triangular elements lead to the following discretized weak formulation:
with wx and wy denote the weight functions, Ωe denote the area in each triangular element, Γint,e the line element at the interface between different materials, and Γe the line element at the computational boundaries.
The derivatives of the fields occurring in the boundary term in Eq. (2) will be handled through the 1st-order BGT-like  TBC to reflect the properties of the fields in the exterior domain properly. We use a vector radiation function
along the computational boundary Γ, which leads to a 1st-order operator on the boundary fields as follows:
In Eqs. (3) and (4) r and θ are the polar coordinates of the cross-section whereby the center of the waveguide has been taken as the origin, and k r,x and k r,y are the complex transverse wavenumbers associated with the x and y components of the field. Solving the wave-equation at the elementwise homogeneous anisotropic exterior domain leads to
with Re(k r)>0 for the leaky-mode case (the one used in this paper) and Im(k r)<0 for the guided-mode case. Using Eqs. (4)–(6) and neglecting the angular dependence of the field, the derivative terms within the boundary term of Eq. (2) can be approximated as
where the caret (^) notation denotes the unit vector. Note that the approximation has one extra term (the 1/(2r) term) compared to the Sommerfeld-like TBC . For this reason, it is one order more accurate than the latter one in order to satisfy the radiation function (3).
Expanding the fields in quadratic nodal-based basis functions will lead to a sparse generalized matrix eigenvalue equation, which can be solved using an eigenvalue solver to obtain the eigenvalues related to the modal indices and eigenvectors associated with the transverse components of the magnetic field of the corresponding modes. Since the eigenvalue to be found also appears in the boundary term, a simple iteration scheme has been used to linearize the resulting non-linear eigenvalue equation. In our implementation, we have used Matlab and ARPACK solver encapsulated by the eigs function.
3. Results and discussions
In order to illustrate the ability of the mode solver to model the microstructured waveguides, we choose several types of structures as our samples, including those that employ solid-material and air as the core, with either circular or non-circular shape of microstructured holes. When available, we will compare our results with the results of other methods. We will also show computational results of the air-silica Bragg fiber recently demonstrated experimentally by Vienne et al. .
3.1. Solid-core structure
As the first sample, we choose a photonic crystal fiber with 6 circular holes arranged in a hexagonal setting as shown in Fig. 1. In order to enable comparison of our results with the results of other methods, we choose the same structure parameters as in the references [13,14]. The diameter of the holes is d=5µm with pitch length of Λ=6.75µm. The refractive index of the background material is n bg=1.45, while the refractive index of the holes is n hole=1. The vacuum wavelength used in the calculation is 1.45µm. Taking advantage of the symmetry of the structure; we only use a quarter circle with a radius of rb =10µm as the computational domain with the curved boundary located just slightly after the hole. The 1st-order BGT-like TBC is applied to the curved boundary while symmetry boundary conditions which consist of a perfect electric conductor (PEC) and/or a perfect magnetic conductor (PMC) are applied at the two boundaries coinciding with the structure symmetry planes. The computational domain is discretized into 1648 triangular elements as shown in Fig. 1.
Table 1 shows our computed results for the first-ten modes and their comparison to the results of other methods. The mode labeling in the table makes use of the similarity of the mode profiles to those of ordinary step-index fiber with additional superscripts a and b denoting the results obtained using PMC and PEC at the horizontal symmetry plane, respectively. The table shows that our results using rather modest mesh size and small computational domain are in good agreement with the results of the multipole method  and the vector FDM-ABC . The agreement in the real part of the effective indices is better than that of the imaginary part. By varying the mesh size and the computational boundary position, we observe larger absolute changes in Re(n eff) than Im(n eff), indicating that the absolute error in the latter is a few order of magnitude better than that in the former. This is attributed in the following: the error in Re(n eff) is due to the contributions from the interior and the boundary of the domain, while the error in the Im(n eff) is dominated by the contribution from the boundary. The boundary terms themselves occupy only a small fraction of the non-zero entries within the FEM matrices and have relatively small values. However, although their absolute errors are small, due to their very small values, it is very hard to get high relative accuracy for the Im(n eff). Since only modest mesh size and a relatively small computational window were used in the results shown in Table 1, the convergence test indicates that in general only the order of magnitude (for some modes also the first non-zero digit) of the computed Im(n eff) is significant, while for the Re(n eff), at least 4 decimal digits for the high order modes and 5 decimal digits for the low order modes are significant. Of course, better accuracy can be expected by using a larger computational window and finer mesh, but with the price of more expensive computational effort.
The transverse fields of modes labeled with the same subscript but different superscript (and also TE- and TM-like modes) are perpendicular to each other at every point in the structure cross-section (see, e.g., Fig. 5 below), hence they can be regarded as pairs. From the computed results (see Table 1, and also later on Figs. 2 and 3) and group theoretical arguments [13,26], we can see that - and -like, - and -like, - and -like modes are degenerate pairs, while - and -like, and also TE 01 - and TM 01 -like modes are non-degenerate pairs. The table also shows that using different symmetry boundary conditions at the two symmetry planes, the computed effective indices of the degenerate modes agree to each other up to 5 decimal digits in Re(n eff). The small numerical birefringence comes from the fact that the mesh used in the discretization has broken the structure symmetry; the effect of which will show up in the discrepancies of the discretization error of the computed degenerate modes. By reducing the discretization error, e.g. using finer mesh, one can expect smaller numerical birefringence .
We also take a similar structure and perform a spectral-scan characterization. We assume that pure silica has been used as the host material and take the more realistic value of the refractive index from its Sellmeier’s equation  for each wavelength. In this way, the material dispersion effect is rigorously taken into account. For the holes, we take a constant value of refractive index of 1. Figure 2 shows the real part of the effective indices and the corresponding dispersion parameter as function of wavelength. The dispersion parameter is calculated using 
The curves show three groups of modes, which correspond to the first-three LP modes of scalar analysis. By using this vectorial solver, it is possible to distinguish modes within the same group, which is more pronounced for longer wavelength region, where the dimension of the structure becomes more comparable to the wavelength. The divergence between curves for - and -like modes also confirms their non-degeneracy as discussed earlier. The figure also shows that the HE 11 -like mode has zero-dispersion wavelength at shorter wavelength than the ordinary fiber. Figure 3(a) shows the confinement loss as function of wavelength. The difference between modes within the same groups is even more pronounced here. Since this loss is one of the important properties of the microstructured waveguides, these results emphasize that a vectorial method is necessary for such structure, especially in the long wavelength regime (and also at large d/Λ ratio , which determines the local air filling fraction in the ring of holes). As the wavelength increases, the modes become less quasi-confined; consequently, the confinement loss also increases. As the HE 11 -like mode is the one with lowest loss, the structure will be effectively single-moded after a certain propagation distance. Figure 3(b) shows the effect of adding more rings of holes around the central core, with the holes in the cladding arranged in a triangular-lattice-like setting. In this case; 2-ring, 3-ring, and 4-ring structures will have 18, 36, and 60 holes in the cladding, respectively. Adding rings of holes will reduce the confinement loss, but will not change the group velocity dispersion (which related to the dispersion parameter) very much. Some small fluctuations in the curve of the structure with 60 holes for short wavelength region are due to the very small value of the loss, which is beyond the machine accuracy of the computing platform.
In order to demonstrate the ability to model waveguides with non-circular microstructured holes, we take a structure with 3 annular-shaped holes as shown in Fig. 4. The host material is pure silica with a refractive index of 1.44402362 at a wavelength of 1.55µm, while the refractive index of the holes is 1. The annular-shaped holes have an inner radius r 1=1µm and outer radius r 2=2µm and angular width of 108°. We take a half-circle with radius r b=2.5µm as our computational domain with the 1st-order BGT-like TBC at the curved boundary, and either PEC or PMC at the boundary located at the structure symmetry plane. The computational domain is discretized into 1937 triangular elements as shown in Fig. 4. The results of the first-six modes of this structure are given in Table 2 with comparison to the results obtained using the vector FDM-ABC . The vector FDM-ABC used in this calculation is the one that uses Fourier decomposition only in the azimuthal direction, while in the radial direction it uses a finite difference discretization. Due to the vectorial character of the modes, scalar method  is simply not enough, hence it will not be considered in the comparison. The table shows good agreement between results of our FEM scheme and the vector FDM-ABC. As in previous example, the - and -like, and also - and -like modes are degenerate pairs, while TE 01 - and TM 01 -like modes are non-degenerate pairs. The vector plots of the transverse magnetic fields of these modes are given in Fig. 5. The real part of the effective indices and the dispersion parameter of this structure are depicted in Fig. 6, while the imaginary part of the effective indices and the confinement loss are shown in Fig. 7. Again, to obtain these spectral plots, we used the Sellmeier’s equation  of pure silica at each wavelength. Figure 6(a) shows that at around 1.483µm, the real part of n eff of HE 21 - and TM 01 -like modes cross over, which indicates their rather dissimilar dispersion properties as shown in Fig. 6(b). Since this structure has much smaller core size and larger (local) air filling fraction than the previous sample, the effect of the air holes is stronger, hence the zero-dispersion wavelength of the HE 11 -like mode is located at a shorter wavelength. Besides, the vectorial character is also more pronounced as indicated by more divergent curves of TE 01 -, TM 01 -, and HE 21 -like modes, both in their real and imaginary part of the effective indices.
3.2. Air-core structure
To demonstrate the ability to model an air-core structure, we pick up a structure recently demonstrated experimentally by Vienne et al. , which is an air-core structure with 3 rings of annular-shaped holes around it. We use the model depicted in Fig. 8, which parameters are taken from the fabricated structure but slightly simplified for computational efficiency. The radius of the air core is r core=10µm. We assume that the holes are of annular-shaped with uniform thickness of t annular=2.3µm. The corners of the holes are rounded with circles tangent to the side of the annular-shaped holes with radius of r corner=t annular/4. The thin rings of the host material located in-between two neighboring rings of holes have thickness of t ring=0.2µm. These thin rings of solid material are supported by thin bridges with a mid position thickness of t bridge=45nm. For computational efficiency, we assume that the number of holes are 24, 34, and 44 in the first, second, and third rings of holes, respectively, which are slightly different than the fabricated one that has 24, 35, and 46 holes at the corresponding rings. We use a computational window of a quarter circle with a radius of r b=19µm, with the curved boundary located not too far from the last ring of holes. The 1st-order BGT-like TBC is applied on this curved boundary, while the PEC and/or PMC are used at the boundaries coinciding with the structure symmetry planes. The computational domain is discretized into 17638 triangular elements. This large number of elements is induced by the very thin bridges located between two adjacent holes in the same ring. In the calculation, we used the wavelength of 1.06µm, which corresponds to the laser used in the experiments to measure the modal near field intensity pattern. The refractive index of the silica as the host material is taken from the Sellmeier’s equation , while the refractive index of air is 1.
The results of the first-six modes are given in Table 3, which shows that TE 01 -like mode exhibits lowest confinement loss; a property that used to be found in the Bragg fibers [8,9]. This is as expected by the designer(s) of this structure, since the structure will be a Bragg fiber if we simply neglect the existence of the thin bridges. The value of the calculated confinement loss of the dominant TE 01-like mode agrees very well to the measured results in the experiment  which is 0.015 dB/cm at λ=1.04µm. It should be noted that since our model has only 2-fold rotational symmetry, all the modes are non-degenerate. Figure 9 shows that the modal fields are very well quasi-confined into the air-core due to the quasi-bandgap effect of the finite alternating solid-air rings in the cladding. Figure 10 shows the mode profiles of the HE 11 - and TE 01-like modes, while their longitudinal component of time averaged Poynting vectors are given in Fig. 11. The profile of the HE 11 -like mode shows the existence of bright resonance spots at the cladding layers that were also observed in the nearfield patterns measured in the experiments of Vienne et al. . To the best of our knowledge, these computational results are the first rigorous results for the realistic model of this novel structure.
4. Concluding remarks
A finite-element-based vectorial optical mode solver furnished with a 1st-order BGT-like TBC has been applied to analyze the microstructured optical waveguides. The boundary conditions allow the calculation of both the real and imaginary part of the modal indices in a relatively small computational domain. Both the solid- and air-core structures with either circular or non-circular holes have been used as the samples. The computed results agree to the results obtained using other methods or experiments.
The authors acknowledge the help of N.A. Issa from Univ. Sydney by providing the results of the vector FDM-ABC for structure with 3 annular-shaped holes for comparison. Discussions with E. van Groesen, M. Hammer, and R. Stoffer are also acknowledged. This work is supported by STW Technology Foundation.
References and links
3. Y. Fink, D.J. Ripin, S. Fan, C. Chen, J.D. Joannopoulos, and E.L. Thomas, “Guiding optical light in air using an all-dielectric structure,” J. Lightwave Technol. 17, 2039–2041 (1999). [CrossRef]
4. A. Ferrando, E. Silvestre, J.J. Miret, and P. Andres, “Nearly zero ultraflattened dispersion in photonic crystal fibers,” Opt. Lett. 25, 790–792 (2000). [CrossRef]
5. J.C. Knight, J. Arriaga, T.A. Birks, A. Ortigosa-Blanch, W.J. Wadsworth, and P.S.J. Russell, “Anomalous dispersion in photonic crystal fiber,” Photonics Technol. Lett. 12, 807–809 (2000). [CrossRef]
6. T.P. White, R.C. McPhedran, C.M. de Sterke, L.C. Botten, and M.J. Steel, “Confinement losses in microstructured optical fibers,” Opt. Lett. 26, 1660–1662 (2001). [CrossRef]
8. S.G. Johnson, M. Ibanescu, M. Skorobogatiy, O. Weisberg, T.D. Engeness, M. Soljacic, S.A. Jacobs, J.D. Joannopoulos, and Y. Fink, “Low-loss asymptotically single-mode propagation in large-core OmniGuide fibers,” Opt. Express 9, 748–779 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-13-748. [CrossRef] [PubMed]
9. I.M. Basset and A. Argyros, “Elimination of polarization degeneracy in round waveguides,” Opt. Express 10, 1342–1346 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-23-1342. [CrossRef]
11. A. Ferrando, E. Solvestre, J.J. Miret, P. Andres, and M.V. Andres, “Full-vector analysis of a realistic photonic crystal fiber,” Opt. Lett. 24, 276–278 (1999). [CrossRef]
12. W. Zhi, R. Guobin, L. Shuqin, and J. Shuisheng, “Supercell lattice method for photonic crystal fibers,” Opt. Express 11, 980–991 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-9-980. [CrossRef] [PubMed]
13. T.P. White, B.T. Kuhlmey, R.C. McPhedran, D. Maystre, R. Renversez, C.M. de Sterke, and L.C. Botten, “Multipole method for microstructured optical fibers. I. Formulation,” J. Opt. Soc. Am B 19, 2322–2330 (2002). [CrossRef]
14. N.A. Issa and L. Poladian, “Vector wave expansion method for leaky modes of microstructured optical fibers,” J. Lightwave Technol. 21, 1005–1012 (2003). [CrossRef]
15. F. Fogli, L. Saccomandi, P. Bassi, G. Bellanca, and S. Trillo, “Full vectorial BPM modeling of index-guiding photonic crystal fibers and couplers,” Opt. Express 10, 54–59 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-1-54. [CrossRef] [PubMed]
16. K. Saitoh and M. Koshiba, “Full-vectorial imaginary-distance beam propagation method based on a finite element scheme: application to photonic crystal fibers,” J. Quantum Electron. 38, 927–933 (2002). [CrossRef]
17. C.P. Yu and H.C. Chang, “Applications of the finite difference mode solution method to photonic crystal structures,” Opt. Quantum Electron. 36, 145–163 (2004). [CrossRef]
18. F. Brechet, J. Marcou, D. Pagnoux, and P. Roy, “Complete analysis of the characteristics of propagation into photonic crystal fibers by finite element method,” Optical Fiber Technol. 6, 181–191 (2000). [CrossRef]
19. Cucinotta, S. Selleri, L. Vincetti, and M. Zoboli, “Holey fiber analysis through the finite element method,” Photonics Tech. Lett. 14, 1530–1532 (2002). [CrossRef]
20. M. Koshiba, “Full-vector analysis of photonic crystal fibers using the finite element method,” IEICE Trans. Electron. E85-C, 881–888 (2002).
21. D. Ferrarini, L. Vincetti, M. Zoboli, A. Cucinotta, and S. Selleri, “Leakage properties of photonic crystal fibers,” Opt. Express 10, 1314–1319 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-23-1314. [CrossRef] [PubMed]
22. K. Saitoh and M. Koshiba, “Leakage loss and group velocity dispersion in air-core photonic bandgap fibers,” Opt. Express 11, 3100–3109 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-23-3100 [CrossRef] [PubMed]
23. H.P. Uranus, H.J.W.M. Hoekstra, and E. van Groesen, “Galerkin finite element scheme with Bayliss-Gunzburger-Turkel-like boundary conditions for vectorial optical mode solver,” J. Nonlinear Optical Phys. Materials (to be published).
24. A. Bayliss, M. Gunzburger, and E. Turkel, “Boundary conditions for the numerical solution of elliptic equations in exterior regions,” SIAM J. Appl. Math. 42, 430–451 (1982). [CrossRef]
25. H.E. Hernandez-Figueroa, F.A. Fernandez, Y. Lu, and J.B. Davies, “Vectorial finite element modeling of 2D leaky waveguides,” Trans. Magnetics 31, 1710–1713 (1995). [CrossRef]
26. R. Guobin, W. Zhi, L. Shuqin, and J. Shuisheng, “Mode classification and degeneracy in photonic crystal fibers,” Opt. Express 11, 1310–1321 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-11-1310. [CrossRef] [PubMed]
27. M. Koshiba and K. Saitoh, “Numerical verification of degeneracy in hexagonal photonic crystal fibers,” Photonics Technol. Lett. 13, 1313–1315 (2001). [CrossRef]
28. I.H. Malitson, “Interspecimen comparison of the refractive index of fused silica,” J. Opt. Soc. Am. 55, 1205–1209 (1965). [CrossRef]
29. L. Poladian, N.A. Issa, and T.M. Monro, “Fourier decomposition algorithm for leaky modes of fibres with arbitrary geometry,” Opt. Express 10, 449–454 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-10-449. [CrossRef] [PubMed]
30. T.M. Monro, D.J. Richardson, N.G.R. Broderick, and P.J. Bennett, “Modeling large air fraction holey optical fibers,” J. Lightwave Technol. 18, 50–56 (2000). [CrossRef]
31. G. Vienne, Y. Xu, C. Jakobsen, H.J. Deyerl, T.P. Hansen, B.H. Larsen, J.B. Jensen, T. Sorensen, M. Terrel, Y. Huang, R. Lee, N.A. Mortensen, J. Broeng, H. Simonsen, A. Bjarklev, and A. Yariv, “First demonstration of air-silica Bragg fiber,” Post Deadline Paper PDP25, Optical Fiber Conference 2004, Los Angeles, 22–27 Feb. 2004.