## Abstract

A relatively simple and efficient numerical method is developed for analyzing the scattering of light by a layered cylindrical structure of arbitrary cross section surrounded by a layered background. The method significantly extends an existing vertical mode expansion method (VMEM) for circular or elliptic cylindrical structures. The original VMEM and its extension give rise to effective two-dimensional formulations for the three-dimensional scattering problems of layered cylindrical structures. The extended VMEM developed in this paper uses boundary integral equations to handle the two-dimensional Helmholtz equations that appear in the vertical mode expansion process. The method is applied to analyze the transmission of light through subwavelength apertures in metallic films and the scattering of light by metallic nanoparticles.

© 2015 Optical Society of America

## 1. Introduction

In many applications, it is necessary to analyze the scattering of light by three-dimensional (3D) objects that are one-dimensional (1D), i.e., layered, in different parts of the structure. We are particularly interested in a single layered cylindrical object surrounded by an infinite layered background. Simple examples of such doubly-layered structures include a metallic film with a cylindrical hole and a cylindrical metallic nanoparticle on a substrate. Existing numerical methods for these scattering problems include general numerical methods, such as the finite-difference time-domain (FDTD) method [1], the finite element method [2], the volume and surface integral equation methods [3–9], as well as special numerical or semi-analytic methods, such as the mode matching method (also called modal method or mode expansion method) [10–16]. While the general numerical methods are widely used, the special methods tailor-made for a class of structures could be simpler or more efficient.

The usual mode matching method is applicable to structures that are piecewise invariant in a spatial variable *z*, and it expands the electromagnetic field in each *z*-invariant layer using the eigenmodes of the layer. For 3D structures with high index-contrast and metallic components, the mode matching method is not very efficient, since a large number of eigenmodes are needed, and these modes are two-dimensional, full vectorial, and expensive to calculate. Recently, a vertical mode expansion method (VMEM) was developed for analyzing circular or elliptic cylindrical objects in a layered background [17, 18]. The method is related to the earlier work of Boscolo and Midrio [19] and its further extensions [20, 21] for photonic crystal slabs, and it has been used to analyze light transmission through apertures in metallic films and scattering of light by metallic nanoparticles.

The structures considered in previous works [17, 18] are separable in each layered region. For circular and elliptic objects, the polar and elliptic coordinates have been used in the horizontal plane perpendicular to the axis of the cylinder, and the solutions in each region are obtained by analytic or numerical separation of variables. In this paper, we develop a general VMEM for cylindrical objects with arbitrary cross sections. In the layered regions, the structures are in general not separable in the horizontal variables. Our approach is to use boundary integral equations (BIEs) to handle the 2D Helmholtz equations that appear in the vertical mode expansion process. We illustrate the new VMEM by numerical examples for plasmonics applications [22].

## 2. The VMEM

As in [17, 18], we consider a cylindrical region *S*_{1}, where the relative permittivity and relative permeability are *ε*^{(1)} (*z*) and *μ*^{(1)} (*z*), respectively. A Cartesian coordinate system is chosen so that *z* is identified as the “vertical” variable, and the *xy* plane is the horizontal plane. Outside *S*_{1} is the infinite region *S* where the material parameters are *ε*^{(0)} (*z*) and *μ*^{(0)} (*z*). As simple examples, we show a cylindrical aperture in a metallic film in Fig. 1(a) and a cylindrical particle on a substrate in Fig. 1(b). The actual 3D part of the structure is limited to 0 < *z* < *D*, so that the media in the top (*z* > *D*) and bottom (*z* < 0) are homogeneous, i.e., *ε*^{(}^{l}^{)} and *μ*^{(}^{l}^{)} are constants and independent of *l* for either *z* < 0 or *z* > *D*. The cross sections of *S*_{1} and *S*_{0} are 2D domains Ω_{1} and Ω_{0} in the *xy* plane. We consider the scattering problem for an incident wave {**E**^{(}^{i}^{)},**H**^{(}^{i}^{)}} given in the top homogeneous medium, where **E** and **H** denote the electric field and the magnetic field multiplied by the free space impedance, respectively.

After a truncation of the *z* variable by the perfectly matched layer (PML) technique [23–28], the *z* components of the total field in region *S _{l}* can be expanded as

*p*=

*e*or

*p*=

*h*, is the profile of the

*j*th transverse electric (TE) or transverse magnetic (TM) mode in region

*S*, respectively, ${\eta}_{j}^{(l,p)}$ is the corresponding propagation constant, {

_{l}**E**

^{(}

^{l}^{)},

**H**

^{(}

^{l}^{)}} is the solution excited by the incident wave for an infinite 1D structure with material profiles

*ε*

^{(}

^{l}^{)}and

*μ*

^{(}

^{l}^{)}, and ${V}_{j}^{(l,p)}$ satisfies the following 2D Helmholtz equation

*. Other components of the electromagnetic field can also be expanded. In particular, it is necessary to expand the horizontal tangential components*

_{l}*H*and

_{τ}*E*, where

_{τ}*τ*is a unit tangential vector along Γ (the common boundary of Ω

_{0}and Ω

_{1}in the

*xy*plane). These expansions involve the normal and tangential derivatives of ${V}_{j}^{(l,p)}$, i.e., ${\partial}_{v}{V}_{j}^{(l,p)}$ and ${\partial}_{\tau}{V}_{j}^{(l,p)}$, on Γ, where

*ν*is a unit normal vector of Γ. More details about these expansions can be found in [17].

Notice that the vertical modes
$\{{\varphi}_{j}^{(l,p)},{\eta}_{j}^{(l,p)}\}$ are first calculated by solving some simple 1D eigenvalue problems [17, 20, 29]. If *z* is truncated and discretized by *N* points, we obtain *N* numerical TE modes and *N* numerical TM modes. The functions
${V}_{j}^{(l,p)}$ are unknown, but each of them satisfies a constant-coefficient 2D Helmholtz equation. The approach used in [17, 18] is to set up a linear system for all
${V}_{j}^{(l,p)}$ on Γ by matching *H _{z}*,

*E*,

_{z}*H*and

_{τ}*E*on the vertical boundary of

_{τ}*S*

_{1}and

*S*

_{0}, and then construct ${V}_{j}^{(l,p)}$ in Ω

*based on separation of variables in the polar or elliptic coordinates. In order to set up the linear system, it is necessary to express ${\partial}_{v}{V}_{j}^{(l,p)}$ on Γ in terms of ${V}_{j}^{(l,p)}$ on Γ. This requires the so-called Dirichlet-to-Neumann (DtN) map ${\mathbf{\Lambda}}_{j}^{(l,p)}$ satisfying*

_{l}The DtN map is a linear operator for functions defined on Γ. If Γ is discretized by *M* points, then
${\mathbf{\Lambda}}_{j}^{(l,p)}$ can be approximated by an *M × M* matrix.

In our case, the Helmholtz equation (3) is in general not separable in Ω* _{l}*, but a BIE can be used to establish a relation between
${V}_{j}^{(l,p)}$ and
${\partial}_{v}{V}_{j}^{(l,p)}$ on Γ. It turns out that the BIE producing the inverse of the DtN map, i.e., the Neumann-to-Dirichlet (NtD) map denoted as
${\mathbf{N}}_{j}^{(l,p)}$, is easier to implement. Therefore, we choose to solve
${\partial}_{v}{V}_{j}^{(l,p)}$ on Γ. Matching the four components

*H*,

_{z}*E*,

_{z}*H*and

_{τ}*E*on the vertical boundary of

_{τ}*S*

_{1}(and

*S*

_{0}), we obtain the following linear system:

Assuming the truncated *z* variable is discretized by *N* points (*z _{i}* for 1 ≤

*i*≤

*N*) and Γ is discretized by

*M*points (

*r**for 1 ≤*

_{m}*m*≤

*M*), then the above is a (4

*NM*) × (4

*NM*) linear system. The four elements of the unknown vector

**are given by**

*x**q*= 1, 2, 3, 4 for (

*l*,

*p*) = (0,

*e*), (0,

*h*), (1,

*e*) and (1,

*h*), respectively. The coefficient matrix is given in 4

*×*4 superblocks, where each superblock is an

*N × N*block matrix or an (

*NM*)

*×*(

*NM*) matrix. The (

*i*,

*j*) blocks of the superblocks are

*M×M*matrices given by

**T**is a matrix (the differentiation matrix) that approximates the tangential derivative operator on Γ. The vectors

*b*_{1},

*b*_{2},

*b*_{3}and

*b*_{4}in the right hand side of Eq. (5) are related to ${H}_{z}^{(1)}-{H}_{z}^{(0)}$, ${E}_{z}^{(1)}-{E}_{z}^{(0)}$, ${H}_{\tau}^{(1)}-{H}_{{}^{\tau}}^{(0)}$ and ${E}_{\tau}^{(1)}-{E}_{\tau}^{(0)}$, respectively, and they are exactly the same as those given in [17].

The matrix **T** can be constructed using interpolations by trigonometric functions [17, 30]. In the next section, we present a BIE method for computing the NtD maps.

## 3. NtD maps by boundary integral equations

For Eq. (3) in 2D domain Ω* _{l}*, the Green’s function is

**= (**

*r**x,y*) and

**= (**

*r*′*x′*,

*y′*) are two points in the

*xy*plane, and ${H}_{m}^{(1)}$ is the

*m*th order Hankel function of first kind. To simplify the notations, we drop the subscript

*j*and the superscripts (

*l, p*) for ${V}_{j}^{(l,p)}$ and ${\eta}_{j}^{(l,p)}$, and denote the Green’s function by

*G*instead of ${G}_{j}^{(l,p)}$.

In the bounded domain Ω_{1}, the Green’s representation theorem is

**∈ Ω**

*r*_{1}, where

*ν*(

**) is the outward unit normal vector of Γ at**

*r*′**∈ Γ, pointing into the exterior domain Ω**

*r*′_{0}. Assuming the boundary Γ is smooth, when

**tends to**

*r*

*r*_{*}∈ Γ, the first term in the right hand side of Eq. (7) is continuous, but the second term (with the minus sign) has a jump discontinuity which contributes an extra

*V*(

*r*_{∗})/2 [31,32]. This gives rise to the BIE

*r*_{∗}∈ Γ. Let

**S**and

**K**(which depend on

*j*,

*l*and

*p*) be the boundary integral operators given by

*ϕ*defined on Γ, then Eq. (8) can be written as

**I**is the identity operator. Therefore, the NtD map is given by

Although the DtN map is formally given by **Λ** = **S**^{−1}(**I**+**K**), numerically it is not very stable to evaluate **S**^{−1}, since when Γ is smooth, **S** is a compact operator in the vector space of continuous functions defined on Γ, and its eigenvalues accumulate at zero [31, 32].

For the exterior domain Ω_{0}, we use the same unit normal vector *ν* on Γ, then the Green’s representation theorem is

**∈ Ω**

*r*_{0}, and the BIE is

*r*_{∗}∈ Γ. The above can be written as

_{0}

When Γ is discretized by *M* points, the operators **S**, **K** and **N** can be approximated by *M × M* matrices. If Γ is given in a parametric form with period 1 as

*t*, that is,

*r**=*

_{i}**(**

*r**t*) for

_{i}*t*=

_{i}*i/M*and 1

*≤ i ≤ M*. If

**B**is a boundary integral operator that maps a function

*ϕ*to another function

*ψ*, i.e.,

*B*(

*r**,*

_{*}**) is the kernel of the operator, we can approximate**

*r*′**B**by an

*M × M*matrix (which is also denoted by

**B**), such that

From Eq. (18), we have

*w*(

*t*) = |d

**(**

*r**t*)/d

*t*|. The function

*ϕ*(

**(**

*r**t*)) can be approximated by its trigonometric interpolation [30]

*L*(

*t*) = sin(

*Mπt*)/[

*M*tan(

*πt*)] if

*M*is an even integer. Notice that

*L*(0) = 1 and

*L*(

*t*) = 0 for 1 ≤

_{j}*j*<

*M*. Therefore, the (

*i, j*) entry of matrix

**B**is

Since the integral operators **S** and **K** have logarithmic singularities at *r*_{*} = *r**′*, special quadrature methods are needed to evaluate *B _{ij}*. For that purpose, it is convenient to use

*ζ*=

*t − t*so that the logarithmic singularity is located at

_{i}*ζ*= 0 and

*ζ*= 1. Using the periodicity of

**(**

*r**t*)and

*L*(

*t*), we have

The above can be evaluated by Alpert’s hybrid Gauss-trapezoidal quadrature rules [33]. For a function *g*(*ζ*) with logarithmic singularities at *ζ* = 0 and *ζ* = 1, Alpert’s quadrature rules give

*h*= 1

*/M*,

*K*

_{1}and

*K*

_{2}are positive integers,

*γ*and

_{k}*δ*for 1

_{k}*≤ k ≤ K*

_{1}, are the weights and nodes of the method. We use a 10-th order method for which

*K*

_{1}= 10,

*K*

_{2}= 6, and

*γ*and

_{k}*δ*are given in Appendix A.

_{k}When the matrix approximations of **S** and **K** are constructed, we can find **N** for domains Ω_{1} and Ω_{0} by Eq. (12) and Eq. (16), respectively. Notice that the identity operator in these equations becomes the *M × M* identity matrix. When the vertical variable *z* is discretized by *N* points, there are 2*N* vertical modes in each cylindrical region. Therefore, it is necessary to calculate 4*N* NtD maps. Since each NtD map can be calculated efficiently using only *O*(*M*^{3}) operations, the CPU time needed for computing all NtD maps is still negligible compared with the time needed for solving the final linear system (5).

## 4. Apertures in metallic films

To demonstrate the new VMEM based on BIEs and NtD maps, we analyze a few subwavelength apertures in a silver film. First, we validate our method by analyzing the transmission of light through an elliptic hole in a silver film [18, 34]. The structure was first studied by Zakharian *et al.* [34] to explore the extraordinary optical transmission (EOT) phenomenon [35, 36]. The semi-axes of the elliptic hole are 400 nm and 50 nm, respectively. The thickness of the silver film is *D* = 124 nm. For incident waves with a free space wavelength 1 *μ*m, the refractive index of the silver is assumed to be 0.226 + 6.99i. The new VMEM and the method of [18] produce indistinguishable numerical results, if they use the same discretizations for *z* and Γ. Since the VMEM of [18] solves a linear system for
${V}_{j}^{(l,p)}$ on Γ, uses numerical separation of variables in elliptic coordinates to construct the DtN maps
${\mathbf{\Lambda}}_{j}^{(l,p)}$, it is significantly different from the new VMEM based on BIEs and NtD maps.

Next, we consider a C-shaped aperture in a silver film of thickness *D* = 120 nm. The top view of the structure is shown in the top-left panel of Fig. 2, where *d* = 100 nm. Our numerical results are shown in the other three panels of Fig. 2. The top-right panel of Fig. 2 shows the wavelength dependence of normalized transmission through a C-shaped air-hole for a normal incident plane wave with an electric field in the *y* direction. The refractive index of silver is taken from [37]. The normalized transmission is defined as the ratio between the total transmitted power (or additional transmitted power if the film without the hole is not opaque) and the power of the incident wave impinging on the aperture. We observe that the normalized transmission reaches a maximum larger than 3 for a free space wavelength slightly larger than 1.3 *μ*m. In particular, if the wavelength is 1 *μ*m, the normalized transmission is 0.186. This agrees well with the result 0.20 of Shi and Hesselink [38] for a similar C-shaped hole with right angle corners.

We also consider the effect of filling the hole with a dielectric material. Earlier studies indicate that transmission could be enhanced if the hole is filled with a proper dielectric material [39, 40]. The bottom-left panel of Fig. 2 shows the dependence of the normalized transmission on the refractive index *n _{c}* of the medium in the hole. The wavelength is fixed at 1

*μ*m and the incident electric field is in the

*y*direction. Notice that a maximum transmission is obtained for

*n*1.97. In the bottom-right panel of Fig. 2, we show the dependence of the normalized transmission on the horizontal angle of the incident electric field. For the fixed wavelength 1

_{c}≈*μ*m and

*n*= 1.97, the normalized transmission decreases as the angle between the electric field and the

_{c}*y*axis increases from 0 to

*π/*2, and reaches a minimum when the electric field is in the

*x*direction. The results of Fig. 2 are obtained using

*N*= 87 and

*M*= 44 for discretizing

*z*and Γ, respectively. The linear system (5) involves 4

*NM*= 15312 unknowns.

For an example with less symmetry, we consider the aperture shown in the top-left panel of Fig. 3, also in a silver film with thickness *D* = 120 nm. The boundary Γ of this aperture is given explicitly (in *μ*m) by

*≤ t ≤*1. Similar to the C-shaped aperture, we consider normal incident plane waves and calculate the normalized transmission for different wavelengths, for different refractive indices of the dielectric material filling the aperture, and for different angles of the incident electric field in the horizontal plane. The results are shown in Fig. 3. Notice that the transmission curve shown in the bottom-right panel of Fig. 3 does not have a reflection symmetry with respect to

*θ*=

*π/*2. These results are also obtained with

*N*= 87 and

*M*= 44.

## 5. Metallic nanoparticles

Metallic nanoparticles are extremely important for nanophotonics applications, mainly due to the localized surface plasmon resonances that can dramatically enhance the field around the particles [22]. To demonstrate the VMEM developed in previous sections, we first consider a C-shaped cylindrical gold particle with the same cross section as the aperture shown in the top-left panel of Fig. 2 and *d* = 100 nm. The particle is surrounded by air and its height is *D* = 120 nm. The refractive index of gold is taken from [37]. In Fig. 4, we show the normalized scattering cross sections as functions of the wavelength, for a few different normal incident plane waves. The normalized scattering cross section is defined as the total scattered power divided by the intensity of the incident light and the top area of the particle. The different curves in Fig. 4 correspond to incident waves with different horizontal angles between the electric field and the *y* axis. Notice that the scattering cross section is independent of the horizontal angle for a wavelength around 1.1 *μ*m. If the wavelength is less than 0.5 *μ*m, the scattering cross section does not have a significant angle dependence.

In practice, the metallic nanoparticles are placed on a substrate. We consider a gold cylindrical particle on a substrate with a frequency-independent refractive index 1.5. The cross section of the particle is the same as the aperture shown in the top-left panel of Fig. 3, and its height is *D* = 120 nm. As before, we consider a few normal incident plane waves with their electric field pointing to different directions. Using the new VMEM, we obtain the results shown in Fig. 5. These results are obtained using the same discretizations for *z* and Γ, i.e. *N* = 87 and *M* = 44.

## 6. Conclusion

The VMEM developed in this paper is a special numerical method for analyzing the scattering of electromagnetic waves by a layered cylindrical object in a layered background. Different from the previous version developed for circular and elliptic cylindrical objects [17, 18], the current VMEM is applicable to cylindrical structures with arbitrary cross sections. It uses BIEs to construct the NtD maps, and establishes a linear system on the vertical boundary of the cylindrical regions. The method is highly competitive, since it gives an effective 2D formulation of the original 3D scattering problems, and it is relatively simple to implement. We have used a highly accurate Chebyshev pseudospectral method [30] to discretize the vertical variable *z*, and used Alpert’s 10-th order hybrid Gauss-trapezoidal rule [33] to solve the BIEs, so that the final linear system is not too large and can be solved by a direct method.

The VMEM should be useful for simulating plasmonic structures, since nanoparticles and apertures are often cylindrical due to the existing fabrication techniques. The method is currently being further improved to handle more complicated multiple cylindrical structures, periodic structures, and multi-layered structures.

## Appendix A

The hybrid Gauss-trapezoidal quadrature rules developed by Alpert [33] can be written as

*S*

_{middle}corresponds to the trapezoidal rule,

*S*

_{left}and

*S*

_{right}depend on the singularities (if any) of

*g*at

*ζ*= 0 and

*ζ*= 1, respectively. A logarithmic singularity at

*ζ*= 0 implies that for

*ζ*close to 0, where

*g*

_{0}and

*g*

_{1}are regular functions of

*ζ*. It is possible to have different singularities at 0 and 1, but if the singularities at the two ends are of the same type, then

*S*

_{left}and

*S*

_{right}can essentially be the same formula. If that is the case, a hybrid Gauss-trapezoidal quadrature rule involves two related positive integers

*K*

_{1}and

*K*

_{2},

*M*> 2

*K*

_{2}is a positive integer,

*h*= 1

*/M*,

*γ*and

_{k}*δ*depend on

_{k}*K*

_{1},

*K*

_{2}, and the singularities. The 10-th order formula for the logarithmic singularity has

*K*

_{1}= 10,

*K*

_{2}= 6, and the weights

*γ*and nodes

_{k}*δ*are given in Table 1.

_{k}## Acknowledgments

The authors would like to thank Dr. Wangtao Lu for some valuable suggestions. This work was partially supported by the Research Grants Council of Hong Kong Special Administrative Region, China, under Project CityU 11301914.

## References and links

**1. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: the Finite-Difference Time-Domain Method*, 2nd ed. (Artech House, 2000).

**2. **J. M. Jin, *The Finite Element Method in Electromagnetics*, 2nd ed. (John Wiley & Sons, 2002).

**3. **W. C. Chew, M. S. Tong, and B. Hu, *Integral Equation Methods for Electromagnetic and Elastic Waves* (Morgan & Claypool, 2008).

**4. **A. M. Kern and O. J. F. Martin, “Surface integral formulation for 3D simulations of plasmonic and high permittivity nanostructures,” J. Opt. Soc. Am. A **26**, 732–740 (2009). [CrossRef]

**5. **B. Gallinet, A. M. Kern, and O. J. F. Martin, “Accurate and versatile modeling of electromagnetic scattering on periodic nanostructures with a surface integral approach,” J. Opt. Soc. Am. A **27**, 2261–2271 (2010). [CrossRef]

**6. **J. M. Taboada, J. Rivero, F. Obelleiro, M. G. Araújo, and L. Landesa, “Method-of-moments formulation for the analysis of plasmonic nano-optical antennas,” J. Opt. Soc. Am. A **28**, 1341–1348 (2011). [CrossRef]

**7. **M. G. Araújo, J. M. Taboada, D. M. Solís, J. Rivero, L. Landesa, and F. Obelleiro, “Comparison of surface integral equation formulations for electromagnetic analysis of plasmonic nanoscatterers,” Opt. Express **20**, 9161–9171 (2012). [CrossRef] [PubMed]

**8. **C. Forestiere, G. Iadarola, G. Rubinacci, A. Tamburrino, L. Dal Negro, and G. Miano, “Surface integral formulations for the design of plasmonic nanostructures,” J. Opt. Soc. Am. A **29**, 2314–2327 (2012). [CrossRef]

**9. **T. V. Raziman, W. R. C. Somerville, O. J. F. Martin, and E. C. Le Ru, “Accuracy of surface integral equation matrix elements in plasmonic calculations,” J. Opt. Soc. Am. B **32**, 485–492 (2015). [CrossRef]

**10. **L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A **14**, 2758–2767 (1997). [CrossRef]

**11. **E. Silberstein, P. Lalanne, J.-P. Hugonin, and Q. Cao, “Use of grating theories in integrated optics,” J. Opt. Soc. Am. A **18**, 2865–2875 (2001). [CrossRef]

**12. **G. Granet and J. P. Plumey, “Parametric formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. A **4**, S145–S149 (2002). [CrossRef]

**13. **M. Hammer, “Hybrid analytical/numerical coupled-mode modeling of guided wave devices,” J. Lightw. Technol. **25**(9), 2287–2298 (2007). [CrossRef]

**14. **J. Mu and W. P. Huang, “Simulation of three-dimensional waveguide discontinuities by a full-vector mode-matching method based on finite-difference schemes,” Opt. Express **16**, 18152–18163 (2008). [CrossRef] [PubMed]

**15. **R. Wang, L. Han, J. Mu, and W. P. Huang, “Simulation of waveguide crossings and corners with complex mode-matching method,” J. Lightw. Technol. **30**, 1795–1801 (2012). [CrossRef]

**16. **V. Liu and S. Fan, “S^{4}: A free electromagnetic solver for layered periodic structures,” Computer Physics Communications **183**, 2233–2244 (2012). [CrossRef]

**17. **X. Lu, H. Shi, and Y. Y. Lu, “Vertical mode expansion method for transmission of light through a single circular hole in a slab,” J. Opt. Soc. Am. A **31**, 293–300 (2014). [CrossRef]

**18. **H. Shi and Y. Y. Lu, “Vertical mode expansion method for analyzing elliptic cylindrical objects in a layered background,” J. Opt. Soc. Am. A **32**, 630–636 (2015). [CrossRef]

**19. **S. Boscolo and M. Midrio, “Three-dimensional multiple-scattering technique for the analysis of photonic-crystal slabs,” J. Lightw. Technol. **22**, 2778–2786 (2004). [CrossRef]

**20. **L. Yuan and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing hole arrays in a slab,” J. Opt. Soc. Am. B **27**, 2568–2579 (2010). [CrossRef]

**21. **L. Yuan and Y. Y. Lu, “An efficient numerical method for analyzing photonic crystal slab waveguides,” J. Opt. Soc. Am. B **28**, 2265–2270 (2011). [CrossRef]

**22. **S. A. Maier, *Plasmonics: Fundamentals and Applications* (Springer, 2007).

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

**24. **W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modified Maxwells equations with stretched coordinates,” Microwave and Optical Technology Letters **7**, 599–604 (1994). [CrossRef]

**25. **H. Derudder, D. De Zutter, and F. Olyslager, “Analysis of waveguide discontinuities using perfectly matched layers,” Electron. Lett. **34**(22), 2138–2140 (1998). [CrossRef]

**26. **P. Bienstman, H. Derudder, R. Baets, F. Olyslager, and D. De Zutter, “Analysis of cylindrical waveguide discontinuities using vectorial eigenmodes and perfectly matched layers,” IEEE Trans. Microwave Theory Tech. **49**(2), 349–354 (2001). [CrossRef]

**27. **Y. Y. Lu and J. Zhu, “Propagating modes in optical waveguides terminated by perfectly matched layers,” IEEE Photon. Technol. Lett. **17**, 2601–2603 (2005). [CrossRef]

**28. **J. Zhu and Y. Y. Lu, “Asymptotic solutions of eigenmodes in slab waveguides terminated by perfectly matched layers,” J. Opt. Soc. Am. A **30**, 2090–2095 (2013). [CrossRef]

**29. **D. Song, L. Yuan, and Y. Y. Lu, “Fourier-matching pseudospectral modal method for diffraction gratings,” J. Opt. Soc. Am. A **28**, 613–620 (2011). [CrossRef]

**30. **L. N. Trefethen, *Spectral Methods in MATLAB* (Society for Industrial and Applied Mathematics, 2000). [CrossRef]

**31. **D. Colton and R. Kress, *Integral Equation Methods in Scattering Theory* (John Wiley & Sons, 1983).

**32. **D. Colton and R. Kress, *Inverse Acoustic and Electromagnetic Scattering Theory*, 3rd ed. (Springer-Verlag, 2013). [CrossRef]

**33. **B. K. Alpert, “Hybrid Gauss-trapezoidal quadrature rules,” SIAM J. Sci. Comput. **20**(5), 1551–1584 (1999). [CrossRef]

**34. **A. R. Zakharian, M. Mansuripur, and J. V. Moloney, “Transmission of light through small elliptical apertures,” Opt. Express **12**, 2631–2648 (2004). [CrossRef] [PubMed]

**35. **T. W. Ebbesen, H. J. Lezec, G. F. Ghaemi, T. Thio, and P. A. Wolff, “Extraordinary optical transmission through sub-wavelength hole arrays,” Nature **391**, 667–669 (1998). [CrossRef]

**36. **F. J. García-Vidal, L. Martin-Moreno, T. W. Ebbesen, and L. Kuipers, “Light passing through subwavelength apertures,” Rev. Mod. Phys. **82**, 729–787 (2010). [CrossRef]

**37. **A. D. Raki, A. B. Djuriic, J. M. Elazar, and M. L. Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,” Appl. Opt. **37**, 5271–5283 (1998). [CrossRef]

**38. **X. Shi and L. Hesselink, “Design of a C aperture to achieve λ/10 resolution and resonant transmission,” J. Opt. Soc. Am. B **21**, 1305–1317 (2004). [CrossRef]

**39. **J. Olkkonen, K. Kataja, and D. Howe, “Light transmission through a high index dielectric-filled sub-wavelength hole in a metal film,” Opt. Express **13**, 6980–6989 (2005). [CrossRef] [PubMed]

**40. **H. Xu, P. Zhu, H. G. Craighead, and W. W. Webb, “Resonantly enhanced transmission of light through subwave-length apertures with dielectric filling,” Opt. Commun. **282**, 1467–1471 (2009). [CrossRef]