Single plasmon hot carrier generation in metallic nanoparticles

Hot carriers produced from the decay of localized surface plasmons in metallic nanoparticles are intensely studied because of their optoelectronic, photovoltaic and photocatalytic applications. From a classical perspective, plasmons are coherent oscillations of the electrons in the nanoparticle, but their quantized nature comes to the fore in the novel field of quantum plasmonics. In this work, we introduce a quantum-mechanical material-specific approach for describing the decay of single quantized plasmons into hot electrons and holes. We find that hot carrier generation rates differ significantly from semiclassical predictions. We also investigate the decay of excitations without plasmonic character and show that their hot carrier rates are comparable to those from the decay of plasmonic excitations for small nanoparticles. Our study provides a rigorous and general foundation for further development of plasmonic hot carrier studies in the plasmonic regime required for the design of ultrasmall devices.


Introduction
Localized surface plasmons (LSPs) in metallic nanoparticles facilitate drastic electric field enhancements and large light absorption cross sections that can be harnessed in nanophotonic applications, such as plasmon-enhanced biosensing 1 , surface-enhanced Raman scattering 2 , data storage 3 or nanoheaters 4,5 . Recently, there has been significant interest in the decay of the LSPs into electrons and holes. The resulting carriers are energetic or "hot" and can be used in solar energy conversion applications, including solar cells 6,7 or photocatalysts [8][9][10] . For example, Mukherjee and coworkers demonstrated that hot electrons can induce challenging chemical reactions, such as the dissociation of hydrogen molecules on gold surfaces 11 . Moreover, the fast decay of LSPs can be used for new quantum information devices and in nanocircuitry 12,13 .
To provide insight and guidance in this rapidly evolving field, a detailed theoretical understanding of hot electron processes, including plasmon decay, hot carrier thermalization and recombination dynamics, is needed. Using semiclassical approaches, which combine a classical description of the LSP with a quantum-mechanical description of hot carriers, several groups analyzed the distribution of hot carriers resulting from the plasmon decay and studied its dependence on the nanoparticle size, material and environment [14][15][16][17] . Providing general insight, the semiclassical approach is frequently based on a bulk dielectric function (such as a Drude model) and therefore cannot be used to describe small nanoparticles where quantum confinement effects play an important role 18,19 . In addition, it has been observed that non-plasmonic excitations, such as electron-hole pairs, can take place at similar energies as the plasmon resonance, but such exci-tations are not described accurately on the basis of a the semiclassical approach [20][21][22] . Finally, a classical description of plasmons is less appropriate in the limit of low plasmon densities, where their quantized character must be captured [23][24][25] . Other groups have employed first-principles realtime time-dependent density-functional theory (TDDFT) to study plasmon 26,27 and hot carrier properties in small metallic nanoparticles 21,[28][29][30] . While this method allows the study of nonlinear properties and scales favorably with the system size, it does not include a quantized treatment of the plasmon. Few attempts have been made to describe the effect of electron-plasmon interactions using quantized plasmons in metallic nanoparticles. Notably, Gerchikov and coworkers 31 and Weick et al. 32 used a separation of centre-of-mass motion and relative motion of the electrons to derive a quantized electron-plasmon Hamiltonian. However, their approach cannot be used to study the decay of other neutral excitations, such as electron-hole pairs.
In this paper, we thus present a fully quantum-mechanical approach to calculating the properties of hot carriers resulting from the decay of neutral excitations, such as LSPs or electronhole pairs. In particular, we employ an effective Hamiltonian which describes the interaction of fermionic quasiparticles with bosonic neutral excitations and determine its parameters, including quasiparticle and plasmon energies and electron-plasmon coupling constants, using quantummechanical calculations. Most importantly, the electron-plasmon coupling strength is derived by comparing the electronic self energy of the effective Hamiltonian with the first-principles self energy within the GW approximation (where the electron self energy is approximated as the product of the electron Green's function G and the screened Coulomb interaction W). After identifying the dominant plasmonic and non-plasmonic neutral excitations, we first calculate the hot carrier generation rates for spherical nanoparticles with different radii and study the effect of quantum confinement on the hot carrier distributions. We find that a larger fraction of the plasmon energy is distributed to the electrons rather than the holes. Secondly, we compare hot carrier rates from the decay of plasmonic and non-plasmonic excitations. We also compare our quantum-mechanical results with semiclassical calculations and show that there is a significant discrepancy in the hot carriers rates for small nanoparticles.

Electron-plasmon coupling
To describe the interaction between charged fermionic quasiparticles and neutral bosonic excitations, such as plasmons and electron-hole pairs, in metallic nanostructures, we employ the following effective Hamiltonian whereb † I (b I ) creates (annihilates) a neutral excitation in state I with energy ω I andĉ † i (ĉ i ) creates (destroys) a quasiparticle in state i with energy i . Also, g I ij is the electron-plasmon coupling, i.e. the matrix element that describes the scattering of quasiparticles from state i into state j via the emission or absorption of a collective excitation in state I. See Supplementary note 1 for the definition of the single plasmon operator within the second quantisation formalism.
A key challenge in using H eff to describe metallic nanoparticles is the accurate determination of the various parameters, including energies of neutral and quasiparticle excitations and their coupling. Quasiparticle energies are formally defined as the poles of the one-electron Green's function and can be measured in photoemission and inverse photoemission experiments 33 . Calculating the Green's function, for example via the GW method, is very challenging for metallic nanoparticles 34 and therefore quasiparticle energies are often approximated using Kohn-Sham energies obtained from density-functional theory (DFT) 35 . Similarly, energies of neutral excitations are defined as poles of a two-particle Green's function and can be obtained by solving the Bethe-Salpeter equation or from time-dependent density-functional theory (TDDFT). Plasmons arise because of the long-ranged nature of the Coulomb interaction. This is captured by the Hartree contribution to the total energy 36 , while exchange-correlation effects often play a minor role for plasmon properties.
Neglecting exchange-correlation effects results in the well-known random-phase approximation (RPA) for neutral excitations 37 .
To determine the electron-plasmon coupling, we follow Lundqvist's approach 38 for the homogeneous electron gas and compare the second-order electron self-energy of H eff with the correlation contribution of the ab initio GW self energy, see Figure 1a. The latter is given by 39 where the coefficient η j has the value +1 for unoccupied orbitals and −1 for occupied orbitals, δ represents a positive infinitesimal and V I mj denotes the fluctuation potential given by with the transition density characterizing the I-th neutral excitation. Note that φ v (r) (φ c (r)) denote single-particle wavefunctions of occupied (empty) states and the coefficients F I vc are obtained by solving the Casida equation, see Methods section.
For the second-order electron self energy ofĤ eff , we find 40 Comparing Eq. (2) and Eq. (5) yields the final expression for the electron-plasmon coupling which is given by This result has the intuitive interpretation that the coupling between the neutral excitation I and the quasiparticles states m and j is due to the electric potential induced by the transition density of the neutral excitation, see Eq. (3). Note that the expression is completely general and does not depend on the dimensionality, material or size of the system under consideration.

Hot carrier generation
To model the decay of a neutral excitations, such as a plasmon, into hot electrons and holes, we calculate the self energy of the neutral excitation and evaluate the Feynman diagram shown in In our numerical calculations, the delta function was replaced by a Gaussian with a standard deviation of 0.12 eV.
The generation rate N I e (E) of hot electrons with energy E from the decay of the I-th neutral excitation is given by A similar formula describes the generation rate of hot holes N I h (E). In our numerical calculations, the second delta function was replaced by a Gaussian with a standard deviation of 0.05 eV.
Finally, the total generation rate N I tot of hot electrons resulting from the decay of the I-th collective excitation is given by 14 where E F denotes the Fermi energy which for metallic nanoparticles is defined as the middle of the gap between the highest occupied state and the lowest unoccupied one and δE is a threshold energy which is typically chosen larger than the available thermal energy. As we are dealing with small nanoparticles, the energy gaps between occupied and unoccupied states are always larger than the thermal energy at room temperature and we therefore set δE to zero. Note that the total rate of excited electrons is equal to the total rate of excited holes.

Semiclassical approach
To model plasmon decay and hot carrier generation in nanoplasmonic systems, a semiclassical approach is usually employed 14,15 . In this approach, the metallic nanoparticle of radius R is assumed to be exposed to an incident electric field (along the z-direction) with strength E 0 . The total potential due to the perturbing field and the induced polarization in the nanoparticle is calculated using the quasistatic approximation according to Here, (ω) denotes the dielectric function of the bulk material which is often described using a Drude model where γ P denotes the plasmon width and ω 0 is the bulk plasmon frequency ω 0 = 4πne 2 /m where e, m and n denote the electron charge, mass and density, respectively. Improved results can be obtained by using non-local approximations to the dielectric function 41,42 .
The hot-carrier generation rates are then obtained by evaluating Fermi's Golden Rule in Eq. (7) with the electron-plasmon coupling strength In order to meaningfully compare the results of the two approaches, we derive an expression for the electric field strength E 1P L 0 which is required to excite a single plasmon, see Supplementary note 2. We find that where µ P denotes the dipole moment of the plasmonic state. The corresponding semiclassical transition dipole moment is given by . When evaluated at the classical Mie frequency ω cl = ω 0 / √ 3, we obtain µ SC (ω cl ) = R 3 ω cl /µ P , which is independent of the plasmon linewidth γ P .

Numerical results
We present results for three sodium nanoparticles: Na 40 , Na 58 and Na 92 , which consist of 40, 58 and 92 Na atoms, respectively. These systems are modelled as jellium spheres with diameters of 1.4 nm, 1.6 nm and 1.9 nm, respectively. We first analyze the neutral excitations of these systems and identify those with plasmonic character and then calculate the hot carrier distributions resulting from their decay. We also compare our results to semiclassical calculations and obtain hot carrier rates for nanoparticles in different dielectric environments.  ticles under consideration. In these graphs, each circle represents an excited state with its radius being proportional to C I . For Na 40 , we find one state (denoted PL for plasmon in Figure2a) with a high oscillator strength, large collectivity and energy close to ω cl . This state gives rise to a strong peak in the absorption spectrum (see inset) and its transition density has a dipolar shape, see Figure 3. We therefore assign this state a plasmonic character. For Na 58 , we find two states with large oscillator strengths and high collectivities. The energies of both states are within 0.3 eV from ω cl and their transition densities are both plasmon-like. We therefore assign both excitations plasmonic character (and denote them by PL1 and PL2 in in Figure2b)). For Na 92 , we find two states with high oscillator strengths. The state with energy ω SP = 2.87 eV has a low collectivity and therefore electron-hole pair character (denoted SP for single pair in in Figure 2c)), while the state with energy ω P = 3.07 eV has a high collectivity and plasmonic character (denoted PL in the figure). Table 1 shows an overall redshift of the plasmonic excitation energy as the nanoparticle size increases. This is caused by quantum confinement effects and the spill-out of the electron wave functions beyond the geometrical radius of the nanoparticle 32 . However, the redshift is not monotonic which has also been observed experimentally for nanoparticles with a diameter smaller than 3.5 nm 48,49 and was predicted theoretically for small nanoparticles 32 . Moreover, the presence of more than one plasmonic resonance for a given nanoparticle could explain the experimentally observed scattered distribution of energies in this size regime 48 . Figure 3 shows the transition densities of the plasmonic excitations in Na 40 , Na 58 and Na 92 .

Identifying plasmonic excitations
The transition densities exhibit a clear dipolar shape, but are significantly more complex than predicted by classical electrodynamics. In particular, the transition densities are not only localized at the surface of the nanoparticles and show multiple regions of positive and negative charge.
While the transition density of the smallest system (Na 40 ) is extended throughout the nanoparticle,    Table 1 shows the total number of hot carriers generated by the decay of a single excitation in the three Na nanoparticles. The total generation rate for the dominant plasmonic excitation (PL in Na 40 and Na 92 and PL2 in Na 58 ) does not show a strong dependence on the nanoparticle radius.

Hot carrier distributions
It is important to recall, however, that these rates are calculated for a single excitation quantum.
As the transition dipole moment of these excitations increases with system size, see Table 1, more excitation quanta are excited per incoming photon and this gives rise to a larger number of hot carriers in larger nanoparticles.  As discussed in the Methods section, our approach can also be used to study the hot-carrier distribution resulting from the decay of neutral excitations without plasmonic character. Figure 4 shows the hot-carrier distributions resulting from decay of the prominent single-pair excitations in Na 92 (denoted SP in Figure 2 c). Similar to the plasmonic excitation in this system, the larger fraction of the excitation energy is transferred to the hot electrons. Interestingly, the total rate of hot carriers from the decay of the electron-hole pair excitation is larger than from the decay of the plasmonic excitation in this system, see Table 1.
The energy of the localized surface plasmon can be easily modified by placing the nanoparticle in different dielectric environments. Specifically, the plasmon energy is reduced as the dielectric constant of the environment increases. To study the effect of the dielectric environment on the hot-carrier generation rates resulting from the decay of plasmonic excitations, we evaluated Eq. (7) with a reduced plasmon energyω P (but keeping the orbital energies and coupling strengths unchanged). Figure 6 shows the resulting total hot-carrier generation rates for three Na nanoparticles as function of the environment-screened plasmon energy. We observe that the reduction of the plasmon energy can lead to significant enhancements in the total hot-carrier rates when the reduced plasmon energy matches the energy of multiple transitions from occupied to empty states in the quasiparticle spectrum. A particularly strong increase is found for plasmon energies smaller than half of the unscreened value ω P . This is caused by a strong increase of the electron-plasmon coupling strength for low-energy transitions, see inset of Figure 6.

Comparison to the semiclassical approach
To compare the hot-carrier rates from our quantum approach to the semiclassical approximation, we carry out semiclassical calculations using the electric field strength E 1P L 0 which generates a single plasmon in the nanoparticle. Table 1 shows the resulting total hot-carrier rate for Na 92 for the classical plasmon energy. Surprisingly, we find that the semiclassical rate is more than an order of magnitude larger than the quantum result.
To understand this discrepancy, we compare the transition dipole moment of the LSP obtained from the quantum calculation with the results from the semiclassical approach. Table 1 shows that the semiclassical transition dipole moment of the plasmon is six times larger than the quantum mechanical one. This is consequence of two factors: (i) as discussed above, the transition density of the plasmons in small nanoparticles deviates significantly from the perfect dipolar shape predicted by classical electrodynamics (see Figure 3) and (ii) in small nanoparticles, the plasmon is one of many excitations which share the total available oscillator strength (according to the f-sum rule), while in the semiclassical approach the total oscillator strength is concentrated in the plasmon, see Supplementary note 3. Because of their smaller transition dipole moments, the quantum plasmons couple less strongly to the transition dipole moments of the hot electron-hole pairs resulting in significantly reduced generation rates compared to the semiclassical calculations.

Discussion
We have presented a new approach for studying hot carrier generation in metallic nanoparticles which takes quantum plasmonic effects, such as the bosonic nature of the plasmon, fully into account. We employ an effective fermion-boson Hamiltonian and determine its parameters for specific nanoparticles using (time-dependent) density-functional theory and many-body perturbation theory within the GW approximation. In particular, an expression for the coupling strength of quasiparticles to neutral excitations is obtained by comparing the self energy of the effective Hamiltonian with the first-principles GW self energy. We have used this approach to study the decay of single plasmons into hot carriers in small sodium nanoparticles with different radii. We find that some systems exhibit multiple plasmonic excitations, while others have electron-hole pair excitations with large oscillator strengths. The hot carrier distributions from the decay of these excitations exhibit a molecular character with discrete peaks that are more closely spaced as the nanoparticle size increases. Interestingly, we find that in all systems a large fraction of the plasmon energy is transferred to the hot electrons. We also compare our results to semiclassical calculations and find that the semiclassical results provide qualitative insights but overestimate hot carrier rates for the small nanoparticles under consideration. Our approach opens the possibility to study hot carrier generation in the quantum plasmonic regime with potential application to quantumcontrolled devices, including single-photon sources, transistors and ultra-compact circuitry at the nanoscale.

Methods
In this section, we describe how the various parameters of the effective electron-plasmon Hamiltonian H eff , Eq. (1), are obtained.

Quasiparticle energies
To model the electronic structure of spherical metallic nanoparticles of radius R = r s N 1/3 (where r s denotes the Wigner-Seitz radius of the bulk material and N is the number of electrons in the nanoparticle), we employ the jellium approach which has been widely used to describe metallic clusters [50][51][52][53] . In this parameter-free method, the positive charge of the atomic nuclei is smeared out homogeneously throughout the volume of the nanoparticle.The jellium approach often yields accurate electronic properties for metals with s-or p-electron bands, but at a significantly reduced computational cost compared to full atomistic descriptions. Other authors 14,54 have used phenomenological spherical well models to describe the electronic structure of metallic nanoparticles, but we have found that such models do not reproduce the experimentally observed ordering of energy levels 55 when electron-electron interactions are included. For example, photoelectron spectroscopy of sodium clusters reveals the following ordering of energy levels (in order of increasing energy): 1s, 1p, 1d, 2s, 1f and 2p 55 . This ordering is reproduced by the jellium approach, but not by the interacting spherical well approach.
Ground state properties, such as the total energy or the electron density, are obtained by solving the Kohn-Sham equations where φ i (r), V KS (r) and i denote the Kohn-Sham orbitals, Kohn-Sham potential and Kohn-Sham energies, respectively. The Kohn-Sham potential consists of a nuclear, a Hartree and an exchangecorrelation contribution, for which we employ the Perdew-Zunger parametrization of the local density approximation (LDA) 56 .
We only carry out calculations for nanoparticles with closed electronic shells. For such systems, the Kohn-Sham potential exhibits spherical symmetry and the Kohn-Sham orbitals can be expressed as the product of a spherical harmonic and a radial function which is obtained by direct integration on a real-space grid. To approximate quasiparticle energies which correspond to electron addition or removal energies, we have calculated the ionization potential of the nanoparticles using the ∆-SCF approach (i.e., by calculating the total energy difference of the neutral and ionized nanoparticles) and shifted all Kohn-Sham energies such that the energy of the highest occupied orbital agrees with the calculated ionization potential. Alternatively, the quasiparticle energies can be obtained from GW calculations, but at a significantly larger computational cost.

Neutral excitations
Neutral excitations are obtained using the random phase approximation (RPA). In particular, we where ω I denotes the energy of the neutral excitations and F I vc is the corresponding eigenvector, which determines the transition density, see Eq. (4). The Casida matrix is given by where K vc,v c denote Coulomb matrix elements The Coulomb integrals were computed using the LIBERI library 58 which we modified to perform integrals using real spherical harmonics. Note that we have chosen to work within a linear-response framework because it allows exploitation of spherical symmetry which is broken by the perturbing where µ vc = e drzφ v (r)φ c (r) denotes the dipole moment matrix element for the vc-transition in the z-direction, which is parallel to the perturbing electric field, µ I is the transition dipole moment of the excitation I and f vc is the occupation number difference.
To define the collectivity of a neutral excitation, we first discuss the two extreme cases. When an excitation is perfectly collective, we expect that all components of |F

Supplementary Information
Supplementary note 1 The plasmon operator and its commutator The operators b † I in the effective Hamiltonian, see Eq. 1 of the main text, generate neutral excitations in state I according to 1 where |GS denotes the correlated ground state of the many-electron system that satisfies b I |GS = 0 for all I. In the random phase approximation, these operators can be expressed as linear combinations of electron-hole pairs where X and Y are the eigenvectors of the RPA pseudo-eigenvalue problem. For real orbitals, this is given by  In the above expression, we used where K vcv c is the Coulomb integral defined in the main text. The pseudo-eigenvectors satisfy the following normalization condition 2

Supplementary note 2
Derivation of the expression for the electric field required to excite one plasmon In time-dependent linear response theory for a spin unpolarized system, the change in the electron density due to the external potential δV (r, ω) = −E 0 z is given by δρ(r, ω) = dr χ(r, r , ω)δV (r , ω), where χ(r, r , ω) denotes the interacting susceptibility. The (retarded) susceptibility can be written as a sum over normal modes where ρ I (r), ω I and γ I denote the transition density, the transition energy and the linewidth of the I-th neutral excitation, respectively.
Close to the plasmon energy ω ≈ ω P , the contribution from the plasmon mode dominates the susceptibility 4 , i.e. χ(r, r , ω) ≈ ρ P (r)ρ P (r ) 1 iγ P and the induced density at the plasmon energy is given by δρ(r, ω P ) = − E 0 γ P ρ P (r) dr z ρ P (r ) ≡ − E 0 γ p µ P ρ P (r), where µ P denotes the plasmon transition dipole matrix element. According to Eq. (4), ρ P (r) is the transition density of a single plasmon excitation. Therefore, we find that a field strength hot holes (red lines) in Na 92 obtained from the semiclassical approach (dashed lines) and the new quantum-mechanical method (solid lines). Note that the latter has been scaled by a factor of 10.
The inset shows the relevant electronic transitions from occupied (blue) to empty (red) levels for the semiclassical model.
To understand the discrepancy between the hot-carrier generation rates from the quantummechanical and the semiclassical calculations, we have calculated the plasmon potentials for the two approaches. Supplementary Figures 3 and 4 (left panels) show the plasmon potential for a Na 92 nanoparticle from the semiclassical approach and quantum-mechanical theory. While the semiclassical result has a purely dipolar shape, the quantum-mechanical result has additional structure resulting from the more complicated plasmon charge distribution, see discussion in the main text.
Morever, the semiclassical potential is significantly stronger than the quantum potential. We have also calculated the product of the plasmon potential and the electron-hole pair state ψ * c (r)ψ v (r) for the transition that gives rise to the main peak of the semiclassical hot-carrier distribution in Sup-plementary Figure 2. Supplementary Figures 3 and 4 show that the overlap is significantly larger for the semiclassical potential explaining larger electron-plasmon coupling strengths and increased hot-carrier generation rates. Supplementary