Open Access
Add to BibSonomyAbstract
The previously developed full-vectorial optical waveguide eigenmode solvers using pseudospectral frequency-domain (PSFD) formulations for optical waveguides with arbitrary step-index profile is further implemented with the uniaxial perfectly matched layer (UPML) absorption boundary conditions for treating leaky waveguides and calculating their complex modal effective indices. The role of the UPML reflection coefficient in achieving high-accuracy mode solution results is particularly investigated. A six-air-hole microstructured fiber is analyzed as an example to compare with published high-accuracy multipole method results for both the real and imaginary parts of the effective indices. It is shown that by setting the UPML reflection coefficient values as small as on the order of 10−40 ∼ 10−70, relative errors in the calculated complex effective indices can be as small as on the order of 10−12.
© 2011 Optical Society of America
1. Introduction
High-accuracy calculation of full-vectorial optical waveguide modes can be demanded in basic research of waveguide problems, such as understanding the corner effect of rectangular dielectric waveguides [1–3], or in the study and design of some special waveguide device structures, such as in precise determination of form birefringence of the fused fiber coupler [4]. As an alternative to such most often employed numerical methods as the finite difference method (FDM) [1, 2, 4] and the finite element method (FEM) [5], we have recently developed the full-vectorial pseudospectral frequency-domain mode solver (PSFD-MS or PSMS) [6, 7] which is based on multidomain pseudospectral methods, and demonstrated its extremely high numerical accuracy. For example, the analysis of the fundamental mode of the step-index optical fiber showed the relative error in the effective index can be as small as on the order of 10−12 [7]. The effective index (neff) is defined as the modal propagation constant divided by the free-space wavenumber.
The application of the pseudospectral method to solve electromagnetics problems is a relatively new approach compared to in the area of fluid dynamics. The method is attractive owing to its high-order accuracy and fast convergence behavior over traditional methods. It has been successfully employed to solve scattering problems in time domain [8–13]. As for frequency-domain problems, the pseudospectral frequency-domain (PSFD) method was formulated in [14] for solving the nonhomogeneous (nonzero-source) Helmholtz equation in a rectangualr-shape problem and in [6] and [7] into eigenvalue equations for eigenmode analysis of optical waveguides with arbitrary step-index profile. Full-vectorial formulations with the eigenvectors containing two transverse field components was developed in [6] and [7] for obtaining vector waveguide modes. A simpler scalar formulation has also been established for solving the Helmholtz equation in the analysis of two-dimensional (2-D) photonic-crystal band diagrams along with the required periodic boundary conditions [15]. In [6], [7], and [15], the multidomain Legendre or Chebyshev collocation method were employed. Spectral polynomials such as Legendre or Chebyshev polynomials were used to develop differential matrices for approximating differential operators in the Helmholtz equation over nonuniform Legendre-Gauss-Lobatto or Chebyshev-Gauss-Lobatto grid points and forming the matrix eigenvalue problem. We have found that the Chebyshev collocation method performs more efficiently than the Legendre one. The multidomain approach combined with the curvilinear mapping technique [16] can efficiently treat general multiple material structures having curved dielectric interfaces and properly fulfill field continuity conditions across the interface, i.e., suitable Dirichlet-type and Neumann-type boundary conditions (DBCs and NBCs), which is essential for achieving high-accuracy solutions along with the inherent high-order scheme of the pseudospectral method. Recently Huang et al. [17] have also employed the multidomain pseudospectral method to formulate a dielectric-waveguide mode solver but only for waveguides with rectangular-shaped cross-section.
We have further generalized the PSMS for treating waveguides with leakage properties, such as microstructured optical fibers [18, 19] and photonic crystal fibers [20] for which the mode power could leak into the cladding region and the mode index or neff would become complex, by applying the perfectly matched layer (PML) [21, 22] absorbing boundary conditions (ABCs) to the boundaries of the computational domain [23]. The detailed formulation will be given in this paper. The anisotropic or uniaxial PML (APML or UPML) proposed in [22] will be used and implemented into the ABC layer to absorb the outward leaky waves with the DBCs and NBCs rederived, forming the more general PSMS-UPML solver. In the implementation, the reflection coefficient and thickness of the UPML are important parameters for controlling the numerical accuracy. PMLs have been popularly employed in waveguide mode solvers using conventional numerical methods, e.g., the FDM [24] and the FEM [25]. One characteristic observed in these conventional methods was that the calculated neff would show oscillatory variation in the typically used PML reflection coefficient range of 10−7 ∼ 10−12. In this paper, we particularly investigate the role of the UPML reflection coefficient in achieving high-accuracy effective-index results. We will consider a six-air-hole microstructured fiber as an example, which was studied in [18] and [19] using a high-accuracy multipole method analysis, and compare our PSMS-UPML calculation with that of [18] and [19] for both the real and imaginary parts of the complex neff ’s. We will show that by setting the UPML reflection coefficient values as small as on the order of 10−40 ∼ 10−70, taking the advantage of the unprecedented accuracy of the PSMS, relative errors in the calculated complex neff ’s can be as small as on the order of 10−12. This issue regarding using small PML reflection coefficients was recently briefly discussed in [26] for one mode of the six-air-hole microstructured fiber, where PMLs based on stretched complex coordinates [27] were employed. In this paper, detailed examinations will be given with different PML parameters and for more modes.
The rest of this paper is outlined as follows. The physical equations and the required DBCs and NBCs in an anisotropic medium are given in Section 2. The PSMS-UPML formulation is derived in Section 3. Numerical analysis of a six-air-hole microstructured fiber is discussed in Section 4 along with the investigation of the role of the UPML reflection coefficient. The conclusion is drawn in Section 5.
2. The H-formulation equations in an anisotropic Medium
We first consider the propagation of the time-harmonic electromagnetic wave in an inhomogeneous anisotropic waveguide uniform in the z direction at the angular frequency ω. The source-free time-harmonic Maxwell’s equations are written as
where ɛ0 and μ0 are permittivity and permeability of free space, respectively, and the relative permittivity tensor [ɛ] and the relative permeability tensor [μ] are defined respectively as By taking the curl of Eq. (4) and using Eq. (3), we obtain the vectorial wave equation where k0 is the wavenumber in free space and . For guided waves propagating in the z direction, H⃗ can be expressed as where β is the modal propagation constant. From Eq. (2), we have By substituting Eq. (6) into Eq. (5) and using Eq. (7), we obtain where the differential operators in the matrix are defined asNow considering an interface of arbitrary shape between two different material regions, a and b, as shown in Fig. 1. With the x-y coordinate system defined in the figure, the available boundary conditions (DBCs and NBCs) are, respectively, expressed as
and where nx and ny are the x- and y-components of the normal unit vector at the interface in Fig. 1 and the superscripts a and b refer to regions a and b, respectively. Of course, all quantities in Eqs. (13)–(15) are evaluated at the interface.
Fig. 1 Interface between two homogeneous regions, Regions a and b, with relative permittivity and permeability tensors, ([ɛa], [μa]) and ([ɛb], [μb]), in the waveguide cross-section. n̂ is a unit vector normal to the interface.
3. The PSMS-UPML Formulation
We discuss the required conditions and materials for the UPML. Assume the material in region a has the relative permeability and permittivity tensors:
where μx, μy, μz, ɛx, ɛy, and ɛz are constants, and that in region b (considered as the UPML region) has the corresponding tensors as [μ]UPML and [ɛ]UPML.In order to impedance match region b to region a such that no outward going wave is reflected at the interface, we select [μ]UPML and [ɛ]UPML such that [28]
Consequently, [μ] and [ɛ] in the whole computational domain including UPML regions can be expressed as with where and are assigned appropriate values to control the field attenuation variables along the x and y directions, respectively, in UPML regions. Assuming the appropriate values of theUPML medium under impedance matching condition have m-power profiles as
where dx and dy are the thicknesses of the UPML as shown in Fig. 2 which shows the cross-section of an arbitrary leaky waveguide with the computational domain surrounded by UPML regions, and ρx and ρy are the distances from the beginning of the UPML along the x and y directions, respectively. Using the theoretical reflection coefficient R [21] (note that the optimal value of R approaches to zero) at the inner interface of the UPML, the maximum values of αx and αy can be determined as where λ is the wavelength corresponding to ω and with min(ɛi) being the smallest of ɛx, ɛy, and ɛz of the anisotropic material region interfaced to the UPML. Note that is a good choice to guarantee that the values of both αmax(x) and αmax(y) are large enough corresponding to R, as for the isotropic material region case [21], for unknown polarized directions in the anisotropic material. With m = 2, sx and sy in Eqs. (18) and (19) become By substituting Eqs. (16) and (17) into Eqs. (9)–(12), we obtain where s′x = ∂(1/sx)/∂x, s′y = ∂(1/sy)/∂y, s″x = ∂2(1/sx)/∂x2, and s″y = ∂2(1/sy)/∂y2.Utilizing the pseudospectral Legendre or Chebyshev method, as formulated in [7], Eq. (8) will lead to a matrix eigenvalue equation of the form, A̿H̄ = β2H̄, with , where the superscript T denotes transpose, H̄x and H̄y are vectors composed of Hx and Hy values, respectively, at Legendre or Chebyshev collocation points, and A̿ is the operator matrix. The eigenvalues can be solved by applying the shift inverse power method (SIPM). The computational domain will be divided into several subdomains of curvilinear quadrilateral shape. Each subdomain will be mapped into a square region [−1, 1] × [−1, 1] in the general curvilinear coordinates (ξ, η). Each subdomain, including the four interfaces, would contain (M + 1) × (N + 1) grid points if Legendre or Chebyshev polynomial of degree M is employed in the ξ direction and that of degree N is employed in the η direction. With both Hx and Hy values to be determined at each grid point, there are 2(M + 1)(N + 1) unknowns in a single subdomain. So for a single-subdomain situation, the matrix A̿ would be a 2k × 2k one, with k = (M + 1)(N + 1), as given in Eq. (26) of [7]. With Eqs. (26)–(29), the A̿ matrix in Eq. (26) of [7] is modified to AUPML in the UPML regions of the computational domain as
where with and sxjj, syjj (j = 1, ⋯ ,k), , , D̿x00, and D̿y00 as defined in Eq. (26) of [7], and I̿ is a k × k identity matrix. And the related DBCs and NBCs can be written, respectively, asIn the general multidomain situation, if the computational domain is divided into p subdomains, the matrix eigenvalue equation before imposing boundary conditions across the subdomain interfaces would be similar to that of Eq. (27) of [7], with the operator matrix being a 2pk × 2pk one. The adjacent subdomains are connected by imposing DBCs and NBCs on the matrix elements, corresponding to the interface grids of the subdomains, in this 2pk × 2pk matrix eigenvalue equation. How such modification of the matrix equation is conducted can be referred to the Appendix in [15]. When we implement the UPML into an arbitrary leaky waveguide problem of Fig. 2, 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. The sx and sy values for different regions are listed in Table 1.
Fig. 2 Cross-section of an arbitrary leaky waveguide problem with the computational domain surrounded by UPML regions.
Table 1. Definition of sx and sy values in the UPML and non-UPML regions.
4. Numerical example: analysis of a six-air-hole microstructured fiber
As a numerical example, we consider a six-air-hole microstructured fiber, or a triangular holey fiber with one ring of six air holes symmetrically arranged around the core, to assess our PSMS-UPML. This structure has been analyzed with high numerical accuracy using the multipole method in [18, 19]. It should be noted that with the multipole method, the complex neff ’s of the modes would be obtained by searching zeros of a complex function resulting from the determinant of a large complex matrix corresponding to a derived homogeneous system of equations. Achievement of the solutions of high accuracy requires special techniques and careful procedures [18, 19]. Our formulation, however, gives a standard eigenvalue matrix equation and the eigen solutions are easy to obtain. The first ten modes for the six-air-hole fiber were solved in [18]. Then, in [19], after discussing in more detail numerical aspects of the formulation, guaranteeing accurate results, and doing numerical verification, more accurate results for the (nondegenerate) sixth (p = 1) mode were provided, with the imaginary part of neff pushed down to 10−11. We thus compare first our analysis of this sixth mode with that of [19].
Owing to the structure symmetry, only a quarter of the fiber cross-section with three UPML regions is included in the computational domain, as shown in Fig. 3(a). As in [18, 19], we assume the hole diameter d = 5μm, the hole pitch Λ = 6.75μm, and the silica background with index n = 1.45. Figure 3(b) shows the domain and mesh division profile including the UPML regions employing the Chebyshev-Gauss-Lobatto grid points, where 32 subdomains are adopted and the mesh pattern of each subdomain depicted is for polynomials of degree M = N = 8, corresponding to (8 + 1)2 × 2 × 32 = 5184 total unknowns in the eigenvector for the transverse magnetic field components. According to [18], the accurate neff for the sixth mode (with class p = 1 as designated in [18]) at wavelength 1.45μm would be 1.438364934 – j1.41647 ×10−6 and the real and imaginary parts converge to 10−9 and 10−11, respectively. In our calculations, we obtain the corresponding neff values as 1.43836493417887 – j1.41647611 × 10−6 when M = N = 20 is used and 1.43836493417887 – j1.41647598 × 10−6 when M = N = 22 is used, which agree with that of [19] very well. In the following, we will use neff, ref = 1.438364934178 – j1.416476 ×10−6 (with 10−12 precision in both the real and the imaginary parts) as a reference to discuss our calculation results using the PSMS-UPML under different UPML parameters. In Fig. 3(a), the computational window sizes are Wx = Wy = 15.75μm, the UMPL thicknesses are dx = dy = 2μm, and the operating wavelength is λ = 1.45μm. Figures 4(a) and (b), respectively, show the relative errors (RErs) in the real and the imagniary parts of neff of the sixth mode relative to the above-mentioned neff, ref with respect to the number of unknowns for different reflection coefficient values (R = 10−6, 10−7, 10−9, 10−10, 10−15, 10−17, 10−20, 10−40, 10−50, 10−60, and 10−70). Note that the reference value is obtained with R = 10−70. It is seen that in both figures the RErs with R ≤ 10−40 can drop to as small as on the order of 10−10 but not for R > 10−40. Note that in Fig. 4(b), the REr of Im[neff] for R = 10−7 appears to be smaller than those using 10−9 ≥ R ≥ 10−20. However, Fig. 4(a) shows that the REr of Re[neff] for R = 10−7 is only better than that for R = 10−6. We thus examine the total relative error of the complex neff, expressed as TREr(neff) and defined as TREr(neff) = |neff – neff, ref|. The TREr(neff)’s versus the number of unknowns are shown in Fig. 5(a), where it is seen that the TREr(neff) for large number of unknowns is better with the smaller R, i.e., TREr(neff)R=10−6 > TREr(neff)R=10−7 > TREr(neff)R=10−9 > TREr(neff)R=10−10 > ⋯ > TREr(neff)R=10−70. However, in Fig. 5(a), the dependence on R of the convergence behavior does not follow this trend when the number of unknowns is smaller than 6400 for TREr(neff) < 10−7. This fact could lead to the wrong conclusion that the optimal R values are within 10−6 ∼ 10−12 when analysis methods of limited accuracy are employed since they could not achieve sufficiently low TREr values to reach better numerical convergency. The maximum number of unknowns used in Figs. 4 and 5(a) is 28224, corresponding to M = N = 20. For this enough high polynomial degree, the TREr(neff) versus R is plotted in Fig. 5(b), which shows the decrease of TREr(neff) with decreasing R until R = 10−40, in consistent with the PML principle that R → 0 for the ideal PML.
Fig. 3 (a) The cross-section of the triangular holey fiber with one ring of six air holes. The computational domain with UPML regions I, II, and III is shown to contain a quarter of the cross-section. (b) Mesh division profile of the computational domain.
Fig. 4 (a) Relative errors in the real part of the effective index for the sixth mode in the holey fiber of Fig. 3 versus the number of unknowns for different UPML reflection coefficient values. (b) Relative errors in the imaginary part of the effective index for the sixth mode in the holey fiber of Fig. 3 versus the number of unknowns for different UPML reflection coefficient values.
Fig. 5 (a) Total relative errors in the effective index for the sixth mode in the holey fiber of Fig. 3 versus the number of unknowns for different UPML reflection coefficient values. (b) Total relative errors in the effective index for the sixth mode in the holey fiber of Fig. 3 versus the UPML reflection coefficient using 28224 unknowns.
The calculated Re[neff]’s and Im[neff]’s corresponding to the R = 10−70 results in Fig. 5(a) are listed in Table 2 for different polynomial degrees and numbers of unknowns. In Table 2 and the following Tables, the digits of each effective index are shown down to the place of 10−14. It can be seen that our results with M = N = 18 have already agreed with the obtained value of [19], i.e., the agreement is good to 10−9 for Re[neff] and to 10−11 for Im[neff]. The calculated Re[neff]’s and Im[neff]’s with M = N = 20 using different Rs from 10−6 to 10−70, corresponding to Fig. 5(b), are listed in Table 3. We also examine the effect of the distance between the six air-holes and the UPML layers, or the size of the computational window. With fixed dx = dy = 2μm, Wx = Wy in Fig. 3 is varied from 13 to 15.75μm. The number of unknowns is taken to be 28224 and M = N = 20. The obtained (Re[neff] – 1.438364934) and Im[neff] are shown in Figs. 6(a) and (b), respectively, versus Wx = Wy for R = 10−7 and 10−70. It is seen that both Re[neff] and Im[neff ] for R = 10−70 almost do not vary with Wx = Wy under the scales shown while those for R = 10−7 show significant variations. The ranges of variation for both Re[neff] and Im[neff] values for R = 10−7 are almost on the same order and are about 2 × 10−7. And it is interesting to observe that for R = 10−7, when Re[neff] reaches the most accurate value (near zero in Fig. 6(a)), Im[neff] in Fig. 6(b) would be with the largest (positive or negative) error, and vice versa. Finally, using R = 10−70, TREr(neff)’s versus the number of unknowns for Wx = Wy = 13.75μm and 15.75μm are plotted in Fig. 7. It is seen that for both window sizes, the TREr(neff) reaches the lowest value when the number of unknowns increases to 14400. No better accuracy can be revealed for more unknowns due to reaching the accuracy limit of the reference value. The results with Wx = Wy = 13.75μm are seen to be a little better than those with Wx = Wy = 15.75μm because smaller window size gives better mesh resolution for the same number of unknowns. According to the above discussions, we can conclude that the reflection coefficient R is the key UPML parameter for achieving high-accuracy neff and its value of as small as 10−70 is suggested.
Fig. 6 (a) Real part and (b) imaginary part of the effective index of the sixth mode of the holey fiber of Fig. 3 versus the width (Wx = Wy) of the computational domain for R =10−7 and 10−70.
Fig. 7 Total relative errors in the effective index of the sixth mode in the holey fiber of Fig. 3 versus the number of unknowns for two different computational domain sizes, Wx = Wy = 13.75 μm and Wx = Wy = 15.75 μm.
Table 2. Real and imaginary parts of the effective index of the sixth mode of the holey fiber of Fig. 3 obtained using R = 10−70 and different numbers of unknowns.
Table 3. Real and imaginary parts of the effective index of the sixth mode of the holey fiber of Fig. 3 obtained using 28224 unknowns and different values for the UPML reflection coefficient.
In the following, we analyze other modes on the same six-air-hole microstructured fiber by taking Wx = Wy = 13.75μm. Only the first five modes are considered. The seventh to the tenth modes studied in [18] are not of much interest due to their much higher leaky losses. For the fundamental (class p = 3) mode, the calculated Re[neff]’s and Im[neff]’s with M = N = 20 using R = 10−50, 10−60, and 10−70, respectively, are listed in Table 4. Then, using R = 10−70, the calculated Re[neff]’s and Im[neff]’s are listed in Table 5 for different polynomial degrees and numbers of unknowns. We can safely conclude that the converged neff is 1.445395256948 – j3.1947 × 10−8. We then list in Table 6 the calculated Re[neff]’s and Im[neff]’s, with the corresponding loss values in dB/m, of the first six modes which are the p = 3, 4, 2, 5, 6, and 1 modes in sequence designated in [18] using M = N = 20 and R = 10−70. Note that p = 3 and p = 4 modes are twofold degenerate, and so are p = 5 and p = 6 modes. The loss in dB/m is obtained by multiplying Im[neff] by (20/[ln(10)])(2π/λ) × 106 [18], where λ is in micrometers. Although the numerical accuracy of the results in [18] is not as high as those in [19], the loss values in dB/m for the first six modes listed in Table 1 of [18], i.e., 1.2, 20, 37, and 52, are fairly close to the corresponding ones in Table 6 since the difference in Im[neff] between [18] and this work is at most 6.5%.
Table 4. Real and imaginary parts of the effective index of the fundamental mode of the holey fiber of Fig. 3 obtained using 28224 unknowns and three different values for the UPML reflection coefficient.
Table 5. Real and imaginary parts of the effective index of the fundamental mode of the holey fiber of Fig. 3 obtained with R = 10−70 and different numbers of unknowns.
Table 6. Real and imaginary parts of the effective indices and losses of the first six modes of the holey fiber of Fig. 3 obtained using 28224 unknowns and R = 10−70.
Finally, normalized field profiles for these six modes are presented in Fig. 8. In [18], normalized fields |Ez| and |ZHz| were plotted for p = 3, p = 4, and p = 2 modes, where Z denotes the impedance of free space. In this PSMS-UPML, which is based on the transverse-magnetic-field formulation, Hx and Hy are first obtained. Then, Ez and Hz can be calculated with high accuracy, although involving the derivatives of Hx and Hy, since the derivative is simply performed by the differential matrix operation under the pseudospectral theory. In Fig. 8, for each mode, the |Ez|, |Hz|, |Ex|, and |Hy| mode-field profiles are shown, with the maximum of each profile set to be unity. The |Ez| and |Hz| profiles for the first three modes (p = 3, p = 4, and p = 2) are found to be consistent with those shown in Figs. 4 and 5 of [18]. One interesting observation is worth mentioning in the following. Note that |Hx| and |Hy| profiles of the p = 2 and p = 6 modes appear to have very little differences but their |Ez| and |Hz| profiles have obvious different features. The same can be observed between the p = 5 and p = 1 modes. Moreover, it can be seen that the |Hx| (|Hy|) profile of the p = 2 mode resembes the |Hy| (|Hx|) profile of the p = 5 mode, and the same feature exists between the p = 6 and p = 1 modes. Such fact should not be surprising since these four modes actually possess very close Re[neff]’s with the same leading digits, 1.438, as seen in Table 6. And more noticeable differences appear in the |Ez| and |Hz| profiles.
Fig. 8 Normalized |Ex|, |Ez|, |Hx|, and |Hy| field profiles of the first six modes of the holey fiber of Fig. 3.
5. Conclusion
We have reported the implementation of UPML absorbing boundary conditions into our previously developed full-vectorial optical waveguide eigenmode solver using pseudospectral frequency-domain (PSFD) formulations for waveguides with arbitrary step-index profile. With UPMLs around the computational domain in this PSMS-UPML solver, leaky waveguides can be treated and their complex modal effective indices determined. The solver is derived under the Hx – Hy formulation by considering the Helmholtz equations. The spatial derivatives in the equations are approximated at collocation points by utilizing Chebyshev-Lagrange interpolating polynomials, leading to a matrix eigenvalue equation with the squared complex propagation constant as the eigenvalue. We have shown other field components can be determined with high accuracy from the obtained Hx and Hy (eigen vector) with the required spatial derivatives treated by the same approximation at collocation points. The high-accuracy performance of the PSMS-UPML is demonstrated by using a six-air-hole microstructured fiber as an analysis example, with the comparison made with published high-accuracy multipole method results [19] for both the real and imaginary parts of the effective indices. In particular, the role of the UPML reflection coefficient, R, in achieving high-accuracy mode solution results is investigated in detail. It is shown that by setting the values of R as small as on the order of 10−40 ∼ 10−70, relative errors in the calculated complex effective indices can be as small as on the order of 10−12.
Acknowledgments
This work was supported in part by the National Science Council of the Republic of China under grants NSC97-2221-E-002-043-MY2 and NSC97-2221-E-151-055-MY2, in part by the Excellent Research Projects of National Taiwan University under grant 98R0062-07, and in part by the Ministry of Education of the Republic of China under “The Aim of Top University Plan” grant. The authors would like to acknowledge the National Center for High-Performance Computing in Hsinchu, Taiwan for providing useful computing resources.
References and links
1. 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]
2. G. R. Hadley, “High-accuracy finite-difference equations for dielectric waveguide analysis II: Dielectric corners,” J. Lightwave Technol. 20, 1219–1231 (2002). [CrossRef]
3. N. Thomas, P. Sewell, and T. M. Benson, “A new full-vectorial higher order finite-difference scheme for the modal analysis of rectangular dielectric waveguides,” J. Lightwave Technol. 25, 2563–2570 (2002). [CrossRef]
4. 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]
5. 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]
6. P. J. Chiang, C. S. Yang, C. L. Wu, C. H. Teng, and H. C. Chang, “Application of pseudospectral methods to optical waveguide mode solvers,” OSA 2005 Integrated Photonics Research and Applications (IPRA ’05) Technical Digest (Optical Society of America, 2005), paper IMG4.
7. P. J. Chiang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44, 56–66 (2008). [CrossRef]
8. B. Yang, D. Gottlieb, and J. S. Hesthaven, “Spectral simulations of electromagnetic wave scattering,” J. Comput. Phys. 134, 216–230 (1997). [CrossRef]
9. B. Yang and J. S. Hesthaven, “A pseudospectral method for time-domain computation of electromagnetic scattering by bodies of revolution,” IEEE Trans. Antennas Propagat. 47, 132–141 (1999). [CrossRef]
10. J. S. Hesthaven, P. G. Dinesen, and J. P. Lynov, “Spectral collocation time-domain modeling of diffractive optical elements,” J. Comput. Phys. 155, 287–306 (1999). [CrossRef]
11. G. Zhao and Q. H. Liu, “The 3-D multidomain pseudospectral time-domain algorithm for inhomogeneous conductive media,” IEEE Trans. Antennas Propagat. 52, 742–749 (2004). [CrossRef]
12. C. H. Teng, B. Y. Lin, H. C. Chang, H. C. Hsu, C. N. Lin, and K. A. Feng, “A Legendre pseudospectral penalty scheme for solving time-domain Maxwell’s equations,” J. Sci. Comput. 36, 351–390 (2008). [CrossRef]
13. B. Y. Lin, H. C. Hsu, C. H. Teng, H. C. Chang, J. K. Wang, and Y. L. Wang, “Unraveling near-field origin of electromagnetic waves scattered from silver nanorod arrays using pseudo-spectral time-domain calculation,” Opt. Express17, 14211–14228 (2009). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-16-14211. [CrossRef] [PubMed]
14. Q. H. Liu, “A pseudospectral frequency-domain (PSFD) method for computational electromagnetics,” IEEE Antennas Wireless Propagat. Lett. 1, 131–134 (2002). [CrossRef]
15. P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E 75, 026703 (2007). [CrossRef]
16. W. J. Gordon and C. A. Hall, “Transfinite element methods: blending-function interpolation over arbitrary curved element domains,” Numer. Math. 21, 109–129 (1973). [CrossRef]
17. C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron. 11, 457–465 (2005). [CrossRef]
18. T. P. White, B. T. Kuhlmey, R. C. Mcphedran, D. Maystre, G. Renversez, C. M. de Sterke, and L. C. Botten, “Multi-pole method for microstructured optical fibers. I. Formulation,” J. Opt. Soc. Am. B 19, 2322–2330 (2002). [CrossRef]
19. B. T. Kuhlmey, T. P. White, G. Renversez, D. Maystre, L. C. Botten, C. M. de Sterke, and R. C. Mcphedran, “Multi-pole method for microstructured optical fibers. II. Implementation and results,” J. Opt. Soc. Am. B 19, 2331–2340 (2002). [CrossRef]
20. P. Russell, “Photonic crystal fibers,” Science 299, 358–362 (2003). [CrossRef] [PubMed]
21. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comp. Phys. 114, 185–200 (1994). [CrossRef]
22. Z. S. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Antennas Propagat. 43, 1460–1463 (1995). [CrossRef]
23. P. J. Chiang and H. C. Chang, “Analysis of leaky optical waveguides using pseudospectral methods,” OSA 2006 Integrated Photonics Research and Applications (IPRA ’06) Technical Digest (Optical Society of America, 2005), paper ITuA3.
24. C. P. Yu and H. C. Chang, “Yee-mesh-based finite difference eigenmode solver with PML absorbing boundary conditions for optical waveguides and photonic crystal fibers,” Opt. Express12, 6165–6177 (2004). http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-25-6165. [CrossRef] [PubMed]
25. Y. Tsuji and M. Koshiba, “Guided-mode and leaky-mode analysis by imaginary distance beam propagation method based on finite element scheme,” J. Lightwave Technol. 18, 618–623 (2000). [CrossRef]
26. P. J. Chiang and Y. C. Chiang, “Pseudospectral frequency-domain formulae based on modified perfectly matched layers for calculating both guided and leaky modes,” IEEE Photon. Technol. Lett. 22, 908–910 (2010). [CrossRef]
27. 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]
28. 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]
