## Abstract

The dispersion and loss in microstructured fibers are studied using a full-vectorial compact-2D finite-difference method in frequency-domain. This method solves a standard eigen-value problem from the Maxwell’s equations directly and obtains complex propagation constants of the modes using anisotropic perfectly matched layers. A dielectric constant averaging technique using Ampere’s law across the curved media interface is presented. Both the real and the imaginary parts of the complex propagation constant can be obtained with a high accuracy and fast convergence. Material loss, dispersion and spurious modes are also discussed.

© 2004 Optical Society of America

## 1. Introduction

Microstructured fibers have been under intensive study as they offer design flexibility in controlling the mode propagation properties. These fibers have some extraordinary properties, such as endless single mode propagation, special dispersion, and high or low nonlinear effects. Air-guiding fibers such as the air-core photonic band gap fibers or Bragg fibers are also of considerable interest. Reviews of photonic crystal fibers can be found in [1–3] and the references therein.

In general, the modeling of these fibers takes advantage of the fact that the E or H field can be decomposed into longitudinal and transverse components in waveguides with invariant index profiles along z-direction. The field can be written as:

where ξ denotes E or H field and the subscripts *t* and *z* denote respectively the transverse and longitudinal components.

The full-vectorial Helmholtz equations can be obtained by substituting Eq. (1) into Maxwell’s equations:

The Helmholtz equations give several important conclusions: it is possible to form an eigen-value problem using the transverse E or H components since Eqs. (2a–b) only contain the transverse components. However, this is not true for the z-components since they are coupled to the transverse components. When the index change is small and the coupled items on the right-hand side are omitted, (which is the scalar approximation), all four equations become the same, and an eigen-value problem can be formed for any E or H components.

In the holey fibers or microstructured fibers, the index contrast of the materials is generally high, hence the scalar wave analysis methods are not accurate to predict their propagation properties; a full-vector approach is required. So far, a few full-vector methods have been used to characterize microstructured fibers, such as the plane wave expansion method (PWM) [4–8], localized function method [9–12], beam propagation method [13–16], finite-element method (FEM) [17–22], and finite-difference method (FDM) in time domain [23–28] or frequency domain [29–31]. Specifically, a highly accurate semi-analytical multipole method [32] has been developed to model fibers with circular air hole inclusions. A brief review of these methods is given below.

The PWM is an extension of modeling for photonic crystals. It assumes an infinite periodic index profile and treats the unit cell or supercell by applying the Bloch boundary conditions. The eigen-matrix is a full matrix and the complexity is the same as the PWM for photonic crystals. It is adequate for index-guiding fibers with many periodic air holes in the cladding, however, the artificial periodic boundary condition and supercell approach are not very suitable for many real fiber structures with a finite number of air holes.

The localized function method based on Galerkin method has been widely used for waveguide analysis, both scalar and vectorial problems. This method applies a set of localized orthogonal functions, such as Sine [33], Laguerre-Gauss [34, 35] (for 1D waveguide), and Hermite-Gauss [6, 9, 10, 12] (for 2D problem), to approximate the unknown mode fields of the localized modes. When the mode is far away from cut-off or well confined, the mode fields can be approximated using tens or hundreds of functions. These methods generally involve integrations, which are computation intensive, and the convergence is generally a problem.

Finite-element method [17–22] is a powerful numerical tool for waveguide problems. It often combines the beam propagation method or simply solves the Helmholtz equation in frequency domain by discretizing the region of interest into triangular cells, which is able to represent fine curved structures. Hence, FEM can provide high accuracy, but the complexity of the algorithm implementation and the computation is generally high.

Finite-difference method using Yee’s mesh [36] is popular for electromagnetic problems, and a compact-2D scheme is often used for waveguides. Compared to FEM, FDM is much easier to implement; yet, it is able to offer a comparable accuracy. The FDM approaches include time domain (FDTD) and frequency domain (FDFD) methods. The compact-2D FDTD [23–28] method solves the eigen-frequency for a given propagation constant, and therefore is unable to process material dispersion. The FDFD approaches include those based on solving the Helmholtz equation [29, 30] or Maxwell’s equations [31] directly. Both material dispersion and material loss are easy to incorporate in FDFD. The latter approach is appealing due to its many merits. There are no second order derivatives; the method is fast and accurate, mathematically simple and straightforward; all six fields are obtained from the transverse E or H field; different boundaries such as Bloch boundary for photonic band gap calculations [37–39] can be readily applied; and sparse techniques can be used to reduce computation.

In practical holey fibers with air hole inclusions, the confinement (by either index guiding or PBG guiding) is not perfect due to the finite cladding; hence, the confinement loss is a significant characteristic of the microstrucured fiber. It has been calculated using several methods, including the semi-analytical multipole method [32], Fourier expansion method [40], and FEM with anisotropic perfectly matched layers (PMLs) [19]. PWM and Hermite-Gauss methods have not been used for leakage loss calculations.

In this paper, we use the compact-2D FDFD approach described in [31] for optical waveguides based on solving the Maxwell’s equations directly to include the calculation of mode leakage loss due to the finite cladding. By applying the anisotropic PML layers, both the dispersion and loss properties can be evaluated in a single run. The method preserves all the advantages discussed above and is simple and straightforward. The curved profile is studied carefully to increase both the accuracy and convergence of the complex propagation constant.

The analysis method is given in section 2; numerical results for a PCF with a single ring of air holes are analyzed in section 3 along with the averaging technique at the media interface, and finally conclusions are given in section 4.

## 2. Analysis method

The leakage loss of a mode can be represented by the imaginary part of its complex propagation constant. To model the leakage, an open boundary condition has to be used, which produces no reflection at the boundary. The PMLs are so far the most efficient absorption boundary condition for this purpose. The split-field PML proposed by Berenger [41] which is often used in FDTD algorithm cannot be applied in frequency domain methods such as the FEM and FDFD methods since it introduces non-Maxwellian equations. The equivalent nonsplit-field anisotropic PML [13, 42–47] has been proposed instead to simulate the open environment in these applications. The Maxwell equations for optical waveguide with anisotropic PML boundaries are expressed as:

For exp(-jωt) convention which is used in this paper:

where σ is the conductivity profile.

In the compact-2D scheme for waveguides, the z-derivatives are replaced by analytical expressions using Eq. (1), and other derivatives are replaced by finite differences in Yee’s mesh. Therefore, the curl equations (3) can be expressed in a matrix form:

where the U and V are sparse matrices which are obtained in the same way as in [31, 37] using a zero boundary outside of the PML layers. There is no need to treat PML in a special way as in split-field scheme.

For simplicity, we assume:

Substituting Eq. (6a) into Eq. (6b) and eliminating H_{z} as in [31], an eigen-value problem can be obtained for H_{t}:

and similarly for E_{t}:

where:

and

In the absence of PML media, the equations given above are reduced to those obtained in [31]. Hence, we can validate our results by setting the thickness of PML layers to zero.

The waveguide is encompassed by PML layers followed by a layer of perfect electric conductor (PEC) or zero boundary. The modes leaking out of the fiber will be absorbed efficiently by the PML with very small reflections; hence, the effect of the artificial boundary on the modes in the PCF will be minimized (especially for those well-localized guided modes).

## 3. Numerical results

The PCF example in [32] is used in this paper since the analytical results using multipole method are available for comparison. The PCF has a single ring of air holes (6 holes) in glass fiber. The parameters used are: the lattice constant *a*=6.75µm, the air hole radius *r*=2.5µm, and the refractive index of the glass is 1.45. The material dispersion is omitted since it is trivial to include it in a compact-2D FDFD scheme. According to the multipole method, the accurate effective mode index at wavelength 1.45µm would be 1.445395345+3.15×10^{-8}i.

This PCF has a symmetry of C_{6ν} (six-fold rotation symmetry and at least one plane of reflection symmetry), and the computation region can be reduced using the symmetry properties by applying a combination of PEC and PMC (perfect magnetic conductor) [48]. The PCFs with such symmetry supports eight mode classes. Figure 1 shows the fiber profile and a quarter of the whole region used for computation of the third and fourth mode classes, which are degenerate pairs with a 90-degree rotation symmetry including the fundamental mode. Glass material is assumed to extend uniformly to infinity, and PML layers are used outside the computation region with a 2^{nd}-order power law profile. The computation region is chosen to be 1.5*a* along both x and y directions and the thickness of the PML layers is 10% of the thickness of the inside area along x or y direction.

The computation region is discretized by a 2D Yee’s mesh. The curved interface crossing a cell is generally approximated using a staircase scheme or averaged using the effective index scheme. The averaging technique is shown to be very effective in increasing convergence and accuracy as in [31]. Considering two different media in the cell, the average dielectric constant in the cell can be evaluated as:

where *f* is the fraction of the first material ε_{a}.

Figure 2 shows the relative error of the calculated complex effective mode index of the fundamental mode with different number of grids. The real part of the effective mode index converges quickly to an accuracy of 10^{-5}~10^{-6} even if a coarse mesh (for example, 30×30) is used; the accuracy is sufficient to obtain group velocity dispersion and other parameters. However, the imaginary part converges rather slowly with a relative error in the range of 10^{-1}~10^{-2} with almost 30% relative error using a 30×30 mesh. Though it converges to the true value with a fine mesh, the slow convergence is still not satisfactory.

We found that the slow convergence of the imaginary part is due to the improper averaging technique at the interface of the two dielectric materials. It can be greatly improved using a more reasonable averaging technique at the interface. Our averaging scheme is shown in Fig. 3. ε_{rx}, ε_{ry} and ε_{rz} are the averaged dielectric constant of the cell located at the same position as E_{x}, E_{y} and E_{z}.

First, we check the averaging technique using Ampere’s Law (the curl equation):

The integration area is the cell formed by the four surrounding H components as shown in Fig. 4. It shows the intersection of the curved interface on the integration surface. In each of these surfaces, the E component is tangential to the interface and will be continuous across it, and the average of the E field is used for the value at the cell center. Taking E_{z} as an example:

For ε_{x} in yz plane and ε_{y} in xz plane, the averaging is easy to do since the boundary is parallel to z direction and the integration cell shrinks to a line, which is shown in Fig. 3 as the dotted cells.

E_{x} and E_{y} in the xy plane as shown in Fig. 3 are not tangential to the dielectric interface, and therefore will not be continuous across the dielectric boundary. As in Fig. 4, when the integration surface for E_{x} moves along x in the cell on the xy plane, the interface will shift, and similar is the case for E_{y}. Since E_{x} and E_{y} are the average field values of the cell in xy plane, another average has to be taken and the averaged dielectric values are [28]:

The average rule described in Eqs. (13, 16–17) could also be derived as the classical rules for the evaluation of effective dielectric constants as described in Milton’s book [49]. The integration can be approximated well using a denser subcell mesh along x or y as the two dashed lines shown in Fig. 3. The same averaging procedure can be applied to magnetic materials with a profile of µ_{r}(r) and our FDFD algorithm is also applicable for these materials.

Since our averaging technique (Eqs. (15–17)) satisfies Ampere’s Law everywhere across the boundary, we introduced it to reduce the possible spurious modes by finite difference scheme. Surprisingly, it is found to be very effective to improve the accuracy and convergence of the imaginary part. The same fiber is calculated again with this averaging technique and the results are shown in Fig. 5 and Table 1.

The convergence and accuracy of the real part shown in the Fig. 5 are similar to those in Fig.2 which was obtained using the previous averaging technique. However, the convergence and the accuracy of the imaginary part are increased at least by one order of magnitude, with an error of 2% for a very coarse 30×30 mesh. The accuracy can be 10^{-3} for a moderately fine mesh.

The converged imaginary part is still 1% larger than the accurate value, and it is due to the finite computation region. We have also varied the computation region to 1.8a and 2.0a, and the systematic error has been found reduced at the cost of increased computation.

Figure 6 shows the calculated fundamental mode and 2^{nd}-order mode of the mode class 3 and 4 in the holey fiber. These modes are well confined by the single ring of the air holes and show the symmetries as discussed above. Once H_{x} and H_{y} (E_{x} and E_{y}) are solved, the other components can be obtained directly using Eqs. (6a–b).

#### 3.1 Effect of dispersive and lossy/gain materials

The effect of dispersive material on the group velocity dispersion can be easily obtained by any algorithm including this method that solves eigen-value problems for the mode propagation constants at a given wavelength. This is achieved by replacing the dielectric constant with a wavelength dependent one. In addition, the total dispersion can be estimated using the sum of material dispersion and waveguide dispersion.

The propagation loss/gain induced by lossy/gain media is also of interest in a waveguide since it is important for long haul transmission or laser applications. The compact-2D FDFD method is also capable of solving it, and no additional work is needed except that the real dielectric constant is replaced with a complex one whose imaginary part represents the loss/gain of the media; the calculated propagation constant would reflect the material loss. To separate the leakage loss from the material loss/gain, the PML layers are removed, leaving only the zero boundaries outside the computation region. The material loss in PCF has been analyzed in [50] using the Hermite-Gauss method and here we just verify their results for a standard fiber (a=2.2µm, cladding index=1.458, core index=1.475+ini and λ=1.55µm). The results are shown in Table 2. Our FDFD is in excellent agreement with the localized function methods.

#### 3.2 Spurious modes

One problem encountered in our method is the introduction of spurious modes. It is known that some spurious modes are from the lack of tangential continuity for the finite difference of the curved dielectric boundaries [28]. However, we did not observe such kind of spurious modes using our proper averaging technique. A major cause of the spurious modes in our method is the cladding modes due to the zero boundary or PML outside. Since the PML and zero boundary still reflect a very small part of energy back to the inside region, an artificial waveguide between the boundary and the air holes is formed. These modes are easily identified as spurious since they have much higher leakage loss. Figure 7 shows the mode patterns of some spurious cladding modes obtained in our calculation. These modes are confined between the boundary and air holes and have very small power portions in the core region.

## 4. Conclusions

In conclusion, we have used the compact-2D FDFD algorithm based on the direct solution of Maxwell’s equations to analyze leakage loss by introducing the anisotropic PML layers. It is found that the curved interface should be averaged according to Ampere’s Law in order to achieve a high accuracy and fast convergence in both the dispersion and leakage loss. Spurious modes are generated in the cladding area due to the artificial waveguide between the absorption boundary and the air holes.

## Acknowledgments

The research at Old Dominion University is supported by NASA Langley Research Center through NASA-University Photonics Education and Research Consortium (NUPERC). Some helpful discussions with Qian Jun at Old Dominion University are acknowledged.

## References and links

**1. **J. Broeng, “Photonic crystal fibers: a new class of optical waveguides,” Opt. Fiber Technol. **5**, 305–330 (1999). [CrossRef]

**2. **P. Russell, “Photonic crystal fibers,” Science **299**, 358–362 (2003). [CrossRef] [PubMed]

**3. **A. Bjarklev, J. Broeng, and A. S. Bjarklev, *Photonic crystal fibres* (Kluwer Academic Publishers, Boston/Dordrecht/London, 2003). [CrossRef]

**4. **Z. Zhu and T. G. Brown, “Analysis of the space filling modes of photonic crystal fibers,” Opt. Express **8**, 547–554 (2001). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-10-547 [CrossRef] [PubMed]

**5. **A. Ferrando, E. Silvestre, J. J. Miret, and P. Andres, “Vector description of higher-order modes in photonic crystal fibers,” J. Opt. Soc. Am. A **17**, 1333–1339 (2000). [CrossRef]

**6. **A. Ferrando, E. Silvestre, J. J. Miret, P. Andres, and M. V. Andres, “Full-vector analysis of a realistic photonic crystal fiber,” Opt. Lett. **24**, 276–278 (1999). [CrossRef]

**7. **J. Broeng, S. E. Barkou, T. Sondergaard, and A. Bjarklev, “Analysis of air-guiding photonic bandgap fibers,” Opt. Lett. **25**, 96–98 (2000). [CrossRef]

**8. **A. Ferrando, E. Silvestre, J. J. Miret, and P. Andres, “Nearly zero ultra-flattened dispersion in photonic crystal fibers,” Opt. Lett. **25**, 79–792 (2000). [CrossRef]

**9. **A. Weisshaar, J. Li, R. L. Gallawa, and I. C. Goyal, “Vector and quasi-vector solutions for optical waveguide modes using efficient Galerkin’s method with Hermite-Gauss basis functions,” J. Lightwave Technol. **13**, 1795–1780 (1995). [CrossRef]

**10. **W. Zhi, R. Guobing, L. Shuqin, and J. Shuisheng, “Supercell lattice method for photonic crystal fibers,” Opt. Express **11**, 980–991 (2003). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-9-980 [CrossRef] [PubMed]

**11. **T. M. Monro, D. J. Richardson, N. G. R. Broderick, and P. J. Bennett, “Holey optical fibers: an efficient modal model,” J Lightwave Technol. **17**, 1093–1102 (1999). [CrossRef]

**12. **D. Mogilevtsev, T. A. Birks, and P. St. J. Russel, “Group-velocity dispersion in photonic crystal fibers,” Opt. Lett. **23**, 1662–1664 (1998). [CrossRef]

**13. **A. Cucinotta, G. Peiosi, S. Selleri, L. Vincetti, and M. Zoboli, “Perfectly matched anisotropic layers for optical waveguide analysis through the finite-element beam-propagation method,” Microwave Opt. Techn. Lett. **23**, 67–69 (1999). [CrossRef]

**14. **K. Saitoh and M. Koshiba, “Full-vectorial finite element beam propagation method with perfectly matched layers for anisotropic optical waveguides,” J. Lightwave Technol. **19**, 405–413 (2001). [CrossRef]

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

**16. **M. Koshiba, Y. Tsuji, and M. Hikari, “Finite element beam propagation method with perfectly matched layer boundary conditions,” IEEE Trans. Magnetics **35**, 1482–1485 (1999). [CrossRef]

**17. **F. Brechet, J. Marcou, D. Pagnoux, and P. Roy, “Complete analysis of the characteristics of propagation into photonic crystal fibers by the finite-element method,” Opt. Fiber Technol. **6**, 181–191 (2000). [CrossRef]

**18. **C. Themistos, B. M. A. Rahman, A. Hadjicharalambous, and K. T. V. Grattan, “Loss/gain characterization of optical waveguides,” J Lightwave Technol. **13**, 1760–1765 (1995). [CrossRef]

**19. **K. Saitoh, M. Koshiba, T. Hasegawa, and E. Sasaoka, “Chromatic dispersion control in photonic crystal fibers: application to ultra-flattened dispersion,” Opt. Express **11**, 843–852 (2003). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-8-843 [CrossRef] [PubMed]

**20. **M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems,” J. Lightwave Technol. **18**, 737–743 (2000). [CrossRef]

**21. **S. Guenneau, S. Lasquellec, A. Nicolet, and F. Zolla, “Design of photonic crystal fibers using finite elements,” International J. Computation and Mathematics in Electrical & Electronics Engineering COMPEL **21**, 534–539 (2002). [CrossRef]

**22. **S. Guenneau, A. Nicolet, F. Zolla, and S. Lasquellec, “Numerical and theoretical study of photonic crystal fibers,” Progress in Electromagnetics Research **41**, 271–305 (2003).

**23. **D. H. Choi and W. J. R. Hoeffer, “The finite-difference time-domain method and its application to eigen-value problems,” IEEE Trans. Microwave Theory Tech. **34**, 1464–1470 (1986). [CrossRef]

**24. **A. Asi and L. Shafai, “Dispersion analysis of anisotropic inhomogeneous waveguides using compact 2D-FDTD,” Electron. Lett. **28**, 1451–1452 (1992). [CrossRef]

**25. **A. C. Cangellaris, “Numerical stability and numerical dispersion of a compact 2D FDTD method used for the dispersion analysis of waveguides,” IEEE Microwave Guided Wave Lett. **3**, 3–5 (1993). [CrossRef]

**26. **F. Zepparelli, P. Mezzanotte, F. Alimenti, L. Roselli, R. Sorrentino, G. Tartarini, and P. Bassi, “Rigorous analysis of 3D optical and optoelectronic devices by the Compact-2D-FDTD method,” Opt. and Quantum Electron. **31**, 827–841 (1999). [CrossRef]

**27. **S. Xiao, R. Vahldieck, and H. Jin, “Full-wave analysis of guided wave structures using a novel 2-D FDTD,” IEEE Microwave Guided Wave Lett. **2**, 165–167 (1992). [CrossRef]

**28. **N. Kaneda, B. Houshmand, and T. Itoh, “FDTD analysis of dielectric resonantors with curved surfaces,” IEEE Trans. Microwave Theory Tech. **45**, 1645–1649 (1997). [CrossRef]

**29. **K. Bierwirth, N. Schulz, and F. Arndt, “Finite-difference analysis of rectangular dielectric waveguide structures,” IEEE Trans. Microwave Theory Tech. **34**, 1104–1113 (1986). [CrossRef]

**30. **P. Lusse, P. Stuwe, J. Schule, 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]

**31. **Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructuredd optical fibers,” Opt. Express **10**, 853–864 (2002). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853 [CrossRef] [PubMed]

**32. **T. P. White, B. T. Kuhlmey, R. C. Mcphedran, D. Maystre, G. Renversez, C. M. de Sterke, and L. C. Botten, “Multipole method for microstructuredd optical fibers. I. Formulation,” J. Opt. Soc. Am. B **19**, 2322 (2002). [CrossRef]

**33. **D. Marcuse, “Solution of the vector wave equation for general dielectric waveguides by the Galerkin method,” IEEE J Quantum Electron. **28**, 459–465 (1992). [CrossRef]

**34. **S. Guo, F. Wu, K. Ikram, and S. Albin, “Analysis of circular fibers with an arbitrary index profile by the Galerkin method,” Opt. Lett. **29**, 32–34 (2004). [CrossRef] [PubMed]

**35. **S. Guo, S. Albin, and R. S. Rogowski, “Comparative analysis of Bragg fibers,” Opt. Express **12**, 198–207 (2004). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-1-207 [CrossRef] [PubMed]

**36. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. **14**, 302–307 (1966). [CrossRef]

**37. **S. Guo, F. Wu, S. Albin, and R. S Rogowski, “Photonic band gap analysis using finite-difference frequency-domain method,” Opt. Express **12**, 1741–1746 (2004). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-8-1741 [CrossRef] [PubMed]

**38. **C. P. Yu and H. C. Chang, “Compact finite-difference frequency-domain method for the analysis of two-dimensional photonic crystals,” Opt. Express **12**, 1397–1408 (2004). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-7-1397 [CrossRef] [PubMed]

**39. **H. Y. D. Yang, “Finite difference analysis of 2D photonic crystals,” IEEE Trans. Microwave Theory Technol. **44**, 2688–2695 (1996). [CrossRef]

**40. **N. A. Issa and L. Poladian, “Vector wave expansion method for leaky modes of microstructuredd optical fibers,” J. Lightwave Technol. **21**, 1005–1012 (2003). [CrossRef]

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

**42. **E. A. Marengo, C. M. Rappaport, and E. L. Miller, “Optimum PML ABC conductivity profile in FDFD,” IEEE Trans. Magnetics **35**, 1506–1509 (1999). [CrossRef]

**43. **F. L. Teixeira and W. C. Chew, “Systematic derivation of anisotropic PML absorbing media in cylindrical and spherical coordinates,” IEEE Microwave Guided Wave Lett. **7**, 371–373 (1997). [CrossRef]

**44. **F. L. Teixeira and W. C. Chew, “Unified analysis of perfectly matched layers using differential forms,” Microwave Opt. Technol. Lett. **20**, 124–126 (1999). [CrossRef]

**45. **T. Tischler and W. Heinrich, “Accuracy limitations of perfectly matched layers in 3D Finite-difference frequency domain method,” IEEE Microwave Theory Tech. **50**, 1885–1888 (2002).

**46. **U. Pekel and R. Mittra, “An application of the perfectly matched layer (PML) concept to the finite element method frequency domain analysis of scattering problems,” IEEE Microwave Guided Wave Lett. **5**, 258–260 (1995). [CrossRef]

**47. **W. C. Chew, J. M. Jin, and E. Michielssen, “Complex coordinate stretching as a generalized absorbing boundary condition,” Microwave & Opt. Technol. Lett. **7**, 363–369 (1997). [CrossRef]

**48. **P. R. McIsaac, “Symmetry-induced modal characteristics of uniform waveguides-II:Theory,” IEEE Trans. Microwave Theory Tech. **23**, 429–433 (1975). [CrossRef]

**49. **G. W. Milton, *The theory of composites* (Cambridge University Press, Cambridge, UK, 2002). [CrossRef]

**50. **R. Guobing, W. Zhi, L. Shuqin, and J. Shuisheng, “Full-vectorial analysis of complex refractive index photonic crystal fibers,” Opt. Express **12**, 1126–1135 (2004). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-6-1135 [CrossRef]