1. Introduction

Accurate numerical solutions of three-dimensional electromagnetic scattering problems have become indispensable in the face of recent advances in nano-optics. When dealing with time-harmonic solutions for homogeneous or piecewise homogeneous scatterers, the surface integral equation (SIE) formulation [1] solved via the method of moments (MoM) [2] stands out as one of the most efficient and versatile techniques. Within this formulation, the unknowns are defined only over the surfaces and interfaces of the scatterers. Consequently, the resulting scheme benefits greatly from a strongly reduced number of discretization elements, superior scaling properties, flexible meshing strategies, and the fact that no artificial truncation of the background medium is required. This has led to a wide range of studies and applications of the SIE-MoM method in nano-optics, notably for plasmonic scatterers embedded in homogeneous media [3–6]. Recently, it has been successfully cultivated into large-scale plasmonic calculations on dense assemblies of nanoparticles [7]. Its flexibility in modeling realistic scatterer shapes has also been exploited and further extended, for instance, through the use of super-shape geometric transformations [8].

Another important feature of the SIE formulation, shared within the class of integral equation (IE)-based formalism [9], is the inclusion of dyadic Green’s functions (DGFs) as its kernels. The electromagnetic (EM) properties of the background medium are fully accounted for in the IEs through these DGFs and simple explicit functions of the material constants. Therefore, generalization of the IEs for homogeneous background media to those governing piecewise homogeneous ones basically amounts to usage of more general DGFs. Taking advantage of this, the volume integral equation (VIE) formulation has been successfully implemented for stratified media [10, 11]. In view of the ever-increasing relevance of light-matter interaction in multilayered environments, as made evident by current nanofabrication technology, it is highly desirable that the more favorable properties of the SIE-MoM approach are also brought to the same level of versatility.

In fact, the SIE-MoM scheme for stratified media has been extensively used in the RF and microwave regimes, either in the popular mixed-potential forms (see [12] and references therein) or in the new and more intuitive combined-field form [13]. The latter is devoid of invoking potential-type DGFs and the associated non-unique definitions of scalar potentials in the IEs. Instead, the hypersingular behavior of the field-type DGFs is dealt with through an alternative rigorous derivation using the pilot-vector approach [14], followed by integration by parts within the MoM representations of the IEs.

Very recently, this combined-field SIE-MoM was implemented for rigorous studies on local density of photonic states near plasmonic structures residing in multilayered substrates [15,16], where the problem is treated as the mathematically-equivalent scattering of fields generated by a Hertzian dipole. To accelerate computationally-expensive evaluations of the IE kernels, the authors employed closed-form approximations of the layered-medium DGFs, based on a well-tuned variant of the discrete complex image method (DCIM) [17]. While being competitive in terms of efficiency, robustness of the DCIM, in general, relies heavily on judicious choices of the sampling paths and the forms of the approximating functions [18], which tends to render its implementation laborious. Moreover, its capability in estimating the incurred error in the spatial domain is still under active research [19]. Its accuracy therefore needs to be ascertained externally, still requiring direct numerical integrations of the approximated integrals.

Another well-known difficulty in SIE-MoM implementations for stratified media is the evaluations of double surface integrals involving spatial-domain singularities of the layered-medium DGFs. Without careful handling of these singularities, the surface integrals cannot be (accurately) evaluated for overlapping or neighboring pairs of discretization elements. When these pairs are relatively far from layer interfaces, the relevant singularities are only due to direct-propagation (primary) terms in the DGFs which comprise no transmission and/or reflection at interfaces. Hence, they can be handled through singularity treatments developed for homogeneous-medium DGFs. These include, among others, the singularity extraction method [20–22], which has widespread usage in nanoplasmonics. On the other hand, when the overlapping or neighboring pairs are situated on or very close to layer interfaces, singularities associated with the indirect-propagation (secondary) terms in the DGFs also become important. Indeed, this does not occur for the systems studied in [15, 16] since the scatterers are placed sufficiently far from layer boundaries. A broad spectrum of scattering problems of practical interest, however, include situations where the scatterers reside on or straddle layer interfaces. A parallel treatment of these secondary-term singularities is therefore necessary for a versatile and well-rounded SIE-MoM scheme.

In this paper, we present an efficient and robust implementation of the SIE-MoM scheme for stratified media. The singularity extraction method is extended by introducing extra terms to also account for the secondary-term singularities of the layered-medium DGFs. The occurring spectral-domain integrals are efficiently and accurately evaluated by direct numerical integrations based on the generalized weighted averages method [23]. Furthermore, to dramatically accelerate the laborious matrix-filling process, we devise a simple yet robust spatial tabulation-interpolation scheme. The associated internal resonance problems of the SIE formulation are coped with by employing the popular Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) formulation [24] which has been shown to be very reliable for plasmonic problems [6]. A brief overview of the underlying formalism and the resulting MoM representations of the dyadic kernels are presented in Section 2. In Section 3, we discuss the details of the numerical implementations. The proposed scheme is further demonstrated by treating several representative problems in Section 4. Finally, several important aspects regarding generalization to more complex situations are elaborated in Section 5.

2. Formulation

We consider a homogeneous scatterer, with its surface denoted by S, embedded in a stratified medium as illustrated in Fig. 1(a). The regions inside and outside S are denoted by subscripts i and o, respectively. Cases with piecewise homogeneous scatterers can be similarly treated by properly unifying all the inner regions under the subscript i and their surfaces into S (see Section 5). The quantities of interest are the regional fields, E_ν(r) and H_ν(r), where ν = i, o. Time dependence of the form exp(−iωt) is assumed and suppressed throughout.

 

Fig. 1 (a) A homogeneous scatterer with EM properties ε_i and µ_i is embedded in a medium with stratified properties ε_o(r) and µ_o(r). (b) A schematic representation of an RWG function f(r) = ±L(r − p^±)/2A^±, r ∈ T^±. T⁺ and T⁻ denote the two triangles that form its supported area, outside of which f(r) vanishes. L is the length of the common edge and A^± is the area of T^±. Implied is the relation ∇ · f(r) = ±L/A^± for r ∈ T^±.

Download Full Size | PDF

Within the PMCHWT formulation, the surface integral equations to solve are given by

To numerically solve Eq. (1)J(r) and M(r) are first approximated as linear combinations of a set of basis functions, {f₁(r),…,f_N(r)}, which hereafter is assumed to consist of the standard Rao-Wilton-Glisson (RWG) functions [25]. S is therefore triangulated into a surface mesh containing a total of N distinct common edges of two triangles. To a given common edge indexed by j ∈ {1,…,N}, one assigns the RWG function f_j(r), which is defined on its supported area, $S_j=T_j^{+}\cup T_j^{-}$ (see Fig. 1(b)). The resulting set of basis functions, {f_j(r)}, furthermore, is also used as the set of testing functions, {f_j(r)}, in discretizing Eq. (1).

It can then be shown that, when the primary terms of the layered-medium DGFs are separately handled in closed form (see Subsection 3.3), the MoM-friendly representation of ${\widehat{D}}_{\mathrm{o}}^{\mathrm{E}}(\mathbf{r})$ takes the form [26]

Here we have used the shorthand notation $〈\dots 〉={\int }_{S_i}\mathrm{d}S{\int }_{S_j}\mathrm{d}{S}^{\prime }\dots$ and the common notations $k_n^2={\omega }^2{\epsilon }_n{\mu }_n$ and $k_{mn}^2={\omega }^2{\epsilon }_m{\mu }_n$. The transverse-magnetic and transverse-electric scalar Green’s functions, g^TM(r,r′) and g^TE(r,r′), can be written in cylindrical coordinates as semi-infinite integrals involving the Bessel function of order zero. They read

In a similar fashion, the MoM representation of ${\widehat{K}}_{\mathrm{o}}^{\mathrm{H}}(\mathbf{r})$ can be written as [9]

Unlike Eq. (4), Eq. (11) might include occurrence of additional line integrals when S_i and/or Sj straddle layer interfaces [13]. These line integrals are invoked when the interaction between f_i(r) and f_j(r′) involves both the terms with m = n and the terms with m ≠ n (see also Section 5).

Lastly, the matrix representations of ${\widehat{D}}_{\mathrm{o}}^{\mathrm{H}}(\mathbf{r})$ and ${\widehat{K}}_{\mathrm{o}}^{\mathrm{E}}(\mathbf{r})$ can readily be obtained by invoking duality. The resulting expressions are equivalent to the ones obtained by performing the simultaneous replacements, ^E ↔ ^H, µ_i ↔ ε_i, and g^TM ↔ g^TE, on Eqs. (4)–(9) and Eqs. (11)–(15).

3. Numerical implementations

3.1. Evaluation of Sommerfeld-type integrals

In view of the propagation factor F^α(k_ρ,z,z′) in Eq. (10), the scalar Green’s function g^α(r,r′) can be written as a sum of Sommerfeld-type integrals (SIs) [28] of the form

Numerical evaluation of the occurring SIs, i.e. Eq. (16), Eq. (17), and their spatial derivatives, is well-known to be challenging and time-consuming. This fact is mainly attributed to the oscillatory and slowly convergent nature of the integrands. Moreover, these integrands possess singularities, associated with guided-wave poles and branch points, on and/or near the integration path [29].

These difficulties can be mitigated by first dividing the integration range into two intervals delimited by a = k₀(n_max +1), where k₀ is the vacuum wavenumber and $n_{\max }^2=c^2\max [\Re ({\epsilon }_i{\mu }_i)]$ In the latter, i runs through all the layer indices of region o. This subdivision ensures that the integrand is free of singularities throughout the second interval, (a,∞). The remaining integration path is further deformed in the complex k_ρ-plane into a half-sine segment of amplitude −b, as illustrated in Fig. 2. Since Bessel functions exhibit exponential growth in the (positive and negative) direction of the imaginary axis, the value of b needs to be carefully chosen for an efficient integration along the deformed path. To this end, we found that the relation

 

Fig. 2 Deformed integration path in the complex k_ρ-plane. {ξ_n} denote a set of break points involved in evaluating the tail integral within the generalized weighted-averages scheme.

Download Full Size | PDF

The integral over (0,a) is evaluated using an adaptive quadrature based on Piessens et al.’s error estimator [30]. The tail integral, on the other hand, requires a more elaborate handling. When r_s = 0, the adaptive quadrature is also employed for the tail integral through the usual variable transformation into a finite integration range. On the other hand, when r_s is nonzero, the tail integrals are evaluated using the generalized weighted-averages (WA) method [23], which has been shown to be highly accurate and efficient in [31]. For completion, a typical implementation of this method is summarized in Appendix A.

It should be noted that, when r_s = Z = 0, some of the SIs are singular. For these SIs, the integrands without quasi-static contributions (see Subsection 3.2) are used instead. Importantly, the resulting nonsingular SIs are strictly required in mediating the tabulation process and the singularity extraction scheme.

3.2. Tabulation-interpolation (TI) scheme

The matrix-filling process in an SIE-MoM procedure unavoidably involves a very large number of SI evaluations. Although usage of the generalized WA method leads to a dramatic increase in efficiency, the resulting framework is generally still unacceptable in terms of computational cost. To remedy this, we employ a spatial tabulation scheme which is specifically tailored for efficient and accurate interpolations of the SIs.

The main challenges in a TI-based acceleration method are associated with the spatial-domain singularities and/or discontinuities of the SIs at r_s = Z = 0. For an accurate and effi-cient TI scheme, the terms containing these irregularities need to be excluded from the tabulated values, to be added back following every interpolation. It is therefore highly desirable that these terms are also expressible in closed form. Since the irregularities originate from the quasi-static (QS) image contributions in the SIs, attaining these desired terms (hereafter referred to as the QS terms) basically amounts to finding functions with the same asymptotic behavior as the QS image contributions. In fact, these QS terms have been derived in [32] for the occurring SIs in Eq. (4) when r, r′ ∈ m. The terms obtained therein correspond to

With the spectral-domain singularities properly canceled out, $g_m+{\tilde{g}}_m^{\mathrm{qs}},g_{mn}+{\tilde{g}}_{mn}^{\mathrm{qs}}$, and their derivatives (see Appendix B) are free of spatial-domain singularities and discontinuities, rendering them suitable for a TI procedure. A look-up table is thus allocated for each (combination) of the SI(s) required by Eqs. (5)–(9) and Eqs. (12)–(15), where each propagation term (indexed by k in Eqs. (16) and (17)) is also tabulated separately. This leads to the use of two-dimensional grids (with independent variables r_s and Z) for the tabulated functions involving g_m and three-dimensional grids (with independent variables r_s, Z_m, and Z_n) for those involving g_mn.

The multivariate interpolations are performed using the tensor product scheme [33], in which the problem is reduced into a series of one-dimensional (1D) interpolations. For these 1D interpolations, a piecewise polynomial interpolant is used to allow a local interpolation within each subgrid, involving a relatively low number of interpolation nodes per interpolation. Another benefit is that the placement of break points (knots) and the polynomial degrees of the piecewise interpolant can be adapted to the local behavior of the approximated function. In our implementation, the polynomial degrees are kept constant across different subgrids, whereas the distances of the knots to the axes origin increase quadratically along each axis. The latter is to account for the fact that the approximated functions are generally more rapidly-varying for smaller values of the independent variables. Meanwhile, the former is very beneficial in terms of efficiency when the interpolating polynomials are evaluated via the barycentric Lagrange formula [34]. The weights involved in this formula are independent of the tabulated values. Thus, they can be precomputed, preceding each series of 1D interpolations with a common interpolation point, even for multiple look-up tables.

We finally note that, for robustness, the interpolation nodes within each subgrid are distributed in a Chebyshev-type configuration [33, 34]. Besides avoiding the well-known wild oscillations of polynomial interpolants on equally-spaced nodes, the use of Chebyshev nodes is important in reducing the sensitivity of the interpolation accuracy to the errors incurred in the tabulation process. Furthermore, its use is required in ensuring the numerical stability of the barycentric Lagrange formula [35].

3.3. Singularity extraction

For r, r′ ∈ m, the propagation factor F^α(k_ρ,z,z′) of Eq. (10) contains a direct-propagation (primary) term which involves no reflection at layer interfaces. In evaluating the layered-medium DGFs, contributions associated with this term can be grouped together (here denoted by ${\overline{\mathbf{G}}}_{\mathrm{pri}}^{\gamma }(\mathbf{r},{\mathbf{r}}^{\prime })$ and written in closed form, i.e. [14],

In evaluating the secondary DGFs, ${\overline{\mathbf{G}}}_{\sec }^{\gamma }(\mathbf{r},{\mathbf{r}}^{\prime })$, singularities are also encountered when the limits r → r′ and r⁰ → I (I denotes a layer interface) are taken simultaneously. A parallel treatment of these secondary-term singularities is therefore required when the overlapping or neighboring pair of S_i and S_j are located on or very close to a layer boundary. Within the MoM representations, these singularities translate into to the singularities of the SIs at r_s = Z = 0, which are fully captured by the QS terms (see Subsection 3.2). Thus, they can be equivalently extracted from these terms. Since the secondary DGFs account for the presence of distinct layers in the background medium, they naturally possess different spatial symmetry and dependence as compared to the primary counterparts. Consequently, a different set of singularity extraction integrals is required in handling these singularities.

For ease of unification, we closely follow the notation and numbering system in [22] and refer the reader to the results therein whenever possible. Furthermore, since the types of singularities encountered are independent of which case (m = n or m ≠ n) is considered, we illustrate the treatments by considering the case of T_i,T_j ∈ m. When combined with the TI scheme discussed in Subsection 3.2, this extraction procedure amounts to adding back only the singularity-free parts of the QS terms to the interpolation results (recall that, by construction, the interpolation results contain no contribution from the QS terms). In effect, this excludes the singular parts of the QS terms from the double surface integrations, which therefore requires them to be integrated separately (via the singularity extraction integrals) and added to the matrix elements.

The singularities contained by the first and last RHS terms of Eq. (4), which are captured by the QS term of Eq. (42), can be handled via the familiar expansion of $\exp (i{k}_m\tilde{R})$ around $\tilde{R}=0$ (see e.g. [3]). As typically suggested, the first two asymmetric terms in the series are excluded from the double surface integrations over T_i and T_j. It can then be shown that the net result of this procedure is equivalent to extracting the terms

They involve the basic integrals over a triangle (with q = −1, +1),

To extract the singularities in the first two RHS terms of Eq. (11), we need to consider the QS term of Eq. (44). Since the singularity of the Hankel function is canceled by the term 1/r_s, the relevant term is $r_s/\tilde{R}(Z+\tilde{R})$, which exhibits a strongly singular behavior around r_s = Z = 0. This singular behavior can be extracted by rewriting this term into $r_s/{\tilde{R}}^2-r_{s}Z/{\tilde{R}}^2(Z+\tilde{R})$ and excluding the first term from the double surface integrations. Note that although the second term is still singular at r_s = Z =0, it is not problematic for the occurring surface integrations since its evaluation is not required when T_i and T_j lie on the same layer interface, either overlapping or neighboring. The latter case is obvious since this term vanishes when Z = 0, whereas the former is more subtle and will be discussed below. Other situations are also not problematic since the relative contribution of this term decreases with smaller values of Z and Gaussian quadratures typically only involve points on the interior of a triangle. As before, it can be shown that this procedure is equivalent to extracting

They can be evaluated semi-analytically using the formulae derived in Appendix C.

Several subtle occurrences in evaluating contributions of the secondary DGFs to the matrix elements are also worth noting. When considering a fully-overlapping pair of T_i and T_j which is parallel to layer interfaces, the first two RHS terms of Eq. (11) vanish since their integrands (of the double surface integrations) are antisymmetric under the exchange r↔r′. By the same token, when the fully-overlapping pair is not parallel to layer interfaces, the same situation still holds but only for the single-reflection terms. Thus, together with the fact that Gaussian quadratures typically only involve interior points, they imply that the singularity treatments for Eq. (11) do not require function evaluations at r_s = Z = 0. Note that similar situations also hold for the weak logarithmic singularity of Eq. (41) and the jump discontinuity of Eq. (43). For these irregularities, however, the double surface integrals involving the point r_s = Z = 0 vanish due to the terms $\widehat{\mathbf{z}}\cdot {\mathbf{f}}_i(\mathbf{r})$ and $\widehat{\mathbf{z}}\cdot {\mathbf{f}}_j({\mathbf{r}}^{\prime })$ in the integrands. Finally, the procedures outline here (and in Appendix C) can be similarly applied in obtaining additional singularity extraction integrals involved in near-field evaluations of Eq. (3), on or very close to a layer interface.

4. Numerical examples

The results presented in this section are calculated using comparatively high numbers of mesh elements, corresponding to roughly 21,000 basis functions, to ensure low discretization errors. In our C implementation, this results in Eq. (1) being solved within 2 h to 6 h on a single thread of an Intel Haswell-E chipset at 3.6 GHz. The difference in computation time is mainly due to the different number of terms in the propagation factor of Eq. (10), which, in turn, is determined by the number of interfaces delimiting the layer(s) in which the scatterer resides.

4.1. Diffraction by single holes

To validate the proposed scheme, we first treat numerically-challenging problems of diffraction by a subwavelength aperture in an optically-thick real-metal film. A schematic representation of the system is depicted in Fig. 3(a). The physical parameters are chosen based on the systems studied theoretically and experimentally in [36]. They correspond to a plane wave of wavelength λ = 660 nm impinging at normal incidence $(\mathbf{k}=|\mathbf{k}|\widehat{\mathbf{z}})$ upon a circular hole in a freely-suspended Ag film of thickness h = 300 nm. The hole diameters considered are also chosen to be within the validity region of the theoretical treatment therein, i.e., the coupled-mode method (CMM) [37] within the fundamental-mode approximation. Seen as a defect in the metal layer, the hole defines the relevant surface S (Fig. 3(b)), the top and bottom parts of which coincide with the metal-air interfaces (see Section 5 for more details). This implies that treatments of secondary-term singularities, discussed in Subsection 3.3, are strictly required. For illustration, a much coarser version of the resulting surface mesh is also shown in Fig. 3(c). The complex-valued refractive index of Ag, denoted by n_M, is taken from [38].

 

Fig. 3 (a) A schematic representation of the single-hole geometry. (b) A side view showing the corresponding scatterer surface S and its relative position. (c) A 3D view of the triangulated S, shown with a reduced number of triangles for visibility.

Download Full Size | PDF

In the CMM, the finite conductivity of the real metal is approximately accounted for through the surface impedance boundary conditions. Additionally, a phenomenological effective hole radius of the form r+qδ, where δ = λ/[2πIm(n_M)] is the skin depth, is used in its formulation to account for field penetration into the metal. By comparing the CMM predictions against rigorous VIE calculations for a hole diameter d = 200 nm and various film thicknesses (up to 150 nm), the value q = 1 was shown to be adequate [36]. The compared quantity was the normalized radial intensity $〈\mathbf{S}(\mathbf{r})〉\cdot \widehat{\mathbf{r}}/{\left|〈\mathbf{S}(\mathbf{r})〉\cdot \widehat{\mathbf{r}}\right|}_{\theta =0}$ as a function of polar θ, where both |r| (taken to be in the far-field zone, k|r|≫1) and the azimuth angle ϕ (taken to be either 0 or π/2) are kept constant. Here 〈S(r)〉 denotes the time-averaged Poynting vector.

A set of corresponding comparisons between the CMM and the SIE-MoM calculations are shown in Figs. 4(a)–4(c), among which only the hole diameter is varied. For far-field evaluation of Eq. (3), we utilize the asymptotic form of layered-medium DGFs in their angular spectrum representation [39]. Figure 4(a) confirms the findings reported in [36] for the same hole diameter. The comparatively better agreement in our case is expected since the CMM inherently assumes negligible direct tunneling across the metal film, which is better satisfied by the larger film thickness (h = 300 nm) considered here. For larger hole diameters, depicted in Figs. 4(b) and 4(c), the comparisons show that the value q = 1 is no longer adequate, while the far-field behavior is still very well described within the fundamental-mode approximation.

 

Fig. 4 Normalized intensities as functions of θ obtained from the CMM (gray lines) and the SIE-MoM (black lines) for ϕ = 0 (║) and ϕ = π/2 (⊥), for hole diameters (a) d = 200 nm, (b) 300 nm, and (c) 600 nm. The incident plane wave E_pl(r) is specified by λ = 660 nm, $\mathbf{k}=\left|\mathbf{k}\right|\widehat{\mathbf{z}}$ and |E_pl(r)| = 1, linearly polarized in the x-direction. The corresponding spatial distributions of |E(r)| at z = 10 nm are shown in (d), (e), and (f). For reference and scaling, the holes are drawn as white circles in the color plots.

Download Full Size | PDF

To justify the different choices of q and to demonstrate the near-field accuracy of the SIE-MoM scheme, the electric field distributions are evaluated at a distance of 10 nm from the exit plane of each system and plotted in Figs. 4(d)–4(f). They clearly show different depths of field penetration into the metal, which are captured by the different values of q. Moreover, the difference in near-field profiles, when compared to the fundamental waveguide mode TE₁₁ (see [36]) assumed in the CMM, does not translate to the far fields as this difference corresponds to contributions from evanescent-wave components.

4.2. Hybrid plasmonic waveguide for nanofocusing

An important advantage of the SIE and other IE-based formulations is the ease of introducing well-specified incident fields into the system. Information regarding the EM sources only enters the IEs through the quantities ${\mathbf{E}}_{\nu }^{\mathrm{inc}}(\mathbf{r})$ and ${\mathbf{H}}_{\nu }^{\mathrm{inc}}(\mathbf{r})$ (see Eqs (1) and (3)), which do not depend on the geometry (for ν = o, any property) of the scatterer. The typical procedure of surface or volume integration over the domain of the EM sources (often placed at infinity) can therefore be bypassed. This is especially beneficial when the incident fields require complicated descriptions of EM sources, e.g., in studies of surface plasmon polariton (SPP) scattering [40]. Moreover, these fields can also be obtained in conjunction with other preestablished numerical schemes.

To illustrate this point, we consider a system of hybrid plasmonic waveguide in the form of a metal-loaded silicon-on-insulator (SOI) slab waveguide [41]. A strip of Ag is patterned onto the SOI, separated by a thin spacer of silica, as depicted in Fig. 5(a). This results in bound hybrid modes in the system, exhibiting strong lateral confinement. These modes are formed through the coupling between the guided slab modes and the SPPs of the metal. By tapering the width of the metal strip (Fig. 5(b)), these bound modes can be adiabatically transformed into predominantly-plasmonic modes, allowing light nanofocusing down to extremely small mode areas. Figure 5(c) shows the taper angle θ_t limited by the adiabatic-focusing criterion and the minimum strip width a limited by current fabrication technology.

 

Fig. 5 (a) Cross-sectional geometry of the hybrid plasmonic waveguide. (b) A 3D view of the tapered metal strip. The surface mesh is shown with a reduced number of triangles for visibility. (c) A top view of the apex region.

Download Full Size | PDF

The geometry considered here closely follows one of systems studied systematically in [41], where the parameters are optimized for low-loss focusing and compatibility with silicon-based photonics. They correspond to a wavelength λ = 1,550 nm, a strip with thickness t = 50 nm and taper angle θ_t = 30°, and a slab waveguide with slab thickness h−s = 160 nm and spacer thickness s = 10 nm. We further choose the base width w and the apex width a of the strip to be 300 nm and 20 nm, respectively. The refractive index of Ag is taken from [38], and those of silicon and silica are taken to be 3.48 and 1.53, respectively. With these parameters, the four-layer slab waveguide only supports the fundamental TE and TM guided modes. We note that the performance of this system was not fully examined in [41]. It was shown that although this system exhibits superior lateral confinement, it does not sustain bound modes for a certain range of strip widths.

In the SIE-MoM scheme, ${\mathbf{E}}_{\mathrm{o}}^{\mathrm{inc}}(\mathbf{r})$ and ${\mathbf{H}}_{\mathrm{o}}^{\mathrm{inc}}(\mathbf{r})$ thus correspond to the fundamental guided slab modes. They are obtained by numerically solving the associated transcendental equations [42] of the four-layer waveguide. Moreover, we consider only the focusing of single modes propagating along the tapering axis (in the positive x-direction). Figures 6(a) and 6(b) show the near-field profiles at 5 nm above the metal strip, for TM and TE incident waves, respectively. As reported in [41], the hybrid waveguide system only supports bound TM-like modes, resulting in profoundly different nanofocusing behaviors of the two slab modes. The expected extreme-subwavelength light confinement at the apex is also clearly exhibited in Fig. 6(a).

 

Fig. 6 Spatial distribution of |E(r)| at z = 55 nm, i.e., 5 nm above the metal strip for the fundamental (a) TM and (b) TE incident waves satisfying max(|E_inc(r)|) = 1, with λ = 1,550 nm. For reference and scaling, the lateral boundaries of the metal strip are drawn in white lines. For the TM case, the same quantity is plotted as a function of z at (x, y) = (x_apex + Δ, y_apex), with Δ = 5 nm, where (c) the lateral and (d) transverse components of the field are shown along with (e) the total field. For comparison, the corresponding components of the TM incident field are shown as gray lines. The layer interfaces are indicated by solid red lines, and the transverse extent of the strip by dashed black lines. In (d), ε_r(r) denotes the relative permittivity of the background medium.

Download Full Size | PDF

The light confinement along the z direction, at 5 nm from the apex, is depicted in Figs. 6(c)–6(e), where the lateral and transverse components of the electric field are shown separately. The spatial profiles of the TM incident electric field are also plotted for comparison. Repulsion of the z-component of the electric field by the metal apex due to coupling to the SPPs can be clearly seen. Figures 6(c) and 6(d) also serve to show that the boundary conditions at layer interfaces are automatically fulfilled. Furthermore, a very strong field enhancement corresponding to max(|E(r)|²)/max(|E_inc(r)|²) > 500 is observed (Fig. 6(e)).

5. Remarks on generalization

Our discussion thus far still carries an important question regarding the appropriate mathematical treatment when the scatterer surface S overlaps with a layer interface I. That is, when a portion of S coincides (Fig. 3) or is in contact (Fig. 5) with I. While J(r) and M(r) are continuous across I, ${\overline{\mathbf{G}}}_{\mathrm{o}}(\mathbf{r},{\mathbf{r}}^{\prime })$ exhibits discontinuities when either r or r′ crosses I. This therefore seems to imply that, in evaluating Eqs. (4) and (11) for S_j ∈ I (recall that S_j is the supported area of f_j(r′)), one has to decide beforehand on the exact manner in which S→I, with different limits representing different situations. For instance, let us consider the top surface of S in Fig. 3(b), replacing the upper half-space with glass. Taking this surface to be infinitesimally above or below the metal-glass interface seems to suggest two different limiting cases, one with a slight protrusion of the hole material into the glass and the other with a thin excess of metal at the interface. On the contrary, as will be shown below, a closer look at the continuity properties of the SIEs (not only of the DGFs) reveals that the two limits are mathematically equivalent, as one would intuitively expect. The case of S_j ∈ I is therefore well-defined and it further determines which limit should be taken for S_i ∈ I.

For the discussion to follow, it is illuminating to first recall that Eq. (1) stems from two SIEs of the form

The last remark also implies that, when both S_i, S_j ∈ I, Eqs. (4) and (11) are evaluated within the limit z → z′, the manner with which (z − z′) → 0 is well-defined. This further entails that the primary term of F^α(k_ρ,z,z′) in Eq. (10), which, if present, takes the form exp(ik_m|z – z′|), has well-defined derivatives with respect to z and z′. As a consequence, the MoM manipulation (involving ∂_z and ∂_z′) performed on the secondary terms, which leads to the convenient forms of Eqs. (4) and (11), is also applicable to the primary term. The relevant continuity properties of the SIEs can thus be examined by considering Eqs. (4) and (11) (with the primary term included) and the continuity relations of F^α with respect to z′.

To this end, let us consider an interface I_m between layer m and layer m + 1, positioned at z = d_m. For a source point r′ satisfying z′ = d_m and an observation point r ∈ layer n, we can evaluate two different propagation factors, $F_{n,m}^{\alpha }(z,d_m)$ and $F_{n,m+1}^{\alpha }(z,d_m)$, where the former (latter) assumes that z′ → d_m from within layer m (m + 1). From the reciprocity relations and the boundary conditions satisfied by F^α, it can be shown that they are related by [13]

Cases with a portion of S overlapping (coinciding or in contact) with I_m can thus be treated without ambiguity, with the attached freedom in choosing how S → I_m is to be interpreted. Moreover, since the two limits for S → I_m are indistinguishable, consistency with the aforementioned meaning of r → S dictates that S_i ∈ I_m should be interpreted as S_i → I_m from within the layer in which V_i resides, or more generally, from within the layer in which the part of V_i associated with S_i resides. An S_i associated with the top or bottom surface of Fig. 3(c) is thus defined within the metal layer, while those associated with the bottom surface of Fig. 5(b), are defined within the air half-space.

In passing, we note that the equivalence of the two limits for S → I_m can be exploited for efficiency. For instance, by taking both the top and bottom parts of S in Fig. 3(c), and thus any Sj defined on them, to be within the metal layer, one avoids any evaluation of terms involving g_mn (Eq. (17)) in the matrix-filling process, along with the three-dimensional tabulation-interpolation scheme associated with it. Another benefit is that, since no S_i nor Sj straddles a layer interface, occurrences of additional line integrals in evaluating Eq. (11) are completely dispensed with. Note that, for S_i or Sj, the term straddling always implies that its common edge coincides with the interface, a situation that greatly simplifies further mathematical treatments. This also serves to show that the line integrals do not need to be invoked unless S truly straddles a layer interface, in which case, the integrals

Finally, cases with piecewise homogeneous scatterers can be treated by again referring to the validity regions of Eq. (31). The treatment parallels how multiple scatterers are handled in homogeneous media since the SIEs for such media are merely special cases of Eq. (31). To briefly illustrate this, we consider two localized scatterers with general relative positions with respect to the layered background and to one another. They are indexed by k = 1,2 and specified by (ε_k, µ_k,V_k,S_k,{m_k}), with {m_k} denoting the set of layer indices associated with scatterer k. For evaluation of Eq. (1), we thus define S = S₁ ∪ S₂ and V_i = V₁ ∪ V₂. Let us now examine the contribution of a triangle pair (T_i,T_j) to the matrix elements, where T_i (T_j) is the testing (basis) triangle. When dealing with ${\widehat{O}}_o(\mathbf{r})$ the surface integral over S stems from a volume integral over V_o. Thus, the contribution of (T_i,T_j) can be evaluated provided only that each triangle has been assigned a layer index m, requiring no special treatment when compared to a single-scatterer case. Concerning ${\widehat{O}}_{\mathrm{i}}(\mathbf{r})$however, the integral over S stems from two independent volume integrals, i.e., over V₁ and V₂. Thus, evaluating the contribution of (T_i,T_j) only requires that T_j has been assigned (ε_k, µ_k), which is chosen based on which S_k contains T_j. In other words, for T_j ∈ S_k, T_i is always seen as being outside V_k and therefore within the validity region of Eq. (31) for S = S_k. Importantly, the symmetry between the contributions of (T_i,T_j) and (T_j,T_i) is generally no longer retained when T_i and T_j belong to different S_k’s. Furthermore, ${\mathbf{E}}_{\mathrm{i}}^{\mathrm{inc}}(\mathbf{r})$ and ${\mathbf{H}}_{\mathrm{i}}^{\mathrm{inc}}(\mathbf{r})$ (if any) should be consistently evaluated, involving only the EM sources inside V_k for T_j ∈ S_k. For the special case of S₁ and S₂ being in direct contact, the continuity of J(r) and M(r) across the scatterers’ interface additionally needs to be enforced. This is typically achieved by representing the two contact surfaces with a single surface mesh, and properly double-counting the contribution of each Sj defined on it. Lastly, it should be noted that when the contact or proximity surfaces (between S₁ and S₂ or between S_k and I) are relatively sharp, one often has to take nonlocal effects into account [43]. These effects, however, cannot be readily accounted for within the framework developed here.

6. Conclusion

We have presented a robust numerical implementation of the surface integral approach to light scattering in stratified media. The limiting difficulties for efficient and well-rounded use of this method in nano-optics are addressed. In particular, those related to singular spatial integrands at layer interfaces and the highly laborious spectral-domain integrations in the matrix-filling process are mitigated. The proposed scheme is thus not limited by the relative positioning of the scatterer(s) in the layered medium, while allowing for a comparatively very high number of mesh elements without the need of invoking parallel computing. Furthermore, it can be readily built on top of a working implementation of the SIE-MoM scheme for homogeneous media, facilitated by the use of available well-tested numerical libraries and modifications thereof. Through the case examples treated here, we have also illustrated the versatility and accuracy of this method for various relevant applications in nano-optics.

Appendix A: The generalized weighted-averages (WA) method

Within the WA framework, one starts by transforming the tail integral (see Subsection 3.1), T, into an infinite sum of partial integrals according to T = u₀ + I, where $I={\displaystyle {\sum }_{n=1}^{\infty }{u}_n,}$

(37)$$u_n={\displaystyle {\int }_{{\xi }_{n-1}}^{{\xi }_n}t\left(k_{\rho }\right)\mathrm{d}{k}_{\rho }},$$

and t(k_ρ) denotes a generic Sommerfeld-type integrand. Here a set of break points, {ξ_n}, with ξ₋₁ = a is introduced to subdivide the integration range into finite subintervals (see Fig 2). The first partial integral, u₀, is evaluated using the aforementioned adaptive quadrature, where ξ₀ is chosen to be an exact zero crossing of the Bessel function satisfying ξ₀ > a. The break points with nonnegative indices are further chosen to be equidistant according to ξ_n = ξ₀ + nh, where h = π/max(r_s,Z). With these choices of break points, each of the remaining partial integrals, u_n, is adequately computed with a 16-point Gaussian quadrature rule [31].

The procedure described thus far effectively reduces the evaluation of T to computing the limit of a sequence of partial sums, i.e. computing I = lim_N→∞ I_N, with $I_N={\displaystyle {\sum }_{n=1}^N{u}_n}$. Given a finite sequence of partial sums I₁,…,I_N associated with a set of equidistant break points, the best estimation of I is obtained from the generalized WA formula [23],

(38)$$I_N^{*}=\frac{{\sum }_{n=1}^N{w}_n{I}_n}{{\sum }_{n=1}^N{w}_n},$$

where the set of generalized weights, {w_n}, is given by

(39)$$w_n={\left(-1\right)}^{n+1}\exp (\gamma {\xi }_n)\left(\begin{array}{c}N-1\\ n-1\end{array}\right){\xi }_n^{N-2-q}.$$

The asymptotic behavior of the integrand is taken into account in calculating the weights through the complex parameter γ and the real constant q. Here γ represents both the exponential decay and the oscillatory nature of the integrand in the form of exp(−γk_ρ), which implies the definition γ = Z + ir_s in our case. On the other hand, q represents the algebraic behavior of the integrand, i.e. it satisfies ${\lim }_{k_{\rho }\to \infty }[t(k_{\rho })\exp (\gamma k_{\rho })-C{k}_{\rho }^q]=0$ with C being a complex constant. The algebraic decay of the Bessel function (of the form $1/\sqrt{k_{\rho }}$) is therefore included in q. Lastly, as also shown in [31], we found that the value of N between 10 and 20 is sufficient to ensure highly accurate evaluation of the tail integral.

Appendix B: Quasi-static terms

Starting from the analytical solution of $g_m^{\mathrm{qs}}(r_s,Z)$ of Eq. (19), a set of integral-free quasi-static (QS) terms (see Subsection 3.2) can be derived to give (omitting the prefactor Ũ_m for brevity)

(40)$${\tilde{g}}_m^{\mathrm{qs}}(r_s,Z)=Z[i{\tilde{I}}_{1m}(r_s,Z)+\pi H_0^{(1)}(k_m{r}_s)/2]-e^{i{k}_m\tilde{R}}/k_m,$$

(41)$${\partial }_Z{\tilde{g}}_m^{\mathrm{qs}}(r_s,Z)=i{\tilde{I}}_{1m}(r_s,Z)+\pi H_0^{(1)}(k_m{r}_s)/2,$$

(42)$${\partial }_Z^2{\tilde{g}}_m^{\mathrm{qs}}(r_s,Z)=i{e}^{i{k}_m\tilde{R}}/\tilde{R},$$

(43)$${\partial }_{r_s}{\tilde{g}}_m^{\mathrm{qs}}(r_s,Z)=Z[i{\tilde{I}}_{2m}(r_s,Z)-\pi k_m{H}_1^{(1)}(k_m{r}_s)/2]-i{r}_s{e}^{i{k}_m\tilde{R}}/\tilde{R},$$

(44)$${\partial }_{r_s}{\partial }_Z{\tilde{g}}_m^{\mathrm{qs}}(r_s,Z)=i{\tilde{I}}_{2m}(r_s,Z)-\pi k_m{H}_1^{(1)}(k_m{r}_s)/2,$$

where $H_0^{(1)}$ and $H_1^{(1)}$ are Hankel functions of the first kind, $\tilde{R}={\left(r_s^2+Z^2\right)}^{1/2}$ is the effective distance, and ${\tilde{I}}_{1m}$ and ${\tilde{I}}_{2m}$ are the aforementioned one-dimensional integrals, which after the modifications discussed in Subsection 3.2, take the forms

(45)$${\tilde{I}}_{1m}(r_s,Z)=(1-k_m^2{r}_s^2/4)\log [(Z+\tilde{R})/r_s]-k_m^{2}Z\tilde{R}/4,$$

(46)$${\tilde{I}}_{2m}(r_s,Z)=\frac{r_s}{\tilde{R}(Z+\tilde{R})}-\frac{1}{r_s}-\frac{k_m^2{r}_s}{2}\left[\left(1-\frac{k_m^2{r}_s^2}{8}\right)\log \left(\frac{Z+\tilde{R}}{r_s}\right)-\frac{k_m^2}{8}Z\tilde{R}\right].$$

Note that when these modified versions of ${\tilde{I}}_{1m}$ and ${\tilde{I}}_{2m}$ are used, the differential relations among Eqs. (40)–(44) are no longer satisfied. Thus the ∂_Z and ${\partial }_{r_s}$ operators therein only serve to indicate the correspondence of these QS terms to the occurring SIs in Eq. (4) and Eq. (11). It is also worth noting that although more terms in Eqs. (40)–(44) can be excluded without affecting the irregularities, they are retained to minimize rounding errors during the tabulation process.

For r ∈ n, r′ ∈ m, the QS terms can be obtained by making use of the relations [44]

(47)$${\partial }_Z^2{\tilde{g}}_{mn}^{\mathrm{qs}}(r_s,Z)=i{\displaystyle {\int }_0^{\infty }{J}_0}(k_{\rho }{r}_s)e^{-k_{\rho }Z}\mathrm{d}{k}_{\rho }=i/\tilde{R},$$

(48)$${\partial }_{r_s}{\tilde{g}}_{mn}^{\mathrm{qs}}(r_s,Z)=-i{\displaystyle {\int }_0^{\infty }{J}_1}(k_{\rho }{r}_s)e^{-k_{\rho }Z}\mathrm{d}{k}_{\rho }/k_{\rho }=i(Z-\tilde{R})/r_s,$$

where the prefactor ${\tilde{U}}_{mn}=\frac{i}{4\pi }{\lim }_{k_{\rho }\to \infty }{U}_{mn}(k_{\rho })$ has been omitted and Z denotes Z_m(z′) + Z_n(z). While higher-order terms are obtained through spatial derivations, Eq. (47) also implies that

(49)$${\partial }_Z{\tilde{g}}_{mn}^{\mathrm{qs}}(r_s,Z)=i\log (Z+\tilde{R}),$$

the integral representation of which is not available. Its evaluation at r_s = Z = 0, however, can be mediated by the integral representation of Eq. (41). We note that only the derivatives of ${\tilde{g}}_m^{\mathrm{qs}}$ and ${\tilde{g}}_{mn}^{\mathrm{qs}}$ are strictly required in the TI scheme since g_m and g_mn do not contain irregularities.

Appendix C: Formulae for surface integrals K₈ and K₉

The surface integral K₈(r;T) of Eq. (29) can be evaluated by decomposition into parallel and perpendicular components with respect to the plane of triangle T, i.e.

(50)$${\mathbf{K}}_8^{\parallel }(\mathbf{r};T)={\displaystyle \sum _{k=1}^3{\mathbf{m}}_k{I}_{\log }^L(\mathbf{r};\partial T_k)},$$

(51)$${\mathbf{K}}_8^{\mathrm{n}}(\mathbf{r};T)=-\mathbf{n}(T)h{I}_{-2}^S(\mathbf{r};T).$$

By Gauss’ theorem, the in-plane component, ${\mathbf{K}}_8^{\parallel }$, has been rewritten as a line integral over the three edges, {∂T_k}, of T. The resulting expression naturally involves the outward-directed unit normal vectors, {m_k}, of {∂T_k} (see Fig. 1(b)) and the logarithmic integral over a line segment,

(52)$$I_{\log }^L(\mathbf{r};\mathrm{\Delta }L)={\displaystyle {\int }_{\mathrm{\Delta }L}\log (|\mathbf{r}-{\mathbf{r}}^{\prime }|)}\mathrm{d}{l}^{\prime }.$$

To obtain a closed-form formula for $I_{\log }^L$, we consider a line segment, ΔL, between two points p₁ and p₂. A sketch of the situation and the relevant parameters are depicted in Fig. 7. Here we closely follow the notations in [22]. Given the free parameters r, p₁, and p₂, the following parameters can be immediately deduced : s = (p₂ −p₁)/|p₂ −p₁|, R⁺ = |p2 −r|, R⁻ = |p₁ −r|, s⁺ = s · (p₂ −r), and s⁻ = s · (p₁ −r). They also imply that s⁺ > s⁻ and $R_0^2={\left(R^{+}\right)}^2-{\left(s^{+}\right)}^2={\left(R^{-}\right)}^2-{\left(s^{-}\right)}^2$. The line integral can then be expressed in terms of these parameters by writing

(53)$$I_{\log }^L(\mathbf{r};\mathrm{\Delta }L)=\frac{1}{2}{\displaystyle {\int }_{s^{-}}^{s^{+}}\log (R_0^2+s^2)\mathrm{d}s},$$

and solving the RHS to give

(54)$$\{\begin{array}{ll}{s}^{+}\log (R^{+})-s^{-}\log (R^{-})-\mathrm{\Delta }L+R_0\left[\arctan \left(\frac{s^{+}}{R_0}\right)-\arctan \left(\frac{s^{-}}{R_0}\right)\right]:\hfill & R_0\ne 0\hfill \\ s^{+}\log (|s^{+}|)-s^{-}\log (|s^{-}|)-\mathrm{\Delta }L:\hfill & R_0=0\hfill \end{array}.$$

 

Fig. 7 A line segment bounded by p₁ and p₂, along which r′ is located and the integrand is evaluated.

Download Full Size | PDF

In Eq. (51)n(T) denotes the unit normal vector of T and h = n(T)·(r − r′) is the plane-to-point normal displacement which remains constant for all r′ in the plane of T. The surface integral $I_{-2}^S(\mathbf{r};T)$ is defined according to Eq. (26) with q = −2. Although a closed-form expression for this integral is not readily available, it can be efficiently evaluated using high-order quadrature rules [45]. Moreover, its evaluation in difficult cases where T_i and/or T_j lie on a layer interface can be avoided. As mentioned in Subsection 3.3, the terms involving $I_{-2}^S(\mathbf{r};T)$ do not need to be evaluated when the triangle pair fully overlap on or very close to a layer interface. For a neighboring pair, the contribution associated with it also vanishes if T_j is parallel to layer interfaces since n(T_j) is (anti)-parallel to $\widehat{\mathbf{z}}$. Moreover, when only T_i is parallel to interfaces, its evaluations can be dispensed with by reversing the roles of T_i and T_j (as testing and basis triangles) and making use of the reciprocity relations satisfied by the DGFs [14]. Finally, by using of the relation ∇′ f (|r − r′|) × (r − r′) = 0 for a differentiable function f, the integral K₉ of Eq. (30) can be rewritten as

(55)$${\mathbf{K}}_9(\mathbf{r};T^{\pm })=\pm \frac{L}{2A^{\pm }}{\mathbf{K}}_8(\mathbf{r};T^{\pm })\times (\mathbf{r}-{\mathbf{p}}^{\pm }),$$

where L, A^±, and p^± are associated with f(r′) (see Fig. 1(b)).

Acknowledgments

This work was supported by VIDI grant from the Netherlands Organization for Scientific Research (NWO).

References and links

1. R. F. Harrington, “Boundary integral formulations for homogeneous material bodies,” J. Electromag. Waves Appl. 3, 1–15 (1989). [CrossRef]  

2. R. F. Harrington, Field Computation by Moment Methods (R. E. Krieger, 1968).

3. 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]  

4. 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]  

5. 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]  

6. M. G. Araújo, J. M. Taboada, D. M. Sols, 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]  

7. D. M. Solís, J. M. Taboada, F. Obelleiro, L. M. Liz-Marzán, and F. J. García de Abajo, “Toward ultimate nanoplas-monics modeling,” ACS Nano 8, 7559–7570 (2014). [CrossRef]  

8. R. Rodríguez-Oliveros and J. A. Sánchez-Gil, “Localized surface-plasmon resonances on single and coupled nanoparticles through surface integral equations for flexible surfaces,” Opt. Express 19, 12208–12219 (2011). [CrossRef]   [PubMed]  

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

10. M. Paulus and O. J. F. Martin, “Light propagation and scattering in stratified media: a Green’s tensor approach,” J. Opt. Soc. Am. A 18, 854–861 (2001). [CrossRef]  

11. J. Alegret, P. Johansson, and M. Käll, “Green’s tensor calculations of plasmon resonances of single holes and hole pairs in thin gold films,” New J. Phys. 10, 105004 (2008). [CrossRef]  

12. K. A. Michalski and J. R. Mosig, “Multilayered media green’s functions in integral equation formulations,” IEEE Trans. Antennas Propag. 45, 508–519 (1997). [CrossRef]  

13. Y. P. Chen, W. C. Chew, and L. Jiang, “A new Green’s function formulation for modeling homogeneous objects in layered medium,” IEEE Trans. Antennas Propag. 60, 4766–4776 (2012). [CrossRef]  

14. W. C. Chew, Waves and Fields in Inhomogeneous Media (Van Nostrand Reinhold, 1990).

15. Y. P. Chen, W. E. I. Sha, W. C. H. Choy, L. Jiang, and W. C. Chew, “Study on spontaneous emission in complex multilayered plasmonic system via surface integral equation approach with layered medium Green’s function,” Opt. Express 20, 20210–20221 (2012). [CrossRef]   [PubMed]  

16. Y. P. Chen, W. E. I. Sha, L. Jiang, and J. Hu, “Graphene plasmonics for tuning photon decay rate near metallic split-ring resonator in a multilayered substrate,” Opt. Express 23, 2798–2807 (2015). [CrossRef]   [PubMed]  

17. Y. P. Chen, W. C. Chew, and L. Jiang, “A novel implementation of discrete complex image method for layered medium Green’s function,” IEEE Antennas Wirel. Propag. Lett. 10, 419–422 (2011). [CrossRef]  

18. A. Alparslan, M. I. Aksun, and K. A. Michalski, “Closed-form green’s functions in planar layered media for all ranges and materials,” IEEE Trans. Microwave Theory Tech. 58, 602–613 (2010). [CrossRef]  

19. E. P. Karabulut, A. T. Erdogan, and M. I. Aksun, “Discrete complex image method with spatial error criterion,” IEEE Trans. Microwave Theory Tech. 59, 793–802 (2011). [CrossRef]  

20. R. D. Graglia, “On the numerical integration of the linear shape functions times the 3-D Green’s function or its gradient on a plane triangle,” IEEE Trans. Antennas Propag. 41, 1448–1455 (1993). [CrossRef]  

21. P. Yla-Oijala and M. Taskinen, “Calculation of CFIE impedance matrix elements with RWG and n × RWG functions,” IEEE Trans. Antennas Propag. 51, 1837–1846 (2003). [CrossRef]  

22. I. Hanninen, M. Taskinen, and J. Sarvas, “Singularity subtraction integral formulae for surface integral equations with RWG, rooftop and hybrid basis functions,” Prog. Electromag. Res. 63, 243–278 (2006). [CrossRef]  

23. J. Mosig, “The weighted averages algorithm revisited,” IEEE Trans. Antennas Propag. 60, 2011–2018 (2012). [CrossRef]  

24. L. N. Medgyesi-Mitschang, J. M. Putnam, and M. B. Gedera, “Generalized method of moments for three-dimensional penetrable scatterers,” J. Opt. Soc. Am. A 11, 1383–1398 (1994). [CrossRef]  

25. S. M. Rao, D. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antennas Propag. 30, 409–418 (1982). [CrossRef]  

26. W. C. Chew, J. L. Xiong, and M. A. Saville, “A matrix-friendly formulation of layered medium Green’s function,” IEEE Antennas Wireless Propag. Lett. 5, 490–494 (2006). [CrossRef]  

27. W. C. Chew and S.-Y. Chen, “Response of a point source embedded in a layered medium,” IEEE Antennas Wireless Propag. Lett. 2, 254–258 (2003).

28. A. Sommerfeld, Partial Differential Equations in Physics (Academic, 1949).

29. M. Paulus, P. Gay-Balmaz, and O. J. F. Martin, “Accurate and efficient computation of the Green’s tensor for stratified media,” Phys. Rev. E 62, 5797–5807 (2000). [CrossRef]  

30. R. Piessens, E. de Doncker-Kapenga, C. W. Überhuber, and D. K. Kahaner, Quadpack : A Subroutine Package for Automatic Integration (Springer-Verlag, 1983). [CrossRef]  

31. R. Golubovic, A. G. Polimeridis, and J. R. Mosig, “Efficient algorithms for computing Sommerfeld integral tails,” IEEE Trans. Antennas Propag. 60, 2409–2417 (2012). [CrossRef]  

32. P. R. Atkins and W. C. Chew, “Fast computation of the dyadic Green’s function for layered media via interpolation,” IEEE Antennas Wireless Propag. Lett. 9, 493–496 (2010). [CrossRef]  

33. C. W. Ueberhuber, Numerical Computation 1: Methods, Software, and Analysis (Springer-Verlag, 1997).

34. J.-P. Berrut and L. N. Trefethen, “Barycentric Lagrange interpolation,” SIAM Review 46, 501–517 (2004). [CrossRef]  

35. N. J. Higham, “The numerical stability of barycentric Lagrange interpolation,” IMA J. Numer. Anal. 24, 547–556 (2004). [CrossRef]  

36. J.-M. Yi, A. Cuche, F. de León-Pérez, A. Degiron, E. Laux, E. Devaux, C. Genet, J. Alegret, L. Martín-Moreno, and T. W. Ebbesen, “Diffraction regimes of single holes,” Phys. Rev. Lett. 109, 023901 (2012). [CrossRef]   [PubMed]  

37. F. de León-Pérez, G. Brucoli, F. J. García-Vidal, and L. Martín-Moreno, “Theory on the scattering of light and surface plasmon polaritons by arrays of holes and dimples in a metal film,” New J. Phys. 10, 105017 (2008). [CrossRef]  

38. P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6, 4370–4379 (1972). [CrossRef]  

39. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University, 2006). [CrossRef]  

40. I. Chremmos, “Magnetic field integral equation analysis of surface plasmon scattering by rectangular dielectric channel discontinuities,” J. Opt. Soc. Am. A 27, 85–94 (2010). [CrossRef]  

41. L. Lafone, T. P. H. Sidiropoulos, and R. F. Oulton, “Silicon-based metal-loaded plasmonic waveguides for low-loss nanofocusing,” Opt. Lett. 39, 4356–4359 (2014). [CrossRef]   [PubMed]  

42. K. Okamoto, Fundamentals of Optical Waveguides (Academic, 2006).

43. F. J. García de Abajo, “Nonlocal effects in the plasmons of strongly interacting nanoparticles, dimers, and waveguides,” J. Phys. Chem. C 112, 17983–17987 (2008). [CrossRef]  

44. I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products (Elsevier Academic, 2007).

45. L. Zhang, T. Cui, and H. Liu, “A set of symmetric quadrature rules on triangles and tetrahedra,” J. Comput. Math. 27, 89–96 (2009).