## Abstract

A new technique of numerical analysis of microstructured optical fibers is presented. The technique combines a standard 2D finite difference equations with the discrete function expansion. By doing this one gets a matrix eigenvalue problem of a smaller size and a simple formulation of radiation boundary condition. The new algorithm was tested for the microstructures of different types and excellent agreement of the obtained results with other methods was achieved.

© 2005 Optical Society of America

## 1. Introduction

The explorations of microstructured optical fibers (MOFs) are one of the most exciting recent developments in fibers optics. MOFs usually have a complex transverse structure (large variety of shapes and arrangement of air holes or thin layers) and have attracted a great deal of interest both in theory and experiment [1, 2]. The development of the tools and the analysis methods makes it possible to explore many remarkable properties of these structures. However, due to large complexity of the structures, analytical or semi-analytical solutions can be found for only a few of them [3, 4]. For more complex structures rigorous numerical methods have been proposed including finite element method (FEM) [5, 6] and finite difference (FD) method [7, 8].

Due to the finite number of holes or layers in the cladding, a realistic structure is usually leaky. A leaky character of MOFs is accounted for by an appropriate implementation of radiation boundary condition. While the radiation boundary condition can be formulated analytically, its numerical implementation is difficult. As a result, FEM and FD often make use of approximate radiation conditions.

Among many different techniques of implementation of radiation boundary conditions perfectly matched layer (PML) [9] seems to be one of the most popular approach. The idea of PML is based on the enlarging of the computational domain by some extra cells loaded with anisotropic lossy media of the various conductivity profiles. This approach has been successfully applied in the FEM [5] and the FD method [8, 10]. However, since the performance of PML depends on the choice of parameters, such as the layer thickness and the conductivity profile, some experience is needed to determine their optimal values. The next problem, which one has to cope with using PML in the eigenvalue formuation, is the separation of leaky modes from the Berenger modes [11]. Finally, PML increases the dimension of the considered problem, which not only entails higher memory requirements but also results in longer CPU time.

Another popular technique of implementation of radiation boundary conditions is known as the Bayliss-Gunzburger-Turkel radiation boundary condition (BGT RBC) [12]. It can, however, be applied only in an approximated form [6] and the accuracy of the calculations may be insufficient.

The approach proposed in this work is based on FD with Yee’s mesh and discrete function expansion (DFE) method. Compared to pure FD, DFE gives reduction of the matrix size (the number of variables) of the considered problem and enables straightforward implementation of radiation boundary conditions.

Since the considered microstructures have regions of cylindrical symmetry the implementation of the FD algorithm in cylindrical coordinates is natural. This approach facilitates taking advantage of the symmetry of the structure and reduces the geometrical dimensions of the computational domain, compared to rectangular grid.

## 2. Formulation

#### 2.1. Yee’s mesh

Since the structures analysed in this paper are homogeneous along the *z* direction, we can assume that *E⃗*(*ρ*,*ϕ*,*z*) = *E⃗*(*ρ*,*ϕ*)*e*^{-γz}
and *H⃗*(*ρ*,*ϕ*,*z*) = *H⃗*(*ρ*,*ϕ*)*e*^{-γz}
. Accordingly, Maxwell’s curl equations have the following form:

and

In the presented approach Yee’s mesh in cylindrical coordinate system is used. The computation domain (in general, a segment of a circle) is divided into *M* cells along the *ρ* direction and *N* along *ϕ* direction (the total number of the cells is equal *K* = *M*∙*N*). Details explaining the notation applied for component designation are presented in Fig. 1. Hence, Eqs. (1,2) can be expressed in a discrete matrix form, as:

and

where Δ*ρ* and Δ*ρ* are the space discretization in the *ρ* and *ϕ* direction, respectively and

From Eqs. (3,4), after some algebra, we can easily obtain the following relation:

where ${k}_{0}^{2}$ = *ω*
^{2}
*μ*
_{0}
*ε*
_{0}.

From the discretization of Gauss law ∇⃗∙ (*ε*
_{0}
*εE⃗*) = 0, we get:

Substituting relations (12) and ${T}_{t}^{\left(h\right)}$
${R}_{\mathit{\text{tt}}}^{\left(h\right)}$
${T}_{t}^{\left(e\right)}$
${R}_{\mathit{\text{tt}}}^{\left(\mathrm{e}\right)}$ = *I* to Eq. (11) we obtain an eigenvalue problem in terms of transverse electric field:

where

Solution of this eigenvalue problem provides us with the effective modal index ${n}_{\mathit{eff}}=\frac{-\mathit{j\gamma}}{{k}_{0}}$ and modal fields distributions of the guided modes. The dimension of the matrix *A* is rather large (*dimension* = 2× *the total number of cells* - *number of the cells of imposed values on the boundary*). However, the matrix is very sparse and as a result the solution of eigenvalue problem is fast.

## 2.2. Effective dielectric constants

One of the fundamental problems in most of discrete methods, which one has to cope with, is a curved interface between different media. An improper technique of refractive index averaging can significantly deteriorate accuracy and convergence of the considered problem. Many different techniques of tackling a media interface problem have been proposed [13, 14, 7]. In the current work averaging method based on scheme presented by Kaneda, Houshmand and Itoh [14] is used. The method was originally formulated for the cartesian coordinate system. In the case of cylindrical coordinate system modifications of effective permittivity formulae are required. Following the reasoning proposed in [14], we get (Fig. 2):

and

where functions *ρ* (*ϕ*) and *ϕ*(ρ) describe the interface between the two media types in the cell (dashed line) and *p* is the fraction of the material of dielectric constant *ε*_{b}
in the cross section of a cell. It is easy to check that for the large values of *ρ* relations (16-18) become identical to the ones presented in [14].

## 2.3. Discrete function expansion method

In standard FD algorithms the entire domain is evenly covered with Yee’s mesh. For some homogeneous regions (where low harmonic content of the field distribution can be assumed) discrete representation of the field often implies oversampling. It is well known that the Discrete Fourier Transform (DFT) allows one to represent a continuous function by amplitudes of its harmonics instead of the samples. As shown in [15, 16, 17], it is beneficial to apply this alternative representation to fields in the homogeneous regions.

In order to illustrate the discrete function expansion method, proposed in this section, let us consider a one dimensional scalar function *f*(*x*) *x* ∈ [0,*L*] satisfying the Dirichlet boundary conditions. Discrete representation of that function can be expressed by the set of the function’s samples {*f*_{n}
= *f*(*x*_{n}
)${\}}_{n=0}^{n=N}$ where *x*_{n}
= *n*Δ*x* or (if the function is smooth) by another (smaller) set of the values {*g*_{n}
${\}}_{n=1}^{n=q}$ where

or, in matrix notation, *f* = *Qg* where *Q* is an orthonormal matrix. This kind of transformation is especially advantageous if the function has low harmonic content. In that case a function can be described by *g* vector (of the dimension significantly smaller than the dimension of *f*) without losing an accuracy. Similarly, for the Neumann boundary conditions, cosine function can be used (for other cases both types of the functions must be used). It is possible to apply the method only for the chosen regions of the domain. In that case matrix *Q* consists of sine and/or cosine functions for the chosen regions and the ones on the main diagonal for the others. In such a case, vector *g* consists of the samples of the function *f* and the coefficients of the expansion for the chosen regions.

The idea presented above can by applied to reduce the size of the considered eigenvalue problem in a very simple way. Let us consider an example presented in Fig. 3 - the quarter circle structure loaded with the different media and the boundary consisting of a perfect electric conductor (PEC) and perfect magnetic conductor (PMC). Since any homogeneous region (for all ϕ) can be selected, we can assume that the field in that region can be easily described by sine or cosine functions (and the expansion can be efficient). Let us apply the expansion method to the *l*-th ring (quarter) with the functions sin ϕ, sin3ϕ, …, sin(2*q* + 1)ϕ for *E*
_{ρ} component and cos ϕ, cos 3ϕ,…, cos (2*q* + 1)ϕ for *E*
_{ρ}, component. This choice of the functions corresponds to the imposed boundary conditions (in general sine and cosine functions must be used for both components). In that case, we can introduce *Ẽ*_{t}
of the following form:

$${\rho}_{1}^{\left(e\right)}{E}_{\varphi 1}\phantom{\rule{.2em}{0ex}}\cdots \phantom{\rule{.2em}{0ex}}{\rho}_{\left(l-1\right)\bullet N}^{\left(e\right)}{E}_{\varphi \left(l-1\right)\bullet N}\phantom{\rule{.2em}{0ex}}{\rho}_{\left(l-1\right)\bullet N+1}^{\left(e\right)}{\tilde{E}}_{\varphi 1}\phantom{\rule{.2em}{0ex}}\cdots \phantom{\rule{.2em}{0ex}}{\rho}_{l\bullet N}^{\left(e\right)}{\tilde{E}}_{\mathit{\varphi q}}$$

$${{\rho}_{l\bullet N+1}^{\left(e\right)}{E}_{\varphi \left(l+1\right)\bullet N}\phantom{\rule{.2em}{0ex}}\cdots \phantom{\rule{.2em}{0ex}}{\rho}_{K}^{\left(e\right)}{E}_{\mathit{\varphi K}}]}^{T}$$

where the tilded terms are the coefficients of the expansion and all other terms are the field samples. The relation between *E*_{t}
and *Ẽ*_{t}
has a form

where

and

$${\tilde{Q}}_{\varphi}=\left[\begin{array}{cccccccccc}1& & & & & & & & & \\ & \ddots & & & & & & & & \\ & & 1& & & & & & & \\ & & & \mathrm{cos}\frac{\pi}{2\left(2N-1\right)}& \mathrm{cos}\frac{3\pi}{2\left(2N-1\right)}& \dots & \mathrm{cos}\frac{\left(2q-1\right)\pi}{2\left(2N-1\right)}& & & \\ & & & \vdots & \vdots & & \vdots & & & \\ & & & \mathrm{cos}\frac{\pi \left(2N-2\right)}{2\left(2N-1\right)}& \mathrm{cos}\frac{3\pi \left(2N-2\right)}{2\left(2N-1\right)}& \dots & \mathrm{cos}\frac{\left(2q-1\right)\pi \left(2N-2\right)}{2\left(2N-1\right)}& & & \\ & & & 0& 0& \dots & 0& & & \\ & & & & & & & 1& & \\ & & & & & & & & \ddots & \\ & & & & & & & & & 1\end{array}\right]$$

Since columns of matrix *Q̃* are orthogonal [18], after normalization *Q̃*^{T}*Q̃* = *I*. Then by substituting (21) into Eq. (14), we obtain a new form of the considered eigenvalue problem of a smaller size

where

The expansion is especially efficient for regions which are homogeneous for angular variation in whole domain. In other cases the implementation of DFE can be more complex. Examples of the construction of the *Q* matrix for various problems can by found in [19]. Additionally, as it is reported in [19], harmonic expansion (applied to the region closest to the origin) reduces the norm of the matrix *Ã*.

## 2.4. Radiation condition

According to the Lorentz gauge and the assumption about variation along the *z* direction, the outgoing electric field can be, in general, expressed by

where ${H}_{m}^{\left(2\right)}$ (∙) is a Hankel function, $\kappa =\sqrt{{\gamma}^{2}+{k}_{0}^{2}{n}_{bg}}$ and *n*
_{bg} is the background refractive index.

Let us apply the expansion method presented above for two outermost rings of the cells of the domain. In that case vector *Ẽ*_{t} (in eigenvalue problem (24)) consists of the samples of the field for the inner region (cells) and the coefficients of the expansion (corresponding to the respective angular variation) for the outer region. The expansion involves both components of the field *E*_{ρ}
and *E*_{ϕ}
.

It is easy to see that the coefficients of the expansion correspond to the terms *a*_{m}${H}_{m}^{\left(2\right)}$ (*κρ*) and *b*_{m}${H}_{m}^{\left(2\right)}$(*κρ*) from Eqs. (26,27). For these elements of the vector *Ẽ*_{t}
, we can create the relation between two adjacent rings. If *Ẽ*_{ρm} denotes the *m*-th element of the *E*_{ρ}
component expansion of the *l*-th ring of the inner radius *ρ*_{n}
and *E*_{ρm+q}
of the following ring of the inner radius *ρ*
_{n+N}, we have

Analogically, for the *E*_{ϕ}
component

Relations (28) and (29) represent a boundary condition which can be easily implement in the considered problem (24). However, since *κ* is a function of *γ*, which is unknown, direct solution of the eigenvalue problem is impossible. In the current work an iterative technique is used. At first, we assume an initial value of the γ and then we solve the eigenvalue problem with the (inexact) boundary condition. The solution of the problem provides us with a new value of *γ*, which is used to construct the boundary condition in the next step. The process ends when changes of the value of *γ* are smaller than a selected threshold.

## 3. Numerical results

In order to demonstrate the efficiency of the presented technique three types of the microstructured fibers were analysed (Fig. 4).

The first sample is a photonic crystal fiber with six circular holes arranged in a hexagonal setting shown in Fig. 4(a). The diameter of the holes is *d* = 5*μm* and the pitch length is Λ = 6.75*μm*. The refractive index of the background material is *n*_{bg}
= 1.45 and *n*_{a}
= 1 for holes. The vacuum wavelength used in calculation is *λ* = 1.45*μm*. The radius of the computation domain is *R* = 10*μ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 PEC and/or PMC (depending on a symmetry class of a considered mode [4]) at the structure’s symmetry planes and the radiation condition, presented in section 2.4, at the curved boundary. The discretization in the *ρ* direction is *M* = 100 and *N* = 150 in the *ρ* direction. For 3 outer rings and 30 inner rings (see Fig. 5) *q*= 10 expansion functions are used, which, compared to the pure FD formulation, reduces the number of variables by about 31%.

Table 1 shows the comparison of the results obtained by the presented technique with the results obtained by other methods. Despite a rather coarse mesh the agreement with the other methods is very good. Especially with the Multipole Method (MM) [3, 4] (10 multipole moments) and the Vector FDM-ABC Scheme (70 Azimuthal, 54 radial terms) [20]. An important advantage of the presented technique is short CPU time. The analysis of a single mode requires only a few iteration (usually fewer than 5), each of which lasts approximately 60s on a standard 2.8GHz desktop (e.g., the Vector FDM-ABC scheme requires approximately 6h on a 1GHz desktop [20]). The agreement with FEM (1648 triangular elements) is better for the real part than for the imaginary part of effective index, which is probably caused by the approximate radiation condition used in [6].

In the scheme presented in [8], the FD method for the cartesian grid with PML was used. The results were given for the fundamental mode only *n*_{eff}
= 1.445394 – 3.18 ∙ 10^{-8}
*j*. Assuming that the accurate value is obtained by the Multipole Method [3, 4], we can easily compare the error of the effective index of fundamental mode in the FD scheme [8] with the error of presented approach as a function of the matrix dimension (see Fig. 6). The error is smaller and the convergence is better for the technique proposed in the current work.

To demonstrate the ability of the proposed technique to model microstructures with non-circular shapes, a structure presented in Fig. 4(b) was analysed. 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*. The computational domain is a half-circle with a radius *R* = 2.5*μm*. The discretization in the *ρ* direction was *M* = 100 and *N* = 300 in the *ϕ* direction. The DFE was used for 3 outer rings and 35 inner rings with *q* = 20 functions, which gives the reduction in the number of variables by about 35%.

In Table 2 the comparison of the obtained results with other methods is shown. The results of the FEM analysis were obtained with 1937 triangular elements discretization [6] and for the Vector FDM-ABC scheme, with the resolution of 800 in radial and the 180 in angular direction.

The results for the analysed structure remain in excellent agreement for the non-circular microstructures, especially for the real part of the effective index.

The last example is depicted in Fig. 4(c). The diameter of the holes and the core is *d* = 4*μm* and the pitch length is Λ = 5*μm*. The radius of the computation domain is *R* = 7.5*μm* and the vacuum wavelength is *λ* = 1.5*μm*. The refractive index of the background is *n*_{bg}
= 1.42, core *n*_{c}
= 1.45 and air *n*_{a}
= 1. The discretization in the *ρ* direction was *M* = 130 and in the (*ϕ*) direction *N* = 200. The DFE was used for 3 outer rings and 30 inner rings with *q* = 10 functions (the reduction of the number of variables is about 25%.

The numerical results are presented in Table 3. The imaginary part of the index of the fundamental mode is smaller than machine precision (1E-16) and since ℜ(*n*_{eff}
) > *n*_{bg}
, ℑ(*n*_{eff}
) = 0 can be assumed.

The FD method (cartesian coordinates, 120×120 discretization points) was also used in [7] to calculate the real part of the effective index of the fundamental mode (the analysis did not include the radiation) yielding *n*_{eff}
= 1.435360. Since satisfying agreement between the result of [7] and the current work is achieved, we can expect that the results obtained for the higher modes are valid.

## 4. Conclusions

We have presented a new technique of the size reduction of the standard 2D FDFD eigenvalue problem. To this end the 2D FDFD problem is formulated in cylindrical coordinate system and the discrete function expansion technique for homogeneous rings is applied. This approach enables straightforward implementation of the radiation boundary conditions. The new algorithm was tested for the different types of the microstructures. Advantages of the proposed method were demonstrated (e.g. short CPU time, fast convergence and simple implementation). Excellent agreement of the obtained results with other methods was achieved.

## References and links

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

**
2
. **
C.
Kerbage
and
B. J.
Eggleton
, ”
Numerical analysis and experimental design of tunable birefringence in microstruc-tured optical fiber
,”
Opt. Express
**
10
**
,
246
–
255
(
2002
),
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-5-246
. [PubMed]

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

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

**
5
. **
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
. [PubMed]

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

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

**
8
. **
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] [PubMed]

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

**
10
. **
J.P.
Berenger
, “
Perfectly matched layer for the FDTD solution of wave-structure. Interaction problems
,”
IEEE Trans. Microwave Theory Tech.
**
44
**
,
110
–
117
(
1996
).

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

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

**
13
. **
T.
Weiland
, “
Verlustbehaftete Wellenleiter mit beliebiger. Randkontur und Materialverteilung
,”
AEU
**
33
**
,
170
–
174
(
1979
).

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

**
15
. **
M.
Mrozowski
, “
A Hybrid PEE-FDTD Algorithm for Accelerated Time Domain Analysis of Electromagnetic Waves
,”
IEEE Microwave and Guided Wave Letters
**
4
**
,
323
–
325
(
1994
). [CrossRef]

**
16
. **
M.
Mrozowski
,
M.
Okoniewski
, and
M.A.
Stuchly
, “
A hybrid PEE-FDTD method for efficient field modeling in cyllindrical coordinates
,”
Electronics Letters
**
32
**
,
194
–
195
(
1996
). [CrossRef]

**
17
. **
M.
Wiktor
and
M.
Mrozowski
,“
Efficient Analysis of Waveguide Components Using a Hybrid PEE-FDFD Algorithm
,”
IEEE Microwave and Wireless Components Letters
**
13
**
,
396
–
398
(
2003
). [CrossRef]

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

**
19
. **
M.
Wiktor
and
M.
Mrozowski
,“
Discrete Projection for Finite Difference Methods
,” 20th Annual Review or Progress in Applied Computational Electromagnetics, Syracuse NY
2004
, Conf. proceedings, S04P08.

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