Toward precise simulations of the coupled ultrafast dynamics of electrons and atomic vibrations in materials

Ultrafast spectroscopies can access the dynamics of electrons and nuclei at short timescales, shedding light on nonequilibrium phenomena in materials. However, development of accurate calculations to interpret these experiments has lagged behind as widely adopted simulation schemes are limited to subpicosecond timescales or employ simpliﬁed interactions lacking quantitative accuracy. Here we show a precise approach to obtain the time-dependent populations of nonequilibrium electrons and atomic vibrations (phonons) up to tens of picoseconds, with a femtosecond time resolution. Combining ﬁrst-principles electron-phonon and phonon-phonon interactions with a parallel numerical scheme to time-step the coupled electron and phonon Boltzmann equations, our method provides microscopic insight into scattering mechanisms in excited materials. Focusing on graphene as a case study, we demonstrate calculations of ultrafast electron and phonon dynamics, transient optical absorption, structural snapshots, and diffuse x-ray scattering. Our ﬁrst-principles approach paves the way for quantitative atomistic simulations of ultrafast dynamics in materials.


I. INTRODUCTION
Time-domain spectroscopies open a window on the dynamics of electrons, phonons, and various elementary excitations, probing matter on timescales characteristic of the interactions among its constituents. In particular, recently developed experimental techniques enable the characterization of the coupled dynamics of electronic and nuclear degrees of freedom with high spatial and time resolutions [1][2][3][4], as well as manipulation of ultrafast structural dynamics [5][6][7][8][9][10] and phase transitions [11][12][13].
The vast information encoded in ultrafast spectroscopy signals underscores a critical need for quantitative theoretical tools. The latter should ideally be able to predict or interpret ultrafast measurements and shed light on the underlying microscopic processes and coupling between different degrees of freedom. A common scenario is the coupled dynamics of excited electrons and phonons, which governs a wide range of phenomena such as carrier and lattice equilibration, photoemission and the associated linewidths, structural transitions, and superconductivity. Heuristic approaches such as twotemperature models (TTMs) [14] or kinetic equations with adjustable interactions [15] are routinely adopted to study aspects of nonequilibrium electron and phonon dynamics. In particular, TTMs can successfully model the dynamics of materials excited above the damage threshold [16,17] and treat electron interactions in that regime in conjunction with molecular dynamics [18,19]. However, at lower excitations and in the context of modeling pump-probe experiments, TTMs and related approaches are not geared toward quantitative predictions [20] and are mainly useful as tools to qualitatively interpret experimental results.
First-principles computational methods based on density functional theory (DFT) and related techniques have made great strides in modeling electron and phonon interactions and nonequilibrium dynamics [21][22][23][24][25][26]. An important approach is real-time time-dependent DFT (rt-TDDFT) [27][28][29], which propagates in time the electronic Kohn-Sham equations while treating atomic motions using Ehrenfest forces. While rt-TDDFT has been widely successful for studying electron dynamics in molecules [30,31], important open challenges remain, including reaching simulation times longer than ∼1 ps, extracting phonon-related quantities when treating periodic systems (crystals) [32], and better understanding, improving, and validating the accuracy of exchange-correlation functionals governing the interactions of electrons and nuclei [33][34][35]. Nonadiabatic molecular dynamics is also a valuable method to model ultrafast dynamics in molecules [36] and solids [37].
A second family of approaches-on which this work hinges-employs perturbation theory to compute the matrix elements for electron and phonon interactions [38], and combines them with the real-time Boltzmann transport equation (rt-BTE) [24,39] or quantum kinetic equations [40] to investigate ultrafast electron dynamics. Owing to its momentum-space formulation, this framework is ideally suited for studies of crystals and periodic systems. However, in spite of recent progress, propagating in time the Boltzmann transport equations (BTEs) for coupled electrons and phonons while fully taking into account their interactions has remained an elusive goal due to computational cost and the need for complex algorithms and workflows.
Here we show explicit simulations of the coupled dynamics of nonequilibrium electrons and phonons for timescales up to tens of picoseconds, using accurate electron-phonon (e-ph) and phonon-phonon (ph-ph) interactions validated through transport calculations. We time-step the electron and phonon BTEs-a large set of coupled integro-differential equationsusing a parallel numerical approach to efficiently compute the relevant collision integrals and obtain time-dependent electron and phonon populations. Focusing on graphene, a material with unconventional physical properties and unique promise for ultrafast devices, we investigate the coupled dynamics of excited electrons and phonons, compute spectroscopic signals such as transient absorption and diffuse x-ray scattering, and analyze the dominant scattering channels, identifying the slow rise of flexural phonons as a bottleneck to equilibration. Taken together, our work demonstrates a leap forward in simulations of ultrafast dynamics in solids and provides a quantitative tool to predict and interpret time-domain spectroscopies.

A. Numerical approach
We solve the coupled electron and phonon BTEs [41] in a homogeneously excited region of a material using a fourthorder Runge-Kutta algorithm. At each time t, we propagate the electronic populations f nk (t ), where n is the band index and k the crystal momentum, and the phonon populations N νq (t ), where ν is the mode index and q the phonon wavevector. Using an in-house modified version of our PERTURBO code [42], we solve the BTEs as a set of integro-differential equations: where I e-ph is the collision integral for electron-phonon (eph) interactions in the electron BTE, while I ph-e and I ph-ph are, respectively, the collision integrals for phonon-electron (ph-e) and phonon-phonon (ph-ph) interactions in the phonon BTE. The collision integrals are obtained using e-ph matrix elements computed in equilibrium at 0 K and ph-ph matrix elements computed in equilibrium at 300 K. Detailed expressions for these integrals are given in Appendix B.
The collision integrals depend explicitly on the electron and phonon populations at the current time step, and thus Eq. (1) is a system of integro-differential equations, with size N e = N b × N k for electrons and N ph = N ν × N q for phonons (N b and N k are the number of electronic bands and k-points, while N ν and N q are the number of phonon modes and q-points). A typical calculation with our scheme involves time-stepping ∼10 6 coupled integro-differential equations, each with 10 7 -10 9 scattering processes in the collision integrals, a formidable computational challenge addressed with an efficient algorithm combining message-passing interface (MPI) and multi-threading (OpenMP) parallelizations. As memory effects are not included, our BTEs amount to a Markovian dynamics of the density matrix in which offdiagonal coherences, which are typically short lived (< 5 fs), are neglected [15]. Using this scheme, we can reach simulation times up to ∼100 ps, with a femtosecond time step.

B. Relaxation of excited carriers
In a first set of simulations, we model photoexcited electron and hole carriers in graphene by choosing appropriate carrier populations at time zero as the initial condition of the BTEs. Pump-probe studies on graphene have shown that fast electron-electron interactions rapidly thermalize the excited carriers within 10-50 fs [43,44]. Therefore, shortly after excitation, electrons and holes are well described by a hot Fermi-Dirac distribution with a temperature of a few thousand degrees for typical experimental settings [45][46][47]. We model the carrier dynamics in graphene after this initial thermalization by setting the electron and hole populations at time zero to Fermi-Dirac distributions at 4000 K temperature, while setting the initial phonon populations to their equilibrium value at 300 K. Figure 1 shows the time evolution of the electron, hole, and phonon populations in the first picosecond after the initial excitation and carrier thermalization. For representative times, we analyze both the carrier and phonon populations ( f nk for electrons, 1 − f nk for holes, and N νq for phonons) along highsymmetry Brillouin zone (BZ) lines, as well as BZ-averaged carrier and phonon concentrations as a function of energy (see Appendix B). In the first 100 fs, the hot electron and hole distributions become narrower in energy as the carriers cool to the Dirac point in about 300 fs. This rapid cooling process is achieved by emitting optical phonons [48] through intravalley and intervalley e-ph scattering. As a result, the optical phonons generated within 1 ps possess wavevectors close to the BZ point for intravalley and K point for intervalley scattering.
Although most of the energy stored in the excited carriers is dissipated in the first 300 fs, surprisingly it takes more than 5 ps for the carriers to fully equilibrate, with negligible population fluctuations. This slow and lingering cooling process is due to two factors-one is the slow rise of out-of-plane flexural phonons due to their weak e-ph coupling (see below), and the other is the weak anharmonic ph-ph coupling allowing optical phonons generated during the rapid carrier cooling to linger for > 10 ps, heating back the carriers via phonon absorption.

C. Phonon equilibration
Accessing the 10-100 ps timescale characteristic of phonon dynamics is an open challenge for existing firstprinciples simulations of ultrafast dynamics. In Fig. 2, we analyze the long-lived equilibration of optical phonons generated by the excited carriers. The optical phonons first decay within ∼5 ps to transverse acoustic (TA) and longitudinal acoustic (LA) phonons, which then emit lower-energy acoustic modes through ph-ph scattering processes, ultimately equilibrating on a 20-100 ps timescale (not shown).
Our calculations can directly identify the main phonon modes and relaxation pathways from the time-dependent 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, 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. excess phonon populations in Fig. 3(a). This result shows that a two-temperature model would fail dramatically in graphene, both because the phonon temperatures are not well defined before 2 ps and because they are mode dependent at longer times.
Our calculations shed light on the entire equilibration cascade, with timescales spanning several orders of magnitudes, from femtosecond for carrier cooling through e-ph interactions to picosecond for phonon down-conversion to equilibration via ph-ph processes over tens of picoseconds. The microscopic details of the scattering processes, accessed easily in our approach due to its formulation in momentum space, are out of reach for existing first-principles simulations. As shown next, the time-dependent electron and phonon populations further allow us to reliably simulate various ultrafast spectroscopies employed as probes of electron and nuclear dynamics.

D. Pump-probe transient absorption
Time-resolved differential transmission measurements are widely employed to shed light on nonequilibrium carrier dynamics and investigate the timescales associated with electronic processes. We build on the methods shown above to compute transient absorption at optical wavelengths. In graphene, at low energies within 1-2 eV of the Dirac cone, the differential transmissivity T /T 0 is proportional to the transient carrier populations; in a single-particle picture and neglecting the state dependence of the optical transition dipoles [49], where E is the carrier energy, a 0 is the absorption coefficient of single-layer graphene, f e,h are BZ-averaged electron or hole populations, and f e,h ( are the corresponding excess carrier populations relative to equilibrium. We carried out calculations of T /T 0 and compared them to two sets of experimental results. In the first simulation, we set the initial electron temperature to 3600 K and probe the T /T 0 over the energy window of relevance (0.6-0.9 eV) in the experiment by Breusing et al. [49]. In the second simulation, the initial electron temperature is set to 4000 K and T /T 0 is computed at 900-nm wavelength to compare with experiments by Brida et al. [50]. For both simulations, the initial phonon temperature is 300 K and the calculated T /T 0 is shown as a function of pump-probe delay time in Fig. 4(a).
In both cases, we obtain very good agreement between theory and experiment.
The experimental result exhibits two distinct temporal regions, which are more clearly seen in the experiments by Brida et al. [50] [red data points in Fig. 4(a)]. At sub-10-fs time delay, T /T 0 increases due to the coupled effects of photoexcitation and electron-electron (e-e) interactions, both of which are not treated explicitly in our calculations but are taken into account in our choice of an initial hot Fermi-Dirac distribution. The second time window at t > 10 fs, modeled here explicitly, shows a slow time decay of T /T 0 with a time constant of a few picoseconds. Our simulated differential transmission spectrum agrees well with experiment [50] in this time window: following a faster decay of T /T 0 due to carrier cooling in the first 300 fs, we correctly predict the onset of a slower cooling of the excited carriers associated with the bottleneck from slowly rising flexural phonons discussed above. While previous work focused on the role of e-e processes [50], our results show the importance of e-ph interactions at subpicosecond times.

E. Diffuse scattering and structural dynamics
Intense experimental efforts are focused on investigating atomic motions and structural dynamics in the time domain [51]. Among other techniques, ultrafast electron [52] and x-ray diffraction [53] have made great strides and are currently able to infer atomic displacement patterns and dominant phonon modes governing the nonequilibrium structural dynamics [51,54]. Using the time-resolved phonon populations for each vibrational mode and wavevector, we can simulate and reconstruct the structural dynamics fully from first principles.
In a classical picture, the time-dependent displacement of atom a can be decomposed into lattice waves with wavevector q and mode ν [53], where A νq is the wave amplitude, a νq its polarization vector projected on atom a, and δ νq a phase factor. For quantized lattice vibrations (phonons), νq is the phonon eigenvector obtained by diagonalizing the dynamical matrix and the amplitude A νq can be expressed in terms of phonon populations N νq as where m a is the atomic mass. The phase factors δ νq are chosen randomly at time zero and do not affect the dynamics. The computed classical atomic displacements u a at time delays of 0.5 and 5 ps are visualized in Fig. 4(b). At 0.5 ps, the excitation of optical phonons from carrier cooling drives a dominant unidirectional vibration of atoms. At longer times, as shown in the 5-ps panel, a larger number of phonon modes are incoherently excited through ph-ph processes and the atomic vibrations become randomized, consistent with the thermal phonon distributions we find at 5 ps. The quantity measured experimentally in x-ray diffraction is the averaged square atomic displacement, u 2 a = νq 1 2 A 2 νq a , which we show in the Supplemental Material [55], where the radius of the blue sphere around each atom represents u 2 a . The averaged atomic displacements increase monotonically with time, reaching greater values as the system approaches thermal equilibrium.
Thermal diffuse x-ray scattering has proven successful for extracting phonon dispersions and investigating momentumdependent phonon dynamics [56,57]. However, interpreting measurements of the diffuse scattering intensity I (q, t ) is challenging due to its rich information content. For example, mode-resolved phonon contributions cannot be recovered straightforwardly from the measured signal, hindering the identification of the dominant phonon modes governing ultrafast structural dynamics and anharmonic phonon processes. Here we demonstrate the inverse process of reconstructing the experimental diffuse scattering signal from the computed mode-resolved phonon populations.
We write the scattering intensity at wavevector q as [53] where M is the Debye-Waller factor, and compute the ultrafast diffuse scattering signal as the difference between the time-dependent and equilibrium scattering intensities, I (q, t )/I = I (q,t )−I (q,−∞) . Our computed diffuse scattering results for graphene are plotted in reciprocal space in Fig. 4(c) at 0.5-and 5-ps time delays. At 0.5 ps, we find peaks of I/I near the BZ corners due to optical phonons generated from intervalley scattering during the initial subpicosecond carrier cooling. Small-q optical phonons are also generated extensively, as discussed above, but are not visible in Fig. 4(c) as the denominator in the differential diffuse scattering I/I [the thermal equilibrium intensity I (q, −∞)] is large at small q due to the large equilibrium population of acoustic phonons. The diffuse scattering result at 5 ps captures the dissipation of the optical phonons' excess energy to the lower acoustic branches, whereby phonons are generated throughout the BZ. While the measured diffuse scattering signal cannot provide information on the mode-resolved dynamics, our simulations can fill this gap, providing detailed information on modedependent scattering processes in momentum space and their contributions to diffuse scattering.

F. Ultrafast dynamics following phonon excitation
Coherent excitation of selected phonon modes has become an important approach for investigating and manipulating material properties [58][59][60]. To demonstrate the flexibility of our rt-BTE numerical framework, we carry out a simulated experiment in which phonons are excited at time zero while electrons and holes are initially kept in their equilibrium room-temperature distributions. After populating a large excess of LO phonons near the BZ center at time zero, we simulate the subsequent phonon dynamics by solving the coupled electron and phonon BTEs, as above but for different initial conditions. The computed time-domain excess phonon populations N (E , t ) are plotted in Fig. 5. After 1 ps, LA phonon modes with ∼100 meV energy are generated through ph-ph processes, and within 2 ps both LA and TA modes with 80-110 meV energy absorb most of the excess energy from the optical modes. The flexural ZA phonons come into 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 above the initial 300 K temperature. Regimes involving intense excitation or phase transitions, atomic diffusion, or material damage require different treatments. Future work will address explicitly the issue of updating the e-ph and ph-ph interactions during the nonequilibrium dynamics.

IV. CONCLUSION
In summary, we have shown a versatile numerical framework, rt-BTE, for modeling the ultrafast coupled dynamics of electrons and phonons. We demonstrated its accuracy and broad applicability through simulations of pump-probe spectroscopy, x-ray diffuse scattering, and structural and phonon dynamics. Our results shed light on e-ph and ph-ph scattering processes governing ultrafast dynamics in graphene, and provide valuable information for the design of ultrafast electronic and optical devices. We plan to make the numerical method available in a future release of PERTURBO [42] to equip the community with a tool for simulating and interpreting ultrafast time-domain experiments. Future extensions will aim to treat explicitly the light excitation pulse and the electron spin to unravel the intertwined nonequilibrium dynamics of electronic, structural, and spin degrees of freedom. Taken together, our first-principles approach demonstrates a paradigm shift in computing the ultrafast electron and atomic vibrational dynamics, bridging the gap between theory and experiment and enabling quantitative predictions of ultrafast phenomena in materials.

ACKNOWLEDGMENTS
The authors thank Jin-Jian Zhou for fruitful discussions.

Density functional theory
We carry out first-principles DFT calculations on graphene with a relaxed in-plane lattice constant of 2.45 Å. The graphene sheet is separated from its periodic replicas by a 9-Å vacuum. The ground-state electronic structure is computed using the QUANTUM ESPRESSO code in the local density approximation (LDA) of DFT. We use a norm-conserving pseudopotential, a 90-Ry plane-wave kinetic energy cutoff, and a Methfessel-Paxton smearing of 0.02 Ry. The groundstate charge density is obtained using a 36 × 36 × 1 k-point grid, following which a non-self-consistent calculation is employed to obtain the Kohn-Sham eigenvalues and wave functions on a 12 × 12 × 1 k-point grid. To construct maximally localized Wannier functions (WFs) with the WANNIER90 code [64], the Kohn-Sham wave functions are first projected onto atomic p z orbitals on each atom and sp 2 orbitals on every other atom [65,66], for a total of 5 wannierized bands. The WF spread is then minimized, and the relevant energy windows are adjusted until the interpolated band structure can smoothly reproduce the LDA result within ∼10 meV throughout the BZ.

Electron-phonon scattering and correction to the phonon dispersions
We use density functional perturbation theory (DFPT) [67] to compute lattice dynamical properties and the e-ph perturbation potentials, and then form the e-ph matrix elements g nn ν (k, q) on coarse 12 × 12 × 1 k-point and q-point BZ grids using PERTURBO [42]. Here and below, the e-ph matrix elements g nn ν (k, q) quantify the probability amplitude for an electron in Bloch state |nk with energy E nk to scatter into a final state |n k + q with energy E n k+q due to emission or absorption of a phonon with branch index ν, wavevector q, and energyhω νq . The electron and phonon energies and the e-ph matrix elements are interpolated on fine grids using WFs with PERTURBO [42]. The average e-ph coupling strengths g 2 ν used in Fig. 3(a) are obtained as g 2 ν = nn kq |g nn ν (k, q)| 2 , where the summation includes electronic states near the Dirac point and phonon wavevectors in the entire BZ.
To account for the Kohn anomaly near q = K [68], we interpolate E nk and g nn ν (k, q) on a random BZ grid of 10 4 points, and obtain graphene phonon dispersions ω νq that include a previously proposed GW correction [69], where π (π * ) labels the occupied (empty) π band, N k = 10 4 is the total number of k-points in the random grid, and B GW q is a parameter employed to converge the LO and TO phonon energies at K. The rescaling factor γ GW q is computed using [70] γ GW with a 0 the graphene lattice constant, γ GW a constant equal to 1.61, and K n the nearest vector to q among those equivalent to K.
To verify the convergence of the q-point grid and compute the electrical mobility, we use PERTURBO [42] to calculate the state-dependent e-ph scattering rates nk within lowest-order perturbation theory [25], where f nk and N νq are the electron and phonon equilibrium occupations at 300 K. Convergence is achieved for a grid of 420 × 420 × 1 q-points (see the Supplemental Material [55]) using a Gaussian broadening of 20 meV to approximate the δ functions in Eq. (A3). The same grid is employed for all ultrafast dynamics calculations.

Phonon-phonon scattering
We use the temperature-dependent effective potential (TDEP) method [71][72][73] to obtain the interatomic force constants and ph-ph matrix elements νν ν (q, q ). The latter are the probability amplitudes for three-phonon scattering processes in which the phonon state |νq with energyhω νq scatters to the states |ν q + q and |ν q , or the inverse process in which the last two phonons combine to generate |νq . For the TDEP calculations, we prepare a number of 12 × 12 × 1 (288-atom) supercells with random thermal displacements corresponding to a canonical ensemble at 300 K. To validate the ph-ph interactions and compute the thermal conductivity, we compute the phonon scattering rates νq as [74] using phonon equilibrium occupations at 300 K. We compute and converge the ph-ph scattering rates (see the Supplemental Material [55]) with an in-house version of PERTURBO [42]. The ph-ph scattering rates converge for grids equal to or finer than 210 × 210 × 1 q-points using a 1-meV Gaussian broadening to approximate the δ functions in Eq. (A4).

Transport properties
The electron mobility is computed from a full solution of the linearized electron BTE with the PERTURBO code [42]. The thermal conductivity is computed using ph-ph matrix elements from TDEP with an in-house implementation of the single-mode RTA formula [63]. Both calculations use the converged BZ grids given above.

APPENDIX B: COUPLED ELECTRON AND PHONON ULTRAFAST DYNAMICS
We simulate the ultrafast dynamics of coupled electrons and phonons due to e-ph and ph-ph scattering processes in graphene in the absence of external fields. The time evolution of the carrier and phonon distributions is obtained by solving the coupled carrier and phonon BTEs: which account for e-ph scattering for electrons, e-ph scattering for phonons (denoted above as ph-e scattering), and ph-ph scattering for phonons. In Eq. (B1), f nk (t ) are time-dependent electron populations and N νq (t ) are time-dependent phonon populations. The e-ph matrix elements g nn ν (k, q) and ph-ph matrix elements νν ν (q , q ) are computed using DFPT in the ground state plus WF interpolation and TDEP, respectively, as specified above. We numerically solve Eq. (B1) using the fourth-order Runge-Kutta method with a time step of 2 fs, using uniform fine BZ grids of 420 × 420 × 1 kand q-points for all e-ph and ph-ph scattering processes. The ph-ph scattering matrix elements are computed on a uniform 210 × 210 × 1 q-point BZ grid and then Fourier interpolated to a 420 × 420 × 1 qpoint grid. These grids are sufficient to converge the scattering processes and ultrafast dynamics. To reduce the number of scattering processes entering the BTEs, we use only electronic states in the energy window of relevance for our calculations, which restricts the number of k-points to a few thousand, while the q-points span the entire BZ. We additionally select the relevant ph-ph processes by imposing energy conservation (within a few times the broadening value), thereby dramatically reducing the number of ph-ph scattering processes, which would otherwise be unmanageable for the fine q-point grid employed. The BZ-averaged energy-dependent carrier populations f (E , t ) and phonon populations N (E , t ) at energy E , first used in Fig. 1, are obtained respectively as f (E , t ) = nk f nk (t )δ(E nk − −E ) and N (E , t ) = νq N νq (t )δ(hω νq − −E ) via tetrahedron integration. Mode crossing has been taken into account when computing all mode-dependent quantities in this work. To reliably label the phonon modes, we ensure that the first and second derivatives of the phonon dispersions are smooth for each phonon mode, and label the modes according to their character at the BZ center.
The average temperature of each phonon mode at time t (see Fig. 3) is computed as T ν (t ) = q T νq (t ), where the state-dependent temperature T νq (t ) is obtained by inverting the Bose-Einstein occupation formula, N νq (t ) = (e¯h ω νq /k B T νq (t ) − 1) −1 , using the phonon populations at each time step. The phonon populations follow a Bose-Einstein thermal distribution at the computed temperature T ν (t ) for times greater than 2 ps for all acoustic modes and the ZO mode, and at times greater than 4-5 ps for the LO and TO optical modes.
We develop a scheme combining MPI and OpenMP parallelization to efficiently compute the e-ph, ph-e, and ph-ph collision integrals at each time step. Briefly, (k, q) pairs for e-ph scattering and (q, q , q ) triplets for ph-ph scattering are distributed among MPI processes; using OpenMP parallelization, rapid on-node operations are employed to compute the collision integrals within each MPI process; the results are then added together to form the collision integrals. This scheme leverages an algorithm we previously employed to time-step the electron BTE [24,42]; however, including ph-ph scattering is dramatically more difficult due to the need to take into account all q-point triplets involved in the ph-ph processes. To accomplish this task, we developed a parallel implementation of the phonon BTE in this work.