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

 

Fig. 1 Summary of the proposed approach.

Download Full Size | PDF

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).

 

Fig. 2 Auxiliary tunnel with (a) single layer and (b) three layers.

Download Full Size | PDF

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]

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]

 

Fig. 3 Supercell with gap.

Download Full Size | PDF

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, µ₀ 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 ∂_RG(R, ω). It is assumed that in a given frequency band F(ω) can be approximated using rational functions as

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₀. The permittivity and inverse permittivity of V_p are denoted by ε_p (t) and ${\overline{\epsilon }}_p(t)$, respectively. For p = 0, permittivity is constant, i.e., ε₀(t) = ε₀ and ${\overline{\epsilon }}_0(t)=1/{\epsilon }_0$. All volumes are non-magnetic with the constant permeability µ_p = µ₀, p = 0, …, P. The surface between the two volumes V_q and V_p is represented by 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_l pointing towards V_p. 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 V_p, 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 ∀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:

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 ε₀(t) = ε₀, ${\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_l are discretized into triangular patches and unknowns 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), ∂_RG_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_mtK_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_n 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 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⁻⁴ Ry.

In all of the examples considered here, it is assumed that the nanostructures reside in free space with permittivity ϵ₀ = 8.8541878176 × 10⁻¹² F/m and permeability μ₀ = 4π × 10⁻⁷ H/m and are excited by a plane wave with electric field

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.

 

Fig. 4 (a) Real and (b) imaginary parts of ε_gap(ω) for gold.

Download Full Size | PDF

 

Fig. 5 (a) Real and (b) imaginary parts of ε_gap(ω) for silver.

Download Full Size | PDF

 

Fig. 6 (a) Real and (b) imaginary parts of ε(ω) for bulk gold.

Download Full Size | PDF

 

Fig. 7 (a) Real and (b) imaginary parts of ε(ω) for bulk silver.

Download Full Size | PDF

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

 

Fig. 8 C_ext(ω) computed using the Mie series solution for (a) a gold sphere of radius 25 nm and (b) a silver sphere of radius 10 nm.

Download Full Size | PDF

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₀ = 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.

 

Fig. 9 C_ext(ω) computed using the TD-SIE solver and the Mie series solution for (a) a gold sphere of radius 25 nm and (b) a silver sphere of radius 10 nm.

Download Full Size | PDF

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₀ = 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).

 

Fig. 10 Electric and magnetic current densities computed using the TD-SIE solver on the surface of the gold dimer (a) without the auxiliary tunnel and (b) with the auxiliary tunnel.

Download Full Size | PDF

 

Fig. 11 C_ext(ω) computed using the TD-SIE solver for the gold dimer (a) without the auxiliary tunnel and (b) with the auxiliary tunnel.

Download Full Size | PDF

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 ϵ₀ 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.

 

Fig. 12 C_ext(ω) computed using the TD-SIE solver for the gold dimer with the auxiliary tunnel constructed using different number of layers.

Download Full Size | PDF

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₀ = 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].

 

Fig. 13 C_ext(ω) computed using the TD-SIE solver for the silver dimer (a) without the auxiliary tunnel and (b) with the auxiliary tunnel.

Download Full Size | PDF

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.