## Abstract

A highly accurate radiation boundary condition for finite difference analysis of open waveguides is introduced. The boundary condition is applicable to the structures embedded in a homogeneous medium and fitted to the cross section of the structure. The numerical tests carried out for a few types of waveguides including microstructured fibers showed that the proposed approach improves the accuracy by about an order of magnitude in comparison with the PML technique and eliminates all its disadvantages.

© 2007 Optical Society of America

## 1. Introduction

The finite difference (FD) technique is one of the most versatile methods of computational electromagnetics and as such it is often used to investigate propagation characteristics of open waveguides [1, 2, 3]. In order for the method to deal with the radiation effects two main categories of techniques can be distinguished. The first one involves transparent boudary conditions (TBC) which are based on an approximation of the Sommerfeld radiation condition (Bayliss-Turkel, Mur, Higdon schemes [4]). The second category (proposed by Berenger [5]) is based on the termination of a computational domain with the region (filled with a lossy anisotropic medium) that absorbs outgoing waves. Since then many variants of an artificial medium, known as a perfectly matched layer (PML), have been proposed [4]. Nowadays, the PML is the method of choice for the finite difference analysis of open problems. The popularity of PML comes from the fact that it gives more accurate results. Numerical tests show that the reflection coefficient on the boundary (obtained from finite difference time domain simulations [4]) is about -40*dB* for TBC and about -120*dB* for PML. However, TBC are still applied in other techniques (beam propagation method [6] or finite element method [7]).

Despite its immense popularity, PML has several obvious disadvantages when it comes to the analysis of guidance properties of radiating structures. Firstly, PML enlarges the computational domain, secondly the results are sensitive to the selection of parameters such as conductivity profile, permittivity permeability, thickness of the layer and a distance to the structure. Finally, artificial Berenger modes appear in the numerical solution and their identification can be very difficult [8, 9].

In this paper we propose a new radiation boundary condition that eliminates all of the above disadvantages. The new approach involves only points located on the outer structures boundary and therefore its application leads to smaller matrix eigenvalue problem than PML technique. The algorithm is based on the idea introduced in our earlier paper [3], but it is applicable to Yee’s mesh in rectangular coordinates and relies on rigorous rather than approximate relationships for the outgoing waves.

## 2. Formulation

In order to present the main idea of the proposed algorithm one-dimensional structures are considered first. For most of them analytical solutions are known and can be used for the first verification of the accuracy of the new technique.

#### 2.1. One dimensional propagation problems

Let us consider a planar slab dielectric waveguide, parallel to *yz* plane, characterized by permittivity *ε*(*x*). Assuming that the fields variation along the *z* direction is represented by a term e^{-γz} and there is no variation in the *y* direction, Maxwell’s equations can be separated into two independent systems for *TE ^{x}* and

*TM*modes.

^{x}Let us consider the equation for *TE ^{x}* modes

According to the assumptions of FD method [4], differential equations (1) are approximated by

where *e _{m}* =

*E*(

_{y}*m*Δ

*x*), ${h}_{m}={H}_{z}\left(\left(m+\frac{1}{2}\right)\Delta x\right)$,

*ε*=

_{m}*ε*(

*m*Δ

*x*) and Δ

*x*denotes the discretization step (see fig. 1).

In a matrix notation, from the above relations we get an eigenvalue problem for a propagation constant:

where *R _{zy}*

^{(e)},

*R*

_{yz}^{(H)}and

*ε*are square

*M*×

*M*matrices

and *e* = [*e*
_{1},*e*
_{2},⋯,*e _{M}*]

^{T},

*h*= [

*h*

_{1},

*h*

_{2},⋯,

*h*]

_{M}^{T}.

#### 2.2. Boundary conditions for ID analysis

Without losing generality, we can assume that a three layered structure (see fig. 1) represents any planar slab dielectric waveguide (*ε*
_{2} = *ε*
_{2}(*x*)).

The field shape in region 1 and 3 is known

$${E}_{y}\left(x\right)={A}_{+}{e}^{-{\kappa}_{3}x},\phantom{\rule{.5em}{0ex}}x>\frac{b}{2}$$

where *κ _{i}* = [-

*γ*

^{2}- k

^{2}

_{0}

*ε*]

_{i}^{1/2},

*k*

_{0}=

*ω*(

*μ*

_{0}

*ε*

_{0})

^{1/2}and

*A*

_{-},

*A*

_{+}are unknown coefficients. For evanescent modes

*κ*, is real and grater than zero, however for leaky modes

_{i}*κ*, is complex and the proper sheet of the Riemann surface must be chosen to satisfy the Sommerfeld radiation condition [10].

_{i}Expressions (7) determine the relations between the field samples at the boundary of the computational domain

where *a*
_{-} = e^{-κ1Δx} and *a*
_{+} = e^{-κ3Δx}.

Let us rewrite problem (3) in a simple form

where *A* =*R _{yz}*

^{(h)}

*R*

_{zy}^{(e)}-

*ωμ*

_{0}

*ε*, then

Applying relations (8) to equations system (10) we get a new eigen problem

However, since *a*
_{-} and *a*
_{+} are the functions of *γ* above problem becomes nonlinear

The most straightforward technique of solving nonlinear problems is a simple iteration method. In the first step we have to solve a regular eigenvalue problem for the initial value *γ*
_{0} and denote the solution as *γ*
_{1}:

In the next step we have to solve the linear problem (13) for *γ*
_{2}, substituting *γ*
_{1} into operator *Ã*. The procedure should be repeated until the difference between two sequential values of *γ* is less than the assumed error.

In general, there is no proof of convergence of that algorithm, so it may happens that the process starts to diverge. If this occurs the convergence can be restored by considering a function *F*(*γ*) =*γ*-*$\widehat{\gamma}$*, where *$\widehat{\gamma}$* is an eigenvalue of the problem *Ã*(*γ*)*ẽ* = *$\widehat{\gamma}$*
^{2}
*ẽ*. A zero of this function is simultaneously the solution of eigenproblem (12). Since complex plane of *γ* can be replaced by two dimensional real space and function *F*(*γ*) by log|*F*(*γ*)| solving of (12) transforms to a problem of finding the minimum of the function.

For this purpose we have selected 3-point simplex method, but any other technique of finding minimum or zero can be applied. Since minimalization is time consuming it is used only to restore the convergence of (13).

#### 2.3. Boundary conditions for 2D analysis

In practice one has to use numerical analysis for more complex waveguides. For this purpose let us consider a cross section of open waveguide in a computational domain covered with rectangular Yee’s mesh. The finite difference method gives a matrix eigenvalue problem in terms of transverse electric field [1, 2]:

The boundary condition is derived based on an analytical representation of outgoing waves. In a homogeneous region, assuming variation along the *z* direction in a form e^{-γz} the electric and magnetic field can be expressed by only two components *E _{z}*(

*ρ*,

*φ*,

*z*) and

*H*(

_{z}*ρ*,

*φ*,

*z*). The other components are unambiguously defined by the relations:

According to the Sommerfeld radiation condition, any outgoing electromagnetic field, can be expressed by the series

where *H _{m}*

^{(2)}(∙) is a Hankel function of the second kind,

*κ*= [

*γ*

^{2}+

*k*

^{2}

_{0}

*με*]

^{1/2},

*η*= [

*μ*/

*ε*)

^{1/2}and

*A*,

_{m}*B*,

_{m}*C*,

_{m}*D*are arbitrary coefficients. In practice we have to reduce the series to a finite number of terms, so let us assume that

_{m}*m*= 0,1, ...,

*Q*.

Because the mesh is rectangular, the fields should be expressed in the Cartesian coordinate system. Only *E _{x}*(

*x*,

*y*,

*z*) and

*E*(

_{y}*x*,

*y*,

*z*) are needed, hence

and

where *x* = *ρ* cos(*φ*) and *y* = *ρ* sin(*φ*).

Substituting the above relations and taking *z* = 0, we get

$$\phantom{\rule{.2em}{0ex}}=\sum _{m=0}^{Q}\{{A}_{m}\frac{\gamma}{\kappa}\phantom{\rule{.2em}{0ex}}\left[\frac{m}{\rho \kappa}{H}_{m}^{\left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(\phi \right)-{H}_{m}^{\prime \left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(\phi \right)\right]$$

$$\phantom{\rule{.2em}{0ex}}-{B}_{m}\frac{\gamma}{\kappa}\phantom{\rule{.2em}{0ex}}\left[\frac{m}{\rho \kappa}{H}_{m}^{\left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(\phi \right)+{H}_{m}^{\prime \left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(\phi \right)\right]$$

$$\phantom{\rule{.2em}{0ex}}+{C}_{m}\frac{\omega \mu}{\kappa \eta}\phantom{\rule{.2em}{0ex}}\left[\frac{m}{\rho \kappa}{H}_{m}^{\left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(\phi \right)+{H}_{m}^{\prime \left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(\phi \right)\right]$$

$$\phantom{\rule{.2em}{0ex}}-{D}_{m}\frac{\omega \mu}{\kappa \eta}\phantom{\rule{.2em}{0ex}}\left[\frac{m}{\rho \kappa}{H}_{m}^{\left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(\phi \right)+{H}_{m}^{\prime \left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(\phi \right)\right]\}$$

and

$$\phantom{\rule{.2em}{0ex}}=\sum _{m=0}^{Q}\{{A}_{m}\frac{\gamma}{\kappa}\phantom{\rule{.2em}{0ex}}[-\frac{m}{\rho \kappa}{H}_{m}^{\left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(\phi \right)-{H}_{m}^{\prime \left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(\phi \right)]$$

$$\phantom{\rule{.2em}{0ex}}-{B}_{m}\frac{\gamma}{\kappa}\phantom{\rule{.2em}{0ex}}[-\frac{m}{\rho \kappa}{H}_{m}^{\left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(\phi \right)+{H}_{m}^{\prime \left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(\phi \right)]$$

$$\phantom{\rule{.2em}{0ex}}+{C}_{m}\frac{\omega \mu}{\kappa \eta}\phantom{\rule{.2em}{0ex}}\left[\frac{m}{\rho \kappa}{H}_{m}^{\left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(\phi \right)-{H}_{m}^{\prime \left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(\phi \right)\right]$$

$$\phantom{\rule{.2em}{0ex}}-{D}_{m}\frac{\omega \mu}{\kappa \eta}\phantom{\rule{.2em}{0ex}}\left[\frac{m}{\rho \kappa}{H}_{m}^{\left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{sin}\left(\phi \right)+{H}_{m}^{\prime \left(2\right)}\left(\kappa \rho \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(m\phi \right)\phantom{\rule{.2em}{0ex}}\mathrm{cos}\left(\phi \right)\right]\}.$$

Using the above series one can find the relation between the field samples outside and inside the computational domain.

Let *S ^{C}* be a set of points of Yee’s mesh (where tangential electric field is defined) inside the computational domain close to the boundary (white circles in fig. 2) and

*S*be a set of points on the boundary (crosses in fig. 2). Let us denote by

^{B}*E*and

^{C}_{t}*E*the field values from set

^{B}_{t}*S*and

^{C}*S*, respectively. From relations (23) and (24) we get

_{B}and

where vector *C* = [*A*
_{0},... ,*A _{Q}*,

*B*

_{0},... ,

*B*,

_{Q}*C*

_{0},... ,

*C*,

_{Q}*D*

_{0},...,

*D*]

_{Q}^{T}. The relation between field values on the boundary and inside the computational domain is

where *M _{C}*

^{(inv)}=

*M*

_{C}^{-1}if the number of the

*S*elements is equal to the size of

^{C}*C*vector, otherwise to invert

*M*the singular value decomposition (SVD) algorithm must be applied [11] (

_{C}*M*

_{C}^{(inv)}

*M*=

_{C}*I*).

A relation (27) can be used to modify operator *A* in eigenvalue problem (14). The new problem involves only fields located inside the structure - all points outside (and on) the boundary can be neglected.

Since the boundary condition matrices *M _{B}* and

*M*

_{C}^{(inv)}are functions of

*γ*, we get a nonlinear eigenvalue problem similar to (12). Again simple iteration is in most cases enough to reach the solution, but if this fails the convergence is restored by applying minimalization procedure described in the previous paragraph.

## 3. Numerical results

For numerical tests we have selected a few types of structures, for which analytical or quasi analytical solutions are known. To emphasize the differences between the proposed approach and standard techniques all of the structures were also analysed using an anisotropic PML with 10 layers (parameters set to the optimal values:
$\sigma =\frac{0.8\left(m+1\right)}{120\pi \Delta n}$ where *m* = 4, *n* is a refractive background index and Δ is a step of the discretization [2, 4]). To refractive index averaging, in cells with an interface between different media, the method presented by Kaneda, Houshmand and Itoh [12] was used.

The effective indices

for the first three modes of a 1D structure with parameters *ε*
_{r1} = *ε*
_{r3} = 1.21, *ε*
_{r2} = 1, *b* = 1*nm*, *λ*
_{0} = 0.2*nm* are presented in Table 1 (corresponding fields distributions are shown in fig. 3 – 5).

Also in structures with larger contrast good agreement was achieved. The results for *ε*
_{r1} = *ε*
_{r3} = 9, *ε*
_{r2} = 1, *b* = 1*nm*, *λ*
_{0} = 1.5*nm* are collected in Table 2.

Since PML is designed for absorbing plane waves, the results obtained for 1D are in good agreement with the theoretical values. However, the convergence and accuracy of PML algorithm is worse than for our formulation especially for high contrast structure. Note that for coarse mesh the percentage error in the value of imaginary part of *n _{eff}* exceeds 10% and is at least two orders of manitude greater than for the new technique.

It must be emphasized that contrary to the PML the proposed technique does not enlarge the computational domain. The results remain unchanged even when the distance between the boundary and the structure is reduced to two cells (for guided as well as for leaky modes).

Moreover, the numerical solution obtained in the simulation gives the possibility of the field evaluation in any point outside the computational domain: (7) in 1D and (19–24) in 2D.

To investigate the performance of the new boundary condition for 2D structures we shall consider a simple optical fiber [9]. The radius of the fiber is *R* = 0.5*μm*, the core index *n _{c}* = 2.9, the background index

*n*= 1. 55 and

_{bg}*λ*

_{0}= 1

*μm*. The discretization assumed in PML simulation was 200 × 200 cells (80,000 of variables) and the domain size was 2

*μm*× 2

*μm*. The PML region was placed 40 cells from the core. This was necessary since a smaller distance to PML results in a significant increase of the simulation error. The alternative analysis was carried out with the same discretization and the proposed boundary condition (

*Q*= 20). The number of variables was only about 20,000, because the boundary condition was imposed very close to the structure

*R*= 0.55

_{B}*μm*(see fig. 6). All numerical results are collected and compared with the theoretical values in Table 3. A significant inaccuracy of the imaginary part of the effective index obtained from PML technique is shown, especially for guided modes.

The numerical tests were also carried out for two different types of the microstructures shown in fig. 7. The first example (fig. 7a) is a photonic crystal fiber with six circular holes arranged in a hexagonal setting. The radius of the holes is *r* = 2.5*μm* and the pitch length is Λ = 6.75*μm*. The refractive index of the background material is *n _{bg}* = 1.45 and

*n*= 1 for holes. The vacuum wavelength used in calculation is

_{a}*λ*= 1.45

*μm*. The radius of the computation domain is

*R*= 9.5

*μm*, but due to the symmetry of the structure only a quarter of the circle has to be analysed. The boundary conditions consist of a perfect electric and/or magnetic conductor at the structure symmetry planes and the radiation condition presented in section 2 at the curved boundary. The discretization assumed in simulation is 150 × 150 cells. The PML approach requires 45,000 variables (the domain is a square 11

*μm*× 11

*μm*). The new boundary condition allows one to fit a boundary to the contour, which reduces the number of variables by 33% to 30,000.

The comparison of the results obtained by the presented technique with the results of other methods is shown in Table 4 and 5. The agreement with the other methods is very good. The results for quasi analytical Multipole Method (MM) [13] are used as a reference.

Since MM provides very accurate values, it is seen that the new boundary condition is more accurate in comparison to PML (an order of magnitude accuracy improvement for the imaginary part).

The last example is a microstructured fiber with non-circular shapes, presented in fig. 7b. The angular-shaped holes have the inner radius *r*
_{1} = 1*μm*, the outer radius *r*
_{2} = 2*μm* and the angular width of 108°. The refractive index of the background material is *n _{bg}* = 1.44402362. The vacuum wavelength is

*λ*= 1.55

*μm*. This time, the computational domain is a half-circle with a radius

*R*= 2.25

*μm*. The discretization was 300×150 cells and the size of the operator was reduced to about 60,000 from 90,000 used in simulation with PML technique (with a rectangular domain 5

*μm*× 2.5

*μm*).

Since the Multipole Method can not handle the fiber under investigation table 6 and 7 compare the obtained results with the Vector FDM-ABC scheme [14] (with the resolution: 800 in radial and the 180 in angular direction).

The results obtained from the new method remain in good agreement with the reference data both for the real and imaginary part. Note that for PML technique the error of the imaginary part of the effective index is extremely high (54.4%) for the *HE*
_{21} - like mode.

The proposed approach leads to the nonlinear eigenvalue problem so it may appear that the method is much more time consuming than PML technique. However, a significant reduction of the problem size makes them comparable. For example, a single eigen solution for the optical fiber lasts about 120*s* for PML and 30*s* for the proposed algorithm (on a standard 2.8GHz desktop). For the microstructured fibers (only partial use of boundary conditions) the difference was smaller (70*s* versus 40*s* and 140*s* vs. 90*s*).

All the results (except a few modes of optical fibers) were obtained in simple iterations with a fast convergence. For example, the convergence process of the fundamental mode of the first microstructure is presented in Table 8

Finally, the presented approach is by far much more effective than some other numerical methods offering high accuracy results (e.g. the Vector FDM-ABC scheme requires approximately 6h on a 1GHz desktop [14]).

It has to be underlined that in proposed technique each mode is analysed separately. However, in practice also in standard PML approach situation is similar. The huge number of artificial (Berenger) modes obtained in simulations with PML causes that solvers are often not able to find more than one proper solution at a single run. Another important problem of PML is a separation of the artificial modes (such modes do not appear in proposed algorithm). The problem occurs even for a homogeneous circular fiber, especially for lower modes [9]. The part of the spectrum obtained in the PML analysis for the first of the microstructures is depicted in fig. 8 (Berenger and leaky modes have similar values and may become indistinguishable in some cases).

## 4. Conclusions

We have presented a new radiation boundary condition for one- and two-dimensional FD eigenvalue problem. The new algorithm was tested for different types of the structures and high accuracy of the results was achieved for all cases. Most of the disadvantages of PML techniques were eliminated.

## Acknowledgment

This work is supported by the Polish Ministry of Science and Higher Education under contract 3 T11F 037 30. We also thank M. Liskin for discussions and providing the PML code for this study.

## References and links

**1. **Z. Zhu and T.G. Brown, “Full-vectorial finite-difference analysis of microstructured optical fibers,” Opt. Express **10**, 853–864 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853. [PubMed]

**2. **Shangping Guo, Feng Wu, Sacharia Albin, Hsiang Tai, and Robert S. Rogowski,“Loss and dispersion analysis of microstructured fibers by finite-difference method,” Opt. Express **12**, 3341–3352 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-15-3341. [CrossRef]

**3. **P. Kowalczyk, M. Wiktor, and M. Mrozowski, “Efficient finite difference analysis of microstructured optical fibers,” Opt. Express **13**, 10349–10359 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-25-10349. [CrossRef]

**4. **A. Taflove and S.C. Hagness, “Computational electrodynamics: the finite-difference time-domain method,” Artech House, Boston (2005), 3rd edn.

**5. **J.P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]

**6. **G.R. Hadley, “Transparent boundary condition for the beam propagation method,” IEEE J. Quantum Electron., **28**, 363–370 (1992). [CrossRef]

**7. **H.P. Uranus and H.J.W.M. Hoekstra, “Modeling of microstructured waveguides using a finite-element-based vectorial mode solver with transparent boundary conditions,” Opt. Express **12**, 2795–2809 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-12-2795. [CrossRef]

**8. **H. Rogier and D. De Zutter, “Berenger and Leaky Modes in Microstrip Substrates Terminated by a Perfectly Matched Layer,” IEEE Trans. Microwave Theory Tech. **49**, 712–715 (2001). [CrossRef]

**9. **H. Rogier and D. De Zutter, “Berenger and Leaky Modes in Optical Fibers Terminated with a Perfectly Matched Layer,” J. Lightw. Technol., **20**, 1141 – 1148 (2002). [CrossRef]

**10. **E.M. Kartchevski, A.I. Nosich, and G.W. Hanson, “Mathematical Analysis of the Generalized Natural Modes of an Inhomogeneous Optical Fiber,” J. Appl. Math. **65**, 2033 – 2048 (2005).

**11. **C.D. Meyer, “Matrix analysis and applied linear algebra”, SIAM, Philadelphia (2000).

**12. **N. Kaneda, B. Houshmand, and T. Itoh, “FDTD analysis of dielectric resonators with curved surfaces,” IEEE Trans. Microwave Theory Tech. **45**, 1645–1649 (1997). [CrossRef]

**13. **T.P. White, B.T. Kuhlmey, R.C. McPhedran, D. Maystre, R. Ranversez, C.M. de Sterke, L.C. Botten, and M.J. Steel, “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]