Quantum plasmonics: from jellium models to ab initio calculations

1 Introduction

Localized surface plasmons (LSPs) are collective oscillations of conduction electrons around the surface of a metallic object that arise as a result of a resonant coupling with an external incident electromagnetic (EM) field [1–4]. Such oscillations are driven by the force exerted by the incident field, whereas the restoring force depends very sensitively on electron-electron interactions and on the curvature of the surface. Then, the frequency of the LSP mode will be determined by the strength of electrostatic interactions and by the geometrical details of the metallic particle. It is thus possible to design and fabricate plasmonic nano-structures with a range of LSP resonance frequencies that can be determined virtually at will [5]. This tunability is at the heart of plasmonic-based nanoantennas [6, 7] or light-harvesting [8] and sensing devices [9]. In addition, the strong induced charge oscillations emit EM fields that are enhanced and concentrated in subwavelength regions [10, 11], which is the basis of plasmon-enhanced spectrosco-pies [12, 13]. Plasmon oscillations are also damped through non-radiative channels and, for instance, the generated heat upon LSP excitation can be used in cancer therapy [14, 15]. Furthermore, LSPs can decay by exciting energetic or “hot” electrons [16, 17], which is a useful mechanism in a number of applications such as photocatalysis or the design of novel photovoltaic devices [18].

In many cases of interest, the plasmonic properties of a nanostructure can be theoretically characterized just by solving the classical Maxwell’s equations with appropriate boundary conditions. As long as the sizes of the different constituents of the system are large enough, the EM response of each particle can be described by a local permittivity, ϵ(ω). The conceptual simplicity of this local-optics approximation together with the existence of a number of efficient numerical implementations [3, 19, 20] explains its widespread use in theoretical classical plasmonics.

For nanodevices where the separation between metallic constituents [21, 22] and/or their radii of curvature is smaller than a few nanometers [23], the local-optics approximation breaks down. Then, the intrinsic spatially non-local nature of the optical response has to be incorporated into the classical theory. This can be done either by using simplified hydrodynamic models [24, 25] or by defining space-dependent dielectric functions in the metal-dielectric interfaces [26]. Even if such a non-locality is included in a very crude manner, many deficiencies of the standard local-optics approximation are remedied. For instance, unlike the local prescription, these classicalnon-local optics approaches are able to describe the well-known [27] confined longitudinal plasmon resonances [28, 29] and do not show any pathological behavior in the response of two metallic objects in the touching-contact limit [22]. However, the unrealistic sharp boundaries between a metallic object and the surrounding dielectric medium are still kept. That is, classical non-local optics neglect the inhomogeneity of the electron density at the surface of a metal (the so-called electron-density spill-out). As a consequence, features like size-dependent shifts of the main LSP of isolated small nanoparticles [30] or the existence of further surface collective modes that cannot be supported on stepped conduction-electron distributions [31] are beyond the scope of classical non-local theories. Nonetheless, in many practical applications, the above-mentioned limitations are not relevant, so nonlocal optics can be considered as state of the art in classical plasmonics.

However, there are recent experiments on the optical properties of two metallic objects at subnanometric separation that dramatically differ from the predictions of both local- and non-local optics [32–38]. In this regime, the optical response is affected by the overlap of electron densities of the two nanoparticles and by the establishment of a photoinduced tunnel current between the metals [39–41]. The latter can be treated, even in the realm of local optics, by including effective dielectric media in what are called quantum-corrected classical methods [41]; however, the lack of electron-density spill-out is a serious drawback. For such systems, classical methods are thus inapplicable, which has contributed to the emergence of the field of quantum nanoplasmonics [42]. As a consequence, methods like time-dependent density functional theory (TDDFT) [43], hitherto restricted to condensed matter physics and quantum chemistry, are gaining more and more importance in the field of plasmonics.

The outline of this review is as follows. In the next section, we shall briefly describe the theoretical foundations of TDDFT and illustrate its applicability by analyzing the optical absorption of simple model systems. Then, we will summarize our recent research about the impact of the atomic structure in the optical properties of plasmonic nanostructures with subnanometric gaps. The corresponding summary and outlook will close the review.

2 TDDFT and simple-metal nanostructures

2.2 Spherical nanoparticles

Focusing on nanoparticles with spherical symmetry, it is fairly easy to calculate and analyze the optical absorption of simple sp-metal nanoparticles described by the spherical jellium model [71–75]. In general, good agreements with experiments are obtained after replacing the term i0⁺ of Eq. (15) by a finite contribution iη with ηη0.1 eV aimed at mimicking non-electronic decay channels (electron-phonon interaction, for instance) and shape/ size dispersion in experimental samples [1, 88, 89]. Moreover, the dynamical screening by d electrons, which plays a prominent role in the optical absorption of noble-metal particles [30, 54, 90, 91], can be simulated in spherical-jellium calculations by a dielectric background with an appropriate permittivity ϵ_b(ω) [cfr. Eq. (26)] [77]. Finally, the combination of effective dielectric media and suitable confining potentials can be used to get insights on the optical properties of doped semiconductor nanocrystals while keeping the simplicity of the spherical jellium model [92].

The TDDFT evaluation of optical absorption spectra of metal nanoparticles containing several thousands of conduction electrons can also be addressed at an affordable computational cost as long as the spherical symmetry is preserved [76, 78–82]. This allows us, for instance, to compare the optical absorption of isolated metal nano-particles at different size regimes. An example is given in Figure 1, where we depict the ALDA dipole photoabsorption spectra of Na_N spherical jellium clusters with N=92, 1000, and 5000 atoms. The mean density of the jellium positive background is nb=3/(4πrs3) with r_s=3.93 Bohr. Therefore, the jellium sphere radii are R=0.94, 2.08, and 3.55=nm for the three clusters, respectively. We have set the damping frequency η to 10 meV for the evaluation of the Green functions. This value is small enough to observe the contribution of individual single-electron excitations to the optical absorption. Finally, an artificial temperature smearing of 50 K has been included in the ground-state KS calculations of the two largest clusters to speed up their convergence. For each frequency of the external E-field, E_ext(r, t)=E₀ exp(-iωt)e_x [i.e. δv_ext(r, ω)=xE₀], we evaluate the induced density using Eq. (16). The dynamical polariz-ability is then given by the induced electric dipole divided by the amplitude of the incident E-field:

(19)α(ω)=−1E0∫xδn(r,ω)dr,

whereas the photoabsorption cross-section is given by

(20)σabs(ω)=4πωcℑα(ω).

Figure 1

ALDA absorption cross section in arbitrary units for spherical NaN jellium clusters of radius R. Red line: N=92 (R=0.94 nm); blue line: N=1000 (R=2.08 nm); black line: N= 5000 (R=3.55 nm). The bulk plasmon frequency, ωp, and classical Mie dipole surface plasmon frequency ωM=ωp/3, are indicated by vertical dashed lines. Note the logarithmic scale of the y-axis.

The general features of the spectrum are very well known [71–73, 88, 89], although it is worthwhile to summarize them. For the smallest cluster (N=92), the presence of the LSP around ω=3 eV is evident, albeit fragmented due to the interaction with surrounding single-electron excitations [93]. In the region around ω=4 eV, a second spectral feature starts to develop, which is a signature of a second surface mode (a Bennett plasmon) resulting from the inhomogeneity of the density profile at the interface [31]. Finally, above the Na volume plasmon frequency ω_p≃6.04 eV, a Landau-damped mode clearly emerges in the spectrum: this is a confined longitudinal volume plasmon. The LSP main resonance is red-shifted with respect to the classical value ωM=ωP/3≃3.49 eV simply because the confining potential at the interface is softer than in the classical picture, thus slowing down the oscillations of the electrons. On the contrary, the confined volume mode is blue-shifted with respect to ω_p due to level quantization. That is, for simple sp-metal clusters, there are two competing mechanisms: diffuseness of the confining potential (which dominates for the surface mode) and discreteness of the one-electron energy levels (which is the most important feature for bulk-like modes) [94]. For the medium-size cluster (N=1000), the LSP frequency approaches to the classical value ω_M and is less fragmented; the Bennett mode is slightly more defined; and not only one but a series of confined volume plasmons appears in the spectrum. Finally, for the largest nanoparticle studied here (N=5000), the one-electron energy levels are very close to be continuous. The LSP comprises almost all the spectral weight, although it still shows a certain structure and its frequency (~3.3 eV) is smaller than ω_M; the Bennett mode appears as a clear shoulder around ω=4.5 eV; and the frequency separation between two consecutive confined volume plasmons is smaller than in Na₁₀₀₀, as expected from the larger size of the cluster. Interestingly, the first confined volume mode appears just belowω_p. As for this size the level quantization should not play any role, the diffuseness of the confining potential is the most plausible explanation of this red shift with respect to the classical value.

2.4 Nanowire dimers

As a first example of the application of time-propagation TDDFT in model plasmonic nanostructures, let us consider the dipole absorption spectrum for a jellium-Na nanowire of radius R=2 nm and for a nanowire dimer where the separation between jellium boundaries is d=0.5 nm [103]. The incident E-field is perpendicular to the nanowire axis and, in the case of the dimer, it is also parallel to the plane determined by the two cylinder axes (see Figure 2). Unlike other similar studies [97, 98], the only symmetry that is taken into account is the translational invariance along the cylinder axis, which is implemented by imposing Born-von Kármán periodic boundary conditions. The TDDFT results are compared with local and non-local classical predictions, both evaluated using the finite element COMSOL package [104]. For the local calculations, ϵ(ω) is a Drude dielectric function

(26)ϵ(ω)=1−ωP2ω(ω+iγω),

where ω_p=6.04 eV is the bulk Na plasmon frequency and γ is a phenomenological damping frequency that accounts for non-radiative decay channels for the conduction electrons (for instance, surface scattering) [1]. For the present case, we have chosen γ=0.38 eV. The non-local hydrodynamics calculations are implemented under the prescriptions by Raza et al. and Hire-math et al. [29, 105, 106]. A damping η=0.1 eV is used for the TDDFT calculations. We must note that the distance d=0.5 nm between wires is close to the onset of the establishment of an induced tunnel current between the wires. As mentioned in Section 1, such a current drastically changes the absorption properties of the system. Therefore, we are at the limit of applicability of classical-optics prescriptions.

Figure 2

Left panel: Optical absorption spectrum of a sodium nanowire (R=2 nm). Right panel: Optical absorption of a sodium nanowire dimer (R=2 nm, d=0.5 nm). Thick black line: TDDFT; thin blue line: non-local optics; dashed red line: local optics. Both classical calculations have been performed with a damping frequency γ=0.38 eV. Figure adapted from Ref. [103].

As can be seen in the left panel of Figure 2, the TDDFT absorption of the single nanowire is dominated by a LSP appearing at ω≃4.09 eV. As expected, this value is red-shifted with respect to the classical counterpart ωs=ωP/2≃4.17 eV. The Bennett plasmon appears as a shoulder around ω≃4.6 eV, although it is not fully developed in line with the results previously shown for the isolated jellium-Na nanoparticle. Finally, albeit not presented in Figure 2, a series of confined bulk longitudinal plasmons can be observed in the range ω=6 eV [103]. The local-optics absorption consists in a single peak centered at ω_S=4.17 eV corresponding to the main LSP. Such a peak is blue-shifted when non-local effects are included, which is the expected result as this standard non-local optics calculation does not describe the electron-density spill-out.

Therefore, the local optics spectrum is closer to the TDDFT one due to a cancellation of errors: lack of non-locality and of density spill-out.

The case of the nanowire dimer is more interesting. First, the geometry of the system allows for a very intense enhancement of the E-field in the region between the wires. Second, the absorption spectrum reflects the well-known hybridization of plasmon modes [11, 107, 108]. Higher multipole-order modes in one wire, which cannot be excited by a homogeneous electric field, are excited by the near induced E-field emitted by the other wire and vice versa. This new mode hybridizes with the main dipole plasmon, thus giving a new resonance, which we label as “Q”. In addition, the dipole surface plasmons of both wires are coupled and red-shifted with respect to the frequency for an isolated wire, a mode that is noted as “D”. It must be mentioned that according to the classical optics picture, the multipole modes of a cylinder are degenerate at a frequency ω_S [47]; that is, all the hybridized modes “depart” from this frequency. Therefore, the spectral feature “Q” might actually contain contributions from higher-order hybridized modes that cannot be resolved in the spectra [97].

As can be observed in the right panel of Figure 2, the description of the absorption spectrum under the three approaches agrees well with the hybridization picture; however, quantitative differences are apparent. Although the distance between the D and Q resonances is practically the same for the three approaches, the relative spectral weight of each resonance as well as the absolute position are very different. Furthermore, the local spectrum shows the prospect of a second hybridization splitting. These discrepancies can be easily traced back to the very different shapes of the induced electron density δn(r, ω). Whereas the quantum TDDFT predicts that it is mainly located outside the jellium boundary, the induced density is concentrated on the surface under the local optics prescription, and δn(r, ω) is spread toward the bulk in the non-local optics approximation. Therefore, the effective electrostatic interaction between the wires is strongest in the quantum picture, whereas it is weakest under non-local optics. Due to the overall dipole symmetry of the induced density along the direction of the incident field, the near E-field emitted by one wire is against the restoring force acting over the induced density in the other wire. Then, the larger the interaction between wires, the less intense the effective restoring force. As a consequence, the red shift of the whole spectrum is strongest in the quantum TDDFT picture. Although the red-shift amount can be easily correlated to the distance between the effective surfaces in both dimers [97], it is evident that in plasmonic nanostructures with narrow gaps we can only achieve full predictive agreement if the inhomogeneity of the induced density is described properly.

The different shapes of the induced densities will be also reflected on the spatial distribution of the induced electric field enhancement (EFE), defined as the ratio between the amplitudes of the total and the external E-fields. In the D mode, where the induced E-field is concentrated in the gap region, the overall shape predicted by TDDFT is rather different than the calculated using classical methods (see Figure 3). This is due to the influence of the electron-density spill out. The differences are even larger when considering the Q mode. In this case, there is not even a qualitative agreement between the classical and quantum calculations. This is not a surprise because the higher-order terms in the multipole expansion of the induced density at each wire are more important in the Q mode than in the D resonance. Therefore, the near E-field in the Q mode is more sensitive to variations in the shape of the corresponding source. This further confirms that in the limit of subnanometric separations, the inhomogeneity of the density has to be determined accurately when analyzing the plasmonic response of these systems.

Figure 3

EFE |E_tot(r)|/|E_ext|corresponding to the dipole (top line) and quadrupole (lower line) resonances of a Na nanowire dimer (d=2 nm). The depicted E-fields have been obtained at the frequencies indicated in each panel. Figure adapted from Ref. [103].

3 Ab initio nanoplasmonics: Na-cluster dimers

As we have seen in the previous section, the plasmonic response of nanodevices with subnanometric gaps depends very sensitively on the amount of electron-density spill-out in the gap region. Furthermore, as commented in Section 1, the establishment of a photoinduced alternating tunnel current between the metal objects can change even more the absorption properties. Using TDDFT and the jellium model, these effects have been studied in detail in the last years, so the main trends are presently well understood [39, 40, 97, 98, 103]. However, as subnanometric length scales are relevant in these nanodevices, it is expected that the atomic structure around the gap region will have at least quantitative impact. Bearing in mind that one of the ultimate goals in theoretical nanoplasmonic should be the development of robust multiscale methods amenable for addressing both classical and quantum aspects of the optical response [109–111], the quantification of the relevance of the atomic structure in prototypical plasmonic nanodevices is of major importance.

Of course, the atomic structure is routinely taken into account in the analysis and prediction of the properties of metallic clusters either from a quantum-chemistry or a condensed matter perspective [112–115]. Moreover, in the last decade, there has been an increasing interest in the ab initio characterization of plasmonic properties of clusters and macromolecules containing up to several hundred atoms. For instance, the plasmon tunability of metallic nanorods and nanowires [116–120], the size-dependent development of collective excitations [121–129], alloys [130–132], core-shell and hybrid nanostructures [133–138], molecular plasmons [55, 139], and small-cluster dimers [140], have been analyzed using TDDFT-based methods. While these studies provide very useful insights about optical absorption properties, only very recently the ab initio characterization of other relevant aspects in nanoplasmonics has started to be addressed [141–145].

Sodium-cluster dimers constitute an appropriate play field to address the question about the impact of the atomic structure in nanoplasmonics. Unlike noble metals, the core electrons do not need to be explicitly considered to evaluate the optical response. Therefore, standard norm-conserving pseudopotentials [146] together with the LDA or GGA level of theory can be safely used. By varying the cluster-cluster distance, d, different plasmonic regimes can be reached and, therefore, it is possible to carry out a systematic study of the importance of the atomic structure at an acceptable computational cost. Specifically, the calculations that we present in this section [141, 142] were also performed with the Octopus package [99–102]. Typical convergence parameters are spatial-grid spacing of 0.026 nm, total propagation time T_max≃30 fs, time step Δt=0.002 fs, and damping frequency η=0.1 eV.

Before analyzing the optical properties of a Na-cluster plasmonic dimer, it is illustrative to discuss very briefly the optical absorption spectra of isolated Na clusters. We consider Na clusters with two different atomic structures: unrelaxed bcc-Na₃₃₁, where the atoms are arranged symmetrically around a central atom on the same bcc lattice as in bulk Na (lattice constant a=0.423 nm), and relaxed ico-Na₂₉₇ cluster, which exhibits icosahedral (I_h) symmetry and whose atomic positions are taken from a previous work [147]. Both clusters are almost spherical, and in principle the spherical jellium model is suitable for evaluating their optical properties. Nonetheless, one cannot expect a perfect agreement between the ab initio and the jellium-model-based results. The optical properties for the jellium systems are obtained from Octopus as well, although we have checked that for isolated clusters there are minor differences with respect to a homemade code that exploits the full spherical symmetry of those systems. A small temperature smearing has been included in all the calculations.

The comparison of the absorption spectra for the atomistic and the jellium clusters is depicted in Figure 4. As it can be seen, the atom- and jellium-TDDFT spectra are fairly similar, although the latter are slightly blue-shifted (<0.2 eV). Also note that the whole spectrum of the ico-Na₂₉₇ is blue-shifted with respect to the bcc-Na₃₃₁. Although a direct comparison is not possible due to the different atomic arrangements, the icosahedral cluster is slightly compressed with respect to the unrelaxed bcc cluster, an effect that we have taken into account in the corresponding jellium calculation. Then, the blue shift of the smallest cluster can be primarily attributed to an overall compression resulting from the relaxation of the structure. As expected from the discussion in the previous section, the four clusters exhibit a spectral feature in the region 4.0–4.5 eV, albeit almost unresolved for the jellium nanoparticles. Nevertheless, none of these qualitative discrepancies are a surprise. In fact, the importance of the atomic structure is manifest when comparing the EFEs. As we may see in Figure 5, the ionic distribution affects the shape of the induced E-field and also its absolute value. Under the jellium model (in what follows, we will only consider the jellium-Na₃₃₁ case), the EFE is overestimated and the tiny structures that reflect the surface corrugation at atomic scale are obviously smoothened out. A more detailed discussion about the inhomogeneities of the induced E-field at atomic scale can be found in the recent work by Barbry et al. [143].

Figure 4

Normalized TDDFT absorption cross section for Na isolated cluster. Solid black line: BCC Na₃₃₁; solid red line: ICO Na₂₉₇; black dotted line: jellium Na₃₃₁; red dotted line: jellium Na₂₉₇. The classical Mie plasmon frequency for the mean density of bulk Na is indicated by a vertical line. Figure adapted from Ref. [141].

Figure 5

Contour plots of the EFE for the BCC Na₃₃₁, ICO Na₂₉₇, and jellium Na₃₃₁ clusters at the corresponding Mie LSP resonance. The E-field amplitudes are normalized to their maximum value, which is indicated in each panel. The incident E-field is directed along the X axis.

Concerning the response properties of Na-clusters dimers, previous TDDFT studies on spherical-jellium dimers clearly suggest the existence of three distinct separation regimes when the incident E-field is parallel to the dimer axis [39, 40, 148]: (a) classical: at large enough separations, the interaction between clusters is mainly EM as there is a negligible photoinduced charge transfer (CT) between the clusters; (b) intermediate: a tunnelling CT between the clusters appears as a result of a photoinduced alternate voltage bias between the clusters (also note that the ground-state densities of the clusters slightly overlap in this regime); and (c) direct contact: the ground-state densities of the clusters overlap because the Fermi level of the system is above the potential barrier in the dimer junction.

The optical absorption for the three cluster dimers are presented in Figure 6 where we normalize the frequency ω of the incident EM field to the frequency of the Mie LSP ωM(t) for each prescription (see Figure 4). By doing so, we will focus on the differences arising from the interaction between the dimers. For d≳0.4 nm, the three dimers are in the regime (a), and the spectra reflects the hybridization of modes that we described in the previous section. However, there are some differences in the spectra, as the hybridization is strongest in the jellium-Na₃₃₁ dimer and weakest in the ico-Na₂₉₇, as can be seen in the upper panels of Figure 6. This can be attributed to the different surface shapes of the clusters that affects the EM near field, which, as we know, drives the hybridization process.

Figure 6

Upper panels: Waterfall representation of the normalized optical absorption cross section for the bcc-Na₃₃₁, ico-Na₂₉₇, and jellium-Na₃₃₁ dimers. The spectra have been shifted to ease the comparison between the different regimes. The separation d (in Å) is indicated in the left panel. Lower panels: Corresponding contour plots of the normalized absorption spectrum.

In the intermediate regime (b), the classical description breaks down; however, the transition from the classical to this intermediate regime is different in the three dimers (this can be clearly seen in the contour plots displayed the lower panels of Figure 6). For the jellium dimer, the regime (b) is reached at d≃0.3 nm. In this case, the D mode frequency stabilizes and the spectral weight of the mode decreases as the nanoparticles approach each other. The hybridized Q mode becomes dominant for d=0.2 nm and, finally, when d=0.1 nm, the mode D is completely quenched and a very weak signature of a CT mode appears instead [41]. In contrast, the intermediate regime is reached by the atomistic clusters at a larger separation (d>0.4 nm), and for d≃0.1 nm the CT mode appears as a strong spectral feature in the absorption spectra, this being the signature of the regime (c). It seems that this regime is not fully reached by the jellium cluster in the range of separations that we are considering.

The nanoparticle dimer is one of the prototypical systems to test light-harvesting properties of nanoplasmonic devices. Therefore, it is worth analyzing how the atomistic description affects both the modal shape and EFE associated with the different resonances of the cluster dimer. As we have shown in Figure 6, for d=0.3 nm, the optical response of the dimer is dominated by the dipolar mode. In the left panels of Figure 7, we render the electric field amplitudes for this D mode calculated at d=0.4 nm for ico-Na₂₉₇ and jellium-Na₃₃₁ dimers (the EFE for the bccNa₃₃₁ is very similar to that for the icosahedral-cluster dimer for the four modes depicted in the figure, and it is not included there for the sake of simplicity). As expected, the modal shape is very similar in the three cases, although the D mode for the atomistic cluster is a bit more delocalized than the jellium one. As a consequence, the EFE is reduced by a factor of around 1.5 when considering the atomic structure. For the quadrupole (Q) modes, the relative impact of the atomic structure is greater, and the overestimation of the EFE by the jellium mode persists. Finally, as can be seen in the right panels of Figure 7, the CT mode is not very strongly affected by the atomic structure.

Figure 7

Contour plots of the total electric field amplitude evaluated on the XY plane of the dimer at four selected resonances. Upper row: ico-Na₂₉₇; lower row: jellium-Na₃₃₁. From left to right: dipolar mode D (d=0.4 nm); quadrupole mode Q (d=0.2 nm), Q mode (d=0.1 nm), and CT mode (d=0.1 nm). In each panel, the E-field amplitude is normalized to its maximum value, which is also shown as a legend. Figure adapted from Ref. [141].

The very different behaviors in the intermediate regime can be easily traced back to the influence that the atomic structure has on the intensity of the photoinduced current between the particles. However, a direct comparison between the three cases is not entirely possible as there is no way to define a separation d in an unambiguous manner for the three geometries at the same time. A better scenario to quantify the impact of the atomic structure in this model nanoplasmonic structure is to consider only one dimer cluster but changing the relative orientation between the clusters while keeping unchanged the distance between the centers of the nanoparticles. Namely, if 2R is the effective diameter of one cluster (2.61 nm for the icosahedral 297-atom nanoparticle), the distance d between the clusters is d=b-2R, b being the distance between the center of the clusters (see Figure 8). As the most realistic individual structure is the ico-Na₂₉₇, in what follows we will restrict our analysis to this cluster. In addition, we will perform some simulations of the plasmon dynamics of the nanostructure [142].

Figure 8

Representation of the two different ico-Na₂₉₇ dimers considered in the present analysis. The applied electric field is orientated along the dimer axis. Figure adapted from Ref. [142].

We will address two main features in this analysis: the relative orientation of the clusters and the atomic relaxation due to the mutual interaction. Concerning the former, we treat the two cases represented in Figure 6. The first one leads to a spatial gap that is limited by two 12-atom faces (F2F orientation) and is the geometry that we used in the previous analysis. The second one corresponds to a spatial gap between two three-atom edges (E2E orientation). The first relative orientation maximized the width of the junction but corresponds to a larger separation between atoms. These two effects tend to cancel each other, thus having an stringent but well-defined scenario to quantify in detail the importance of the atomic structure. In addition, the position of the atoms in the junction region is relaxed in the E2E geometry due to the smaller distance between atoms. The energy optimization is performed using the FIRE algorithm [149], recently implemented in Octopus [102]. Once the structure is relaxed, the absorption spectra are calculated as in the previous analysis.

The dipole optical absorption in the range of separations 0≤d≤0.6 nm is presented in Figure 9. We have also included the optical absorption for the E2E overlapping clusters (d=-0.1 nm), where a major reconstruction of the junction occurs. As can be seen, there are marginal differences between the E2E and F2F optical spectra for the separation regimes (a) and (c), which correspond to the classical and overlapping limits. The irrelevance of the relative orientation in the classical limit is not really surprising due to the overall I_h symmetry of each cluster. In addition, in the overlapping regime (c), the two clusters cannot be considered as individual entities anymore and the system is actually a single “peanut-like” nanoparticle. In fair correspondence with the comparison that we have made for individual clusters, it is not expected that the different arrangements in the junction (which now is inside the nanoparticle) and at the surface will have a major impact in the absorption. As we have already seen, in this overlapping regime, the spectrum is dominated by two LSPs, labeled now as D(CT) and Q(CT), in order to distinguish their shape and emphasize that there is a CT between the two lobes of the peanut-shaped nanoparticle in both modes. Also note that these modes can be addressed using classical optics if the shape of the coalesced nanoparticle is properly modeled [150, 151].

Figure 9

Optical photoabsorption cross sections, σ_abs(ω), of Na₂₉₇ dimers as a function of the distance d as defined in the text. Upper panels: contour plots. The vertical dashed line marks the approximate distance where the hybridized Q mode becomes the main spectral feature in the spectrum. Lower panel: waterfall plot of the absorption cross sections. The red and blue lines correspond to the F2F and E2E orientation, respectively. For each distance, the separation in nm between nearest faces or edges is indicated between parenthesis. Figure adapted from Ref. [142].

On the contrary, the relative orientation in the intermediate regime (b) is crucial. For this range of distances, the optical absorption of the dimer results from a competition between near induced E-fields and the establishment of a tunneling CT between the nanoparticles. It is evident from Figure 9 that the F2F dimer is not yet in this regime when d~0.4 nm. However, for the E2E dimer, the value of the frequency of the D mode is already locked and there is a net transfer of spectral weight from the D mode to the Q mode, which is a clear signature of the breakdown of the classical picture. In fact, the quantum effects are only discernible in the F2F spectrum when d<0.4 nm. At d<0.3 nm, the two spectra are qualitatively different, as in the E2E dimer the D mode has been already quenched, whereas in the F2F dimer the dipole mode still dominates the spectrum.

Let us remember that the atomic positions of the E2E geometry are relaxed. In fact, the cluster-cluster interaction leading to a further relaxation of the structures with respect to the original energy-optimized ones is only significant if d≲0.4 nm for the E2E geometry and if d≲0.25 nm for the F2F configuration. For shorter distances, the atomic reconstruction will be important; however, in this regime, the specific ionic arrangement is of secondary importance. Therefore, the atomic relaxation will play a role only in the range of those critical distances where the spectrum changes its overall structure; that is, when the Q mode starts to be dominant. By comparing the absorption spectra depicted in Figure 9 with that obtained without relaxing the atoms in the gap region, we indeed observe a change in the optical spectrum [142]. However, the changes are merely quantitative and less important than those related to the relative orientation.

In Ref. [142], a detailed study of the time evolution of driven induced current intensities and EFEs was presented. In this review, we present some new results that complement those already discussed in Ref. [142]. Specifically, we consider the action over the system of a weak uniform quasi-monochromatic laser pulse of mean frequency ω_ext and duration τ=20π/ω_ext (i.e. 10 periods) given by

where E₀=0.51×10⁶ V m^-1. Note that the maximum intensity of this laser pulse is 34.5 kW cm^-2, well below the onset of non-linear effects, already analyzed under the jellium model prescription by Marinica and coworkers [40].

As an example of time-resolved dynamics that can be addressed using TDDFT, in Figure 10 we depict the evolution of the x component of the total E-field in the center of the dimer junction driven by the external field [Eq. (27)] and the driven electric dipole,

for selected resonant couplings in the E2E nanodimer. For all the cases, as corresponds to the resonance condition, there is a quarter-period phase shift between the external field (represented in the top panel of Figure 10) and the dipole. Moreover, the maximum amplitude of the dipole oscillations are delayed with respect to the maximum intensity of the external field. On the other hand, the different phase shifts between the driven dipole and the E-field can be readily observed. Note that the D mode at d=0.5 nm is the narrowest resonance among all those presented in Figure 10 (see the absorption spectra in Figure 9).

Figure 10

Time evolution of total electric fields (solid blue lines) at the middle of the dimer junction (r=0) and induced dipoles (dashed green lines) due to the resonant coupling of selected modes of the E2E nanoparticle dimer with an applied external monochromatic field (depicted in the upper panel) of frequency ω_ext and maximum amplitude E₀=10^-6 a.u.

Figure 11

Time evolution of induced currents (dashed green lines) and total electric fields (solid blue lines) at the middle of the dimer junction (r=0) due to the resonant coupling of selected modes of the F2F nanoparticle dimer with an applied external monochromatic field (depicted in the upper panel) of frequency ω_ext and maximum amplitude E₀=10^-6 a.u. The subsequent free evolution during a time interval ΔT=4π/ω_ext after the disappearance of the external field is also displayed.

Therefore, when the driving external field starts to fade out, the phase shifts and the relative amplitude between the total E-field and the dipole remain practically the same. However, the relative evolution of these two quantities changes appreciably for the rest of the modes when the intensity of the external perturbation decays.

This effect is easier to observe after the total disappearance of the external field, as the motion of the electrons has to be self-sustained (up to the decay due to the energy transfer from the coherent plasmon mode to an incoherent set of electron-hole excitations), and thus is closer to the actual mode. In Figure 11, we represent the time-resolved evolution of the E-field in the center of the dimer and the current intensity between the nanoparticles for the D modes at d=0.5 nm and d=0.3 nm of the F2F dimer. In both cases, the CT is established via tunneling through the potential barrier between the clusters. Then, the junction behaves in principle as an ohmic resistor [40], which explains the negligible phase shift between the current and the total E-field as well as the practically constant ratio between those quantities. However, upon the disappearance of the driving field, the free evolutions of the induced current and the E-field are not exactly the same when d=0.3 nm, whereas the free evolution seems to be unaffected when d=0.5 nm. The former simply indicates that the self-sustained photoinduced current when d=0.3 nm has to be characterized by an inductive impedance, as it is the case of the modes in the direct-contact regime (c).

4 Outlook and conclusions

TDDFT has been a method of choice in ab initio condensed matter and cluster physics for more than two decades. However, its application in nanoplasmonics is far more recent. Quantum simulations not only allow us to explore physical phenomena that cannot be addressed under the classical description of light-matter interaction, but are also an invaluable tool to assess approximate methods. Model calculations, like those discussed in the first part of this review, can be employed for a critical judgment of new theoretical developments with a broader scope of applicability than TDDFT.

We cannot forget that the nanodevices with actual technological interest are much larger than those systems that can be treated in a fully quantum manner. Then, simplified descriptions of the electron motion will be certainly required. For instance, as the limitations of the HA are well known [54], the mere reproduction of the main spectral features of an isolated metallic object in the dipole quasistatic limit is not enough. It would be interesting to perform more stringent assessments of its predictive accuracy, and the model calculations discussed in this article are indeed a good test. Also note that the HA reduces the many-body problem to a single-body Schrödinger equation, the wavefunction being the square root of the density, so its usefulness in the simulation of dynamical phenomena can be assessed as well.

We have also presented a detailed study of the importance of the atomic structure, so far overlooked in many quantum simulations of the plasmonic response of nanodevices. We have focused on the optical response of a prototypical plasmonic nanostructure and seen that the atomic structure becomes critical in the tunneling CT regime between the constituent nanoparticles. In addition, effects like atomic relaxation and different relative orientations have been also explicitly addressed. On the other hand, our findings confirm the good qualitative accuracy of the jellium model, either combined with quantum mechanical calculations or classical theories, except when there are very narrow gaps in the plasmonic nanostructure. As we have seen, this is precisely the regime where the most complicated processes occur. Thus, we have shown that a fair account of the atomic structure of the junction for very close nanoparticles needs to be taken into account to obtain theoretical results with enough predictive power.

Finally, our calculations of photoinduced currents between nanoparticles constitute a first step toward ab initio optoelectronics but at the TDDFT level of theory, which incorporates exchange and correlation quantum many-body effects [152, 153]. In addition, the explicit time propagation of the KS equations enables the time-resolved analysis of the dynamics of the collective motion of electrons driven by an arbitrary external E-field. Furthermore, it is possible to simulate the transient regime in the resonant coupling between the nanoparticle and an external quasimonochromatic perturbation, as well as the free evolution of the electron system upon the fade out of the external perturbation. The latter will give us, for instance, first-hand information about plasmon decay via generation of electron-hole pairs. Finally, the frozenion approximation, which we have adopted throughout this study, can be relaxed at an affordable computational cost. This paves the way to the semiclassical simulation of other decay channels that are absent in the, until now, few TDDFT simulations in nanoplasmonic devices.