Atomistic Theory of Hot-Carrier Relaxation in Large Plasmonic Nanoparticles

Recently, there has been significant interest in harnessing hot-carriers generated from the decay of localized surface plasmons in metallic nanoparticles for applications in photocatalysis, photovoltaics, and sensing. In this work, we develop an atomistic method that makes it possible to predict the population of hot-carriers under continuous wave illumination for large nanoparticles of relevance to experimental studies. For this, we solve the equation of motion of the density matrix, taking into account both the excitation of hot-carriers and subsequent relaxation effects. We present results for spherical Au and Ag nanoparticles with up to 250,000 atoms. We find that the population of highly energetic carriers depends on both the material and the nanoparticle size. We also study the increase in the electronic temperature upon illumination and find that Ag nanoparticles exhibit a much larger temperature increase than Au nanoparticles. Finally, we investigate the effect of using different models for the relaxation matrix but find that the qualitative features of the hot-carrier population are robust. These insights can be harnessed for the design of improved hot-carrier devices.


Introduction
Metallic nanoparticles (NPs) feature strong and highly tunable light-matter interactions which makes them attractive for applications in photovoltaics and photocatalysis [1][2][3][4][5] .In particular, such NPs exhibit collective charge oscillations known as localized surface plasmons (LSPs) 6 which give rise to a pronounced peak in the optical absorption spectrum at the localized surface plasmon resonance (LSPR) frequency.LSPs are strongly damped oscillations that can decay via the Landau damping mechanism into energetic electron-hole pairs, known as hot carriers, which can then be harnessed for different applications.
After their generation 7 , hot carriers lose their energy as a result of interactions with other electrons, lattice vibrations 3,8,9 or defects.In addition, they can also interact with molecular species adsorbed to the surface of the NP which in turn can induce adsorbate vibrations, desorption events or chemical transformations.A detailed understanding of how these various processes influence and depend on the population of hot carriers is important for optimizing devices.
To gain microscopic insights into hot-carrier relaxation processes, various theoretical approaches have been developed.For example, Dubi and Sivan 10 solved the Boltzmann equation to calculate the steady-state population of hot carriers taking both electron-electron and electron-phonon scattering effects into account.A fully quantum-mechanical approach based on the equation of motion for the electronic density matrix was developed Govorov and coworkers 11,12 who employed a relaxation time approximation.Another way to include relaxation effects is by solving a master equation for the electronic occupations of the NP states.Electron-phonon, electron-electron and electron-photon scattering mechanisms can be included, and their relative importance estimated 13,14 .
While the aforementioned modelling efforts allowed important insights into the relaxation dynamics of hot carriers, it is important to note that they all employed highly simplified models of the NP electronic structure: they often use simple spherical well models or assume that the nanoparticle's density of states is the same as that of a bulk free electron gas.As a consequence, such models cannot capture hot carriers generated in d bands 15 , describe the dependence of the electronic structure on the exposed facets of the NP 16 , or describe plasmonic heterostructures atomistically 17 , all of which can play a critical role [18][19][20][21][22] .
In this paper, we start from the density matrix formalism developed in 11 and combine it with an atomistic model of the NP's electronic structure based on the tight-binding approach.
To enable the simulation of large NPs, we cast the expression of the hot-carrier population in a basis-independent way and use the Kernel Polynomial Method (KPM) 23,24 to evaluate it.
In recent years, KPM and other spectral methods have been established as promising tools to study quantum-mechanical phenomena in mesoscopic systems.For example, they have been used to extract spectral properties [25][26][27] and transport properties [28][29][30][31][32][33][34] in multi-billion atom lattices.Recently, they have also been leveraged to obtain the rate of hot-carrier generation in million-atom plasmonic nanoparticles 35 .Here, we use this approach to study the hot-carrier population in large spherical Au and Ag nanoparticles.
The paper is organized as follows.In Section 2 we explain how the optical properties of the NP, the electronic Hamiltonian and the density matrix are constructed before using these quantities to determine the non-equilibrium steady-state hot-carrier population.In Section 3, we discuss the sensitivity of the hot-carrier population on the values of the relaxation times and then study the dependence of the hot-carrier population on the diameter of spherical Au and Ag NPs.This section ends with a discussion of the temperature increase in illuminated NPs.

Theory and computational details
In this section, we describe our theoretical and computational methodology.For further details, see Ref. 35 .We start with the optical properties of the NP, discuss the electronic Hamiltonian and density matrix and finally show how the steady-state population can be obtained in a basis-independent way.

Optical properties
The interaction of the NP with light is modeled by considering the effect of a monochromatic external electric field of angular frequency ω given by E ext (r, t) = E 0 cos (ωt − kz) ẑ, where kẑ denotes the wave vector and E 0 is the amplitude.In this paper we are interested in NPs whose diameters range from ∼ 2 nm to ∼ 20 nm, which are much smaller than the typical wavelength of the external electric fields.In this case E ext can be considered uniform inside the NP.This is the quasistatic approximation (QSA).To calculate the total electric field experienced by the electrons inside the NP, it is then sufficient to solve Laplace's equation ∇•(ε(r, ω)∇ϕ(r, ω)) = 0 for the total electric potential ϕ, with boundary condition at infinity ϕ (r, t) = −E 0 z cos (ωt) and appropriate dielectric functions ε(ω) in each region.Here we are only considering spherical nanoparticles for which the analytical solution of Laplace's equation is given by where ε m is the dielectric constant of the environment and the measured dielectric function of the bulk material (either Ag or Au) were taken from Ref. 36 .Here, we consider NPs in vacuum and set ε m = ε 0 , the vacuum permittivity.The factor reflects the enhancement of the electric field inside the NP.The potential energy Φ (r, t) of an electron at position r and time t is given by Φ (r, t) = −eϕ (r, t), where −e < 0 is the electron's charge.This optical perturbation excites the electron-hole pairs which then interact with phonons and other electrons and relax to a non-equilibrium population under continuous wave illumination.

Electronic Hamiltonian
The dynamics of electrons is described by a tight-binding model H = H 0 + Φ (t), where represents the Hamiltonian in the absence of illumination and describes the optical perturbation.Here, c † i,α (c i,α ) denotes the creation (annihilation) operator of an electron in orbital α of an atom at position R i , t αβ ij is the hopping and onsite matrix and the double sum over i and j is typically constrained to a small number of neighbours, making t αβ ij a sparse matrix.The tight-binding models used in this paper are obtained through an orthogonal two-center Slater-Koster parameterization 37 .For Au, we use a parameterization consisting of 5d, 6s and 6p atomic orbitals and for Ag, we use 4d, 5s and 5p orbitals.

Density matrix
To determine the non-equilibrium population of electrons inside the NP, we solve the equation of motion for the electronic density matrix ρ (t) given by 11 where |n⟩ denotes an eigenstate of H 0 with corresponding eigenenergy ε n and ρ mn = ⟨m|ρ|n⟩.
Relaxation effects are described phenomenologically through the matrix Γ mn .Also, δρ mn = ρ mn − ρ 0 mn , where ρ 0 mn = δ nm f (ε n ) is the density matrix before the perturbation has been turned on, given by the Fermi-Dirac distribution f (ε) = 1 + e β(ε−µ) −1 , with chemical potential µ and β = 1/(k B T ) (with k B and T denoting the Boltzmann constant and the temperature, respectively).For all our calculations, we use room temperature (T = 298K).
The equation of motion for the density matrix can be solved using perturbation theory assuming a weak monochromatic perturbation Φ (t) = φ ω e iωt + φ † ω e −iωt 11 .The result can be used to calculate the excess density of electrons with energy E, δN relative to the equilibrium value.This yields Evaluating this expression requires knowledge of the eigenfunctions and eigenvalues of the Hamiltonian H 0 .Calculating these quantities using a diagonalization procedure becomes numerically challenging for large NPs of relevance to nanoplasmonic experiments.

Basis-independent population
To study large NPs, it is useful to express Eq. 1 in a basis-independent form.To achieve this, we first assume that Γ can be expressed as a function of the eigenenergies, i.e.Γ nm = Γ (ε n , ε m ).Then, following 35 and introducing a double integral over ε and ε ′ weighted by with Note that Λ (ε, ε ′ ) contains all information about relaxation processes through Γ.The Lorentzian form of the term inside the square brackets ensures energy conservation in optical transitions from ε to ε ′ with a finite linewidth Γ (ε, ε ′ ).The difference of Fermi functions ensures that only transitions from occupied states to unoccupied state are allowed.Ω (ε, ε ′ ) is the energy-resolved optical transition matrix.Importantly, Eq. 2 factorizes the information about statistics and relaxation processes in Λ (ε, ε ′ ) and the information about the singleparticle electron dynamics in Ω (ε, ε ′ ).Typically, Ω (ε, ε ′ ) is the most difficult object to calculate, but once it is obtained, any relaxation matrix can be studied without additional computational cost as a consequence of the factorization.Ω (ε, ε ′ ) is calculated efficiently using the Kernel Polynomial Method 23 in the same way as in Ref 35 .

Results and Discussion
In this section, we calculate the steady-state electron population in spherical Au and Ag NPs of different sizes.We also explore different simple models for the relaxation matrix, but find that the gross features of the electron population do not depend on the model for the relaxation matrix.

Relaxation matrix model
The relaxation matrix Γ (ε, ε ′ ) is one of the central objects in Eq. 2, and its functional form depends on the details of electron-electron and electron-phonon scattering mechanisms.
While our approach allows a considerable amount of flexibility in choosing Γ, here we explore a few simple models to understand their effect on the steady-state electron population of spherical Au and Ag NPs.Following 11 , we parameterize the matrix by the energy relaxation time τ ε (τ ε = 1.3 meV for both Au and Ag) and the momentum relaxation time τ p (τ p = 78 meV for Au and τ p = 20 meV for Ag).The momentum relaxation time is obtained from the Drude fit to the optical constants of Au and Ag obtained experimentally 38 .The energy relaxation time estimates the electronic thermalization time and has considerable contributions coming from both electron-electron and electron-phonon scattering.Pumpprobe experiments provide information about both these processes and place the corresponding times in the order of magnitude of 0.1 − 1 ps [39][40][41] , depending on the experimental apparatus, NP sizes and frequency at which the measurement was made (among other factors).For simplicity, we consider their collective effect to yield τ ε ≈ 0.5 ps.In 11 , Γ nn = τ ε and To mimic this form for the relaxation matrix, we assign τ ε to transitions which have very similar energies, that is |ε n − ε m | < 20 meV, and assign τ p otherwise.This relaxation matrix is referred to as Γ 1 and given by To test the sensitivity of our results to the numerical values of the relaxation times, we also present results for a relaxation matrix in which the value of τ p is doubled (denoted by Γ 2 ) and a relaxation matrix in which the value of τ ε is doubled (denoted by Γ 3 ).Comparing the results of the three different relaxation matrix we observe that the energy relaxation time mostly determines the overall magnitude of δN : the populations obtained with Γ 1 and Γ 3 (which have the same momentum relaxation time) have the same shape, but differ in their magnitude.This can be understood from Eq. 1, since 1/Γ nn = 1/τ ε enters as an overall prefactor.In contrast, the shape of δN is affected by the momentum relaxation time, especially around the Fermi level.However, we find that the qualitative features of the electron population are the same for all relaxation models and we therefore present results only for Γ 1 in the remainder of this paper.

Size dependence
In this section, we explore the effect of NP size on the steady-state population.Fig. 2 shows the rate of carrier generation in Au NPs, obtained using the method from 35 , as well as the electron and hole population inside the NP after relaxation, at the LSPR frequency ℏω = 2.4 eV.For very small nanoparticles (see results for 2 nm diameter), the generation rate exhibits a molecule-like behaviour with many discrete peaks arising from quantum confinement effects.
The steady-state hot-carrier population exhibits a similar peak structure as the generation rate, but the relative height of the peaks is different compared to the generation rate: after Figure 3: Band structure and density of states of Au and Ag, using the tight-binding model from Ref. 37 .The state that can participate in optical transitions at the LSPR frequency lie within the blue rectangle, and the Fermi energy is set to zero.
the carriers are generated, they relax via electron-electron and electron-phonon interactions, losing energy and approaching the Fermi level.For example, electrons are constantly generated at E = 2.2 eV contributing to = 2.2 eV peak in the electron population (see bottom panel of Fig. 2).However, the electrons constantly relax back to the Fermi energy and contribute to the peak near E = 0 eV.A similar picture holds for the holes generated at To quantitatively study the number of high-energy electrons (δN HE e ) and holes (δN HE h ) in the steady-state population, we evaluate δN where E thr = 1.5 eV denotes a threshold energy.The value of the threshold energy is chosen such that the contribution from the high-energy peaks in δN (E) is captured.
Fig. 4 (top left) shows that significantly more energetic holes are available in the steadystate population than energetic electrons, in particular for larger NPs.This asymmetry can Next, we repeat the analysis for Ag NPs.Ag has a qualitatively similar band structure to Au (see bottom panel of Fig. 3), but both the d bands and the unoccupied sp bands are further away from the Fermi level.The LSPR frequency is ℏω = 3.5 eV instead of ℏω = 2.4 eV.Compared to Au, the LSPR region has a smaller overlap with these bands in Ag (see Fig. 3 ), so we expect the high-energy peaks to be sharper in both the generation rate and the steady-state population due to the smaller number of available transitions.This is indeed what is observed in Fig. 5.The pronounced peak in the hot-hole generation at −3 eV persists in the hole population, but the sharp peak at 3.5 eV in the hot-electron generation rate almost disappears after relaxation effects are taken into account.In its place, a relatively uniform and broad distribution of hot electrons appears.
Because of this sensitivity to the availability of compatible interband transitions in the hot-electron sector, there is a strong size dependency in the number of hot electrons δN HE e to NP size, see Fig. 4 (top right).In contrast, the number of high-energy holes is almost constant.Note that an energy threshold E thr = 2.0 eV was used in Eqs. 3 and 4 for Ag NPs.

Electronic temperature
The nanoparticle's temperature is a key ingredient in explaining many chemical phenomena, and the electronic temperature can be estimated from energy dissipation with the heat equation 42 , or by finding the solution to Boltzmann's equation in steady state 10 .Our formalism also provides a way to compute the temperature of the electrons inside the nanoparticle via its connection to the Fermi function.In equilibrium, the electronic temperature T can be determined from the slope of the Fermi-Dirac population at the Fermi level (which we set to zero here), i.e. f ′ (0) = −1/(4k B T ) (with k B being the Boltzmann constant).Upon illumination, the electron population is no longer determined by a Fermi-Dirac distribution, but instead given by N the density of states which varies slowly near the Fermi level, see Fig. 3.
To extract an electronic temperature from N (E), we assume that it behaves like an equilibrium population with an elevated temperature T ′ near the Fermi level, i.e.N (E) ≈ f T ′ (E)D(0) (where we now explicitly write the temperature as a subscript to the Fermi-Dirac function).Near the Fermi level, we thus have The new temperature of the NP is then found to be .
The increase in temperature as function of the illumination intensity is shown in Fig. 4 (bottom).The temperature increase is mostly proportional to the field intensity.In larger NPs (most pronounced for Ag), a crossover to a parabolic dependence can be observed at larger illumination intensities.Larger NPs exhibit a larger increase in temperature compared to smaller NPs as there is a larger portion of electrons and holes close to the Fermi energy (see Figs. 2 and 5).Ag NPs show a much higher temperature increase than Au NPs for the same field strength.This is explained mainly by the observation that the field enhancement (given in Eq. 2.1) in Ag is much larger than in Au.Using the measured dielectric functions from 36 , we estimate α Au = 1.35 for Au and α Ag = 5.02 for Ag at the LSPR frequency.
The squared field enhancement enters δN (E) explaining the large difference in temperature increase in Ag and Au NPs.

Conclusion
We have developed an atomistic approach to calculate the steady-state hot-carrier population in spherical Au and Ag NPs containing hundreds of thousands of atoms.This was achieved by expressing the population in a basis-independent way and using the Kernel Polynomial Method.Relaxation processes are included in this formalism through an energy-dependent relaxation matrix Γ (ε, ε ′ ).We explored different simple models for the relaxation matrix, but found that the qualitative predictions do not depend on this choice.The hot-carrier populations of Ag and Au nanoparticles exhibit similar features, but the number of highly energetic electrons in Ag NPs found to be highly sensitive to the NP size.We explain this observation in terms of Ag's band structure which features an occupied band whose energy relative to the Fermi level is almost exactly the same as Ag LSP energy.Furthermore, the increase of the electron temperature in the NP is estimated and it is found that Ag NPs exhibit a much larger increase as a consequence of strong electric field enhancement inside the NP.The approach developed in this work can be used to study the hot-carrier population in NPs of complex shapes and compositions.This will be the focus of future work.

Figure 1
Figure 1 compares the carrier population in Au and Ag NPs at their LSPR frequencies for the three relaxation matrices under an external electric field |E ext | = 8.7 × 10 −4 V/nm, corresponding to an illumination power of 1mW/µm 2 .Carriers with positive (negative) energy correspond to electrons (holes).The sign of the hole contribution has been inverted for clarity.Transitions from the d band to the sp band generate hot holes and cold electrons in both metals, giving rise to a hole peak at the position of the d bands (−2 eV for Au, −3.5 eV for Ag) and a corresponding electron peak close to the Fermi level.Transitions in which both the initial and the final states derive from the sp band of the bulk material (also known as intraband transitions) are responsible for the electron peak at 2.2 eV in Au and the smaller electron peak 3.5 eV in Ag.The large number of electrons and holes near the

Figure 1 :
Figure1: Steady-state hot-carrier population in spherical Au (left) and Ag (right) NPs with a diameter of 20 nm for three different relaxation matrices (Γ 1 , Γ 2 , Γ 3 ) at their corresponding LSPR (2.4 eV for Au, 3.5 eV for Ag) at room temperature (T = 298K).The Fermi energy is set to zero, and the hole contribution to the hot-carrier population has been flipped in sign for clarity.The inset shows the high energy peak of the hot-electron population.

Figure 2 :
Figure 2: Hot-carrier generation rate (top) and population (bottom) of spherical Au NPs at the LSPR frequency (ℏω = 2.4 eV) for three different diameters at room temperature (T = 298K).

Figure 4 :
Figure 4: Top panels: Number of high-energy electrons and holes in the steady-state population.In Au, high-energy carriers are defined to have an energy larger than 1.5 eV relative to the Fermi level.In Ag, the corresponding energy threshold is 2.0 eV.Bottom panels: Increase in electronic temperature as function of external illumination strength for different nanoparticle diameters.

Figure 5 :
Figure 5: Hot-carrier generation rates (top) and populations (bottom) for three different sizes of Ag NPs at the LSPR frequency (ℏω = 3.5 eV) at room temperature (T = 298K).