Modeling of plasmonic properties of nanostructures for next generation solar cells and beyond

ABSTRACT

Plasmonic particles and nanostructures are widely used in photovoltaic and photonics. Surface plasmons were found to enhance different types of solar cells including plasmonic DSSCs, plasmonic solid semiconductor solar cells, plasmonic organic solar cells, and plasmonic perovskite solar cell. Size, composition, and shape of plasmonic nanoparticles as well as nanometer-distance control between particles are key design factors of plasmonic nanostructures. Modeling is rapidly gaining in importance for mechanistic understanding and rational design of plasmonic nanostructures. We review the modeling approaches used to model plasmon resonance features of nanostructures, from classical approaches that can routinely handle most particle sizes used in solar cells to approaches beyond classical electrodynamics such as ab initio approaches based on time-dependent density functional theory (TD-DFT). We highlight recently emerging approaches which have the potential to significantly enhance modeling capabilities in the coming years, in particular, by allowing atomistic (ab initio) modeling at realistic length scales, i.e. of particle sizes beyond 10 nm which are of most interest to plasmonic solar cells but remain problematic with traditional DFT-based techniques, such as density functional tight binding (DFTB) based approaches, time-dependent orbital-free DFT, and machine learning-based approaches, as well as many-body perturbation theory which is expected to gain usage with advances in computing power.

Graphical Abstract

KEYWORDS:

• Surface plasmon
• plasmonic solar cell
• time-dependent density functional theory
• time-dependent density functional tight binding
• time-dependent orbital-free dft

Introduction

When metal nanoparticles (NPs), such as silver (Ag) and gold (Au), as well as some doped semiconductor nanoparticles are irradiated by light, their localized surface plasmons, i.e. collective oscillations of charge density, are excited, thus causing light scattering and light absorption. The classical theory established early on that for metal particles with a diameter below about 40 nm, the energy absorbed is used for the excitation of surface plasmons (SPs), which is a predominant process [1,2]. On the other hand, for particles with a diameter above 100 nm, the scattering effects dominate. Particles of such large sizes are therefore often used for light management e.g. in dye-sensitized solar cells (DSSC) and other types of solar cells. Both absorption and scattering are present in NPs of intermediate sizes. The effects of interest to us here are produced by particles of sizes up to about 100 nm, which, as we will see below, are already too large for direct atomistic (ab initio) modeling, even though the smallest sized NPs used in applications (on the order of a nm) can be directly modeled with the mainstream ab initio methods and tools. The number and wavelength of surface plasmon resonances (SPRs) determines the optical spectrum peaks and their intensity and width. The peaks can be tuned by modifying the size, morphology, chemical composition, and other properties.

In recent years, plasmonic phenomena found increasing use in photovoltaics (PV) and other applications. In 1997, the concept of a DSSC that utilizes localized surface plasmons was introduced [3]. In this DSSC, near-field light was used to enhance both the photoabsorption of Ru-dye and the generation of photoelectrons. By introducing AgNPs into a film of N3 dye, the absorbance of the dye was increased by a factor of about 150. Absorption enhancement has also been reported when the composite films are formed with AuNPs and black dye (BD). A key factor for enhancing the photoabsorption of a plasmonic DSSC with such a composite film is an overlap in the photoabsorption spectra between the metal NPs and the dye [4]. The NP spectrum, specifically the plasmonic resonance peak, is thus a key design criterion. Already first DSSCs with Ag NPs modified with polyacrylate-based comb-shaped block copolymer showed improved photoelectric conversion efficiency (PCE) of 2.5%, compared with about 1.5% for a DSSC without Ag NPs[5]. In the following years, plasmonic solar cells have been extensively researched and efficiency improvements shown in different types of cells beyond DSSC were reported [3–18]: plasmonic solid semiconductor solar cells [19–38], plasmonic organic solar cells [39–45], plasmonic perovskite solar cells [46–53] and other types of plasmonic solar cells [54–59] based on the features of their localized surface plasmons and surface plasmon polaritons [60].

Plasmons are the collective oscillation of free electrons in a material. As shown in Figure 1, plasmons can be roughly classified into surface plasmons and plasmons in bulk. Surface plasmons have been extensively researched, including the basic theory of plasmons and the various approaches to applications such as surface-enhanced Raman scattering, optical near-field optical microscopy, and nano-manufacturing [61,62]. Surface plasmons can be classified further into surface plasmon polaritons (SPPs) and localized surface plasmons (LSPs), as schematically shown in Figures 2 and Figure 3, respectively.

Figure 1. Classification of plasmons

Figure 2. Concept of surface plasmon polariton (SPP)

Figure 3. Concept of localized surface plasmon

A plasmon is a charge density wave, which is a type of longitudinal wave. In general, when a wave is excited by another wave, both waves must match in both frequency and wave number vector. Because light such as sunlight is a type of transverse electromagnetic wave, the longitudinal wave of a plasmon in bulk usually cannot be excited by using the transverse wave of light due to the dispersion relations.

However, on a metal surface in the Kretschmann configuration in which a prism is used, SPPs can be excited by such light due to the different boundary conditions at the surface from those at the bulk [63]. This configuration enables a special dispersion relation of the light. The excitation of SPP produces a transverse magnetic wave (TM wave), which can propagate along the surface of a metal. A feature of a TM wave is that its oscillation is localized near the surface and its amplitude is attenuated exponentially with increasing distance from the surface. Therefore, the strong electric field called an evanescent wave, sometimes a hundred times stronger than that of irradiated light, can be generated near the surface of metals by the resonant excitation of SPPs.

In contrast, LSPs can be excited by light such as sunlight without a special configuration such as the Kretschmann configuration. Compared with SPPs, this is an advantage when LSPs are applied to solar cells. LSPs are excited in metal NPs whose diameter, usually ranging from several nanometers to several tens of nanometers, is much smaller than the wavelength of irradiating light. In that case, NPs are effectively exposed to irradiation of a longitudinal wave due to negligibly small effect of the retardation originating from the finite velocity of light, although a transverse wave of sunlight actually irradiates NPs. Therefore, oscillations of an electric dipole can be induced at metal NPs. These oscillations produce a strong electric field at the surface of the metal NPs, called optical near-field whose electric field, which is sometimes a hundred times stronger than that of the irradiating light. The optical near-field is not propagated, but is localized within an area about an NP diameter from the surface of metal NP.

Two types of applications of surface plasmons to photoelectric conversion have been reported. The first application uses surface plasmons to enhance the photoabsorption of original photoabsorbing materials of solar cells by the two effects described below. This type of solar cell is usually called a plasmonic solar cell. In contrast, the second application is that the photoelectric conversion does not have original photoabsorbing materials without metal NPs, but collects photocarriers directly from metal NPs using plasmon resonance absorption [64,65].

There are two main LSP effects that are expected to contribute to the improvement in PCE of plasmonic solar cells.

1. Enhancement effect by the electric field: improvement in photocurrent density J by enhancement of photocarrier generation rate using the strong electric field of optical near-field.

2. Effect of light scattering: improvement in J by anti-reflection and/or by increasing the optical length including light trapping.

Both the electric field enhancement effect and the light scattering effect can be revealed by the initial light absorption of metal NPs and by light scattering at metal NPs, respectively. As mentioned in the previous section, because the light absorption/scattering ratio depends on factors such as the size, shape, and configuration of NPs, the above two effects can be controlled mainly by these factors when NPs are inserted into solar cells.

There are also two collateral LSP effects that can be expected:

1. Improvement of the open circuit voltage V_OC

2. Improvement of the fill factor, FF.

By using the electric field enhancement and light scattering effects, the minimum thickness of the photoabsorbing layer to absorb sunlight completely can be reduced. The open circuit voltage (V_OC) in Si solar cells can be improved using these effects by reducing the total thickness of solar cells[58]. Furthermore, the FF related to the internal serial resistance of solar cells can also increase by reducing the total thickness of the solar cells using effects 1 and 2 for LSPs.

As mentioned above, light absorption that induces a strong electric field and/or the light scattering by LSPs at metal NPs can be observed by the irradiation of light with certain wavelengths. The ratio of absorption to scattering depends on the size of the metal NPs, generally increasing with decreasing size. This size-dependence of the absorption/scattering ratio has been reported for Ag NPs in Ref. [66]. In accordance with the theory, the absorption of Ag NPs whose diameter was smaller than about 50 nm was dominant over the extinction of light. The optical properties of NPs and the optical near-field near NPs are influenced by the following factors: 1) dielectric function of the NPs, 2) wavelength of incident light, 3) size and shape of metal NPs, 4) dielectric constant of the surrounding medium, and 5) configuration of metal NPs. Therefore, an appropriately chosen array of metal NPs can induce an interaction between LSPs and thus generate the so-called ‘hot spots’ – areas of greatly enhanced electromagnetic field compared with field strengths induced by single metal NP.

Rational design of plasmonic nanostructures therefore requires knowing the optical properties as a function of composition and morphology. The ability to compute optical properties is therefore critical for rational design. Ideally, one would like to be able to compute not only the optical spectrum for a given composition and geometry but also to design nanostructures that possess desired optical properties, i.e. to solve the inverse design problem. Modeling is rapidly gaining in importance for mechanistic understanding and rational design of plasmonic nanostructures [67–70], although atomistic electronic structure methods are still relatively less used. The difficulties in the modeling of plasmonic structures stem, in particular, from the length scale. Experimentally used NPs range from a few nm for hundreds of nm; in particular, NPs of sizes on the order of 1–10 nm are quite common. At this scale, the classic theory [1] already can fail while quantum mechanics based, ab initio approaches are quite expensive to use. As we will see below, time-dependent (TD) (Kohn-Sham) DFT is the most widely used ab initio approach to compute the plasmonic resonance peak. However, even ground state Kohn-Sham DFT [71,72] is very costly already for NPs of a few nm, which already imply models with hundreds of atoms. To model directly NPs of dozens of nm, models with many thousands of atoms are required which is impractical. Additional complications arise when computing spectra: linear response (Casida) TD-DFT [73,74] which is widely used in smaller-scale applications, is best suited for systems with a bandgap (non-metallic) and also becomes impractical even with a few nm due to the rapid explosion of the number of excitations that have to be included to compute the spectrum through the plasmonic peak [75,76]. For this reason, real-time TD-DFT [77] has found wider use in plasmonics but still suffers from the high cost of Kohn-Sham DFT. The scale also poses serious cost hurdles for the application of approaches including many body effects [78–81]. It is therefore important to consider novel avenues for atomistic modeling which are suited for the modeling of plasmonics.

In this review, we survey the computational modeling approaches used to model plasmon resonance features. We describe the classical approaches which are still widely used and can handle most particle sizes of interest to photovoltaics and, as we will show, based on comparisons with ab initio calculations, can remain valid even for relatively small NPs, but we place particular focus on approaches beyond classical electrodynamics such as ab initio approaches based on TD-DFT and related methods. While some more well-known approaches such as classical and TD-DFT were previously reviewed [82–84], the bulk of DFT-based literature deals with system sizes limited to a few nm, which is on the lower part of the range of particle sizes used in PV. We therefore specifically highlight emerging approaches which have the potential to significantly enhance modeling capabilities in the coming years, in particular, by allowing modeling at realistic length scales (of particles beyond 10 nm which are of most interest to plasmonic solar cells but remain problematic with Kohn-Sham DFT), such as time-dependent orbital-free DFT [85,86], linear response and real-time time-dependent density functional tight binding [87–90] and machine learning-based approaches, and by allowing modeling beyond the single-particle based (DFT) picture with many-body perturbation theory [78–81]. The latter approach, while being computationally expensive, is also finding increasing use in the modeling of plasmonics as access to computing power expands, and is therefore explicitly included in this review.

When using the terms ‘nanoparticles’, ‘nanoclusters’, ‘nanorods’, we usually follow the word choices of the original works. What should be called a nanoparticle or rather a cluster is somewhat debatable; especially as many DFT models use relatively small structures on the order of 1 nm, calling them clusters might be more appropriate. However, besides following the language of the reviewed works, we note that in many DFT works, systems with hundreds to thousands of atoms were used whose core preserves (which is not the case of molecular clusters) the structure (such as the crystal structure) of the respective bulk phase. This, and the fact that they are nanosized makes us adopt these terms.

Modeling of plasmonic properties

Classical approaches

When metal NPs and nanostructures are large enough, their electronic structure is metallic, i.e. a continuous band of electronic states across the Fermi energy is formed, the band structure, and it does not allow to distinguish between discretized states. Under this condition of bulk-like properties, their optical response can be modeled with classical electrodynamics methods, because oscillations of electron density behave classically at large length scales (which is of course in agreement with the assumption of the Drude model).

Classically, in the approximation of a spherical nanoparticle, the localized surface plasmonic response (LSPR) frequency of the main absorption peak in metals and sufficiently n-doped semiconductors NP (i.e. when they behave as metallic) with spherical shape can be estimated with the classical Drude model [91,92] as

(2.1.1) ${\omega }_{LSPR}={\left(\frac{{\omega }_p^2}{{ϵ}_{\mathrm{\infty }}+2ϵ}\right)}^{\frac{1}{2}}$(2.1.1)

where ${ϵ}_{\mathrm{\infty }}$ is the high-frequency dielectric constant and $ϵ$ the dielectric function of the environment. This relation can be derived from the plasmon resonance condition, the permittivity of bulk material including the response of inner electrons and neglecting the shift from the imaginary axis. The LSPR frequency assumes a static electric field approximation with particle sizes much smaller than the wavelengths (i.e., the dipole approximation). The plasma frequency ${\omega }_p$ of a volume plasmon [93] is given by

(2.1.2) ${\omega }_p=\sqrt{\frac{4\pi n_e}{m_e^{eff}}}$(2.1.2)

with $n_e$ the free electron density and $m_e^{eff}$ the electron’s effective mass (we use atomic units). There is no general solution for particles of arbitrary sizes using the original Drude model. Note that size-dependence can be included in a classical description of surface plasmons through consideration of restoring forces, i.e., a damping term depending on the size of the NP in the equation of motion of electrons, that are often neglected in other the formal derivation [94]. The latter correction is particularly important for small NPs.

Common numerical approaches to compute optical response of plasmonic materials with relatively large geometric features, generally more than a few 100 nm, include, for instance, the finite element method (FEM) and the finite difference time domain (FDTD) approach [95–97]. In both, the Maxwell equations of macroscopic electromagnetism are solved for a specific geometry and the obtainable dielectric constant within the Drude model, which provides an accurate description for free electron metals, determines the optical absorption for the near- and far-field. FDTD and FEM are only two examples of established numerical methods at this length scale. For a comprehensive review for classical methods and their applications, the reader is referred to Refs. [98,99]. Already with this rather simple approach, many plasmonic phenomena can be modeled, and it proved itself as a powerful technique capable of generating significant insight. For instance, Shen et al. [47] used FDTD to calculate the power conversion efficiency enhancement in Au-coated MAPbI₃ by plasmonic absorption enhancement [100]. Lee et at. also used FDTD to simulate the effect of nanoprisms and nanospheres (sizes of several 100 nm) of Au on the flux amplification of hot electrons in Schottky nanodiodes [101].

It was, however, early realized that the classical model and methods fail to describe intricate plasmonic features of smaller structures including NPs. The classic model fails to describe plasmons in geometries of small dimensions due to electron spill-out, nonlocal response and single-particle transitions [100,102–107]. The latter can be caused by discrete electronic states which are a result of finite sizes of NPs. Besides, the Drude theory is a single-pole model. Then, extensions of the theory were developed. For instance, McMahon demonstrated the importance of nonlocal Drude theory, which is a hydrodynamic model in real space represented through an extension of the dielectric function in the classical continuum limit allowing for conduction band excitations of electrons [108,109]. Therein, the corresponding optical response is defined with an additional spatially dispersive permittivity term. The nonlocality in the hydrodynamic Drude model (HDM) allows separating the longitudinal and transversal dielectric functions [110]. Because it is an extension of the local Drude model, HDM can also be implemented in existing FDTD codes [111,112]. Then, the nonlocal response of, e.g., systems with confined electrons gas in much smaller systems or those with more complex geometries can be simulated more accurately [109]. McMahon applied the method on nanowires with different diameters. Comparing classical calculations with local and nonlocal character for nanostructures with sizes up to 40 nm, they found significant differences in spectral shape for small NPs (5 nm) and blue shift with reduced intensity for larger ones. Fernandez-Dominguez et al. studied the plasmonic properties of two touching nanowires with a nonlocal classical theory employing the transformation-optics approach [113] and detected significant plasmonic field enhancements. Later, two separated nanowires forming a junction were studied also with a nonlocal classical theory [114]. Also for this case, plasmonic field enhancement was observed. It should be noted that charge transfer effects due to tunneling are not considered in the theory.

Further developments of classical methods aimed to reduce computational complexity and demands or try to include missing effects, as the above-mentioned tunneling effect, to increase the accuracy and hence predictive power of the method, besides revealing the nature of plasmonic effects. Luo et al. [104] suggested a local analogue model that maps the nonlocal plasmon properties in a metal to those of a composite material consisting of a thin dielectric shell and a metal core with the aim to obtain the same optical activity in the composite materials using a local theory as in the pure metal employing a nonlocal description. This allows describing the optical responses of angstrom-sized features with classical electrodynamics[104]. Because nonlocal features are rendered to a local model of the composite, rather simple calculations can be carried out using a standard FDTD solver for Maxwell equations. As mentioned above, tunneling effects between nanostructures with close vicinity (of 1 nm or less, a regime of interest for quantum plasmonics) are not considered in a classical approach. Esteban et al. [115] incorporated these effects via an fictitious effective medium that fills the vacuum of the junction and mimics the effect of tunneling, which motivated various studies to compare this and corrections/improvements to the classical Drude model [105,116,117].

Often, classical approaches using analytic models prevail in their computational efficiency. Within their validity, analytic models can obtain striking agreement to experimental results. Restrictions apply regarding the details of the nanoparticle such as shape, size, distributions in the host materials, etc. (as consequence of the made approximations). One frequently employed ansatz to compute the plasmonic properties, in particular the absorption coefficient, is the effective medium theory (EMT)[118]. In the EMT, the macroscopic properties of a composite material are modeled under certain assumptions. In the simplest form, the EMT sets the optical response of a composite material (that is a host medium with embedded nanoparticles) to be equal to the response of a homogeneous material with the same dielectric polarization in the visible region. This crude formulation, however, does not consider the particle size and details of the intrinsic optical response of the nanoparticle (states, scattering, etc.).

For dilute nanoparticle distributions, where the particles have a low fractional volume (or small particle sizes) and interactions between nanoparticles can be neglected, the Maxwell-Garnett (MG) EMT [119] describes the effective complex dielectric function of a composite with embedded spherical nanoparticles. The dipole approximation restricts particle sizes to be much smaller than the wavelength of the radiation, often only a few nanometers. However, details on the electronic structure, e.g. widening band gaps with decreasing particle size, are not considered. Further extensions of the MGEMT incorporate confinement effect through use of Mie theory [1]. Another extension of the EMT was introduced by Hinsen and Felderhof based on inclusion of higher-order multipoles in the representation of the effective dielectric function, Clausius-Mossotti relation [120]. Recently, Arefinia demonstrated the capabilities of classical approaches using these analytic models to reproduce experimental absorption coefficient of Au and Ag nanoparticles in hybrid composites, that are mixtures of organic and inorganic materials where the latter are the metal clusters [121]. However, the limitation of nanoparticle sizes remains, although remarkable agreement between experimental and simulations of metallic nanoparticles (Au, Ag with average diameters of about 20 to 40 nm with spherical shapes) were obtained.

Already the classical equations (2.1.1–2) imply possibility of control of the plasmon peak by controlling the electron concentration (which can be achieved by static doping or injection) as well as the effective mass (by choosing a material with a corresponding conduction band shape) but, being derived in the classical ansatz, ignore the effects of bandstructure, the inhomogeneity of the electron density and the nonlocality of the electron density response to the field as well charge transfer phenomena (including tunneling between NPs). They are therefore suited for relatively large nanoparticles with sizes larger than 10–100 nm, while for smaller nanoparticles there is need to explicitly account for discretization of energy levels and confinement effects. It has, however, been pointed out that for some types of systems the classical regime can hold well down to a couple of nm: Sinha-Roy et al. [122] concluded, by comparing classical and real-time TD-DFT calculations on Au and Ag rod between 19 and 145 atoms of different lengths and aspect ratios that the classical plasmonic frequencies are in surprisingly good agreement with those predicted by ab initio theory for the rods and quasi-one-dimensional chains, as long as collective surface-plasmon resonances lie far below the onset of d-electron excitations. Della Sala et al. [123] have modeled externally doped silicon nanocrystals classically as well as with TD-DFT and concluded that the classical approach can be effective for particles larger than 2 nm. Earlier, Pi and Delerue [124] observed convergence between tight-binding calculations and the Drude model from about 4 nm for P-doped Si. Other quantitative comparisons of the classical ansatz with quantum mechanical treatments, in particular, with time-dependent DFT (see below), have been made in the literature [105,125]. For example, Figure 4 shows a comparison of the absorption cross-section of a dimer of Au nanoparticles as a function of inter-particle separation obtained with the classical model (lines) and real-time TD-DFT (using the jellium approximation) [125].

Figure 4. Comparison of the absorption cross-section of a dimer of jellium Au nanoparticles as a function of inter-particle separation obtained with the classical model (lines) and real-time TD-DFT (color map). CTP and CTP1 are the charge transfer plasmon modes of the coupled dimer. Bonding dimer plasmon (BDP) originating from the hybridization of the dipolar plasmon modes of the individual nanospheres and the higher frequency hybridized bonding quadrupolar plasmon (BQP) modes are also indicated. Adapted with permission from Ref. [125]. Copyright 2013 American Chemical Society

Corrections to the Drude model have been proposed [126–130] to effectively capture the resonance frequency shift due to nanosizing and quantization effects. Corrected Drude models were used to compute plasmonic resonances in Al-doped ZnO [131], B-doped Si [132], and naturally n-doped HgS [100].

Such corrections effectively provide a better estimation of the peak position but they do not overcome in substance the intrinsic limitations of classical models. For that, one needs to solve a multidimensional Schrödinger equation. This is computationally costly but can be done with ab initio approaches. Alternatively, one can solve a simplified quantum problem. Liu et al. [133] modeled a NP as a sphere with a uniform charge density and solved the Schrödinger equation for a particle in a sphere using the DVR (discrete variable representation) [134] approach. They applied the model to a near-spherical ZnO NP and obtained good agreement with DFT for distribution of excess electrons. The spherical models ignore the corrugation of a real NP which can have as important an effect as other quantum corrections [135]. In particular, accounting for the atomistic structure of the NP has been shown to lead to an increased peak width (Figure 5). Quantum approaches involving the solution of the Schrödinger equation may still ignore important atomistic details such as the ionic potential, which may be replaced by a uniform (jellium) background or by an effective potential. Effective potentials allow nevertheless capturing some important phenomena such as density accumulation near the surface [133].

Figure 5. Illustration of the effects of inclusion of quantum effects and atomistic structure of the nanoparticle on the example of a spherical (BCC) Na₃₃₁ and icosahedral (ICO) Na₂₉₇ nanoparticles. The solid lines are absorption spectra including the main LSP peak computed with real-time TD-DFT. Dashed blue line: jellium model. The classical local-optics results are shown as dashed-dotted green line. Adapted with permission from Ref. [135]. Copyright 2014 American Physical Society

Linear response TD-DFT and TD-DFTB based on Casida formalism and related approaches

Computing the plasmonic resonant peak in a truly quantum approximation, however, requires computation of electron excitations. This can be done with time-dependent DFT [136,137]. Linear response (LR) time-dependent TD-DFT based on Casida formalism [73,74] is widely used for computing optical properties of various types of materials and is implemented in most popular molecular modeling codes such as Gaussian [138] and has recently been implemented in major periodic codes [139,140]. In LR TD-DFT, the excitation spectrum ω is obtained from the eigenvalue problem

(2.2.1) $\left[\begin{array}{cc}A& B\\ B^{\ast }& A^{\ast }\end{array}\right]\left[\begin{array}{c}X\\ Y\end{array}\right]=\omega \left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\left[\begin{array}{c}X\\ Y\end{array}\right]$(2.2.1)

where matrices A and B depend on the integrals over Kohn-Sham orbitals [74]

(2.2.2) $K_{ia\mu ,jb\nu }=\int \int {\varphi }_{i\mu }^{\ast }\left(\mathit{r}\right){\varphi }_{a\mu }\left(\mathit{r}\right)\left(\frac{1}{\left|\mathit{r}-{\mathit{r}}^{\prime }\right|}+\frac{{\delta }^2{E}_{XC}}{\delta {\rho }_{\mu }\left(\mathit{r}\right)\delta {\rho }_{\nu }\left({\mathit{r}}^{\mathbf{\prime }}\right)}\right){\varphi }_{j\nu }^{\ast }\left({\mathit{r}}^{\prime }\right){\varphi }_{b\nu }^{\hspace{0.17em}}\left({\mathit{r}}^{\prime }\right)d\mathit{r}d{\mathit{r}}^{\prime }$(2.2.2)

where indices i, j and a, b label occupied and virtual orbitals $\varphi$, respectively, and indices μ and $v$ denote spin, ρ is the density, and E_XC the exchange-correlation energy. This approach is simpler and often less demanding computationally and in terms of human effort than real-time (RT) propagation (see below) when the number of involved orbitals is not high. In realistically-sized nanoparticles, however, it is very high, and the cost grows rapidly, notably, due to the necessity of including very large numbers of unoccupied as well as occupied orbitals in the calculation to achieve excitation energies up to and beyond the plasmonic resonance peak [75]. LR TD-DFT can be implemented using a density-response function in a way well suited to periodic systems [141,142]; Yan et al. [143] used such a TD-DFT approach to investigate surface plasmons of Ag(111) and H/Ag(111).

The Casida formulation is best suited for systems with a bandgap and a relatively low density of states but has also been used to model plasmonic spectra. In Ref. [144], LR TD-DFT calculations were performed on Ag tetrahedral clusters with up to 120 atoms. In Ref. [145] LR TD-DFT was used to compute spectra including plasmonic peaks of small (up to 165 atoms) silver, gold, and bimetallic Au-Ag particles of different shapes and sizes. The accuracy can be judged from the comparison with the experimental spectrum for Ag₅₅ shown in Figure 6 and is more qualitative than quantitative. One should consider however the compounded influences of the Casida approximation, the approximation to the exchange-correlation functional (LDA for optimization and LB94 [146] for properties), the relatively small (double-ζ) basis used, and the scalar relativistic ZORA (zeroth order regular approximation) approximation. LSPR of gold nanoclusters (up to 314 Au atoms) covered by a thiolate monolayer were studied by LR TD-DFT by Malola et al. [147], which is an example of combining metal NP with molecular surface modifications.

Figure 6. Absorption spectrum of a Ag₅₅ silver NP with TD-DFT (black line), TD-DFT+TB (red line) with DZ basis set and TD-DFTB (blue line) with DFTB.org/Hyb-0-2 for Ag and auorg-1-1 for Au parameter. The isosurface plot is for the TD-DFT+TB transition density of the most intense peak. Adapted with permission from Ref. [145] under CC-BY-NC-ND license

The issue of the cost associated with the need to include a very large number of excited states can be palliated with TD-DFTB (time-dependent density functional tight binding) [89,90], with which thousands of excited states can be included routinely and CPU time requirements can be lowered by three orders of magnitude vs. DFT [144]. The self-consistent-charge density-functional-tight-binding (SCC-DFTB) approach achieves computational efficiency by using an expansion of the exact DFT energy functional around a reference density. The reader is referred to Refs. [148–150] for a description of DFTB. The method needs to be parameterized with the so-called Slater-Koster (SK) parameters and allows for the formulation of TD-DFTB based on the Casida formalism [89,90]. Asadi-Aghbolaghi et al. [145] compared (LR) TD-DFTB with (LR) TD-DFT calculations with their proposed method which combines a full DFT ground state calculation with tight-binding approximations in the linear response calculation (called TD-DFT+TB), when computing optical properties of Ag, Au and bimetallic Ag−Au nanoparticles with tetrahedral icosahedral structures with up to 165 atoms. The number of included excited states ranged 3,000 − 30,000 depending on the cluster. An example of a comparison between the three methods and the experiment for Ag₅₅ is shown in Figure 6 where spectra including plasmonic peaks are shown. DFTB [149,151] is usually parameterized to DFT calculations using GGA (generalized gradient approximation) functionals, and as a result tends to underestimate bandgaps and excitation energies (as can be appreciated from Figure 6), in addition to errors due to the DFTB approximation itself such as accuracy of parameterization of integrals and limitations due to the use of pair-wise repulsive potentials. It is therefore important to benchmark TD-DFTB calculations to TD-DFT or experiments. The TD-DFT+TB approach provided an accuracy close to that of TD-DFT while lowering the computational cost. Especially as DFTB parameterizations suitable for metallic nanoparticles and customized for plasmonics become available, one can expect an increase in TD-DFTB applications in plasmonics. D’Agostino et al. [144] showed, also computing Ag nanoparticles (4,000 included excitations), that with appropriate (tuned) parameterization, (LR) TD-DFTB can be about as accurate as (LR) TD-DFT when computing spectra of plasmonic systems.

Liu et al. [152] studied with TD-DFTB the absorption spectra for nanoparticle (up to Ag₁₆₄) silver dimers at various inter-particle distances. 4,000 excited states were included. Both homodimers and heterodimers were considered. They found that as the inter-particle distance decreases, a red shift arises from contributions of the transition dipole moment that are aligned along the z (inter-dimer) axis; blue shifts occur for peaks that originate from transition dipole moment components in the x and y directions. When the nanoparticles were similar in size, the features in the absorption spectra were more sensitive to the inter-particle distance. Alkan and Aikens studied with TD-DFTB silver nanorods and nanorod dimers with up to 2000 atoms in side-by-side and end-to-end configurations[153]. To reach excitation energies up to 4 eV, up to 18,000 excited states had to be included. TD-DFTB was compared to TD-DFT for small clusters and dimers and it was found, unsurprisingly, that TD-DFTB excitation energies are underestimated. Nevertheless, similar spectral shapes and features were obtained with the two methods. Douglas-Gallardo et al. [154] studied with TD-DFTB icosahedral aluminum nanoclusters (up to Al₉₂₃) with different degrees of oxidation (by oxygen), and aluminum nanoclusters with a disordered structure to study how the loss of crystallinity affects the optical properties. It was found that oxygen induces a redshift and a broadening of the dipole surface plasmon resonance band, and that loss of crystallinity leads to broadening of spectral features. In Ref. [155], the interactions, in an ultra-near-field regime, between a localized surface plasmon excitable in a tetrahedral Ag₂₀ cluster and a molecular exciton with excitation energy in the same range were studied with TD-DFT for smaller molecules and TD-DFTB for larger. For metal–molecule distances below 5 Å, the optical response of the system resulted in the appearance of a double peak structure.

Another approach developed to alleviate the issue of the computational cost with a large number of excited states computes the spectrum from the calculation of the complex dynamical polarizability independently at each frequency; this ensures good parallelizability. Baseggio et al. introduced the method for the calculation of plasmonic properties [156,157] and applied it to model Ag and Au clusters with up to 309 atoms. They showed near-quantitative agreement with Casida TD-DFT; e.g. for Au clusters deviations within 0.2 eV of the plasmon resonance peak position were obtained with the two methods [156]. Recently, calculations of spectra from real frequency-dependent polarizability computed separately at each frequency have also been proposed [158,159] but are yet to be explored for plasmonics.

Real-time TD-DFT

The most widely used TD-DFT approach in modeling of plasmonic spectra appears to be time-propagation in real space [77]. In it, the absorption spectrum is computed from the imaginary part of the complex frequency dependent dynamic polarizability $\alpha \left(\omega \right)$ as

(2.3.1) ${\sigma }_{abs}\left(\omega \right)=\frac{4\pi \omega }{c}\mathrm{I}\mathrm{m}\left[\alpha \left(\omega \right)\right]$(2.3.1)

where ${\sigma }_{abs}\left(\omega \right)$ is the absorption cross-section. The dynamic polarizability is computed from perturbation-induced time-dependent electron density variation $\delta n\left(\mathit{r},t\right)$ which obtains following a delta-like perturbation with an electric field $\mathbf{E}\left(\mathit{r},t\right)=k_0\delta \left(t\right){\mathit{e}}_{\mathit{x}}$, where ${\mathit{e}}_{\mathit{x}}$ is a unit vector of the electric field direction (here assumed to be along the x axis) and $k_0$ is the amplitude, which is typically small enough (on the order of ${10}^{-3}$ a.u.) to ensure a linear response regime. Harmonic perturbations can also be used. The density variation computed over the propagation time T is transformed into its spectral representation as

(2.3.2) $\delta n\left(\mathit{r},\omega \right)=\underset{0}{\overset{T}{\int }}\delta n\left(\mathit{r},t\right)e^{\left(i\omega -\gamma \right)t}dt$(2.3.2)

where the damping frequency γ is introduced to account for the spectral broadening due to various losses. It is easy to see that $\delta n\left(\mathit{r},\omega \right)$ is a complex electron density induced by the field ${\mathbf{E}}_{\omega }\left(\mathit{r},t\right)=k_0{e}^{-i\omega t}{\mathit{e}}_{\mathit{x}}$ (i.e. by the corresponding Fourier component of the delta-like perturbation). The dynamic polarizability is then computed as

(2.3.3) $\alpha \left(\omega \right)=-\frac{1}{k_0}\int x\delta n\left(\mathit{r},\omega \right)d\mathit{r}$(2.3.3)

In applications, one most often uses the adiabatic (i.e. the exchange-correlation potential depending only on the instantaneous density – the approximation widely used with TD-DFT for other applications as well) local density approximation (ALDA) with the exchange-correlation functional of Gunnarson and Lundqvist [160]. ALDA provides in many cases a rather accurate description with relative errors within 10% (which is small compared to the well-known band-gap issue of DFT with relative errors of up to 50%) [142,161], while representing a rather simple formalism. Other exchange-correlation approximations (e.g. AGGA) are also used [84,122,157,162]. It should be mentioned the failure of TD-ALDA for periodic systems has motivated developments of many-body techniques for these systems such as GW and BSE formalisms (see below), while other works aim at increasing TD-DFT accuracy beyond that of typical LDA and GGA functionals with, for instance, the adiabatic Gritsenko-van Leeuwen–van Lenthe–Baerends – solid-correlation (GLLB-SC) exchange-correlation functional which reduced self-interaction error by the introduction of an orbital energy-dependent localization of the exchange hole and has the advantage of improved d-states description [162–164].

The time propagation of the density $n\left(\mathit{r},t\right)=\sum _i^{\hspace{0.17em}}{f}_i{\left|{\mathrm{\Psi }}_i\left(\mathit{r},t\right)\right|}^2$ (where f_i are occupation numbers of the KS orbitals ${\mathrm{\Psi }}_{\mathrm{i}}\left(\mathrm{r},\mathrm{t}\right)$) can be done via time propagation of the basis coefficients of orbitals:

${\mathrm{\Psi }}_i\left(\mathit{r},t\right)=\sum _n^{\hspace{0.17em}}{C}_{ni}\left(t\right){\varphi }_n\left(\mathit{r}\right)$

(2.3.4) $i\mathit{S}\frac{\mathrm{\partial }\mathit{C}\left(t\right)}{\mathrm{\partial }t}=\mathit{H}\left(t\right)\mathit{C}\left(t\right)$(2.3.4)

where C is the matrix of expansion coefficients, H is the Hamiltonian matrix $H_{nm}=⟨{\varphi }_n|{\hat{H}}_{KS}|{\varphi }_m⟩$, and S the overlap matrix $S_{nm}=⟨{\varphi }_n|{\varphi }_m⟩$.

Zhang et al. employed and compared RT TD-DFT and FDTD calculations on Al NPs with tetragonal shape and their size dependence on optical properties [165]. For clusters containing up to 560 atoms, corresponding to cluster sizes of less than 3.7 nm, FDTD and RT TD-DFT obtained noticibly different spectra (besides, RT TD-DFT and LR TD-DFT obtained agreement in the low energy excitations). The detected discrepancy can, of course, be explained with the missing reproduction of confinement effects by the FDTD method. At very large particle sizes, agreement can be expected to be much better although, as mentioned above, KS-based ab initio methods are just not yet capable of handleing the required system sizes. Instead, the developments of RT TD-OFDFT and RT TD-DFTB are more promising (see below). Kuisma et al. [162] studied with RT TD-DFT silver clusters of size up to Ag₅₆₁. They noted that already at about 2 nm LSPR peaks agreed with the classical limit (of 3.4 eV), while higher peak energies were obtained with smaller NPs. Lopez-Lozano et al [166] used RT TD-DFT to compute surface plasmon resonances of gold nanorods of length up to 7 nm. They studied the size and aspect ratio dependence and showed that by changing the length and width, the energy and the character of the resonance can be tuned. Li et al. [167] modeled Na nanoclusters up to Na₃₃₁; in particular, they used different cluster shapes including some rather jagged boundaries and obtained good agreement with experiment (peak position accurate on the order of 0.1 eV) for cases where comparable experimental data were found. Iida et al. [168] computed LSPR for Au clusters with up to 1414 atoms, which is one of the largest scale calculations at the full RT TD-DFT.

In plasmonic systems, when nanoparticles are in close proximity to each other, there is inter-particle interaction due in particular to electron density spillover beyond the NP surface. Charge transfer excitations and electron tunneling can occur. LSP can hybridize in such systems. These are practically important phenomena, in particular, for the realization of nanodevices with sub-nm gaps [169,170]. A critical phenomenon is plasmonic field enhancement, including in the region between nanoparticles. The enhancement is strongly affected by LSP hybridization [125,135,171]. A dimer of two interacting nanoparticle serves as a prototypical system to study such interactions. Marinica et al. [125] computed optical properties and field enhancement as a function of inter-particle distance of a dimer of spherical nanoparticles, with real-time TD-DFT and the jellium model (of gold) (Figure 7). They showed, by comparing real time TD-DFT results with both classical as well as a linear quantum mechanical description of the system, that the latter approximations fail to describe field enhancement even for moderate incident light intensities. The charge transfer current between the NPs is also underestimated in the linear approximation. Teperik at al. [105] modeled with real time TD-DFT coupled Na nanowire dimers of radii up to 10 nm also using the jellium approximation. They computed optical spectra, induced charge densities, and near fields, revealing nonlocal quantum effects (consisting of electron tunneling between the nanowires at small junction widths and dynamical screening) determining the plasmonic modes and field enhancement. They compared TD-DFT results to the Drude model (see section 2.1).

Figure 7. Plasmonic field enhancement at the point between the nanoparticles as a function of separation between jellium Au nanoparticles computed in Ref. [125] with RT TD-DFT, for different light intensities. Adapted with permission from Ref. [125]. Copyright 2013 American Chemical Society

In Figure 8, the optical responses of Na nanowire dimers computed with different methods in Ref. [172] are shown. In particular, they compare accurate TD-DFT calculations to classical models with and without additional corrections such as nonlocal and tunneling terms (as in the nonlocal HDM and the quantum corrected model [115], respectively). Most obvious differences in the plasmonic response to the TD-DFT can be seen in the Drude and NLHD calculations were intensity ratios and relative peak positions are mismatched. For the Drude model, even the number of modes assigned to charge transfer (for separation less than 0 which was also considered, i.e. overlap of nanowires) is overestimated. The comparison between QCM and TD-DFT shows good agreement besides a mismatch in the intensity ratio between dipole and quadrupole excitations. Also, it should be noted that qualitative differences of the changes of the type of plasmon excitation over different separation distances arise. For instance, the QP peak in TD-DFT continuously changes to a charge transfer plasmon (C2) with decreasing separation while in NLHD, this is somewhat abrupt.

Figure 8. Comparison of length-normalized cross sections computed with different methods for a dime of Na nanowires with diameters of 9.8 nm at different separation distances S. From left to right (jellium) TD-DFT, a quantum corrected model (QCM), Drude model and the nonlocal hydrodynamic model (NLHD). Labeled plasmon poles include dipole plasmon (DP), quadrupole plasmon (QP), charge transfer plasmon (C1 and C2) and higher-order mode (HM). Adapted with permission from Ref. [172]. Copyright 2013 American Physical Society

With TD-DFT, the enhancement of the local field strength between the nanowire dimers is quenched compared to classical models. Moreover, the characteristics of the field enhancement significantly differ. In the classical calculation, the point of largest enhancement is on the symmetry axis between the wires whereas in the quantum theory calculations it is slighty shifted. As pointed out by Teperik, the classical models do not consider the effect of quantum tunneling at small separation distances (on the subnanometer scale) between the wires causing these disagreements. In the quantum corrected model developed by Esteban et al., tunneling effects were effectively included by a fictious medium and reproducing the quenching [115]. In addition, accumulation of opposite effective charges across the junction in the wires should be considered. Considering both of these effects improves the physical description going beyond earlier reports of plasmonic properties with (quasi-)local and nonlocal classical methods applied to similar systems (i.e. couple nanowires) as reported by Luo et al. [104] and Fernandez-Dominguez et al. [114] At the same time, the comparison to these results employing classical models helps reveal the nature of plasmonic effects. It was concluded that the classical description overestimates the number of resonances and the local field enhancement due to the coupling effect across the vacuum junction. Furthermore, there is no significant difference except shifts in frequencies among the local and nonlocal model. Obviously the TD-DFT results have a more accurate description of the electronic effects than the classical approaches including electron spill-out, dynamic screening and transitions. For instance, Andersen et al. found agreement between classical and TD-DFT results in a Na nanowire dimer with nanometer-scale separations, while at subnanometer distances tunneling yields a blue-shift in the plasmonic peaks caused by a decrease of induced charges [173]. In general, however, the nature of the screening effects depend on the involved electrons. In transition metal elements, the local character of d-electrons and their interband transitions can result in variations of plasmonic response and often manifest in broader features at higher excitation energies [142,143,162,166,174]. Also, they tend to opposite fields leading to dynamical screening, whereas the s- and p-electrons, e.g. in Na, behave as quasi-free electrons [122,175]. Because the treatment of d-electrons in DFT is non-trivial and several schemes such as DFT+U and long-range corrected functional were proposed to accuractly describe these states, this issue remains challenging and of importance [122,162].

Apart from all benefits that TD-DFT offers and the success that classical models achieved, many surface plasmon-related phenomena are rather complex requiring explicit atomistic modeling of realistic sizes and shapes of particles without fallback on simplifications such as the jellium model, while many of the even small nanoparticle sizes are still outside the scope of TD-DFT. However, the continuum electrodynamics-based approaches remain their validity only for the large NPs and structures. Between these lengths scales, a wide theoretical gap opens that either of these theoretical frameworks is incapable or simply inadequate to describe plasmonic properties (although in some cases classical approaches continue to work well at surprisingly small scale, as mentioned above). The recent developments in RT TD-DFTB and RT TD-OFDFT might be able to fill this gap briding the atomic scale with the mesoscale.

Real-time TD-DFTB

The use of DFTB promises significant computational cost savings and therefore ability to treat larger systems also with the real-time time-dependent approaches. RT TD-DFTB has been gaining use in plasmonics [87,88,176,177]. In RT TD-DFTB, the Liouville-von Neumann equation

(2.3.1.4) $\frac{\mathrm{\partial }\hat{\rho }}{\mathrm{\partial }t}=-i\left(S^{-1}\hat{H}\left[\hat{\rho }\right]\hat{\rho }-\hat{\rho }\hat{H}\left[\hat{\rho }\right]S\right)$(2.3.1.4)

where S is the overlap matrix (which is part of the DFTB parameterization) is numerically integrated to propagate the reduced one-electron density matrix $\hat{\rho }$ in the presence of a perturbation, typically delta-like or harmonic. Here $\hat{H}={\hat{H}}_{GS}+V_{ext}\left(t\right)$, where ${\hat{H}}_{GS}$ is the ground state Hamiltonian and $V_{ext}\left(t\right)$ the external potential due to the perturbing field. Optical absorption spectra are then computed similar to Eq. (2.3.1).

RT TD-DFTB was used by Bonafé et al. [87] to compute optical response of isocahedral silver nanoparticles with up to 309 atoms. They saw evidence of breathing oscillations of the NP in the sub-picosecond regime via surface plasmon excitations whereby excited electrons populate excited states of antibonding character, weakening the bonds and causing an expansion of the nanoparticle, which then vibrates around a new equilibrium radius. Lu et al. [176] used RT TD-DFTB to study optical properties of gold nanoparticles with 3-mercapto-1,2-propanediol ligands on the surface which exhibit plasmonic circular dischroism (PCD). Absorption spectra were computed with the method for different directions of the external field and the average orientation of the transition dipole moment was used to help understand the PCD behavior.

Ilawe and Wong [177] showed that RT TD-DFTB produces a similar LSPR peak of Na₅₅ as RT TD-DFT [167] and computed the response of a plasmonic nanoantenna made of four Na₅₅ NPs. Douglas-Gallardo et al. [88] computed optical properties of different icosahedral silver and gold nanoclusters of diameters between 1 and 4 nanometers (up to 1415 atoms, see Figure 9) and studied the plasmon-induced hot-carrier generation process under laser illumination. RT TD-DFTB is getting recognized as a promising modeling method for many types of plasmonic structures [177] and is expected to gain in importance in the coming years.

Figure 9. RT TD-DFTB optical absorption spectra of gold (left) and silver (right) NPs of different sized. Adapted from from ChemRxiv. 2019; doi: 10.26434/chemrxiv.7905344.v1under CC BY NC ND License

Real-time TD-OFDFT

Practically used plasmonic nanoparticles often have sizes (10¹–10³ nm) which are not feasible for routine DFT calculations. With LR TD-DFT, the number of included excitations (related to the number of considered pairs of occupied and virtual orbitals in Eq. (2.2.2)) needed to reach the plasmonic peak rapidly balloons beyond feasibility, in particular, due to a high density of states in larger systems. Such calculations are still infeasible even with TD-DFTB where thousands of excited states can routinely be included [76]. However, more than tens of thousands excitations would have to be included for NPs of more than a few nm. Real-time TD-DFT circumvents this issue, but the underlying Kohn-Sham (KS) DFT calculations are still very costly even for particle sizes on the order of 1–10 nm and are prohibitively costly for 100–1000 nm. This ultimately has to do with the near-cubic headline scaling of KS-DFT. Linear-scaling versions of KS-DFT [178–180] could be useful but come at a cost of large prefactors and additional approximations; specifically, their basic ansatz of forming a sparse Hamiltonian matrix (typically by the use of finite-support basis functions) is not well suited for metallic systems.

Orbital-free (OF) DFT [86,181] offers a genuinely linear-scaling approach to DFT with small prefactors allowing routine calculations of million-atom systems. Even though not all kinetic energy functionals (KEF) are linear-scaling, semi-local KEFs [182,183] are, and non-local KEFs [184–188] in general results in a quadratic to n log(n) scaling which is still much advantageous over KS-DFT. The use of OFDFT in applications has until recently been held back for many years by the lack of sufficiently accurate KEFs, which aim to replace the non-interacting kinetic energy (KE) E_kin or kinetic energy density (KED) τ(r) of KS-DFT,

(2.3.2.1) $E_{kin}^{KS}=\int d\mathit{r}{\tau }_{KS}\left(\mathit{r}\right)$(2.3.2.1)

where

(2.3.2.2) ${\tau }_{\mathrm{K}\mathrm{S}}\left(\mathit{r}\right)=-\frac{1}{2}\sum _i{\psi }_i\left(\mathit{r}\right)\mathrm{\Delta }{\psi }_i\left(\mathit{r}\right)$(2.3.2.2)

or its positively definite version

(2.3.2.3) ${\tau }_{KS}^{+}\left(\mathit{r}\right)=\frac{1}{2}\sum _i{\left|\mathrm{\nabla }{\psi }_i\left(\mathit{r}\right)\right|}^2={\tau }_{\mathrm{K}\mathrm{S}}\left(\mathit{r}\right)+\frac{1}{4}\mathrm{\Delta }\rho \left(\mathit{r}\right)$(2.3.2.3)

with functionals of density (or any set of density-dependent variables):

(2.3.2.4) $\tau \left(\mathit{r}\right)=\tau \left[\rho \left(\mathit{r}\right)\right]\approx {\tau }_{KS}\left(\mathit{r}\right)E_{kin}=E_{kin}\left[\rho \left(\mathit{r}\right)\right]\approx E_{kin}^{KS}$(2.3.2.4)

However, for light metals such as Mg, Al, or Na considered in a previous section, existing KEFs are already quantitatively accurate [189–193]. Moreover, in recent years, significant advances in KEF development, including advances in semi-local (i.e. linear-scaling) KEFs [182,183] as well as machine-learned KEFs [194–200], are being made. The newest functionals appearing in recent years can treat systems with more non-uniform densities and give hope that in relatively near future many materials used in plasmonics (metals as well as non-metals) can be modeled with OFDFT with at least semi-quantitative accuracy.

The absence of the orbitals permits much better scaling but also makes impossible the use of LR TD-DFT (as well as DOS analysis and other analyses typically done with KS DFT which rely on KS states). However, real-time TD-DFT is possible with OFDFT. TD-OFDFT can be rationalized by considering the OFDFT formulation using a collective orbital $\psi \left(\mathit{r}\right)=\sqrt{\rho \left(\mathit{r}\right)}$. In this case one solves the equation

(2.3.2.5) $\hat{H}\psi \left(\mathit{r}\right)=\mu \psi \left(\mathit{r}\right)$(2.3.2.5)

where μ is the chemical potential, with

(2.3.2.6) $\hat{H}=-\frac{1}{2}\mathrm{\Delta }+\frac{\delta E^{Pauli}}{\delta \rho \left(\mathit{r}\right)}+v_s\left(\mathit{r}\right)$(2.3.2.6)

Here $E^{Pauli}=E_{kin}-E_{vW}$ with $E_{vW}=\frac{1}{8}\int \frac{{\left|\mathrm{\nabla }\rho \left(\mathit{r}\right)\right|}^2}{\rho \left(\mathit{r}\right)}d\mathit{r}$ being the von Weizsäcker kinetic energy. $v_s\left(\mathit{r}\right)=\frac{\delta E_H\left[\rho \left(\mathit{r}\right)\right]}{\delta \rho \left(\mathit{r}\right)}+\frac{\delta E_{XC}\left[\rho \left(\mathit{r}\right)\right]}{\delta \rho \left(\mathit{r}\right)}+V_{ext}\left(\mathit{r}\right)$, where E_H is the Hartree, E_XC the exchange-correlation, and V_ext the external potential. The TD-OFDFT is then introduced by considering a time-dependent collective orbital defined based on the time-dependent density $\psi \left(\mathit{r},t\right)=\sqrt{\rho \left(\mathit{r},t\right)}{e}^{iS\left(\mathit{r},t\right)}$. One then solves the equation

(2.3.2.7) $i\frac{\mathrm{\partial }\psi \left(\mathit{r},t\right)}{\mathrm{\partial }t}=\hat{H}\psi \left(\mathit{r},t\right)$(2.3.2.7)

with

(2.3.2.8) $\hat{H}=-\frac{1}{2}\mathrm{\Delta }+\frac{\delta E^{Pauli}}{\delta \rho \left(\mathit{r},t\right)}+v_s\left(\mathit{r},t\right)$(2.3.2.8)

These equations can be numerically solved with techniques similar to those used in RT KS-DFT. RT TD-OFDFT has recently been implemented in the DFTpy code [85]. In Ref. [85], optical spectra of Mg clusters up to Mg₅₀ have been computed, using a simple Thomas-Fermi-von Weizsäcker KEF and an LDA exchange-correlation functional, with RT TD-OFDFT and were compared with KS-DFT results with decent agreement.

Given its computational cost advantage, RT TD-OFDFT has the potential to be a game changer and is expected to gain significant ground in the modeling of plasmonics in the coming years, in particular, by modeling systems of realistic sizes. The development of KEFs suited for optical calculations will remain an active area of improvement of RT TD-OFDFT, in particular for the correct description of localized d-states. Another challenging area is the development of local pseudopotentials (LPP) which are required by OFDFT for the atoms constituting most widely used plasmonic materials such as Ag or Au. A good LPP for Ag has been proposed [201]; the Ag case is simpler as its valence properties are largely determined by a single s electron. For other atoms beyond light metals new ideas are needed. Machine learning holds promise in this area [202].

Beyond the single-particle Hamiltonian: including many-body and excitonic effects

DFT remains a theory exact by construction only for the description of the ground state (GS) properties and can fail when applied to excited states[161]. This means it is not suitable in its native form for the prediction of the electronic bandgap, the key property in optoelectronics, to compare with experimentally detected gaps, e.g., via (direct/indirect) photoemission spectroscopy (PES) results[161], and other optical properties including excitonic effects. The neglected effects of self-interaction and the well-known derivative discontinuity [203,204] represent the two main shortcomings that hinder the use of DFT for high-accuracy excited state calculations.

To solve these single-particle derived issues, the usage of hybrid functionals has been suggested or, alternatively, topic of the present section, the replacement of the exchange-correlation potential contribution with the so-called self-energy (Σ) either in a ‘one-shot’ approach or in an iterative, self-consistent way, i.e. updating the Green function G or both the Green function and the screening potential, W. The time-ordered Green’s function (governed by the T, the time-ordering operator) assumes the general form

(2.4.1) $G\left(\mathit{x},t;{\mathit{x}}^{\mathbf{\prime }},t{ }^{\mathrm{\prime }}\right)=\frac{⟨{\mathrm{\Psi }}_0|T\left(\hat{\psi }\left(\mathit{x},\mathrm{t}\right){\hat{\psi }}^{†}\left({\mathit{x}}^{\prime },{\mathrm{t}}^{\prime }\right)\right)|{\mathrm{\Psi }}_0⟩}{⟨{\mathrm{\Psi }}_0|{\mathrm{\Psi }}_0⟩}$(2.4.1)

where $\hat{\psi }\left(\mathit{x},t\right)$ and ${\hat{\psi }}^{†}\left({\mathit{x}}^{\mathbf{\prime }},t^{\prime }\right)$ are the Heisenberg creation and annihilation field operators, respectively, and ${\psi }_0$ the N-particle ground state. After a Fourier transform, in the frequency domain the Eq. (2.4.1) assumes the so-called Lehmann (spectral) representation

(2.4.2) $G\left(x,x\text{'};\omega \right)={\displaystyle \sum _k\frac{g_k\left(x\right)g_k^{*}(x\text{'})}{\omega -\left(E_k^{N+1}-E_0^N\right)+i\eta }}+{\displaystyle \sum _k\frac{f_k\left(x\right)f_k^{*}\left(x\prime \right)}{\omega -\left(E_k^{N-1}-E_0^N\right)-i\eta }}$(2.4.2)

where $f_k\left(\mathit{x}\right)=⟨{\mathrm{\Psi }}_k^{N-1}|\psi \left(\mathit{x}\right){|\hat{\mathrm{\Psi }}}_0^N⟩$ and $g_k\left(\mathit{x}\right)=⟨{\mathrm{\Psi }}_0^N|\hat{\psi }\left(\mathit{x}\right){|\mathrm{\Psi }}_k^{N+1}⟩$. In eq. (2.4.2) $E_k^{N\pm 1}$ is the energy associated with the k^th excited states of the N ± 1-particle systems, $|{\mathrm{\Psi }}_k^{N\pm 1}\text{rhigt}⟩$. The relevance of Eq. (2.4.2) is that its poles represent the energy difference that can be compared, with experimentally detected energy difference obtained by direct and inverse PES, i.e. the energies of the system upon addition and removal of a particle.

The imaginary part of Eq. (2.4.2), the spectral representation, takes the form (similar relation holds for the advanced Green function, G^A)

(2.4.3) $A\left(\mathit{x},{\mathit{x}}^{{ }^{\mathbf{\prime }}};\omega \right)=-\frac{1}{\pi }Im{G}^R\left(\mathit{x},{\mathit{x}}^{{ }^{\mathbf{\prime }}};\omega \right)=\sum _k{f}_k\left(\mathit{x}\right)f_k^{\ast }\left({\mathit{x}}^{\prime }\right)\delta \left(\omega +E_k^{N-1}-E_0^N\right)+\sum _k{g}_k\left(\mathit{x}\right)g_k^{\ast }\left({\mathit{x}}^{\prime }\right)\delta \left(\omega +E_0^N-E_k^{N+1}\right)$(2.4.3)

and assumes a physical meaning of localized density of states. Here G^A is the advanced and G^R the retarded Green’s function, corresponding to different signs (direction of time) in the time propagator. The concept of quasiparticle (QP) is a direct result of Eq. (2.4.3). In terms of non-interacting electrons, the spectral function consists of several delta function peaks, each of them corresponding to the excitation of a single non-interacting particle. Turning on the interaction, the matrix elements of the spectral function will consist of several elements with finite amplitudes: the resulting peak associated with this merging process is the quasiparticle associated to the excitation. The peak may oscillate before being damped, giving similarly rise to a higher energy satellite peak which is the result of the plasmon-QP interaction.

By means of the Dyson equation ansatz, the Green function of an interacting system is expressed as a function of G₀, i.e. that of the non-interacting system [205,206] as

(2.4.4) $G^{-1}\left(\mathit{r},{\mathit{r}}^{\prime },\omega \right)=G_0^{-1}\left(\mathit{r},{\mathit{r}}^{{ }^{\mathbf{\prime }}}\omega \right)+\mathrm{\Sigma }\left(\mathit{r},{\mathit{r}}^{\prime },\omega \right)-{\mathrm{\Sigma }}_{xc}^0\left(\mathit{r},{\mathit{r}}^{\prime }\right)$(2.4.4)

where the previously introduced self-energy, Σ is described as

(2.4.5) $\mathrm{\Sigma }\left(\mathit{r},{\mathit{r}}^{\prime },\omega \right)=\frac{i}{2\pi }\int e^{-i\omega { }^{\mathrm{\prime }}\delta }G\left(\mathit{r},{\mathit{r}}^{{ }^{\mathbf{\prime }}};\omega +\omega { }^{\mathrm{\prime }}\right)W\left(\mathit{r},{\mathit{r}}^{{ }^{\mathbf{\prime }}};\omega \right)d{\omega }^{\prime }$(2.4.5)

and where the screened Coulomb interaction is expressed as a function of the bare Coulomb interaction (V) and the microscopic dielectric function

(2.4.6) $W\left(\mathit{r},{\mathit{r}}^{{ }^{\mathbf{\prime }}};\omega \right)=\int V\left(\mathit{r},{\mathit{r}}^{{ }^{\mathbf{\prime }}}\right){\epsilon }^{-1}\left(\mathit{r},{\mathit{r}}^{{ }^{\mathbf{\prime }}};\omega \right)dr{ }^{\mathrm{\prime }}$(2.4.6)

On the other hand, the diagonal ansatz of the Green’s function (Lehmann) allows for a formulation where poles (E_n) and eigenstates ${\mathrm{\Psi }}_n$ satisfy the general equation:

(2.4.7) $\left(-\frac{{\mathrm{\nabla }}^2}{2}+V_{ext}+V_{Hartree}\left(\mathbf{r}\right)\right){\mathrm{\Psi }}_{\mathrm{n}}\left(\mathbf{r}\right)+\int dr{ }^{\mathrm{\prime }}\mathrm{\Sigma }\left(\mathbf{r},\mathbf{r}{ }^{\mathrm{\prime }},\omega \right){\mathrm{\Psi }}_{\mathrm{n}}\left(\mathbf{r}\right)=E_n\left(\omega \right){\mathrm{\Psi }}_{\mathrm{n}}\left(\mathbf{r}\right)$(2.4.7)

This equation motivates the choice of DFT as the best initial guess for excited state as well, since when the self-energy approaches the exchange-correlation potential of KS ($\mathrm{\Sigma }=V_{xc}^{KS}\left(\mathit{r}\right)\delta \left(\mathit{r}-{\mathit{r}}^{\prime }\right)$) it becomes equivalent to the single-particle Kohn-Sham equations.

We now briefly introduce a theoretical tool of which applications are discussed in the following and which enables the calculation of optical spectra including local-field and excitonic effects, namely the Bethe-Salpeter equation (BSE) which exploits an excitonic two-particle (electron and hole) Hamiltonian whose resonant part, Hermitian, includes QP energies calculated at the so-called GW level (named from the Green functions G and screened Coulombic interaction W), whereby an exchange term (V) that includes local field effects, and a screened Coulomb interaction term (W) enter as

(2.4.8) ${\widehat{H}}_{vck,v^{\prime }{c}^{\prime }{k}^{\prime }}=\left({ϵ}_{ck}^{GW}-{ϵ}_{vk}^{GW}\right){\delta }_{kk\prime }{\delta }_{vv\prime }{\delta }_{cc\prime }+2V_{\left(vck\right)\left(v\prime c\prime k\prime \right)}-W_{\left(vck\right)\left(v\text{'}c\prime k\prime \right)}$(2.4.8)

It is worth stressing that several flavors of the GW approach are used in applications. G₀W₀ is considered the best approach in terms of computational cost – result benefit ratio, actually requiring only the diagonal elements of the self-energy to provide G₀ (more complex is the procedure, see Eq. (2.4.4), to get the fully interacting function, G). On the other hand, it depends strongly on the initial guess, which usually is a DFT or Hartree Fock GS calculation (energy and wavefunction), which in turn is used as starting point for the Green’s function G₀ and the screened Coulomb potential, W₀. Several studies have focused on the possibility of addressing this issue (see, among the others, Ref. [207]), trying to eliminate the dependence [208] which may lead to large discrepancies between calculated and experimental values, in terms of bandgap/band edge (and relative deeper states) values. An alternative procedure to the one-shot G₀W₀ is the eigenvalues update in the Green’s function G only, keeping the screening Coulomb part as obtained from DFT (GW₀ approximation) or the full update, i.e. the Green’s function and screening matrix eigenvalues (self-consistent GW, GW or scGW). The difference between the latter approach and the simple G₀W₀ mainly consists in the Coulomb correlation that, once included in the calculation, massively reduces the difference between the two approaches [209]. An alternative way to get good results already at the G₀W₀ level is to use hybrid functionals as initial guess [210,211].

With the development of increasingly efficient computational resources, plasmonics has become subject of investigation at the many-body, beyond single-particle, level of theory. As previously discussed, the method of choice for the study of NPs remains the real-time TD-DFT. There are, however, issues related to the use of TD-DFT. The peaks due to volume plasmons in TD-DFT are often quite small/vanishing, a feature ascribed to the fact that the plasma modes are excited differently by the external perturbation. In other words, localized surface plasmons and volume plasmons will ‘communicate’ differently with the perturbation, with the latter activated via the communication with the former. GW in this sense helps in solving this shortcoming [78].

The idea of investigating in depth nanoclusters and their plasmonic features in nanoplasmonic applications is motivated by the fact that NP shape/size can be tailored, tailoring accordingly the resonance of the associated devices, further paving the way to nanoobjects with enhanced performances in optoelectronics because of an improved work frequency. Interesting in this sense is the work of Matsko [79] who at first has investigated, with the G₀W₀ method, Si clusters with size below 30 Å: such size is critical since it is believed to be intermediate between the bulk-like behavior (>30 Å) from that of clusters (<30 Å) in terms of surface and volume plasmon distinguishability. In the same work particular attention is similarly paid to the effect of passivation (with H) of Si clusters. To evaluate the plasmonic properties trend in Si NPs the spectral function of the system is considered which at the GW level reproduces two satellites for the QP peak, one per band edge, above (CB) and below (VB) the Fermi energy, respectively. The difference between the position of the two peaks is the plasma frequency which – at the GW level – is overestimated with respect to the experimental data because of the Green function poles affect the spectral function (plasmaron, i.e., a plasmon-hole interaction quasiparticle), while it does but to a lesser extent in the imaginary part of the Green function which is thus investigated and compared as a benchmark to experimental plasma frequency. An impressive improvement is found (from 28% to <0.4% deviation) when comparing the experimental plasma frequency to the theoretically calculated as imaginary part of Σ. The effect of passivation (see Figure 10 for a direct comparison between a bare and a passivated cluster in terms of spectral function, A(ω), imaginary part of the self-energy, Im(Σ), and inverse of the Green’s function, G⁻¹(ω)) is to induce a change in the frequency and damping of the plasmons with the final result of a net differentiation between surface and volume plasmons even in the smallest passivated cluster. The same feature is not observed in the bare Si NP.

Figure 10. A(ω) (solid line), Im(Σ) (dotted line), G⁻¹(ω) (dashed line) functions for bottom valence electron state (a and c) and total A_tot(ω) (b and d) for the Si₃₀ and Si₃₀H₃₄ nanoclusters. Arrows indicate the peaks corresponding to the volume and surface plasmons (VP and SP, respectively). Reproduced from arXiv:1806.08993 under Creative Commons license

More recently, Matsko has investigated with GW the formation of normal surface plasmon modes in sodium NPs [78]. The work aims to provide theoretical insight which is especially important in view of the difficulties with experimental studies of normal SPP data for small NPs. Indeed, for NPs smaller than the incident light wavefunction, it is unfeasible to match wavevector and frequency of light with those of the SPP on the object interface. This shortcoming limits the experimental knowledge of such systems, in spite of the potential significance of normal surface plasmons for applications of such nano-objects in plasmonic devices. The primary role in this theoretical analysis is played by the calculation of the loss function, Im[ɛ⁻¹(ω)], where the function $ϵ$⁻¹(r, r’,ω) that accounts for the system’s response at r to a perturbation at r’ is represented by the equation [78]

(2.4.9) $E^{-1}(r,r\prime ,\omega )={\displaystyle \sum _{q,G,G\prime }{e}^{i\left(q+G\right)r}{\epsilon }_{GG\prime }^{-1}\left(q,\omega \right)e-i{(q+G\prime )}^{r^{\prime }}}$(2.4.9)

where the inverse of the dielectric function is constructed from the response function which in turn includes the mean-field QP energy. Several clusters (from 20 to 300 atoms) were included in the analysis, while the G₀W₀ calculations were performed to solve Eq. (2.4.9)[212]. Results of such analysis are reproduced in Figure 11 where the loss function for Na clusters of increasing size is shown. In the small Na NPs ($\le$ 2 nm) only the localized surface plasmon resonance is observed, while larger NPs show the raise of the SPP with the intensity of the LSPR decreasing, indicating 2 nm as the lower limit for the devices to properly work with SPPs.

Figure 11. The loss function Im[ε⁻¹(ω)] for sodium clusters containing 20–300 atoms. Peaks corresponding to LSPR, SPP and VP excitations are marked for the upper curve, i.e. for the cluster Na₃₀₀. Reproduced from arXiv:1912.08604 under Creative Commons license

Several other papers have exploited GW to properly evaluate plasmonic properties of metallic clusters. Pavlyukh et al. [80] have calculated, by means of the combination of GW approach and Gaussian-type basis functions, lifetimes of plasmons in ${\mathrm{N}\mathrm{a}}_9^{+}$ and Pt₃ clusters (calculated from the plasmon peak’s width), finding a value of 0.7 fs, while a fast damping of the quasiparticles (11.3 and 3.4 fs) is found for the latter species, which clearly shows a more complex electronic structure where relativistic and electron correlation effects play a dominant role. Guerrini et al. [213] investigated J-aggregates of molecular crystals made of push-pull organic dyes. They used the concepts of plasmonicity index and generalized plasmonicity index (GPI), to provide thus a more ‘manageable’ definition of plasmonics. The GPI (η) is here defined as the product between the excited state lifetime (Γ⁻¹) and the plasmonic energy associated with the transition density (ρ_ξ(r))that is

$\eta ={\mathrm{\Gamma }}^{-1}{E}_{\xi }^{plas}$

(2.4.10) $E_{\xi }^{plas}=\int \frac{{\rho }_{\xi }^{\ast }\left(\mathit{r}\right){\rho }_{\xi }\left({\mathit{r}}^{\prime }\right)}{\left|\mathit{r}-{\mathit{r}}^{\prime }\right|}d\mathit{r}d{\mathit{r}}^{\prime }$(2.4.10)

Here the many-body approach is employed, first to correct, by means of the GW approach, the DFT bandgaps, while optical spectra of monomers and J-aggregates are further improved by including local-field and excitonic effects, nominally solving the BSE. From the BSE Hamiltonian, information about transition densities is obtained, thus Eq. (2.4.7) is solved and the GPI is calculated and compared between monomers and aggregates. Authors report an overall GPI lowering induced by the aggregation of molecules in the crystal.

Mowbray et al. [81] have on the other hand applied a self-developed G₀(W₀+ ΔW)-BSE theoretical setup [214,215] to the study of the interaction of light with molecules (nominally, benzene C₆H₆, terrylene C₃₀H₁₆, and fullerene C₆₀) physisorbed on a metallic substrate. The methodology employed is multi-step, treating at first the isolated molecules (DFT first and then G₀W₀ to get a reasonably good convergence for the bandgap of the molecules) and then their interactions with the surface. The core of the latter part is based on the introduction of the surface screening, process performed by correcting the quasiparticle G₀W₀ energies by means of the G₀(W₀+ ΔW) method [214,215], where $\mathrm{\Delta }W$ is the surface-induced dynamical Coulomb interaction. Aiming at including the role played by the metal surface, the bare Coulomb interaction is replaced by the screened $V+\mathrm{\Delta }W$ while the screening W in BSE is replaced with $W_0+\mathrm{\Delta }W$: the surface screening will be then included in the exchange (implicitly) and in the correlation (explicitly) terms of the self-energy. Among the cases investigated, fullerene deserves particular attention. Its π plasmon resonance is characterized by an intense peak at 6.4 eV, energy corresponding to the highest bright exciton energy in C₆₀. Such exciton (the others are at 3.9 and 5.2 eV, respectively) is mostly affected by the interaction with the metal surface. Additionally, the hybridization between C₆₀ and the metal surface seems to turn the dark exciton of the fullerene (higher in energy, at ~6.8 eV) into a bright, optically active, one. This results from several conditions that the system must fulfill. Among them the condition that the surface-plasmon frequency must be smaller, or at most comparable, with the fullerene dark exciton frequency [214]. Figure 12 illustrates the mechanism of interaction between the metal and the fullerene. After the hybridization, the surface-plasmon branch evolves as another, higher in energy, branch (ω₃) and does not continue as a π plasmon. Another branch (ω₂) emerges at the π plasmon energy. Thus, when molecular π plasmon and the surface plasmon (this is mostly effective for Ag and Au) resonate, the interaction with the surface induces three separate optically active modes.

Figure 12. Fullerene’s G₀[W₀ + ΔW]-BSE optical absorption intensity as a function of the incident x-polarized light’s energy ℏω in eV and the substrate Wigner–Seitz radius r_s in a₀ at a minimum height z₀ = 6, a₀ ≈ 3.18 Å above the substrate. The upper panel shows the excited electron (blue) and hole (red) densities schematically (left) and from G₀W₀-BSE calculations (right) for fullerene’s third bright π–π* exciton at 6.0 eV in gas phase. The exciton frequency ω_π, surface plasmon frequency ω_s, their hybridized modes ω_1,2 from Eq. (2.4.11) and the third quadrupolar exciton mode ω₃ are marked. Reproduced with permission from Ref. [81]. Copyright 2019 American Chemical Society

Formulaically, the hybridization between the exciton and the surface plasmon is obtained from a quantum system whose Hamiltonian eigenenergies are (V is the coupling energy; ω_s and ω_ex are surface plasmon and exciton frequencies, respectively):

(2.4.11) ${\omega }_{1,2}=\frac{{\omega }_s+{\omega }_{ex}}{2}\mp \sqrt{\left[\frac{{({\omega }_s-{\omega }_{ex})}^2}{4}\right]+{\left|V\right|}^2}$(2.4.11)

Topological insulators (T.I.) have been similarly extensively investigated in view of their usage in nanoplasmonic devices [216,217]. For this class of materials, indeed, a strong plasmonic response has been detected with features that seem to improve over those of Au and Ag in the UV- and blue-region[216]. Accordingly, Sb₂Te₃ [217] but also alloys like Bi_xSb_1-xTe_ySe_1-y [216] have been investigated at the MBPT level of calculation. For the former, Nechaev et al. [217] provided a theoretical comparison between DFT and GW computed Electron Energy Loss Spectrum (EELS) reporting (see above) the imaginary part of ɛ⁻¹(ω), i.e. the loss function. While results seem to be quite insensitive to the exchange-correlation potential used for the GS property calculations (DFT level), the inclusion of local field effects plays undoubtedly a major role in reproducing the experimental result. Including them, indeed, has the effect of splitting the main EELS peak, feature not detected experimentally and here ascribed to the presence of semicore states in antimonium, an effect previously reported for other systems[218]. A clear linear (parabolic) dispersion is observed for the low (high)-energy bulk plasmon that, agreeing with results for other T.I.s (i.e. Bi₂Se₃)[219], seems to be fingerprint not only of Sb₂Te₃, but of the entire class of 3D layered T.I.

Non-metal plasmonic structures

For the most commonly used metallic nanoparticles, such as gold and silver NPs, the range of realized plasmon resonance frequencies is limited by their electron density, and they typically show an LSPR resonance in the visible region. For applications, including applications in solar cells, where resonances in other parts of the spectrum are desired, one may have to go beyond metallic systems. Besides, non-metal materials have attracted some attention for their high abundance and low cost compared to noble metals. In particular, doped semiconductors offer a convenient way to control the free electron density and thereby shift the resonance to the infrared. With semiconductors, one can control the free electron density by either the dopant concentration using conventional n-doping or by injecting the electrons in the conduction band. The latter approach has the advantage of allowing dynamic control of electron concentration during device operation. Doping by surface species is also possible. The concentration of free injected electrons for the case of a near-spherical nanoparticle, which can be expressed as

(2.5.1) $n=\frac{3}{4\pi }\frac{N_e}{R^3}$(2.5.1)

where R is the effective radius, directly modifies the plasma frequency and therefore the plasmon resonance peak position as per Eqs. (2.1.1–2).

Della Sala et al. [123] studied with LR TD-DFT externally doped hydrogenated silicon clusters of different sizes up to Si₂₈₁H₁₇₂ with different free electron densities and computed absorption spectra and GPI, using the jellium model approximation. They showed that LSPR peak energy and GPI increase with the number of electrons. Tight-binding calculations in the random-phase approximation were performed in Ref. [124] to compute LSPR of silicon nanocrystals of 1.8 to 4 nm doped with large concentrations of phosphorus atoms[124]. The computed plasmon energy was higher than that predicted with the Drude model (converging to it for clusters of 4 nm), and increase in LSPR peak excitation energy with doping (free electron) concentration was clearly demonstrated. Hydrogen-terminated Si clusters of different sizes were studied with the GW approach by Matsko [79].

Another prominent example of non-metallic plasmonic materials with great potential, besides doped semiconductors, are two-dimensional materials like graphene due to remarkable confinement and light absorption properties, but due to severe energy loss in plasmonic excitations, practical applications are not yet achieved. To understand and then control energy dissipative decays of plasmonic excitations, different decay channels that can be proper to them by, for instance, electron dephasing, hot carrier mobility and phonon assisted dissipation and their underlying mechanisms must be revealed. Kuda-Singappulige et al. performed non-adiabatic ab initio Ehrenfest [220] molecular dynamic simulations and RT TD-DFT on a naphthalene molecule, a building block of graphene, to elucidate the dynamics of plasmon-like excitations and its dependence on vibrational modes [221]. In contrast to adiabatic molecular dynamics simulations based on the Born-Oppenheimer approximation [222], the non-adiabatic simulations allow non-ground state potentials to govern the nuclear motion (that is described in the classical mean-field approximation) and to include transitions between excited states [223,224]. This combination of computational techniques allowed analyzing the effects and interaction of different mechanisms through analyzing the time evolution of dipole moments, charge density differences and vibrational modes. Interestingly, plasmon-like excitations that could couple to dark electronic states (i.e. dipole forbidden) had much faster decay rates (through higher order effects).

Machine learning for modeling of surface plasmons

Machine learning (ML) has been actively penetrating many areas of materials modeling [225,226], and is only now starting to get used for the modeling of plasmonics. ML is also making inroads into experimental plasmonics [227–229]. Several methods such as neural networks (NN) [230,231], linear regression methods of various flavors[232], and K-nearest neighbors [233] have been employed to model the optical properties of plasmonic systems.

Nelson and Di Vece used a NN to find the core-shells configurations of Ag NPs with optimal optical absorption spectra for use in perovskite solar cell[231]. The NN was trained on FDTD simulations of a limited number of configurations and used to predict spectra of other configurations. Malkiel et al. [234] used multilayer (with up to 8 layers and hundreds of neurons in each layer) NNs with large numbers of neurons to map on one hand the nanostructure geometry (which was of a complex shape) to the spectrum and on the other hand (by separate NNs) the optical spectra (under light of different polarizations) and other material properties (such as permittivity) to the geometry. In this way, they machine-learned not only the direct design problem of predicting the optical response in both polarizations of a nanostructure based on its geometry, but also the important inverse design problem of finding geometries capable of providing required optical properties. The NNs were trained to finite-element based solutions of the Maxwell equations. He et al. [235] fitted a multilayer NN (with hundreds of neurons per layer) to finite-difference time domain (FDTD) solutions of the Maxell’s equations for nanospheres and nanorods and their dimers and also achieved both the forward prediction of the optical properties and the inverse prediction of dimensional parameters of NPs. Arzola-Flores and González [236] compared ridge regression (Tikhonov-regularized linear regression), a multilayer perceptron neural network, and K-nearest neighbors when machined-learning the wavelength of the dipole surface plasmon resonance (SPR) of gold concave nanocubes as a function of geometries features (edge lengths and depth radius) to data computed with discrete dipole approximation (DDA) [237], which is a simple model employing dipoles arranged in a cubic lattice that yield the response to an applied electric field. The best approach in this case was found to be the NN (with two hidden layers) with an accuracy of 94% (which was defined as the R² value between the predicted and exact values of the peak wavelength), which is roughly in line with Pearson coefficients obtained by He et al. [235] A large NN was also used by Li et al. [238] to map the geometry of a nanostructure to the spectrum by training on FDTD data.

It is expected that the application of more powerful ML techniques to datasets computed with more accurate methods, such the ab initio methods described above, will significantly enhance the modeling of plasmonics going forward. As ab initio methods are much costlier (than the rather simplistic DDA model as well as than the classical electrodynamics-based models) and therefore training datasets need to being relatively small, more powerful ML approaches such as Gaussian process regression [239] (the only method shown to outperform NN in controlled comparisons [240]) and ML methods specifically developed for working with sparse data [241,242] are expected to be helpful. A key value of ML is precisely in cutting the calculation cost which would have been incurred with more accurate but more expensive methods and be able to determine optical properties of NPs within a very large design space from a small number of reference calculations.

Conclusions

Modeling of plasmonic phenomena, and in particular of plasmonic features in the absorption spectra of nanoparticles and their arrays is important for rational design of plasmonic materials for applications in plasmonic solar cells and beyond. As most widely used plasmonic particles are metallic or heavily doped semiconductors, modeling their spectral features faces specific challenges not necessarily present in the modeling of organic materials and (weakly doped) inorganic semiconductors. While classics-based and TD-DFT approaches were previously reviewed [82–84], new applications as well as new methodologic developments surged in recent years. We therefore specifically highlighted, besides these well-known methods, emerging approaches which have the potential to significantly enhance modeling capabilities for plasmonics in the coming years, in particular, by allowing modeling at realistic length scales. These methods include time-dependent orbital-free DFT [85,86], linear response and real-time time-dependent density functional tight binding [87–90] and machine learning based approaches, as well as modeling beyond the single-particle based (DFT) picture with many-body perturbation theory [78–81].

The Casida TD-DFT formalism has been and continues to be used for the modeling of plasmonics but has not been the most popular approach to model the spectra of plasmonic NPs, in particular, because of the necessity of inclusion of a very large number of unoccupied orbitals and associated cost, as well as accuracy, issues. Recent progress in TD-DFTB partially alleviates the cost issue, making possible to include routinely tens of thousands of excited states even on a desktop workstation and more with appropriate HPC resources. Real-time approaches such as RT TD-DFT and more recently TD-DFTB found most use in plasmonics as they in principle alleviate some of the issues associated with the Casida formalism application in plasmonics. KS-DFT-based RT-DFT calculations remain, however, costly and difficult to apply to systems of realistic size. RT TD-DFTB allows routine modeling of larger (more realistic) systems and is likely to enjoy increased use especially as new DFTB parameterizations, including those tailored to plasmonics, are produced.

While DFTB allows achieving a three orders of magnitude computational cost saving vs DFT, its scaling is still non-linear with system size and the DFTB approximation imposes its own limits on accuracy. Another area which has the potential to revolutionize the modeling of plasmonics is time-dependent orbital-free DFT which brings orders of magnitude advantage in cost vs. KS-DFT-based modeling and, most importantly, in principle allows for near-linear scaling opening the way for routine simulations will millions of atoms. For some classes of materials, OFDFT is already quantitatively accurate and first results with RT TD-OFDFT suggest that it is a viable route for the modeling of plasmonic NPs of sizes too large for KS TD-DFT. As more reliable kinetic energy functionals become available, TD-OFDFT is deemed to rapidly grow in importance.

Besides advances in first principles computational method development such as TD-OF DFT and TD-DFTB that can help to bridge the lengths scales of nanoparticles between classical and TD-DFT calculations, new approaches combining computational efficiency and accuracy in robust multiscale implementations for plasmonic simulations will be essential to fully understand and reproduce experimental findings. Machine learning is only recently getting applied to the modeling of plasmonic properties, with very scarce literature to date. ML for plasmonics has the potential to experience rapid growth in the coming years, in particular in the area of rational design and accelerated discovery of plasmonic materials with desired properties.

While approaches such as TD-DFTB and TD-OFDFT have the potential to facilitate routine modeling of larger systems, the continued improvements in computing power also allows wider use of many-body approaches such as GW which palliate the critical deficiency of DFT-based approaches (i.e. their single-electron ansatz).

Acknowledgments

This work was in part supported by National Sciences and Engineering Research Council of Canada (NSERC). J.L. was financially supported by the Ministry of Science and Technology (MOST) Taiwan under Grant No. MOST109-2112-M-110-009.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada [RGPIN-2020-05725]; Ministry of Science and Technology (MOST) Taiwan [MOST109-2112-M-110-009].