Transient analysis of electromagnetic wave interactions on plasmonic nanostructures using a surface integral equation solver

Ismail E. Uysal; H. Arda Ülkü; Hakan Bağci

doi:10.1364/JOSAA.33.001747

1. INTRODUCTION

Metallic nanostructures support surface plasmon modes upon illumination by electromagnetic fields at optical frequencies [1]. These modes significantly enhance scattering from nanostructures in the near- and far-field regions [1], enabling their use in various applications, such as the design of nanoantennas for biochemical sensing [2], subwavelength waveguides for telecommunication [3], and ring resonators for integrated compact optical components [4]. Characteristics of these plasmon modes, such as resonance frequency, propagation distance, and amplitude decay rate, depend on the nanostructure’s material properties, shape, and size as well as the material properties of the background medium where the nanostructure resides [1]. Therefore, simulation tools, which are capable of characterizing plasmonic field interactions on arbitrarily shaped geometries with certain material properties, have to be utilized to design a plasmonic nanostructure with desired field patterns in the near- and far-field regions.

Simulation tools implemented to analyze plasmonic field interactions make use of finite-difference time-domain (FDTD) schemes [5–7], finite element method (FEM) [8], discrete dipole approximation (DDA) [9], or surface integral equation (SIE) solvers [10–14]. Among these methods, FDTD is the most frequently used one due to its simplicity in formulation and implementation. However, FDTD (like any other differential equation-based solver such as FEM) requires volumetric discretizations of the plasmonic nanostructure and the surrounding medium. This becomes inefficient considering the fact that exponentially decaying behavior of the plasmonic fields could only be accurately captured using very small discretization elements, especially in the vicinity of metal/dielectric interfaces. Also, the unbounded surrounding medium has to be truncated into a (bounded) computational domain [15]. This is achieved by using absorbing boundary conditions, which approximate field behavior at infinity. This approximation introduces significant errors, especially while analyzing resonant phenomena requiring long durations of simulations [16]. Additionally, FDTD suffers from numerical dispersion introduced by finite difference approximations used for discretizing time and space derivatives. This produces additional error in the solution for electrically large objects and possibly leads to unphysical effects [17].

On the contrary, SIE solvers [10–14] only discretize the surfaces of the nanostructures (the dimension of the discretization domain is reduced to two from three—in comparison with FDTD and FEM) and the solution automatically satisfies the radiation condition (the field behavior at infinity is accurately captured without the need for approximate absorbing boundary conditions). However, the use of SIE solvers for characterizing plasmonic field interactions has been limited to the frequency domain. Consequently, these solvers cannot be used when material properties are functions of fields, such Kerr nonlinear materials [18,19], and they have to be executed many times (one execution per frequency) to obtain broadband results, which are typically needed for designing plasmonic devices/systems expected to operate over the support of multiple resonances.

These disadvantages can be overcome through formulation and implementation of SIE solvers for characterizing transient plasmonic field interactions. Time-domain SIE (TD-SIE) solvers have long been used for analyzing electromagnetic field interactions on perfect electrically conducting scatterers and dielectric objects with nondispersive material properties [20–31]. However, despite their advantages, TD-SIE solvers have not been used for analyzing scattering from objects with dispersive material properties, with only one exception [32], which makes use of finite-difference delay modeling (FDDM) or convolution quadrature (CQ) techniques [33]. This can be attributed to the following challenges: (i) For dispersive media, the dielectric permittivity does not have a closed-form expression in the time domain unless one of the simplified models, such as Drude or Drude–Lorentz [19], is used. Consequently, the time-domain Green function also does not have a closed-form expression. (ii) Because both the permittivity and the Green function are dispersive, scattered field computations require evaluation of additional convolutions between the permittivity, the Green function, and the equivalent currents. Schemes for discretizing these additional convolutions have never been formulated or implemented. (iii) The dispersive Green function has an infinite tail in the time domain, which significantly increases the computational cost of the convolutions mentioned in (ii) [26–30]. But this cost can be significantly reduced by using the recently developed plane-wave time-domain (PWTD) method [26] or fast Fourier transform (FFT)-based schemes [27–30,34].

The TD-SIE solvers [31,32], which make use of FDDM or CQ techniques [33], formulate and discretize the SIEs in the Laplace domain, approximate the Laplace domain parameter $s$ in the $z$ -domain using a finite-difference scheme, and finally convert the resulting system into the time domain through the inverse $z$ -transform. Since the SIE is constructed in the Laplace domain, these TD-SIE solvers do not face the challenges (i) and (ii) listed above. But they still require a fast method to efficiently compute (discretized) temporal convolutions after conversion into the time domain. This could similarly be remedied by FFT-based schemes [27–30,34].

In this work, a different approach is used. The TD-SIE is constructed directly in the time domain and temporal samples of the dispersive Green function and permittivity are obtained using separate inverse Fourier transforms applied to their frequency-domain representations. More specifically, the well-known frequency-domain Poggio-Miller-Chan-Harrington-Wu-Tsai surface integral equation (PMCHWT-SIE) [35,36] is extended in the time domain to account for scatterers with dispersive material properties, more specifically plasmonic structures. The time-domain PMCHWT-SIE (TD-PMCHWT-SIE) is solved using the well-known marching on-in-time (MOT) scheme [20–30]. The MOT scheme expands the electric and magnetic current densities, which are introduced on the surfaces of the nanostructures, using Rao-Wilton-Glisson (RWG) [37] and polynomial basis functions [21,28] in space and time, respectively. Inserting this expansion into the TD-PMCHWT-SIE and Galerkin testing the resulting equation at discrete times yield a system of equations that is solved for the current expansion coefficients by marching in time. The additional convolutions are converted into a double discrete summation, where the inner summation corresponds to the discretized convolution of the Green function and temporal basis function [38] [addressing the challenge (ii) above]. The temporal samples of the time-domain dielectric permittivity and Green function, which are required by this MOT discretization procedure, are obtained numerically from their frequency-domain samples [addressing the challenge (i) above]. This is achieved by representing the frequency-domain Green function and permittivity in terms of a summation of weighted rational functions. The weighting coefficients are found by applying the fast relaxed vector fitting (FRVF) scheme to the frequency-domain samples [39–41]. Time-domain functions are then obtained by analytically computing the inverse Fourier transform of the summation. The resulting MOT-PMCHWT solver is used in characterizing transient field interactions on several different geometries. Accuracy of the solver is demonstrated via comparison of the results with those obtained by frequency-domain solvers or analytical solutions obtained using the Mie series expansion.

2. FORMULATION

A. Time-Domain PMCHWT-SIE

Let $V_{1}$ denote the volume of the plasmonic nanostructure with time-dependent permittivity $ε_{1} (t)$ and constant permeability $μ_{1}$ . Assume that the nanostructure resides in a nondispersive unbounded background medium denoted by $V_{0}$ . Constant permittivity and permeability in $V_{0}$ are represented by $ε_{0}$ and $μ_{0}$ , respectively. Let $S$ denote the surface enclosing $V_{1}$ . An electromagnetic wave with electric and magnetic field intensities, $E^{inc} (r, t)$ and $H^{inc} (r, t)$ , excites the nanostructure. It is assumed that $E^{inc} (r, t) \approx 0$ and $H^{inc} (r, t) \approx 0$ , $\forall r \in V_{1}$ and $t < 0$ , and they are essentially band limited to $f_{\max}$ . Using the surface equivalence principle [42] and the boundary conditions on $S$ , the time derivative of the TD-PMCHWT-SIE is constructed as [25]

\hat{n} (r) \times \partial_{t} E^{inc} (r, t) = - \hat{n} (r) \times [\partial_{t} E_{0}^{sca} (r, t) - \partial_{t} E_{1}^{sca} (r, t)],

\hat{n} (r) \times \partial_{t} H^{inc} (r, t) = - \hat{n} (r) \times [\partial_{t} H_{0}^{sca} (r, t) - \partial_{t} H_{1}^{sca} (r, t)] .

Here,

\hat{n} (r)

is the unit normal vector at

r \in S

, which is pointing toward

V_{0}

, and

E_{p}^{sca} (r, t)

and

H_{p}^{sca} (r, t)

are the scattered electric and magnetic field intensities in

V_{p}

,

p \in {0, 1}

.

E_{p}^{sca} (r, t)

and

H_{p}^{sca} (r, t)

are expressed in terms of equivalent electric and magnetic current densities,

J (r, t)

and

M (r, t)

, which are introduced on

S

using the surface equivalence principle:

\partial_{t} E_{0}^{sca} (r, t) = - L_{0} {μ_{0} J (r, t)} + Q_{0} {ε_{0}^{- 1} J (r, t)} - K_{0} {M (r, t)},

\partial_{t} H_{0}^{sca} (r, t) = - L_{0} {ε_{0} M (r, t)} + Q_{0} {μ_{0}^{- 1} M (r, t)} + K_{0} {J (r, t)},

\partial_{t} E_{1}^{sca} (r, t) = L_{1} {μ_{1} J (r, t)} - Q_{1} {{\bar{ε}}_{1} (t) * J (r, t)} + K_{1} {M (r, t)},

\partial_{t} H_{1}^{sca} (r, t) = L_{1} {ε_{1} (t) * M (r, t)} - Q_{1} {μ_{1}^{- 1} M (r, t)} - K_{1} {J (r, t)} .

In Eqs. (3)–(6), the integral operators

L_{p} {\cdot}

,

Q_{p} {\cdot}

, and

K_{p} {\cdot}

are given by

L_{p} {X (r, t)} = \int_{S} g_{p} (R, t) * \partial_{t}^{2} X (r^{'}, t) d s^{'},

Q_{p} {X (r, t)} = \nabla \int_{S} g_{p} (R, t) * \nabla^{'} \cdot X (r^{'}, t) d s^{'},

K_{p} {X (r, t)} = \nabla \times \int_{S} g_{p} (R, t) * \partial_{t} X (r^{'}, t) d s^{'} .

Here,

g_{p} (R, t)

is the Green function of the unbounded medium that has the same permittivity and permeability as

V_{p}

,

R = | r - r^{'} |

is the distance between source and observation points,

r^{'}

and

r

, “

*

” denotes temporal convolution, and

{\bar{ε}}_{1} (t)

is the time-domain inverse permittivity in

V_{1}

.

B. Time-Domain Permittivity and Inverse Permittivity

Time-domain permittivity and inverse permittivity of the plasmonic nanostructure, $ε_{1} (t)$ and ${\bar{ε}}_{1} (t)$ , are obtained by inverse Fourier transforming their frequency-domain counterparts, $ε_{1} (ω)$ and $1 / ε_{1} (ω)$ , respectively:

ε_{1} (t) = F^{- 1} {ε_{1} (ω)},

{\bar{ε}}_{1} (t) = F^{- 1} {{\bar{ε}}_{1} (ω)} = F^{- 1} {1 / ε_{1} (ω)} .

Here, it is assumed that frequency samples of

ε_{1} (ω)

and

{\bar{ε}}_{1} (ω)

, i.e.,

ε_{1} (ω_{i})

and

{\bar{ε}}_{1} (ω_{i})

,

i = 1, \dots, N_{ω}^{VF}

, where

ω_{i} \in 2 π [f_{\min}^{VF}, f_{\max}^{VF}]

, are available from tabulated experimental data, such as those provided for silver and gold in [43] and [44]. Using the FRVF scheme [39–41] (see Appendix A) on

ε_{1} (ω_{i})

and

{\bar{ε}}_{1} (ω_{i})

, one can approximate

ε_{1} (ω)

and

{\bar{ε}}_{1} (ω)

as

ε_{1} (ω) - ε_{1}^{\infty} \approx \sum_{m = 1}^{N_{p}^{VF}} \frac{b_{m}^{ϵ}}{j ω + a_{m}^{ϵ}},

{\bar{ε}}_{1} (ω) - \frac{1}{ε_{1}^{\infty}} \approx \sum_{m = 1}^{N_{p}^{VF}} \frac{b_{m}^{\bar{ϵ}}}{j ω + a_{m}^{\bar{ϵ}}},

where

j = \sqrt{- 1}

is the imaginary number,

ε_{1}^{\infty} = 1

is the high frequency limit of

ε_{1} (ω)

as

ω \to \infty

, and

{a_{m}^{ϵ}, b_{m}^{ϵ}}

and

{a_{m}^{\bar{ϵ}}, b_{m}^{\bar{ϵ}}}

are the pole/residue pairs obtained by applying the FRVF scheme to

ε_{1} (ω_{i})

and

{\bar{ε}}_{1} (ω_{i})

. Inverse Fourier transforming Eqs. (12) and (13) yields

ε_{1} (t) \approx ε_{1}^{\infty} δ (t) + \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{ϵ} u (t) e^{- a_{m}^{ϵ} t} = ε_{1}^{\infty} δ (t) + γ_{1} (t),

{\bar{ε}}_{1} (t) \approx \frac{1}{ε_{1}^{\infty}} δ (t) + \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{\bar{ϵ}} u (t) e^{- a_{m}^{\bar{ϵ}} t} = \frac{1}{ε_{1}^{\infty}} δ (t) + {\bar{γ}}_{1} (t),

where

δ (\cdot)

and

u (\cdot)

are the Dirac delta and the unit step functions, respectively. It should be noted here that the causality of

ε_{1} (t)

and

{\bar{ε}}_{1} (t)

is ensured by enforcing

a_{m}^{ϵ} > 0

and

a_{m}^{\bar{ϵ}} > 0

during the execution of the FRVF scheme.

C. Time-Domain Green Function

The integral operators $L_{p} {\cdot}$ , $Q_{p} {\cdot}$ , and $K_{p} {\cdot}$ present in Eqs. (3)–(6) and given by Eqs. (7)–(9) require $g_{p} (R, t)$ and its derivative with respect to $R$ , $\partial_{R} g_{p} (R, t)$ [due to the presence of the curl operation in Eq. (9)], to be known. For $p = 0$ (in the background medium), since $ε_{0}$ and $μ_{0}$ are simply constants, $g_{0} (R, t) = δ (t - R / c_{0}) / (4 π R)$ and $\partial_{R} g_{0} (R, t) = - δ^{'} (t - R / c_{0}) / (4 π c_{0} R) + δ (t - R / c_{0}) / (4 π R^{2})$ [45]. Here, $c_{0} = 1 / \sqrt{ε_{0} μ_{0}}$ is the speed of light in the background medium and $δ^{'} (\cdot)$ is the derivative of the Dirac delta function with respect to its argument. On the other hand, for $p = 1$ , general analytical expressions for $g_{1} (R, t)$ and $\partial_{R} g_{1} (R, t)$ are not immediately available due to the dependence of $ε_{1} (t)$ on time. This can be explained by the fact that, for the dispersive $ε_{1} (ω)$ , the inverse Fourier transform of the frequency-domain Green function $g_{1} (R, ω) = e^{- j k_{1} R} / (4 π R)$ , may not exist in closed-form. Note that, here, $k_{1} = ω \sqrt{ε_{1} (ω) μ_{1}}$ is the complex wave number in $V_{1}$ . To overcome this problem and obtain expressions for $g_{1} (R, t)$ and $\partial_{R} g_{1} (R, t)$ , the FRVF scheme is used in the same way as it is used to obtain the expressions of $ε_{1} (t)$ and ${\bar{ε}}_{1} (t)$ in Section 2.B. However, a naive application of the FRVF scheme to the samples of $g_{1} (R, ω)$ and $\partial_{R} g_{1} (R, ω)$ does not yield accurate representations for $g_{1} (R, t)$ and $\partial_{R} g_{1} (R, t)$ . Singular term $1 / (4 π R)$ (corresponding to a Dirac delta function in time) and phase $e^{j ω R / c_{1}}$ , $c_{1} = 1 / \sqrt{ε_{1}^{\infty} μ_{1}}$ (corresponding to a delay in time) are first extracted from $g_{1} (R, ω)$ . Then, the FRVF scheme (see Appendix A) is applied to the samples of the remaining function at $ω_{i} \in 2 π [f_{\min}^{VF}, f_{\max}^{VF}]$ , $i = 1, \dots, N_{ω}^{VF}$ to yield

e^{j ω R / c_{1}} [g_{1} (R, ω) - \frac{e^{- j ω R / c_{1}}}{4 π R}] = \frac{e^{- j (k_{1} R - ω R / c_{1})} - 1}{4 π R} \approx d^{g} (R) + \sum_{m = 1}^{N_{p}^{VF}} \frac{b_{m}^{g} (R)}{j ω + a_{m}^{g} (R)} .

Here,

{a_{m}^{g} (R), b_{m}^{g} (R)}

are the pole/residue pairs and

d^{g} (R)

is the constant term (with respect to frequency dependence) generated by the FRVF scheme. Note that

a_{m}^{g} (R)

,

b_{m}^{g} (R)

, and

d^{g} (R)

are functions of

R

. Using Eq. (16), one can write

g_{1} (R, ω) \approx e^{- j ω R / c_{1}} [\frac{1}{4 π R} + d^{g} (R) + \sum_{m = 1}^{N_{p}^{VF}} \frac{b_{m}^{g} (R)}{j ω + a_{m}^{g} (R)}] .

Inverse Fourier transforming Eq. (17) yields

g_{1} (R, t) \approx \frac{δ (τ_{1})}{4 π R} + d^{g} (R) δ (τ_{1}) + \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{g} (R) u (τ_{1}) e^{- a_{m}^{g} (R) τ_{1}},

where

τ_{1} = t - R / c_{1}

is the delay representing the time retardation due to the finite speed of light. A similar method is used to obtain an expression for

\partial_{R} g_{1} (R, t)

from the samples of

\partial_{R} g_{1} (R, ω)

. Note that

\partial_{R} g_{1} (R, ω) = - e^{- j k_{1} R} / (4 π R) [j k_{1} + 1 / R]

. Extracting the singularity and the phase, and applying the FRVF scheme (see Appendix A) to the samples of the remaining function at

ω_{i} \in 2 π [f_{\min}^{VF}, f_{\max}^{VF}]

,

i = 1, \dots, N_{ω}^{VF}

yield

e^{j ω R / c_{1}} [\partial_{R} g_{1} (R, ω) + (\frac{j ω}{c_{1} R} + \frac{1}{R^{2}}) \frac{e^{- j ω R / c_{1}}}{4 π}] = - \frac{e^{- j (k_{1} R - ω R / c_{1})}}{4 π} (\frac{j k_{1}}{R} + \frac{1}{R^{2}}) + \frac{1}{4 π R^{2}} + \frac{j ω}{4 π c_{1} R} \approx d^{g_{R}} (R) + j ω f^{g_{R}} (R) + \sum_{m = 1}^{N_{p}^{VF}} \frac{b_{m}^{g_{R}} (R)}{j ω + a_{m}^{g_{R}} (R)} .

Here,

{a_{m}^{g_{R}} (R), b_{m}^{g_{R}} (R)}

are the pole/residue pairs,

d^{g_{R}} (R)

is the constant term (with respect to frequency dependence), and

f^{g_{R}} (R)

is the coefficient of the linear term generated by the FRVF scheme. Using Eq. (19), one can write

\partial_{R} g_{1} (R, ω) \approx e^{- j ω R / c_{1}} [j ω (f^{g_{R}} (R) - \frac{1}{4 π c_{1} R}) + d^{g_{R}} (R) - \frac{1}{4 π R^{2}} + \sum_{m = 1}^{N_{p}^{VF}} \frac{b_{m}^{g_{R}} (R)}{j ω + a_{m}^{g_{R}} (R)}] .

Inverse Fourier transforming Eq. (20) yields

\partial_{R} g_{1} (R, t) \approx δ^{'} (τ_{1}) (f^{g_{R}} (R) - \frac{1}{4 π c_{1} R}) + δ (τ_{1}) (d^{g_{R}} (R) - \frac{1}{4 π R^{2}}) + \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{g_{R}} (R) u (τ_{1}) e^{- a_{m}^{g_{R}} (R) τ_{1}} .

One could assume

\partial_{R} g_{1} (R, t)

could be directly obtained by taking the derivative of Eq. (18) with respect to

R

, but this is not possible since the coefficients

a_{m}^{g_{R}} (R)

,

b_{m}^{g_{R}} (R)

,

d^{g_{R}} (R)

, and

f^{g_{R}} (R)

depend on

R

and are obtained numerically using the FRVF (i.e., they do not have closed-form expressions). Also, the causality of

g_{1} (R, t)

and

\partial_{R} g_{1} (R, t)

is ensured by enforcing

a_{m}^{g} (R) > 0

and

a_{m}^{g_{R}} (R) > 0

during the execution of the FRVF scheme.

D. Discretization and MOT Scheme

To numerically solve the TD-PMCHWT-SIE in Eqs. (1) and (2), first $S$ is discretized into triangular patches. Then $J (r, t)$ and $M (r, t)$ are expanded in space and time as

J (r, t) = \sum_{j^{'} = 1}^{N_{t}} \sum_{n^{'} = 1}^{N_{s}} J_{n^{'}} (j^{'} Δ t) f_{n^{'}} (r) T_{j^{'}} (t) = \sum_{j^{'} = 1}^{N_{t}} \sum_{n^{'} = 1}^{N_{s}} {{\bar{I}}_{j^{'}}^{J}}_{n^{'}} f_{n^{'}} (r) T_{j^{'}} (t),

M (r, t) = \sum_{j^{'} = 1}^{N_{t}} \sum_{n^{'} = 1}^{N_{s}} M_{n^{'}} (j^{'} Δ t) f_{n^{'}} (r) T_{j^{'}} (t) = \sum_{j^{'} = 1}^{N_{t}} \sum_{n^{'} = 1}^{N_{s}} {{\bar{I}}_{j^{'}}^{M}}_{n^{'}} f_{n^{'}} (r) T_{j^{'}} (t) .

Here,

f_{n^{'}} (r)

represent the well-known RWG functions [37],

T_{j^{'}} (t) = T (t - j^{'} Δ t)

are the shifted Lagrange interpolation functions [21,28],

Δ t

is the time step size,

{{\bar{I}}_{j^{'}}^{J}}_{n^{'}} = J_{n^{'}} (j^{'} Δ t)

and

{{\bar{I}}_{j^{'}}^{M}}_{n^{'}} = M_{n^{'}} (j^{'} Δ t)

are the unknown coefficients, and

N_{s}

and

N_{t}

are the total numbers of spatial basis functions on

S

and the time steps, respectively. It should be mentioned here that each one of

f_{n^{'}} (r)

is defined on a pair of triangular patches with support

S_{n^{'}}

and the Lagrange interpolation functions

T (t)

are piece-wise continuous polynomial functions defined on support

[- Δ t, d Δ t]

, where

d

is the order of

T (t)

. Inserting Eqs. (22) and (23) into Eqs. (1) and (2) and testing the resulting equations with

\hat{n} (r) \times f_{n} (r) δ (t - j Δ t)

,

n = 1, \dots, N_{s}

, yield a system of equations:

[\begin{matrix} {\bar{\bar{Z}}}_{0}^{JJ} & {\bar{\bar{Z}}}_{0}^{JM} \\ - {\bar{\bar{Z}}}_{0}^{JM} & {\bar{\bar{Z}}}_{0}^{MM} \end{matrix}] [\begin{matrix} {\bar{I}}_{j}^{J} \\ {\bar{I}}_{j}^{M} \end{matrix}] = [\begin{matrix} {\bar{V}}_{j}^{J} \\ {\bar{V}}_{j}^{M} \end{matrix}] - \sum_{j^{'} = 1}^{j - 1} [\begin{matrix} {\bar{\bar{Z}}}_{j - j^{'}}^{JJ} & {\bar{\bar{Z}}}_{j - j^{'}}^{JM} \\ - {\bar{\bar{Z}}}_{j - j^{'}}^{JM} & {\bar{\bar{Z}}}_{j - j^{'}}^{MM} \end{matrix}] [\begin{matrix} {\bar{I}}_{j^{'}}^{J} \\ {\bar{I}}_{j^{'}}^{M} \end{matrix}] .

Here, vectors

{\bar{V}}_{j}^{J}

and

{\bar{V}}_{j}^{M}

of size

N_{s} \times 1

store the tested incident fields and their elements are given by

{{\bar{V}}_{j}^{J}}_{n} = - {⟨ f_{n} (r), \partial_{t} E_{0}^{inc} (r, t) ⟩}_{t = j Δ t},

{{\bar{V}}_{j}^{M}}_{n} = - {⟨ f_{n} (r), \partial_{t} H_{0}^{inc} (r, t) ⟩}_{t = j Δ t} .

Note that here

{⟨ f_{n} (r), X (r, t) ⟩}_{t = j Δ t} = \int_{S_{n}} {f_{n} (r) \cdot X (r, t) d s |}_{t = j Δ t}

represents the inner product operation due to spatial testing. The elements of the matrices

{\bar{\bar{Z}}}_{j - j^{'}}^{JJ}

,

{\bar{\bar{Z}}}_{j - j^{'}}^{JM}

, and

{\bar{\bar{Z}}}_{j - j^{'}}^{MM}

of size

N_{s} \times N_{s}

are given by

{{\bar{\bar{Z}}}_{j - j^{'}}^{JJ}}_{n, n^{'}} = {{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 0}}_{n, n^{'}} + {{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 1}}_{n, n^{'}},

{{\bar{\bar{Z}}}_{j - j^{'}}^{JM}}_{n, n^{'}} = {{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 0}}_{n, n^{'}} + {{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 1}}_{n, n^{'}},

{{\bar{\bar{Z}}}_{j - j^{'}}^{MM}}_{n, n^{'}} = {{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 0}}_{n, n^{'}} + {{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 1}}_{n, n^{'}},

where

{{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 0}}_{n, n^{'}} = - μ_{0} {⟨ f_{n} (r), L_{0} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} + \frac{1}{ε_{0}} {⟨ f_{n} (r), Q_{0} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t},

{{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 0}}_{n, n^{'}} = - {⟨ f_{n} (r), K_{0} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t},

{{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 0}}_{n, n^{'}} = \frac{ε_{0}}{μ_{0}} {{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 0}}_{n, n^{'}},

{{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 1}}_{n, n^{'}} = - μ_{1} {⟨ f_{n} (r), L_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} + \frac{1}{ε_{1}^{\infty}} {⟨ f_{n} (r), Q_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} + {⟨ f_{n} (r), Q_{1} {{\bar{γ}}_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t},

{{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 1}}_{n, n^{'}} = - {⟨ f_{n} (r), K_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t},

{{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 1}}_{n, n^{'}} = - ε_{1}^{\infty} {⟨ f_{n} (r), L_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} - {⟨ f_{n} (r), L_{1} {γ_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} + \frac{1}{μ_{1}} {⟨ f_{n} (r), Q_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} .

Here,

{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, p}

,

{\bar{\bar{Z}}}_{j - j^{'}}^{JM, p}

, and

{\bar{\bar{Z}}}_{j - j^{'}}^{MM, p}

,

p \in {0, 1}

, correspond to contributions from the scattered fields generated in the background medium (

p = 0

) and inside the nanostructure (

p = 1

). Their entries, which are given by Eqs. (30)–(35), are computed using the schemes described in Section 2.E. Once

{\bar{\bar{Z}}}_{j - j^{'}}^{JJ}

,

{\bar{\bar{Z}}}_{j - j^{'}}^{JM}

, and

{\bar{\bar{Z}}}_{j - j^{'}}^{MM}

are constructed, unknown vectors

{[{\bar{I}}_{j}^{J} {\bar{I}}_{j}^{M}]}^{T}

,

j = 1, \dots, N_{t}

, are obtained recursively by time marching, as described next [20–30]. First,

{[{\bar{I}}_{1}^{J} {\bar{I}}_{1}^{M}]}^{T}

at time

t = Δ t

is found by solving Eq. (24) with right-hand side

{[{\bar{V}}_{1}^{J} {\bar{V}}_{1}^{M}]}^{T}

(

j = 1

).

{[{\bar{I}}_{1}^{J} {\bar{I}}_{1}^{M}]}^{T}

are then used to compute the scattered fields at time

t = 2 Δ t

, which are added to the tested incident fields

{[{\bar{V}}_{2}^{J} {\bar{V}}_{2}^{M}]}^{T}

to yield the right-hand side of Eq. (24).

{[{\bar{I}}_{2}^{J} {\bar{I}}_{2}^{M}]}^{T}

is found by solving Eq. (24) with this right-hand side (

j = 2

). Then,

{[{\bar{I}}_{1}^{J} {\bar{I}}_{1}^{M}]}^{T}

and

{[{\bar{I}}_{2}^{J} {\bar{I}}_{2}^{M}]}^{T}

are used to compute the scattered fields at time

t = 3 Δ t

, which together with

{[{\bar{V}}_{3}^{J} {\bar{V}}_{3}^{M}]}^{T}

form the right-hand side and permit the computation of

{[{\bar{I}}_{3}^{J} {\bar{I}}_{3}^{M}]}^{T}

, and so on. The computational cost of this time-marching scheme is dominated by that of computing the scattered fields, i.e., the discrete summation on the right-hand side of Eq. (24) (see Section 2.F).

E. Computation of the Matrix Entries

1. Contributions from the Background Medium ( $p = 0$ )

As described in Section 2.C, the Green function in the background medium, i.e., for $p = 0$ , $g_{0} (R, t) = δ (t - R / c_{0}) / (4 π R)$ . Inserting this in the expressions of ${{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 0}}_{n, n^{'}}$ , ${{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 0}}_{n, n^{'}}$ , and ${{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 0}}_{n, n^{'}}$ in Eqs. (30)–(32) and using the definitions of integral operators $L_{0} {\cdot}$ , $Q_{0} {\cdot}$ , and $K_{0} {\cdot}$ in Eqs. (7)–(9), one can obtain

{{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 0}}_{n, n^{'}} = \frac{ε_{0}}{μ_{0}} {{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 0}}_{n, n^{'}} = - \frac{ε_{0}}{4 π} \int_{S_{n}} f_{n} (r) \cdot \int_{S_{n^{'}}} \frac{f_{n^{'}} (r^{'}) \partial_{t}^{2} T_{j^{'}} (t) |_{t = j Δ t - R / c_{0}}}{R} d s^{'} d s - \frac{1}{4 π μ_{0}} \int_{S_{n}} \nabla \cdot f_{n} (r) \int_{S_{n^{'}}} \frac{\nabla^{'} \cdot f_{n^{'}} (r^{'}) T_{j^{'}} (j Δ t - R / c_{0})}{R} d s^{'} d s,

{{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 0}}_{n, n^{'}} = \frac{1}{4 π} \int_{S_{n}} f_{n} (r) \cdot \int_{S_{n^{'}}} \frac{\hat{R} \times f_{n^{'}} (r^{'})}{R} {[\frac{\partial_{t^{'}}^{2} T_{j^{'}} (t^{'})}{c_{0}} + \frac{\partial_{t} T_{j^{'}} (t^{'})}{R}]}_{t^{'} = j Δ t - R / c_{0}} d s^{'} d s .

In Eq. (37),

\hat{R} = (r - r^{'}) / R

. The surface integrals present in Eqs. (36) and (37) are computed using numerical quadrature [46] with proper singularity treatment [47]. One can also use semi-analytical expressions provided in [25,48–50].

2. Contributions from the Plasmonic Medium ( $p = 1$ )

These require computation of two types of convolutions.

Type-1 Convolutions: Type-1 refers to convolutions of $g_{1} (R, t)$ with $\nabla \cdot f_{n^{'}} (r) T_{j^{'}} (t)$ and $f_{n^{'}} (r) \partial_{t}^{2} T_{j^{'}} (t)$ and the convolution of $\partial_{R} g_{1} (R, t)$ with $f_{n^{'}} (r) \partial_{t} T_{j^{'}} (t)$ . Let $C_{T} (R, t)$ , $C (R, t)$ , and $C_{R} (R, t)$ denote the following temporal convolutions:

C_{T} (R, t) = g_{1} (R, t) * \partial_{t}^{2} T_{j^{'}} (t),

C (R, t) = g_{1} (R, t) * T_{j^{'}} (t),

C_{R} (R, t) = \partial_{R} g_{1} (R, t) * \partial_{t} T_{j^{'}} (t) .

Using the definitions in Eqs. (38)–(40), one can rewrite the expressions of

{⟨ f_{n} (r), L_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

,

{⟨ f_{n} (r), Q_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

, and

{⟨ f_{n} (r), K_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

, which are needed to compute

{{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 1}}_{n, n^{'}}

,

{{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 1}}_{n, n^{'}}

, and

{{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 1}}_{n, n^{'}}

:

{⟨ f_{n} (r), L_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} = \int_{S_{n}} f_{n} (r) \cdot \int_{S_{n^{'}}} C_{T} (R, t) f_{n^{'}} (r^{'}) {d s^{'} d s |}_{t = j Δ t},

{⟨ f_{n} (r), Q_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} = - \int_{S_{n}} \nabla \cdot f_{n} (r) \int_{S_{n^{'}}} C (R, t) \nabla^{'} \cdot f_{n^{'}} (r^{'}) {d s^{'} d s |}_{t = j Δ t},

{⟨ f_{n} (r), K_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} = \int_{S_{n}} f_{n} (r) \cdot \int_{S_{n^{'}}} \hat{R} \times f_{n^{'}} (r^{'}) C_{R} (R, t) {d s^{'} d s |}_{t = j Δ t} .

The surface integrals in Eqs. (41)–(43) are computed using numerical quadrature [46]. For every quadrature point pair, $C_{T} (R, t)$ , $C (R, t)$ , and $C_{R} (R, t)$ have to be computed. Their explicit expressions are obtained by inserting Eqs. (18) and (21) into Eqs. (41)–(43):

C_{T} (R, t) = [d^{g} (R) + \frac{1}{4 π R}] {\partial_{t^{'}}^{2} T (t^{'}) |}_{t^{'} = t - j^{'} Δ t - R / c_{1}} + \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{g} (R) \int_{- Δ t}^{\min (t - j^{'} Δ t - R / c_{1}, d Δ t)} u (t - t^{'} - j^{'} Δ t - R / c_{1}) e^{- a_{m}^{g} (R) (t - t^{'} - j^{'} Δ t - R / c_{1})} \partial_{t^{'}}^{2} T (t^{'}) d t^{'},

C (R, t) = [d^{g} (R) + \frac{1}{4 π R}] T (t - j^{'} Δ t - R / c_{1}) + \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{g} (R) \int_{- Δ t}^{\min (t - j^{'} Δ t - R / c_{1}, d Δ t)} u (t - t^{'} - j^{'} Δ t - R / c_{1}) e^{- a_{m}^{g} (R) (t - t^{'} - j^{'} Δ t - R / c_{1})} T (t^{'}) d t^{'},

C_{R} (R, t) = [f^{g_{R}} (R) - \frac{1}{4 π R c_{1}}] \partial_{t^{'}}^{2} {T (t^{'}) |}_{t^{'} = t - j^{'} Δ t - R / c_{1}} + [d^{g_{R}} (R) - \frac{1}{4 π R^{2}}] \partial_{t^{'}} {T (t^{'}) |}_{t^{'} = t - j^{'} Δ t - R / c_{1}} + \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{g_{R}} (R) \int_{- Δ t}^{\min (t - j^{'} Δ t - R / c_{1}, d Δ t)} u (t - t^{'} - j^{'} Δ t - R / c_{1}) e^{- a_{m}^{g_{R}} (R) (t - t^{'} - j^{'} Δ t - R / c_{1})} \partial_{t^{'}} T (t^{'}) d t^{'} .

For every quadrature point pair,

C_{T} (R, t)

,

C (R, t)

, and

C_{R} (R, t)

are computed using Eqs. (44)–(46). For these computations, terms involving

1 / R

and

1 / R^{2}

are treated carefully using the proper singularity extraction schemes [47]. The boundaries of the time integrals are determined using the support of

T (t)

,

\partial_{t} T (t)

, and

\partial_{t}^{2} T (t)

and the unit step function

u (t - t^{'} - j^{'} Δ t - R / c_{1})

. It should also be noted here that these integrals are computed using closed-form expressions since

T (t)

,

\partial_{t} T (t)

, and

\partial_{t}^{2} T (t)

consist of polynomial functions [38].

Type-2 Convolutions: Type-2 refers to the convolutions of $g_{1} (R, t)$ with $\nabla \cdot f_{n^{'}} (r) {\bar{γ}}_{1} (t) * T_{j^{'}} (t)$ and $f_{n^{'}} (r) γ_{1} (t) * \partial_{t}^{2} T_{j^{'}} (t)$ . Note the presence of the extra temporal convolution (in comparison with type-1 convolutions). Let $D (R, t)$ and $D_{T} (R, t)$ denote the following temporal convolutions:

D (R, t) = g_{1} (R, t) * {\bar{γ}}_{1} (t) * T_{j^{'}} (t),

D_{T} (R, t) = g_{1} (R, t) * γ_{1} (t) * \partial_{t}^{2} T_{j^{'}} (t) = \partial_{t}^{2} g_{1} (R, t) * γ_{1} (t) * T_{j^{'}} (t) .

Using the definitions in Eqs. (47) and (48), one can rewrite the expressions of

{⟨ f_{n} (r), Q_{1} {{\bar{γ}}_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

and

{⟨ f_{n} (r), L_{1} {γ_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

, which are needed to compute

{{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 1}}_{n, n^{'}}

and

{{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 1}}_{n, n^{'}}

:

{⟨ f_{n} (r), Q_{1} {{\bar{γ}}_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} = - \int_{S_{n}} \nabla \cdot f_{n} (r) \int_{S_{n^{'}}} D (R, t) \nabla^{'} \cdot f_{n^{'}} (r^{'}) d s^{'} d s |_{t = j Δ t},

{⟨ f_{n} (r), L_{1} {γ_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} = \int_{S_{n}} f_{n} (r) \cdot \int_{S_{n^{'}}} D_{T} (R, t) f_{n^{'}} (r^{'}) {d s^{'} d s |}_{t = j Δ t} .

The surface integrals in Eqs. (49) and (50) are computed using numerical quadrature [46]. But this computation can be done very efficiently by making use of type-1 convolutions that are already computed using the same quadrature points. In other words, one can find a way to represent the inner products

{⟨ f_{n} (r), Q_{1} {{\bar{γ}}_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

and

{⟨ f_{n} (r), L_{1} {γ_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

(type-2 convolutions) in terms of the inner products

{⟨ f_{n} (r), Q_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

and

{⟨ f_{n} (r), L_{1} {f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

(type-1 convolutions). This is described next. Note that convolutions

F_{j^{'}} (t) = {\bar{γ}}_{1} (t) * T_{j^{'}} (t)

and

F_{j^{'}}^{T} (t) = γ_{1} (t) * T_{j^{'}} (t)

do not depend on space but only on the material properties. Let

F_{j^{'}} (t)

and

F_{j^{'}}^{T} (t)

be discretized using the temporal basis function

T_{l^{'}} (t)

as

F_{j^{'}} (t) = {\bar{γ}}_{1} (t) * T_{j^{'}} (t) \approx \sum_{l^{'} = 1}^{N_{t}} F_{l^{'} - j^{'}} T_{l^{'}} (t),

F_{j^{'}}^{T} (t) = γ_{1} (t) * T_{j^{'}} (t) \approx \sum_{l^{'} = 1}^{N_{t}} F_{l^{'} - j^{'}}^{T} T_{l^{'}} (t),

where

F_{l^{'} - j^{'}}

and

F_{l^{'} - j^{'}}^{T}

are the expansion coefficients. To find them, Eqs. (51) and (52) are tested at times

s Δ t

, where

s

is an integer. Using the fact that

T ([s - l^{'}] Δ t) = 1

only when

s = l^{'}

and

T ([s - l^{'}] Δ t) = 0

otherwise, and inserting the expressions for

γ (t)

and

\bar{γ} (t)

from Eqs. (14) and (15) into the resulting equations, one can obtain

F_{l^{'} - j^{'}} = F_{j^{'}} (l^{'} Δ t) = {\bar{γ}}_{1} (t) * {T_{j^{'}} (t) |}_{t = l^{'} Δ t} = \int_{- Δ t}^{\min (l^{'} - j^{'}, d) Δ t} {\bar{γ}}_{1} ([l^{'} - j^{'}] Δ t - t^{'}) T (t^{'}) d t^{'} = \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{\bar{ϵ}} \int_{- Δ t}^{\min (l^{'} - j^{'}, d) Δ t} u ([l^{'} - j^{'}] Δ t - t^{'}) e^{- a_{m}^{\bar{ϵ}} ([l^{'} - j^{'}] Δ t - t^{'})} T (t^{'}) d t^{'},

F_{l^{'} - j^{'}}^{T} = F_{j^{'}}^{T} (l^{'} Δ t) = γ_{1} (t) * T_{j^{'}} (t) |_{t = l^{'} Δ t} = \int_{- Δ t}^{\min (l^{'} - j^{'}, d) Δ t} γ_{1} ([l^{'} - j^{'}] Δ t - t^{'}) T (t^{'}) d t^{'} = \sum_{m = 1}^{N_{p}^{VF}} b_{m}^{ϵ} \int_{- Δ t}^{\min (l^{'} - j^{'}, d) Δ t} u ([l^{'} - j^{'}] Δ t - t^{'}) e^{- a_{m}^{ϵ} ([l^{'} - j^{'}] Δ t - t^{'})} T (t^{'}) d t^{'} .

Note that the boundaries of the time integrals in Eqs. (53) and (54) are determined using the supports of

T (t)

and the unit step function

u ([l^{'} - j^{'}] Δ t - t^{'})

. These integrals are computed using closed-form expressions since

T (t)

consists of polynomial functions, respectively [38]. It should also be noted here

F_{l^{'} - j^{'}} = 0

and

F_{l^{'} - j^{'}}^{T} = 0

for

j^{'} \geq l^{'} + 1

. Inserting Eqs. (51) and (52) into Eqs. (47) and (48), respectively, yields

D (R, t) = g_{1} (R, t) * F_{l^{'} - j^{'}} (t) \approx \sum_{l^{'} = 1}^{N_{t}} F_{l^{'} - j^{'}} g_{1} (R, t) * T_{l^{'}} (t),

D_{T} (R, t) = \partial_{t}^{2} g_{1} (R, t) * F_{l^{'} - j^{'}}^{T} (t) \approx \sum_{l^{'} = 1}^{N_{t}} F_{l^{'} - j^{'}}^{T} \partial_{t}^{2} g_{1} (R, t) * T_{l^{'}} (t) = \sum_{l^{'} = 1}^{N_{t}} F_{l^{'} - j^{'}}^{T} g_{1} (R, t) * \partial_{t}^{2} T_{l^{'}} (t) .

Finally, inserting Eqs. (55) and (56) into Eqs. (49) and (50), respectively, yields

{⟨ f_{n} (r), Q_{1} {{\bar{γ}}_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} = - \int_{S_{n}} \nabla \cdot f_{n} (r) \int_{S_{n^{'}}} [\sum_{l^{'} = 1}^{N_{t}} F_{l^{'} - j^{'}} g_{1} (R, t) * T_{l^{'}} (t)] \nabla^{'} \cdot f_{n^{'}} (r^{'}) {d s^{'} d s |}_{t = j Δ t} = \sum_{l^{'} = 1}^{N_{t}} {⟨ f_{n} (r), Q_{1} {f_{n^{'}} (r) T_{l^{'}} (t)} ⟩}_{t = j Δ t} F_{l^{'} - j^{'}},

{⟨ f_{n} (r), L_{1} {γ_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t} = \int_{S_{n}} f_{n} (r) \cdot \int_{S_{n^{'}}} [\sum_{l^{'} = 1}^{N_{t}} F_{l^{'} - j^{'}}^{T} g_{1} (R, t) * \partial_{t}^{2} T_{l^{'}} (t)] f_{n^{'}} (r^{'}) {d s^{'} d s |}_{t = j Δ t} = \sum_{l^{'} = 1}^{N_{t}} {⟨ f_{n} (r), L_{1} {f_{n^{'}} (r) T_{l^{'}} (t)} ⟩}_{t = j Δ t} F_{l^{'} - j^{'}}^{T} .

Equations (57) and (58) show that

{⟨ f_{n} (r), Q_{1} {{\bar{γ}}_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

and

{⟨ f_{n} (r), L_{1} {γ_{1} (t) * f_{n^{'}} (r) T_{j^{'}} (t)} ⟩}_{t = j Δ t}

(type-2 convolutions) can be easily computed using

{⟨ f_{n} (r), Q_{1} {f_{n^{'}} (r) T_{l^{'}} (t)} ⟩}_{t = j Δ t}

and

{⟨ f_{n} (r), L_{1} {f_{n^{'}} (r) T_{l^{'}} (t)} ⟩}_{t = j Δ t}

(type-1 convolutions).

F. Computational Complexity

In this section, computational complexity of the MOT-PMCHWT-SIE solver is described step by step following the derivation of the matrix entries (Section 2.E) and the MOT scheme (Section 2.D):

(i) Computation of ${{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 0}}_{n, n^{'}}$ , ${{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 0}}_{n, n^{'}}$ , and ${{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 0}}_{n, n^{'}}$ in Eqs. (36) and (37) scales with $O (N_{s}^{2})$ because $g_{0} (R, t)$ is simply a Dirac function in time and space (no temporal tail).
(ii) Computation of type-1 convolutions in Eqs. (41)–(43) has a complexity of $O (N_{t} N_{s}^{2})$ because $g_{1} (R, t)$ has a temporal tail.
(iii) Convolutions in Eqs. (53) and (54) depend only on the permittivity of the plasmonic medium. Therefore $F_{l^{'} - j^{'}}$ and $F_{l^{'} - j^{'}}^{T}$ are computed once and stored in the memory. Cost of computing $F_{l^{'} - j^{'}}$ and $F_{l^{'} - j^{'}}^{T}$ scales with $O (N_{t})$ because $γ_{1} (t)$ and ${\bar{γ}}_{1} (t)$ have temporal tails.
(iv) Computation of type-2 convolutions in Eqs. (57) and (58) scales with $O (N_{t} N_{s}^{2})$ as a consequence of (ii) and (iii). As a result, computation of ${{\bar{\bar{Z}}}_{j - j^{'}}^{JJ, 1}}_{n, n^{'}}$ , ${{\bar{\bar{Z}}}_{j - j^{'}}^{JM, 1}}_{n, n^{'}}$ , and ${{\bar{\bar{Z}}}_{j - j^{'}}^{MM, 1}}_{n, n^{'}}$ scales with $O (N_{t} N_{s}^{2})$ . This assumes that the results of discrete summations in Eqs. (57) and (58) are directly incorporated into the matrix elements.

Finally, the computational complexity of the MOT scheme given in Eq. (24) scales with $O (N_{t}^{2} N_{s}^{2})$ . This cost can be reduced to $O (N_{t} N_{s} \log^{2} (N_{t}) \log (N_{s}))$ using the PWTD method [26] or $O (N_{t} N_{s}^{3 / 2} \log (N_{s}) \log (N_{t} N_{s}) \log (N_{t}))$ using FFT-based schemes [27–30,34]. In such cases, the discrete summations in Eqs. (57) and (58) should not be precomputed into the MOT matrices but should be incorporated into time marching through the use of $\sum_{j^{'} = 1}^{l^{'} - 1} F_{l^{'} - j^{'}} {\bar{I}}_{j^{'}}^{J}$ and $\sum_{j^{'} = 1}^{l^{'} - 1} F_{l^{'} - j^{'}}^{T} {\bar{I}}_{j^{'}}^{M}$ . This computation can be done very efficiently using a blocked FFT scheme without impacting the computational complexity of the accelerated time marching [29]. The use of these accelerated MOT schemes reduce the computational complexity of the proposed solver to essentially that of FDTD schemes.

3. NUMERICAL RESULTS

In this section, accuracy and applicability of the proposed MOT-PMCHWT-SIE solver are demonstrated through its application to the analysis of scattering from several nanostructures. In all examples, it is assumed that the nanostructure is residing in free space and excited by a plane wave with electric field

E_{0}^{inc} (r, t) = \hat{p} E_{0}^{inc} G (t - \hat{k} \cdot r / c_{0}),

where

E_{0}^{inc} = 1 V / m

is the amplitude,

\hat{p}

is the polarization unit vector,

\hat{k} = \hat{x} \sin θ^{inc} \cos φ^{inc} + \hat{y} \sin θ^{inc} \sin φ^{inc} + \hat{z} \cos θ^{inc}

is the unit vector along the direction of propagation defined by the angles

θ^{inc}

and

φ^{inc}

,

G (t) = \cos (2 π f_{0} [t - t_{0}]) e^{- {(t - t_{0})}^{2} / 2 σ^{2}}

is a Gaussian pulse with modulation frequency

f_{0}

, duration

σ

, and delay

t_{0}

. In all examples,

t_{0} = 8 σ

and

σ = 3 / (2 π f_{bw})

, where

f_{bw}

denotes an effective band. This specific selection of

σ

ensures that 99.998% of the incident energy is within the frequency band

[f_{0} - f_{bw}, f_{0} + f_{bw}]

. The free-space wavelengths at the minimum and maximum frequency of the excitation are

λ_{\max} = c_{0} / [f_{0} - f_{bw}]

and

λ_{\min} = c_{0} / [f_{0} + f_{bw}]

.

In all examples, the unknown coefficients of the equivalent electric and magnetic current densities, $J_{n} (k Δ t)$ and $M_{n} (k Δ t)$ , $k = 1, \dots, N_{t}$ , $n = 1, \dots, N_{s}$ , are computed by the proposed MOT-PMCHWT-SIE solver under this excitation. After the time-domain simulation is completed, the Fourier transforms of $J_{n} (k Δ t)$ and $M_{n} (k Δ t)$ are computed using the discrete time Fourier transform (DTFT). Normalizing the results with the DTFT of the samples of the Gaussian pulse, $G (k Δ t)$ , yields coefficients of the frequency-domain (i.e., time harmonic) electric and magnetic current densities denoted by $J_{n} (l Δ ω)$ and $M_{n} (l Δ ω)$ , $l = 1, \dots, N_{f}$ , where $Δ ω = 2 π Δ f$ , $Δ f$ is the frequency step, and $N_{f}$ is the number of frequency samples. This is followed by the computation of the frequency-domain scattering and extinction cross sections, $C_{sca} (l Δ ω)$ and $C_{ext} (l Δ ω)$ [51]:

C_{sca} (l Δ ω) = \frac{1}{16 π^{2} {| E_{0}^{inc} |}^{2}} \int_{Ω} {| F (\hat{r}, l Δ ω) |}^{2} d Ω,

C_{ext} (l Δ ω) = \frac{1}{k_{0} {| E_{0}^{inc} |}^{2}} Im {E_{0}^{inc} \hat{p} \cdot F (\hat{k}, l Δ ω)} .

Here,

\hat{r} = \hat{x} \sin θ \cos φ + \hat{y} \sin θ \sin φ + \hat{z} \cos θ

is the unit vector along the direction defined by angles

θ

and

φ

,

d Ω = \sin θ d θ d φ

is the differential solid angle, and

F (\hat{r}, l Δ ω)

represents the scattered electric field pattern in the far field and computed by inserting

J_{n} (l Δ ω)

and

M_{n} (l Δ ω)

into

F (\hat{r}, l Δ ω) = - j (l Δ ω) μ_{0} N (\hat{r}, l Δ ω) + j k_{0} \hat{r} \times L (\hat{r}, l Δ ω),

where

N (\hat{r}, l Δ ω) = \sum_{n = 1}^{N_{s}} J_{n} (l Δ ω) \int_{S_{n}} f_{n} (r^{'}) e^{j k_{0} r^{'} \cdot \hat{r}} d s^{'},

L (\hat{r}, l Δ ω) = \sum_{n = 1}^{N_{s}} M_{n} (l Δ ω) \int_{S_{n}} f_{n} (r^{'}) e^{j k_{0} r^{'} \cdot \hat{r}} d s^{'} .

Note that here

k_{0} = (l Δ ω) \sqrt{ε_{0} μ_{0}}

, it is assumed that

r \to \infty

during the derivation of Eqs. (62)–(64), and

\hat{r} \cdot F (\hat{r}, l Δ ω) = 0

is explicitly enforced during the computation of

C_{sca} (l Δ ω)

. Extinction efficiency of a scatterer,

Q_{ext} (l Δ ω)

, is defined as the ratio of

C_{ext} (l Δ ω)

to its geometrical cross section on the plane perpendicular to the direction of propagation of the incident field

\hat{k}

. In all examples presented in this section,

Q_{ext} (l Δ ω)

,

C_{ext} (l Δ ω)

, and

C_{sca} (l Δ ω)

computed by the proposed MOT-PMCHWT-SIE solver are compared to those obtained by the frequency-domain (FD) PMCHWT-SIE solver of [12], which computes

J_{n} (l Δ ω)

and

M_{n} (l Δ ω)

on the same discretization but directly in the frequency domain. Note that FD-PMCHWT-SIE solver has to be executed

N_{f}

times for this operation.

A. Accuracy of FRVF Scheme

In this section, the accuracy of the expansions in Eqs. (12) and (13) is demonstrated. The samples of the scatterer’s permittivity, $ε_{1} (ω_{i})$ , $i = 1, \dots, N_{ω}^{VF}$ , which are used in the FRVF scheme to generate the terms in these expansions, are obtained from the experimental data of Johnson–Christy tabulated for gold and silver in [43]. It should be noted here that experiments were carried out at only 49 frequency points; therefore, spline interpolation is used to generate samples at $ω_{i} \in 2 π [f_{\min}^{VF}, f_{\max}^{VF}]$ , $i = 1, \dots, N_{ω}^{VF}$ , $N_{ω}^{VF} = 1000$ required by Eqs. (12) and (13). Here, $f_{\min}^{VF} = 155$ , and $f_{\max}^{VF} = 1595 THz$ (corresponding to free-space wavelengths of $λ_{\max}^{VF} = c_{0} / f_{\min}^{VF} = 1937$ and $λ_{\min}^{VF} = c_{0} / f_{\max}^{VF} = 188 nm$ ), and the number of terms in the expansions $N_{p}^{VF} = 100$ . Figures 1 and 2 compare $ε_{1} (ω)$ in Eq. (12) obtained using the FRVF to that provided by the four-term Lorentz model described in [7], for gold and silver, respectively. The figures clearly show that both the FRVF expansion and the Lorentz model generate the real part of $ε_{1} (ω)$ accurately. However, the imaginary part generated by the Lorentz model does not match the experimental data. For the sake of completeness, Figs. 3(a) and 3(b) plot amplitude of $γ_{1} (t)$ and ${\bar{γ}}_{1} (t)$ , which correspond to the Fourier transform of the summation in Eqs. (12) and (13), i.e., $ε_{1} (t)$ and ${\bar{ε}}_{1} (t)$ without the Dirac delta term [see Eqs. (14) and (15)].

Fig. 1. (a) Real and (b) imaginary parts of $ε_{1} (ω)$ for gold.

Download Full Size | PDF

Fig. 2. (a) Real and (b) imaginary parts of $ε_{1} (ω)$ for silver.

Download Full Size | PDF

Fig. 3. Amplitude of (a) $γ_{1} (t)$ and (b) ${\bar{γ}}_{1} (t)$ for gold and silver.

Download Full Size | PDF

B. Gold and Silver Spheres

In this example, the scatterer is a silver or gold sphere of radius 50 nm. The permittivities of gold and silver are obtained using the Johnson–Christy [43] experimental data. For the FRVF scheme, $N_{p}^{VF} = 100$ , $N_{ω}^{VF} = 1000$ , $λ_{\min}^{VF} = 188$ , and $λ_{\max}^{VF} = 1937 nm$ . The excitation parameters are $\hat{p} = \hat{x}$ , $\hat{k} = \hat{z}$ , $f_{0} = 900$ , $f_{bw} = 600 THz$ , $λ_{\min} = 200$ , and $λ_{\max} = 1000 nm$ . The currents induced on the sphere surface are discretized using $N_{s} = 1926$ RWG basis functions and the simulation is executed for $N_{t} = 2500$ time steps with step size $Δ t = 0.0333 fs$ . Figures 4(a) and 4(b) plot the amplitudes of $J_{1} (k Δ t)$ and $M_{1} (k Δ t)$ , $k = 1, \dots, N_{t}$ computed during the simulations of the gold and silver spheres, respectively. Figures clearly show the stability of the results. Figures 5(a) and 5(b) compare $Q_{ext} (l Δ ω)$ , $l = 1, \dots, N_{f}$ , $N_{f} = 100$ , $Δ f = 12 THz$ , computed by the MOT-PMCHWT-SIE and FD-PMCHWT-SIE solvers, and those obtained from Mie series solution. Results agree well with each other.

Fig. 4. Amplitude of $J_{1} (t)$ and $M_{1} (t)$ induced on (a) the gold sphere and (b) the silver sphere.

Download Full Size | PDF

Fig. 5. $Q_{ext} (ω)$ computed for (a) the gold sphere and (b) the silver sphere.

Download Full Size | PDF

C. Gold Rounded Cube

Next, scattering from a gold rounded cube [52] is analyzed using the proposed MOT-PMCHWT-SIE solver. The dimension of the cube is 200 nm and the radius of fillets on the edges and corners is 20 nm. The permittivity of gold is obtained using the Johnson–Christy [43] experimental data. For the FRVF scheme, $N_{p}^{VF} = 20$ , $N_{ω}^{VF} = 1000$ , $λ_{\min}^{VF} = 292$ , and $λ_{\max}^{VF} = 1088 nm$ . The excitation parameters are $\hat{p} = \hat{x}$ , $\hat{k} = \hat{z}$ , $f_{0} = 685$ , $f_{bw} = 315 THz$ , $λ_{\min} = 300$ , and $λ_{\max} = 810 nm$ . The currents induced on the cube surface are discretized using $N_{s} = 6543$ RWG basis functions and the simulation is executed for $N_{t} = 900$ time steps with step size $Δ t = 0.05 fs$ . Figure 6 compares $Q_{ext} (l Δ ω)$ , $l = 1, \dots, N_{f}$ , $N_{f} = 200$ , $Δ f = 3.15 THz$ , obtained using the proposed MOT-PMCHWT-SIE solver, to that computed by the null-field method [52]. Results are in good agreement with each other.

Fig. 6. $C_{ext} (ω)$ computed for the gold rounded cube.

Download Full Size | PDF

D. Gold Rounded Triangular Prism

In this example, the scatterer is a gold rounded triangular prism of height 40 nm and edge length 200 nm [52]. The radius of the fillets on the edges and the corners is 10 nm. The permittivity of gold is obtained using the Johnson–Christy [43] experimental data. For the FRVF scheme, $N_{p}^{VF} = 20$ , $N_{ω}^{VF} = 1000$ , $λ_{\min}^{VF} = 292$ , and $λ_{\max}^{VF} = 1088 nm$ . The excitation parameters are $\hat{p} = \hat{x}$ , $\hat{k} = \hat{z}$ , $f_{0} = 685$ , $f_{bw} = 315 THz$ , $λ_{\min} = 300$ , and $λ_{\max} = 810 nm$ . The currents induced on the cube surface are discretized using $N_{s} = 6543$ RWG basis functions and the simulation is executed for $N_{t} = 900$ time steps with step size $Δ t = 0.05 fs$ . Figure 7 compares $C_{ext} (l Δ ω)$ , $l = 1, \dots, N_{f}$ , $N_{f} = 200$ , $Δ f = 3.15 THz$ obtained using the MOT-PMCHWT-SIE and FD-PMCHWT-SIE solvers to that computed by the DDA [52]. Results obtained by the MOT-PMCHWT-SIE and FD-PMCHWT-SIE solvers match well with each other, but the result computed by the DDA differs from those, especially around the peak values. The same type of mismatch is also observed in [52] and explained by the fact that the DDA loses accuracy as the permittivity of the scatterer increases [53].

Fig. 7. $C_{ext} (ω)$ computed for the gold rounded triangular prism.

Download Full Size | PDF

E. Gold-Coated Silica Sphere

Scattering from a silica sphere coated with gold [7] is analyzed using the proposed MOT-PMCHWT-SIE solver. The radius of the silica sphere and the thickness of the gold layer are 40 nm and 10 nm, respectively. The relative permittivity of the silica is 2.04 while the permittivity of the gold is obtained using the Johnson–Christy [43] experimental data. For the FRVF scheme, $N_{p}^{VF} = 100$ , $N_{ω}^{VF} = 1000$ , $λ_{\min}^{VF} = 188$ , and $λ_{\max}^{VF} = 1937 nm$ . The excitation parameters are $\hat{p} = \hat{x}$ , $\hat{k} = \hat{z}$ , $f_{0} = 940$ , $f_{bw} = 560 THz$ , $λ_{\min} = 200$ , and $λ_{\max} = 788 nm$ . The currents induced on the sphere surface are discretized using $N_{s} = 2130$ RWG basis functions and the simulation is executed for $N_{t} = 2000$ time steps with step size $Δ t = 0.0333 fs$ . Figure 8 plots $C_{ext} (l Δ ω)$ , $l = 1, \dots, N_{f}$ , $N_{f} = 100$ , $Δ f = 11.2 THz$ computed by the MOT-PMCHWT-SIE solver to that obtained from the Mie series solution. Results are in good agreement.

Fig. 8. $Q_{ext} (ω)$ computed for the gold shell with silica core.

Download Full Size | PDF

F. Silver Dimer

In this example, the scatterer is an $x$ -directed dimer consisting of two silver spheres residing on the $x y$ -plane [54]. The radius of the spheres and the shortest distance between them are 15 nm and 1.5 nm, respectively. The permittivity of silver is obtained using the Palik experimental data [44]. For the FRVF scheme, $N_{p}^{VF} = 100$ , $N_{ω}^{VF} = 1000$ , $λ_{\min}^{VF} = 196$ , and $λ_{\max}^{VF} = 1033 nm$ . Two simulations are carried out for two polarizations of the incident field, $\hat{p} = \hat{x}$ and $\hat{p} = \hat{y}$ . The remaining excitation parameters are $\hat{k} = \hat{z}$ , $f_{0} = 750$ , $f_{bw} = 250 THz$ , $λ_{\min} = 300$ , and $λ_{\max} = 600 nm$ . The currents induced on the sphere surfaces are discretized using $N_{s} = 1428$ RWG basis functions. Both simulations are executed for $N_{t} = 2000$ time steps with step size $Δ t = 0.0333 fs$ . Figure 9 compares $C_{ext} (l Δ ω)$ , $l = 1, \dots, N_{f}$ , $N_{f} = 100$ , $Δ f = 5 THz$ , computed by the MOT-PMCHWT-SIE and FD-PMCHWT-SIE solvers. Results agree well with each other and also with those provided in [54].

Fig. 9. $C_{ext} (ω)$ computed for the silver dimer under two excitations with different polarizations.

Download Full Size | PDF

G. Gold Disk with a Nonconcentric Cavity

In the last example, the scatterer is a gold disk embedding a nonconcentric cavity and residing on the $x y$ -plane [55]. The radius and thickness of the disk are 50 nm and 10 nm, respectively. The radius of the circular cavity is 20 nm and its center is offset from the center of the disk by 29 nm in the $x$ -direction. The permittivity of gold is obtained using the Johnson–Christy [43] experimental data. For the FRVF scheme, $N_{p}^{VF} = 20$ , $N_{ω}^{VF} = 1000$ , $λ_{\min}^{VF} = 496$ , and $λ_{\max}^{VF} = 1216 nm$ . The excitation parameters are $\hat{p} = \hat{y}$ , $θ^{inc} = 78$ and $φ^{inc} = 0$ , $f_{0} = 450$ , $f_{bw} = 150 THz$ , $λ_{\min} = 500$ , and $λ_{\max} = 1000 nm$ . The currents induced on the surface of the scatterer are discretized using $N_{s} = 4380$ RWG basis functions. The simulation is executed for $N_{t} = 1500$ time steps with step size $Δ t = 0.083 fs$ . Figure 10 compares $C_{sca} (l Δ ω)$ , $l = 1, \dots, N_{f}$ , $N_{f} = 100$ , $Δ f = 3 THz$ , obtained using the proposed MOT-PMCHWT-SIE solver, to that computed by the frequency-domain FEM [55]. The results agree with each other reasonably well. The positions of the two main peaks in the results are close enough yet the peak at 660 nm is less pronounced and slightly blueshifted in the result obtained by the MOT-PMCHWT-SIE solver.

Fig. 10. $C_{sca} (ω)$ computed for the gold disk with the nonconcentric cavity.

Download Full Size | PDF

4. CONCLUSION

An MOT scheme for solving the PMCHWT-SIE enforced on surfaces of plasmonic nanostructures is described. The unknown equivalent electric and magnetic current densities introduced on these surfaces are expanded by RWG and polynomial basis functions in space and time, respectively. Inserting this expansion into the PMCHWT-SIE and Galerkin testing the resulting equation at discrete times yield a system of equations. This system is then solved for the unknown expansion coefficients using the MOT scheme.

The PMCHWT-SIE requires evaluation of additional convolutions involving the plasmonic medium’s permittivity and Green function and the temporal basis function to compute the scattered fields. This convolution is discretized in a way that is fully consistent with the MOT scheme. The computation of the discretized convolution is carried out with almost no additional cost and without changing the computational complexity of the MOT scheme. Time-domain samples of the permittivity and the Green function required by this computation are obtained from their frequency-domain samples using the FRVF algorithm. Frequency samples are generated using tabulated data obtained from experiments.

Numerical results that demonstrate the accuracy of the MOT-PMCHWT-SIE solver are presented. The quantum-corrected version of the same solver, which can accurately be applied to scatterers separated by subnanometer distances, is under development.

APPENDIX A

This section briefly describes the fundamental idea behind the FRVF scheme. The details can be found in [39–41]. Let $F (ω)$ represent a function defined in the frequency domain. FRVF approximates $F (ω)$ as a sum of rational functions:

F (ω) \approx j ω f + d + \sum_{m = 1}^{N} \frac{b_{m}}{j ω + a_{m}} .

Here,

N

is termed fitting order,

d

and

f

are optional parameters, and

a_{m}

and

b_{m}

are poles and residues, respectively, while

d

and

f

are real or can be forced to be zero, and

a_{m}

and

b_{m}

are either real or have to come in complex conjugate pairs.

The FRVF scheme is iteratively executed. Each iteration has two stages: identification of poles and residues of the right-hand side of Eq. (A1). First, poles are identified as described next. Let $S (ω)$ and $S (ω) F (ω)$ represent two auxiliary functions:

S (ω) = 1 + \sum_{m = 1}^{N} \frac{b_{m}^{S}}{j ω + {\tilde{a}}_{m}},

S (ω) F (ω) \approx j ω f^{SF} + d^{SF} + \sum_{m = 1}^{N} \frac{b_{m}^{SF}}{j ω + {\tilde{a}}_{m}} .

Here,

{\tilde{a}}_{m}

are “guessed” poles of

S (ω)

and

S (ω) F (ω)

for the given iteration. Multiplying Eq. (A2) with

F (ω)

and equating the resulting expression to the right-hand side of Eq. (A3), and evaluating the final equation at frequency samples

ω_{i}

,

i = 1, \dots, N_{samp}

, where

N < N_{samp}

, yield an overdetermined system:

\bar{\bar{A}} \bar{X} = \bar{Y} .

Here,

{[\bar{\bar{A}}]}_{i} = [\begin{matrix} \frac{1}{j ω_{i} + {\tilde{a}}_{1}} \dots \frac{1}{j ω_{i} + {\tilde{a}}_{N}} 1 & j ω_{i} \frac{- F (ω_{i})}{j ω_{i} + {\tilde{a}}_{1}} \dots \frac{- F (ω_{i})}{j ω_{i} + {\tilde{a}}_{N}} \end{matrix}]

represents the

i

th row of

\bar{\bar{A}}

,

\bar{X} = {[\begin{matrix} b_{1}^{SF} \dots b_{N}^{SF} & d^{SF} & f^{SF} & b_{1}^{S} \dots b_{N}^{S} \end{matrix}]}^{T}

is the unknown vector, and

{\bar{Y}}_{i} = F (ω_{i})

. Unknown

\bar{X}

is obtained by solving Eq. (A4) in the least-squares sense. One can show that zeros of the

S (ω)

are equal to the poles of

F (ω)

by rewriting Eqs. (A1)–(A3) in terms of sum of partial fractions [39]. These zeros are equal to the eigenvalues of the matrix

\bar{\bar{D}} - \bar{H} {\bar{B}}^{T}

, where

\bar{\bar{D}}

is a diagonal matrix with entries

{\bar{\bar{D}}}_{m, m} = {\tilde{a}}_{m}

,

{\bar{H}}_{m} = 1

, and

{\bar{B}}_{m} = b_{m}^{S}

. This completes the identification of the poles of the right-hand side of Eq. (A1) for the given iteration. Once these poles are known, the residues are identified as described next. Eq. (A1) is evaluated at

ω_{i}

,

i = 1, \dots, N_{samp}

, yielding an overdetermined system as in Eq. (A4), where now

{[\bar{\bar{A}}]}_{i} = [\begin{matrix} \frac{1}{j ω_{i} + a_{1}} \dots \frac{1}{j ω_{i} + a_{N}} & 1 & j ω_{i} \end{matrix}]

and

\bar{X} = {[\begin{matrix} b_{1} \dots b_{N} & d & f \end{matrix}]}^{T}

. Solving this system in the least-squares sense yields the residues of the right-hand side of Eq. (A1) for the given iteration and completes the iteration. Poles obtained in one iteration are used as an initial guess in the next one and the iterations are repeated until a preset level of accuracy is obtained, i.e., an accurate representation for

F (ω)

is obtained.

REFERENCES

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

2. J. N. Anker, W. P. Hall, O. Lyandres, N. C. Shah, J. Zhao, and R. P. V. Duyne, “Biosensing with plasmonic nanosensors,” Nat. Mater. 7, 442–453 (2008). [CrossRef]

3. D. J. de Aberasturi, A. B. Serrano-Montes, and L. M. Liz-Marzán, “Modern applications of plasmonic nanoparticles: from energy to health,” Adv. Opt. Mater. 3, 602–617 (2015). [CrossRef]

4. S. I. Bozhevolnyi, V. S. Volkov, E. Devaux, J.-Y. Laluet, and T. W. Ebbesen, “Channel plasmon subwavelength waveguide components including interferometers and ring resonators,” Nature 440, 508–511 (2006). [CrossRef]

5. C. Oubre and P. Nordlander, “Optical properties of metallodielectric nanostructures calculated using the finite difference time domain method,” J. Phys. Chem. B 108, 17740–17747 (2004). [CrossRef]

6. A. Vial, A.-S. Grimault, D. Macas, D. Barchiesi, and M. L. de La Chapelle, “Improved analytical fit of gold dispersion: application to the modeling of extinction spectra with a finite-difference time-domain method,” Phys. Rev. B 71, 085416 (2005). [CrossRef]

7. F. Hao and P. Nordlander, “Efficient dielectric function for FDTD simulation of the optical properties of silver and gold nanoparticles,” Chem. Phys. Lett. 446, 115–118 (2007). [CrossRef]

8. Y. Hu, S. J. Noelck, and R. A. Drezek, “Symmetry breaking in gold-silica-gold multilayer nanoshells,” ACS Nano 4, 1521–1528 (2010). [CrossRef]

9. K. L. Kelly, E. Coronado, L. L. Zhao, and G. C. Schatz, “The optical properties of metal nanoparticles: the influence of size, shape, and dielectric environment,” J. Phys. Chem. B 107, 668–677 (2003). [CrossRef]

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

11. M. Araújo, J. Taboada, J. Rivero, D. Sols, and F. Obelleiro, “Solution of large-scale plasmonic problems with the multilevel fast multipole algorithm,” Opt. Lett. 37, 416–418 (2012). [CrossRef]

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

13. B. Gallinet, A. M. Kern, and O. J. 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]

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

15. A. Taflove and S. C. Hagness, Computational Electrodynamics (Artech House, 2005).

16. K. Sirenko, V. Pazynin, Y. K. Sirenko, and H. Bagci, “An FFT-accelerated FDTD scheme with exact absorbing conditions for characterizing axially symmetric resonant structures,” Prog. Electromagn. Res. 111, 331–364 (2011). [CrossRef]

17. P. G. Petropoulos, “Stability and phase error analysis of FD-TD in dispersive dielectrics,” IEEE Trans. Antennas Propag. 42, 62–69 (1994). [CrossRef]

18. J. H. Greene and A. Taflove, “General vector auxiliary differential equation finite-difference time-domain method for nonlinear optics,” Opt. Express 14, 8305–8310 (2006). [CrossRef]

19. A. Karlsson and S. Rikte, “Time-domain theory of forerunners,” J. Opt. Soc. Am. A 15, 487–502 (1998). [CrossRef]

20. D. Vechinski and S. M. Rao, “Transient scattering from two-dimensional dielectric cylinders of arbitrary shape,” IEEE Trans. Antennas Propag. 40, 1054–1060 (1992). [CrossRef]

21. G. Manara, A. Monorchio, and R. Reggiannini, “A space-time discretization criterion for a stable time-marching solution of the electric field integral equation,” IEEE Trans. Antennas Propag. 45, 527–532 (1997). [CrossRef]

22. H. Bagci, A. Yilmaz, J.-M. Jin, and E. Michielssen, “Fast and rigorous analysis of EMC/EMI phenomena on electrically large and complex cable-loaded structures,” IEEE Trans. Electromagn. Compat. 49, 361–381 (2007). [CrossRef]

23. A. A. Ergin, B. Shanker, and E. Michielssen, “The plane-wave time-domain algorithm for the fast analysis of transient wave phenomena,” IEEE Trans. Antennas Propag. 41, 39–52 (1999). [CrossRef]

24. B. Shanker, A. A. Ergin, M. Lu, and E. Michielssen, “Fast analysis of transient electromagnetic scattering phenomena using the multilevel plane wave time domain algorithm,” IEEE Trans. Antennas Propag. 51, 628–641 (2003). [CrossRef]

25. B. Shanker, M. Lu, J. Yuan, and E. Michielssen, “Time domain integral equation analysis of scattering from composite bodies via exact evaluation of radiation fields,” IEEE Trans. Antennas Propag. 57, 1506–1520 (2009). [CrossRef]

26. P.-L. Jiang and E. Michielssen, “Multilevel plane wave time domain-enhanced MOT solver for analyzing electromagnetic scattering from objects residing in lossy media,” Proc. IEEE 3B, 447–450 (2005). [CrossRef]

27. A. E. Yilmaz, J.-M. Jin, and E. Michielssen, “Time domain adaptive integral method for surface integral equations,” IEEE Trans. Antennas Propag. 52, 2692–2708 (2004). [CrossRef]

28. H. Bagci, A. E. Yilmaz, V. Lomakin, and E. Michielssen, “Fast solution of mixed-potential time-domain integral equations for half-space environments,” IEEE Trans. Geosci. Remote Sens. 43, 269–279 (2005). [CrossRef]

29. H. Bagci, A. Yilmaz, and E. Michielssen, “An FFT-accelerated time-domain multiconductor transmission line simulator,” IEEE Trans. Electromagn. Compat. 52, 199–214 (2010). [CrossRef]

30. A. Yilmaz, D. S. Weile, B. Shanker, J.-M. Jin, and E. Michielssen, “Fast analysis of transient scattering in lossy media,” IEEE Antennas Wireless Propag. Lett. 1, 14–17 (2002). [CrossRef]

31. X. Wang and D. S. Weile, “Electromagnetic scattering from dispersive dielectric scatterers using the finite difference delay modeling method,” IEEE Trans. Antennas Propag. 58, 1720–1730 (2010). [CrossRef]

32. X. Wang and D. S. Weile, “Implicit Runge-Kutta methods for the discretization of time domain integral equations,” IEEE Trans. Antennas Propag. 59, 4651–4663 (2011). [CrossRef]

33. C. Lubich, “Convolution quadrature and discretized operational calculus. I,” Numer. Math. 52, 129–145 (1988). [CrossRef]

34. C. Lubich and A. Ostermann, “Linearly implicit time discretization of non-linear parabolic equations,” IMA J. Numer. Anal. 15, 555–583 (1995). [CrossRef]

35. J. Putnam and L. Medgyesi-Mitschang, “Combined field integral equation formulation for inhomogeneous two and three-dimensional bodies: the junction problem,” IEEE Trans. Antennas Propag. 39, 667–672 (1991). [CrossRef]

36. K. Donepudi, J.-M. Jin, and W. C. Chew, “A higher order multilevel fast multipole algorithm for scattering from mixed conducting/dielectric bodies,” IEEE Trans. Antennas Propag. 51, 2814–2821 (2003). [CrossRef]

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

38. I. E. Uysal, H. A. Ulku, and H. Bagci, “MOT solution of the PMCHWT equation for analyzing transient scattering from conductive dielectrics,” IEEE Antennas Wireless Propag. Lett. 14, 507–510 (2015). [CrossRef]

39. B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Trans. Power Del. 14, 1052–1061 (1999). [CrossRef]

40. B. Gustavsen, “Improving the pole relocating properties of vector fitting,” IEEE Trans. Power Del. 21, 1587–1592 (2006). [CrossRef]

41. D. Deschrijver, M. Mrozowski, T. Dhaene, and D. De Zutter, “Macromodeling of multiport systems using a fast implementation of the vector fitting method,” IEEE Microwave Wireless Compon. Lett. 18, 383–385 (2008). [CrossRef]

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

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

44. E. D. Palik, Handbook of Optical Constants of Solids (Academic, 1998), Vol. 3.

45. P. M. Morse and H. Feshbach, Methods of Theoretical Physics (McGraw-Hill, 1953).

46. S. Rao, Time Domain Electromagnetics (Elsevier, 1999).

47. D. Wilton, S. Rao, A. Glisson, D. Schaubert, O. Al-Bundak, and C. Butler, “Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains,” IEEE Trans. Antennas Propag. 32, 276–281 (1984). [CrossRef]

48. H. A. Ulku and A. A. Ergin, “Analytical evaluation of transient magnetic fields due to RWG current bases,” IEEE Trans. Antennas Propag. 55, 3565–3575 (2007). [CrossRef]

49. H. A. Ulku and A. A. Ergin, “On the singularity of the closed-form expression of the magnetic field in time domain,” IEEE Trans. Antennas Propag. 59, 691–694 (2011). [CrossRef]

50. H. A. Ulku and A. A. Ergin, “Application of analytical retarded-time potential expressions to the solution of time domain integral equations,” IEEE Trans. Antennas Propag. 59, 4123–4131 (2011). [CrossRef]

51. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

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

53. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11, 1491–1499 (1994). [CrossRef]

54. E. R. Encina and E. A. Coronado, “Plasmon coupling in silver nanosphere pairs,” J. Phys. Chem. C 114, 3918–3923 (2010). [CrossRef]

55. M. Amin and H. Bagci, “Scattering properties of vein induced localized surface plasmon resonances on a gold disk,” in Proceedings of High Capacity Optical Networks and Enabling Technologies (HONET), Riyadh, Saudi Arabia (2011), pp. 237–240.

Transient analysis of electromagnetic wave interactions on plasmonic nanostructures using a surface integral equation solver

Abstract

1. INTRODUCTION

2. FORMULATION

A. Time-Domain PMCHWT-SIE

B. Time-Domain Permittivity and Inverse Permittivity

C. Time-Domain Green Function

D. Discretization and MOT Scheme

E. Computation of the Matrix Entries

1. Contributions from the Background Medium ( $p = 0$ )

2. Contributions from the Plasmonic Medium ( $p = 1$ )

F. Computational Complexity

3. NUMERICAL RESULTS

A. Accuracy of FRVF Scheme

B. Gold and Silver Spheres

C. Gold Rounded Cube

D. Gold Rounded Triangular Prism

E. Gold-Coated Silica Sphere

F. Silver Dimer

G. Gold Disk with a Nonconcentric Cavity

4. CONCLUSION

APPENDIX A

REFERENCES

Cited By

Figures (10)

Equations (68)

Journal of the Optical Society of America A

Abstract

1. INTRODUCTION

2. FORMULATION

A. Time-Domain PMCHWT-SIE

B. Time-Domain Permittivity and Inverse Permittivity

C. Time-Domain Green Function

D. Discretization and MOT Scheme

E. Computation of the Matrix Entries

1. Contributions from the Background Medium (p=0)

2. Contributions from the Plasmonic Medium (p=1)

F. Computational Complexity

3. NUMERICAL RESULTS

A. Accuracy of FRVF Scheme

B. Gold and Silver Spheres

C. Gold Rounded Cube

D. Gold Rounded Triangular Prism

E. Gold-Coated Silica Sphere

F. Silver Dimer

G. Gold Disk with a Nonconcentric Cavity

4. CONCLUSION

APPENDIX A

REFERENCES

Cited By

Figures (10)

Equations (68)

Journal of the Optical Society of America A

1. Contributions from the Background Medium ( $p = 0$ )

2. Contributions from the Plasmonic Medium ( $p = 1$ )