II. THEORETICAL DETAILS

We consider two approaches to better understand the nature of optical excitations in nanoscale systems. We first use real-space, real-time time-dependent density functional theory (TDDFT) to obtain a quantum description of MNPs.^{Reference Castro, Appel, Oliveira, Rozzi, Andrade, Lorenzen, Marques, Gross and Rubio51} DFT theory reduces the full many-body problem to a calculation where the electronic charge density is the key material quantity that defines material properties. In DFT, a single-electron Schrodinger-like equation is solved for the Kohn–Sham states that define the electronic charge density. Although a single-electron-like equation is solved, each electron interacts self-consistently with the other electrons via the Hartree interaction, which is determined from the total electronic charge density, and via the exchange and correlation potentials, also defined in terms of the electronic charge density.^{Reference Castro, Appel, Oliveira, Rozzi, Andrade, Lorenzen, Marques, Gross and Rubio51} Here, we treat the MNPs as Au jellium nanospheres with a uniform positive charge background to model the ions.^{Reference Townsend and Bryant21,Reference Townsend and Bryant22} The jellium model has been studied extensively.^{Reference Prodan and Nordlander52,Reference Zuloaga, Prodan and Nordlander53} Most recently, the jellium background has been replaced by the lattice of ionic cores to account for atomistic effects as well.^{Reference Zhang, Feist, Rubio, García-González and García-Vidal54} In this study, our goal is to show how to characterize the optical excitations and to distinguish single-particle and collective excitations. The simpler jellium model is sufficient for this purpose. We consider the limit of small MNPs where size quantization and the finite number of electrons play a key role and show how these effects evolve as the nanoparticle size increases.

The TDDFT calculations proceed in three steps. First, DFT is used to determine the electronic ground state of the jellium MNP.^{Reference Castro, Appel, Oliveira, Rozzi, Andrade, Lorenzen, Marques, Gross and Rubio51} The Schrodinger-like equation is solved to determine the single-particle-like Kohn–Sham orbitals (${\phi _i}$, with orbital energy ε_i) that are occupied for ε_i less than the Fermi level ε_F to determine the ground-state electron charge density. Then, TDDFT is used to determine the response of the MNP to an instantaneous electromagnetic delta-impulse by evolving each occupied Kohn–Sham orbital in time after the impulse. The time-dependent dipole moment of the MNP is determined. Its Fourier transform gives the frequency-dependent response of the system. A typical spectrum for a 100-electron MNP is shown in Fig. 1. The width of each peak is determined by the length of the time evolution. Once the excitation frequencies are determined from the spectrum, the TDDFT time evolution is repeated for electromagnetic driving fields at the excitation frequencies identified from the spectrum. This allows us to explicitly investigate the spatial character of each excitation as it evolves, define which single-electron transitions between Kohn–Sham orbitals contribute to the excitation, and determine how these contributions evolve in time. As we will demonstrate, this information allows us to distinguish single-particle and plasmonic (i.e., collective) excitations. However, TDDFT does not describe how these material excitations should be quantized, nor does it describe whether the excitations are fermions or bosons, or whether multiple plasmonic excitations are independent, like bosons, or coupled.

FIG. 1. Optical response of spherical Au MNPs. (a) Frequency-dependent response for a 100-electron Au MNP. Frequency on the horizontal axis is displayed in units of the classical surface plasmon frequency. (b) Time-dependent induced charge density for the classical surface plasmon [ω = 0.90ω_sp, as indicated by the circle in (a)] of a 100-electron MNP. The color shows the induced charge along the axis (x) parallel to the driving field through the center of the MNP. R is the MNP radius. The bottom (blue) line shows the phase of the driving electric field. The solid line with circles is the ground-state electron density as a function of position x. The ground-state electron density is found in a DFT calculation and then used as the starting density for the TDDFT time simulations. The scale for the ground-state electron density, normalized by the density of the jellium background, is shown at the top. T is the period of the driving field. (c) Time-dependent induced charge density for the quantum core plasmon (ω = 0.74ω_sp) of a 100-electron MNP. (d) Frequency-dependent response for a 600-electron Au MNP.

To address these quantization issues without imposing the form of the quantization as constraints the theory must obey, we need a full quantum mechanical theory of these excitations that finds the actual many-body eigenstates of the system rather than just the charge density, as in DFT. With the full eigenstates, one can determine how they are quantized. To gain some insight into what such a theory might reveal, we have studied some simple models for interacting electrons in small nanoscale systems. The goal is to study systems small enough that they can be analyzed exactly. This limits the models to systems with a small number of electrons. Here, we study a linear chain of atoms, with electrons hopping from atom to atom along the chain. Such a model describes atomic chains on surfaces^{Reference Wallis, Nilius and Ho4,Reference Yan, Yuan and Gao5} and linear molecules.^{Reference Bernadotte, Evers and Jacob2,Reference Gao, Ruud and Luo3} Most importantly it should provide insight on how plasmonic excitations arise in these systems.

Here, we consider a single-band model for the linear chain of atoms with one electronic state on each atomic site and short range, next-neighbor hopping t between atoms. If there are N _s sites in the chain and m _e electrons in the system, then there are $\left( {\matrix{ {{\it N}_{\rm{s}} } \cr {{\it m}_{\rm{e}} } \cr } } \right) = {{{\it N}_{\rm{s}} !} \over {{\it m}_{\rm{e}} !\left( {{\it N}_{\rm{s}} - {\it m}_{\rm{e}} } \right)!}}$ many-electron states in the system. In a charge neutral system, each site has a nuclear charge Z = m _e/N _s. An electron at site i interacts with the nuclear charge Z at site j via the attractive Coulomb interaction

(1)$${V_{{\rm{nuc}}}}\left( {i,j} \right) = {{ - {\lambda _{{\rm{nuc}}}}Z} \over {\left( {\left| {i - j} \right| + {\xi _{{\rm{nuc}}}}} \right)}}\quad ,$$

where λ_nuc is a scale factor that includes any dielectric screening and the length scale for the site separation and ξ_nuc is a cutoff that accounts for the orbital spread of the electron orbital on a site. An electron at site i interacts with an electron at site j via the repulsive Coulomb interaction

(2)$${V_{{\rm{ee}}}}\left( {i,j} \right) = {{{\lambda _{{\rm{ee}}}}} \over {\left( {\left| {i - j} \right| + {\xi _{{\rm{ee}}}}} \right)}}\quad ,$$

where λ_ee is the scale factor that includes any dielectric screening and the length scale for the site separation and ξ_ee is the cutoff that accounts for the spread of the electron orbital on a site. The effects of electron–electron interaction are balanced by the electron–nuclear interaction when λ_ee = λ_nuc. Optical response in this model is proportional to the dipole moment of the electrons along the chain.

All electronic states for this model can be easily determined by standard numerical diagonalization techniques when the system has less than ten thousand states. Here, we consider a half-filled band of spinless electrons (m _e = N _s/2). In that case, we can consider chains with up to 16–20 atoms. For longer chains or when electron spin is included, it becomes impractical to determine all the many-electron states. Instead, one must use diagonalization techniques that find limited parts of the spectrum. Here, we focus on the smaller chains where all excitations can be found and analyzed. This allows us to investigate the characteristics that will distinguish plasmonic-like, collective excitations from single-particle excitations. Developing this insight will allows us to study larger systems where we must start with educated choices about which selected parts of the spectrum to investigate to find the plasmon-like excitations.

III. DENSITY FUNCTIONAL THEORY OF METAL NANOPARTICLES

Collective/plasmonic excitations can be strongly mixed with single-particle excitations in nanoscale systems, making it difficult to characterize the excitations. We consider a 2-nm-diameter Au MNP with 100 electrons to show how the character of its excitations can be defined.^{Reference Townsend and Bryant21,Reference Townsend and Bryant22} Figure 1(a) shows the spectral response of the 100-electron gold MNP found using TDDFT and the jellium model. The spectrum has well separated, narrow peaks, several with large oscillator strengths, and the rest substantially weaker. There are no losses in the calculation so the width of each peak is determined by the length of the time simulation. The time evolution is long enough to ensure that peaks are resolved and that each resonance frequency can be accurately determined. Typically, the simulations are run for tens of periods of the driving field. This is comparable to the plasmon lifetime, although the sharpest spectra shown here would be broadened by typical losses and inhomogeneous effects. To reveal the nature of each excitation, we drive the MNP at the excitation resonance frequency. The time evolution of the induced charge density when driven at the two indicated frequencies is shown in Figs. 1(b) and 1(c). The induced charge density oscillation for the most prominent peak is localized mostly at the MNP surface, behaving as one would expect for the classical surface plasmon. There is structure in the induced charge density near the surface that arises due to the quantum confinement. The excitation illustrated in Fig. 1(c) also shows an oscillating induced charge density, but it is localized to the core of the MNP. We have previously referred to this as a quantum core plasmon, because it is a plasmon-like internal excitation that arises because of confinement effects. The dependence on size is shown by comparing the response of a 600-electron MNP [Fig. 1(d)] with the response for the 100-electron MNP. In the larger MNP, the response of the core plasmon is much weaker and not easily identifiable. For the 100-electron MNP, there are two additional strong excitations between the core and surface plasmons that have mixed character of both the core and surface plasmons. In addition, there are a large number of weak resonances, especially for the smaller MNP. These resonances are sensitive to the boundary conditions used to define the computational region for the DFT. This indicates that these excitations have higher-energy electrons that can escape the MNP and are sensitive to the details of the computational boundary. The energy and spatial character of the classical surface and quantum core plasmons do not depend on the boundary conditions, indicating that these are internal excitations of the MNP.

It is tempting to conclude that this analysis has determined that the prominent excitations are plasmonic. However, further analysis will show that the core plasmon is really predominately single-particle-like. This is already hinted at by the size dependence. The core plasmon is much weaker than the surface plasmon in the larger dot, indicating that the strength of the core plasmon does not scale with the size of the system, as it would if it were a plasmon-like, collective, many-electron excitation.

To get a clearer picture of which excitations are collective and which are single-particle-like, one must determine which single-electron transitions contribute to each excitation and how these contributions evolve as the MNP is driven on resonance. In a time simulation of a driven MNP, the MNP starts at t = 0 in its DFT ground state with the DFT Kohn–Sham orbitals for ε_i < ε_F occupied. These occupied orbitals change in time when the driving field is applied for t > 0. The projection of the time-evolving occupied orbital, ${\phi _i}\left( t \right)$, onto its initial, t = 0, orbital gives the probability ${\left| {\left\langle {{\phi _i}\left( t \right)} \right|\left. {{\phi _i}\left( 0 \right)} \right\rangle } \right|^2}$ that ith state has evolved without changing its initial state. The probability that the state undergoes a transition to another state j is given by the off-diagonal projection ${\left| {\left\langle {{\phi _i}\left( t \right)} \right|\left. {{\phi _j}\left( 0 \right)} \right\rangle } \right|^2}$. Many of these single-electron transitions could contribute to an excitation. How they contribute will define the character of the excitation. Figure 2 shows these transition probabilities at one time when the 100-electron MNP is driven continuously at the surface plasmon resonance. The off-diagonal transitions are small, in the 10⁻⁴ range. The decrease in the probabilities to remain unchanged in the ground-state Kohn–Sham state is in the same range. The total induced charge density is small and the response is well within the linear regime. The orbital symmetries of the states in the transitions are indicated on the axes in Fig. 2. The Fermi level is in the 2f level. As expected for optical excitation, the allowed transitions have a change of orbital angular momentum of ±1. Four types of transitions are highlighted. Strong transitions from occupied states just below the Fermi level to empty states just above the Fermi level have the appropriate change in angular momentum. However, the change in energy due to these single-electron transitions is much smaller than the excitation energy, indicating that these transitions must contribute collectively to the excitation. As the total occupation of the ground-state Fermi sea is depleted by the driving field, a second class of transitions between states inside the Fermi sea develops and agitates the Fermi sea. At the same time, charge accumulates above the ground-state Fermi sea and a third class of transitions between states above the Fermi level develops. These two classes of transitions that arise because the ground-state Fermi sea is depleted also have energies well below the excitation energy and must contribute collectively to the excitation. Finally, a fourth class of transitions couples states well below the Fermi level to states well above the Fermi level. For these transitions, the energy changes are comparable to the excitation energy indicating that these transitions add single-particle-like character to the excitation. These are transitions that can also contribute to the hot electron generation mentioned previously.

FIG. 2. Projections ${\left| {\left\langle {{\phi _n}\left( t \right)} \right|\left. {{\phi _m}\left( 0 \right)} \right\rangle } \right|^2}$ of time-evolved Kohn–Sham states for the classical surface plasmon onto the ground-state Kohn–Sham states in the 100-electron MNP. The snapshot is shown for an arbitrarily chosen time. The circles highlight the four classes of transitions discussed in the text.

The transition probabilities for collective and single-particle contributions have distinctly different time dependences. For the class of transitions that makes single-particle-like contributions, the transition probabilities increase monotonically as the MNP is driven for the time-scale of the simulations. This monotonic increase is expected from Fermi's golden rule for single-particle transitions weakly driven well below saturation. For the three classes of transitions that contribute collectively, the transition probabilities initially increase on average as the MNP is driven, but this average increase is strongly sinusoidally modulated at the resonance frequency and by low frequency beating. These differences in the time dependence for the single-particle and collective transitions are also apparent in the time dependence of the occupation of each ground-state level.^{Reference Townsend and Bryant22} Levels just below the Fermi level deplete and then refill. This oscillation in level occupation occurs at the resonance frequency. In addition, charge accumulates in states just above the ground-state Fermi level and then empties from those states, in phase with the depletion and refilling of the states below the Fermi level. This collective sloshing of charge back and forth from just below the Fermi level to just above the Fermi level is the hallmark of a plasmon-like excitation. At the same time, single-particle-like transitions from states well below the Fermi level to well above the Fermi level deplete the levels well inside the Fermi sea monotonically, leading to an inversion of charge well above the Fermi level. For the 100-electron MNP, all four classes of transitions contribute to each resonance because the resonances are strongly mixed. Nonetheless, the surface plasmon resonance is dominated by transitions that slosh charge back and forth from just below to just about the Fermi level. In contrast, the resonance labeled a quantum core plasmon is dominated by the transitions that invert the charge from well below to well above the Fermi level. From this analysis, it becomes clear that the quantum core “plasmon” is predominantly a single-particle excitation.

The spatial distribution of induced charge associated with each class of transitions further defines the character of the excitations. The sloshing transitions induce charge density oscillations at the surface. Transitions between states both inside the Fermi sea or both above the Fermi sea make little contribution to the induced charge oscillation of an excitation. The single-particle-like transitions contribute to the induced charge density oscillations inside the MNP, either to structure the charge density oscillation near the surface or to define the oscillations deeper inside the MNP that characterize the core plasmon.

The TDDFT analysis shows that the single-particle and plasmonic transitions in nanoscale systems can be identified and distinguished even though they can be strongly mixed. To determine the nature of an excitation, both its time-dependent induced charge oscillation and the time dependence of the transitions that contribute to it must be characterized. Single-particle-like excitations have contributions from single-particle transitions with energy changes comparable to the excitation energy. These contributions increase monotonically with time and occur mostly inside the MNP. The collective surface plasmon-like excitations arise from charge sloshing back and forth from just below to just above the Fermi level. Many low-energy transitions contribute collectively to this excitation.

IV. ONE-DIMENSIONAL CHAIN MODELS FOR NANOPLASMONICS

The TDDFT analysis says nothing about how to quantize the excitations in nanoscale and atomic-scale systems, nothing about whether the excitations are bosonic or fermionic, and nothing about whether the excitations are independent or coupled nonlinearly. Exact solutions possible for small systems provide a framework to address these questions. However, the first step in this approach is still to see if the excitations have well-defined characteristics that can distinguish them as plasmonic or single-particle. In this section, we discuss initial exact results for short linear chains of atoms to see what insight can be gained. We only consider here chains with fewer than 20 atoms where exact results for the entire spectrum can be obtained. This allows us to examine the full spectrum to find the signs of plasmonic excitations.

As we have already discussed, excitations in nanoscale systems can be strongly mixed with both single-particle and plasmonic character. Here, we follow an approach recently developed by Bernadotte et al.^{Reference Bernadotte, Evers and Jacob2} to identify single-particle and plasmonic excitations in TDDFT calculations for short linear molecules. They found that many excitations could be identified as single-particle or plasmonic. To make this identification, they varied the scale of the Coulomb interaction in the TDDFT using the same scale factor λ_ee that we use to define the electron–electron repulsion in our model. In their TDDFT calculations, single-particle-like excitations had excitation energies that depended only weakly on the Coulomb scale factor λ_ee. In contrast, other excitations had excitation energies that scaled with the square root of the Coulomb interaction strength. This square root scaling is expected for plasmonic excitations. When these excitations crossed, they had mixed character.

We have followed this approach for our exact calculations. Here, we illustrate the results by discussing chains with 12 atoms and 6 electrons. This corresponds to a half-filled band of spinless electrons and should have metallic behavior. The full excitation spectrum, with nearly a thousand states, is shown as a function of λ_ee in Fig. 3. All excitation energies are shown relative to the ground-state energy. Qualitatively similar spectra are obtained for eight-atom chains with 4 electrons and 70 excitations and for 16-atom chains with 8 electrons and ten thousand excitations. The spectra for different chain lengths have a similar energy span. The main difference is the density of states. At λ_ee = 0, there is no electron–electron repulsion and no collective response. All excitations are single or multiple independent-electron excitations. As λ_ee increases, the electron–electron interaction is turned on, the excitations become correlated and collective response is possible.

FIG. 3. The spectra of a 12 atom, 6 electron linear atomic chain as a function of the electron–electron interaction strength λ_ee. The excitation energies ΔE are shown relative to the ground-state energy. The inset shows spectra for small λ_ee.

The inset in Fig. 3 shows the excitation spectra for λ_ee < 0.5t, where the electron–electron repulsion is comparable to the hopping. A few low-energy and high-energy excitations are well separated from the quasicontinuum of excitations at intermediate energies. At low energies, these excitations correspond to single-electron excitations in a short chain. The quantum confinement of a finite chain length ensures that the single-electron states have discrete energies. At high energies, the discrete excitations are those formed by removing a single-electron excitation from the fully excited state. Figure 4(a) shows a blowup of the low-energy spectra. The excitation energies increase weakly as the electron–electron repulsion is added. The excitation energies for a few of the states increase faster than others. However, there is no clear separation of states, as predicted by TDDFT, into a set of single-particle states which are nearly independent of λ_ee and another set of states with excitation energies that scale as the square root of λ_ee. An analysis of the low-energy excitations shows that they are single independent-electron excitations. In this range of λ_ee, interaction is too weak to have much effect.

FIG. 4. The low-energy spectra of a 12 atom, 6 electron linear atomic chain as a function of the electron–electron interaction strength λ_ee for (a) weak electron–electron interaction and (b) stronger interaction where collective effects develop. The metallic and Wigner crystal regimes and the regimes of weak and strong Coulomb effects are indicated.

This weak dependence on λ_ee might seem counterintuitive because the electron–electron repulsion is comparable to the hopping. However, the full Coulomb interaction includes both the electron–electron repulsion which leads to collective response and the electron–nucleus attraction which contributes to the independent single-electron energies. Because the chains are chosen to be charge neutral, these interactions compensate each other and the net Coulomb interaction is weak. The spectra for λ_ee increasing up to 25t are shown in Fig. 3. The low-energy excitations for this range of λ_ee are shown in Fig. 4(b). For low λ_ee, the ground-state electron density is uniform along the chain, as would be expected for a metallic state. For very large λ_ee (λ_ee ≈ 25t), the first excited state and the ground state become degenerate. For these λ_ee, the ground-state electron density becomes localized to every other site, as would be expected for a Wigner crystal where the electrons separate and localize into an array due to their mutual repulsion. For a chain with an even number of sites, there are two ways to realize this ground-state density and the ground state becomes doubly degenerate. This suggests that metallic behavior extends up to λ_ee ≈ 25t and the onset of Wigner crystallization. In addition, the highest two excitations become degenerate for λ_ee > 1.5t. The highest energy state becomes fully localized to one half of the chain. Because there are two ways for this to happen, the highest excitation becomes doubly degenerate for λ_ee > 1.5t. For smaller λ_ee, even the highest excitation has a uniform electron density along the chain. These results suggest that strong collective response begins to appear for λ_ee > 1.5t but metallic behavior persists until λ_ee ≈ 25t. These regimes are indicated in Fig. 4. In this range of λ_ee, the low-energy excitations remain well separated from the quasicontinuum and show a clear strong linear dependence on λ_ee. An analysis of these states shows that they are formed from many single and multiple independent-electron states, indicating a strong collective character.

To better identify the optically active single-particle-like and plasmonic/collective excitations, we eliminate from the spectra all excitations that are dipole forbidden. This eliminates all excitations that have the same parity as the ground state and eliminates all excitations made up only from independent-electron states with multiple-electron excitations. This dramatically reduces the density of states in the spectra, as shown in Fig. 5. The largest (blue) circles indicate excitations with the largest oscillator strength. Circles of smaller size (respectively green, red, and gray) have oscillator strengths reduced by one order of magnitude for each size step. The smallest (black) dots indicate all other excitations with finite oscillator strength.

FIG. 5. Optically active excitations in the low-energy spectra of a 12 atom, 6 electron linear atomic chain as a function of the electron–electron interaction strength λ_ee for (a) weak electron–electron interaction and (b) stronger interaction where collective effects develop. The largest (blue) circles indicate excitations with the largest oscillator strength. Circles of smaller size (respectively green, red, and gray) have oscillator strengths reduced by one order of magnitude for each size step. The smallest (black) dots indicate all other excitations with finite oscillator strength.

At λ_ee = 0, only the single-electron excitations are optically active. They occur for excitation energies less than 4t. In the low interaction regime for 0 < λ_ee < 0.5t, these single-particle-like excitations lose oscillator strength as λ_ee increases while higher-energy excitations above the single-electron transitions become optically active. Several of these higher-energy states in the quasicontinuum have higher oscillator strength than nearby states, suggesting that they might be plasmonic/collective states. However, the character of these states changes rapidly as λ_ee increases because of multiple level crossings and strong mixing with the other nearby states in the high density of states. These higher-energy regions in the spectra with strong oscillator strength might correspond to plasmonic-like resonances but they will be short lived. As λ_ee increases further and correlated response becomes important [Fig. 5(b)], the low-energy excitations which are single-particle-like for λ_ee < 0.5t become correlated, collective excitations made from many single- and multiple-electron transitions. The excitation energies increase linearly with λ_ee both at small λ_ee and at large λ_ee with a crossover at intermediate λ_ee. At higher λ_ee, these appear to be plasmon-like collective excitations. There is no indication that collective and single-particle-like excitations can both be present in the spectra at the same λ_ee as predicted by TDDFT.

The predictions of the exact calculations are in striking contrast to the TDDFT predictions for linear chains.^{Reference Bernadotte, Evers and Jacob2} TDDFT predicts that single-particle-like and collective/plasmonic excitations can coexist at the same λ_ee. Excitation energies for single-particle excitations vary weakly with λ_ee. The TDDFT^{Reference Bernadotte, Evers and Jacob2} also predicts that plasmonic excitation energies have a square root dependence on λ_ee. The exact calculations presented here show that the single-particle-like excitations evolve into the collective excitations as λ_ee increases with a crossover region at intermediate λ_ee. At low and high λ_ee the excitation energies depend linearly on λ_ee. There is no indication that well-defined single-particle-like excitations and collective, plasmon-like excitations coexist at the same interaction strength. These different predictions clearly present a challenge that must be resolved to develop an understanding for quantum excitations in these systems.