Open Access
Add to BibSonomyAbstract
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.
©2004 Optical Society of America
1. 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 time-domain (FDTD) method [3], and derived directly from Maxwell’s equations have been proposed [4–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–14]. To obtain high-accuracy results, proper treatment of the fields near the dielectric interface is essential [15–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 two-core HF, and the air-guiding PCF will be analyzed using the established method in Section 3. The conclusion is drawn in Section 4.
2. 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–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 Ex, Ey, and Ez, respectively. The values of sx and sy 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]
Table 1. The values of sx and sy for the PML region I, II, and III.
and the maximum conductivity σ max can then be determined as
where c is the speed of light in free space. For the case of m=2, s in Eq. (4) becomes
Fig. 1. The cross-section of an arbitrary waveguide problem with the PMLs placed at the edges of the computing domain.
By applying the central difference scheme for the differential operators in Eqs. (3), the time-harmonic Maxwell curl equations in Eqs. (1) can be converted into the matrix form
where
Δx and Δy are grid sizes, I is an identity matrix, and ε x, ε y, and εz are diagonal matrices representing relative permittivities εx, εy, and εz, respectively, at the corresponding grid points. A x, A y, B x, B y, C x, C y, D x, 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. 2. 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. (10) will appear with the same form as in [5].
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 EL and ER 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 Ey at point (i+1/2, j+1/2) with respect to x becomes
where rxΔ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 Ey,L is the y component of EL. Considering the continuity BCs at the dielectric interface, we can have
where nx and ny are the x and y components, respectively, of the normal unit vector n, and Ey,R and Ex,L can be approximated by interpolation and extrapolation as
Similarly, for ∂Hz/∂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.
3. Numerical results
3.1 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 ng=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.
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.
3.2 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 neff=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.
3.3 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 Wx×Wy 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; Wx=8.0 µm; Wy=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 x-polarized 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 one-ring HF (with 6 air holes) and the two-ring HF (with 18 air holes) and the same good agreement with [20] is achieved.
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.
Fig. 8. Modal birefringence of the one-ring HF using the present method with the proper BC matching scheme and the stair-case approximation.
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-8 j [8] which was obtained by using the multipole method.
Table 2. The calculated modal index as a function of the number of grid points along the x-axis and the reference result is neff,ref=1.445395345-3.15×10-8j [8].
3.4 Two-core holey fibers
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 x-polarized 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 effective-index 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. 10(b) with Fig. 7(b), we observe that the confinement losses of the two-core HF are a little smaller than those of the one-core HF.
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.
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.
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.
3.5 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.
Fig. 12. The computing window with PMLs for a hollow PCF with six rings of air holes in the cladding.
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 modes 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].
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. The field distributions of the surface modes and fundamental core modes at points A, B, C, and D indicated in Fig. 13(a).
4. 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–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.
Acknowledgments
This work was supported in part by the National Science Council of the Republic of China under grants NSC92-2215-E-002-004 and NSC92-2215-E-002-008 and in part by the Ministry of Economic Affairs of the Republic of China under grant 91-EC-17-B-08-S1-0015.
References and links
1. C. L. Xu, W. P. Huang, M. S. Stern, and S. K. Chaudhuri, “Full-vectorial mode calculations by finite difference method,” Inst. Elec. Eng. Proc.-J 141, 281–286 (1994). [CrossRef]
2. P. Lüsse, P. Stuwe, J. Schüle, and H.-G. Unger, “Analysis of vectorial mode fields in optical waveguides by a new finite difference method,” J. Lightwave Technol. 12, 487–494 (1994). [CrossRef]
3. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. AP-14, 302–307 (1966).
4. T. Ando, H. Nakayama, S. Numata, J. Yamauchi, and H. Nakano, “Eigenmode analysis of optical waveguides by a Yee-mesh-based imaginary-distance propagation method for an arbitrary dielectric interface,” J. Lightwave Technol. 20, 1627–1634 (2002). [CrossRef]
5. 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 [CrossRef] [PubMed]
6. H. C. Chang and C. P. Yu, “Yee-mesh-based finite difference eigenmode analysis algorithms for optical waveguides and photonic crystals,” in OSA Integr. Photon. Res. Dig., 2004, paper IFE4.
7. J. C. Knight, T. A. Birks, P. St. J. Russell, and D. M. Atkin, “All-silica single-mode optical fiber with photonic crystal cladding,” Opt. Lett. 21, 1547–1549 (1996). [CrossRef] [PubMed]
8. T. P. White, B. T. Kuhlmey, R. C. Mcphedran, D. Maystre, G. Renversez, C. M. de Sterke, and L. C. Botten, “Multipole method for microstructured optical fibers. I. Formulation,” J. Opt. Soc. Am. B 19, 2322–2330 (2002). [CrossRef]
9. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Computat. Phys. 114, 185–200 (1994). [CrossRef]
10. H. C. Chang and C. P. Yu, “Yee-mesh-based finite difference optical waveguide eigenmode solver with perfectly matched layer absorbing boundaries conditions,” in Proceedings of 5th Asia-Pacific Engineering Research Forum on Microwaves and Electromagnetic Theory, pp. 158–164, Kyushu University, Fukuoka, Japan, July 29–30, 2004. [PubMed]
11. S. Guo, F. Wu, and S. Albin, “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]
12. W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave Opt. Technol. Lett. 7, 599–604 (1994). [CrossRef]
13. R. Mittra and U. Pekel, “A new look at the perfectly matched layer (PML) concept for the reflectionless absorption of electromagnetic waves,” IEEE Microwave Guided Wave Lett. 5, 84–86 (1995). [CrossRef]
14. C. M. Rappaport, “Perfectly matched absorbing boundary conditions based on anisotropic lossy mapping of space,” IEEE Microwave Guided Wave Lett. 5, 90–92 (1995). [CrossRef]
15. Y. P. Chiou, Y. C. Chiang, and H. C. Chang, “Improved three-point formulas considering the interface conditions in the finite-difference analysis of step-index optical devices,” J. Lightwave Technol. 18, 243–251 (2000). [CrossRef]
16. Y. C. Chiang, Y. P. Chiou, and H. C. Chang, “Improved full-vectorial finite-difference mode solver for optical waveguides with step-index profiles,” J. Lightwave Technol. 20, 1609–1618 (2002). [CrossRef]
17. G. R. Hadley, “High-accuracy finite-difference equations for dielectric waveguide analysis I: Uniform regions and dielectric interfaces,” J. Lightwave Technol. 20, 1210–1218 (2002). [CrossRef]
18. G. R. Hadley, “High-accuracy finite-difference equations for dielectric waveguide analysis II: Dielectric corners,” J. Lightwave Technol. 20, 1219–1231 (2002). [CrossRef]
19. Ditkowski, J. S. Hesthaven, and C. H. Teng, “Modeling dielectric interfaces in the FDTD-method: A comparative study,” in 2000 Progress in Electromagnetics Research (PIERS 2000) Proceedings, Cambridge, Massachusetts, July 5–14, 2000.
20. K. Saitoh and M. Koshiba, “Full-vectorial imaginary-distance beam propagation method based on a finite element scheme: Application to photonic crystal fibers,” IEEE J. Quantum Electron. 38, 927–933 (2002). [CrossRef]
21. M. Koshiba and K. Saitoh, “Numerical verification of degeneracy in hexagonal photonic crystal fibers,” IEEE Photon. Technol. Lett. 13, 1313–1315 (2001). [CrossRef]
22. R. F. Cregan, B. J. Mangan, J. C. Knight, T. A. Birks, P. St. J. Russell, P. J. Roberts, and D. C. Allan, “Single-mode photonic band gap guidance of light in air,” Science 285, 1537–1539 (1999). [CrossRef] [PubMed]
23. K. Saitoh, N. A. Mortensen, and M. Koshiba, “Air-core photonic band-gap fibers: the impact of surface modes,” Opt. Express 12, 394–400 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-3-394 [CrossRef] [PubMed]
24. H. K. Kim, J. Shin, S. Fan, M. J. F. Digonnet, and G. S. Kino, “Designing air-core photonic-bandgap fibers free of surface modes,” IEEE J. Quantum Electron. 40, 551–556 (2004). [CrossRef]
