Efficient vertical mode expansion method for scattering by arbitrary layered cylindrical structures

Hualiang Shi; Ya Yan Lu

doi:10.1364/OE.23.014618

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₁, where the relative permittivity and relative permeability are ε⁽¹⁾ (z) and μ⁽¹⁾ (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₁ is the infinite region S where the material parameters are ε⁽⁰⁾ (z) and μ⁽⁰⁾ (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₁ and S₀ are 2D domains Ω₁ and Ω₀ in the xy plane. We consider the scattering problem for an incident wave {E⁽ⁱ⁾,H⁽ⁱ⁾} 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.

Fig. 1 (a): A cylindrical aperture in a metallic film; (b): A cylindrical particle on a substrate.

Download Full Size | PDF

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

H_{z} = H_{z}^{(l)} + \frac{1}{μ^{(l)}} \sum_{j = 1}^{\infty} {[η_{j}^{(l, e)}]}^{2} ϕ_{j}^{(l, e)} V_{j}^{(l, e)},

E_{z} = E_{z}^{(l)} + \frac{1}{ε^{(l)}} \sum_{j = 1}^{\infty} {[η_{j}^{(l, h)}]}^{2} ϕ_{j}^{(l, h)} V_{j}^{(l, h)},

where

ϕ_{j}^{(l, p)} (l = 0, 1; j = 1, 2, \dots)

, for p = e or p =h, is the profile of the jth transverse electric (TE) or transverse magnetic (TM) mode in region S_l, respectively,

η_{j}^{(l, p)}

is the corresponding propagation constant, {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

\frac{\partial^{2} V_{j}^{(l, p)}}{\partial x^{2}} + \frac{\partial^{2} V_{j}^{(l, p)}}{\partial y^{2}} + {[η_{j}^{(l, p)}]}^{2} V_{j}^{(l, p)} = 0

in Ω_l. Other components of the electromagnetic field can also be expanded. In particular, it is necessary to expand the horizontal tangential components H_τ and E_τ, where τ is a unit tangential vector along Γ (the common boundary of Ω₀ and Ω₁ 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_{τ} 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 ${ϕ_{j}^{(l, p)}, η_{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₁ and S₀, and then construct $V_{j}^{(l, p)}$ in Ω_l 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 $Λ_{j}^{(l, p)}$ satisfying

Λ_{j}^{(l, p)} V_{j}^{(l, p)} = \partial_{v} V_{j}^{(l, p)} on Γ .

The DtN map is a linear operator for functions defined on Γ. If Γ is discretized by M points, then $Λ_{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 $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₁ (and S₀), we obtain the following linear system:

[\begin{matrix} A^{(11)} & 0 & A^{(13)} & 0 \\ 0 & A^{(22)} & 0 & A^{(24)} \\ A^{(31)} & A^{(32)} & A^{33} & A^{34} \\ A^{(41)} & A^{(42)} & A^{(43)} & A^{(44)} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{matrix}] = [\begin{matrix} b_{1} \\ b_{2} \\ b_{3} \\ b_{4} \end{matrix}] .

Assuming the truncated z variable is discretized by N points (z_i for 1 ≤ i ≤ N) and Γ is discretized by M points (r_m for 1 ≤ m ≤ M), then the above is a (4NM) × (4NM) linear system. The four elements of the unknown vector x are given by

x_{q} = [\begin{matrix} \partial_{v} v_{1}^{(l, p)} \\ \partial_{v} v_{2}^{(l, p)} \\ ⋮ \\ \partial_{v} v_{N}^{(l, p)} \end{matrix}], \partial_{v} v_{j}^{l, p} = [\begin{matrix} \partial_{v} V_{j}^{(l, p)} (r_{1}) \\ \partial_{v} V_{j}^{(l, p)} (r_{2}) \\ ⋮ \\ \partial_{v} V_{j}^{(l, p)} (r_{M}) \end{matrix}],

where 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

\begin{array}{l} A_{i j}^{(11)} = \frac{1}{μ^{(0)} (z_{i})} {[η_{j}^{(0, e)}]}^{2} ϕ_{j}^{(0, e)} (z_{i}) N_{j}^{(0, e)}, \\ A_{i j}^{(13)} = - \frac{1}{μ^{(1)} (z_{i})} {[η_{j}^{(1, e)}]}^{2} ϕ_{j}^{(1, e)} (z_{i}) N_{j}^{(1, e)}, \\ A_{i j}^{(22)} = \frac{1}{ε^{(0)} (z_{i})} {[η_{j}^{(0, h)}]}^{2} ϕ_{j}^{(0, h)} (z_{i}) N_{j}^{(0, h)}, \\ A_{i j}^{(24)} = - \frac{1}{ε^{(1)} (z_{i})} {[η_{j}^{(1, h)}]}^{2} ϕ_{j}^{(1, h)} (z_{i}) N_{j}^{(1, h)}, \\ A_{i j}^{(31)} = \frac{1}{μ^{(0)} (z_{i})} \partial_{z} ϕ_{j}^{(0, e)} (z_{i}) T N_{j}^{(0, e)}, \\ A_{i j}^{(32)} = i k_{0} ϕ_{j}^{(0, h)} (z_{i}) I, \\ A_{i j}^{(33)} = - \frac{1}{μ^{(1)} (z_{i})} \partial_{z} ϕ_{j}^{(1, e)} (z_{i}) T N_{j}^{(1, e)}, \\ A_{i j}^{(34)} = - i k_{0} ϕ_{j}^{(1, h)} (z_{i}) I, \\ A_{i j}^{(41)} = - i k_{0} ϕ_{j}^{(0, e)} (z_{i}) I, \\ A_{i j}^{(42)} = \frac{1}{ε^{(0)} (z_{i})} \partial_{z} ϕ_{j}^{(0, h)} (z_{i}) T N_{j}^{(0, h)}, \\ A_{i j}^{(43)} = i k_{0} ϕ_{j}^{(1, e)} (z_{i}) I, \\ A_{i j}^{(44)} = - \frac{1}{ε^{(1)} (z_{i})} \partial_{z} ϕ_{j}^{(1, h)} (z_{i}) T N_{j}^{(1, h)}, \end{array}

where i is the imaginary unit, and T is a matrix (the differentiation matrix) that approximates the tangential derivative operator on Γ. The vectors b₁, b₂, b₃ and b₄ in the right hand side of Eq. (5) are related to

H_{z}^{(1)} - H_{z}^{(0)}

,

E_{z}^{(1)} - E_{z}^{(0)}

,

H_{τ}^{(1)} - H_{^{τ}}^{(0)}

and

E_{τ}^{(1)} - E_{τ}^{(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

G (r, r^{'}) = \frac{i}{4} H_{0}^{(1)} (η | r - r^{'} |), r \neq r^{'},

where r = (x,y) and r′ = (x′,y′) are two points in the xy plane, and

H_{m}^{(1)}

is the mth 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

η_{j}^{(l, p)}

, and denote the Green’s function by G instead of

G_{j}^{(l, p)}

.

In the bounded domain Ω₁, the Green’s representation theorem is

V (r) = \int_{Γ} [G (r, r^{'}) \frac{\partial V}{\partial V} (r^{'}) - V (r^{'}) \frac{\partial G (r, r^{'})}{\partial v (r^{'})}] d s (r^{'})

for r ∈ Ω₁, where ν(r′) is the outward unit normal vector of Γ at r′ ∈ Γ, pointing into the exterior domain Ω₀. Assuming the boundary Γ is smooth, when r tends to 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

\frac{V (r_{*})}{2} = \int_{Γ} [G (r_{*}, r^{'}) \frac{\partial V}{\partial V} (r^{'}) - V (r^{'}) \frac{\partial G (r_{*}, r^{'})}{\partial v (r^{'})}] d s (r^{'})

for r_∗ ∈ Γ. Let S and K (which depend on j, l and p) be the boundary integral operators given by

(S ϕ) (r_{*}) = 2 \int_{Γ} G (r_{*}, r^{'}) ϕ (r^{'}) d s (r^{'}),

(K ϕ) (r_{*}) = 2 \int_{Γ} \frac{\partial G (r_{*}, r^{'})}{\partial v (r^{'})} ϕ (r^{'}) d s (r^{'})

for any function ϕ defined on Γ, then Eq. (8) can be written as

(I + K) V = S \partial_{v} V on Γ,

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

N = {(I + K)}^{- 1} S .

Although the DtN map is formally given by Λ = S⁻¹(I+K), numerically it is not very stable to evaluate S⁻¹, 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 Ω₀, we use the same unit normal vector ν on Γ, then the Green’s representation theorem is

V (r) = \int_{Γ} [V (r^{'}) \frac{\partial G (r, r^{'})}{\partial v (r^{'})} - G (r, r^{'}) \frac{\partial V}{\partial v} (r^{'})] d s (r^{'})

for r ∈ Ω₀, and the BIE is

\frac{V (r_{*})}{2} = \int_{Γ} [V (r^{'}) \frac{\partial G (r_{*}, r^{'})}{\partial v (r^{'})} - G (r_{*}, r^{'}) \frac{\partial V}{\partial v} (r^{'})] d s (r^{'})

for r_∗ ∈ Γ. The above can be written as

(K - I) V = S \partial_{v} V on Γ,

and it gives rise to the following NtD map for domain Ω₀

N = {(K - I)}^{- 1} S .

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

r = r (t), 0 \leq t \leq 1,

then it can be discretized by a uniform sampling of t, that is, r_i = r(t_i) for t_i = i/M and 1 ≤ i ≤ M. If B is a boundary integral operator that maps a function ϕ to another function ψ, i.e.,

ψ (r_{*}) = \int_{Γ} B (r_{*}, r^{'}) ϕ (r^{'}) d s (r^{'}), r_{*} \in Γ,

where B(r_*, r′) is the kernel of the operator, we can approximate B by an M × M matrix (which is also denoted by B), such that

[\begin{matrix} ψ (r_{1}) \\ ψ (r_{2}) \\ ⋮ \\ ψ (r_{M}) \end{matrix}] \approx B [\begin{matrix} ϕ (r_{1}) \\ ϕ (r_{2}) \\ ⋮ \\ ϕ (r_{M}) \end{matrix}] .

From Eq. (18), we have

ψ (r_{i}) = \int_{0}^{1} B (r_{i}, r (t)) ϕ (r (t)) w (t) d t,

where w(t) = |dr(t)/dt|. The function ϕ (r(t)) can be approximated by its trigonometric interpolation [30]

ϕ (r (t)) \approx \sum_{j = 1}^{M} ϕ (r_{j}) L (t - t_{j}),

where L(t) = sin(Mπt)/[M tan(πt)] if M is an even integer. Notice that L(0) = 1 and L(t_j) = 0 for 1 ≤ j < M. Therefore, the (i, j) entry of matrix B is

B_{i j} = \int_{0}^{1} B (r_{i}, r (t)) L (t - t_{j}) w (t) d t .

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_i so that the logarithmic singularity is located at ζ = 0 and ζ = 1. Using the periodicity of r(t)and L(t), we have

B_{i j} = \int_{0}^{1} B (r_{i}, r (ζ + t_{i})) L (ζ + t_{i} - t_{j}) w (ζ + t_{i}) d ζ .

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

\int_{0}^{1} g (ζ) d ζ \approx h \sum_{k = 1}^{K_{1}} γ_{k} [g (δ_{k} h) + g (1 - δ_{k} h)] + h \sum_{k = K_{2}}^{M - K_{2}} g (t_{k}),

where h = 1/M, K₁ and K₂ are positive integers, γ_k and δ_k for 1 ≤ k ≤ K₁, are the weights and nodes of the method. We use a 10-th order method for which K₁ = 10, K₂ = 6, and γ_k and δ_k are given in Appendix A.

When the matrix approximations of S and K are constructed, we can find N for domains Ω₁ and Ω₀ 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 2N vertical modes in each cylindrical region. Therefore, it is necessary to calculate 4N NtD maps. Since each NtD map can be calculated efficiently using only O(M³) 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 $Λ_{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.

Fig. 2 Top-left: a C-shaped aperture in a silver film; Top-right: transmission spectrum for a C-shaped hole; Bottom-left: transmission versus refractive index n_c; Bottom-right: transmission versus horizontal angle of incident electric field.

Download Full Size | PDF

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_c ≈ 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 μm and n_c = 1.97, the normalized transmission decreases as the angle between the electric field and the 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 4NM = 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

\begin{array}{l} x (t) = 0.15 \cos (2 π t) + 0.05 \sin (4 π t + 0.8), \\ y (t) = 0.1 \sin (2 π t) + 0.02 \cos (4 π t), \end{array}

for 0 ≤ 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.

Fig. 3 Top-left: A “triangular” aperture in a metallic film; Top-right: the dependence of normalized transmission on wavelength; Bottom-left: the dependence on the refractive index n_c of the material filling the aperture; Bottom-right: the dependence on the horizontal angle of the incident electric field.

Download Full Size | PDF

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.

Fig. 4 Normalized scattering cross sections of a gold C-shaped cylindrical particle, for normal incident waves with different horizontal angles between the electric field and the y axis.

Download Full Size | PDF

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.

Fig. 5 Normalized scattering cross sections of a gold cylindrical particle on a substrate with refractive index 1.5 for normal incident plane waves with different angles between the electric field and the y axis.

Download Full Size | PDF

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

\int_{0}^{1} g (ζ) d ζ \approx S_{left} + S_{middle} + S_{right},

where 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

g (ζ) = g_{0} (ζ) \log (ζ) + g_{1} (ζ)

for ζ close to 0, where g₀ and g₁ 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₁ and K₂,

S_{left} = \sum_{k = 1}^{K_{1}} γ_{k} g (δ_{k} h), S_{right} = \sum_{k = 1}^{k_{1}} γ_{k} g (1 - δ_{k} h), S_{middle} = \sum_{k = K_{2}}^{M - K_{2}} g (k h),

and where M > 2K₂ is a positive integer, h = 1/M, γ_k and δ_k depend on K₁, K₂, and the singularities. The 10-th order formula for the logarithmic singularity has K₁ = 10, K₂ = 6, and the weights γ_k and nodes δ_k are given in Table 1.

Table 1. Nodes δ_k and weights γ_k of the 10-th order Alpert’s hybrid Gauss-trapezoidal rule for functions with a logarithmic singularity.

View Table

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⁴: 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]

δ_k	γ_k
1.17508 93812 27308E−03	4.56074 68820 84207E−03
1.87703 41298 31289E−02	3.81060 63223 84757E−02
9.68646 83914 26860E−02	1.29386 49972 89512E−01
3.00481 86680 02884E−01	2.88436 03814 08835E−01
6.90133 15571 73356E−01	4.95811 19143 44961E−01
1.29369 57380 83659E+00	7.07715 46005 94529E−01
2.09018 77297 98780E+00	8.74192 43652 85083E−01
3.01671 93131 49212E+00	9.66136 19865 15218E−01
4.00136 97478 72486E+00	9.95788 78660 78700E−01
5.00002 56617 93423E+00	9.99866 57874 23845E−01

Efficient vertical mode expansion method for scattering by arbitrary layered cylindrical structures

Abstract

1. Introduction

2. The VMEM

3. NtD maps by boundary integral equations

4. Apertures in metallic films

5. Metallic nanoparticles

6. Conclusion

Appendix A

Acknowledgments

References and links

Cited By

Figures (5)

Tables (1)

Equations (30)

Optics Express