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 -40dB for TBC and about -120dB 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^x modes.

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

 

Fig. 1. A cross section of the slab dielectric waveguide covered with Yee’s mesh.

Download Full Size | PDF

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 ₁,e ₂,⋯,e_M]^T, h = [h ₁,h ₂,⋯, 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 (ε ₂ = ε ₂(x)).

The field shape in region 1 and 3 is known

where κ_i = [-γ ²- k² ₀ ε_i]^1/2, k ₀ = ω(μ ₀ ε ₀)^1/2 and A _-,A ₊ are unknown coefficients. For evanescent modes κ_i, 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].

Expressions (7) determine the relations between the field samples at the boundary of the computational domain

where a _- = e^-κ₁Δx and a ₊ = e^-κ₃Δx.

Let us rewrite problem (3) in a simple form

where A =R_yz ^(h) R_zy ^(e) - ωμ ₀ ε, 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 γ ₀ and denote the solution as γ ₁:

In the next step we have to solve the linear problem (13) for γ ₂, substituting γ ₁ 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 }$ ² ẽ. 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 ⁽²⁾ (∙) is a Hankel function of the second kind, κ = [γ ² + k ² ₀ με]^1/2, η = [μ/ε)^1/2 and A_m, B_m, C_m, D_m are arbitrary coefficients. In practice we have to reduce the series to a finite number of terms, so let us assume that 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

and

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

 

Fig. 2. Tangential electric field of Yee’s mesh. A fragment of the computational domain in a neighborhood of the boundary.

Download Full Size | PDF

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^B be a set of points on the boundary (crosses in fig. 2). Let us denote by E^C_t and E^B_t the field values from set S^C and S_B, respectively. From relations (23) and (24) we get

and

where vector C = [A ₀,... ,A_Q,B ₀,... ,B_Q,C ₀,... ,C_Q,D ₀,...,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^C elements is equal to the size of C vector, otherwise to invert M_C the singular value decomposition (SVD) algorithm must be applied [11] (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 = 1nm, λ ₀ = 0.2nm are presented in Table 1 (corresponding fields distributions are shown in fig. 3 – 5).

 

Fig. 3. The field distribution for n_eff = 0.99526-0.00134j.

Download Full Size | PDF

Also in structures with larger contrast good agreement was achieved. The results for ε _r1 = ε _r3 = 9, ε _r2 = 1, b = 1nm, λ ₀ = 1.5nm 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).

Table 1. The results obtained for a structure with parameters ε _r1 = ε _r3 = 1.21, ε _r2 = 1, b = 1nm, λ₀ = 0.2nm. The percentage error of the real and imaginary part is given in brackets.

View Table | View all tables in this article

 

Fig. 4. The field distribution for n_eff = 0.98086-0.00533j.

Download Full Size | PDF

 

Fig. 5. The field distribution for n_eff = 0.95621-0.01196j.

Download Full Size | PDF

Table 2. The results obtained for a structure with parameters ε _r1 = ε _r3 = 9, ε _r2 = 1, b= 1nm for λ ₀ = 1.5nm. The percentage error of the real and imaginary part is given in brackets.

View Table | View all tables in this article

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_bg = 1. 55 and λ ₀ = 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_B = 0.55μ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.

 

Fig. 6. The computational domain for PML technique - left hand side and for analytical boundary conditions (ABC) - right hand side.

Download Full Size | PDF

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_a = 1 for holes. The vacuum wavelength used in calculation is λ = 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.

 

Fig. 7. Two different types of a microstructured optical fibers.

Download Full Size | PDF

Table 3. The results obtained for an analysis of the optical fiber. The percentage error of the real and imaginary part is given in brackets.

View Table | View all tables in this article

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

Table 4. The real part of the effective indices of the structure with 6 circular holes obtained by different methods. The percentage error (in brackets) is calculated in reference to Multipole Method (MM) [13] (10 multipole moments).

View Table | View all tables in this article

Table 5. The imaginary part of the effective indices of the structure with 6 circular holes obtained by different methods. The percentage error (in brackets) is calculated in reference to Multipole Method.

View Table | View all tables in this article

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 ₂₁ - like mode.

Table 6. The real part of the effective indices of the structure with 3 angular-shaped holes obtained by different methods. The percentage error (in brackets) is calculated in reference to Vector FDM-ABC Scheme (70 Azimuthal, 54 radial terms) [14].

View Table | View all tables in this article

Table 7. The imaginary part of the effective indices of the structure with 3 angular-shaped holes obtained by different methods. The percentage error (in brackets) is calculated in reference to FDM-ABC method.

View Table | View all tables in this article

Table 8. The convergence of the fundamental mode for structure presented in fig. 7a.

View Table | View all tables in this article

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 120s for PML and 30s for the proposed algorithm (on a standard 2.8GHz desktop). For the microstructured fibers (only partial use of boundary conditions) the difference was smaller (70s versus 40s and 140s vs. 90s).

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

 

Fig. 8. Berenger and leaky modes for structure presented in fig. 7a. Leaky modes are surrounded by a circle.

Download Full Size | PDF

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.