Ultrafast dynamics of electrons and phonons: from the two-temperature model to the time-dependent Boltzmann equation

The advent of pump-probe spectroscopy techniques paved the way to the exploration of the ultrafast dynamics of electrons and phonons in crystalline solids. Following photo-absorption of a pump pulse and the initial electronic thermalization, the dynamics of electronic and vibrational degrees of freedom is dominated by electron-phonon and phonon-phonon scattering processes.The two-temperature model (TTM) and its generalizations -- as, e.g., the non-thermal lattice model (NLM) -- provide valuable tools to describe these phenomena and the ensuing coupled electron-phonon dynamics over timescales ranging between few hundreds of femtoseconds and tens of picoseconds. While more sophisticated theoretical approaches are nowadays available, the conceptual and computational simplicity of the TTM makes it the method of choice to model thermalization processes in pump-probe spectroscopy, and it keeps being widely applied in both experimental and theoretical studies. In the domain of ab-initio methods, the time-dependent Boltzmann equation (TDBE) ameliorates many of the shortcomings of the TTM and it enables a realistic and parameter-free description of ultrafast phenomena with full momentum resolution. After a pedagogical introduction to the TTM and TDBE, in this manuscript we review their application to the description of ultrafast process in solid-state physics and materials science as well as their theoretical foundation.

Pump-probe experiments provide a versatile tool to directly probe these phenomena CONTACT F. Caruso. Email: caruso@physik.uni-kiel.de arXiv:2202.07237v1 [cond-mat.mtrl-sci] 15 Feb 2022 [34,35]. Time-and angle-resolved photoemission spectroscopy, for example, enables to directly probe the dynamics of photoexcited electrons and holes with energy and momentum resolution [36,37]. Additionally, pump-probe scattering techniques as, e.g., ultrafast electron diffuse scattering complement optical and photoemission measurements by providing direct insight into the dynamics of the crystalline lattice and electron-phonon scattering processes with time and momentum resolution, and they are suitable to attain a detailed understanding of the energy flow between hot electrons and the crystalline lattice [38][39][40].
A rigorous description of electron-phonon coupling and its influence on the dynamics inevitably requires to account for phase-space constraints in the phonon-assisted scattering processes as well as for the anisotropies of the electron-phonon coupling. A suitable approach to attain this requirement is the time-dependent Boltzmann equation (TDBE). Ultrafast dynamics simulations based on the TDBE formalism have been widely employed to described non-equilibrium processes involving electrons and phonons in condensed matter. The TDBE has been applied to investigate the thermalization of electrons and phonons following photoirradiation in metals [126][127][128][129][130][131][132][133][134], semiconductors [135][136][137][138], and two-dimensional [139][140][141][142][143][144] and layered [145] materials, and it has further been extended to investigate the ultrafast magnetization dynamics in photo-excited ferromagnets [146,147]. In comparison to simplified models -as, e.g., the TTM and NLM -a major advantage of the TDBE approach consists in enabling the investigation of the coupled non-equilibrium dynamics of electrons and phonon populations with a full momentum resolution. This aspect is particularly important to capture the anisotropic population of the electronic and vibrational states in reciprocal space which may be established out of equilibrium. For example, the absorption of circularly polarized light in transition metal dichalcogenide monolayers is governed by valley-selective circular dichroism [148][149][150]: carriers are photoexcited in either the K or K valleys in the Brillouin zone depending on the light helicity, leading to an anisotropic electronic distribution which cannot be captured via a density-of-state approximations [151]. Similarly, highly anisotropic phonon populations in reciprocal space can be established upon the preferential emission of strongly-coupled phonons [143,144], leading to the characteristic spectral signatures in ultrafast diffuse scattering techniques [83,145,[152][153][154].
This manuscript aims at introducing the fundamentals of model approaches and of the time-dependent Boltzmann equation in the context of the non-equilibrium dynamics of coupled electrons and phonons. In particular, in Sec. 2 we introduce the TTM and in Sec. 2.1 its generalization to several temperature reservoirs, i.e., to the NLM. In Sec. 3, we introduce the TDBE and discuss its theoretical foundation (Sec. 3.1) as well as its relation to the TTM. In Sec. 4, after a concise overview of the experimental signatures of non-equilibrium processes in graphene (Sec. 4.1) we proceed to discuss the theoretical description of the ultrafast carrier and lattice dynamics based on the NLM (Sec. 4.2) and the TDBE (Sec. 4.3) approaches. Summary and outlook are discussed in Sec. 5.

The two-temperature model
The underlying ideas of the TTM [41][42][43][44] are most easily illustrated by considering the thermalization dynamics of two systems S 1 and S 2 . Quantum mechanics is not required at this point; S 1 , S 2 , and their interactions can be assumed to be governed by the laws of thermodynamics. If S 1 and S 2 are initially at thermal equilibrium at the temperatures T 1 and T 2 , respectively, in absence of interactions [ Fig. 1 (a)] the temperature of each subsystem will remain unaltered through time [ Fig. 1 (b)]. If heat can be exchanged [ Fig. 1 (c)], on the other hand, one can expect that for a sufficiently small temperature difference T 2 − T 1 the energy transferred from S 1 to S 2 (from S 2 to S 1 ) during the interval ∆t will be linear in T 2 − T 1 (T 1 − T 2 ), leading to [41]: Here, the energy of S 1 is given by E 1 = c 1 T 1 , where c 1 is the specific heat, and similarly for S 2 . g 1 and g 2 are proportionality constants with units of energy × (temperature × time) −1 . Energy conservation throughout the thermalization dynamics requires ∆E 1 = −∆E 2 , which in turn promptly leads to condition g 1 = g 2 = g. In the limit of infinitesimal time interval, the thermalization dynamics of S 1 and S 2 can thus be formulated by the coupled first-order differential equations: Equations (3) and (4) are the central equations of the two-temperature model. If the temperature dependence of c 1 and c 2 are neglected, Eqs. (3) and (4) admit analytical solutions in the form of decaying exponentials (see Appendix). A detailed discussion of the analytical solution of the TTM can be found in Ref. [57]. The resulting time dependence of the temperatures T 1 and T 2 is illustrated in Fig. 1 (d). In short, the presence of a coupling constant g tends to restore a regime of thermal equilibrium where T 2 = T 1 . More generally, owing to the non-trivial temperature dependence of the specific heat, the TTM must be solved numerically via time-stepping algorithms (e.g., the Euler or Runge-Kutta algorithms). The TTM is often modified to introduce a time-dependent source (driving) term coupled to S 1 or S 2 [50,155]. This scenario is clearly reminiscent of pump-probe experiments [ Fig. 1 (e)] whereby either electrons or lattice are driven out of equilibrium through the coupling to ultrashort pulses. If a coupling in the form S(t) = αe − |t| τ θ(t) is added to Eq. (4), i.e., a laser pulse with the amplitude of α that lasts for time period of τ and starts at t = 0, the model still admits analytic solution (Appendix), leading to the trend illustrated in Fig. 1 (f). For a more detailed discussion on the time-dependent source term and its different functional forms we refer the reader to the Refs. [49,60,118]. The initial increase of T 2 reflects the raise of temperature of S 2 induced by the interaction with a source, whereas at a later stage, the thermalization follows a similar trend to that of Fig. 1 (d).
The TTM can be straightforwardly employed to model the coupled electron-phonon phonon dynamics in condensed matter via the following steps [44,50]: one of the subsystems (S 1 ) is identified with lattice, whereas the other (S 2 ) with the electrons. The temperatures T 1 and T 2 are identified with the effective temperatures of the lattice (T ph ) and electrons (T el ) and the source term S(t) is employed to model the coupling to an external light source. c 1 and c 2 are replaced by the phonon (C ph ) and electron (C el ) heat capacities, which can be expressed as [50]: where D el and D ph are the electron and phonon density of states. Similarly, the coupling constant g can be expressed as [1,50]: where the thermalization rate of electron-lattice system is governed by the electronphonon coupling strength and second moment of the phonon spectrum, which are related to the Eliashberg function α 2 F (Ω) as λ ω 2 = 2 dΩΩα 2 F (Ω). The equations of the TTM, i.e., Eqs. (3) and (4), can thus be expressed as: The theoretical foundation of the TTM, and its relation to the TDBE is discussed in Sec. 3.2. While the electron and phonon heat capacities C el and C ph can be im- mediately obtained from calculations based on density-functional theory and densityfunctional perturbation theory [156], respectively, the parameters g can be estimated from first-principles calculations of the Eliashberg function via well-established simulation packages [157]. This procedure enables the solution of the TTM entirely ab-initio, without resorting to free parameters [50,55,62,68,79]. Alternatively, the g can be deduced by experimental data, e.g., by fitting Eq. (9) to pump-probe photoemission measurements (see, e.g., Sec. 4) [76].
As discussed in Sec. 3.2, the application of the TTM to the non-equilibrium dynamics of electrons and phonons in solids can be justified through its formal derivation from the time-dependent Boltzmann equation [44]. The description of ultrafast processes via Eqs. (3) and (4), however, entails two main approximations: (i) at each time steps throughout the dynamics, electrons are assumed to populate electronic bands according to a Fermi-Dirac function at the effective temperature T el ; (ii) the lattice is assumed to be at thermal equilibrium throughout the dynamics, i.e., all bosonic occupations are described by the Bose-Einstein statistics at the effective temperature T ph . These approximations limit the domain of applicability of the TTM. Because of the approximation (i), the TTM is unsuitable to describe the early stages of the electron dynamics (t < 100 fs), which can be characterized by population inversion, the anisotropic excitation of electron-hole pairs in the Brillouin zone, and electronelectron scatterings. The domain of application of the TTM is thus restricted to metals, semimetals, and doped semiconductors with short electron thermalizaton times, since, on the other hand, electronic excitations in gapped systems (e.g., semiconductors) are inherently linked to a regime of population inversion that cannot be properly modeled via a Fermi-Dirac function. Additionally, the approximation (ii) makes the TTM unsuitable to describe the non-equilibrium dynamics of the lattice. The TTM is sometimes extended to model the population inversion and other forms of nascent non-equilibrium distributions by defining separately electron and hole thermal baths [158], i.e., electron and hole temperatures [159], or by dividing the electronic bath into a majority of thermal and a small portion of non-thermal carriers [49,88,115,116]. However, in these extensions, the issue of thermalized lattice bath (ii) is still present. In the following, we discuss how this limitation can be overcome by extending the TTM to account for anisotropic coupling to different phonon modes.

The non-thermal lattice model
Ultrafast diffuse-scattering experiments and first-principles calculations provide strong evidence that non-thermal regimes of the lattice -i.e., vibrational states characterized by bosonic occupations which deviate significantly from the Bose-Einstein statisticscan be established upon photo-excitation in both semiconducting and metallic layered compounds, such as, black phosphorus [145], MoS 2 [144], graphite [83], graphene [143], and TiS 2 [152]. Even for simple metals such as Al, the anisotropic coupling between electrons and acoustic phonons can trigger the emergence of a non-equilibrium vibrational states persisting for several picoseconds [38]. Generally, whenever the electronphonon interaction is dominated by one or several strongly-coupled modes, these lattice vibrations may provide a preferential decay channel for the relaxation of photo-excited electrons and holes [79]. As mentioned in the introduction, such a scenario can lead to the formation of hot phonons, i.e., a non-thermal state of the lattice [74,76,80,124,125]. Because of the assumption that the lattice can be described by a Bose-Einstein distribution at temperature T ph , the TTM is unsuitable to describe these phenomena [38,83,145] To enable the description of hot phonons and non-thermal states of the lattice, a generalization of the TTM to account for anisotropies in the coupling with different subsets of lattice vibrations -referred to as non-thermal lattice model (NLM) [38,60,120,121] or three-temperature model [79,80], depending on the level of approximation -has recently been proposed. In short, while in the TTM the electrons are coupled to the lattice via a single coupling constant g [ Fig. 2 (a)], in the NLM the lattice is partitioned into phonon groups characterized by distinct coupling constants g µ [ Fig. 2 (b)]. Phonons characterized by higher coupling play a primary role in the electronic relaxation, and therefore exhibit a larger effective temperature on short timescales, whereas on longer timescales thermal equilibrium is restored by phononphonon coupling. This behaviour is schematically illustrated in Fig. 2 (c) for aluminum, where only three acoustic phonons are present.
In general, by partitioning the N ph normal modes of the lattice into N g groups, where 1 < N g ≤ N ph (for N g = 1, the TTM is recovered), the total energy can be expressed as E ph = Ng µ=1 E µ ph . The rate of change of the energy of the µ-th phonon group can be expressed as ∂ t E µ ph = C µ ph ∂ t T µ ph , where C µ ph is the specific heat of the µ-th phonon group, which is obtained by replacing the group density of states in Eq. (6). The NLM can thus be formulated as a set of coupled first-order differential equations which relates the electronic temperature T el to the temperatures T µ ph of the N g phonon groups: Here, the coupling constant g µ is defined as in Eq. (7) by restricting the sum to all phonons in the µ-th group, τ µµ denotes the time constant for the decay of the µ-th subgroup due to phonon-phonon interaction with phonons in the µ -th group, and it can be obtained from first-principles via calculations of the phonon-phonon scattering matrix elements [60]. As illustrated in Fig. 2 (c) for Al [38], at variance with TTM, the NLM enables to account for the establishment of a non-thermal regime in the lattice. In particular, following photo-excitation of the electrons, energy is transferred to each phonon group µ proportionally to the coupling constant g µ , thus, leading to a discrepancy among the vibrational temperatures T µ ph . On longer timescales, the onset of phonon-phonon scattering (via the second term in Eq. (10)), drives the lattice back to a thermalized regime, where all vibrational temperatures coincides.
With the exception of elemental metals, where the identification of different phonon groups is straightforward due to the reduced vibrational degrees of freedoms (see, e.g., Fig. 2), in compounds with several atoms in the unit cell the definition of phonon groups is to some extent arbitrary, and it represents a shortcoming of the NLM. In some limiting cases (e.g., cuprates, graphene, and MgB 2 [74,79,80]), it is possible to distinguish strongly and weakly coupled modes and, correspondingly, phonons can be grouped according to their coupling strength. The arbitrariness in the definition of phonon groups can be lifted by resorting to a first-principles description of the ultrafast dynamics of electrons and phonons based on the time-dependent Boltzmann equation.

The time-dependent Boltzmann equation
The time-dependent Boltzmann equation (TDBE) constitutes an optimal compromise between accuracy and efficiency to investigate the ultrafast dynamics of coupled electron-phonon systems. In the TDBE, the dynamics of electronic and vibrational excitations are described by changes of the electron and phonon distribution functions f nk (t) and n qν (t), respectively, whereas electron and phonon energies are left unchanged throughout the dynamics. At thermal equilibrium, f nk and n qν are time independent and they coincide with the Fermi-Dirac and the Bose-Einstein occupations f 0 nk and n 0 qν : Here, ε F is the Fermi energy, ε nk is the single-particle energy of a Bloch electron, and ω qν the phonon energy. This case is exemplified by the left panel of Fig. 3, where the Fermi-Dirac occupations are superimposed to the band structure of monolayer MoS 2 , with yellow (blue) denoting fully occupied (empty) states with f nk = 1 (f nk = 0). In this framework, a regime of non-equilibrium requires either f nk or n qν (or both) to differ from the equilibrium Fermi-Dirac and the Bose-Einstein occupations, as illustrated in the right panel of Fig. 3. The non-equilibrium distributions change over time, and their dynamics is determined by the TDBE: where ∂ t = ∂/∂t and Γ nk and Γ qν denote the collision integrals for electrons and phonons. The numerical solution of Eqs. (14) and (15) requires the development of suitable approximations for the evaluation of the collision integrals. In short, Γ nk and Γ qν account for the several scattering mechanisms which may lead to changes of the distributions functions as, e.g., electron-electron, electron-phonon, phonon-phonon, and impurity scattering as well as the coupling to external fields. The recent development of electronic structure codes for the study of electron-phonon and phonon-phonon coupling has enabled to estimate the contribution of these scattering processes to collision integrals, enabling the investigation of the coupled electron-phonon dynamics entirely from first principles [33,59,134,143,144].

First-principles expressions for the collision integrals
In the following, we outline the derivation of the collision integral due to the electronphonon and phonon-phonon interactions from Fermi's golden rule. A similar treatment can be generalized to other contributions to the collision integral as, e.g., electronelectron and impurity scattering or coupling to external fields. The Hamiltonian for an anharmonic crystal in presence of electron-phonon and phonon-phonon interactions can be expressed as:Ĥ HereĤ e = nk ε nkĉ † nkĉ nk is the electronic Hamiltonian, whereĉ † nk andĉ nk are fermionic creation and annihilation operators, respectively.Ĥ p = qν ω qν [â qνâ † qν + 1/2] is the Hamiltonian of the lattice in the harmonic approximation.â † qν andâ qν are phonon creation and annihilation operators. The eigenstates ofĤ p can be expressed as |χ s = qν |n s qν , where |n s qν are the eigenstates of the number operator N qν =â † qνâqν with eigenvalue n s qν . The quantity n s qν is the occupation number of the bosonic mode characterized by quantum numbers ν and q. The superscript s-th indicates that the occupations are relative to the s-th eigenstate |χ s of the harmonic Hamiltonian. With this notation, the energy of the harmonic lattice can be expressed as E s = χ s |Ĥ p |χ s = qν ω qν [n s qν + 1/2]. The third term in Eq. (16) is the electron-phonon coupling Hamiltonian: where N p is the number of unit cells in the Born-von Kármán (BvK) supercell [1]. The electron-phonon coupling matrix elements g ν mn (k, q) can be derived from first principles within the framework of density-functional perturbation theory and they are defined as g ν mn (k, q) = ψ mk+q |∆ qνv KS |ψ nk , where ∆ qνv KS is the linear change of the effective Kohn-Sham potential v KS due to a phonon perturbation, and |ψ nk are single-particular orbitals, solutions of the single-particle Kohn-Sham equations [1,156].
Finally, the last term in Eq. (16) is the phonon-phonon coupling Hamiltonian, which arises from anharmonicities of the lattice and it can be expressed as [160]: Ψ νν ν qq q denotes the phonon-phonon scattering matrix elements, which is related to the probability amplitude of three-phonon scattering processes. The derivation of the collision integrals for the Hamiltonian in Eq. (16) begins with the observation that electron-phonon and phonon-phonon interactions in solids are typically weak and, thus, the termsĤ ep andĤ pp in the Hamiltonian in Eq. (16) can be treated as perturbations. Correspondingly, the rate Γ i→f of transitions from an initial state |i to a final state |f can be obtained via Fermi's golden rule: where E tot i and E tot f are the total energies of the initial and final states, respectively, andV is an arbitrary perturbation. To proceed further, we focus on the electronphonon interaction and we consider initial and final states in the form of a Born-Oppenheimer ansatz as |i = |Ψ i |χ i , where |χ i and |Ψ i are eigenstates of the harmonic HamiltonianĤ p and of the electronic HamiltonianĤ e , respectively. The matrix elements of the electron-phonon coupling Hamiltonian can be expressed as: The matrix elements of fermionic operators in Eq. (20) differs from zero only if |Ψ i and |Ψ f differ only in the occupation of the nk and mk + q states. In such case, it yields: Similarly, from elementary considerations on the action of bosonic operators on the vibrational eigenstates |χ s , one can deduce that the matrix element χ f |â qν |χ i = √ n qν (and similarly, χ f |â † −qν |χ i = n −qν + 1) if the state |χ f differs from |χ i only in the occupation of the phonon qν (−qν) such that n f qν = n i qν −1 (and n f qν = n i qν +1), and it vanishes otherwise. The total rate (number of events per unit time) for the absorption of a phonon with quantum numbers qν can be directly derived from Eq. (19) making use of Eqs. (20) and (21) and summing over all initial and final electronic states: The energy difference between the initial and final states has been expressed as The diagrammatic representation of a phonon absorption process of this kind is reported in Fig. 4 (b). The same procedure can be repeated by considering phonon emission processes [ Fig. 4 (a)], yielding: The total rate of change in the phonon occupation n qν due to the electron-phonon interactions Γ ep qν can thus be defined as the difference between the rates of phonon (36) represents the difference of the rate for an electron in state |nk to scatter out of the state ( rst two terms) and the rate for an electron to scatter into the state |nk (last two terms). Both processes can be mediated either by phonon absorption ( rst and third term) or phonon emisssion (second and forth term). We note that we let q → −q in the terms involving phonon emission to write them also in terms of f mk+q instead of f mk−q , making use of ω qν = ω −qν and the fact that the matrix elements for phonon emission and absorption are related by complex conjugation. The four scattering processes included in Γ (co) nk are illustrated in gure 2.
Equation (28) is solved iteratively to obtain the Eelddependent occupancies f nk . Then the experimentally accessible macroscopic average of the current density J(r) can be obtained via where we made use of equations (5), (21), and (29) and where V and V uc denote the crystal and unit cell volume, respectively. In equation (38) we introduced the diagonal velocity the E-eld dependence of all quantities for clarity. In the case of weak electric elds, we can restrict ourselves to the linear response of the current density, which de nes the conductivity tensor: Here α, β run over the three Cartesian directions and we introduced the short-hand notation ∂ E β f nk = (∂f nk /∂E β )| E=0 . From equation (28), we can obtain an expression for the linear response coef cients ∂ E β f nk by taking derivatives on both sides with respect to the electric eld: where we introduced the partial decay rate and its analog τ −1 mk+q→nk with the indices nk and mk + q swapped. Here, f 0 nk denotes the equilibrium occupancies in the absence of an electric eld, which are given by the Fermi-Dirac distribution evaluated at the band energies, where µ is the chemical potential. We also used the fact that, ignoring the Berry curvature [61], the diagonal matrix elements of the velocity operator are simply given by v α nk = −1 ∂ε nk /∂k α . Equation (40) is known in the literature [60] as the Boltzmann transport equation. Its solution yields the linear response coef cients ∂ E β f nk , which are needed in equation (39) to obtain the conductivity tensor.
The electrical conductivity in equation (39) scales with the density of carriers. This is generally not an issue when studying metals, for which temperature, bias voltage, and defects do not alter the carrier density near the Fermi energy. However, in semiconductors the carrier density can change by many orders of magnitude with doping, temperature, and applied voltage. In these cases, in order to single out the intrinsic transport properties of the material, it is convenient to introduce the carrier drift mobility, which is de ned as the ratio between conductivity and carrier density: The charge carrier density entering the electron mobility tensor, n c = n el , is de ned as Figure 2. The four processes included in the collision rate in equation (25) derived from the Fan-Migdal self-energy: Scattering of an electron out of state |nk via phonon absorption (green, rst term) and emission (purple, second term) and scattering of an electron into state |nk via phonon absorption (brown, third term) and emission (orange, fourth term). (36) represents the difnce of the rate for an electron in state |nk to scatter out e state ( rst two terms) and the rate for an electron to ter into the state |nk (last two terms). Both processes can ediated either by phonon absorption ( rst and third term) honon emisssion (second and forth term). We note that we → −q in the terms involving phonon emission to write also in terms of f mk+q instead of f mk−q , making use of = ω −qν and the fact that the matrix elements for phonon sion and absorption are related by complex conjugation. four scattering processes included in Γ (co) nk are illustrated gure 2. quation (28) is solved iteratively to obtain the Eeldndent occupancies f nk . Then the experimentally accesmacroscopic average of the current density J(r) can be ined via re we made use of equations (5), (21), and (29) and where nd V uc denote the crystal and unit cell volume, respecy. In equation (38) we introduced the diagonal velocity matrix elements v nk = nk|p/m|nk and explicitly indicated the Eeld dependence of all quantities for clarity. In the case of weak electric elds, we can restrict ourselves to the linear response of the current density, which de nes the conductivity tensor: Here α, β run over the three Cartesian directions and we introduced the short-hand notation From equation (28), we can obtain an expression for the linear response coef cients ∂ E β f nk by taking derivatives on both sides with respect to the electric eld: where we introduced the partial decay rate and its analog τ −1 mk+q→nk with the indices nk and mk + q swapped. Here, f 0 nk denotes the equilibrium occupancies in the absence of an electric eld, which are given by the Fermi-Dirac distribution evaluated at the band energies, where µ is the chemical potential. We also used the fact that, ignoring the Berry curvature [61], the diagonal matrix elements of the velocity operator are simply given by v α nk = −1 ∂ε nk /∂k α . Equation (40) is known in the literature [60] as the Boltzmann transport equation. Its solution yields the linear response coef cients ∂ E β f nk , which are needed in equation (39) to obtain the conductivity tensor.
The electrical conductivity in equation (39) scales with the density of carriers. This is generally not an issue when studying metals, for which temperature, bias voltage, and defects do not alter the carrier density near the Fermi energy. However, in semiconductors the carrier density can change by many orders of magnitude with doping, temperature, and applied voltage. In these cases, in order to single out the intrinsic transport properties of the material, it is convenient to introduce the carrier drift mobility, which is de ned as the ratio between conductivity and carrier density: The charge carrier density entering the electron mobility tensor, n c = n el , is de ned as re 2. The four processes included in the collision rate in tion (25) derived from the Fan-Migdal self-energy: Scattering electron out of state |nk via phonon absorption (green, rst ) and emission (purple, second term) and scattering of an ron into state |nk via phonon absorption (brown, third term) mission (orange, fourth term).
Prog. Phys. 83 (2020) 036501 (36) represents the d ference of the rate for an electron in state |nk to scatter o of the state ( rst two terms) and the rate for an electron scatter into the state |nk (last two terms). Both processes c be mediated either by phonon absorption ( rst and third ter or phonon emisssion (second and forth term). We note that let q → −q in the terms involving phonon emission to wr them also in terms of f mk+q instead of f mk−q , making use ω qν = ω −qν and the fact that the matrix elements for phon emission and absorption are related by complex conjugati The four scattering processes included in Γ (co) nk are illustra in gure 2.
Equation (28) is solved iteratively to obtain the Ee dependent occupancies f nk . Then the experimentally acc sible macroscopic average of the current density J(r) can obtained via where we made use of equations (5), (21), and (29) and wh V and V uc denote the crystal and unit cell volume, resp tively. In equation (38) we introduced the diagonal veloc  (25) derived from the Fan-Migdal self-energy: Scatterin of an electron out of state |nk via phonon absorption (green, rst term) and emission (purple, second term) and scattering of an electron into state |nk via phonon absorption (brown, third term) and emission (orange, fourth term).
Rep. Prog. Phys. 83 (2020) 036501 (36) represents the difference of the rate for an electron in state |nk to scatter out of the state ( rst two terms) and the rate for an electron to scatter into the state |nk (last two terms). Both processes can be mediated either by phonon absorption ( rst and third term) or phonon emisssion (second and forth term). We note that we let q → −q in the terms involving phonon emission to write them also in terms of f mk+q instead of f mk−q , making use of ω qν = ω −qν and the fact that the matrix elements for phonon emission and absorption are related by complex conjugation. The four scattering processes included in Γ (co) nk are illustrated in gure 2.
Equation (28) is solved iteratively to obtain the Eelddependent occupancies f nk . Then the experimentally accessible macroscopic average of the current density J(r) can be obtained via where we made use of equations (5), (21), and (29) and where V and V uc denote the crystal and unit cell volume, respectively. In equation (38) we introduced the diagonal velocity the E-eld depend In the case of w to the linear respon conductivity tensor Here α, β run ove introduced the shor From equation ( The charge carrier sor, n c = n el , is de  (25) derived from the Fan-Migdal self-energy: Scattering of an electron out of state |nk via phonon absorption (green, rst term) and emission (purple, second term) and scattering of an electron into state |nk via phonon absorption (brown, third term) and emission (orange, fourth term). Figure 5. Schematic representation of the four phonon-assisted scattering processes included in the electron collision integral due to the electron-phonon interaction (Eq. (25)). Adapted from Ref. [8].
emission (Γ em qν ) and absorption (Γ abs qν ) processes: Equation (24) is the phonon collision integral due to the electron-phonon interaction. Its time dependence arises from the changes of the electron and phonon distribution functions (f nk and n qν ) over time. In conditions of thermal equilibrium between electrons and the lattice, as for instance in an ideal situation in which electron and phonon occupations are described by Fermi-Dirac and Bose-Einstein statistics at a given temperature, the rates Γ abs qν and Γ em qν are equal and opposite in sign, indicating that the total change in the phonon number n qν vanishes, since the emission and absorption of phonons are perfectly compensated.
Following similar steps, the electronic collision integral due to the electron-phonon interaction can be derived as: Each term in this expression is directly related to phonon-assisted electronic transitions involving states in the vicinity of the Fermi surface, as depicted in Fig. 5. In particular, the first and second terms in Eq. (25) arise from processes in which an electron scatter from mk + q in to nk via the emission and absorption of a phonon, respectively. The third and fourth terms arise from the scattering from nk into mk + q due to phononemission and absorption processes. In analogy to Eq. (24), at thermal equilibrium, scattering processes in and out of nk balance each other, leading to Γ ep nk = 0. Note that the rate Γ ep nk of electron population relaxation defined with Eq. (25) and the corresponding lifetime τ ep nk = /Γ ep nk should not be confused with the quasiparticle decay rate and lifetime as obtained from the imaginary part of the electron self-energy [161]. It is the former and not the latter that is usually extracted from pump-probe measurements.
The derivation of the scattering rate due to the phonon-phonon scattering involves the tedious (but otherwise straightforward) evaluation of several matrix elements of bosonic operators. The result can be recast in the form: where the modified Kronecker's δ G q equals unity if q = 0 or q = G, where G is a reciprocal-lattice vector, and it is zero otherwise. Equations (26) vanishes identically if the lattice is at thermal equilibrium.
Combining Eqs. (25)- (26), the time-dependent Boltzmann equation can be rewritten as: A numerical procedure to solve Eqs. (27) and (28) consists in: (i) defining an initial electronic (or vibrational) excited state characterized by electronic (vibrational) occupations which differs from the equilibrium ones (as e.g. in Fig. 3); (ii) solve the differential equation using iterative methods (as, e.g., the Euler of Runge-Kutta methods); (iii) update the collision integrals via Eqs. (25)-(26) at each time step. Besides the aforesaid microscopic scattering process, electron-electron interaction plays a major role in thermalization of electron-lattice system, especially in the 10fs timescale. The corresponding collision integral Γ ee nk [f nk (t)] enters Eq. (27) and it is usually defined in terms of electron scatterings {k, k } → {k + q, k − q} (and vice versa) mediated by the statically screened Coulomb interaction. For a detailed description of electron-electron scattering processes we refer the reader to Refs. [121,132,134,162]. Note that the electron-plasmon scatterings, as a part of the dynamical electron-electron interaction, is as well considered to be essential in the early stage of electron thermalization process [126,163,164], however, these dynamical effects are rearly taken into account when studying hot carrier thermalization.

Relation between the TTM and the time-dependent Boltzmann equation
As demonstrated by Allen [44], the thermalization dynamics of electrons and phonons in presence of the electron-phonon interaction can be recast in the form of a TTM, as formulated via Eqs. (3) and (4), starting entirely from first principles. In short, by expressing the coupling parameter g in terms of the Eliashberg function α 2 F , all free parameters of the models are fixed. A detailed derivation of this scheme can be found, for instance, in Refs. [60,121].
In the following, we report an alternative derivation of the TTM with scope of emphasizing the link with the phonon lifetime, as obtained from phonon self-energy due to the electron-phonon interaction. We begin by considering the total energy of a set of non-interacting electrons (E el ) and phonons (E ph ): The rate of change of the energies can be expressed as: where ∂ t = ∂/∂t, the left-hand side has been rewritten making use of the chain rule ∂ t E = ∂E ∂T ∂T ∂t , and the heat capacity C el(ph) = ∂E el(ph) /∂T el(ph) has been introduced. Equations (31) and (32) rely on the assumption that electron energies ε nk and phonon frequencies ω qν do not depend on time. The time derivative of the fermionic and bosonic distribution functions entering Eqs. (31) and (32) can be obtained via the time-dependent Boltzmann equation.
A central assumption of the non-thermal lattice model is that, at each time t, the electronic occupation f nk can be approximated by a Fermi-Dirac function at temperature T el : Similarly, it is assumed that the lattice remains at thermal equilibrium throughout the dynamics, that is, all vibrational modes are at the same temperature T ph (a condition that may be violated in the non-equilibrium dynamics of the lattice following photoexcitation [144]). Under these assumptions, the arguments of the collision integrals can be simplified. For instance, the argument of Eq. (24) can be rewritten as [60]: where we introduced the Bose-Einstein function: Additionally, we set ε mk+q − ε nk = ω qν because of the Dirac-δ in Eq. (24), and we used n qν = n BE (ω qν , T ph ). A Tailor expansion of n BE to first order yields: In the last equality, we introduced the specific heat of a single harmonic oscillator c qν = ω qν ∂n qν /∂T .

13
Making use of Eqs. (34), the integrand within the phonon-electron collision integral in Eq. (24) can be rewritten as: Where we introduced the phonon lifetime due to the electron-phonon interaction [1]: Via the identity τ −1 qν = −2 Im Π qν =, Eq. (36) further establishes a simple relation between the phonon self-energy due to the electron-phonon interaction Π qν and the corresponding collision integral. The limit of validity of this identity are defined by the assumptions made thus far. Combining (37) with Eq. (32): A similar identity can be derived by applying the same procedure for the electronic term. After introducing the effective electron-phonon coupling constant g: the equations of the TTM, Eq. (3) and (4) are recovered.

Experimental signature of hot-carrier dynamics in graphene
The ultrafast dynamics of Dirac carriers in pump-probe experiments is usually approximated as a four-step process consisting of (i) photo-excitation, (ii) formation of a quasi-equilibrium state, (iii) full electron thermalization, and (iv) energy transfer to the lattice (acoustic phonons). In the step (i), the interaction with a pump pulse drives the system into an electronic excited state, with carriers photo-excited above the Fermi level. With this a non-equilibrium electron distribution is created. The step (ii) involves electron-electron scatterings, impact ionization, and Auger processes (timescale of ∼ 10 fs) [162,191], which promote photo-excited electrons and holes towards the Fermi level. Some experiments indicate that under suitable conditions a regime of population inversion can be established, i.e., two electron distributions with separate chemical potentials (and temperatures) [77,192]. Interestingly, some theoretical works suggested that strong scatterings between non-equilibrium electrons and phonon are active already at this stage [132,193]. Even more, this early electronphonon scattering process might be responsible for majority of energy flow, leaving a small portion of excess energy for the later scatterings between equilibrium electrons and phonons [132]. In step (iii), Auger recombination [172,194], scatterings with optical phonons [76,122], as well as plasmon emission [163,164] bring the electrons to a thermalized regime (∼ 100 fs), where the electronic distribution function is described by a high-temperature Fermi-Dirac function. Finally, in step (iv), the photo-excited electrons and holes dissipate their energy via phonon-assisted scattering processes, mostly via the acoustic modes (∼ 1 − 10 ps). Supercollisions [76] and the formation of hot phonons [195] have been proposed as underlying physical processes to explain the energy transfer to the lattice [122,196]. Further experimental investigations on the parent compound graphite corroborated the hot-phonon picture [37,197]. Overall, these studies revealed a complex non-equilibrium dynamics characterized by the coexistence of several scattering mechanisms which are challenging to decipher based on purely experimental investigations.  Refs. [76,77]. Simulations based on the NLM (more precisely, three-temperature model) with first-principles parameters are reported as continuous lines. Panels (a-b) and (c-d) are adapted from Ref. [68] and Ref. [79], respectively.
[77]. As a result of substrate-induced hole-doping, the Fermi energy is located 0.2 eV below the Dirac point [ Fig. 6 (a)]. Figures 6 (b-d) illustrate several time snapshots of the tr-ARPES spectral function in the vicinity of the Dirac point. Before excitation [ Fig. 6 (b)], the spectral intensity reflects occupied electronic states at equilibrium, with broadening arising from finite temperature and experimental resolution. Owing to the optical selection rules governing the coupling with linearly-polarized light, the spectral intensity is dominated by contributions arising from the left sub-band [198]. indicates the recovery of an equilibrium regime. High-quality tr-ARPES measurements can be analysed to extract the effective electronic temperatures T el and its time dependence [71,72,76,77,123,[199][200][201], thus, establishing a direct link with the TTM and NLM discussed in Secs. 2 and 2.1. For a given time delay, T el can be determined by fitting the normalized momentum-integrated photoemission intensity with a Fermi-Dirac function. This procedure is illustrated in Fig. 6 (e) where the fit Fermi-Dirac function (continuous line) is superimposed to the measured energy distribution curves (dots). The time dependence of T el , shown in Fig. 6 (f), is obtained by repeating this procedure for several measured time delays and it closely resembles the trend reported for the light-driven TTM illustrated in Fig. 1 (c): the photo-excitations of the electrons by the pump pulse manifests itself through the initial increase of electronic temperature, whereas the subsequent cooling reflects the thermalization dynamics due electron-phonon scattering processes.

Theoretical modelling of carrier thermalization in graphene
First-principles calculations of the electron-phonon interaction can be combined with the time-propagation algorithms discussed in Secs. 2-3 to investigate the origin of the hot-carrier dynamics and its fingerprints in spectroscopy. The band structure and and probe beams: before ( t < 0) and after ( t = 0) pump, during ( t = 0.5 ps) and after ( t = 2.5 ps) the carrier thermalization. Calculated di↵erential tr-ARPES intensity relative to equilibrium for t = 0 (f) and 0.5 ps (g). Di↵erential pump-probe signal N (E, t) (integrated over momentum) as obtained from our calculations (h) and from tr-ARPES measurements (i), reproduced from Ref. [2]. (j) Electron linewidths due to electron-phonon interaction at various time delays.
of thermal temperatures. Conversely, these features reveal a time-dependent renormalization of the quasiparticle energies, which we attribute to the increase of the population of the acoustic phonons resulting from the higher e↵ective temperatures. This finding is compatible with recent observations in transient optical measurements [54,55] of narrowing energy gap between the ⇡ and ⇡ ⇤ bands.
This e↵ect can be quantified by inspecting the change in the real part of the electron self-energy, ⌃ = Re[⌃ nk (" nk , t) ⌃ eq nk (" nk )], where ⌃ eq denotes the Fan-Migdal self-energy of the system at equilibrium. ⌃ provides information pertaining to the time-dependent change of the quasiparticle energies due to electronphonon interaction. For " nk = 1 eV, we obtain ⌃ = 10 meV at t = 0.5 ps. On the other hand, at t = 0 ps, when T sc approaches its maximum value, the quasiparticle energies are almost unchanged ( ⌃ < 1 meV). This suggests that the energy renormalization (at these lower energies) stems from the electronic coupling to acoustic modes. To further validate this argument, we note that, while the population of acoustic phonons may be easily enhanced by thermal excitation due to their low energy, the opposite is true for the strongly coupled modes. In particular, owing to the high energy of the longitudinal optical phonons (~! LO ' 200 meV), their population is not significantly enhanced by the increase of the e↵ective temperature T sc , and therefore it is not expected to introduce visible changes in the spectral function.
To enable the comparison with tr-ARPES measurements across the full range of time delays, we evaluate the di↵erential momentum-integrated pump-probe signal N (!, t) = R dk[A k (!, t) A eq k (!)], where A eq k denotes the spectral function at equilibrium. Our calculations for N (!, t) are reported in Fig. 2 (h), whereas experimental data from [2] are shown in Fig. 2 (i). The comparison with experiment demonstrates good agreement across the full range of time delays.
Finally, we discuss the time dependence of the electron linewidths due to electron-phonon interaction, illustrated in Fig. 2 (j). Close to the Fermi energy, the linewidths exhibit a fast increase after the pump, followed by a slow decrease. After 2.5 ps, the linewidth coincides with the one preceding photo-excitation ( t < 0), marking the return to equilibrium. This behavior can be ascribed to the time-dependent change of the phase phonon dispersion of hole-doped graphene is illustrated in Fig. 7 (a) and (b), respectively [68], whereas the phonon density of state (F ) and the Eliashberg function (α 2 F ) are shown in the right panel of Fig. 7 (b). The Eliashberg function reflects the weighted phonon density of states, to which individual phonons contribute according to their electron-phonon coupling strength [202]. The sharp peaks in α 2 F at 160 and 190 meV indicate that the electron-phonon coupling arises primarily from longitudinal optical (LO or E 2g ) phonons at Γ and transverse optical (TO or A 1 ) modes at K [green dots in Fig. 7 (b)], whereas remaining phonons couple relatively weakly to the electrons. This finding suggests that a few phonons are likely to provide a preferential decay channel for the relaxation of excited carriers, leading to the emergence of hot phonons -namely, phonons characterized by a higher vibrational temperature. The pronounced anisotropy of the electron-phonon interaction in graphene can be explicitly accounted for by formulating a NLM (i.e., a three-temperature model) which discriminates between strongly-(sc) and weakly-coupled (wc) phonons [79], and by determining the model parameters from first-principles calculations [156,157,203]. The effective electronic temperature T el and the vibrational temperatures (T wc and T sc ) obtained from the solution of the NLM [Eqs. (10) and (11)] are illustrated in Figs. 7 (c) and (d) for different excitation conditions [79]. The electronic temperature follows the characteristic trend expected for photoexcited electrons and it is in good agreement with the experimental values extracted from tr-ARPES. The preferential excitation of stronglycoupled modes upon carrier relaxation leads to a pronounced rise of T sc , whereas the vibrational temperatures of weakly-coupled modes T wc is smaller.
By post-processing the effective temperatures T el , T wc , and T sc , one can define a simple procedure to estimate the transient changes of many-body interactions and their signatures in photoemission. The self-energy due to electron-phonon interaction can be modified to account for changes of the electronic and vibrational distribution functions as: The time dependence arise from the dependence of the electron and phonon distri-17 bution functions on the effective temperatures, given by f nk = [e εnk/kBTel + 1] −1 and n νq = [e ωνq/kBTν − 1] −1 , respectively, where T ν = T sc (T ν = T wc ) for the strongly (weakly) coupled modes. The tr-ARPES spectral function can be directly derived from Eq. (40) via A k (ω, ∆t) = π −1 n Im [ ω − ε nk − Σ nk (ω, ∆t)] −1 . The spectral function of graphene obtained from this procedure is illustrated in Fig. 8, and it reproduces the main spectral signatures of photoexcitations revealed by experiments as well as the characteristic time scales of photoexcitations.
Other time-dependent physical features of ultrafast hot carrier cooling visible in pump-probe spectroscopy can be theoretically captured in a similar fashion, i.e., by benefiting from the time dependence of electron and phonon temperatures. For instance, time dynamics of electron excitations (e.g., plasmons and screened interband transitions) immediately after the laser excitation can be simulated via optical conductivity σ(ω; T el , T ph ) or dielectric function (q, ω; T el , T ph ) [88,124]. The spectral function of electron excitations (or the so-called electron energy loss function) is defined as [68,204] −Im where Ω 2 (q, ω; T el , T ph ) = v(q)q 2 Re π αα (q, ω; T el , T ph ) and Γ(q, ω; T el , T ph ) = −v(q)q 2 Im π αα (q, ω; T el , T ph ), define the energy and corresponding damping (e.g., Landau damping and electron-phonon channel) of electronic excitations as a function of time delay (via time dependence of T el and T ph ). v(q) is the Coulomb interaction and π αα is the current-current response function (photon self-energy) for the polarization direction α. Useful optical quantities such as optical absorption σ(ω; T el , T ph ) = iπ αα (q = 0, ω)/ω or photoconductivity ∆σ(ω; T el , T ph ) = σ(ω; T el , T ph ) − σ(ω; T el = T ph = T 0 ) can also be simulated in the same way. The combination of Eq. (41) and the TTM equations was recently utilized in order to track ultrafast dynamics of Dirac plasmon in graphene, including energy loss of plasmon in non-equilibrium regime [68]. It was concluded that although steady state plasmon is mostly damped due to scatterings with acoustic phonons, the non-equilibrium Dirac plasmon dissipates most of its energy on intrinsic strongly coupled optical phonons. The same approach was also applied to study laser-induced band renormalizations close to van Hove singularity point of graphene band structure [205,206], where it was concluded that electron-phonon coupling plays a significant role in photo-induced band modifications in graphene [95]. Note that the approach of incorporating time dependence into Eqs. (40) and (41) via the TTM results is actually general and can be applied to any spectral function or selfenergy defined in terms of electron and phonon distribution functions. For instance, this approach was exploited to monitor hot-phonon dynamics in conventional superconductor MgB 2 via phonon spectral function B(ω; T el , T ph ) and many-body phonon self-energy Π(ω; T el , T ph ) of strongly-coupled E 2g mode [80]. All in all, this way of modelling ultrafast carrier and lattice thermalization can provide many microscopic insights of spectroscopic quantities in time domain, and was proven to have a wellbalanced ratio between accuracy and numerical cost. phonon populations (see Fig. 2). In the first 300 fs, the excess electronic energy rapidly generates modes with strong e-ph coupling, including the A 1 and E2g optical modes, respectively with wavevectors near the K and points of the BZ. A slower change in phonon populations occurs between 1 and 5 ps, when the optical phonons decay into LA and TA modes with wavevectors halfway between and the K or M points of the BZ. Long-wavelength (small-q) acoustic phonons are also generated in the same time window. A second stage of phonon equilibration sets in after ∼20 ps, when energy redistribution occurs primarily within the LA and TA branches (see Fig. 2).
We carry out a mode-by-mode analysis of the phonon populations and effective temperatures, taking into account mode crossing as discussed in Appendix B. In Fig. 3(a), we plot for each mode the excess population relative to equilibrium, where Nν (t ) are BZ-averaged phonon populations for each mode. Also shown in Fig. 3(a) are the mode-resolved e-ph coupling strengths, g 2 ν (see Appendix B). Due to their strong e-ph coupling, the in-plane longitudinal optical (LO) and transverse optical (TO) modes are rapidly emitted during the initial fast carrier cooling, highlighting the key role of e-ph interactions at early times [4].
The LO, TO, and LA modes are populated extensively before 1 ps, while the TA and ZO modes are excited more gradually through ph-ph interactions over timescales longer than 1 ps. Due to their weak e-ph coupling, the out-of-plane flexural phonon modes (ZA and ZO) are generated at a significantly slower rate than in-plane phonons, resulting in a slow rise of their populations in the first 10 ps. The ZA mode, with the weakest e-ph scattering strength, exhibits the slowest rise in population among all modes. In turn, the weak e-ph coupling and slow generation of flexural phonons is responsible for a carrier cooling bottleneck: As electrons and holes interact with hot ZA and ZO populations, they gain energy from the flexural modes for over 10 ps, well beyond the initial subpicosecond carrier cooling.
Although at short times the phonon distributions are still nonthermal, we find that after 2-5 ps all phonon populations are thermal and well approximated by hot Bose-Einstein distributions. The mode-resolved effective phonon temperatures are computed as BZ averages of state-dependent temperatures (see Appendix A) and shown in Fig. 3(b). We find that the effective phonon temperatures are strongly mode dependent throughout the simulation, their trend mirroring the respective FIG. 2. Phonon scattering channels. Optical phonons at and K are rapidly populated in the first 300 fs due to their strong coupling with the excited carriers, but later decay into lower-energy optical and acoustic modes through ph-ph processes. The dominant energy and momentum-conserving phonon decay channels are shown with red arrows at 5-and 20-ps delay times after the initial electronic excitation. play only after 3 ps. Throughout the simulation, low-energy electron-hole pairs are excited near the Dirac cone via ph-e processes, resulting in a modest ultrasonic attenuation of the phonon dynamics. Overall, the excess energy initially imparted to the excited phonons mainly dissipates through slow ph-ph processes, the entire phonon relaxation taking tens of picoseconds. Our approach is uniquely able to shed light on these longer timescales characteristic of phonon dynamics.

III. DISCUSSION
A key advantage of our rt-BTE approach is the possibility, rather unique among existing first-principles methods for ultrafast dynamics, to validate the interactions against experiments, thus guaranteeing the quantitative accuracy of our simulated dynamics. We validate the e-ph and ph-ph interactions employed in our calculations, respectively, by computing electrical and heat transport properties of graphene (see Appendix A). We obtain a room-temperature electron mobility of 186000 cm 2 /V s, in excellent agreement with experimental results in suspended graphene [61,62]. We also compute the thermal conductivity in the single-mode relaxation time approximation (RTA) of the BTE, obtaining a room-temperature value of 482 W/mK in excellent agreement with previous RTA calculations [63]; this result implies that the full solution of the BTE (beyond the RTA) would give a thermal conductivity consistent with experiment [63]. These results show that the interactions employed in our ultrafast dynamics are precise, so our computed timescales are expected to be quantitatively accurate.
Finally, note that our method treats e-ph and ph-ph interactions perturbatively, using values for these interactions computed on the equilibrium structure. Therefore our approach is most suitable to model the low-excitation regime in which the atomic positions remain fairly close to their equilibrium values; for example, in our calculations the phonon temperatures do not rise significantly abov temperature. Regimes involving intense e transitions, atomic diffusion, or material d ferent treatments. Future work will add issue of updating the e-ph and ph-ph inte nonequilibrium dynamics.

IV. CONCLUSION
In summary, we have shown a versatil work, rt-BTE, for modeling the ultrafast of electrons and phonons. We demonstrat broad applicability through simulations of troscopy, x-ray diffuse scattering, and stru dynamics. Our results shed light on e-p tering processes governing ultrafast dyna and provide valuable information for the electronic and optical devices. We plan to m method available in a future release of PERT the community with a tool for simulatin ultrafast time-domain experiments. Future to treat explicitly the light excitation puls spin to unravel the intertwined nonequili electronic, structural, and spin degrees of gether, our first-principles approach demo shift in computing the ultrafast electron and dynamics, bridging the gap between theo and enabling quantitative predictions of u in materials.

Density functional theo
We carry out first-principles DFT calcu with a relaxed in-plane lattice constan graphene sheet is separated from its per 9-Å vacuum. The ground-state electronic puted using the QUANTUM ESPRESSO code approximation (LDA) of DFT. We use pseudopotential, a 90-Ry plane-wave kin and a Methfessel-Paxton smearing of 0.0 state charge density is obtained using a 3 grid, following which a non-self-consis employed to obtain the Kohn-Sham eige functions on a 12 × 12 × 1 k-point grid. phonon populations (see Fig. 2). In the first 300 fs, the excess electronic energy rapidly generates modes with strong e-ph coupling, including the A 1 and E2g optical modes, respectively with wavevectors near the K and points of the BZ. A slower change in phonon populations occurs between 1 and 5 ps, when the optical phonons decay into LA and TA modes with wavevectors halfway between and the K or M points of the BZ. Long-wavelength (small-q) acoustic phonons are also generated in the same time window. A second stage of phonon equilibration sets in after ∼20 ps, when energy redistribution occurs primarily within the LA and TA branches (see Fig. 2).
We carry out a mode-by-mode analysis of the phonon populations and effective temperatures, taking into account mode crossing as discussed in Appendix B. In Fig. 3(a), we plot for each mode the excess population relative to equilibrium, where Nν (t ) are BZ-averaged phonon populations for each mode. Also shown in Fig. 3(a) are the mode-resolved e-ph coupling strengths, g 2 ν (see Appendix B). Due to their strong e-ph coupling, the in-plane longitudinal optical (LO) and transverse optical (TO) modes are rapidly emitted during the initial fast carrier cooling, highlighting the key role of e-ph interactions at early times [4].
The LO, TO, and LA modes are populated extensively before 1 ps, while the TA and ZO modes are excited more gradually through ph-ph interactions over timescales longer than 1 ps. Due to their weak e-ph coupling, the out-of-plane flexural phonon modes (ZA and ZO) are generated at a significantly slower rate than in-plane phonons, resulting in a slow rise of their populations in the first 10 ps. The ZA mode, with the weakest e-ph scattering strength, exhibits the slowest rise in population among all modes. In turn, the weak e-ph coupling and slow generation of flexural phonons is responsible for a carrier cooling bottleneck: As electrons and holes interact with hot ZA and ZO populations, they gain energy from the flexural modes for over 10 ps, well beyond the initial subpicosecond carrier cooling.
Although at short times the phonon distributions are still nonthermal, we find that after 2-5 ps all phonon populations are thermal and well approximated by hot Bose-Einstein distributions. The mode-resolved effective phonon temperatures are computed as BZ averages of state-dependent temperatures (see Appendix A) and shown in Fig. 3(b). We find that the effective phonon temperatures are strongly mode dependent throughout the simulation, their trend mirroring the respective FIG. 2. Phonon scattering channels. Optical phonons at and K are rapidly populated in the first 300 fs due to their strong coupling with the excited carriers, but later decay into lower-energy optical and acoustic modes through ph-ph processes. The dominant energy and momentum-conserving phonon decay channels are shown with red arrows at 5-and 20-ps delay times after the initial electronic excitation.

Non-equilibrium lattice dynamics from the time-dependent Boltzmann equation
To illustrate the application of the TDBE formalism to the ultrafast electron-phonon dynamics of 2D materials and its suitability for the description of these phenomena, we report in Fig. 9 ultrafast dynamics simulations for graphene obtained from the numerical solutions of Eqs. (27) and (28) and reproduced from Ref. [143]. In panels (a)-(c), the electron (f nk ) and hole (1-f nk ) distribution functions are superimposed to the band dispersion of graphene for energies above and below the Fermi level (dashed) for several time delays. The initial (t = 0) distribution of electrons and holes loses its excess energy via emitting phonon in the vicinity of the Γ and K high symmetry points and relaxes back to Fermi level within few picoseconds. The non-equilibrium phonon population is illustrated in Figs. 9 (d)-(f), where the point size is proportional to the phonon occupation n qν (t). The phonon population at 300 fs ( Fig. 9 (e)) is characterized by an enhancement of the number LO and TO phonons at Γ and K, indicating that these modes constitute the primary decay path for excited electrons for the initial phases of the relaxation. These findings are compatible with the Eliashberg function shown in Fig. 7 (b). On longer time scales, the phonon-phonon scattering leads to the thermalization of lattice and a redistribution to the excess energy of these modes to low-energy modes Fig. 9 (g). More generally, phonon-emission processes are constrained by energy and momentum conservation laws, which reduce the phase-space available for phonon-assisted electronic transitions. For example, graphene carriers in the vicinity of the Fermi level at K can scatter to states within the same Dirac cone at K, leading to the emission of phonons with q Γ, alternatively they can scatter to the second Dirac cone at −K, emitting phonons with crystal momentum around K or -K. Transitions to other region of the BZ are forbidden as they violate conservation laws. This stringent momentum selectivity confines the emitted phonons to narrow regions in reciprocal space, leading to the emergence of hot spots in the BZ, i.e., regions characterized by an enhanced vibrational temperature. Figure 10 (a) illustrates the BZ and high-symmetry points of monolayer MoS 2 , reproduced from Ref. [207]. The superimposed color coding denotes the effective vibrational temperature obtained by inverting the Bose-Einstein distribu- phonon populations (see Fig. 2). In the first 300 fs, the excess electronic energy rapidly generates modes with strong e-ph coupling, including the A 1 and E 2g optical modes, respectively with wavevectors near the K and points of the BZ. A slower change in phonon populations occurs between 1 and 5 ps, when the optical phonons decay into LA and TA modes with wavevectors halfway between and the K or M points of the BZ. Long-wavelength (small-q) acoustic phonons are also generated in the same time window. A second stage of phonon equilibration sets in after ∼20 ps, when energy redistribution occurs primarily within the LA and TA branches (see Fig. 2).
We carry out a mode-by-mode analysis of the phonon populations and effective temperatures, taking into account mode crossing as discussed in Appendix B. In Fig. 3(a), we plot for each mode the excess population relative to equilibrium, where N ν (t ) are BZ-averaged phonon populations for each mode. Also shown in Fig. 3(a) are the mode-resolved e-ph coupling strengths, g 2 ν (see Appendix B). Due to their strong e-ph coupling, the in-plane longitudinal optical (LO) and transverse optical (TO) modes are rapidly emitted during the initial fast carrier cooling, highlighting the key role of e-ph interactions at early times [4].
The LO, TO, and LA modes are populated extensively before 1 ps, while the TA and ZO modes are excited more gradually through ph-ph interactions over timescales longer than 1 ps. Due to their weak e-ph coupling, the out-of-plane flexural phonon modes (ZA and ZO) are generated at a significantly slower rate than in-plane phonons, resulting in a slow rise of their populations in the first 10 ps. The ZA mode, with the weakest e-ph scattering strength, exhibits the slowest rise in population among all modes. In turn, the weak e-ph coupling and slow generation of flexural phonons is responsible for a carrier cooling bottleneck: As electrons and holes interact with hot ZA and ZO populations, they gain energy from the flexural modes for over 10 ps, well beyond the initial subpicosecond carrier cooling.
Although at short times the phonon distributions are still nonthermal, we find that after 2-5 ps all phonon populations are thermal and well approximated by hot Bose-Einstein distributions. The mode-resolved effective phonon temperatures are computed as BZ averages of state-dependent temperatures (see Appendix A) and shown in Fig. 3(b). We find that the effective phonon temperatures are strongly mode dependent throughout the simulation, their trend mirroring the respective FIG. 2. Phonon scattering channels. Optical phonons at and K are rapidly populated in the first 300 fs due to their strong coupling with the excited carriers, but later decay into lower-energy optical and acoustic modes through ph-ph processes. The dominant energy and momentum-conserving phonon decay channels are shown with red arrows at 5-and 20-ps delay times after the initial electronic excitation. play only after 3 ps. Throughout the simulation, low-energy electron-hole pairs are excited near the Dirac cone via ph-e processes, resulting in a modest ultrasonic attenuation of the phonon dynamics. Overall, the excess energy initially imparted to the excited phonons mainly dissipates through slow ph-ph processes, the entire phonon relaxation taking tens of picoseconds. Our approach is uniquely able to shed light on these longer timescales characteristic of phonon dynamics.

III. DISCUSSION
A key advantage of our rt-BTE approach is the possibility, rather unique among existing first-principles methods for ultrafast dynamics, to validate the interactions against experiments, thus guaranteeing the quantitative accuracy of our simulated dynamics. We validate the e-ph and ph-ph interactions employed in our calculations, respectively, by computing electrical and heat transport properties of graphene (see Appendix A). We obtain a room-temperature electron mobility of 186000 cm 2 /V s, in excellent agreement with experimental results in suspended graphene [61,62]. We also compute the thermal conductivity in the single-mode relaxation time approximation (RTA) of the BTE, obtaining a room-temperature value of 482 W/mK in excellent agreement with previous RTA calculations [63]; this result implies that the full solution of the BTE (beyond the RTA) would give a thermal conductivity consistent with experiment [63]. These results show that the interactions employed in our ultrafast dynamics are precise, so our computed timescales are expected to be quantitatively accurate.
Finally, note that our method treats e-ph and ph-ph interactions perturbatively, using values for these interactions computed on the equilibrium structure. Therefore our approach is most suitable to model the low-excitation regime in which the atomic positions remain fairly close to their equilibrium values; for example, in our calculations the phonon  1. Simulated carrier and phonon dynamics in graphene. Excited electrons (red) and holes (purple) in graphene are initialized in a Fermi-Dirac distribution at 4000 K electronic temperature at time zero. The subsequent nonequilibrium dynamics is analyzed by mapping the populations f nk (t ) for electrons and 1 − f nk (t ) for holes on the electronic band structure near the Dirac cone (top) and the phonon populations N νq (t ) on the phonon dispersions (bottom), using a point size proportional to the population values. Shown left of each population panel are the corresponding Brillouin-zone-averaged carrier or phonon concentrations as a function of energy (in arbitrary units).
phonon populations (see Fig. 2). In the first 300 fs, the excess electronic energy rapidly generates modes with strong e-ph coupling, including the A 1 and E 2g optical modes, respectively with wavevectors near the K and points of the BZ. A slower change in phonon populations occurs between 1 and 5 ps, when the optical phonons decay into LA and TA modes with wavevectors halfway between and the K or M points of the BZ. Long-wavelength (small-q) acoustic phonons are also generated in the same time window. A second stage of phonon equilibration sets in after ∼20 ps, when energy redistribution occurs primarily within the LA and TA branches (see Fig. 2).
We carry out a mode-by-mode analysis of the phonon populations and effective temperatures, taking into account mode crossing as discussed in Appendix B. In Fig. 3(a), we plot for each mode the excess population relative to equilibrium, where N ν (t ) are BZ-averaged phonon populations for each mode. Also shown in Fig. 3(a) are the mode-resolved e-ph coupling strengths, g 2 ν (see Appendix B). Due to their strong e-ph coupling, the in-plane longitudinal optical (LO) and transverse optical (TO) modes are rapidly emitted during the initial fast carrier cooling, highlighting the key role of e-ph interactions at early times [4].
The LO, TO, and LA modes are populated extensively before 1 ps, while the TA and ZO modes are excited more gradually through ph-ph interactions over timescales longer than 1 ps. Due to their weak e-ph coupling, the out-of-plane flexural phonon modes (ZA and ZO) are generated at a significantly slower rate than in-plane phonons, resulting in a slow rise of their populations in the first 10 ps. The ZA mode, with the weakest e-ph scattering strength, exhibits the slowest rise in population among all modes. In turn, the weak e-ph coupling and slow generation of flexural phonons is responsible for a carrier cooling bottleneck: As electrons and holes interact with hot ZA and ZO populations, they gain energy from the flexural modes for over 10 ps, well beyond the initial subpicosecond carrier cooling.
Although at short times the phonon distributions are still nonthermal, we find that after 2-5 ps all phonon populations are thermal and well approximated by hot Bose-Einstein distributions. The mode-resolved effective phonon temperatures are computed as BZ averages of state-dependent temperatures (see Appendix A) and shown in Fig. 3(b). We find that the effective phonon temperatures are strongly mode dependent throughout the simulation, their trend mirroring the respective FIG. 2. Phonon scattering channels. Optical phonons at and K are rapidly populated in the first 300 fs due to their strong coupling with the excited carriers, but later decay into lower-energy optical and acoustic modes through ph-ph processes. The dominant energy and momentum-conserving phonon decay channels are shown with red arrows at 5-and 20-ps delay times after the initial electronic excitation.

023072-3
(g) , with N ph = 9 being the number of phonon modes of monolayer MoS 2 , for crystal momenta within the first BZ and for selected time steps. The same color bar (shown beside panel i) is used for panels a−i.
At t = 0 (Figure 3a), the lattice is at thermal equilibrium, as reflected by the constant vibrational temperature in the BZ (Tq = T ph 0 = 100 K). As the coupled electron−phonon dynamics begins (t > 0), the excited carriers in the valence and conduction bands relax back to Fermi level by transferring energy to the lattice through the emission of phonons. The influence of these processes on the phonon distribution function is accounted for by the phonon−electron collision integral (Γ qν pe ) in eq 4, which leads to an increase of n qν (and thus of T qν ) as phonons with matching crystal momenta q and index ν are emitted. After t = 100 fs (Figure 3b) the lattice has abandoned the initial thermalized state, as revealed by the emergence of hot-spots in the BZ characterized by a higher average vibrational temperature Tq. In particular, an increase of the vibrational temperature is observed for momenta close to Γ and K, which in turn reflects an enhancement of the phonon population.
To understand the origin of these features, it should be noted that the emission of phonons, and thus the change of T qν , is triggered by electronic transitions within the valence and conduction bands, which are heavily constrained by energy and momentum conservation laws. For the excited electronic distribution of Figure 2a, for instance, phonon-assisted transitions within the valence band primarily involve two types of processes: (i) intravalley transitions, connecting initial and valley also enables the emission of phonons around schematic illustration of the allowed inter-and phonon-assisted transitions is provided in Figure  Supporting Information. Umklapp processes are also these pictures, because transitions connecting differ be folded back to the first BZ via translation by lattice vector. This picture enables us to attribute the increase of vibrational temperature to the preferen of phonons at Γ and K, which is dictated by selectivity in the electronic transitions.
This mechanism leads to a further enhancem temperature anisotropy in the BZ for t = 500 fs Additionally, an increase in vibrational temperature at the M point and, less pronouncedly, at Q, whic transitions involving the Q pocket in the conductio longer time scales, phonon−phonon scattering, ac by the phonon−phonon collision integral (Γ qν pp counterbalances a nonthermal vibrational state by lattice toward a thermalized regime (namely, T qν This behavior is manifested for t = 1.5 and 3 ps (Figu progressive reduction of the temperature anisotrop The mode-and momentum-resolved vibrational T qν , superimposed to the phonon dispersion in Figu = 0.1, 0.5, and 3 ps, may further change significantly phonon branches, because the contribution of each the relaxation process is dictated by its own electr coupling strength. 32,44 In particular, the stronger optical phonons makes them a more likely decay ch At t = 0 (Figure 3a), the lattice is at thermal equilibrium, as reflected by the constant vibrational temperature in the BZ (Tq = T ph 0 = 100 K). As the coupled electron−phonon dynamics begins (t > 0), the excited carriers in the valence and conduction bands relax back to Fermi level by transferring energy to the lattice through the emission of phonons. The influence of these processes on the phonon distribution function is accounted for by the phonon−electron collision integral (Γ qν pe ) in eq 4, which leads to an increase of n qν (and thus of T qν ) as phonons with matching crystal momenta q and index ν are emitted. After t = 100 fs (Figure 3b) the lattice has abandoned the initial thermalized state, as revealed by the emergence of hot-spots in the BZ characterized by a higher average vibrational temperature Tq. In particular, an increase of the vibrational temperature is observed for momenta close to Γ and K, which in turn reflects an enhancement of the phonon population.
To understand the origin of these features, it should be noted that the emission of phonons, and thus the change of T qν , is triggered by electronic transitions within the valence and conduction bands, which are heavily constrained by energy and momentum conservation laws. For the excited electronic distribution of Figure 2a, for instance, phonon-assisted transitions within the valence band primarily involve two types of processes: (i) intravalley transitions, connecting initial and final electronic states both located close to the same highsymmetry point (Γ or K); (ii) intervalley transitions, with the initial and final electronic states located at Γ and K, respectively (or vice versa). Phonon-assisted transitions across the gap are forbidden by energy conservation. Because of momentum conservation, processes of type (i) result in the emission of longwavelength phonons (q ≃ 0) with momenta close to Γ, whereas processes of type (ii) can involve only the emission of phonons with momenta around K. A similar picture applies to transitions in the conduction band. Here, however, the presence of the Q valley also enables the emission of phonons a schematic illustration of the allowed inte phonon-assisted transitions is provided in Supporting Information. Umklapp processes a these pictures, because transitions connecting be folded back to the first BZ via translatio lattice vector. This picture enables us to attrib increase of vibrational temperature to the pre of phonons at Γ and K, which is dictate selectivity in the electronic transitions.
This mechanism leads to a further enh temperature anisotropy in the BZ for t = 5 Additionally, an increase in vibrational tempe at the M point and, less pronouncedly, at Q transitions involving the Q pocket in the con longer time scales, phonon−phonon scatter by the phonon−phonon collision integra counterbalances a nonthermal vibrational st lattice toward a thermalized regime (namely This behavior is manifested for t = 1.5 and 3 ps progressive reduction of the temperature ani The mode-and momentum-resolved vibra T qν , superimposed to the phonon dispersion i = 0.1, 0.5, and 3 ps, may further change signifi phonon branches, because the contribution o the relaxation process is dictated by its own coupling strength. 32,44 In particular, the str optical phonons makes them a more likely dec relaxation of excited carriers, as compared to This trend is reflected in Figures 3g,h by the temperature of these modes throughout the i dynamics, suggesting that the electronic co modes plays a primary role in the emergenc state of the lattice. A comprehensive picture of the formatio nonthermal vibrational state is provided in illustrates the time evolution of the av temperature Tq for momenta along the M-  Here, the phonon distribution function n qν (t) is obtained from the solution of the coupled TDBE for electrons and lattice. The preferential emission of phonons around Γ, K, and -K is reflected by the enhanced vibrational temperature at these high-symmetry points for the initial phases of hot-carrier relaxation. Phonon hot-spots in the BZ can persist for several picoseconds, until the onset of phonon-phonon scattering processes restores a regime of thermal equilibrium.

Summary and Outlook
In this manuscript we introduced parameter-free theoretical approaches for the study of the coupled ultrafast dynamics of electrons and phonons in solids, and we reviewed their application to the prototypical two-dimensional materials, such as graphene and MoS 2 .
The two-temperature and the non-thermal lattice models recast the thermalization dynamics of electrons and phonons into a set of coupled first-order differential equations, which is formally equivalent to the one governing the temperature evolution of coupled thermal baths. While on the one hand these approximations are sufficient to capture the characteristic timescales and dynamics of the electron thermalization with the lattice revealed by pump-probe experiments, on the other hand, the nonequilibrium electron, vibrational and structural dynamics are poorly described, restricting the domain of applicability of these approaches. In particular, the emergence of non-equilibrium states of the lattice often entails phonon populations characterized by strong anisotropy in the Brillouin zone that cannot be captured by a Bose-Einstein function.
The time-dependent Boltzmann equation overcomes these limitations. As compared to models, the advantage of this approach is two-fold: (i) it provides a description of ultrafast electron and phonon dynamics with full momentum resolution (ii) it accounts for phase-space constraints in phonon-assisted electronic transitions. These points are important prerequisites to capture the emergence of phonon hot-spots in the Brillouin zone, which have recently been revealed by ultrafast diffuse scattering experiments [38,39,83,145,152]. Additionally, the collision integrals can be made fully parameter-free by determining electron-phonon and phonon-phonon coupling matrix elements via available first-principles electronic-structure codes [138,157,208].
Fostered by the relentless advances in the experimental characterization of ultrafast phenomena, first-principles approaches for coupled electron-phonon dynamics constitute a rapidly evolving research field. Many non-equilibrium phenomena, however, still represent de facto a challenge for ab-initio techniques. These include, for instance, coherence and dephasing of electronic and vibrational excitations [209], correlation effects on the electron dynamics [210], structural phase transitions [211][212][213], topological Floquet states [214], and the transient formation of quasiparticle states [215]. The development of accurate and efficient first-principles techniques capable of tackling these phenomena remains a key research priority in the field of ultrafast science. In this respect, approaches based on the density-matrix formalism [194] and non-equilibrium Green's functions [210] are promising routes to overcome the inherent limitations of the Boltzmann equation: recent studies explored new directions to reduce the notoriously high computational cost of NEGF methods [216] paving the way towards a rigorous quantum mechanical description of coherent and correlated dynamical phenomena in condensed matter [217].
Below we report the analytical solution of the TTM for the limiting case of temperature-independent heat capacities. The TTM equations in absence of external driving field are reported below: The TTM thus forms a set of two coupled first-order differential equation, which can be solved analytically via standard techniques. Its solution is: T ph (t) = − C el a 1 C ph e −γt + a 2 , where we introduced the constants a 1 = (T 0 el − T 0 ph )C ph /(C ph + C el ), a 2 = T 0 el − c 1 , and γ = g/(1/C el + 1/C ph ). The solution is easily verified by substitution.
By introducing a source term coupled to the electronic temperature, the TTM gets 21 modified as follows: where S(t) describes the interaction with a light pulse, and the second equation remains unchanged. Analytic solution of the TTM is still possible for a simple time dependence of S(t). For example, by considering S(t) = αe − t τ θ(t), the solution can be written as: T ph (t) = T el (t) + with the constants: