Efficient finite difference analysis of microstructured optical fibers

doi:10.1364/OPEX.13.010349

« Show journal navigation

Efficient finite difference analysis of microstructured optical fibers

P. Kowalczyk, M. Wiktor, and M. Mrozowski »View Author Affiliations

Optics Express, Vol. 13, Issue 25, pp. 10349-10359 (2005)
http://dx.doi.org/10.1364/OPEX.13.010349

View Full Text Article

Acrobat PDF (120 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Full Text
Article Info
References (20)
Cited By
Figures (6)
Metrics

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.

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

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]

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

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]

]. For more complex structures rigorous numerical methods have been proposed including finite element method (FEM) [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]

] and finite difference (FD) method [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]

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

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

] 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

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

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

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]

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

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]

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

[\begin{matrix} ρ^{- 1} & 0 & 0 \\ 0 & ρ & 0 \\ 0 & 0 & ρ^{- 1} \end{matrix}] [\begin{matrix} 0 & γ & \frac{\partial}{\partial ϕ} \\ - γ & 0 & - \frac{\partial}{\partial ρ} \\ - \frac{\partial}{\partial ϕ} & \frac{\partial}{\partial ρ} & 0 \end{matrix}] [\begin{matrix} E_{ρ} \\ ρ E_{ϕ} \\ E_{z} \end{matrix}] = - jω μ_{0} [\begin{matrix} H_{ρ} \\ ρ H_{ϕ} \\ H_{z} \end{matrix}]

(1)

and

[\begin{matrix} ρ^{- 1} & 0 & 0 \\ 0 & ρ & 0 \\ 0 & 0 & ρ^{- 1} \end{matrix}] [\begin{matrix} 0 & γ & \frac{\partial}{\partial ϕ} \\ - γ & 0 & - \frac{\partial}{\partial ρ} \\ - \frac{\partial}{\partial ϕ} & \frac{\partial}{\partial ρ} & 0 \end{matrix}] [\begin{matrix} H_{ρ} \\ ρ H_{ϕ} \\ H_{z} \end{matrix}] = jω ε_{0} [\begin{matrix} ε_{ρ} & 0 & 0 \\ 0 & ε_{ϕ} & 0 \\ 0 & 0 & ε_{z} \end{matrix}] [\begin{matrix} E_{ρ} \\ ρ E_{ϕ} \\ E_{z} \end{matrix}]

(2)

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:

Fig. 1. The discrete field components in the n-th cell (n = 1..K).

[\begin{matrix} T_{t}^{(e)} & 0 \\ 0 & T_{z}^{(e)} \end{matrix}] [\begin{matrix} γ R_{tt}^{(e)} & R_{tz}^{(e)} \\ R_{zt}^{(e)} & 0 \end{matrix}] [\begin{matrix} E_{t} \\ E_{z} \end{matrix}] = - jω μ_{0} [\begin{matrix} H_{t} \\ H_{z} \end{matrix}]

(3)

and

[\begin{matrix} T_{t}^{(h)} & 0 \\ 0 & T_{z}^{(h)} \end{matrix}] [\begin{matrix} γ R_{tt}^{(h)} & R_{tz}^{(h)} \\ R_{zt}^{(h)} & 0 \end{matrix}] [\begin{matrix} H_{t} \\ H_{z} \end{matrix}] = jω ε_{0} [\begin{matrix} P_{t} & 0 \\ 0 & P_{z} \end{matrix}] [\begin{matrix} E_{t} \\ E_{z} \end{matrix}]

(4)

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

T_{t}^{(e)} = diag (ρ_{1}^{{(h)}^{- 1}}, \dots, ρ_{K}^{{(h)}^{- 1}}, ρ_{1}^{(e)}, \dots, ρ_{K}^{(e)}), T_{z}^{(e)} = diag (ρ_{1}^{{(e)}^{- 1}}, \dots, ρ_{K}^{{(e)}^{- 1}})

(5)

T_{t}^{(h)} = diag (ρ_{1}^{{(e)}^{- 1}}, \dots, ρ_{K}^{{(e)}^{- 1}}, ρ_{1}^{(h)}, \dots, ρ_{K}^{(h)}), T_{z}^{(h)} = diag (ρ_{1}^{{(h)}^{- 1}}, \dots, ρ_{K}^{{(h)}^{- 1}})

(6)

P_{t} = diag (ε_{ρ 1}, \dots, ε_{ρK}, ε_{ϕ 1}, \dots, ε_{ϕK}), P_{z} = diag (ε_{z 1}, \dots, ε_{zK})

(7)

R_{tt}^{(e)} = [\begin{matrix} 0 & I \\ - I & 0 \end{matrix}] = R_{tt}^{(h)}, R_{tz}^{(e)} = [\begin{matrix} F \\ - R \end{matrix}] = R_{zt}^{{(h)}^{T}}, R_{zt}^{(e)} = [- F R] = R_{tz}^{{(h)}^{T}}

(8)

R = Δ ρ^{- 1} [\begin{matrix} - 1 & 1 \\ ⋱ & ⋱ \\ ⋱ & 1 \\ ⋱ \\ - 1 \end{matrix}], F = Δ ϕ^{- 1} [\begin{matrix} - 1 & 1 \\ ⋱ & ⋱ \\ ⋱ & ⋱ \\ ⋱ & 1 \\ - 1 \end{matrix}]

(9)

E_{t} = [\begin{matrix} E_{ρ 1} \\ ⋮ \\ E_{ρK} \\ ρ_{1}^{(e)} E_{ϕ 1} \\ ⋮ \\ ρ_{K}^{(e)} E_{ϕK} \end{matrix}], E_{z} = [\begin{matrix} E_{z 1} \\ ⋮ \\ E_{zK} \end{matrix}], H_{t} = [\begin{matrix} H_{ρ 1} \\ ⋮ \\ H_{ρK} \\ ρ_{1}^{(h)} H_{ϕ 1} \\ ⋮ \\ ρ_{K}^{(h)} H_{ϕK} \end{matrix}], H_{z} = [\begin{matrix} H_{z 1} \\ ⋮ \\ H_{zK} \end{matrix}]

(10)

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

γ^{2} T_{t}^{(h)} R_{tt}^{(h)} T_{t}^{(e)} R_{tt}^{(e)} E_{t} + γ T_{t}^{(h)} R_{tt}^{(h)} T_{t}^{(e)} R_{tz}^{(e)} E_{z} + T_{t}^{(h)} R_{tz}^{(h)} T_{z}^{(e)} R_{zt}^{(e)} E_{t} = k_{0}^{2} P_{t} E_{t}

(11)

where

k_{0}^{2}

= ω ² μ ₀ ε ₀.

From the discretization of Gauss law ∇⃗∙ (ε ₀ εE⃗) = 0, we get:

γ P_{z} E_{z} = D_{t} E_{t}

(12)

D_{t} = T_{z}^{(e)} ∙ [- R^{T} - F^{T}] ∙ T_{t}^{{(e)}^{- 1}}

(13)

Substituting relations (12) and

T_{t}^{(h)}

R_{tt}^{(h)}

T_{t}^{(e)}

R_{tt}^{(e)}

= I to Eq. (11) we obtain an eigenvalue problem in terms of transverse electric field:

A E_{t} = γ^{2} E_{t}

(14)

where

A = T_{t}^{(h)} R_{tt}^{(h)} T_{t}^{(e)} R_{tz}^{(e)} P_{z}^{- 1} D_{t} P_{t} + T_{t}^{(h)} R_{tz}^{(h)} T_{z}^{(e)} R_{zt}^{(e)} - k_{0}^{2} P_{t}

(15)

Solution of this eigenvalue problem provides us with the effective modal index

n_{eff} = \frac{- jγ}{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.

Fig. 2. The media interface

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

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]

, 7

]. In the current work averaging method based on scheme presented by Kaneda, Houshmand and Itoh [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]

] 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

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]

], we get (Fig. 2):

ε_{ρ} = Δ ϕ^{- 1} \ln \frac{ρ_{1}}{ρ_{2}} {[\int_{ρ_{1}}^{ρ_{2}} \frac{dρ}{ρ [(ϕ (ρ) - ϕ_{1}) ε_{a} + (ϕ_{2} - ϕ (ρ)) ε_{b}]}]}^{- 1}

(16)

ε_{ϕ} = Δϕ {[\ln \frac{ρ_{2}}{ρ_{1}} \int_{ϕ_{1}}^{ϕ_{2}} \frac{dϕ}{ε_{a} \ln \frac{ρ (ϕ)}{ρ_{1}} + ε_{b} \ln \frac{ρ_{2}}{ρ (ϕ)}}]}^{- 1}

(17)

and

ε_{z} = p ε_{b} + (1 - p) ε_{a}

(18)

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

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]

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

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]

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

{\begin{matrix} f_{1} = g_{1} \sin (\frac{π}{N}) + g_{2} \sin (\frac{2 π}{N}) + \dots + g_{q} \sin (\frac{qπ}{N}) \\ f_{2} = g_{1} \sin (\frac{π 2}{N}) + g_{2} \sin (\frac{2 π 2}{N}) + \dots + g_{q} \sin (\frac{qπ 2}{N}) \\ ⋮ \\ f_{N - 1} = g_{1} \sin (\frac{π (N - 1)}{N}) + g_{2} \sin (\frac{2 π (N - 1)}{N}) + \dots + g_{q} \sin (\frac{qπ (N - 1)}{N}) \end{matrix}

(19)

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.

Fig. 3. Example of the structure used to illustrate the scheme of the DFE method implementation (different colors represents different media).

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(2q + 1)ϕ for E _ρ component and cos ϕ, cos 3ϕ,…, cos (2q + 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:

{\tilde{E}}_{t} = [E_{ρ 1} \dots E_{ρ (l - 1) ∙ N} {\tilde{E}}_{ρ 1} \dots {\tilde{E}}_{ρq} E_{ρ (l + 1) ∙ N} \dots E_{ρK}

ρ_{1}^{(e)} E_{ϕ 1} \dots ρ_{(l - 1) ∙ N}^{(e)} E_{ϕ (l - 1) ∙ N} ρ_{(l - 1) ∙ N + 1}^{(e)} {\tilde{E}}_{ϕ 1} \dots ρ_{l ∙ N}^{(e)} {\tilde{E}}_{ϕq}

{ρ_{l ∙ N + 1}^{(e)} E_{ϕ (l + 1) ∙ N} \dots ρ_{K}^{(e)} E_{ϕK}]}^{T}

(20)

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

E_{t} = \tilde{Q} {\tilde{E}}_{t}

(21)

where

\tilde{Q} = [\begin{matrix} {\tilde{Q}}_{ρ} \\ {\tilde{Q}}_{ϕ} \end{matrix}]

(22)

and

{\tilde{Q}}_{ρ} = [\begin{matrix} 1 \\ ⋱ \\ 1 \\ 0 & 0 & \dots & 0 \\ \sin \frac{π}{2 N - 1} & \sin \frac{3 π}{2 N - 1} & \dots & \sin \frac{(2 q - 1) π}{2 N - 1} \\ ⋮ & ⋮ & ⋮ \\ \sin \frac{π (N - 2)}{2 N - 1} & \sin \frac{3 π (N - 2)}{2 N - 1} & \dots & \sin \frac{(2 q - 1) π (N - 2)}{2 N - 1} \\ 1 \\ ⋱ \\ 1 \end{matrix}]

{\tilde{Q}}_{ϕ} = [\begin{matrix} 1 \\ ⋱ \\ 1 \\ \cos \frac{π}{2 (2 N - 1)} & \cos \frac{3 π}{2 (2 N - 1)} & \dots & \cos \frac{(2 q - 1) π}{2 (2 N - 1)} \\ ⋮ & ⋮ & ⋮ \\ \cos \frac{π (2 N - 2)}{2 (2 N - 1)} & \cos \frac{3 π (2 N - 2)}{2 (2 N - 1)} & \dots & \cos \frac{(2 q - 1) π (2 N - 2)}{2 (2 N - 1)} \\ 0 & 0 & \dots & 0 \\ 1 \\ ⋱ \\ 1 \end{matrix}]

(23)

Since columns of matrix Q̃ are orthogonal [18

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

], after normalization Q̃^TQ̃ = I. Then by substituting (21) into Eq. (14), we obtain a new form of the considered eigenvalue problem of a smaller size

\tilde{A} {\tilde{E}}_{t} = γ^{2} {\tilde{E}}_{t}

(24)

where

\tilde{A} = {\tilde{Q}}^{T} A \tilde{Q}

(25)

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

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.

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

E_{ρ} = \sum_{m = 1}^{\infty} a_{m} H_{m}^{(2)} (κ ρ) \sin (m ϕ) e^{- γ z}

(26)

E_{ϕ} = \sum_{m = 0}^{\infty} b_{m} H_{m}^{(2)} (κ ρ) \cos (m ϕ) e^{- γ z}

(27)

where

H_{m}^{(2)}

(∙) is a Hankel function,

κ = \sqrt{γ^{2} + k_{0}^{2} n_{b g}}

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}^{(2)}

(κρ) and b_m

H_{m}^{(2)}

(κρ) 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

{\tilde{E}}_{ρ m} = \frac{H_{m}^{(2)} (κ ρ_{n}^{(e)})}{H_{m}^{(2)} (κ ρ_{n + N}^{(e)})} {\tilde{E}}_{ρ m + q}

(28)

Analogically, for the E_ϕ component

{\tilde{E}}_{ϕ m} = \frac{ρ_{n}^{(2)} H_{m}^{(2)} (κ ρ_{n}^{(h)})}{ρ_{n + N}^{(h)} H_{m}^{(2)} (κ ρ_{n + N}^{(h)})} {\tilde{E}}_{ϕ m + q}

(29)

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

Fig. 4. Three different types of the microstructures used to testing presented technique.

Fig. 5. Illustration of the computation domain - DFE method for some outer and inner rings is used.

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

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

] (10 multipole moments) and the Vector FDM-ABC Scheme (70 Azimuthal, 54 radial terms) [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]

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

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]

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

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

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

Table 1. The effective indices of the structure with 6 circular holes obtained by different methods

MM		FEM		FDM-ABC		FD-DFE
ℜ(n_eff)	-ℑ(n_eff)	ℜ(n_eff)	-ℑ(n_eff)	ℜ(n_eff)	-ℑ(n_eff)	ℜ(n_eff)	-ℑ(n_eff)
1.445395	3.19E-08	1.445394	4.11E-08	1.445395	3.07E-08	1.445391	3.21E-08
1.438584	5.31E-07	1.438576	3.97E-07	1.438589	5.43E-07	1.438585	5.09E-07
1.438445	9.73E-07	1.438442	7.13E-07	1.438444	9.62E-07	1.438437	9.04E-07
1.429957	1.59E-05	1.429952	1.67E-05			1.429954	1.57E-05
1.429248	8.73E-06	1.429261	9.17E-06			1.429264	8.83E-06

Fig. 6. The relative error of the calculated effective index: a)real part b)imaginary part. Current work - solid line, FD with rectangular mesh and PML [8

] - dashed line.

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μm, the outer radius r ₂ = 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.

Table 2. The effective indices of the structure with 3 angular-shaped holes obtained by different methods.

FEM		FDM-ABC		FD-DFE
ℜ(n_eff)	-ℑ(n_eff)	ℜ(n_eff)	-ℑ(n_eff)	ℜ(n_eff)	-ℑ(n_eff)
1.35580	4.95E-05	1.35584	5.00E-05	1.35588	5.11E-05
1.21479	1.11E-03	1.21476	1.25E-03	1.21488	1.25E-03

Table 3. The effective indices of the cored structure from Fig. 4(c).

FD-DFE
ℜ(n_eff)	-ℑ(n_eff)
1.435355	0
1.389851	1.29E-06
1.387892	5.44E-07
1.384163	4.13E-06

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]

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(060.2400) Fiber optics and optical communications : Fiber properties

ToC Category:
Research Papers

Citation
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/oe/abstract.cfm?URI=oe-13-25-10349

Sort: Journal | Reset

References

P. Russell, "Photonic Crystal Fibers," Science 299, 358-362 (2003). [CrossRef] [PubMed]
C. Kerbage, B. J. Eggleton, "Numerical analysis and experimental design of tunable birefringence in microstructured optical fiber," Opt. Express 10, 246-255 (2002), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-5-246">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-5-246</a>. [PubMed]
T.P. White, R.C. McPhedran, C.M. de Sterke, L.C. Botten, M.J. Steel, "Confinement losses in microstructured optical fibers," Opt. Lett. 26, 1660-1662 (2001). [CrossRef]
T.P. White, B.T. Kuhlmey, R.C. McPhedran, D.Maystre, R. Ranversez, C.M. de Sterke, L.C. Botten, M.J. Steel, "Multipole method for microstructured optical fibers. I. Formulation," J. Opt. Soc. Am. B 19, 2322-2330 (2002). [CrossRef]
D. Ferrarini, L. Vincetti, M. Zoboli, A. Cucinotta, and S. Selleri, "Leakage properties of photonic crystal fibers," Opt. Express 10, 1314-1319 (2002), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-23-1314">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-23-1314</a>. [PubMed]
H.P. Uranus, 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), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-12-2795">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-12-2795</a>. [CrossRef] [PubMed]
Z.Zhu, T.G. Brown, "Full-vectorial finite-difference analysis of microstructured optical fibers," Opt. Express 10, 853-864 (2002), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853</a>. [PubMed]
Shangping Guo, Feng Wu, Sacharia Albin, Hsiang Tai, Robert S. Rogowski,"Loss and dispersion analysis of microstructured fibers by finite-difference method," Opt. Express 12, 3341-3352 (2004), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-15-3341">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-15-3341</a>. [CrossRef] [PubMed]
J.P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Comput. Phys. 114, 185-200 (1994). [CrossRef]
J.P. Berenger, "Perfectly matched layer for the FDTD solution of wave-structure. Interaction problems," IEEE Trans. Microwave Theory Tech. 44, 110-117 (1996).
H. Rogier, 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]
A. Bayliss, M. Gunzburger, E. Turkel, "Boundary conditions for the numerical solution of elliptic equations in exterior regions," SIAM J. Appl. Math. 42, 430-451 (1982). [CrossRef]
T.Weiland, "VerlustbehafteteWellenleiter mit beliebiger. Randkontur und Materialverteilung," AEU 33, 170-174 (1979).
N. Kaneda, B. Houshmand, T. Itoh, "FDTD analysis of dielectric resonators with curved surfaces," IEEE Trans. Microwave Theory Tech. 45, 1645-1649 (1997). [CrossRef]
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]
M. Mrozowski, M. Okoniewski, M.A. Stuchly, "A hybrid PEE- FDTD method for efficient field modeling in cyllindrical coordinates," Electronics Letters 32, 194-195 (1996). [CrossRef]
M. Wiktor, M. Mrozowski,"Efficient Analysis of Waveguide Components Using a Hybrid PEE-FDFD Algorithm," IEEE Microwave and Wireless Components Letters 13, 396-398 (2003). [CrossRef]
C.D. Meyer, "Matrix analysis and applied linear algebra", SIAM, Philadelphia (2000).
M.Wiktor, M. Mrozowski,"Discrete Projection for Finite Difference Methods," 20th Annual Review or Progress in Applied Computational Electromagnetics, Syracuse NY 2004, Conf. proceedings, S04P08.
N.A. Issa, L. Poladian, "VectorWave Expansion Method for Leaky Modes of Microstructured Optical Fibers," J. Lightwave Technol. 21, 1005-1012 (2003). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper