Yee-mesh-based finite difference eigenmode solver with PML absorbing boundary conditions for optical waveguides and photonic crystal fibers

Abstract: Eigenvalue equations for solving full-vector modes of optical waveguides are formulated using Yee-mesh-based finite difference algorithms and incorporated with perfectly matched layer absorbing boundary conditions. The established method is thus able to calculate the complex propagation constants and the confinement losses of leaky waveguide modes. Proper matching of dielectric interface conditions through the Taylor series expansion of the fields is adopted in the formulation to achieve high numerical accuracy. The method is applied to the study of the holey fiber with triangular lattice, the two-core holey fiber, and the air-guiding photonic crystal fiber.


Introduction
The finite difference (FD) method has been popularly used for solving full-vector modes of optical waveguides.The method has the advantage of simple formulation and numerical implementation.It can be formulated either through the wave equations [1] or the Helmholtz equation [2].Recently, formulations based on the Yee's mesh, as in the finite-difference timedomain (FDTD) method [3], and derived directly from Maxwell's equations have been proposed [4][5][6].Such schemes are attractive in that the obtained mode fields can easily be incorporated into the FDTD computation.Although most of conventional optical waveguides are designed to propagate well-confined guided modes with ideal real-valued propagation constants, calculation of leaky properties of waveguide modes having complex propagation constants has become more important due to recent intensive study on the confinement losses of photonic crystal fibers (PCFs) [7,8].We have recently considered the Yee-mesh-based FD algorithm for formulating the eigenvalue problem for waveguide mode solutions and incorporate into it perfectly matched layer (PML) absorbing boundary conditions [9] so that leaky waveguide modes can be analyzed [10].Guo et al. reported a similar formulation using anisotropic PMLs [11].We however employed an alternative formulation for the PML by mapping Maxwell's equations into an anisotropic complex stretched coordinate [12][13][14].To obtain high-accuracy results, proper treatment of the fields near the dielectric interface is essential [15][16][17][18].For example, proper matching of interface conditions through the Taylor series expansion of the fields could achieve excellent accuracy [16].Similar interface matching procedure can be performed in the Yee-mesh-based algorithm [4,19].In [11] a dielectric constant averaging technique using Ampere's law across the curved material interface was proposed to increase numerical convergence and accuracy.In this paper we will apply the proper interface matching procedure to our formulation and several different PCF structures will be analyzed.
After presenting the formulation of the Yee-mesh-based FD method with PML absorbing boundary conditions and examining its numerical accuracy in Section 2, several microstructured fibers or PCFs including the holey fiber (HF) with triangular lattice, the twocore HF, and the air-guiding PCF will be analyzed using the established method in Section 3. The conclusion is drawn in Section 4.

Formulation
Figure 1 shows the cross-section of an arbitrary waveguide problem with the computing window surrounded by PML regions I, II, and III, each having thickness of d, with x and y being the transverse directions and z being the direction of propagation.Regions I and II have the normal vectors parallel to the x-axis and y-axis, respectively, and regions III are the four corner regions.In the PML region, using the stretched coordinate transform [12][13][14], Maxwell's equations can be written as where ω is the angular frequency, µ 0 and ε 0 are the permittivity and permeability of free space, respectively, and ε r is the relative permittivity of the medium considered.Using the modified differential operator ∇′ and the two-dimensional (2-D) Yee's mesh shown in Fig. 2 under the exp[-jβ z] field dependence assumption with β being the modal propagation constant, we have the Maxwell's curl equations in terms of the six components of the electric and magnetic fields for the whole problem ) where ε x , ε y , and ε z are the ε r values at the corresponding grid points of E x , E y , and E z , respectively.The values of s x and s y in Eqs.(3) are summarized in Table 1 with the parameter s defined as where σ e and σ m are the electric and magnetic conductivities of the PML, respectively, and n is the refractive index of the adjacent computing domain.In Eq. ( 4), we can see the impedance matching condition of the PML medium is which means that the wave impedance of a PML medium exactly equals to that of the adjacent medium in the computing window regardless of the angle of propagation.Assume that the electric conductivity of the PML medium has an m-power profile as where ρ is the distance from the beginning of the PML.At the interface of the PML and the computing window, the theoretical reflection coefficient for the normal incident wave is [13] max where c is the speed of light in free space.For the case of m = 2, s in Eq. (4) becomes ,( 1/ 2, ) ,( 1/ 2, ) , , ,( , )
After some mathematical work, an eigenvalue matrix equation in terms of the transverse magnetic fields can be obtained as where The eigenvalue equation ( 12) is solved by using the shift inverse power method.
As for the dielectric interface in a mesh as shown in Fig. 2, an improved scheme is proposed by using interpolation and extrapolation to approximate the fields on both sides of the dielectric interface [4,19].By properly matching the boundary conditions (BCs) at the dielectric interface, a modified FD formulation can be obtained.Figure 3 shows an example with a dielectric interface lying between two sampled grid points with E L and E R representing the electric fields on the left-hand and right-hand sides of the interface, respectively, and n being the normal unit vector which points outward from medium 1 to medium 2.
The first-order derivative of E y at point (i+1/2, j+1/2) with respect to x becomes { } where r x ∆x is the distance from the intersection of the dielectric interface along the x-axis to the point (i+1/2, j+1/2) and E y,L is the y component of E L .Considering the continuity BCs at the dielectric interface, we can have where n x and n y are the x and y components, respectively, of the normal unit vector n, and E y,R and E x,L can be approximated by interpolation and extrapolation as Similarly, for ∂H z /∂x at point (i+1, j+1/2), we have Applying the continuity BCs for the magnetic field, we can obtain Following the similar procedure, we can also obtain the derivatives of the electric and magnetic fields with respect to y for the case with a dielectric interface lying in a mesh.

Channel waveguide
We first consider a square channel waveguide discussed by Hadley [18] with the cross-section shown as the inset in Fig. 4. For the symmetric geometry of the waveguide, we only need to consider a quarter of the waveguide with the following parameters: waveguide width d = 1 µm, the indices n g = 2.83 and n ν = 1, and wavelength at 1.5 µm.For the purpose of comparing with the result of [18], instead of PMLs, we use the same BCs for the edges of the computing domain as in [18], i.e., perfect magnetic conductors (PMCs) placed on the top and bottom and perfect electric conductors (PECs) on the left and right sides of the computing domain.The modal index for the TE-like fundamental mode calculated by Hadley [18] was 2.65679692×10 -8 , which was obtained by dealing with the points near the dielectric interface and corners using a basis set series solutions of the Helmholtz equation in the cylindrical coordinate.
Applying the present method, the fundamental TE-like mode can be successfully obtained.As for the dielectric interface in the waveguide, in addition to the proper BC matching scheme discussed above and the simple stair-case approximation, an index averaging (IA) scheme which has been shown to be useful and efficient in increasing the numerical accuracy [5] is also adopted.In the IA scheme the value of the relative permittivity ε at each grid point is determined by an averaging formula.For example, in Fig. 2., ε z (i, j) is defined as where ε 1 and ε 2 are the relative permittivities of media 1 and 2, respectively, in the mesh and f is the filling fraction of medium 1 in the mesh.The relative errors in the modal index of our calculations to the above value reported by Hadley [18] versus the grid size with the three different interface treatments are shown in Fig. 4. It is seen that using the IA scheme and the proper BC matching scheme can achieve very high accuracy and the truncation errors are of higher order than that in the stair-case approximation.Relative errors in the modal index of the fundamental mode and the corresponding computation time for the strongly guiding optical fiber using the present method with the proper BC matching scheme and the stair-case approximation.

Optical fiber
We then consider a strongly guiding optical fiber with core radius being 3.0 µm, core index being 1.45, and air cladding.For the symmetry of the structure, we can only consider a quarter of the fiber.The size of the computing window is 6.0 µm × 6.0 µm, and the zero-value BC is applied at the edges.Figure 5 shows the calculated relative errors in the fundamental modal index with respect to the analytical solution n eff = 1.438604 and the corresponding computation time based on a Pentium IV 3.0-GHz personal computer using different numbers of the grid points along the x-axis.The lines with circular dots are the results using the proper BC matching scheme.One can see that the relative errors rapidly decrease to the order of 10 -6 as the number of grid points is larger than 40.On the contrary, the results adopting stair-case approximation, which is the lines with square dots, possess slower convergent property.By comparing the computation time, we observe that the results obtained using the proper BC matching scheme need more grid points to approximate the modal field, and thus denser characteristic matrix in the formulation and longer computation time are required than when stair-case approximation is used.

Holey fibers with triangular lattice
We now consider the loss properties of a HF with triangular lattice having three-ring air holes in the cladding region, which corresponds to 36 air holes.The cross-section is illustrated in Fig. 6 with the size of the computing window being W x × W y and PMLs with thickness d being placed on the top and the right side of the window.In the following calculation we adopt the parameters used in [20]: a = 2.3 µm; W x = 8.0 µm; W y = 7.3 µm; d = 2 µm; R = 10 -8 ; silica index n = 1.45.(Silica material is assumed in all the PCFs analyzed in this paper.)In [20] Saitoh and Koshiba proposed a finite-element imaginary-distance beam propagation method (FE-ID-BPM) for studying PCFs.Our calculated effective indices and losses for the xpolarized fundamental guided modes of the three-ring HF with r/a = 0.25, 0.3, and 0.35 for λ = 0.1 to 2.0 µm are illustrated in Fig. 7(a) and (b), respectively, by the solid lines.The effective index in this case decreases with the increase of the wavelength or the hole size.From the magnitude of the loss, we can see that larger air holes can provide better confinement resulting in smaller loss.Our results agree very well with those from Saitoh and Koshiba [20] shown in Fig. 7(a) and (b) by the circular dots.We have also analyzed the onering HF (with 6 air holes) and the two-ring HF (with 18 air holes) and the same good agreement with [20] is achieved.We like to remark here that for obtaining the results of Fig. 7, calculation with the staircase approximation is good enough.However, when analysis of modal birefringence defined as the difference between the x-polarized and y-polarized modal effective indices of the fundamental modes, the more rigorous proper BC matching scheme has to be employed in order to obtain better numerical convergence, For example, Koshiba and Saitoh have numerically examined the fundamental mode degeneracy of the one-ring HF using their FEM method [21].We study the same structure (a = 6.75 µm and r = 2.5 µm) and the modal birefringence as a function of the number of grid points along the x-axis is shown in Fig. 8.It is seen that using the proper BC matching scheme, the numerical convergence appears to be much better than only using the stair-case approximation.We also list in Table 2 the calculated modal indices with variant numbers of grid points along the x-axis using the PMLs and the proper BC matching scheme.Our results are seen quite close to the reference result n eff,ref = 1.445395345-3.15×10 - j [8] which was obtained by using the multipole method.We then consider a kind of two-core HFs having three-ring air holes surrounding the two cores with the computing window illustrated in Fig. 9.The PMLs with d = 2 µm and R = 10 -8 are placed on the top and the right side of the computing window to gain the loss information of the two-core HF. Figure 10(a) and (b) shows the effective indices and the losses of the xpolarized even mode for the two-core PCF with the hole radius r being 0.25a, 0.3a, and 0.35a calculated at different wavelengths.It demonstrates that the effective indices and the losses can be successfully obtained by using the present method with PMLs.From Fig. 10(b) one can see that smaller air holes cause weaker confinement of light, resulting in larger loss than larger air holes would do.Also, Fig. 10 shows that at longer wavelengths, smaller effectiveindex values and larger losses in guided modes are obtained.These general trends in the fiber characteristics are similar to those given in Fig. 7 for the one-core HF.By comparing Fig.The x-polarized odd mode has similar properties in the effective index and the loss, as shown in Fig. 11(a) and (b), except that the losses of the odd modes are slightly smaller than those of the even mode.

Air-guiding photonic crystal fFiber
PCFs with air cores have been attractive structures with many promising applications [22].Consider the air-guiding PCF or hollow PCF with the cladding region formed by six rings of air holes with pitch a being 2 µm and air filling fraction f in a unit cell being 0.7, as shown in Fig. 12.Using our FD method with PMLs having d = 2 µm and R = 10 -8 as indicated in Fig. 12, we can find the guided modes within the photonic band gap (PBG) or near the edges of the PBG.The computing window size is 17 µm × 15.2 µm, and for a 170 × 152 grid division, the computation time is 267 seconds for one wavelength using a Pentium IV 3.0-GHz personal computer.The effective indices of the fundamental x-polarized guided modes for this hollow PCF are plotted in Fig. 13(a) with the dashed lines representing the PBG boundaries.It can be observed that there are two modal dispersion curves with flat and steep slopes, respectively, which correspond to the fundamental core modes and surface modes [23] of the PCF. Figure 14 shows the field distributions of the guided modes at points A, B, C, and D along the modal dispersion curves as indicated in Fig. 13(a).At points B and C, which are located on the flat curve representing the fundamental core modes, one can see that stronger field (colored white) is located in the air core with weaker field (colored black) in the silica around the air core.The other kind of guided modes is the surface modes, labeled A and D, with most of the energy concentrated in the silica region around the air core.Figure 13(b) illustrates the losses of these two kinds of modes within the frequency range of the PBG.Outside the PBG or near the edges of the PBG, the core modes have more field penetrating in the cladding, resulting in larger losses than those of the core modes lying in the central part of the PBG.As for the surface modes, the losses remain almost the same magnitude and larger than those of the core modes at most frequencies.Ideally, the hollow PCFs should be lossless due to the much smaller absorption and Rayleigh scattering of the air.However, recent research found that some losses in the hollow PCF may result from the coupling between the core modes, the surface modes, and the leaky cladding modes [23].Although most of the core mode field is concentrated in the air core as shown in Fig. 14, the overlap of field distributions of core modes and surface modes in the silica region around the air core makes the low-loss core modes coupled to the surface modes having larger losses.Compared with the core modes, the surface modes have more overlap field with the leaky cladding m odes in the cladding region, and thus easily couple to these modes, yielding the high loss in the hollow PCF.It has been reported that careful design of the air-core size can help produce a single-mode hollow PCF free of surface modes [24].

Conclusion
We have formulated eigenvalue equations for solving full-vector modes of optical waveguides using Yee-mesh-based finite difference algorithms with PML absorbing boundary conditions.We have made use of the formulation for the PML in which Maxwell's equations are mapped into an anisotropic complex stretched coordinate [12][13][14].Proper matching of dielectric interface conditions through the Taylor series expansion of the fields is adopted in the formulation to obtain better numerical accuracy [4,19].The established method is found to be able to achieve high accuracy through the analysis of a square channel waveguide with reported accurate modal index [18].The method is shown to be capable of providing the complex propagation constants and the confinement losses of leaky waveguide modes.The method has been successfully applied to the analysis of the holey fiber with triangular lattice, the two-core holey fiber, and the air-guiding photonic crystal fiber.The efficient calculation of the propagation loss characteristics and mode field profiles of the core modes and surface modes in the air-guiding PCF is demonstrated.

Fig. 1 .
Fig. 1.The cross-section of an arbitrary waveguide problem with the PMLs placed at the edges of the computing domain.
central difference scheme for the differential operators in Eqs.(3), the timeharmonic Maxwell curl equations in Eqs.(1) can be converted into the matrix form

2 .
and D y are square matrices determined by the central difference scheme and the boundary conditions.Please note that the subscript (i, j) in the left hand side of each equation of Eqs.(11) denotes the mesh number in the 2-D computing domain while those in the right hand side represent the corresponding grid (point) numbers as shown in Fig.If no PML region is applied as the boundary conditions at the edges of the computing window, we can have A x = B x , A y = B y , C x = D x , and C y = D y , and Eqs. (

Fig. 3 .
Fig.3.The situation of a dielectric interface lying between two sampled points.

Fig. 4 .
Fig.4.Relative errors in the modal index of the fundamental TE-like mode for the square channel waveguide using the present method with three different schemes in dealing with the dielectric interface.
Fig.5.Relative errors in the modal index of the fundamental mode and the corresponding computation time for the strongly guiding optical fiber using the present method with the proper BC matching scheme and the stair-case approximation.

Fig. 6 .
Fig.6.The computing window with PMLs for the holey fiber with three-ring air holes.

Fig. 7 .Fig. 8 .
Fig. 7. (a) The effective indices and (b) the losses of the x-polarized fundamental guided modes in the three-ring holey fiber with a = 2.3 µm.

rFig. 10 .
Fig. 10.(a) The effective indices and (b) the losses of the even modes in the two-core holey fiber with r/a = 0.25, 0.30, and 0.35.

r
Fig. 11.(a) The effective indices and (b) the losses of the odd modes in the two-core holey fiber with r/a = 0.25, 0.30, and 0.35.

Fig. 12 .
Fig. 12.The computing window with PMLs for a hollow PCF with six rings of air holes in the cladding.

Fig. 13 .
Fig. 13.(a) The effective indices of the surface modes and fundamental core modes of the hollow PCF of Fig. 11 with air filling fraction f = 0.7.(b) Losses of the core modes and the surface modes.

Fig. 14 .
Fig. 14.The field distributions of the surface modes and fundamental core modes at points A, B, C, and D indicated in Fig. 13(a).

Table 1 .
The values of s x and s y for the PML region I, II, and III.

Table 2 .
[8] calculated modal index as a function of the number of grid points along the x-axis and the reference result is n eff,ref = 1.445395345-3.15×10 -[8].
PMLPML Fig.9.The computing window with PMLs for the two-core holey fiber.