Relaxation of electrons in quantum-confined states in Pb/Si(111) thin films from master equation with first-principles-derived rates

Atomically thin films of Pb on Si(111) provide an experimentally tunable system comprising a highly structured electronic density of states. The lifetime of excited electrons in these states is limited by both electron-electron (e-e) and electron-phonon (e-ph) scattering. We employ the description by a master equation for the electronic occupation numbers to analyze the relative importance of both scattering mechanisms. The electronic and phononic band structures, as well as the matrix elements for electron-phonon coupling within deformation potential theory were obtained from density functional calculations, thus taking into account quantum confinement effects. For the relaxation dynamics, the contribution of impact ionization processes to the lifetime is estimated from the imaginary part of the electronic self-energy calculated in the GW approximation. By numerically solving rate equations for the occupations of the Pb-derived electronic states coupled to a phononic heat bath, we are able to follow the distribution of the electronic excitation energy to the various modes of Pb lattice vibrations. While e-e scattering is the dominant relaxation mechanism, we demonstrate that the e-ph scattering is highly phonon-mode-specific, with a large contribution from surface phonons. At electron energies of about 0.3 eV above the Fermi surface, a 'phonon bottleneck' characteristic of relaxation in nanostructures with well-separated electronic states is observed. The time scales extracted from the simulations are compared to data from pump-probe experiments using time-resolved two-photon photoemission.


I. INTRODUCTION
The thermalization of hot carriers in metals after optical excitation is accomplished both by the Coulomb scattering among the carriers (electron-electron (e-e) interactions) and by the scattering of electrons and holes by lattice vibrations (electron-phonon (e-ph) interaction). In a well-established picture 1,2 , the relaxation can be understood as a two-step process: at early times (t < 0.3 ps), e-e scattering dominates and brings the electrons to a thermal (or possibly non-thermal) distribution. At later times (> 0.3 ps) the e-ph interaction establishes equilibrium between the electronic distribution and the lattice temperature. In this second stage, the high density of excited carriers close to the Fermi energy (within an energy interval corresponding to few phonon quanta) is thought to be responsible for most of the energy flow between the electronic and the phononic system. If so, the e-ph coupling inferred from thermalization experiments should relate directly to the microscopic e-ph coupling constant that governs electric resistivity or the superconducting transition temperature 3 . In this prevailing view, the role of e-ph interactions already in the early stages of relaxation is usually ignored. However, this simple picture is questioned by studies, both experimental and theoretical 4,5 , suggesting overlapping timescales of e-e and e-ph-driven thermalization. Moreover, there is little knowledge how the electrons far above the Fermi level (several tenth of eV) interact with the phonons. For instance, long-lived population of such states, e.g. at the Pb-covered Si(111) surface, has been observed in photoemission experiments 6 . The situation at high energies is in contrast to the e-ph interaction in close vicinity to the Fermi surface, which is crucial for a variety of physical phenomena such as electrical resistivity or supercon-ductivity induced by electron-phonon coupling in thin films 7 , and is quite well explored utilizing the concept of the Eliashberg function (for an overview, see Ref. 8).
In conclusion, there is a need for more studies of the eph interaction at energies further away from the Fermi energy. In this paper, we attempt to obtain a better understanding of the relative importance of e-e and e-ph interaction in highly excited states of a metal and their respective contributions to the early stage of relaxation. To introduce our approach, we have chosen thin multilayer Pb films on Si(111). The fact that this materials system shows a highly structured electronic density of states due to confinement effects 9,10 has been a great advantage for analyzing the energy-dependent lifetime of the excited electrons using time-resolved pumpprobe spectroscopy. 11 The experimental results were rationalized in Ref. 11 in terms of e-e interaction only, and it was concluded that the electronic lifetime closely follows the behavior expected from Landau's theory of Fermi liquids. 12, 13 Yet, a contribution of e-ph scattering to the lifetime cannot be excluded completely based on the achieved level of agreement between experiment and theory. Therefore, we aim at a detailled analysis of the role of e-ph scattering for the features observed in photoemission. Since ample experimental and computational data are available for the Pb/Si(111) films, we consider this system a good test case for quantitative studies of electronic relaxation dynamics. In a previous paper 14 by us, we have worked out a realistic atomistic description for multilayer Pb films on Si(111) and have carried out first-principles calculations of the electronic and phononic band structure and of e-ph coupling in electronic states far away from the Fermi level. While the e-ph interaction in bulk solids has become accessible to first-principles calculations by using density functional perturbation theory together with Wannier interpolation methods to enhance the number of reciprocal-space sampling points [15][16][17][18] , thin films on a substrate are still difficult to treat on a microscopic level because the adequate supercell typically contains tens to hundreds of atoms and computational costs are high. For the Pb films on Si(111), for instance, the complex phase diagram 19,20 results in various reconstructions requiring large supercells for their description 21 . In the present work, we constructed a √ 3 × √ 3 unit cell of Si(111) matched with a (2 × 2) unit cell of Pb(111) to describe the atomic structure consisting of 40 Pb and 30 Si atoms. 14 While the two-dimensional Brillouin zone of reconstructed surface plus interface is smaller than the Brillouin zone of a bulk material, the supercell contains a large number of bands, both in the electronic and phononic spectra. Therefore, a thoughtful selection of bands will be required to arrive at a tractable model for e-ph coupling. The approach via density functional perturbation theory and the calculation of the Eliashberg function would be too cumbersome for large supercells.
In this paper, building upon the knowledge of our previous work 14 , we elaborate on the consequences of these microscopic data for the e-ph scattering rate using a kinetic master equation. The detailed modeling of e-ph scattering is combined with a description of the e-e interaction at the level of Fermi liquid theory. This combination allows us to simulate the temporal evolution of electronic populations on the relevant scales and to make contact with experimental observations.

II. THEORY
The general problem of an excited electronic system coupled to lattice degrees of freedom can be approached from various perspectives. If one is satisfied with a classical description of the atomic positions and velocities and their dynamics can be described in the trajectory approximation, carrying out non-adiabatic molecular dynamics simulations (e.g. with the methodology described in 22,23 ) is the method of choice. As an advantage, this approach can handle large deviations of the atomic positions from their ground state, and the forces acting on the atoms are calculated directly within the first-principles electronic structure framework. Thus, it is suitable for systems with very strong and non-linear electron-phonon coupling, as encountered e.g. in two-dimensional materials 24 . In this work, we emphasize the quantum nature of the phonons, while the weak coupling of the electrons to phonons and to external fields can be treated in first-order perturbation theory. Casting the problem into the form of a model Hamiltonian, it reads H = H 0 + H int with H 0 being the ground-state Hamiltonian with phonons described in the harmonic approximation, The creation and annihilation operators c nk , c † nk and b IQ , b † IQ obey the usual anticommutator relations for fermions and commutator relations for bosons, respectively. The first, integer index n specifies the band, while the second index k describes the crystal momentum in the form of a two-dimensional vector within the Brillouin zone of a thin slab. Capital letters are used to index phonon modes, whereas small letters refer to electronic bands. In contrast to the molecular dynamics approach mentioned at the beginning of this paragraph, the full quantum treatment is best suited when the coupling terms in the interaction Hamiltonian H int are weak, and the model Hamiltonian H 0 provides already a good starting point for the coupled dynamics. Utilizing model Hamiltonians for describing electronic dynamics is a well-established technique in the field of ultrafast soild-state optics, see e.g. 25,26 . In semiconductor bulk materials and quantum wells, the dispersion ε nk entering the Hamiltonian can be approximated as being quadratic (and sometimes as being linear, e.g. for graphene 27 ), and a full solution of the relaxation dynamics for various scattering mechanisms has been achieved in these cases. Here, we are interested in a realistic description of the ground state of a particular system. For this reason, all the band energies and phonon frequencies entering H 0 are determined by density-functional theory calculations. The VASP code 28 with the settings described in Ref. 14 has been employed for this purpose. The electronic single-particle energies ε nk are taken to be equal to the Kohn-Sham eigenvalues obtained with the GGA-PBE exchange-correlation functional 29 . The phonon frequencies Ω IQ and the corresponding eigenmodes are obtained from DFT calculations using the method of finite atomic displacements within a supercell, as detailed in Ref. 14 and 30. In case of the Pb/Si(111) films, such a detailed first-principles description is considered necessary in view of the experimental findings: The two-photon photoemission spectra show peaks at certain intermediate-state energies of the electrons that are referred to as quantum well states (QWS). These are energies where the electronic density of states is high and/or where the excited electrons are long-lived. For a correct prediction of the energetic position of the QWS, the (1 × 1) periodicity of a free-standing Pb(111) films is not sufficiently accurate. 31 . It is required to take the larger ( √ 3 × √ 3) periodicity enforced by the Si(111) substrate into account. As major achievement of the first-principles calculations in Ref. 14, we were able to reproduce quantitatively the dependence of the energetic position of the QWS on the number of Pb layers in the film, as well as the very small dispersion of the occupied QWS in the films with an odd number of Pb layers. On this basis, the present work is addressing the role of the elec-tronic lifetime in the QWS for the experimentally detected peaks. The interaction Hamiltonian H int contains any further interactions required to describe the problem at hand. These interactions could e.g. be the electron-electron interactions beyond the effective mean-field description of density functional theory (see below). Moreover, the interaction with an external electromagnetic field, e.g. of a laser pulse, can be considered as part of H int . Most importantly for the present study, H int contains a term H ep describing in linear order the coupling of the electrons to quantized phonons, The term in parentheses is linear in each phonon coordinate.
In principle it is possible to describe the quantum nonequilibrium dynamics under the action of H exactly by a density matrix. Schemes for evolving the density matrix in time have been worked out 32 , and applications to surfaces and low-dimensional systems can be found in the literature. 27,[33][34][35] However, since the system we want to describe is quite complex, we resort to a simpler description of the dynamics which is appropriate if the coherent excitation by an optical pulse and the subsequent relaxation take place on separable time scales. While quantum coherence is important during the interaction of the system with the light field, electron-electron scattering usually leads to a fast loss of coherence. 36 For Pb films, an example of calculations taking the effects of coherence into account can be found in 6 . In the limit of vanishing coherence, only the diagonal elements of the density matrix, the populations f nk of states indexed by n and the wave vector k, are important. For the investigation of the ultrafast population dynamics in our system, the quantities which we have to look at are the electronic occupation numbers f nk = ĉ † nkĉ nk and the phononic occupation numbers n IQ = b † IQb IQ . For the latter, we employ a bath approximation In the numerical calculations presented below, we will use different baths, one for each high-lying optical mode of the Pb film (Ω IQ ≥ 2 THz) with temperature T I , and a common one for all low-frequency phonons of the Pb film (Ω IQ < 2 THz) with temperature T 0 . More details are given in the appendix. Using the Markov approximation and the second-order Born approximation for the transitions, it is possible to derive from the density-matrix equations a set of coupled differential equations that can be cast into the form of a master equation (cf. Ref. 37).
The expressions for the rates, both for scattering into and out of the state nk, are made up of an electronic and a phononic contribution each, i.e., Γ nk = Γ (ee) nk . This holds for both Γ in and Γ out that both consist two terms owing to electron-electron scattering and electronphonon scattering: Exploiting conservation of crystal momentum parallel to the film, k−k = ±Q with the sign depending on phonon emission or absorption, the electron-phonon scattering rates originating from the Hamiltonian (2) can be expressed according to Fermi's golden rule as These expressions include processes where the electron absorbs a phonon as well as phonon emission processes. This is denoted by the ± signs in the equations, where the minus sign stands for absorption, and the plus sign for both spontaneous and induced emission, proportional to n IQ + 1. It is our goal to calculate the contribution of e-ph scattering to the lifetime of specific quantum well states (QWS) in Pb/Si(111) films. In a very simple picture, the conduction band electrons of Pb with crystal momentum normal to the surface or interface of the Pb(111) films are confined, similar to the quantum-mechanical particle-ina-box problem. In an atomistic picture, these conduction band states are derived from the 6p z orbitals of the Pb atoms and their wavefunctions extend both above the surface and into the Si(111) substrate, see Ref. 14 . As described in the experimental paper 11 , there are significant differences between the lifetimes in films with an even and an odd number of Pb layers. Therefore, we study two representative systems, a Pb film with 4 monolayers (ML) and one with 5 ML on Si(111). Side views of the corresponding slabs are depicted in Fig. 1. Motivated by the experimental focus on excited electrons in unoccupied bands, we include e-ph scattering rates for the electrons excited into QWS. Since the population of the valence bands was not analyzed in these experiments, the hole states are treated in less detail, and and only Coulomb scattering, as described in Section III B, will be considered among the holes. To solve the rate equations, we need explicit expressions for the quantities D mk−Q nk, IQ and n IQ (T ) in eq. (6) entering the decay rates Γ in, (ep) and Γ out, (ep) . Both quantities depend on the phonon . The bands with large Pb 6pz character are highlighted by the red symbols and thick red lines. The dashed ellipses mark the regions of quantum wells states (QWS) whose lifetime under e-ph scattering is shown in Fig. 2(a). For 5 ML Pb, the occupied quantum well resonance at ∼ −0.25 eV has been included in the plot. The electronic eigenvalues are represented on a 32 × 32 Monkhorst-Pack grid. The plots are shown along a diagonal cut K' -Γ -K through the Brillouin zone of the ( branches Ω IQ . Of all phonon modes Ω IQ of the supercell obtained with our first-principles approach 30 , those with Pb character are taken into account, see Fig To keep the number of individual scattering processes at a tractable level, we also restrict ourselves to a subspace of the electronic bands: Since we are interested in the electron-phonon coupling in QWS in Pb, only those electronic bands that have a significant overlap with the Pb 6p z orbitals, as indicated by the VASP calculation, are retained in the Hamiltonian H ep in eq. (2). The electronic states belonging to a specific Pb-derived band are grouped together into subsets indexed by α(k) ∈ {n} of all band indices n. To be specific, we used the five (six) lowest-lying conduction bands with appreciable Pb 6p z character for the 4 ML and 5 ML Pb film, respectively, i.e. α = 1, . . . 5, (6). These 5 (6) bands are displayed in Fig. 1 by the thick red lines and symbols, together with the full band structure (dashed lines) that is also shown (over a wider range of energies and wavevectors) in Fig. 3 and 4 of Ref.
14. Due to the use of a supercell and backfolding of the bands, these bandstructures are different from the bandstructure of Pb(111)(1 × 1) slabs that had been used previously 11 in the experimental data analysis.
For evaluating the electron-phonon scattering rates, eq. (6), we use techniques based on deformation potential theory that allows us to obtain Γ out (ep) nk from firstprinciples calculations of the phonon spectrum, the electronic wavefunctions, and Kohn-Sham eigenvalues, with only few approximations. As the most significant one, we neglect of the Q-dependence of the deformation potential, while keeping its dependence on band index n and crystal momentum k. This is a good approximation for optical phonons and corresponds to keeping the leading (constant) term in an expansion in powers of Q, cf. Ref. 38. In the energy-conserving δ-function in eq. (6), we retain the finite phonon energy Ω IQ ≈ Ω I0 , but neglect the dispersion of the optical phonon branches. This is justified since the dispersion remains small (cf. Fig. 6 in Ref. 14) due to the large real-space unit cell, and hence small Brillouin zone, of the Pb films. Within these approximations (see appendix for the derivation), the matrix element for electron-phonon scattering in eq. (6) can be replaced by 3) supercell used to model the Pb/Si(111) film, M Pb and v Pb atom are the atomic mass and atomic volume of Pb. D nk, I is the deformation potential of the n th electronic band under the phonon mode I. The D nk, I have been obtained from DFT calculations 14 by evaluating the electronic eigenvalue shift under finite displacements of the atomic positions given by the corresponding mode eigenvector of the phonon. The two δ-symbols reflect conservation of crystal momentum in e-ph scattering, and the projection of H ep to the finite electronic subspace, as described above. The matrix elements I mk nk account for the difference between intra-band (n = m) and interband scattering (n = m), and for the dependence on both the initial and final electron momenta k and k = k − Q. They are obtained from the overlap of the corresponding DFT wave functions. More details of the derivation are given in the appendix. In summary, this approach allows us to arrive at a simplified and computationally tractable, yet parameter-free description of e-ph scattering even for such a complex systems as an overlayer on a substrate.

III. RESULTS
In metals, e-ph and e-e scattering are closely intertwined, since the vast majority of phonons is emitted by sec-ondary electrons and holes rather than by the charge carriers initially excited by the light pulse. This is because ee scattering quickly generates an avalanche of secondary electron-hole pairs with small energies around the Fermi level. Since these secondary electrons and holes are produced with high density and their energy still exceeds typical phonon energies, they play a major role in determining the rate at which the energy is dissipated from the electronic system into the lattice. Nevertheless, we start our discussion by considering the contribution of both e-ph and e-e scattering separately.
A. Relaxation due to e-ph scattering First we investigate how the population of a QWS decays under the sole effect of e-ph scattering. For this purpose, we initially populate a single QWS at the Γ-point and let the population evolve according to the master equation (3) using only the rates Γ in, (ep) and Γ out, (ep) . The results are shown in Fig. 2(a). At comparable energies of the QWS of ∼ 0.5eV, the decay is much faster in the 4 ML than in the 5 ML Pb film. This is to be expected from the different size of the deformation potentials in the two films reported in Ref. 14 . The relaxation rate increases with the temperature of the phonon heat bath, which is indicative of the role of stimulated emission of phonons. By decreasing the phonon temperature from 400 K and 100 K, the lifetime of the QWS in the 4 ML film increases from 1.3 to 2.7 ps. For the 5 ML film, the lifetimes fall between 13 and 37 ps.

B. Relaxation due to e-e scattering
The lifetime of hot electrons due to e-e scattering can be described by a self-energy formalism, as discussed in Ref. 39. The loss term Γ out (ee) nk in eq. (4) is given by Γ out, (ee) nk = −2ImΣ(ε nk )/ . The self-energy Σ is obtained from a GW calculation of bulk Pb. Here, G stands for the electronic Green function, and W for the screened Coulomb interaction. These quantities are calculated from the DFT wave functions and Kohn-Sham eigenvalues using the built-in capabilities of VASP 40 . To be specific, a 11 × 11 × 11 k-point mesh is used, and the denominator in G is evaluated with a small shift of the transition energy away from the real axis, η = 0.08 eV, much smaller than typical values used in GW calculations of band structures. The result obtained for −2ImΣ in the conduction band is fitted to the α(ε nk − E F ) 2 dependence expected from Landau's theory of the Fermi liquid. Our result α = 0.022 (eV) −1 is in excellent agreement with earlier GW calculations of bulk Pb 41 . Finally, we obtain the expression Γ out, (ee) nk plotted in Fig. 2(b).
In the relaxation of highly excited electrons and holes, the energy is dissipated to secondary electron-hole pairs. This process is very efficient in metals, since, in contrast to semiconductors, there is no energy gap preventing the generation of secondary particles. These effects are included in the the scattering-in term Γ in (ee) nk of eq. (5). Although this term can be obtained from the master equation 39 as well, we choose for computational convenience a simpler treatment in our present study. The gain term is assumed to factorize into an energy-dependent and a time-dependent factor, Γ in, (ee) nk = Φ(x)N (t).The distribution function Φ describes the secondary electrons and holes produced via impact ionization by a relaxing high-energy electron. Following the work of Baranov and Kabanov 5 , we use for Φ a stationary solution of the Boltzmann equation with a Coulomb scattering kernel, where the electronic temperature T el = 650 K was chosen in accordance with the energy of E dep = 0.1eV = π 2 6 g(E F )(k B T el ) 2 deposited by the laser and the electronic heat capacity of Pb The time-dependent factor N (t) for creation of secondary electrons and holes is determined by energy conservation in the e-e scattering. Our simplified treatment assumes that the initial electron in state nk ends up at the Fermi energy, transferring all its initial energy to secondary electron-hole pairs. This motivates the choice where g(ε) is the electronic density of states. Both g(ε) and N (t) are evaluated numerically using as input the DFT band structure of slab models for Pb/Si(111)( √ 3 × √ 3) multilayer films.

C. Competition between e-e and e-ph scattering
In this Section, we compare simulation results for Pb films of 4 ML and 5 ML thickness as representatives of films with an even and odd number of layers studied by optical pump-probe experiments in Ref. 11. While these experiments measure the total probability for two-photon photoemission, our simulations model the population of the intermediate electronic states that are reached by the electron after applying the pump pulse and subsequent relaxation. The second step of the two-photon photoemission, which kicks the electron into the vacuum, is not modeled. Provided that the probability of ionization by the probe pulse is a smooth function of energy, the measured yield can be considered approximately proportional to the population of the intermediate state.
The initial distribution is chosen such that it describes the response of our specific system, multilayers of Pb on Si(111), to a short optical pulse with frequency centered around hν = 1.9 eV. This corresponds to the photon energy of the pump laser used in the experiment 11 . The polarization of the electric field, denoted by the unit vector e, is chosen parallel to the Pb film surface. Before the laser pulse arrives, the system is described by a Fermi-Dirac distribution with low temperature, T el → 0; hence f nk (t < 0) = Θ(E F − ε nk ). To be specific, we evaluate dipole matrix elements 42 In the numerical evaluation 43 , a broadening of the δfunction by 0.02 eV is used. The proportionality factor A 0 is chosen such that the energy of excited electrons and holes deposited in the Pb films amounts to ∼ 0.1 eV per supercell area, equivalent to 3.7 µJ/cm 2 . By solving the master eq. (3) numerically, we are able to follow the relaxation of the excited electrons in real time.
We define an energy and time dependent population of the intermediate state Fig. 3 shows on a logarithmic scale the energy distribution P (ε, t j ) of the excited electrons for various times t j after the excitation. For plotting the results, the δfunction in eq. (11) has been replaced by a rectangle with a width of 0.06 eV. We start with a discussion of the initial distribution, shown by the thick black line, calculated according to the transition dipole strength, eq. (10). For the 4 ML Pb film ( Fig. 3(a)), the distribution is highly structured with a sharp maximum at 0.58 eV and a broad peak around 1.21 eV. For the 5 ML Pb film (thick black line in Fig. 3(b)), only the peak at 1.21 eV (and possibly a short-lived peak at higher energies) remain visible, while the low-energy peak is much less pronounced. These results are in excellent agreement with the experimental observations of Ref. 11. In this work, a high-energy peak in the range of 1.1 to 1.2 eV was observed for films with an odd number of Pb monolayers, whereas the peak at 0.6 eV was dominant in Pb films with an even number of layers. Note that, due to experimental limitations of the probe laser energy, excited electrons with energies lower than ∼ 0.5 eV could not be detected in Ref. 11. Next, we analyze the relaxation of the energy distributions for later times. From Fig. 3 it is obvious that all distributions develop a low-energy part corresponding a quasi-thermal distribution of secondary electrons, showing up as an exponentially decreasing function of energy. The high-energy part of the initial spectrum decays mainly due to e-e scattering, thereby creating secondary electrons via impact ionization. Therefore, the highenergy tails decay quickly, simultaneously accompanied by an increasing weight of the secondary-electron distribution. Now turning to longer time scales, we observe that the low-energy part in the 4 ML Pb film and the population in the energy window between 0.4 and 0.6 eV decay more slowly. At the same time, a broad shoulder on top of the secondary electron distribution builds up at ∼ 0.3 eV (marked by the arrow in Fig. 3(a) and magnified in the inset). Both phenomena can be traced back to the effect of e-ph scattering: Excited electrons of the initial distribution at energies of 0.4 to 0.6 eV 'glide down' the electronic band structure (cf. Fig. 1(a) ), thereby emitting phonons. The relatively slow rate at which this occurs leads to a 'phonon bottleneck', i.e., to the buildup of the shoulder centered at ∼ 0.3 eV. This effect, related to the discreteness of electronic states, is quite common for the e-ph relaxation in nanostructures, and has been observed e.g. in quantum dots 44 as well as in two-dimensional layered semiconductors 45 . While there is still continuous energy dissipation to the crystal lattice by very low-energy secondary electrons, the 'phonon bottleneck' affects electrons at higher energies and results in additional, but delayed production of phonons emitted by these electrons 'gliding down' the conduction bands. A similar, but weaker 'phonon bottleneck' can be observed on top of the secondary electron distribution in the simulations for the 5 ML Pb film in Fig. 3(b), marked by the arrow and magnified in the inset. In Fig. 4, we analyze the time scales associated with the electronic population decay at selected energies where peaks had been found in the populations in Fig. 3. The lines without symbols in Fig. 4 show the decay according to the full relaxation dynamics, including both e-e and e-ph scattering. All populations show a nearly exponential decay, albeit with different decay rates. The lifetimes for the population maxima at 1.21 eV and 0.58 eV, extracted via exponential fits, are 21 fs and 101 fs, respectively. For the lowest energy of 0.46 eV, we find an initial rise of the population due to scattering-in from electrons at higher energies, followed by a population decay after about 30 fs, corresponding to a lifetime of 183 fs. The rather broad (in time) maximum of the 0.46 eV curve results from a compensation of the rates of incoming electrons from higher energies, mostly originating from e-e scattering at these high energies, and losses due to both e-e and e-ph scattering, the latter one gaining in relative importance as we go to lower energies. Interestingly, the described broad maximum of the population evolution at low energies is also seen in the experimental data 11 . In the computer simulation, we can deliberately turn off the e-e scattering channel after a very short initial time interval. The e-e scattering was permitted only in the very early times, somewhat arbitrarily chosen to be less than 6 fs, since some mechanism is required to establish a realistic smoothened electron distribution including an appropriate low-energy secondary-electron part. The result of these runs are displayed in Fig. 4 by the lines with circular symbols. If the relaxation after 6 fs proceeds by e-ph scattering only, the population at 0.46 eV initially decays on a time scale of 350 fs (green symbols in Fig. 4(a)), which can be taken as an estimate of the e-ph scattering rate in this energy range. This initial decay is followed by a much slower decay over several picoseconds. At the higher electron energies of 0.58, 0.81 and 1.21 eV, the scattering-in events of electrons from higher energies are equally probable or even more frequent than the scattering-out events, and hence a net contribution of e-ph scattering to the decay is not detectable on the time scale shown in Fig. 4. A similar analysis has been carried out for the 5 ML Pb film, see Fig. 4(b). For the peak energies at 1.21, 0.81 and 0.52 eV, overall lifetimes of 21, 47 and 126 fs are obtained from exponential fits to the full relaxation dynamics. Again, it is possible to estimate the relative importance of e-ph scattering by watching the population decay after the e-e scattering has been 'turned off'. From the slopes of the curves marked by the circular symbols in Fig. 4(b), characteristic times of 24 ps and 4.1 ps are obtained for 0.81 and 0.52 eV electron energy, respectively. At the highest electron energy of 1.21 eV, again we find that the contribution of e-ph scattering is too small to be detectable on the time scale shown in Fig. 4. From this analysis, we learn that the contribution of eph scattering to the total lifetime of the peaks at energies larger than 0.5 eV is much smaller compared to the e-e contribution. This finding confirms the original analysis of the experimental data by Kirchmann et al. 11 where e-ph scattering had been disregarded. Analyzing the experimental data for many Pb film thicknesses, they concluded that the low-energy peak is clearly observed in films with an even number of atomic layers and has a lifetime of 115 ± 10 fs, while the high-energy peak is visible only in the odd-layer films and has an energydependent lifetime which turns out to be 10 ± 5 fs for 5 ML Pb. Our simulation results of 101 fs and 21 fs are in reasonable agreement with their experimental findings, in particular if it is taken into account that the e-e scattering rate is very sensitive to the precise energetic position of the peak. The lifetimes extracted from our simulations are summarized by the circular symbols in Fig. 2(b). Despite the additional decay channel of e-ph scattering being taken into account, the simulated lifetimes lie above the lifetime of isolated electrons due to e-e scattering alone. This is because the simulations describe a realistic distribution of excited electrons, and the incoming flux from higher-lying electronic states effectively 'conserves' the population of the lower lying states over longer times.

D. Excitation of lattice vibrations
Although the contribution of e-ph scattering to the lifetime of the quantum well states was found to be small, the low-energy states populated by the secondary electrons couple significantly to the lattice vibrations. At these low energies, the e-ph scattering as loss mechanism even dominates over e-e scattering, since the lifetime due to e-e interactions rises above 300 fs for electrons below 0.33 eV according to eq. (9). With the help of the simulations, it is possible to follow the energy transferred from the electrons to each of the phonon modes separately. Fig. 5 shows the increase in time of the excess vibrational energy (in addition to the thermal energy corresponding to the initial substrate temperature) in the various Pb vibrational modes. Summation over all modes yields a total energy transfer between the electronic and the lattice degrees of freedom of 8.2 meV/ps for the 4 ML Pb and 0.79 meV/ps the 5 ML Pb film, respectively. The smallness of these quantities (a few meV compared to 0.1 eV electronic energy in the film) gives an a posteriori justification for the perturbative expression used for the interaction Hamiltonian H int . It is seen from Fig. 5 that the energy transfer is highly mode-selective 46 . One particular surface phonon mode receives a major part of the energy. The dominance of this single mode depends on film thickness; in the 4 ML film it is clearly more pronounced than in the 5 ML film where several modes participate in the energy uptake. In the 4 ML Pb film, this is a mode with frequency 2.26 THz (labeled 47 in our previous publication 14 ), while in the 5 ML Pb film it is a phonon mode at 2.03 THz. In experiments using a high laser fluence, the strong coupling to a specific mode has indeed been observed; in the 5 ML Pb film, it was detected experimentally by a periodic shift of the quantum well energy with a frequency of 2.0 ± 0.1 THz 47 , which is in excellent agreement with the frequency of 2.03 THz identified in our simulation. For both film thicknesses, the modes at frequencies below 2 THz, which are similar to phonon modes in bulk Pb, receive only a much smaller amount of energy on average. As seen from Fig. 5, the increase of the energy over long time scales, e.g. over 4 ps, is sub-linear and saturation is expected at even longer times. We attribute this slow response to the rather long time required by the electrons to relax from high energies down to E F via multiple phonon emission processes. Due to the 'phonon bottleneck' described above, this takes longer than expected from Allan's formula 3 . This formula requires as sole input the electron-phonon coupling constant of the bulk material that can e.g. be determined from ultrafast reflectivity measurements 48 . Eventually, we comment on the possibility to directly observe the excitation of the lattice experimentally. The excitation of low-lying Pb modes results in irregular displacements of the atomic positions on medium to large length scale that can be followed in the Debye-Waller factor of a time-resolved electron diffraction experiment. Preliminary experimental data for Pb films 49 indicate that the low-energy phonons are indeed getting excited on the time scale of a few picoseconds. Interestingly, recent diffraction experiments were able demonstrate the mode-selective energy transfer to phonons even for bulk materials such as aluminium 50 or nickel 51 . Apparently, phonon excitation by strongly excited charge carriers opens up the exploration of a new class of nonequilibrium phenomena, not only in bulk metals, but also in semiconductors 52 and nanostructures. Due to the widely different timescales of e-ph interaction and phonon-phonon interaction, the non-equilibrium distribution of phonons created by the hot carriers can persist over a time span of several picoseconds. Thus, non-thermal phonon distributions could be a more widespread phenomenon than previously thought, and further research along these lines could be fruitful.

IV. CONCLUSION
Simulations of electronic relaxation have been performed for a realistic system, metallic multilayer Pb films on Si(111), with the help of a parameter-free approach based on density functional theory. Electron-electron (e-e) and electron-phonon (e-ph) scattering were both included in the master equation. Not surprisingly, e-e scattering was found to dominate over e-ph scattering for short times and highly excited electrons more than 0.5 eV above the   ing the fastest e-ph relaxation found in the present simulations with a characteristic time of 350 fs as a marker for the cross-over between e-e and e-ph scattering, we conclude that e-ph scattering significantly contributes to the relaxation of electrons in Pb at energies below 0.3 eV. Indeed, the simulations show a population pile-up around 0.3 eV due to the cross-over of the e-e and e-ph scattering time scales at this energy. After 300 fs, it is fair to describe the electronic population by a thermal distribution, however, only up to an excess population of electrons getting stuck in the 'phonon bottleneck'. Remark-ably, this does by no means imply that the phonon populations could be described by a temperature as well. Contrarily, the simulations show that even up to 4 ps after excitation high-frequency surface vibrational modes are preferentially excited by strongly mode-selective phonon emission.
In summary, our simulations enable us to disentangle the contributions of e-e and e-ph scattering at short times, < 0.3 ps after optical excitation. Although a first glance at the data is compatible with electronic thermalization by e-e scattering, a contribution of e-ph scattering can be observed already in this early stage, in particular at low electron energies. The phonon system requires much longer (several picoseconds) to equilibrate. Additional simulations beyond the scope of the present work are desirable to gain an improved understanding of the energy transfer between the two subsystems in the later stages of relaxation. The derivation of the form of the electron-phonon matrix element, eq. (8), starts from the Hamiltonian H ep brought to its real-space representation, The physical interpretation of ∆V I Q ( r) is the additional perturbing potential in the Kohn-Sham equation due to the presence of a phonon of mode I and wave vector Q. The vector notation refers to vectorial quantities in threedimensional (3D) space. In a very general setting, Fermi's golden rule requires us to calculate matrix elements of the form involving 3D real-space integration. We now switch to the situation of interest, a two-dimensional film with periodic boundary conditions in the x, y coordinates. Within deformation potential theory 38 , the integration over the coordinates x and y in the film plane are carried out to yield a (spatially independent) shift of the Kohn-Sham eigenvalue, denoted by D nk, I Q , multiplied by a δ-symbol of parallel momentum conservation, The bold symbols denote two-dimensional vectors in the film plane; the '+' sign refers to absorption, the '−' to emission of a phonon of wavevector Q. Technically speaking, the deformation potentials are obtained in the following way: A phonon eigenvector labeled I, as obtained from the phonopy code, is scaled with √ M Pb and then added or subtracted from the Cartesian positions of the Pb atoms in the unit cell. For the two geometries obtained in this way (called + and −), static DFT calculations are carried out, yielding Kohn-Sham energy eigenvalues ε + nk I and ε − nk,I . The eigenvalue shift determines the deformation potential via Note that the (optical) deformation potential D nk, I0 has by definition the physical unit of energy/length. Now we return to the evaluation of the matrix element in eq. (A2). We work out the equations for general wavevectors k and k and will enforce momentum conservation in a later step. After squaring the matrix element to obtain the expression for the rate constant, we can write Here, the constant eigenvalue shift has been taken out of the brackets, and φ nk (z) denotes the Bloch-periodic part of the full 3D wavefunction Ψ nk ( r) suitably averaged over x and y. This expression requires to evaluate integral(s) over the remaining z coordinate that denotes the spatial direction normal to the film in which the electrons are confined, where the definition has been used. To obtain the overall scattering rate due to one specific phonon mode I, the summation over Q in eq. (A1) needs to be carried out, or equivalently an integration over the Brillouin zone of the phonon. Therefore the expression for the rate will include Qx,Qy If we now assume that the deformation potential D is only weakly dependent on Q (which is reasonable for optical phonons in a large supercell), we can take it out of the sums and integral, set Q = 0 and D nk,I0 =: D nk,I , and perform the integration over Q z prior to the integration over z and z . Motivated by this procedure, we define a quantity (cf. the analogous case for a 1D wire, rather than a 2D film, in the appendix of Ref. 53) As a consequence of the definition of J n k n k in eq. (A8), it follows that the Q z integration yields δ(z − z ), and the two factors with arguments z and z in the expression (A7) turn out to be complex conjugates of each other. Performing the integral over z together with δ(z − z ) leaves us with a single integral over z, but now with the squared moduli of the wavefunction in the integrand. In the numerical implementation, we not only integrate over z, but simultaneously perform the spatial average over the supercell by integrating over x and y and dividing the result by the area of the supercell A supercell : For this expression to be valid, the squared wavefunctions must be normalized such that in other words, Ψ must have the physical unit length −1/2 . Finally, by inserting the expression (A10) for the double integration in eq. (A7) and properly normalizing to A supercell , we arrive at the result given in eq. (8).
In practice, the squared wavefunctions |Ψ nk (x, y, z)| 2 are obtained from the DFT calculations, using the PARCHG keyword of the VASP code. To reduce the sheer number of calculations, we calculate I n k n k explicitly for all combination of band indices n and n of the bands in the energy range [E F , E F + 2eV], and for the three combinations of the two k-vectors, k x = (1/4, 0, 0) and k y = (0, 1/4, 0). For other combinations of k-vectors, we interpolate using the angle θ between the wave vectors as variable, defined via the scalar product cos θ = kk kk , and obtain Already from inspecting the numerical values of the deformation potential D nk, I , it becomes clear that the coupling of an excited electron to a phonon can be highly phonon-mode specific. In particular, coupling to the highest-lying Pb modes is strong. We therefore use an independent heat bath for each the N opt = 6 highest-lying Pb modes. For the 4ML Pb/Si(111) film, these six phonons (surface and interface modes) of the Pb layer are in the 2.1 to 2.5 THz range (modes labelled 44,45,46,47,48 and 49 in Ref. 14). For the 5ML Pb/Si(111) film, the analogous modes are found at somewhat lower frequencies, in the 2.0 to 2.3 THz range. All the other, lower-lying modes are taken to constitute a common acoustic phonon bath with temperature T 0 . In addition, a constant lattice temperature T sub of the Si substrate acting as a heat sink is part of the description of the vibrational system. As initial condition, all temperatures are set to a base temperature of 100 K at the beginning of the simulation. In general, the population n I (t) of a phonon mode I varies slowly on the scale of electronic relaxation. Yet we allow for energy exchange between the various heat baths. Previous simulations 54 using classical molecular dynamics showed that the surface phonon in a Pb/Si(111) monolayer is damped due to mode conversion on a time scale of τ conv = 30 ps. In the present simulation, we use this time constant to couple each of the baths at temperature T I (t) to the common acoustic bath at T 0 (t). Moreover, the acoustic bath may transfer energy to the Si substrate lattice on an even longer time scale τ sub = 160 ps. This value was adopted from measurements observing the equilibration of Pb films with the Si substrate by Witte et al. 55 . The above considerations lead to the following equations describing the evolution of lattice temperatures: The quantities c acu V and c opt V denote the partial heat capacities of the acoustic phonon bath, and of one optical mode, respectively. In the 'classical' approximation, i.e. when the equipartition theorem holds, we have c opt  by the arc length (in k-space) of the solution, approximated by the cumulated lengths of the line segments, as illustrated in Fig. 6. The right hand side of the Master equation is thus made up by the sum over m, I, all grid points and pertaining line segments, in total ∼ 10 6 terms.