{'fraction_non_alphanumeric': 0.04291260098625538, 'fraction_numerical': 0.013254992480677089, 'mean_word_length': 4.290286771507863, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real-Time Excited Carrier Dynamics in Momentum Space from ab initio Nonadiabatic Molecular Dynamics

The application of nonadiabatic molecular dynamics (NAMD) approach is severely limited to studying the carrier dynamics in momentum space since a supercell is required to sample the phonon excitation and electron-phonon ( e-ph ) interaction at diﬀerent momenta in a molecular dynamics simulation. Here we develop an ab initio NAMD approach for the real-time carrier dynamics in momentum space (NAMD k ) by directly introducing e-ph coupling into the Hamiltonian based on harmonic approximation. The application of NAMD k on the hot carrier dynamics in graphene reveals that an energy threshold of 0.2 eV, which is deﬁned by two strongly coupled optical phonon modes, separates the hot electron relaxation into fast intervalley scattering and slow intravalley scattering regions with lifetimes of pico- and nano-second, respectively. Our results reveal the phonon-speciﬁc hot carrier relaxation mechanism in graphene, and the NAMD k approach provides a powerful tool to understand the real-time carrier dynamics in momentum space for diﬀerent materials.


Introduction
Tracking the quantum dynamics of excited charge carriers in solid materials in multi-dimensions including time and energy domains, and real and momentum spaces is fundamental to the understanding of many dynamical processes in optoelectronics, spin-and valley-tronics, solar energy conversion, and so on. [1][2][3][4][5][6] Different time-resolved experimental techniques including ultrafast time-and angle-resolved photoemission spectroscopy (TR-ARPES) with time, energy, and momentum resolution have been rapidly developed and applied to investigate charge carrier dynamics in various materials. However, without the input of ab initio investigations, it is rather difficult to understand the physical mechanisms behind the experimental spectra. Thus, it is urgent to develop an ab initio simulation approach to achieve a state-of-the-art understanding of multi-dimensional carrier dynamics in solids.
Ab initio approach based on perturbation theory provides useful information to understand carrier lifetimes governed by different scattering mechanisms. However, the real-time dynamics information cannot be obtained straightforwardly. In recent decades, the real-time time-dependent density functional theory (TDDFT) based on Ehrenfest dynamics [7,8] and nonadiabatic molecular dynamics (NAMD) on account of time-dependent Kohn-Sham equation (TDKS) and surface hopping have been applied to investigate the quantum dynamics of excited charge carrier. [6,9,10] The electron-phonon (e-ph) coupling, spin-orbital coupling (SOC), and many-body electron-hole interaction have been included with different theoretical strategies. [6,11,12] However, in all these methods, it is rather difficult to achieve real-time carrier dynamics in momentum space. In previous NAMD simulation, since the phonon excitation is simulated using ab initio molecular dynamics (AIMD) within periodic boundary conditions, only phonons at the Γ point are included in the simulation. Thus, the electron transition from one k to another through e-ph is forbidden even though the electronic states can be simulated using multi-k grids. To sample the phonon excitation and electron-phonon interaction at different momenta, a supercell needs to be used so that the phonons at other q-points can be folded to the Γ point. The q-grid density is determined by the size of the supercell, and due to the computational cost, usually, only a few q-points can be included in the NAMD simulation. As a contrast, the e-ph scattering between different momenta usually needs to be simulated with a very dense of k-and q-grids, especially when the electronic band dispersion is strong. Therefore, an ab initio approach to describe the real-time quantum dynamics of the photoexcited carrier in momentum space is essential.
In this work, by introducing the e-ph coupling elements into the timepropagation Hamiltonian, we have extended the ab initio NAMD approach from the real-space (labeled as NAMD r) to the momentum space (NAMD k). Different from the previous NAMD r approach, where an AIMD simulation using a large supercell is required for q-grid sampling, in NAMD k the k and q sampling is performed by the calculation of e-ph matrix elements using primitive cell. The computational cost is significantly reduced. Moreover, the NAMD k approach provides a straightforward picture not only for the dynamics of excited electrons in momentum space, but also the timedependent phonon excitation of the lattice due to the e-ph scattering. Using this approach, we have investigated the hot carrier dynamics in graphene. It is found that there is an energy threshold at 0.2 eV above the Fermi level (E f ), which separate the hierarchical relaxation dynamics from fast [picosecond(ps)] to slow [nanosecond(ns)] regions, which are dominated by the intervalley and intravalley e-ph scattering, respectively. The energy threshold is determined by strongly coupled optical phonon modes A 1 and E 2g . Our work not only reveals the phonon mode-specific energy threshold for hot electron relaxation in graphene, it also provides a powerful tool which can be widely applied to different solid system for the excited carrier dynamics in momentum space.

Methodology
In the NAMD r approach, the charge carrier (electron or hole) wavefunction is expanded on the basis of instantaneous adiabatic Kohn-Sham (KS) orbitals, which are obtained by solving the KS equation at atomic configuration R(t), Based on the classical path approximation (CPA), R(t) can be obtained by AIMD. The charge carrier wavefunction follows the time-dependent Schrodinger equation (TDSE) i ∂ ∂t |Ψ(r; R(t)) =Ĥ el (r; R(t)) |Ψ(r; R(t)) .
(2) Article Then a set of differential equations for the coefficients c m (t) is produced where n is the energy of the adiabatic KS state and d mn is the NAC between KS states m and n. The NAC can be written as Here m and n are the eigenvalues of the KS orbitals m and n, d mn is the e-ph coupling term, andṘ are the nuclear velocity. NAC is the crucial term in the NAMD r approach, which determines not only the time-dependent coefficient evolution but also the hopping probability in the subsequent surface hopping step. [9,10] According to Bloch's theory, if ψ m and ψ n have different k vectors, the NAC is zero (see more details in the supplementary materials). The essential reason is that in the AIMD simulation with periodic conditions, only phonons modes at the Γ point are included. Therefore the NAMD r approach can not sufficiently simulate the carrier dynamics in momentum space. The NAMD k approach is based on the harmonic approximation. In the NAMD k approach, different with eq. (1) we expand the charge carrier wavefunction using the KS orbitals of the equilibrium atomic configuration R 0 where the KS orbital |ψ nk (r; R 0 ) with band index n and momentum k is the eigenstate of the equilibrium configuration R 0 . Naturally, the charge carrier Hamiltonian is divided into two partŝ H el (r; R(t)) =Ĥ 0 (r; R 0 ) + ∆V (r; R(t)), where ∆V is the variation of potential induced by nuclear displacements ∆R(t) = R(t) − R 0 . Combining the above equations, we get a new coefficient evolution equation here H 0 mk ,nk = ψ mk |Ĥ 0 |ψ nk = nk δ mn,k k , is the diagonal KS energy matrix and is the e-ph coupling Hamiltonian. Here m and n are the notation of the KS orbitals, and k and k are the notation of the momentum. As transforming into momentum space, the e-ph term can be rewritten as where N p is the number of unit cell according to Born-von Kámánn boundary conditions, Q qν (t) is the normal mode coordinate of corresponding vibration mode of phonon with momentum q in the branch ν, l qν is the "zero-point" displacement amplitude, and g mnν (k, q) is e-ph matrix element. One can that in this way, the NAC in NAMD r approach is replaced by e-ph coupling Hamiltonian in eq. (10) , which naturally includes the coupling between electronic states with different momentum k and the scattering with phonons at different momentum q.
In the NAMD k method, to get the real-time carrier dynamics, the e-ph coupling matrix element g mnν (k, q) and the normal mode coordinate Q qν (t) (t) are required. g mnν (k, q) can be calculated by the DFPT method using primitive cell. [13] Q qν (t) can be obtained using different methods. For example, it can be expressed in terms of phonon populations as Q qν (t) = l qν [(2n qν + 1)/2] 1/2 (e −iω qν t + e iω -qν t ), (11) where the initial population of phonons at t = 0 (t 0 ) is given by the Bose-Einstein distribution n qν = 1 e ω qν /k B T −1 . It can also be obtained from the molecular dynamics simulation using normal mode decomposition method. [14] The normal mode coordinate of a specified phonon can be achieved by decomposing nuclear displacement from a molecular dynamics trajectory. More details can be found in the Supplementary Materials. Finally, the fewest switches surface hopping (FSSH) is applied to include the stochastic factor of the carrier dynamics. More details can be found in the Supplementary Materials.

Results and Discussion
To verify the validity of the NAMD k method, we choose graphene as a prototypical system and simulate the hot electron relaxation process, which has been investigated intensively. [15][16][17][18][19][20][21][22][23][24][25] The band structure of graphene has six Dirac cones near the E f at K and K points in the 1 st Brillouin zone (BZ), which is also known as six valleys. As discussed in the Supplementary Materials, the NAMD k approach using 9 × 9 × 1 k-grid can successfully reproduce the NAMD r results using 9 × 9 × 1 supercell. Using NAMD k, it is found the hot electron dynamics behavior reaches convergence when the k-grid is as dense as 150×150×1, which is difficult to be simulated by NAMD r approach. The hot electron relaxation may involve intervalley and intravalley e-ph scattering. We first study the single electron dynamics initially excited at 1.0 eV above the E f (E ini = 1.0 eV) in the K valley. To get statistics on the quantum behavior of an excited electron, we sample 30 trajectories by randomly setting the initially excited electron at different k-points in K valley at 1.0 eV above the E f . Figure 1 (a) shows five snapshots of the hot electron population on the band structure plotting within 2 ps and Figure 1 (b) indicates the corresponding hot electron distribution in the 1 st BZ. It can be seen although the hot electron is initially excited in K valley, the K-K intervalley scattering almost immediately starts. Figure 1 (c) presents the time-dependent electron population in K and K valleys. The valley lifetime (τ K ), which is defined as the timescale when the equilibrium between K and K is reached, is around 0.4 ps. The intervalley scattering suggests that the hot electron couples with phonons with large momentum. Figure 1 (d) shows the hot electron energy relaxation, in which the color bar indicates the electron population and the blue line represents the averaged energy. It can be seen that there is an energy bottleneck for hot electron relaxation, which locates at around 0.2 eV above the E f . Above this energy, the relaxation is a relatively fast process, which corresponds to energy relaxation from 1.0 to around 0.2 eV within 2 ps. Using a gaussian function, the lifetime for this fast energy relaxation (τ E ) can be estimated to be 0.56 ps. Furthermore, a quantize character with an energy difference of around 0.2 eV as indicated in Figure 1 (c). Following the fast process, there happens a much slower relaxation process from 0.2 eV to the Dirac point. The timescale of the slow process is difficult to be estimated with a 2 ps simulation. The dependence of τ E on E ini . When E ini > 0.2 eV, there are a fast and a slow relaxation processes, and τ E f ast is plotted. When E ini < 0.2 eV, there is only slow relaxation process and τ E slow is plotted. The green line represents the reciprocal of density of mk ,nk = ψ mk |∆V |ψ nk between state ψ mk and ψ nk plotted with different E mk and E nk scale.
To understand further the slow relaxation process close to the E f , we perform a 100 ps NAMD k simulation for hot electron relaxation with E ini = 0.1 eV. The energy and valley dynamics are shown in Figure 2 (a-b). Both the energy and valley relaxation become much slower and there is no longer a fast relaxation process for the energy relaxation. τ E and τ K are estimated to be 1.2 and 3.0 ns, respectively. Such a relaxation behavior is completely different with E ini = 1.0 eV. We further study the hot electron relaxation dynamics with different E ini from 0.1 to 1.5 eV. As shown in Figure 2 (c), it is found that E ini = 0.2 eV is a critical point for different relaxation behaviors. If E ini > 0.2 eV, there will be both a fast and a slow relaxation process. The τ E for Article the fast process (τ E f ast ) ranges from 0.3 to 3 ps, which follows the inversely proportional trend to the density of states (DOS). When E ini < 0.2 eV, there is only the slow process and τ E slow dramatically increases by 1-3 orders of magnitude. The correlation between τ K and E ini is shown in Figure 2 (d), and it shows a very similar trend with τ E . It suggests that when E ini > 0.2 eV, the equilibrium between K and K valley can be reached at ps timescale due to frequent intervalley e-ph scattering. Whereas, when E ini < 0.2 eV, the intervalley e-ph scattering becomes rare and τ K reaches ns timescale. The slow relaxation process is dominated by the intravalley scattering.
In the NAMD r approach, NAC expressed in Eq. (4) is the crucial term determining the carrier dynamics. Accordingly, in NAMD k approach, e-ph coupling H ep mk ,nk = ψ mk |∆V |ψ nk between state ψ mk and ψ nk plays the key role. In Figure 2 (e-f) we plot the averaged H ep mk ,nk , where the x and y axes represent the energy of ψ mk and ψ nk (labeled as E mk and E nk , respectively). In Figure 2 (e), where E mk and E nk range within [0.0, 1.5] eV, the largest H ep mk ,nk can be roughly fitted by two lines which are expressed as E mk − E nk = 0.2 eV , suggesting the coupling between two electronic states reaches the largest when their energy difference is around 0.2 eV. Thus, when E ini > 0.2 eV, the hot electron prefers to relax to an electronic state 0.2 eV lower in energy. It perfectly explains the quantize character with an energy difference of 0.2 eV observed in the fast relaxation process shown in Figure 1 (d), suggesting that hot electron relaxation is strongly coupled to phonons with energy around 0.2 eV. When E ini < 0.2 eV, as indicated by the square marked with the white dashed lines in Figure 2 (f), the couplings between E mk and E nk are much smaller. And the matrix elements close to the diagonal line, where E mk − E nk is very small, play a crucial role. It implies that in this case, the coupling to the phonons with small energies is essential.
During the hot electron relaxation, the energy of the electrons transfers to phonons through e-ph coupling. And Figure 3 shows the phonon excitation dynamics along with the hot electron relaxation. Figure 3 (a) shows the four snapshots of phonon excitation within 2 ps with E ini = 1.0 eV. It can be seen that within the hot electron relaxation process, only the optical modes A 1 and E 2g , which belong to LO and TO branches, are notably excited. In addition, there is minor excitation for LA and TA modes. Figure 3 (b-c) shows the timedependent phonon occupations and energies of these four different phonon modes. It can be seen that the excitation of the optical A 1 and E 2g modes with ω A 1 = 0.16 eV and ω E 2g = 0.19 eV are dominant. The excitation of A 1 is responsible for the intervalley electron scattering and both A 1 and E 2g contribute to the quantized character in the energy relaxation process and together define the critical energy threshold around 0.2 eV. Figure 3 (d-f) presents the time-dependent phonon excitation with E ini = 0.1 eV. In this case, the LA and TA phonons around Γ point excitation are dominant as shown in Figure 3 (d) and (e). Due to the thermal expansion of the electronic states, the optical phonon modes A 1 also have minor contribution. But since its energies are much higher than the LA and TA phonons, its contribution to excited phonon energy is still the major part. Since the energies and momentum of LA and TA phonons are both very small, the energy and valley dynamics are much slower. The phonon excitation results perfectly explain the distinct hot electron relaxation behavior when E ini is above or below the energy threshold of 0.2 eV. Finally, we study the multi-hot electron relaxation by simulating electron temperature (T e ) decrease. In the TR-ARPES measurements, after photoexcitation, the hot electron will reach equilibrium with a certain temperature through electron-electron (e-e) scattering, and then relax to a lower temperature through e-ph coupling. Figure 4 (a) shows 5 snapshots in the T e relaxation with initial electron temperature T eini = 3193 K. In this case T e decreases to 639 K at 10 ps. Figure 4 (b) shows the time-dependent relaxation dynamics with T eini = 3193, 2200, and 1060 K. For all three cases, T e converges to around 500 K at 10 ps. The relaxation from 500 K to lower temperature is very slow.
The simulation by NAMD k can explain the most important experimental results. Different experimental groups reported that when E ini < 0.2 eV, slow decay via acoustic phonons plays a role. [26][27][28][29][30] With E ini > 0.2 eV, the hot electron can efficiently scatter with the optical phonons. [19,23,26,30] The timescale is in the same magnitude with our results. Especially, the strong coupling with A 1 and E 2g modes and quantized energy-loss is perfectly in agreement with the report by Na et al. [23]. It also agrees with the theoretical work by Bernardi et al. based on quasi-classical time-dependent Boltzamnn transport equation (rt-BTE). [25] It can be noted that hot electron relaxation mainly excites A 1 , E 2g , LA and TA. All of these phonons are in-plane modes. The out-of-plane phonon modes are not directly excited and thus we propose the buckling of graphene will need the phonon-phonon interaction, which can be studied by rt-BTE method. [25] Comparing with the NAMD r, the NAMD k method has the great advantage treating the solid system. First, only the e-ph coupling based on the calculation of primitive cell is required for NAMD k. By contrast, for NAMD r, the AIMD simulation and the electronic structure calculations for each time step along the AIMD trajectory using a supercell are required. The computational cost is significantly reduced by NAMD k. Second, the |ψ nk (r; R 0 ) basis sets used here can be understood as diabatic basis sets. Each electronic state has a clear notion of band index and momentum. We do not need to re-order the electronic states when they cross each other, thus the well-known state crossing problem in surface hopping can be avoided. [31][32][33][34][35] Third, in the NAMD k method, the e-ph coupling can be understood in the time domain. Especially, the time-dependent phonon excitation induced by electron relaxation can be achieved, which is helpful to understand the photo-induced lattice structure distortion, which can not be achieved by the NAMD r method using CPA, where the lattice distortion is pre-determined by the AIMD simulation. Finally, since the e-ph coupling considering the SOC effects, as well as the exciton-phonon coupling elements can be calculated based on previous studies, the NAMD k method can be easily extended for the investigation of SOC promoted spin dynamics and exciton dynamics using GW + real-time Bethe-Salpeter equation (GW+rtBSE). The major approximation used in the NAMD k method is the harmonic approximation. It will be very important and interesting to include the anharmonic effects in this approach in the future.

Summary
The newly developed NAMD k approach breaks the bottleneck of NAMD application on carrier dynamics in momentum space. Based on harmonic approximation, e-ph coupling matrix can be used to calculate the coupling between electronic states with different momentum, and the phonon excitation dynamics can be simulated simultaneously. NAMD k approach is applied to study the hot carrier dynamics in graphene. It is found the optical phonon mode defines an energy threshold which separates the hot electron carrier dynamics into fast and slow regions, which are dominated by intervalley and intravalley e-ph scattering, respectively. We propose the NAMD k approach to be an indispensable tool in exploring excited carrier dynamics in momentum space, and it will pave a new pathway to the theoretical design of devices and materials for optoelectronics, valleytronics and solar energy conversion. i ∂ ∂t |Ψ(r; R(t)) =Ĥ el (r; R(t)) |Ψ(r; R(t)) .
(S1) Article Here, both electronic HamiltonianĤ el (r; R(t)) and electronic state |Ψ(r; R(t)) are dependent on time through nuclear trajectory R(t).
Begin with the equilibrium atomic structure R 0 , we separate the electronic Hamiltonian into two part H el (r; R(t)) =Ĥ 0 (r; R 0 ) + ∆Ĥ(r; R(t)), where the first partĤ 0 (r; R 0 ) is equilibrium Hamiltonian at equilibrium positions R 0 , which is time-independent. The second part ∆Ĥ is the variation of Hamiltonian, which is induced by nuclear displacements ∆R(t) = R(t)−R 0 . As the nuclear coordinates only relate to potential term V (r; R) in electronic HamiltonianĤ el =T +V , the ∆Ĥ can be written as Then, we expand electronic wavefunction by the eigenstates of equilibrium Hamiltonian where k represents the crystal momentum and n is the band number. the basis set |ψ nk (r; R 0 ) are calculated in the framework of Kohn-Sham (KS) density functional theory (DFT), and satisfŷ where nk is eigen energy of KS state |ψ nk . Here, we should notice that the expansion basis wavefunctions |ψ nk are time-independent. By substituting Eqs. (S2-S5) into Eqs. (S1) and multiplying by ψ mk | on both sides of equation, we get the time-depent equation of expansion coefficient where H 0 mk ,nk is diagonal energy matrix element and H ep mk ,nk is e-ph coupling (EPC) matix element Before the following discussions, let me introduce some useful notations. We consider M nuclei or ions in the unit cell. The position vector and Cartesian coordinates of the nucleus κ in the primitive unit cell are denoted by R κ and R κα , respectively. We describe the infinitely extended solid using Bornvon Kámánn (BvK) boundary conditions. In this approach, periodic boundary conditions are applied to a large supercell which contains N p unit cells, identified by the direct lattice vectors L p , with p = 1, 2, ..., N p . The position of the nucleus κ belonging to the unit cell p is indicated by R pκ = L p + R κ . The Bloch wave vectors q are taken to define a uniform grid of N p points in one unit cell of the reciprocal lattice, and the vectors of the reciprocal lattice are indicated by G.
To calculate the EPC term, we expand the potential energy V (r; R(t)) in terms of nuclear displacements ∆R(t). The potential to first order in displacements is We can prove that ∂ qκ v(r; R 0 ) and ∆ qν v(r; R 0 ) are lattice-periodic functions. Then, Eq. (S11) can be written as ∆V (r; R(t)) = N −1/2 p qν e iq·r ∆ qν v(r; R 0 )Q qν (t)/l qν . (S15) For electronic states, wave function of KS state |ψ nk can be expressed in the following Bloch form: with u nk a lattice-periodic function. The wave function ψ nk is taken to be normalized to one in the supercell, while the periodic part u nk (r) is normalized to one in the crystal unit cell. By combining Eq. (S8), (S15) and (S16), we have (S17) where g mnν (k, q) is e-ph matrix element g mnν (k, q) = u mk+q |∆ qν v(r; R 0 )|u nk uc . (S18) Here the subscript "sc" and "uc" indicate that the integral is carried out within an entire supercell and one unit cell of the crystal, respectively.

Surface Hopping Probability
The FSSH method represents time-evolving electron-nuclear system by an ensemble of trajectories, propagating under the influence of deterministic (via TDSE) and stochastic (via SH) factors. Evolution of the system of interest is defined in a joint space combining classical phase-space for nuclei and discrete quantum states for electrons. The paths are constructed such that the FSSH probabilities for all states at all times, averaged over the trajectory ensemble, are equal to the corresponding probabilities obtained from the TDSE. The latter are given by diagonal elements of the density matrix: The off-diagonal elements determine the probabilities of transitions between electronic states. For simplicity, we take the diagonal elements ρ nk,nk as ρ nk Article where ω nk = nk . Then we consider a two states system, where electron occupies state of nk and hopping to state of mk .
Here, EPC is given as: where Q 0 qν is the amplitude of the Q qν (t). For a real system, number of unit cells N p → +∞, so EPC will be a small value, H ep mk ,nk → 0. As a result, occupation of initial state will change very slow, that meansc nk (t) ≈ 1, then The coefficient of final statẽ (S32) For ω qν → ω f i , we ignore the term e i(ω f i +ω qν )t −1 i(ω f i +ω qν ) . Theñ The probability becomes here, we define As we can see in the figure, the function F (t, ω f i − ω qν ) has a sharp peak at ω f i = ω qν . Do to this behavior we say that F, and thus the probability, shows a resonant nature. The half of distance between the first two zeros are defined as the resonance width ∆ω = 2π/t (S36) As t → ∞, the function F (t, ω f i − ω qν ) → 2πδ(ω f i − ω qν )t. That means we need a long time t to guarantee energy conservation. However, in the SH method we calculate the probability in a short time ∆t, and propagate within CPA, as a result, we can not guarantee energy conservation. For a short time step ∆t, the probability P nk→mk (t; t + ∆t) = 2Im ρ * mk ,nk (t)H ep mk ,nk ∆t ρ nk (t) (S37) TheH ep mk ,nk is EPC with energy conservation correction, N p ν |g mnν (k, q)| 2 [n qν δ( ω f i − ω qν ) + (n qν + 1)δ( ω f i + ω qν )] (S38) Here we use gaussian function to fit δ function, where σ represents the spread of energy conservation.
To choose reasonable energy spread, we take tests for different σ to simulate the electronic dynamics with intitial energy of 1.0 eV. As shown in Figure S1, the electron dynamics get convergent as the energy spread increases. For energy spread larger than 10 meV, the dynamics nearly do not change. But we have Article to point it out that the convergent energy spread is related to the interval between electron states. In graphene, the Dirac conical band structure results in a very small density of electronic states near the Fermi surface, so a sufficiently dense k -point grid is required to accurately model the hot electron relaxation process. We tested the calculation of different k -point densities, as shown in Figure S2. It can be seen that with the increase of k -point density, the electron relaxation process gradually tends to converge. When the density of k -points reaches 90, the relaxation process is no longer changed. However, it should be noted that as the electron energy gets closer to the Dirac point, the density of its electronic states decreases, so the denser mesh is required. In order to simulate the finer electronic relaxation process, we use the k -point grid of 150 × 150 × 1 in simulations for E ini = 0.4-1.5 eV, and k -point grid of 450 × 450 × 1 in simulations for E ini < 0.3 eV.

Camparison with original NAMD simulation
In order to verify the rationality of the new method, we compare the calculation results of the NAMD k with the original method NAMD r. Because the original method requires the supercell to simulate the electronic relaxation process of multiple k -points, we can only simulate the supercell with 9 × 9 × 1 unit cells due to the limit of computational amount. Therefore, we compare this result with the result of NAMD textitk using the k -point grid of 9 × 9 × 1. As shown in Figure S3, the calculation of the new method is basically consistent with that of the original method. This also shows the limitations of the Fig. S3 (a) Original NAMD simulation within 9 × 9 × 1 supercell. (b) NAMD with EPC simulation with k-points grid of 9 × 9 × 1.
NAMD r method for strongly dispersive systems, while the NAMD k method can be simulated using a much denser k -point grid.