## Abstract

A time domain surface integral equation (TD-SIE) solver is developed for quantum-corrected analysis of transient electromagnetic field interactions on plasmonic nanostructures with sub-nanometer gaps. “Quantum correction” introduces an auxiliary tunnel to support the current path that is generated by electrons tunneled between the nanostructures. The permittivity of the auxiliary tunnel and the nanostructures is obtained from density functional theory (DFT) computations. Electromagnetic field interactions on the combined structure (nanostructures plus auxiliary tunnel connecting them) are computed using a TD-SIE solver. Time domain samples of the permittivity and the Green function required by this solver are obtained from their frequency domain samples (generated from DFT computations) using a semi-analytical method. Accuracy and applicability of the resulting quantum-corrected solver scheme are demonstrated via numerical examples.

© 2017 Optical Society of America

## 1. Introduction

Recent advances in nano-fabrication have permitted prototyping and production of plasmonic structures with sub-nanometer holes and gaps [1, 2]. When two nanostructures are separated with a gap shorter than 0.5 nm, effects of quantum (electron) tunneling are expected to become visible on the structures’ scattering spectrum [3–8]. Indeed this has been demonstrated experimentally for a single dimer [9–13] and an ensemble of dimers [14, 15]. Consequently, this means that, simulation tools, which are capable of accurately modeling quantum-tunneling effects, are needed to design geometrically complicated nanoscale plasmonic devices. Unfortunately, numerical schemes developed for solving classical equations of electromagnetics are not intrinsically equipped to account for tunneling of electrons between two nanostructures. One could use a quantum mechanical solver for this purpose, but this approach would be computationally very expensive since the whole device has to be modeled at the atomic level. One way to overcome this bottleneck and account for electron tunneling accurately without sacrificing from efficiency is to incorporate the quantum mechanical solution locally (around and inside the gap) into a classical electromagnetic solver. This approach has been termed “quantum-corrected model” (QCM) in [5, 6]. The QCM proposed in [5–8] replaces the gap between two nanostructures with an auxiliary tunnel whose permittivity is represented using a Drude model [5]. To determine the Drude model’s parameters, first, quantum-tunneling probability is computed by solving the one-dimensional (1D) Schrödinger equation for a single electron between the nanostructures (i.e., on a domain with size equal to the length of the gap). This tunneling probability is used to obtain the tunnel’s conductivity. Finally, the conductivity is incorporated into the Drude model via a damping parameter [5, 6] or by adding it directly to the imaginary part of the permittivity [7, 8]. Once the permittivity of the auxiliary tunnel is known, a frequency domain surface integral equation (FD-SIE) solver [5–7] or a frequency domain finite element method (FD-FEM) [8] is used to characterize the time harmonic electromagnetic field interactions on the combined structure (two nanostructures plus auxiliary tunnel connecting them). It should be noted here that another Drude representation is used for modeling the (bulk) permittivity of the nanostructures [5].

The approach summarized above suffers from three drawbacks: *(i)* During the solution of the 1D Schrödinger equation for a single electron, atomic structure of the material is not accounted for. This might result in inaccuracies: It has been argued in [16] that modeling many-body interactions of electrons and taking the atomic structure into account yield more accurate values for optical parameters of a material. *(ii)* The Drude model used for representing the (bulk) permittivity of the nanostructures does not take into account the inter-band contributions. Including only intra-band contributions produces permittivity values that do not match to those obtained by experiments [17]. Indeed, experimentally generated permittivity models have been used in numerical characterization of gold and silver nanoparticles in many different studies [9, 10, 12]. *(iii)* A classical electromagnetic solver, which relies on a volumetric discretization, like the FD-FEM used in [8], requires a high sampling rate in the vicinity of the nanostructure surfaces to accurately capture the fast decaying plasmonic fields. This significantly increases the number of unknowns to be solved for.

To increase the accuracy of the permittivity model for the auxiliary tunnel, a first-principles approach of density functional theory (DFT) can be used. Indeed, in [18] the conductivity of the tunnel is extracted from DFT computations assuming that the tunnel is connecting two slabs of sodium represented by the jellium model. In [19], time dependent density functional theory (TDDFT) computations are used together jellium approximation (for sodium) to find the polarizability and the permittivity within the tunnel. In this work, the permittivity of both of the auxiliary tunnel and nanostructures is obtained from DFT computations [20, 21]. Unlike the approaches used in [18, 19], jellium approximation/model is avoided. The DFT computations used in this work account for many-body interactions of electrons in a crystal structure of atoms as well as inter- and intra-band contributions to the permittivity [21–23] [avoiding drawbacks *(i)* and *(ii)*]. Then, the permittivity is represented as a sum of rational functions in frequency domain. Weights of the rational functions are extracted using the fast relaxed vector fitting (FRVF) scheme [24–26] on the frequency domain permittivity samples, which are obtained from the DFT computations. To facilitate the simulation of transient electromagnetic field interactions, time domain permittivity is obtained analytically by applying inverse Laplace transform to the sum of rational functions. It should be noted here that the computation of the time domain permittivity using the procedure described here does not assume any specific model but requires only frequency domain samples of the permittivity to be known [avoiding drawback *(ii)*].

At the final stage, transient electromagnetic field interactions on the combined structure (two nanostructures plus the auxiliary tunnel connecting them) are computed using a time domain surface integral equation (TD-SIE) solver. More specifically, the Poggio-Miller-Chan-Harrington-Wu-Tsai surface integral equation (TD-PMCHWT-SIE) formulation, which can account for different dielectric volumes, each of which represents the nanostructures, the tunnel, and the background medium, is used [27–29]. The TD-PMCHWT-SIE is enforced on the interfaces between these dielectric volumes, equivalent electric and magnetic currents are introduced on these interfaces, and they are expanded in terms of spatial and temporal basis functions [30–32]. Inserting this expansion into the TD-PMCHWT-SIE and Galerkin testing the resulting equation at discrete times yield a system of equations. This system is solved for the unknown expansion coefficients of the currents by time marching. It should be noted here that this scheme calls for computation of the time domain samples of the unbounded-medium Green function associated with each dielectric volume. Those are obtained from frequency domain samples of the Green functions using the same procedure used for computing the time domain samples of the permittivity. It should be emphasized here that the TD-PMCHWT-SIE solver calls for only the discretization of the interfaces between different dielectric volumes (instead of a volumetric discretization), enforces the radiation condition implicitly without using approximate absorbing boundary conditions, and uses a time step size that does not depend on the spatial discretization, i.e., time step size is not limited to a Courant-Frederich-Lewy (CFL) like condition [33] [avoiding drawback *(iii)*].

## 2. Formulation

This section describes the steps followed in quantum-corrected analysis of transient electromagnetic field interactions. Figure 1 presents a flow diagram summarizing these steps, each of which is detailed in Sections 2.1–2.4

#### 2.1. Tunnel model

It has been theoretically and experimentally shown that the quantum tunneling is observed when the distance between two structures is smaller than 0.5 nm [5–13]. Additionally, the permittivity of the auxiliary tunnel introduced between the two nanostructures depends on this distance [5]. It should be noted here that, since this distance might vary depending on the shape of the nanostructure surfaces, the auxiliary tunnel is constructed layer by layer. In other words, the distance between the two surfaces is sampled and each layer corresponds to a “sampled” tunneling gap distance. The permittivity of a layer is then obtained from the DFT computations that use the sampled distance associated with that layer (see Section 2.2). As expected higher number of layers translates to higher accuracy in the tunnel model since the spatial dependence of the permittivity is sampled with a higher accuracy. This is especially true for shorter tunneling distances since the variation in the permittivity of the resulting auxiliary tunnel is much faster.

For an example of how the tunnel model is constructed, consider a dimer consisting of two spheres, as shown in Fig. 2(a). In this case, the tunneling gap distance, which is defined as the distance between two spherical surfaces along the direction that is perpendicular to the dimer axis, is clearly not constant and varies between *d*_{min} (shortest distance) and *d*_{max} = 0.55 nm (distance beyond which the tunneling effects are assumed to vanish). It should be noted here the schematics in Figs. 2(a) and 2(b) are not drawn to scale. Figure 2(b) shows three (cylindrical) layers used for representing the auxiliary tunnel in Fig 2(a). End surfaces (caps) of the layers are not flat but curved surfaces that match the surfaces of the spheres. This means that each of these layers must be assigned an effective length as described next. For the inner most layer [shown with red in Fig. 2(b)] the effective layer length is the shortest distance between the two “caps” of the layer. For the other two layers, effective length is the average between the shortest and longest length of the layer. Finally, effective length of a given layer is used in DFT computations to obtain the permittivity of that layer (see Section 2.2).

#### 2.2. Permittivity from DFT computations

Information about structural, electronic, magnetic, and chemical properties of a material can be obtained by studying the quantum mechanical wave function of its electrons. However, computation of the wave functions in a many-electron system from the approximation-free solution of the Schrödinger equation is not easy if not impossible. The DFT overcomes this problem by using density functional of electrons as the quantity of interest instead of the wave function. Indeed, all related information of a material can also be retrieved from the density functional of electrons [34, 35]. The DFT formulates the many-electron problem in a crystal with fixed nuclei (i.e., under Born-Oppenheimer approximation) using a self-consistent density functional. Even after this simplification, electron density functional cannot be solved for using analytical methods and a numerical scheme has to be used [36]. In this work, the DFT computations are carried out using the commercially-available WIEN2k solver [20], which makes use of the full potential linearized augmented plane wave (FP-LAPW) method [37]. Additionally, WIEN2k solver uses the generalized gradient approximation (Perdew-Burke-Ernzerhof functional) to account for the exchange correlation interaction [38]. The rest of this section describes how the permittivity of the bulk material and auxiliary tunnel is computed.

The complex dielectric permittivity tensor ${\overline{\overline{\epsilon}}}^{\text{DFT}}(\omega )={\overline{\overline{\epsilon}}}^{\text{inter}}(\omega )+{\overline{\overline{\epsilon}}}^{\text{intra}}(\omega )$ has two contributions from inter-band and intra-band transitions of electrons. The imaginary part of inter-band contribution to the dielectric permittivity tensor is given by [21, 22]

*n*≠

*n*′) states. In Eq. (1),

*ħ*is the reduced Planck constant,

*ϵ*

_{0}is the permittivity of free space,

*q*is the charge of the electron,

_{e}*m*is the mass of the electron,

_{e}*p*

_{z;n,n}_{′},

**is the momentum matrix element along**

_{k}*z*-direction,

*f*(

*E*

_{n}_{,}

**) is the Fermi distribution,**

_{k}*E*

_{n}_{,}

**is the electron energy,**

_{k}**k**is the crystal wave vector. The real part of the inter-band contribution of the dielectric permittivity tensor is computed using the Kramers-Kronig relation [39]

*P*denotes the principal value integral. The intra-band contribution is given by a Drude model

*γ*is the lifetime broadening (or damping frequency) and ${\omega}_{p,zz}^{2}$ is the plasma frequency computed as

*E*is the Fermi energy and

_{F}*n*runs over all states.

Two different configurations are considered. In the first configuration, for the bulk material, all quantities of interest in Eqs. (1)–(4) are obtained from the DFT computations on a unit cell of the crystal with periodic boundary conditions [21, 22]. It should be noted here that due to the symmetry of the crystal, the permittivity tensor
${\overline{\overline{\epsilon}}}^{\text{DFT}}(\omega )$ obtained under this configuration has only one independent component
${\overline{\overline{\epsilon}}}_{zz}^{\text{DFT}}(\omega )$. This component is equal to the permittivity of the bulk material, i.e.,
${\epsilon}_{\text{bulk}}(\omega )={\overline{\overline{\epsilon}}}_{zz}^{\text{DFT}}(\omega )$. In the second configuration, for the auxiliary tunnel, all quantities of interest in Eqs. (1)–(4) are obtained from the DFT computations on a 1 × 1 × 6 supercell as shown in Fig. 3. The gap of length *d*_{gap}, which represents tunneling gap distance, is introduced between two groups of three unit cells, i.e., the gap is introduced between two metal blocks of three unit cells. This ensures that the resulting permittivity is not affected by the periodic boundary conditions. Due to the symmetry of the problem the dielectric permittivity tensor has two independent components
${\overline{\overline{\epsilon}}}_{xx}^{{}^{\text{DFT}}}(\omega )={\overline{\overline{\epsilon}}}_{yy}^{{}^{\text{DFT}}}(\omega )$ and
${\overline{\overline{\epsilon}}}_{zz}^{\text{DFT}}(\omega )$. However, only
${\overline{\overline{\epsilon}}}_{zz}^{\text{DFT}}(\omega )$ is used here since the auxiliary tunnel dimension along the *z*-direction is much larger than dimensions in other directions and the incident electric field is assumed to be *z*-polarized.

The supercell in Fig. 3 can be treated as a layered medium where the effective permittivity is given by
${\overline{\overline{\epsilon}}}_{zz}^{\text{DFT}}(\omega )$, which is computed using Eqs. (1)–(4). Then the permittivity of the gap *ε*_{gap}(*ω*) is obtained using the expression [40]

*h*is the length of metal blocks, i.e., three lattice constants. It should be noted here

*d*

_{gap}is set to the effective length of the auxiliary tunnel layer as explained in Section 2.1, and ${\overline{\overline{\epsilon}}}_{zz}^{\text{DFT}}(\omega )$ and

*ε*

_{gap}(

*ω*) are computed for every different value of

*d*

_{gap}.

#### 2.3. Weighted rational function fit

The bulk material and auxiliary tunnel permittivities *ε*_{bulk}(*ω*) and *ε*_{gap}(*ω*) are formulated and computed in frequency domain as described in Section 2.2. Let *ε*(*ω*) represent any one of *ε*_{bulk}(*ω*) or *ε*_{gap}(*ω*). To be able to analyze the transient electromagnetic field interactions on the nanostructure, one has to carry out simulations in time domain. Consequently, the TD-SIE solver described in Section 2.4 requires time domain samples of *ε*(*ω*),
$\overline{\epsilon}(\omega )=1/\epsilon (\omega )$, Green function
$G(R,\omega )={e}^{i}{}^{\omega R\sqrt{\epsilon (\omega ){\mu}_{0}}}/(4\pi R)$, and its spatial derivative
${\partial}_{R}G(R,\omega )=(i\omega \sqrt{\epsilon (\omega )\mu}-1/R){e}^{i\omega R\sqrt{\epsilon (\omega ){\mu}_{0}}}/(4\pi R)$, to be computed. Here, *µ*_{0} is the permeability in free space. The transformation from frequency domain to time domain cannot be carried out analytically since the frequency domain samples of *ε*(*ω*) are obtained numerically. Therefore, the numerical scheme initially proposed in [27] is used for computing the time domain samples. This scheme is briefly described next.

Let *F*(*ω*) represent any one of *ε*(*ω*),
$\overline{\epsilon}(\omega )$, *G*(*R, ω*), or *∂ _{R}G*(

*R, ω*). It is assumed that in a given frequency band

*F*(

*ω*) can be approximated using rational functions as

*N*is the number of rational functions,

*d*and

*f*are constants (which can explicitly be enforced to be zero) and

*a*and

_{k}*b*are the poles and the residues associated with the rational functions. To find the unknown coefficients,

_{k}*d*,

*f*,

*a*, and

_{k}*b*, the FRVF scheme [24–26] is used. The FRVF scheme minimizes the difference between the samples of

_{k}*F*(

*ω*) and the right hand side of Eq. (6) (computed in the given frequency band) to find the unknown coefficients. It should be noted here that during this operation, it is enforced that

*d*and

*f*are real valued, and

*a*and

_{k}*b*are real valued or come in complex conjugate pairs. Additionally, the FRVF scheme chooses stable poles, i.e., the resulting

_{k}*a*satisfy Re {

_{k}*a*} < 0. Once the coefficients are obtained, one can find the time domain expression for

_{k}*F*(

*t*) by inverse Fourier transforming the right hand side of Eq. (6):

*δ*(.) is the Dirac delta function,

*δ*′(.) is its first derivative, and

*u*(.) is the unit step function. The expression in Eq. (7) allows

*ε*(

*t*), $\overline{\epsilon}(t)$,

*G*(

*R, t*), and

*∂*(

_{R}G*R, t*) to be computed in closed form in time domain.

#### 2.4. TD-SIE

Let
$V={\cup}_{p=1}^{P}{V}_{p}$ represent the total volume of the combined structure consisting of the nanostructures and the auxiliary tunnel. Each one of the nanostructures and the layers of the tunnel is assigned a volume denoted by *V _{p}*,

*p*= 1, …,

*P*. The combined structure resides in an unbounded non-dispersive background medium represented with

*V*

_{0}. The permittivity and inverse permittivity of

*V*are denoted by

_{p}*ε*(

_{p}*t*) and ${\overline{\epsilon}}_{p}(t)$, respectively. For

*p*= 0, permittivity is constant, i.e.,

*ε*

_{0}(

*t*) =

*ε*

_{0}and ${\overline{\epsilon}}_{0}(t)=1/{\epsilon}_{0}$. All volumes are non-magnetic with the constant permeability

*µ*=

_{p}*µ*

_{0},

*p*= 0, …,

*P*. The surface between the two volumes

*V*and

_{q}*V*is represented by

_{p}*S*,

_{l}*l*= 1, …,

*L*,

*p*= 0, …,

*P*,

*q*= 0, …,

*P*,

*p*≠

*q*, and ${\widehat{\mathbf{n}}}_{l}(\mathbf{r})$ is the unit normal vector on

*S*pointing towards

_{l}*V*. Let $\{{\mathbf{E}}_{p}^{\text{inc}}(\mathbf{r},t),{\mathbf{H}}_{p}^{\text{inc}}(\mathbf{r},t)\}$ and $\{{\mathbf{E}}_{p}^{\text{sca}}(\mathbf{r},t),{\mathbf{H}}_{p}^{\text{sca}}(\mathbf{r},t)\}$ represent the incident and the scattered electromagnetic fields in

_{p}*V*, respectively. It is assumed that ${\mathbf{E}}_{p}^{\text{inc}}(\mathbf{r},t)$ and ${\mathbf{H}}_{p}^{\text{inc}}(\mathbf{r},t)$ are vanishingly small ∀

_{p}**r**∈

*V*,

*t*< 0 and essentially band limited to

*f*

_{max}. The fields $\{{\mathbf{E}}_{p}^{\text{inc}}(\mathbf{r},t),{\mathbf{H}}_{p}^{\text{inc}}(\mathbf{r},t)\}$ and $\{{\mathbf{E}}_{p}^{\text{sca}}(\mathbf{r},t),{\mathbf{H}}_{p}^{\text{sca}}(\mathbf{r},t)\}$ satisfy two boundary conditions for

**r**∈

*S*:

_{l}**J**

*(*

_{l}**r**,

*t*) and

**M**

*(*

_{l}**r**,

*t*) introduced on

*S*[27–29]:

_{l}*l*′ runs over the indices of the surfaces that “touches”

*V*and integral operators are defined as

_{p}Here, *G _{p}* (

*R, t*) is the Green function of the unbounded medium that has the same permittivity and permeability as

*V*(

_{p}*ε*(

_{p}*t*) and

*µ*),

_{p}*R*= |

**r**–

**r**′| is the distance between points

**r**and

**r**′ and “∗” denotes temporal convolution. It should be noted here that for

*p*= 0, temporal convolutions in Eqs. (10)–(13) are simplified using the fact that

*ε*

_{0}(

*t*) =

*ε*

_{0}, ${\overline{\epsilon}}_{0}(t)=1/{\epsilon}_{0}$, and ${G}_{0}(R,t)=\delta (t-R\sqrt{{\epsilon}_{0}{\mu}_{0}})/(4\pi R)$. For

*p*= 1, …,

*P*,

*ε*(

_{p}*t*), ${\overline{\epsilon}}_{p}(t)$,

*G*(

_{p}*R, t*) are computed as described in Sections 2.2 and 2.3. Inserting Eqs. (10) and (11) into Eqs. (8) and (9) yields TD-PMCHWT-SIE [27, 28] in unknown

**J**

*(*

_{l}**r**,

*t*) and

**M**

*(*

_{l}**r**,

*t*). To numerically solve the TD-PMCHWT-SIE,

*S*are discretized into triangular patches and unknowns

_{l}**J**

*(*

_{l}**r**,

*t*) and

**M**

*(*

_{l}**r**,

*t*) are expanded using the Rao-Wilton-Glisson basis functions in space [30] and Lagrange interpolation functions in time [31, 32]. The unknown coefficients of this expansion are computed using the marching on-in-time (MOT) scheme as described in [27–29]. The computation of the spatiotemporal convolutions involving

*ε*(

_{p}*t*), ${\overline{\epsilon}}_{p}(t)$,

*G*(

_{p}*R, t*),

*∂*(

_{R}G_{p}*R, t*), spatial basis functions, and temporal interpolation functions, which is required by this MOT scheme, is carried out efficiently using the method described in [27].

After the time marching is completed, i.e., all expansion coefficients are known, **J*** _{l}* (

**r**,

*t*) and

**M**

*(*

_{l}**r**,

*t*) can be constructed for

**r**∈

*S*,

_{l}*l*= 1, …,

*L*and $\{{\mathbf{E}}_{p}^{\text{sca}}(\mathbf{r},t),{\mathbf{H}}_{p}^{\text{sca}}(\mathbf{r},t)\}$ can be computed for

**r**∈

*V*,

_{p}*p*= 1, …,

*P*using Eqs. (10) and (11). It is worth mentioning here that these time-domain results (obtained with only one execution of the MOT scheme) can easily be converted to frequency-domain by Fourier transform to generate broadband data.

## 3. Numerical results

This section presents numerical examples that demonstrate the accuracy and applicability of the proposed approach. For all examples considered here, DFT computations are carried out by the WIEN2k package [20] that implements the FP-LAPW method [37]. The bulk material and auxiliary tunnel permittivities *ε*_{bulk}(*ω*) and *ε*_{gap}(*ω*) are then obtained from the DFT computations using the method described in Section 2.2. Two materials are considered: gold (Au) and silver (Ag). The parameters of the DFT simulations required by the WIEN2k implementation are selected as described next. Au and Ag have face-centered cubic lattices (belonging to space group Fm-3m, no. 225) with experimental lattice constants of 0.407 nm and 0.408 nm, respectively [41]. First, the crystal structures are fully optimized and the equilibrium lattice constants of 0.408 nm and 0.409 nm are obtained for Au and Ag, respectively. These optimized lattice constants are used for the DFT computations. All computations use 24 atoms and a dense 96 × 96 × 16 k-mesh. For the total energy convergence, the product *R*_{mt}*K*_{max} = 7. Here, *R*_{mt} is the smallest muffin-tin (MT) sphere radius and *K*_{max} is the magnitude of the largest wave vector **k** + **K _{n}** used in the plane wave expansion,

**k**is the wave vector inside the first Brillouin zone, and

**K**are the reciprocal lattice vectors. The core states are separated from the valance states by setting the cutoff energy to -6 Ry. The maximum angular momentum number for the wave function expansion inside the atomic spheres

_{n}*l*

_{max}= 10, while the magnitude of the largest wave vector

**G**used in the Fourier expansion of the charge density

*G*

_{max}= 12(Ry)

^{1/2}. The MT radii for Au and Ag are 2.50 Bohr. Self-consistent calculations are stopped when iterations have an energy difference less than 10

^{−4}Ry.

In all of the examples considered here, it is assumed that the nanostructures reside in free space with permittivity *ϵ*_{0} = 8.8541878176 × 10^{−12} F/m and permeability *μ*_{0} = 4*π* × 10^{−7} H/m and are excited by a plane wave with electric field

*f*

_{0}, duration

*σ*, and delay

*t*

_{0}. In all examples,

*t*

_{0}= 8

*σ*and

*σ*= 3/(2

*π f*

_{bw}), where

*f*

_{bw}denotes an effective bandwidth. With this selection of parameters, less than 99.998% of

*G*(

*t*)’s power resides within the frequency band [

*f*

_{0}−

*f*

_{bw},

*f*

_{0}+

*f*

_{bw}].

After the time domain simulation is completed, frequency domain (i.e., time-harmonic) equivalent electric and magnetic current densities are obtained by dividing the Fourier transform of **J*** _{l}* (

**r**,

*t*) and

**M**

*(*

_{l}**r**,

*t*),

*l*= 1, …,

*L*, by the Fourier transform of

*G*(

*t*). Then, frequency domain current densities are used to compute the extinction cross section

*C*

_{ext}(

*ω*) in frequency domain as explained in [27]. In all examples,

*C*

_{ext}(

*ω*) is plotted against the photon energy

*E*=

*ħω*.

#### 3.1. Permittivity values

The permittivity of the auxiliary tunnel *ε*_{gap}(*ω*) is computed for gold and silver using Eq. (5) with *d*_{gap} ∈ {0.1, 0.2, 0.3, 0.4, 0.5} nm. Figures 4 and 5 plot Re{*ε*_{gap}(*ω*) and Im{*ε*_{gap}(*ω*)} versus the photon energy *E* = *ħω* for gold and silver, respectively. Similarly, the bulk permittivity *ε*_{bulk}(*ω*) is obtained for gold and silver using
${\epsilon}_{\text{bulk}}(\omega )={\overline{\overline{\epsilon}}}_{zz}^{\text{DFT}}(\omega )$, where
${\overline{\overline{\epsilon}}}_{zz}^{\text{DFT}}(\omega )$ is computed as described in Section 2.2. Figures 6 and 7 plot Re{*ε*_{bulk}(*ω*)} and Im{*ε*_{bulk}(*ω*)} versus the photon energy *E* = *ħω* for gold and silver, respectively. Same figures also provide values of *ε*_{bulk}(*ω*) obtained using the Johnson-Christy [42] and the Drude [5] models. It is clearly shown that values of *ε*_{bulk}(*ω*) obtained using the DFT computations are closer to the experimental values obtained using the Johnson-Christy model especially at higher frequencies.

It is worth mentioning here that comparisons of Figs. 4(a) and 4(b) to Figs. 6(a) and 6(b) and Figs. 5(a) and 5(b) to Figs. 7(a) and 7(b) show that *ε*_{gap}(*ω*) (obtained from DFT computations) “converges” to *ε*_{bulk}(*ω*) (obtained from DFT computations) as *d*_{gap} approaches zero. This indicates that the models making use of DFT computations to obtain *ε*_{bulk}(*ω*) and *ε*_{gap}(*ω*) are consistent.

To demonstrate the effect of *ε*_{bulk}(*ω*) on the scattering properties, Mie series solution [43] is used to compute the extinction cross section *C*_{ext}(*ω*) of gold sphere of radius 25 nm and a silver sphere of radius 10 nm with *ε*_{bulk}(*ω*) obtained from the DFT computations and using the Johnson-Christy and Drude models. Figures 8(a) and 8(b) plot *C*_{ext}(*ω*) versus the photon energy *E* = *ħω* for gold and silver spheres, respectively. It is clear that *C*_{ext}(*ω*) computed for the spheres with *ε*_{bulk}(*ω*) obtained from the DFT computations is more accurate than that computed for the spheres with *ε*_{bulk}(*ω*) values obtained using the Drude model. These figures suggest that using *ε*_{bulk}(*ω*) and *ε*_{gap}(*ω*) both of which are obtained from the DFT computations might produce more accurate scattering properties than using those obtained from the Drude model. It should also be noted here one could use the Johnson-Christy model for *ε*_{bulk}(*ω*) and the Drude model for *ε*_{gap}(*ω*) as done in [9]. However, this approach is not consistent since *ε*_{gap}(*ω*) does not converge to *ε*_{bulk}(*ω*) as the length of the tunnel *d*_{gap} approaches zero [5, 6].

#### 3.2. Scattering from gold and silver spheres

In this example, the scatterers are a single gold sphere of radius 25 nm and a single silver sphere of radius 10 nm. The excitation parameters *f*_{0} = 785 THz and *f*_{bw} = 665 THz. The TD-SIE solver uses 3018 spatial basis functions to discretize the currents induced on the sphere surface and the simulation is executed for 1500 time steps with step size 0.03 fs. Figures. 9(a) and 9(b) compare *C*_{ext}(*ω*) obtained using the currents computed by the TD-SIE solver to *C*_{ext}(*ω*) analytically obtained from the Mie series solution for gold and silver spheres, respectively. Very good match is observed between the analytical and the numerical solutions demonstrating the accuracy of the TD-SIE solver.

#### 3.3. Scattering from a gold dimer

In this example, scattering from a dimer consisting of two gold spheres is analyzed. Spheres have a radius of 25 nm and are aligned along *z*-axis. The excitation parameters *f*_{0} = 785 THz and *f*_{bw} = 665 THz. Two groups of simulations are carried out. In both groups, five different simulations are executed while the distance between the spheres (*d*_{min}) is changed from 0.1 nm to 0.5 nm with 0.1 nm steps. In the first group, effects of quantum tunneling are ignored, i.e., no auxiliary tunnel is located between the spheres. The current densities induced on the spheres are discretized using 5352 spatial basis functions. In the second group of simulations, the auxiliary tunnel is constructed using four layers as explained in Section 2.1. The current densities induced on the surfaces of the spheres plus the auxiliary tunnel is discretized using 5616 spatial basis functions. Both groups of simulations are executed for 1500 time steps with step size 0.03 fs.

Figures 10(a) and 10(b) plot the amplitude of electric and magnetic current densities computed at **r** = (2.681, 0.794, 0.197) nm for the first and second groups of simulations with *d*_{min} = 0.1 nm. Figure 10(a) clearly shows that the current densities obtained by the simulation in the first group (i.e. without the auxiliary tunnel) oscillate with a single frequency. This indicates the presence of only one resonance. On the other hand, the oscillations of the current densities obtained by the simulation in the second group (i.e. with the auxiliary tunnel) appear more complicated due to the presence of an additional resonance. These plasmon resonances are identified more clearly in the frequency response of the dimer as demonstrated by Figures 11(a) and 11(b).

These two figures plot *C*_{ext}(*ω*) of the gold dimer obtained by the first and second groups of simulations. It should be noted here that, in both figures, curves representing *C*_{ext}(*ω*) for different values of *d*_{min} are shifted vertically for illustration purposes. In Fig. 11(a), it is clearly shown that the main plasmon resonance around 1.7 eV redshifts as *d*_{min} decreases [5, 6]. Figure 11(b) shows that, unlike the simulations without the quantum correction, the main plasmon resonance is observed around 1.9 eV and blueshifts as *d*_{min} decreases. Additionally, a second resonance emerges around 1.2 eV. This resonance mode is induced due to the electrons tunneled between two spheres and has been termed “charge transfer plasmon mode” in previous studies. It has been theoretically predicted in [3–5] and demonstrated by electron energy loss spectroscopy (EELS) experiments in [12, 13].

It is also worth mentioning here that *C*_{ext}(*ω*) obtained by the first and second groups of simulations converge to each other for the same value of *d*_{min} as it increases. This is expected since quantum tunneling effects become less and less dominant and the auxiliary tunnel starts acting like the background medium with permittivity *ϵ*_{0} as *d*_{min} increases.

Finally, to demonstrate the convergence of accuracy in the layered construction of the auxiliary tunnel, the simulation, which uses an auxiliary tunnel with *d*_{min} = 0.1 nm, is repeated as the number of layers used in the tunnel model is changed from one to three, four, and seven. Figure 12 plots *C*_{ext}(*ω*) of the dimer computed in each of these simulations. It clearly shows that four- and seven-layer models produces almost identical results demonstrating the convergence of tunnel model’s accuracy.

#### 3.4. Scattering from a silver dimer

In the last example, scattering from a dimer consisting of two silver spheres is analyzed. Spheres have a radius of 10 nm and are aligned along *z*-axis. The excitation parameters *f*_{0} = 506 THz and *f*_{bw} = 493 THz. Two groups of simulations are carried out. In both groups, five different simulations are executed while as *d*_{min} is changed from 0.1 nm to 0.6 nm with 0.1 nm steps. In the first group, effects of quantum tunneling are ignored, i.e., no auxiliary tunnel is located between the spheres. The current densities induced on the spheres are discretized using 4584 spatial basis functions. In the second group of simulations, the auxiliary tunnel is constructed using four layers as explained in Section 2.1. The current densities induced on the surfaces of the spheres plus the auxiliary tunnel is discretized using 4908 spatial basis functions. Both groups of simulations are executed for 1500 time steps with step size 0.03 fs.

Figures 13(a) and 13(b) plot *C*_{ext}(*ω*) of the silver dimer obtained by the first and second groups of simulations, respectively. Figure 13(a) clearly shows that the main plasmon resonance observed around 2.55 eV redshifts to 2.1 eV as *d*_{min} is reduced to 0.1 nm. Additionally, a second plasmon resonance, which is observed around 2.85 eV, redshifts to 2.62 eV, and a third resonance starts to appear around 2.81 eV as *d*_{min} is decreased. Figure 13(b) shows that the main plasmon resonance, which is observed around 2.66 eV for *d*_{min} = 0.5 nm, blueshifts to 2.76 eV for *d*_{min} = 0.1 nm. Figure 13(b) also shows that a second plasmon resonance appears around 1.25 eV and blueshifts to 1.48 eV when *d*_{min} = 0.1 nm. Like the results obtained for the gold dimer, this resonance is associated with the charge transfer plasmon mode as theoretically predicted in [3–5] and experimentally demonstrated in [12, 13].

Also like the results obtained for the gold dimer, *C*_{ext}(*ω*) of the silver dimer obtained by the first and second groups of simulations converge to each other for the same value of *d*_{min} as it increases. However, convergence seems to be slower than for the silver dimer and effects of quantum tunneling is still observed for *d*_{min} = 0.6 nm even though they are significantly less dominant than those for *d*_{min} = 0.1 nm.

## 4. Conclusion

A quantum-corrected TD-SIE solver is developed for analyzing electromagnetic field interactions on plasmonic nanostructures with sub-nanometer gaps. The sub-nanometer gap, where the quantum tunneling of electrons happen, is replaced with an auxiliary tunnel. The complex permittivities of the nanostructures and the auxiliary tunnel are obtained from the DFT computations. The electromagnetic field interactions on the combined structure (the nanostructures plus the auxiliary tunnel connecting them) are computed using a TD-SIE solver. The proposed approach is used for computing electromagnetic fields scattered from gold and silver dimers. The results clearly show the effect of quantum tunneling on the scattering spectrum of the dimers.

## References and links

**1. **M. S. Tame, K. McEnery, Ş. Özdemir, J. Lee, S. Maier, and M. Kim, “Quantum plasmonics,” Nat. Phys. **9**, 329–340 (2013). [CrossRef]

**2. **K. Kneipp, H. Kneipp, and J. Kneipp, “Probing plasmonic nanostructures by photons and electrons,” Chem. Sci. **6**, 2721–2726 (2015). [CrossRef]

**3. **I. Romero, J. Aizpurua, G. W. Bryant, and F. J. G. de Abajo, “Plasmons in nearly touching metallic nanoparticles: singular response in the limit of touching dimers,” Opt. Express **14**, 9988–9999 (2006). [CrossRef] [PubMed]

**4. **J. Zuloaga, E. Prodan, and P. Nordlander, “Quantum description of the plasmon resonances of a nanoparticle dimer,” Nano Lett. **9**, 887–891 (2009). [CrossRef] [PubMed]

**5. **R. Esteban, A. G. Borisov, P. Nordlander, and J. Aizpurua, “Bridging quantum and classical plasmonics with a quantum-corrected model,” Nat. Commun. **3**, 825 (2012). [CrossRef] [PubMed]

**6. **R. Esteban, A. Zugarramurdi, P. Zhang, P. Nordlander, F. J. Garcia-Vidal, A. G. Borisov, and J. Aizpurua, “A classical treatment of optical tunneling in plasmonic gaps: extending the quantum corrected model to practical situations,” Faraday Discuss. **178**, 151–183 (2015). [CrossRef] [PubMed]

**7. **U. Hohenester, “Quantum corrected model for plasmonic nanoparticles: A boundary element method implementation,” Phys. Rev. B **91**, 205436 (2015). [CrossRef]

**8. **L. Wu, H. Duan, P. Bai, M. Bosman, J. K. Yang, and E. Li, “Fowler–Nordheim tunneling induced charge transfer plasmons between nearly touching nanoparticles,” ACS Nano **7**, 707–716 (2013). [CrossRef]

**9. **K. J. Savage, M. M. Hawkeye, R. Esteban, A. G. Borisov, J. Aizpurua, and J. J. Baumberg, “Revealing the quantum regime in tunnelling plasmonics,” Nature **491**, 574–577 (2012). [CrossRef] [PubMed]

**10. **J. A. Scholl, A. García-Etxarri, A. L. Koh, and J. A. Dionne, “Observation of quantum tunneling between two plasmonic nanoparticles,” Nano Lett. **13**, 564–569 (2013). [CrossRef]

**11. **W. Zhu and K. B. Crozier, “Quantum mechanical limit to plasmonic enhancement as observed by surface-enhanced Raman scattering,” Nat. Commun. **5**, 5228 (2014). [CrossRef]

**12. **S. F. Tan, L. Wu, J. K. Yang, P. Bai, M. Bosman, and C. A. Nijhuis, “Quantum plasmon resonances controlled by molecular tunnel junctions,” Science **343**, 1496–1499 (2014). [CrossRef] [PubMed]

**13. **S. Kadkhodazadeh, J. B. Wagner, H. Kneipp, and K. Kneipp, “Coexistence of classical and quantum plasmonics in large plasmonic structures with subnanometer gaps,” Appl. Phys. Lett. **103**, 083103 (2013). [CrossRef]

**14. **H. Jung, H. Cha, D. Lee, and S. Yoon, “Bridging the nanogap with light: continuous tuning of plasmon coupling between gold nanoparticles,” ACS Nano **9**, 12292–12300 (2015). [CrossRef] [PubMed]

**15. **H. Cha, D. Lee, J. H. Yoon, and S. Yoon, “Plasmon coupling between silver nanoparticles: Transition from the classical to the quantum regime,” J. Colloid Interface Sci. **464**, 18–24 (2016). [CrossRef]

**16. **P. Zhang, J. Feist, A. Rubio, P. García-González, and F. J. García-Vidal, “Ab initio nanoplasmonics: The impact of atomic structure,” Phys. Rev. B **90**, 161407 (2014). [CrossRef]

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

**18. **U. Hohenester and C. Draxl, “Ab initio approach for gap plasmonics,” Phys. Rev. B **94**, 165418 (2016). [CrossRef]

**19. **W. Yan, M. Wubs, and N. A. Mortensen, “Projected dipole model for quantum plasmonics,” Phys. Rev. Lett. **115**, 137403 (2015). [CrossRef] [PubMed]

**20. **P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, and J. Luitz, *WIEN2k: An Augmented Plane Wave and Local Orbitals Program for Calculating Crystal Properties* (Vienna University of Technology, 2001).

**21. **C. Ambrosch-Draxl and J. O. Sofo, “Linear optical properties of solids within the full-potential linearized augmented planewave method,” Comput. Phys. Commun. **175**, 1–14 (2006). [CrossRef]

**22. **W. S. Werner, K. Glantschnig, and C. Ambrosch-Draxl, “Optical constants and inelastic electron-scattering data for 17 elemental metals,” J. Phys. Chem. Ref. Data , **38**, 1013–1092, (2009). [CrossRef]

**23. **S. Laref, J. Cao, A. Asaduzzaman, K. Runge, P. Deymier, R. W. Ziolkowski, M. Miyawaki, and K. Muralidharan, “Size-dependent permittivity and intrinsic optical anisotropy of nanometric gold thin films: a density functional theory study,” Opt. Express **21**, 11827–11838 (2013). [CrossRef] [PubMed]

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

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

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

**27. **I. E. Uysal, H. A. Ulku, and H. Bagci, “Transient analysis of electromagnetic wave interactions on plasmonic nanostructures using a surface integral equation solver,” J. Opt. Soc. Am. A **33**, 1747–1759 (2016). [CrossRef]

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

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

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

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

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

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

**34. **P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. **136**, B864 (1964). [CrossRef]

**35. **W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. **140**, A1133 (1965). [CrossRef]

**36. **S. Cottenier, “Density functional theory and the family of (L)APW-methods: a step-by-step introduction,” Instituut voor Kern-en Stralingsfysica, KU Leuven, Belgium **4**, 41 (2002).

**37. **P. Blaha, K. Schwarz, P. Sorantin, and S. Trickey, “Full-potential, linearized augmented plane wave programs for crystalline systems,” Comput. Phys. Commun. **59**, 399–415 (1990). [CrossRef]

**38. **J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. **77**, 3865–3868 (1996). [CrossRef] [PubMed]

**39. **F. Wooten, *Optical Properties of Solids* (Academic, 1972).

**40. **L. Brekhovskikh, *Waves in Layered Media* (Academic, 1980), Chap. 1, Sec. 13, pp. 98–99.

**41. **A. T. Fromhold, *Quantum Mechanics for Applied Physics and Engineering* (Academic, 1981).

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

**43. **G.-G. Siu and L. Cheng, “Mie solution of light scattering from spheres of radii up to 80λ with digit-array method,” J. Opt. Soc. Am. B **19**, 1922–1929 (2002). [CrossRef]