Microscopic theory of ultrafast dynamics of carriers photoexcited by THz and near-infrared linearly-polarized laser pulses in graphene

We investigate the dynamics of photoexcited carriers and nonequilibrium phonons in graphene by solving the microscopic kinetic Bloch equations. The pump and drift effects from the laser field as well as the relevant scatterings (including the Coulomb scattering with dynamic screening) are explicitly included. When the pump-photon energy is high enough, the influence of the drift term is shown to be negligible and the isotropic hot-electron Fermi distribution is established under the scattering during the linearly polarized laser pulse investigated here. However, in the case with low pump-phonon energy, the drift term is important and leads to a net momentum transfer from the electric field to electrons. Due to this net momentum and the dominant Coulomb scattering, a drifted Fermi distribution different from the one established under static electric field is found to be established in several hundred femtoseconds. We also show that the Auger process investigated in the literature involving only the diagonal terms of density matrices is forbidden by the dynamic screening. However, we propose an Auger process involving the interband coherence and show that it contributes to the dynamics of carriers when the pump-photon energy is low. In addition, the anisotropically momentum-resolved hot-phonon temperatures due to the linearly polarized light are also investigated, with the underlying physics revealed.


I. INTRODUCTION
Graphene is a promising material which has received intensive attention both theoretically and experimentally due to its excellent transport and optical properties. [1][2][3][4][5][6][7][8][9][10][11][12] Especially, its linear band structure and zero band gap lead to unique optical and electrical properties, making it suitable for various optoelectronic applications. 1,8 A comprehensive understanding of the optical and electrical properties is one of the prerequisites for its potential applications.
Experimentally, the time-resolved pump-probe measurement is widely used to investigate the dynamics of carriers in graphene. [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] With this method, both the differential reflection and the differential transmission (DT) after pump are measured. [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] In the experiments with medium photoexcited electron density N p e (around 10 12 cm 2 ) and probe-photon energy much higher than the equilibrium Fermi energy, the positive DT is often observed with its fast relaxation of several hundred femtoseconds followed by a slower picosecond relaxation. 13,20,21,25 The mechanism leading to these two distinct relaxations is attributed to the rapid carrier-phonon thermalization and the slow decay of hot phonons, respectively. 20,22,31 On the other hand, in recent years, experiments with higher photoexcited electron densities (N p e > 10 13 cm −2 ) 23, 26,29 show that a negative DT appears when the probe-photon energy (> 1500 meV) is much higher than the equilibrium Fermi energy. 23,26 Such negative DT was usually claimed to arise from the renormalization of the single particle energy by the electron-electron Coulomb Hartree-Fock contribution. 23,26 It is noted that although the observed DT with high probe-photon energy and low photoexcited electron density is usually positive, negative DT was also reported with N p e = 3 × 10 11 cm −2 by Plochocka et al. where the sample is 70-100 layer epitaxial graphene exhibiting the Dirac-type electron spectrum. 19 In contrast to the negative DT with probe-photon energy much higher than equilibrium Fermi energy, other types of negative DTs with lower probe-photon energy have already been observed in previous works. 14, 15 Sun et al. reported the appearance of the negative DT under weak pump intensity with the probe-photon energy between 528 to 697 meV, 14 which is lower than twice of the Fermi energy (E F ∼ 350 meV). This negative DT was suggested to come from the weakening of the Pauli blocking due to the heating of the electrons by the pump pulse. 14 However, this claim has not yet been verified by the microscopic theoretical approach. Soon after that, the negative DT was also observed by George et al. in the case when the probe-photon energy is as low as a few tens of meV (with the frequency about several terahertz). 15 With such low probe-photon energy, the increase of the intraband absorption is assumed to be responsible for the negative DT. 15 In addition to these phenomena, with the pump-probe method, the optical gain has also been reported in optically pumped graphene very recently. 28 In addition to these experimental works, many theoretical studies have been devoted to get a clear comprehension on the dynamics of photoexcited carriers and phonons. 20,25,29,[31][32][33][34][35][36][37][38][39] Among these works, the simple rate equation approach 20,25,39 assumes the establishment of the Fermi distribution as a starting point and, therefore, fails to incorporate the pump process and the buildup of the Fermi distribution. The widely used time-dependent Boltzmann equation approach 31,34,38 can fulfill the above aims but neglects the interband coherence, 35 which should be very important during the pump process. More elaborative microscopic investigations with the interband coherence included have been carried out based on manyparticle density-matrix approach. 33,[35][36][37] However, these theoretical investigations are still far away from completion, as stated in the following.
Among those papers with the interband coherence included, 33,[35][36][37] the scattering terms under the Markovian approximation with the Coulomb screening in the static limit are applied. Nevertheless, with the static screening, the electron-electron Coulomb scattering diverges when the initial and final states of the scattered electrons satisfy the condition v F |q| = | ω|, 40,41 where q and ω are the differences in momentum and energy between the initial and final states, respectively. To avoid this divergence, in these papers, the delta function originated from the Markovian approximation and reflecting the energy conservation in the Coulomb scattering term, was replaced by a Lorentzian function with the broadening attributed to the non-Markovian effect. 42 Nevertheless, under the more natural dynamic screening, 31,[43][44][45] this broadening is not needed. This is because that the dynamic screening also diverges in a higher order when v F |q| = | ω|, 46,47 and therefore reduces the original divergent scattering amplitudes to zero.
Another important issue frequently investigated is the influence of the Auger process. In previous works, the interband coherence was neglected and the Auger process was defined as the process that one of the two scattered electrons is transferred from one band to the other while the other electron remains in the same band. 33,34,36 In the works by Winzer et al., 33,36 such process was particularly stressed. However, under the dynamic screening, this Auger process is totally forbidden because it also satisfies the strong restriction v F |q| = | ω| due to the energy and momentum conservation. However, in the case with the optical transition, the interband coherence is generated. With the interband coherence, there are Coulomb scattering processes with three out of the four band indices in the Coulomb interaction Hamiltonian being the same and the left one being different. This is also a kind of Auger process. However, the influence of such process has not been addressed in the literature.
Finally, in the microscopic equations obtained based on the many-particle density-matrix framework by Knorr et al., 33,[35][36][37] the pump term was obtained under the vector potential gauge. However, the drift term which describes the acceleration of the electrons by the electric field is not included. The drift term may be unimportant when the pump-photon energy is very high as investigated in these works. Nevertheless, it can be expected to be important when the pump-photon energy is low. On the other hand, in the time-dependent Boltzmann equation approach, the drift term is included but the pump term in the equa-tion has to be introduced based on a Fermi-golden-rulelike approximation and hence is suitable only for a weak pump intensity. 34 Moreover, the Rabi flopping 48 is absent in this approach. Therefore, a general equation based on the microscopic many-particle density-matrix approach that includes both the drift and pump terms naturally is needed.
In this paper, we perform a microscopic investigation on the dynamics of photoexcited carriers and phonons in graphene by setting up and solving the kinetic Bloch equations with the electron-phonon, electron-impurity and electron-electron Coulomb scatterings explicitly included. 48,49 In our study, the dynamic screening 46,50-52 is adopted in the Coulomb scattering. Moreover, with the gauge invariant approach, 48 both the drift and pump terms are obtained naturally. We look into the dynamics of carriers photoexcited by the linear polarized laser pulses with the pump-photon energy in both nearinfrared and THz regimes under medium pump intensity. When the pump-photon energy is high enough, the influence of the drift term is shown to be negligible compared to that from the pump process. Moreover, the anisotropic photoexcited electrons tend to be isotropic under the scattering and an isotropic hot-electron Fermi distribution is established before the end of the pulse investigated here. On the contrary, when the pump-photon energy is low, the drift term causes a net momentum transfer from the electric field to electrons. Together with the dominant Coulomb scattering, a drifted Fermi distribution is established in several hundred femtoseconds. Moreover, the form of the drifted Fermi distribution differs from the one obtained under static electrics field. 45 Besides, we also show that the temporal evolution of the DT measured by Hale et al. can be well fitted with our microscopic calculation. 25 In addition, although the Auger process involving only the diagonal terms of the density matrices is forbidden by the dynamic screening, we show that the Auger process involving the interband coherence still contributes to the dynamics of carriers. However, its contribution is important only when the pump-photon energy is low. Similarly, it is shown that the terms neglected in the rotation-wave approximation, 48 which is widely accepted in semiconductor optics, also become important in the case of low pump energy. We also investigate the dynamics of phonons. Due to the linearly polarized pump pulse, the q-resolved hot-phonon temperatures are anisotropic. Moreover, the anisotropic hotphonon temperatures under high pump energy are very different from those under low pump energy. Finally, we investigate the negative DT by fitting the temporal evolution of DT in the experimental work by Sun et al.. 14 Our results support their suggestion that the negative DT comes from the weakening of the Pauli blocking due to the heating of the electrons by the pump pulse.
This paper is organized as follows. In Sec. II, we set up the model, lay out the kinetic Bloch equations and then give a simple analysis on the kinetic Bloch equations. In Sec. III the results obtained numerically from the kinetic Bloch equations are presented. We summarize in Sec. IV.

II. MODEL AND FORMALISM
We start our investigation from graphene on SiO 2 and SiC substrates. In the absence of external field, the energy dispersion of graphene near the K and K ′ points can be described by the effective Hamiltonian 53 ( ≡ 1) Here, v F represents the Fermi velocity, µ = 1 (−1) for K (K ′ ) valley and τ are the Pauli matrices in the pseudospin space formed by the A and B sublattices of the honeycomb lattice. The eigenvalues of H µ,AB 0 are ε ηk = ηv F k and the corresponding eigenstates are ψ µ kη = (µηe −iµθ k , 1) T / √ 2 with η = 1 (−1) for conduction (valence) band and θ k representing the polar angle of k. The additional electron-optical-filed interaction Hamiltonian under the vector potential gauge is given by 54 Note that this Hamiltonian is consistent with the one obtained from the tight-binding approach. 55,56 Here, we only investigate the linearly polarized light with the electric field along the x direction. For simplicity, the vector potential is taken to be Our investigation is performed in the base set of the eigenstates of the Hamiltonian [Eq. (1)]. Under this base set, the Hamiltonians are given by with σ representing the Pauli matrices in the conduction and valence band space,k indicating the direction of k andẑ standing for the unit vector along the z direction.
Exploiting the nonequilibrium Green's function approach with the gradient expansion as well as the generalized Kadanoff-Baym ansatz, the gauge invariant kinetic equations 48 of the carriers under the Markovian approximation read 45,49,57,58 (see Appendix A for detailed derivation of the coherent and drift terms) (6) in which ρ µk stand for the density matrices with their diagonal terms ρ µk,ηη = f µk,η representing the electron distribution functions and the off-diagonal terms ρ µk,1,−1 = ρ * µk,−1,1 standing for the interband coherence. The drift terms can be written as which describe the acceleration of electrons by the laser field E. The coherent terms are given by Here, [A, B] ≡ AB − BA is the commutator and H µ Pump = −|e|v F µA x sin θ k σ y /c comes from the laser field describing the pump process.
Σ HF µk = − qŜ µ k,k−q ρ µk−qŜ µ k−q,k V r (q, 0) stand for the Hartree-Fock self-energies, 59 in whichŜ µ k,k−q are the form factors with their elements being S µ k,k−q,ηη ′ = ψ µ † kη ψ µ k−qη ′ and V r (q, ω) = V 0 q /ǫ(q, ω) represent the screened twodimensional Coulomb potentials with V 0 q = 2πv F r s /q being the bare Coulomb potential. Here r s stands for the dimensionless Wigner-Seitz radius. 51,60-63 ǫ(q, ω) is the dielectric function under the dynamic screening, 48,51,60 which is given by 46,[50][51][52] (9) It is noted that the laser field E can not only pump electrons but also accelerate them, as revealed in our gauge invariant kinetic Bloch equations [Eq. (6)] 48 by the pump and the drift terms, respectively. However, in previous works by Knorr et al., 33,[35][36][37] the drift term is not included and hence the case of low pump energy, for which the drift term is expected to be important, can not be investigated with their equations. In contrast, our kinetic Bloch equations can deal with both the high and low pump energies completely.
Under the laser field with pump-photon energy ω L , the off-diagonal terms of the density matrices are generated around the resonant absorption point (where 2v F k ≈ ω L ). They oscillate with the frequency around ω L . Following the approach adopted in semiconductor optics, 48 we introduce the slowly varying interband-polarization component P µk = ρ µk,1,−1 e iωLt . Then, the kinetic equations can be rewritten into with the detuning ν µ σ z ] and the Rabi frequency When the pump-photon energy is high enough, these equations can be further simplified by neglecting the high-frequency oscillating terms (rotation-wave approximation) 64 and the Rabi frequency becomes Then, one finds that, except for the drift and scattering terms, Eq. (10) is of the same form as that in semiconductor optics. 48 However, when the pump energy is low, the rotation-wave approximation is no longer valid and the Rabi frequency must be calculated with Eq. (11). We further define a new set of matricesρ µk witĥ ρ µk,η,η = f µk,η andρ µk,1,−1 =ρ * µk,−1,1 = P µk . Under such definition,ρ µk vary slowly with time and the scattering terms in Eq. (10), which include the electronelectron Coulomb, electron-impurity, electron-optical phonon and electron-remote-interfacial (RI) phonon scatterings, are given as in which [>←→<] stands for the same term in the previous [ ] but interchanging >←→<, and {η ←→ η ′ } * represents the conjugate of the same term in the previous { } except the exchange of η and η ′ . In Eq. (15), N ±λ q = n λ ±q +1/2±1/2 with n λ q standing for the distribution function of the phonon in branch λ; V a (q, ε η ′ k+q − ε ηk ) = [V r (q, ε η ′ k+q −ε ηk )] * ;ρ > µk = 1−ρ µk andρ < µk =ρ µk . The optical phonons include the transverse optical phonons (KO) near the K (K ′ ) point and the longitudinal (LO) as well as transverse optical (TO) phonons near the Γ point. The detailed forms of the scattering matrix elements U (ε η2k − ε η2k−q ) and M µµ ′ λ k,k−q,ηη1 are given in Appendix B.
In this work, the q-resolved dynamics of the RI and optical phonons are also studied, while the acoustic phonons are set to be in equilibrium with the environment at temperature T 0 . The kinetic equations of the hot phonons are in which ∂ t n λ q | ep come from the carrier-phonon scattering and ∂ t n λ q | pp describe the anharmonic decay of hot phonons in relaxation time approximation. They are given by where τ pp is the phenomenological relaxation time from the phonon-phonon scattering and n 0 λ,q is the number of the λ branch phonons at environmental temperature T 0 .
Before we show our numerical results, we first give a simple analysis on the Auger process. If only the diagonal terms of the density matrices (i.e., η = η ′ , From this equation, one finds that it describes the scattering process with one electron being scattered between bands η and η 1 while the other electron between bands η ′ 2 and η ′ 1 under the energy conservation ). Moreover, the vertices of the involved two scattering pro- , are conjugated with each other. For the Auger-type scattering process here, it means that three out of the four involved band indices (η, η 1 , η ′ 1 and η ′ 2 ) are the same and the left one is different, corresponding to the case that one of the two scattered electrons is transferred from one band to the other while the other electron remains in the same band. Then, taking η 1 = η ′ 1 = η ′ 2 = η for example and substituting the energy dispersion, the requirement of energy conservation becomes must be satisfied. Although this condition is obtained in the special case η 1 = η ′ 1 = η ′ 2 = η, one can prove that it is valid as long as the band indices in the delta function is Auger type. On the other hand, it has been shown that, under this condition, the dielectric function under the dynamic screening ǫ(q, ) are reduced to zero and the corresponding processes are forbidden. Nevertheless, if the scatterings involving the off-diagonal terms are taken into consideration, terms like indicates an Auger type scattering process for the scattering with guarantees that the scattering process restricted by the energy conversion [whose vertex ] is non-Auger type, which means that the requirement Eq. (20) is no longer necessary. Therefore, these scattering terms are permitted and are expected to contribute to the dynamic of carriers. Moreover, for the corresponding scattering term [Eq. (21)], there is a prefactor Since the other factors in the scattering term all vary slowly with time, the contribution of such scattering term oscillates very fast when ω L is large and hence it has little effect on the dynamics of carriers. However, when ω L is small enough, this kind of Auger process is expected to have marked influence on the evolution of electron distribution function.

III. RESULTS
The kinetic Bloch equations, Eqs. (6) and (16) with all terms included are numerically solved following the method detailed in Ref. 65, except the drift term. In this paper, the drift term is treated up to the third order to suppress the noise (see Appendix C). Then, the temporal evolutions of the carrier and phonon distribution functions are obtained. From Eq. (10), one finds that f µk,η = f −µk,η and P µk = −P −µk . Therefore, we  only show the results in K valley (µ = 1) with the valley index µ no longer appearing for simplicity. With the evolution of carrier distribution, the temporal evolution of the optical transmission at the probe-photon energy ω pr is calculated from 13,66,67 where n ref is the refractive index of the substrate and N lay is the number of graphene layers. Here, the corresponding interband optical conductivity, which is determined by carrier distributions, is given by 13,54,68,69 ωpr ]/T 0 ωpr , with T 0 ωpr representing the transmission before the pump. We present our numerical results with material parameters given in Table I. The full width at half maximum (FWHM) of the pump pulse is chosen to be 100 fs (corresponding to σ t = 42.5 fs) and τ pp = 2.8 ps. 78 The equilibrium chemical potential is taken to be 18.6 meV (the corresponding electron density is 1.43 × 10 11 cm −2 ), the environmental temperature is set to be 300 K and the impurities are absent, unless otherwise specified.

A. High pump-photon energy
We first investigate the dynamics of carriers and phonons with high pump-photon energy. The substrate is chosen to be SiO 2 . In the calculation, ω L = 1500 meV (corresponding to the near-infrared light), unless otherwise specified.
as function of k with ωL = 150 meV along the x and y directions at different times. Here, the pump process is switched off. The black dotted curve in (c) is (vF|k|−µ1)/(kBT1) with T1 = 333 K and µ1 = 11.4 meV. To make the results clearer, the curve for t = 100 fs plotted in (c) is ln(f −1 k,1 −1)+1. (d) Temporal evolution of the total momentum knet(t) with ωL = 150 meV (solid curve with the scale on the left hand side of the frame). The temporal evolution of the electric field −E(t) is also plotted for reference (dotted curve with the scale on the right hand side of the frame).

Influence of drift term
We first investigate the influence of the drift term on the dynamics of electrons. To make it clearer, the pump process is switched off by excluding H µ Pump in Eq. (8). In the calculation, the maximal electric field in Eq. (3) is taken as E 0 = 150 kV/cm and the pump-photon energy ω L = 150 meV. The k dependences of ln(f −1 k,1 − 1) along the k x and k y directions at different times are plotted in Figs. 1(a) and (b), respectively. The dynamics of holes is similar and is not shown here. From the figures, one finds that before t = −100 fs the influence of the electric field is negligible. Moreover, the electron distribution at that moment is the Fermi distribution, with as shown in the figure. Here, T η and µ η are the electron temperature and chemical potential in corresponding band η, respectively. After t = −100 fs, the electrons are drifted under the influence of the electric field [Eq. (3)] along the x direction. From this oscillating electric field, a net oscillating momentum k net (t) = k kf k,1 (t)/ k f k,1 (t) along the x direction is gained by electrons. As a result, the location of the minimum of ln(f −1 k,1 − 1) along the k x directions oscillates around k = 0, as shown in Fig. 1(a). However, the minimum along the k y directions remains fixed at k = 0 but its magnitude oscillates, as revealed in Fig. 1(b). To further illustrate the influence of the electric field on the net momentum k net (t), we plot the temporal evolution of k net in Fig. 1(d). The electric field −E(t) is also plotted in the same figure for reference. One finds that its phase is always π/2 in advance of k net . Moreover, after the pulse (t > 100 fs, when |E| < 0.1|E 0 |), the net momentum gained from the electric field tends to zero. In addition, it is found that for curves at t = −100, 0 and 100 fs in Fig. 1(a), the magnitudes of slopes of the curves away from the minima decrease with the temporal evolution, indicating that electrons are heated by the drift pulse [see Eq. (24)]. Then, with the negligible net momentum, the isotropic hot-electron Fermi distribution is expected to be established under the Coulomb scattering after the pulse. To show this, we plot the k dependences of ln(f −1 k,1 − 1) along the k x and k y directions at t = 100 and 130 fs in Fig. 1(c). One finds that at t = 100 fs, the curve is almost linear with k for each direction, which indicates the establishment of the Fermi distribution along that direction. To further show the buildup of the isotropic Fermi distribution, one fits the distribution with Eq. (24) along the direction θ k and obtains the corresponding hotelectron temperature T ηθ k and chemical potential µ ηθ k . Then, the establishment of the isotropic Fermi distribution is identified when T ηθ k along different directions are very close. Here, it is found that at t = 130 fs the relative difference of the temperature increase T 1θ k −T 0 along different directions is less than 10 %.
Then, we investigate the influences of the pumpphoton energy and impurity density on the thermalization of electrons by the drift term. We plot the pumpphoton energy and impurity density dependences of the hot-electron temperature T 1 in Fig. 2. With different ω L and N i , the isotropic Fermi distributions are established at different times. Here, T 1 is fitted along the k x direction at t = 230 fs, when the isotropic Fermi distribution has just been established (relative difference between temperature increase T 1θ k − T 0 along the k x and k y directions is less than 10 %) for the slowest case in our calculation. The impurity density dependence of T 1 is obtained with ω L = 300 and 500 meV. It is shown that T 1 decreases monotonically with the increase of the pump-photon energy. Moreover, a peak appears in the impurity density dependence as shown in Fig. 2. These phenomena can be well understood from the Boltzmann equation under the relaxation time approximation: which is the same as the intraband conductivity given in the literature. 39,54,67 From this equation, it is found that the conductivity and hence the energy absorbed from the electric field decreases monotonically with the increase of ω L and a peak exists in its N i (or equivalently 1/τ p ) dependence. These phenomena can also be understood more physically. One finds from Eq. (25) that the conductivity tends to zero in the limit of no scattering or infinitely large scattering. This can be comprehended as the effect of the electric field on the electron is totally cancelled out in one oscillation period in the case without scattering and the electrons can not be drifted by the electric field in the limit of infinitely large scattering strength. Then, in the case with finite scattering strength, a peak must appear as shown in Fig. 2. As for the decrease of T 1 with the increasing ω L , it is under-stood as the electron becomes more and more difficult to follow the laser field with the increase of the oscillation frequency.
To further compare the influence from the drift and the pump terms, we also perform the calculation with the pump process included but the drift term excluded. In the calculation, ω L = 500 meV, N i = 0 and the other parameters remain the same. At the time in Fig. 2, i.e., t = 230 fs, the isotropic Fermi distribution is checked to have already been established and the obtained hotelectron temperature is T 1 = 854 K, much larger than 7 K established by the drift term only. In the sample usually reported in the experiments, 79,80 the mobility is larger than 1000 cm 2 /Vs (N i < 2.5 × 10 12 cm −2 ). Within this impurity density range, the temperature increase due to the drift term is less than 37 K. This is much smaller than that from the pump process (the influence of impurity on the pump process is small with such pump intensity and impurity density range). Moreover, since the thermalization due to the drift effect decreases with increasing ω L , the drift term is negligible in the case with pump-photon energy higher than 500 meV. This is consistent with the investigation in semiconductors where the energy of the pump photon is very high due to the large band gap. 48 However, in the case with very low pump-photon energy as we will investigate in the next subsection, the drift term is expected to be very important, especially for the case where the electric field does not change direction during the pulse and a large net momentum is transferred to the electrons via the drift term.

Dynamics of electrons
Before the investigation on electron dynamics, we first calculate the DT and compare it with the experimental data in single-layer graphene by Hale et al. 25 [ Fig. 3(a)]. It is noted that this experiment has been fitted in our previous work. 31 However, in that paper, the Fermi distribution was assumed to be established at all the time and the hot phonons are averaged with a fitting parameter E max , which describes the upper energy limit of the hot carriers being able to emit phonons. In this paper, with the inclusion of the pump term and the q-resolved hot phonons, the buildup of the Fermi distribution can be studied elaborately and the fitting parameter E max is no longer needed. Then, the influence of the approximations in our previous paper can be checked. In the calculation here, the probe-photon energy ω pr = 1100 meV, the FWHM is taken to be 180 fs and the pump-photon energy ω L = 1500 meV, as indicated in the experiment. It is noted that with N i = 0 and such high pump-photon energy, the drift effect of the laser field is negligible as shown in the previous subsection. In addition to these parameters, the fitting parameter τ pp = 2.8 ps and the maximal electric field E 0 = 300 kV/cm which leads to the maximal electron density being 2.7×10 12 cm −2 . This electron density is comparable with the experimental estimation and that in our previous paper (4.6 × 10 12 cm −2 ). For the relaxation of DT which begins from t = 133 fs, as shown in Fig. 3(a), one finds that our numerical results agree very well with the experimental data with a fast relaxation followed by a slow one with the relaxation times being 0.28 and 1.33 ps, respectively. The discrepancy between our numerical results with the experimental data at the initial increase of the DT is understood to come from the approximation in our model where the interference between the pump and probe pulses is neglected and the Markovian approximation 48 is applied.
Then, we investigate the dynamics of electrons, especially the buildup of the Fermi distribution under the linearly polarized light. We plot ln(f −1 k,1 − 1) as func-tion of k along different directions at different times in Figs. 3(b)-(d). With the strong scattering and the negligible drift effect, the distribution is center symmetric and hence the distribution functions along the −x and −y directions are not shown. Besides, the dynamics of holes is almost the same 31 and is not plotted here. From Fig. 3(b), one finds that the distribution of the pumped electrons are anisotropic and mainly around k = 1.14 nm −1 ≈ ω L /(2 v F ) along the y-axis as shown by the valley there. This comes from the linearly polarized pump light as shown by sin θ k in Eq. (11). During the pump process, the scattering also plays an important role by smearing out the valleys and making the distribution tend to be isotropic as shown in Figs. 3(b) and (c). Moreover, it is found that the slopes of the curves around the Dirac point decrease with the temporal evolution, indicating their thermalizations under the scattering. At t = 131 fs shown in Fig. 3(d), one finds that the curves become almost linear with k for each direction. By fitting the curves with Eq. (24), the obtained hot-electron temperature T 1θ k along the x and y directions are 2160 and 2221 K, respectively with the relative difference for these two directions being less than 3 %. By noticing that the pulse [Eq. (3)] ends after t = 165 fs (when |E| < 0.1|E 0 |), one concludes that in the case with such long pump pulse width, the hot-electron Fermi distribution is established during the pulse and once it is established, it is an isotropic one. Considering that the relaxation of the DT begins at 133 fs which is after the establishment of the Fermi distribution, the approximation of the buildup of the Fermi distribution in our previous work is acceptable. However, this approximation leads to the inaccuracy of the hot-electron temperature (3200 K at the beginning of the relaxation in our previous work 31 ).
To further show the influence of the Auger process involving the interband coherence, we also plot the results without the Auger process in Figs. 3(b)-(d). One finds that the exclusion of the Auger process has little influence on the evolution of the distribution function. This is in agreement with our argument in Sec. II that the Auger process which survives the dynamic screening in the presence of the optical polarization does not contribute due to the high pump-photon energy. Similarly, also due to the high pump energy, the rotation-wave approximation is checked to be valid (not shown in the figure).

Dynamics of phonons
With the inclusion of the q-resolved hot phonons, we are able to give a more detailed investigation on the dynamics of phonons than our previous work. 31 We calculate the temperature T ph (q) for phonons with momentum q at branch λ by T ph (q) = ω λ /[k B ln(1 + 1/n λ q )] and then, plot them in momentum space at t = 131 fs for the KO and RI 1 phonons in Figs. 4(a) and (b). The qresolved temperatures of the other phonons are similar to T ph (q) (K) these two phonon branches. It is shown that due to the anisotropic distribution of electrons, the thermalization of phonons under the electron-phonon scattering is also anisotropic. Moreover, it is found that the thermalization of phonons with small q (q < ω qλ /v F , which are the phonons within the black dotted circle) is negligible.
This can be understood from the energy conservation during the electron-phonon scattering process [Eq. (15) one finds that the phonons with q < ω qλ /v F and q > ω qλ /v F are involved in the interband (η 2 = η 3 ) and intraband (η 2 = η 3 ) electronphonon scattering, separately. Therefore, the negligible thermalization of the phonons with q < ω qλ /v F means that before the establishment of the hot-electron Fermi distribution (t ≤ 131 fs), the recombination of electrons and holes due to the interband electron-phonon scattering is marginal. In order to show this clearer, we plot the q dependence of the angle averaged phonon temperature T ph q = 2π 0 dθ q T ph (q)/(2π) at different times in Figs. 4(c)-(f), as done by Knorr et al.. 32,35 The locations v F q = ω qλ for each branch of phonons are indicated by arrows in Figs. 4(c)-(f). It is found that in the first several hundred femtoseconds, only the thermalization of the phonons with v F q > ω qλ is efficient. After that, the temperatures of the phonons with v F q < ω qλ increase markedly. Especially, for the KO and RI 2 branches, the hottest phonons appear among the phonons with v F q < ω qλ as shown in the inset of Fig. 4(f). This is consistent with the results by Knorr et al. 32,35 and can be understood as the interband electron-phonon scattering happens only among the low energy electrons with its energy |ε ηk | < ω qλ . However, the energy of the photoexcited electrons are very high (ω L /2 > 3.9ω qλ ) and before the low-energy electrons are heated under the scattering, the energy gained by the electrons first relaxes through the intraband electron-phonon scattering process. Then, after the thermalization of the low energy electrons, the heating of the phonons with v F q < ω qλ starts to be important.
One also observes from Figs. 4(a) and (b) that the hottest phonons are along the q y direction. This comes from the angular dependence of the electron-phonon scattering matrices. Since the hottest phonon appears with q > ω qλ /v F , we only consider the intraband scattering process. Then, with only the diagonal terms of the density matrices taken into consideration, the scattering matrices in Eq. (15) |M µµ ′ λ k,k−q,ηη | 2 are propor-tional to 1 − cos(θ k − θ k−q ) for the KO phonons and 1+cos(θ k −θ k−q ) for the RI 1 phonons. 63,70,71 One immediately understands that the parallel or antiparallel scattering are the strongest. Then, since the photoexcited electrons are mainly along the k y direction, the hottest phonons, which come from the electron-phonon scattering, also have to be along the y direction as shown in the figures. In addition, Figs

B. Low pump-photon energy
We then investigate the dynamics of electrons and phonons with the pump-photon energy ω L down to the THz regime. The maximal pump field in Eq. (3) is taken as E 0 = 5 kV/cm. The substrate is again chosen to be SiO 2 .

Dynamics of electrons
In this subsection, we investigate the dynamics of electrons with ω L = 15 meV (3.6 THz). With such low pump-photon energy, the corresponding oscillation period of the electric field is 276 fs, which is even larger than twice of the FWHM. Therefore, during the pump process, the laser field E almost keeps along one direction and a large net momentum is transferred to electrons. In this case, the drift term is very important since it describes the momentum transfer. Furthermore, due to the transferred net momentum, a drifted Fermi distribution after the pulse is expected. In previous works with static electric field, 45,81 it has been shown that in the case with strong Coulomb scattering, the drifted Fermi distribution is expressed as with u η being the drift velocity. However, in graphene, it has been shown by Zhou and Wu that the Coulomb scattering is not strong enough compared with the drift effect of the static electric field and the drifted Fermi distribution is instead described by 45 Although both distributions are the same in semiconductors with parabolic spectrum, they have different properties in graphene due to the linear dispersion. For Eq. (27), one finds that u η = k kf k,η (t)/ k f k,η (t) is nothing but the center-of-mass drift velocity. However, in Eq. (26), u η is only an effective parameter. Moreover, in the k dependence of ln(f −1 k,1 − 1) along the direction of u η , Eq. (26) indicates a minimum at k = 0 and the magnitudes of the slopes for θ k and θ k + π directions are η ln(f −1 k,η − 1) along θ k = 0 and π directions for conduction band electrons (η = 1) at t = −50, 25, 100, 300 and 500 fs, as well as that for valence band electrons (η = −1) at t = 300 fs. Here, ωL = 15 meV and the impurity density Ni = 0. The black dotted curve is calculated from Eq. (26) with T1 = 323.6 K, u1 = 0.12vF and µ1 = 16.3 meV. (b) Temporal evolution of T ηθ k and u ηθ k (with the scale on the right-hand side of the frame) for electrons in conduction (solid curves) and valence (dashed curves) bands obtained by fitting Eq. (26) along directions θ k = 0 and π (red squares), as well as θ k = π/4 and 5π/4 (blue dots). different. However, Eq. (27) shows that the minimum deviates from k = 0 and the magnitudes of the slopes are identical. To reveal the distribution established here and the dynamics of electrons, we plot ln(f −1 k,1 − 1) as functions of k along the direction of the electric field E (i.e., θ k = 0 and π) at different times in Fig. 5(a). One finds that during the pump process, the minimum of the curve is first drifted away from the Dirac point as shown by the curves at t = −50 and 25 fs. However, before the pulse ends (t > 100 fs, when |E| < 0.1|E 0 |), the minimum has already started to tend to the Dirac point and the magni-tudes of the slopes for the curves away from Dirac point become different. Then, at t = 300 fs, the minimum is around the Dirac point and the curves along θ k = 0 and π directions are almost linear with k, but with different magnitudes of slopes. With such different magnitudes of slopes, the corresponding drift of the center-of-mass k net is 0.02 nm −1 . This indicates the establishment of the distribution Eq. (26), as shown by the fitting plotted as the black dotted lines in the figure. This is understood that the drifted Fermi distribution is established after the pulse in our present investigation. Then, with the Coulomb scattering being the dominant scattering, the distribution as Eq. (26) is established.
It is also shown in Fig. 5(a) that at t = 300 fs, the similar drifted hot-electron Fermi distribution in the valence band (η = −1) has also been established since the minimum of the corresponding curve is around k = 0 and the curves along θ k = 0 and π directions are almost linear with k with different magnitudes of slopes. Moreover, the distributions along the other directions θ k are also similar. To further investigate the buildup of the unitary drifted Fermi distribution along all directions, we fit the distributions of electrons along different directions in both the conduction and the valence bands with Eq. (26). The obtained T ηθ k and u ηθ k are plotted in Fig. 5(b). One finds that after t = 300 (380) fs, the relative difference between the drift velocities along different directions for electrons in conduction (valence) band is less than 5 % and the relative difference of the temperature increases T 1θ k − T 0 (T −1θ k − T 0 ) is less than 10 %. This indicates the buildup of the unitary drifted Fermi distribution as Eq. (26) for all directions in conduction (valence) band. It is noted that the case studied here is n-doped, which leads to different buildup times, hot-electron temperatures and drift velocities of the unitary drifted Fermi distributions for electrons in conduction and valence bands. However, under the interband Coulomb scattering, the temperatures and drift velocities in conduction and valence bands tend to be identical. From our calculation, the relative difference of the drift velocities is less than 10 % after t = 650 fs and that of the temperature increases is less than 10 % after t = 350 fs, as shown in the figure.
We also investigate the influence of impurities and pump-photon energy on the drift velocity of electrons. u 1θ k obtained from the distributions along θ k = 0 and π directions under different pump-photon energies and impurity densities are plotted in Fig. 6(a). It is noted that after t = 300 fs, the relative differences of the drift velocities along different directions are less than 5 % for all cases investigated here. One finds that in the case without impurity, the drift velocity u 1θ k and hence the total momentum relax slowly (with the relaxation time being about 1.8 ps). However, when the impurities with N i = 5 × 10 11 cm −2 (the corresponding mobility is about 5000 cm 2 /Vs) are introduced into the system, the relaxation of the total momentum becomes very fast (with the relaxation time being 46 fs) and at t = 300 fs, the corresponding u 1θ k = 0.003v F is negligible. This is because that the relaxation of the drift velocity comes from the electron-impurity and electron-phonon scatterings. However, the electron-phonon scattering strength is weak and only when the impurity density is large (e.g., 5 × 10 11 cm −2 ), can the relaxation of the net momentum become fast. Besides, it is also shown in the figure that when ω L increases from 15 to 30 meV, the value of u 1θ k at t = 300 fs decreases markedly. This is because that with ω L = 30 meV, the laser field E changes direction during the pulse and the net momentum gained by electrons is partially cancelled out.
We further investigate the influence of the pump and Auger processes on the evolution of the photoexcited electron density by switching them off separately. The results are plotted in Fig. 6(b). In the calculation, ω L = 15 meV and N i =0. It is noted that without the pump process, electrons can still be excited from the valence band to conduction band. This is because that the electrons are heated by the drift effect of the laser field, which leads to the interband electron-phonon scattering. However, one finds from the figure that the dominant mechanism for the excitation of the electrons is the pump process in the case investigated here. Moreover, it is also found that the Auger process indeed affects the dynamics of electrons and it leads to the change of the excited electron density at the end of the pulse (t = 100 fs) as large as 22 %, which is consistent with our analysis in Sec. II. Furthermore, with such low pump-photon energy, the widely accepted rotation-wave approximation in semiconductors optics is no longer valid. The temporal evolution of the photoexcited electron density under the rotation wave approximation is also plotted in Fig. 6(b). One finds that the fast oscillation of the electron density during the pump process disappears and the excited electron density at the end of the pulse changes about 17 %.

Dynamics of phonons
We further investigate the dynamics of phonons. We plot the q-resolved temperatures of the hot RI 1 and KO phonons with ω L = 15 meV under the impurity densities N i = 0 and 5 × 10 11 cm −2 in Fig. 7. The other phonons are similar and not shown. The temperatures of these two branches of phonons are obtained at t = 300 fs when the unitary distribution of electrons for all directions has been established. From the figures, one finds that in the case without impurity, the q-resolved temperatures of the hot phonons are highly anisotropic. Moreover, by considering that u η in both conduction and valence bands are along the θ k = 0 direction, it is shown that the hottest phonons appear in the direction parallel to u η and the coldest phonons antiparallel to u η . This phenomenon is also shown in semiconductors under static electric field. 82 It can be understood from the generation rate of phonons due to the electron-phonon scattering [Eq. (17)] after the establishment of the drifted Fermi distribution of electrons. With only the diagonal terms of the density matrices taken into consideration, Eq. (17) is given by (29) Here, to simplify the analysis, the temperatures of phonons and electrons are taken to be the same T and the chemical potentials and drift velocities for electrons in different bands are also set to be unified µ and u, respectively. Then, due to the energy conservation δ(ε η2k+q − ε η1k − ω qλ ), Eq. (29) is simplified to e (η1vFk−u·k−µ+ω qλ )/(kB T ) [1 − e −u·q/(kBT ) ]. One immediately understands that the hottest phonons appear among phonons with u parallel to q and the coldest phonons appear among the phonons with u antiparallel to q as shown in Figs. 7(a) and (b). This phenomenon can also be understood physically from the relaxation of the drift velocity of electrons. During the relaxation, the net momentum gained from the laser field is transferred from the electrons to the phonons, which implies the tendency of the absorption of the phonons with q · u η < 0 and the emission of the phonons with q · u η > 0. On the other hand, in the system where the net momentum relaxes very fast through the electron-impurity scattering, the distributions of phonons become more isotropic as shown in Figs. 7(c) and (d).
C. Negative DT in system with high Fermi energy In this subsection, we further investigate the negative DT due to the weakening of the Pauli blocking. Sun et al. measured the DT with high probe-photon energy (> 500 meV) and showed that a peak exists in the temporal evolution of DT. 14 Moreover, it is also found that there is a range of the probe-photon energy where the negative DT appears during the evolution. To understand this, they suggested that their epitaxial sample consists of a stack of graphene layers 13,83,84 with some layers heavily doped and some layers undoped. Moreover, the Fermi energy of the heavily doped layers is assumed to be around half of the upper boundary of the probe energy range for the appearance of the negative DT. Then, with the energies of the resonant absorption states of the probe pulse smaller than the Fermi energy in the heavily doped layers, the corresponding distribution difference f kω,1 − f kω,−1 in Eq. (23) increases in the undoped layers but decreases in the heavily doped layers after the pulse. This leads to the decrease of the conductivity in undoped layers but increase in heavily doped ones. These two opposite tendencies compete in the evolution and are responsible for the peak and the negative DT. Moreover, with the decrease of the probephoton energy, the variation of the conductivity from the undoped layers is strengthened whereas that from the heavily doped layers is weakened. Furthermore, if the probe-photon energy is higher than twice of the Fermi energy in the heavily doped layers, the conductivities in both the heavily and undoped layers decrease. These properties lead to the fact that the negative DT appears only in a probe-energy range smaller than twice of the Fermi energy.
To check above conjecture, we fit their experimental data with our microscopic approach. The substrate is chosen to be SiC, the environmental temperature is 50 K, the FWHM of the pulse is 100 fs and the pumpphoton energy ω L = 1550 meV as indicated by Sun et al.. 14 The maximal electric field in Eq. (3) is taken to be E 0 = 175 kV/cm. With such pump intensity, the photoexcited electron density is about 4 × 10 11 cm −2 per layer. The probe-photon energy ω pr = 550 meV is the same as that in the experiment. It is noted that with such high probe energy, the intraband absorption is suppressed. 68 To model the multilayer structure in their sample, N lay σ ωpr (t) in Eq. (22) is replaced by N h σ h (t) + N u σ u (t) with σ h (t) and σ u (t) being the conductivity in heavily and undoped graphenes, respectively and N h and N u being the corresponding layer numbers. In our computation, σ h (t) and σ u (t) are calculated separately with the corresponding equilibrium Fermi energy being µ 0 = 350 and 0 meV, respectively, as estimated in the experiment. N h = 2 and N u = 20 are taken as fitting parameters which are very close to those estimated in the experiment. The results are plotted in Fig. 8(a). One finds that our numerical results repeat the peak and the negative DT shown by the experimental data. Especially, when t > 500 fs, the fitting is quite well. The discrepancy between our numerical results with the experimental data for the magnitude of the peak is again from the approximation mentioned in Sec. IIIA2 as well as the simplification of the complex layer configuration to two kinds of doped layers. The electron distributions before and after the pulse are plotted in the inset of Fig. 8(a), which is consistent with the arguments by Sun et al.. 14 Moreover, when the probe-photon energy is small enough or larger than twice of the Fermi energy, the DT keeps positive as shown by the curves with ω pr = 520 and 710 meV in Fig. 8(b). Based on these results, our microscopic investigation supports their argument that the negative DT comes from the weakening of the Pauli blocking in the heavily doped layers.

IV. SUMMARY
In summary, we have investigated the nonequilibrium dynamics of carriers and phonons in graphene under the linearly polarized near-infrared and THz laser pulses via the microscopic kinetic Bloch equation approach. 48,49 The carrier-impurity, carrier-phonon and carrier-carrier Coulomb scatterings are explicitly included and the dynamic screening for the Coulomb scattering is utilized. Moreover, based on the vector potential gauge and the gauge invariant approach, 48 both the drift and pump terms are naturally included in our approach.
We first investigate the thermalization of electrons due to the drift term. It is shown that the thermalization is weakened with the increase of the pump-photon energy and it is even negligible when the pump-photon energy is higher than 500 meV. Moreover, a peak appears in the impurity density dependence of the thermalization. Then, the dynamics of carriers and phonons with the pump-photon energy up to the near-infrared regime is studied. We show that the temporal evolution of the DT reported by Hale et al. 25 can be well fitted. Moreover, under the linearly polarized laser pulse, the electrons are photoexcited anisotropically. However, their distribution tends to be isotropic under the scattering. As a result, the isotropic hot-electron Fermi distribution is found to be established within t = 131 fs, which is before the end of the pulse (the FWHM of the pulse is 180 fs). This is consistent with the previous works. 26,31,35 In contrast to the high pump-photon energy case, in the case with the pump-photon energy down to the terahertz regime, we show that a large net momentum is transferred to electrons through the drift term and a drifted Fermi distribution different from the one under static field is established in several hundred femtoseconds. 45 We further show that the Auger process investigated in the previous works 33,34,36 which only involves the diagonal terms of density matrices is forbidden by the dynamic screening. However, the Auger process involving the interband coherence contributes to the dynamics of carriers. Nevertheless, its influence is important only when the pumpphoton energy is low. In contrast, for the rotation-wave approximation which is widely accepted in semiconductor optics, we show that it fails when the pump energy is low. We also apply our microscopic theory to investigate the negative DT in the case where the equilibrium Fermi energy is high and the probe-photon energy is around twice of the Fermi energy. We show that the negative DT appears due to the weakening of the Pauli blocking under the pump pulse, which supports the suggestion by Sun et al. in their experimental work. 14 In addition, the momentum-resolved hot-phonon temperatures are also studied. When the pump-photon energy is high, the hot phonon temperatures are anisotropic due to the linearly polarized laser field and the hottest phonons appear with their momentum in the direction perpendicular to the laser field E. Moreover, we also show that the electron-hole recombination due to the interband electron-phonon scattering is marginal in the first several hundred femtoseconds. On the other hand, when the pump-photon energy is low and the relaxation of the net momentum is slow, the hot-phonon temperatures are also highly anisotropic with the hottest and the coldest phonons appearing in the directions parallel and antiparallel to the drift velocity, respectively, due to the drifted Fermi distribution of electrons. However, if the relaxation of the net momentum is fast as in the case with high impurity density, the hot-phonon temperatures become more isotropic.
In order to satisfy the requirement of the numerical stability, m ′ in F θ,m ′ µ,nm,ηη ′ is taken to be m if −|e|E(t) · x sin θ m+0.5 > 0 and m + 1 otherwise. Moreover, n ′ in F r,n ′ m ′ µ,nm,ηη ′ is chosen to be n if −|e|E(t) ·x cos θ m > 0 and n + 1 otherwise. Besides, if n ′ − 1 ≥ 0, m ′ is taken to be m. On the other hand, if n ′ − 1 < 0, the term n ′ − 1 in the equation is replaced by 0 and m ′ = m + M/2. It is noted that with this numerical scheme, the temperature variation due to the noise is within 3 K for the pulse investigated in Fig. 2.