## Abstract

A finite-element-based vectorial optical mode solver is used to analyze microstructured optical waveguides. By employing 1^{st}-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

## 1. Introduction

Since the introduction of the holey fiber [1], various waveguiding structures that utilize the arrangement of microstructured holes [2] or thin layers [3] 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 [6]. The confinement loss discriminates one mode from the others; hence, one can have structures with wide effectively-single-mode operation wavelength range [7] 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 [11], supercell lattice method [12], multipole method [13], Fourier decomposition method [14], beam propagation method (BPM) [15,16], and also mode solver based on finite difference method (FDM) [17], 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 1^{st}-order Bayliss-Gunzburger-Turkel-like (BGT-like) transparent
boundary conditions (TBC) [23] 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 [23]. 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, *H*_{x} and *H*_{y} the *x* and *y* components of
the magnetic field, while ${n}_{\mathit{\text{xx}}}^{2}$ , ${n}_{\mathit{\text{yy}}}^{2}$ , and ${n}_{\mathit{\text{zz}}}^{2}$ 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:

$$-\underset{{\Gamma}_{e}}{\int}\frac{1}{{n}_{\mathit{yy}}^{2}}{w}_{x}\left({\partial}_{x}{H}_{x}+{\partial}_{y}{H}_{y}\right)dy+\underset{{\Gamma}_{e}}{\int}\frac{1}{{n}_{\mathit{xx}}^{2}}{w}_{y}\left({\partial}_{x}{H}_{x}+{\partial}_{y}{H}_{y}\right)dx\}$$

$$+\underset{e}{\sum _{\mathit{InterfaceElement}}}\left\{-\underset{{\Gamma}_{\mathrm{int},e}}{\int}\frac{1}{{n}_{\mathit{yy}}^{2}}{w}_{x}\left({\partial}_{x}{H}_{x}+{\partial}_{y}{H}_{y}\right)dy+\underset{{\Gamma}_{\mathrm{int},e}}{\int}\frac{1}{{n}_{\mathit{xx}}^{2}}{w}_{y}\left({\partial}_{x}{H}_{x}+{\partial}_{y}{H}_{y}\right)dx\right\}$$

$$+\underset{e}{\sum _{\mathit{TriangularElement}}}\underset{{\Omega}_{e}}{\iint}\{\frac{1}{{n}_{\mathit{zz}}^{2}}\left({\partial}_{x}{w}_{y}-{\partial}_{y}{w}_{x}\right)\left({\partial}_{x}{H}_{y}-{\partial}_{y}{H}_{x}\right)+\left[{\partial}_{x}\left(\frac{1}{{n}_{\mathit{yy}}^{2}}{w}_{x}\right)+{\partial}_{y}\left(\frac{1}{{n}_{\mathit{xx}}^{2}}{w}_{y}\right)\right]\left({\partial}_{x}{H}_{x}+{\partial}_{y}{H}_{y}\right)$$

$$+{k}_{0}^{2}{n}_{\mathit{eff}}^{2}\left(\frac{1}{{n}_{\mathit{yy}}^{2}}{w}_{x}{H}_{x}+\frac{1}{{n}_{\mathit{xx}}^{2}}{w}_{y}{H}_{y}\right)-{k}_{0}^{2}\left({w}_{x}{H}_{x}+{w}_{y}{H}_{y}\right)\}dx\phantom{\rule{.2em}{0ex}}dy=0$$

with *w*_{x} and *w*_{y} 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 1^{st}-order BGT-like [24] 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 1^{st}-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

and

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/(2*r*) term) compared to the Sommerfeld-like TBC [25]. 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*. [31].

#### 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
*r*_{b} =10µm as the computational domain with the curved boundary located just slightly after the hole. The 1^{st}-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 [13] and the vector FDM-ABC [14]. 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 ${\mathit{\text{HE}}}_{11}^{a}$ - and ${\mathit{\text{HE}}}_{11}^{b}$ -like, ${\mathit{\text{HE}}}_{21}^{a}$ - and ${\mathit{\text{HE}}}_{21}^{b}$ -like, ${\mathit{\text{EH}}}_{11}^{a}$ - and ${\mathit{\text{EH}}}_{11}^{b}$ -like modes are degenerate pairs, while ${\mathit{\text{HE}}}_{31}^{a}$ - and ${\mathit{\text{HE}}}_{31}^{b}$ -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 [27].

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 [28] 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 [4]

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 ${\mathit{\text{HE}}}_{31}^{a}$ - and ${\mathit{\text{HE}}}_{31}^{b}$ -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 [30], 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 1^{st}-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 [14]. 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 [29] 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 ${\mathit{\text{HE}}}_{11}^{a}$ - and ${\mathit{\text{HE}}}_{11}^{b}$ -like, and also ${\mathit{\text{HE}}}_{21}^{a}$ - and ${\mathit{\text{HE}}}_{21}^{b}$ -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 [28] 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*. [31],
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 [28], 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 [31] 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*. [31]. 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 1^{st}-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.

## Acknowledgments

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

**1. **J.C. Knight, T.A. Birks, P.S.J. Russell, and D.M. Atkin, “All-silica single-mode optical fiber with photonic crystal cladding,” Opt. Lett. **21**, 1547–1549 (1996). [CrossRef] [PubMed]

**2. **P. Russell, “Photonic Crystal Fibers,” Science **299**, 358–362 (2003). [CrossRef] [PubMed]

**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]

**7. **T.A. Birks, J.C. Knight, and P.S.J. Russell, “Endlessly single-mode photonic crystal fiber,” Opt. Lett. **22**, 961–963 (1997). [CrossRef] [PubMed]

**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]

**10. **A. Argyros, N. Issa, I. Basset, and M. van Eijkelenborg, “Microstructured optical fiber for single-polarization air guidance,” Opt. Lett. **29**, 20–22 (2004). [CrossRef] [PubMed]

**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.