## Abstract

We propose an accurate and computationally efficient rational Chebyshev multi-domain pseudo-spectral method (RC-MDPSM) for modal analysis of optical waveguides. For the first time, we introduce rational Chebyshev basis functions to efficiently handle semi-infinite computational subdomains. In addition, the efficiency of these basis functions is enhanced by employing an optimized algebraic map; thus, eliminating the use of PML-like absorbing boundary conditions. For leaky modes, we derived a leaky modes boundary condition at the guide-substrate interface providing an efficient technique to accurately model leaky modes with very small refractive index imaginary part. The efficiency and numerical precision of our technique are demonstrated through the analysis of high-index contrast dielectric and plasmonic waveguides, and the highly-leaky ARROW structure; where finding ARROW leaky modes using our technique clearly reflects its robustness.

© 2016 Optical Society of America

## 1. Introduction

Optical waveguides are the fundamental building blocks around which many photonic devices are developed. Prior to the fabrication of these devices, thorough understanding of the light propagation in optical waveguides is essential. Among the methods which provide analytical or semi-analytical solutions are the mode-matching method [1, 2] and the eigenmode expansion method [3]. However, the most commonly used modal solution techniques are based on the finite difference method [4, 5] and the finite element method [6, 7]. Although successful in handling many problems, they produce matrices of huge size and have a relatively slow convergence rate of order *𝒪*(1/*N ^{p}*), where

*N*is the number grid points and

*p*is a positive constant. On the other hand, multi-domain pseudo-spectral methods (MDPSMs) [8, 9, 10] emerged recently and proved to be very accurate and efficient with a significant reduction in memory usage and computation time in addition to the fast convergence rate. The MDPSM handles semi-infinite extensions of computational domains through different approaches. One technique was proposed in [11] using Laguerre-Gauss basis functions to mimic the semi-infinite extension of the computational domain. Although this approach was successful, it required prior estimation of a very important scaling factor which is related to the decay rate of the field in outer subdomains, and significantly affects the convergence of the method. Another alternative approach was proposed in [12] by using Chebyshev functions in all subdomains while truncating the computational domain. Moreover, perfectly-matched layers (PMLs) are also commonly used to truncate the computational domain. However, there are still serious simulation problems in which PMLs collapse, confront inevitable reflections or even exponential growth [13]. Other major drawbacks are the high sensitivity to the selection of the PML parameters such as conductivity profile, thickness of the layer, permittivity, and the distance to the structure. Moreover, non-physical Berenger PML modes appear in the numerical solution that cannot be wiped out [14].

In this paper, we introduce an efficient numerical technique that can alleviate the problems associated with the existing treatments of semi-infinite computation domains in the context of MDPSMs. Our approach is based on using rational Chebyshev functions in semi-infinite subdomains. The rational basis functions are constructed by composing cardinal Chebyshev functions with a carefully selected conformal map. The proposed approach renders itself to have advantages summarized as:

- Unlike Laguerre functions [11], rational Chebyshev functions are bounded i.e. no exponential function with unknown scaling is required to ensure convergence. Hence, convergence is independent of any prior estimates.
- Due to the properties of rational Chebyshev functions, we don’t need to employ Perfectly-Matched Layers (PMLs) or PML-like absorbing boundary conditions.
- In addition to efficient handling of guided modes, the proposed scheme is very efficient in handling leaky modes. This is accomplished by deriving a nonlinear radiation (Robin-like) boundary condition at the guide-substrate interface, which is linearized using a Taylor series expansion. Through this treatment, exact analytical representation of leaky modes is newly proposed.

## 2. Rational Chebyshev multi-domain pseudo-spectral method

The wave equation for 2D optical waveguides is

*n*

^{2}(

*y*), the whole domain is first divided into a finite number of subdomains at interfaces of discontinuity denoted by

*d*,

_{i}*i*= 1, 2,...,

*M*for a structure having

*M*interfaces of different materials as shown in Fig. 1. This leads to a computational domain with (

*M*+ 1) subdomains, Ω

*,*

_{i}*i*= 1, 2,...,

*M*+ 1. The interface boundary conditions at the points of discontinuity,

*d*are given by

_{i}*d*, respectively.

_{i}#### 2.1. Rational Chebyshev functions with algebraic maps

The proposed method is based on the cardinal Chebyshev function expansion [11, 15] defined for *ỹ* ∈ [−1, 1] as follows

*T*(

_{N}*ỹ*) = cos(

*N*cos

^{−1}(

*ỹ*)) is the Chebyshev polynomial of degree

*N*,

*σ*

_{0}=

*σ*= 2, and

_{N}*σ*= 1 for

_{j}*j*= 1,...,

*N*− 1. The basis functions in Eq. (5) are defined based on the Gauss-Lobatto collocation points

*ỹ*= −cos(

_{k}*kπ/N*),

*k*= 0, 1,..,

*N*. For simplicity of reference, we define the subdomain Ω

_{1}= (−∞,

*d*

_{1}), Ω

*= (*

_{i}*d*

_{i−1},

*d*),

_{i}*i*= 2,...,

*M*, and Ω

_{M+1}= (

*d*, ∞).

_{M}Since the cardinal Chebyshev basis functions are defined on the interval [−1, 1], conformal map are needed in order to define suitable basis functions in each subdomain, Ω* _{i}* shown in Fig. 1.

Referring to Fig. 1, we have two types of subdomains. First, bounded inner subdomains, Ω_{2},...,Ω* _{M}* and two semi-infinite subdomains Ω

_{1}and Ω

_{M+1}.

- For the bounded subdomains, the maps are simply linear and take the form and their mapped Gauss-Lobatto collocation points are given by
- For the semi-infinite subdomains Ω
_{1}≡ (−∞,*d*_{1}) and Ω_{M+1}≡ (*d*, ∞), we define the algebraic maps_{M}*ϕ*_{1}(*y*;*N*_{1}) and*ϕ*_{M+1}(*y*;*N*_{M+1}) in Ω_{1}and Ω_{M+1}, respectively. In Ω_{1}the map*ϕ*_{1}(*y*;*N*_{1}) is given by [16] where*τ*(*N*_{1}) =*κ*_{1}*N*_{1}is an linear function with*κ*_{1}> 0 and*N*_{1}+ 1 is the number of basis functions used in Ω_{1}. The new collocation points, ${y}_{j}^{(1)}$ in the subdomain Ω_{1}are given by$${y}_{j}^{(1)}={d}_{1}-\tau ({N}_{1})\frac{1+\text{cos}(j\pi /{N}_{1})}{1-\text{cos}(j\pi /{N}_{1})},$$at which we obtain the physical values of the electric or magnetic field components. Similarly, for the semi-infinite subdomain Ω_{M+1}, we employ the map*ϕ*_{M+1}(*y*;*N*_{M+1}) given by$${\varphi}_{M+1}(y;{N}_{M+1})=\frac{(y-{d}_{m})-\alpha ({N}_{M+1})}{(y-{d}_{m})+\alpha ({N}_{M+1})},$$where*α*(*N*_{M+1}) =*κ*_{2}*N*_{M+1},*κ*_{2}> 0. The corresponding collocation points, ${y}_{j}^{(M+1)}$ in the subdomain Ω_{M+1}are then given by

The solution of the wave Eq. (1) is obtained by expanding the field in terms of mapped basis functions obtained by composing Chebyshev functions by the suitable conformal map for each subdomain Ω* _{i}*, as we have described and imposing Dirichlet zero-boundary conditions at

*y*= ±∞ as shown in Fig. (2).

The field expansion in a subdomain Ω* _{i}* takes the form

*N*+ 1 denotes the number of basis functions in a subdomain Ω

_{i}*and the circle “○” in Eq. (12) denotes a compositions of the functions*

_{i}*C*(

_{j}*ỹ*) with the map

*ϕ*(

_{i}*y*), namely (

*C*○

_{j}*ϕ*)(

_{i}*y*) =

*C*(

_{j}*ϕ*(

_{i}*y*)). The maps

*ϕ*(

_{i}*y*) all satisfy

*ϕ*(Ω

_{i}*) = [−1, 1] i.e. the function*

_{i}*ϕ*maps the subdomain Ω

_{i}*onto the interval [−1, 1]. The unknown coefficients ${u}_{j}^{(i)}$ in Eq. (12) represent the field values in the subdomain Ω*

_{i}*at the mapped Gauss-Lobatto collocation points for that subdomain, given by ${y}_{j}^{(i)}={\varphi}_{i}^{-1}({\tilde{y}}_{j})$. Cardinality of the basis functions in Eq. (12) is preserved i.e. we have $\left({C}_{j}\u25cb{\varphi}_{i}\right)\left({y}_{k}^{(i)}\right)={\delta}_{kj}$, where*

_{i}*δ*is the Kronecker delta.

_{kj}It is worth noting that composing the algebraic maps *ϕ*_{1}(*y*) and *ϕ*_{M+1}(*y*) defined by Eqs. (8) and (10), respectively, leads to rational Chebyshev basis functions defined on the semi-infinite subdomains Ω_{1} and Ω_{M+1}. The proposed treatment accounts for semi-infinite extensions of the outer subdomains, hence non-physical modes have no chance to appear in the solution spectrum. It also leads to an efficient replacement of the non-physical perfectly-matched layers (PMLs).

The multi-domain pseudo-spectral method reported [12] applies a domain truncation approach. In that approach, the outer semi-infinite subdomains Ω_{1} and Ω_{M+1} shown in Fig. 1 are truncated by replacing the original computational domain (−∞, ∞) by a computational window [−*W*, *W*] for *W* > max{|*d*_{1}|, |*d _{M}*|} and imposing Dirichlet zero-boundary conditions at

*y*= ±

*W*. Truncation of the domain assumes rapid decay rate of the field in outer subdomains. Although this approach is successful, the appropriate value of

*W*should vary according to the refractive index profile (RIP) to accurately describe the field in outer regions. Further, the value of

*W*should increase as the number of basis functions in the outermost subdomains increases to maintain high convergence rate (see [16] and references therein).

Another approach to handle the subdomains Ω_{1} and Ω_{M+1} was proposed in [11]. It was based on using Laguerre-Gauss basis functions. However, that approach requires a priori estimate for a scaling parameter whose value significantly affects both the accuracy and rate of convergence.

#### 2.2. Leaky modes treatment

For multi-layer planar waveguides where the substrate has a high refractive index, leaky modes arise and accurate calculation of the effective index imaginary part is challenging. Instead of expanding the field in the outer subdomain representing the substrate with high refractive index, we use the analytic expression of the field *u*(*y*) in the substrate region. From this expression we derive the suitable boundary conditions at the guide-substrate interface. This approach can be described follows

For a substrate with refractive index *n _{s}* in the subdomain Ω

_{M+1}, the field

*u*(

*y*) denoted by

*u*

^{(M+1)}can be expressed as

*y*and evaluating at $y={d}_{M}^{-}$, yield

_{1}is similar and yields the following condition

*p*> 1

#### 2.3. Formation of the linear eigenvalue problem

A linear eigenvalue problem for modal solution is deduced by substituting Eq. (12) into the wave equation Eq. (1), and demanding that Eq. (1) is satisfied in each subdomain Ω* _{i}* at the corresponding set of collocation points
${y}_{j}^{(i)}$ except for the boundary points of each subdomain which are used to impose the conditions Eq. (3) and Eq. (4). For guided modes, zeros-boundary conditions are imposed at

*y*= ±∞, while the conditions Eq. (16) or Eq. (17) hold for leaky modes. Define the first order differentiation matrix

*D*≡

*D*

^{(1)}such that its

*kj*-th entry

*d*=

_{kj}*dC*(

_{j}/dỹ*ỹ*), defined in [20]. The differentiation matrix of the second order derivative is then given by

_{k}*, we define the diagonal matrix Λ*

_{i}*(*

_{i}*f*) for an arbitrary function

*f*by

*are given by and*

_{i}*, we define the restricted matrices*

_{i}*I*in Eq. (23) is an identity matrix of size (

_{i}*N*+ 1) × (

_{i}*N*+ 1). Using the previous notations and definitions, we get for each subdomain Ω

_{i}*the eigenvalue problem where the (*

_{i}*N*− 1) × (

_{i}*N*+ 1) matrix

_{i}*K*is given by

It remains to point out that extension of the proposed approach to handle 3D problems has three main points that should be considered. First, rational Chebyshev basis functions would be incorporated in subdomains with semi-infinite extensions. Second, assuming full-vectorial formulation in terms of the magnetic field components *H _{x}* and

*H*, the form of approximate expansion of the field components is a double sum of products of basis functions of the form

_{y}*x*,

_{i}*y*), and (

_{j}*N*+1) and (

*M*+ 1) denote the number of basis functions in the direction of the

*x*-axis and

*y*-axis, respectively. Finally, the conditions imposed at horizontal and vertical interfaces are the continuity of the field components

*E*and

_{z}*H*(see [10] and references therein).

_{z}## 3. Numerical results

To demonstrate the efficiency of our approach, we used our suggested approach to perform modal solution of some challenging optical waveguides.

#### 3.1. Symmetric step-index slab waveguide

The first example is a three layer symmetric single-mode step-index slab waveguide [21]. It consists of core with refractive index *n*_{c} = 1.5 sandwiched between a cladding with refractive index *n*_{cl} = 1.48. The guiding layer width is 2 *μ*m and the operating wavelength is 1.5 *μ*m. The analytic value of the effective refractive index,
${n}_{\text{eff}}^{\text{exact}}=1.489305274606916$ is obtained by solving the characteristic equation for TM modes using the analytical method described in [22]. The normalized propagation constant *b* is defined by
$b=\left({n}_{\text{eff}}^{2}-{n}_{\text{cl}}^{2}\right)/\left({n}_{\text{c}}^{2}-{n}_{\text{cl}}^{2}\right)$ [23]. Its exact value is *b*^{exact} = 0.455715434937231. The relative error in the normalized propagation constant is given by

*b*are listed Table 1 as the total number of basis function,

*N*is increased.

The proposed method outperforms the Matrix Method [21], where the minimum obtained relative error is only 1.818E–04 using a truncated series of *N* = 400 terms, while the proposed method provides a relative error of 4.909E–14 at *N* = 80 (five times less than) matrix method. The semi-log plot of Δ*b* versus *N* is shown in Fig. (4b) where a very high convergence rate is achieved. This assures that the interface boundary conditions are considered accurately. The normalized fundamental TM field, |*H _{x}*| is shown in Fig. (4a).

#### 3.2. Planar waveguide with continuous Gaussian RIP

We consider a planar waveguide with continuous Gaussian RIP [12, 24] shown in Fig. (6a). It is described by

where*n*= 1.512,

_{s}*C*= 0.0833,

*t*= 2.66

*μ*m and the operating wavelength

*λ*= 0.6328

*μ*m. To apply our technique, we divide the computation domain (−∞, ∞) into two semi-infinite subdomains as shown in Fig. (5).

It is worth noting that the collocation points have high concentration near the center of the waveguide where field is more confined. The treatment in [12] is based on single domain expansion where a domain truncation technique was used to truncate the computational domain from (−∞, ∞) to a finite computational window [−*W*, *W*]. Although this approach was successful, the appropriate value of *L* should vary as different graded-index profiles are considered. This problem is handled here by dividing the computational domain into two semi-infinite subdomains. In each subdomain, the field is expanded using rational Chebyshev functions, taking into account the semi-infinite extension. Convergence of the effective refractive index of the first TE mode is shown in Table 2 with the third column representing an approximate relative error defined by
$\left|{n}_{\text{eff}}^{\text{new}}-{n}_{\text{eff}}^{\text{old}}\right|/{n}_{\text{eff}}^{\text{new}}$ and the first ten TE modes are listed in Table 3 and the first four TE modes are shown in Fig. (6b).

The results reported in Table 3 show a great agreement of the present method with method in [12]. However, the former method based on domain truncation needs appropriate estimation of the size of the truncated domain i.e. the value of the parameter *W*, especially for cases when the field decays slowly as *y* → ±∞.

#### 3.3. Discontinuous planar waveguide with strong discontinuity

The third example studies the case of a discontinuous RIP given by

*n*= 1.512,

_{s}*C*= 0.0833,

*t*= 2.66

*μ*m and

*λ*= 0.6328

*μ*m. This refractive index profile is in Fig. (7a). The field profile of the fundamental TE mode obtained by the proposed method is shown in Fig. 7b. This example was studied in [11] using a pseudo-spectral domain decomposition with Laguerre-Gauss functions used in the outer subdomains and using the same approach with domain truncation in [12]. The results shown in Table 4 depict the convergence behaviour of the proposed method as the total number of basis functions increases. The approximate relative error is also given. The obtained convergent value by the proposed method is in good agreement with that obtained in [11] which is 1.579107.

#### 3.4. Slow plasmon resonant nano-structure

We now turn to studying plasmonic structures. First, we study slow surface plasmon-polariton (SPP) modes for a symmetric insulator-metal-insulator (IMI) structure shown in Fig. 8. The propagation constants of long (short) range L(S)R-SPP modes, *k*_{l(s)rsp} are given, respectively by the exact formula given in [25]. The dielectric functions of the metal and dielectric are given by *ε _{d}* = 1 and

*ε*= −23 + 1.69

_{m}*j*, respectively. The effective refractive index for long (short) range SPP modes

*n*

_{l(s)}and the propagation length

*L*

_{l(s)}are then given by

*δn*

_{l(s)}and the relative error in the propagation length

*δL*

_{l(s)}are defined by

We consider a case of a film with thickness *t* = 40 nm. The effective index and propagation length for the long range SPP mode for are reported in Table 5 and those for short range SPP mode are given in Table 6. The field profiles for the normalized ℜ(*H _{x}*(

*y*)) for the LR-SPP and SR-SPP modes are shown in Fig. 8. The rapidly decaying relative errors for both the effective index and propagation length reflect the excellent performance and accuracy of the proposed method. In Table 7, films with different thickness values are reported, where

*N*is the number of basis functions in each subdomain.

#### 3.5. The gap SPP mode

We now consider the calculation of the gap SPP (G-SPP) mode in the symmetric metal-insulator metal(MIM) structure shown in Fig. 9 with the field plot of the G-SPP mode. This mode is the only surviving mode for all values of gap thickness, *t* (nm). The G-SPP propagation constant *k _{gsp}* is given by the exact expression [26, 27]

*n*and the propagation length

_{gsp}*L*with the corresponding relative errors given, respectively by

_{gsp}*δn*and

_{gsp}*δL*are reported in Table 8 for gap of thicknesses

_{gsp}*t*= 20nm and those for

*t*= 1000nm are reported in Table 9. The exact values were obtained by solving Eq. (34).

In spite of the strong damping of the field inside the metal region in the MIM nano-structure, the convergence of the domain truncation method [12] depends on the size of the computational window [−*W*, *W*]. In addition, the numerical results shown in Fig. (10a) and (10b) confirm that convergence of the domain truncation approach cannot be ensured by just using a large value of the computational window parameter *W*.

#### 3.6. Anti-resonant reflecting optical waveguide (ARROW)

Finally, we consider an anti-resonant reflecting optical waveguide (ARROW) structure [28, 29] as shown in Fig. 11a. The operating wavelength *λ* = 0.6328 *μ*m, the refractive indices are *n _{a}* = 1 for air and

*n*= 3.85 for the substrate. The refractive indices for the first and second claddings are

_{s}*n*

_{1}= 2.3 and

*n*

_{2}= 1.46, respectively.

Their thicknesses are *d*_{1} = 0.142*λ* and *d*_{2} = 3.15*λ*, respectively. The refractive index of the core layer is *n _{c}* = 1.46 and its thickness

*d*= 6.3

_{c}*λ*. The results of the normalized propagation constant,

*β/k*

_{0}are reported in Table 10. The normalized transverse electric field

*E*corresponding to TE

_{x}_{1}and TE

_{18}leaky modes are shown in Fig. 11b. The inset figures demonstrate the effective handling of interface boundary conditions.

## 4. Conclusion

In conclusion, an efficient multi-domain pseudo-spectral method based on rational Chebyshev functions is introduced. Careful selection of the rational maps and employing the rational Chebyshev functions had a significant improvement of the accuracy over traditional domain truncation techniques. The efficiency of our approach was demonstrated through accurate calculation of plasmonic and leaky modes propagation constants. The reported results showed the precision and accuracy of the proposed method with an excellent agreement with other methods in the literature.

## References and links

**1. **Y. C. Shih, “The mode-matching method,” in *Numerical Techniques for Microwave and Millimeter-Wave Passive Structures*, T. Itoh, ed. (Wiley, 1989), pp. 592–621.

**2. **A. S. Sudbo, “Film mode matching: a versatile numerical method for vector mode field calculations in dielectric waveguides,” Pure Appl. Opt. **2**(3), 211–233 (1993). [CrossRef]

**3. **P. Bienstman and R. Baets, “Optical modelling of photonic crystals and VCSELs using eigenmode expansion and perfectly matched layers,” Opt. Quantum Electron. **33**, 327–341 (2001). [CrossRef]

**4. **Y. P. Chiou, Y. C. Chiang, C. H. Lai, C. H. Du, and H. C. Chang, “Finite difference modeling of dielectric waveguides with corners and slanted facets,” J. Lightwave Technol. **27**, 2077–2086 (2009). [CrossRef]

**5. **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 (2007). [CrossRef]

**6. **S. Selleri, L. Vincetti, A. Cucinotta, and M. Zoboli, “Complex FEM modal solver of optical waveguides with PML boundary conditions,” Opt. Quantum Electron. **33**, 359–371 (2001). [CrossRef]

**7. **Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. **11**, 1545–1549 (1993). [CrossRef]

**8. **C. C. Huang, “Solving the full anisotropic liquid crystal waveguides by using an iterative pseudospectral-based eigenvalue method,” Opt. Express **14**, 3363–3378 (2011). [CrossRef]

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

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

**11. **C. C. Huang, C. C. Huang, and J. Y. Yang, “An efficient method for computing optical waveguides with discontinuous refractive index profile using spectral collocation method with domain decomposition,” J. Lightwave Technol. **20**, 2284–2296 (2003). [CrossRef]

**12. **J. Zhu, X. Zhang, and R. Song, “A unified mode solver for optical waveguides based on mapped barycentric rational chebyshev differentiation matrix,” J. Lightwave Technol. **28**, 1802–1810 (2010). [CrossRef]

**13. **A. F. Oskooi, L. Zhang, Y. Avniel, and S. G. Johnson, “The failure of perfectly matched layers, and towards their redemption by adiabatic absorbers,” Opt. Express **16**, 11376–11392 (2008). [CrossRef] [PubMed]

**14. **H. Rogier and D. D. Zutter, “Berenger and leaky modes in optical fibers terminated with a perfectly matched layer,” J. Lightwave Technol. **20**, 1141–1148 (2002). [CrossRef]

**15. **D. B. Haidvogel and T. Zang, “The accurate solution of Poisson equation by expansion in Chebyshev polynomials,” J. Comput. Phys. **30**, 167–180 (1979). [CrossRef]

**16. **J. P. Boyd, “The optimization of convergence for Chebyshev polynomial methods in an unbounded domain,” J. Comput. Phys. **45**, 42–79 (1982). [CrossRef]

**17. **C. C. Huang, “Numerical calculations of ARROW structures by pseudo-spectral approach with Murs absorbing boundary conditions,” Opt. Express **14**, 11631–11652 (2006). [CrossRef] [PubMed]

**18. **F. Tisseur and Karl Meerbergen, “The quadratic eigenvalue problem,” SIAM Rev. **43**, 235–286 (2001). [CrossRef]

**19. **J. -P. Dedieu and F. Tisseur, “Perturbation theory for homogeneous polynomial eigenvalue problems,” Linear Algebra Appl. **358**, 71–94 (2003). [CrossRef]

**20. **G.-X. Fan, Q. H. Liu, and J. S. Hesthaven, “Multidomain pseudospectral time-domain simulations of scattering by objects buried in lossy media,” IEEE Trans. Geosci. Remote Sens. **40**, 1366–1373 (2002). [CrossRef]

**21. **L. Wang and C. S. Hsiao, “A matrix method for studying TM modes of optical planar waveguides with arbitrary index profiles,” IEEE J. Quantum Electron. **37**, 1654–1660 (2001). [CrossRef]

**22. **K. Kawano and T. Kitoh, “Analytical Methods,” in *Introduction to Optical Waveguide Analysis: Solving Maxwell’s Equations and the Schrödinger Equation* (John Wiley & Sons, 2004), pp. 15–18.

**23. **C. Vassallo, “1993–1995 Optical mode solvers,” Opt. Quantum Electron. **29**, 95–114 (1997). [CrossRef]

**24. **S. Banerjee and A. Sharma, “Propagation characteristics of optical waveguiding structure by direct solution of the helmholtz equation for total fields,” J. Opt. Soc. Amer. A **6**, 1884–1894 (1989). [CrossRef]

**25. **S. I. Bozhevolnyi and T. Sondergaard, “General properties of slow-plasmon resonant nanostructures: nano-antennas and resonators,” Opt. Express **15**, 10869–10877 (2007). [CrossRef] [PubMed]

**26. **E. N. Economou, “Surface plasmons in thin films,” Phys. Rev. **182**, 539–554 (1969). [CrossRef]

**27. **R. Zai, M. D. Selker, P. B. Catrysee, and M. L. Brongersma, “Geometries and materials for subwavelength surface plasmon modes,” J. Opt. Soc. Am. A (21), 2442–2446 (2004).

**28. **E. Anemogiannis, E. N. Glytsis, and T. K. Gaylord, “Determination of guided and leaky modes in lossless and lossy planar multilayer optical waveguides: reflection pole method and wavevector density method,” J. Lightwave Technol. **17**, 929–941 (1999). [CrossRef]

**29. **E. I. Golant and K. M. Golant, “New method for calculating the spectra and radiation losses of leaky waves in multilayer optical waveguides,” Tech. Phys. **51**, 1060–1068 (2006). [CrossRef]