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
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 ). The second category (proposed by Berenger ) 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 . 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 ) is about -40dB for TBC and about -120dB for PML. However, TBC are still applied in other techniques (beam propagation method  or finite element method ).
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 , but it is applicable to Yee’s mesh in rectangular coordinates and relies on rigorous rather than approximate relationships for the outgoing waves.
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 TEx and TMx modes.
Let us consider the equation for TEx modes
where em =Ey(mΔx), , ε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 Rzy (e) , Ryz (H) and ε are square M × M matrices
and e = [e 1,e 2,⋯,eM]T, h = [h 1,h 2,⋯, hM]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
where κi = [-γ 2- k2 0 εi]1/2, k 0 = ω(μ 0 ε 0)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 .
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 =Ryz (h) Rzy (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(γ) =γ-, where is an eigenvalue of the problem Ã(γ)ẽ = 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 Ez(ρ,φ,z) and Hz(ρ, φ,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 Hm (2) (∙) is a Hankel function of the second kind, κ = [γ 2 + k 2 0 με]1/2, η = [μ/ε)1/2 and Am, Bm, Cm, Dm 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 Ex(x,y,z) and Ey(x,y,z) are needed, hence
where x = ρ cos(φ) and y = ρ sin(φ).
Substituting the above relations and taking z = 0, we get
Using the above series one can find the relation between the field samples outside and inside the computational domain.
Let SC 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 SB be a set of points on the boundary (crosses in fig. 2). Let us denote by ECt and EBt the field values from set SC and SB, respectively. From relations (23) and (24) we get
where vector C = [A 0,... ,AQ,B 0,... ,BQ,C 0,... ,CQ,D 0,...,DQ]T. The relation between field values on the boundary and inside the computational domain is
where MC (inv) = MC -1 if the number of the SC elements is equal to the size of C vector, otherwise to invert MC the singular value decomposition (SVD) algorithm must be applied  (MC (inv) MC = 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 MB and MC (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: 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  was used.
The effective indices
Also in structures with larger contrast good agreement was achieved. The results for ε r1 = ε r3 = 9, ε r2 = 1, b = 1nm, λ 0 = 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 neff 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 . The radius of the fiber is R = 0.5μm, the core index nc = 2.9, the background index nbg = 1. 55 and λ 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 RB = 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.
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 nbg = 1.45 and na = 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.
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)  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 nbg = 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  (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 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
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 ).
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 . 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).
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.
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]