Intrinsic and extrinsic effects on intraband optical conductivity of hot carriers in photoexcited graphene

We present a numerical study on the intraband optical conductivity of hot carriers at quasiequilibria in photoexcited graphene based on the semiclassical Boltzmann transport equations (BTE) with the aim of understanding the effects of intrinsic optical phonon and extrinsic coulomb scattering caused by charged impurities at the graphene{substrate interface. We employ iterative solutions of the BTE and the comprehensive model for the temporal evolutions of hot-carrier temperature and hot-optical-phonon occupations to reduce computational costs, instead of using full-BTE solutions. Undoped graphene exhibited large positive photoconductivity owing to the increase in thermally excited carriers and the reduction in charged impurity scattering. The frequency dependencies of the photoconductivity in undoped graphene with high concentrations of charged impurities significantly deviated from that observed in the simple Drude model, which is attributed to temporally varying charged impurity scattering during terahertz (THz) probing in the hot-carrier cooling process. Heavily doped graphene exhibited small negative photoconductivity similar to that of the Drude model. In this case, charged impurity scattering is substantially suppressed by the carrier-screening effect, and the temperature dependencies of the Drude weight and optical phonon scattering govern negative photoconductivity. In lightly doped graphene, the photoconductivity changes its sign temporally after the photoexcitation, depending on the carrier and optical phonon temperatures and the pump uence. Moreover, the photoconductivity spectra depend not only on the material property of graphene sample but also on the waveform of the THz-probe pulse. Our approach provides a quantitative understanding of non-Drude behaviors and the temporal evolution of photoconductivity in graphene.

Microscopic models based on semiconductor Bloch equations for graphene have been employed for the quantitative and qualitative understanding of hot-carrier dynamics considering intrinsic effects such as carrier-carrier and carrier-phonon scattering 16,38,39 . However, in addition to the intrinsic effects, extrinsic effects such as charged impurities, surface optical phonons, and dielectric properties on the substrate play essential roles in the electrical and optical conductivities of hot carriers 17,26,28 . Furthermore, most previous studies discussed hot-carrier dynamics based on transient transmission change that determines the response around the center frequencies of the THz-probe pulses. To obtain more quantitative insights into hot-carrier dynamics governed by intrinsic and extrinsic effects, analyses based on the frequency dependence of THz conductivity are necessary. Such effects are essential for understanding the physics underlying various graphene-based device applications. However, it remains a practical numerical challenge to obtain the frequency-dependent optical conductivity for intraband and interband transitions with the full solution of a carrier distribution function by solving the semiconductor Bloch equation in the 2D momentum space of graphene, even if only intrinsic carrier-carrier and carrier-opticalphonon scatterings are considered 16 . Therefore, a suitably approximate method is required for understanding the hot-carrier dynamics affected by extrinsic effects.
In this study, we calculate the frequency-dependent intraband optical conductivity of hot carriers in photoexcited graphene based on the Boltzmann transport equation (BTE), including the intrinsic and extrinsic interactions in the collision term. Because the intraband transition is dominant in the THz-frequency region of 0.1-10 THz, the microscopic polarization for interband transition in the semiconductor Bloch equation may be negligible in the calculation of optical conductivity at the quasi-equilibrium carrier distribution in graphene after the photoexcitation 16,18,22 . Under such conditions, arXiv:2009.10960v1 [cond-mat.mes-hall] 23 Sep 2020 the semiclassical BTE may offer an alternative 40,41 . Under linearly polarized photoexcitations, highly nonequilibrium anisotropic electron and hole distributions are created and rapidly relaxed to the uniform hot Fermi-Dirac (FD) distribution with two different chemical potentials in the conduction and valence bands via carriercarrier and carrier-optical phonon scattering within several tens of femtoseconds 23,38,39 . Following the recombination of the photoexcited electron and hole pairs via Auger recombination and the interband optical phononemission process, the carrier distribution is relaxed to the hot-FD distribution with a single chemical potential 24 . Thereafter, the thermalized carriers and optical phonons at the quasi-equilibrium cool to the equilibrium via the energy transfer from the hot carriers and optical phonons to other types of phonons through phonon-phonon interactions caused by lattice anharmonicity 42 and supercollision (SC) carrier cooling process which are disordermediated electronacoustic phonon scatterings [43][44][45] . Because the energy relaxation caused by the optical and acoustic phonon emission is inefficient for low-energy carriers near the Dirac point, the SC carrier cooling process becomes important. To consider the cooling processes of carriers and optical phonon modes, at least five coupled BTEs expressed as nonlinear integrodifferential equations must be solved for carriers in conduction and valence bands and in three dominant optical phonon modes 46 of graphene. Further reductions in computational costs are still required.
The BTE solution using relaxation-time approximation (RTA) has been used extensively because it makes it easy to solve BTE 47 (see Appendix B). The RTA is valid in the spherical energy band for elastic scattering under low-field conditions. For the inelastic scattering process, RTA is only valid for isotropic scattering in nondegenerate semiconductors, where the FD distribution can be approximated with the Boltzmann distribution. However, in the case of graphene, the RTA solution underestimates or overestimates the relaxation rate for the inelastic scattering because the Boltzmann distribution is not valid. Furthermore, the sub-picosecond temporal resolution of the OPTP experiments cannot instantaneously capture the carrier dynamics. The carrier distribution and momentum relaxation rate in the cooling process change significantly within a probe time that is approximately equal to the THz-pulse duration, and therefore, the THz conductivity cannot be adequately analyzed by the calculation based on the RTA. To ensure calculation accuracy, we employ an iterative method [48][49][50] . This provides the BTE solution for the intraband complex conductivity in graphene near the (quasi-)equilibrium under a weak electric field with appropriate accuracy. Moreover, to calculate the hot-carrier intraband conductivity of graphene, we perform analysis using the BTE combined with a comprehensive temperature model based on rate equations to describe the temporal evolution of hot-carrier temperature and hot-optical-phonon occupations 15,37,42,45,51 . We present numerical simulations of the intraband optical conductivity of hot carriers in photoexcited graphene with different Fermi energies, considering the intrinsic and extrinsic carrier scatterings and the temporal variations of the carrier distribution and momentum relaxation rate during THz probing.
The remainder of this paper is organized as follows. In Section II, we present the iterative solutions of the BTEs for graphene at near-equilibrium under a weak electric field to calculate the direct current (DC) and intraband optical conductivity. Furthermore, the temperature model based on coupled rate equations is introduced to apply the iterative solutions to the quasi-equilibrium hot-carrier state in the cooling process following photoexcitation. Section III presents our main numerical results for the THz conductivity of graphene using different carrier-doping concentrations. Finally, the major conclusions are summarized in Section IV.

A. Iterative solutions of BTE in steady state
The iterative solution of the BTE for obtaining the steady-state conductivity of graphene was introduced in Ref. 50 . The BTE in a homogenous system under a timedependent electric field describes the temporal evolution of the carrier distribution 47,52 , which is given as (1) Here, f λ (k, t) is the electron distribution function for the conduction band (λ = 1) and valence band (λ = −1), k is the wave vector of the carriers, e is the elementary charge, E(t) is the electric field, and ∂f λ (k, t)/ ∂t| collision is the collision term describing the change in the distribution function via carrier scattering. The BTE in a steady state under a constant electric field E , (df λ (k, t)/dt = 0), is given by For spherical bands under a low field, E, the general solution of Eq. (2) is provided approximately by the first two terms of the zone spherical expansion.
Here, f 0 (ε λk ) = 1/ [exp {(ε λk − µ (T e )) /k B T e } + 1] is the FD distribution for the corresponding equilibrium electron distribution at the electron temperature, T e , ε λk = ± v F |k| (ε 1k ≥ 0, for the conduction band. ε −1k ≤ 0 for the valence band) is the electron energy within the Dirac approximation of the graphene energy-band structure 53 , and µ (T e ) is the temperaturedependent chemical potential of the 2D-MDF (see Appendix A and Ref. 26,54,55 ). Moreover, g (ε λk ) is the perturbation part of the distribution, and α is the angle between E and k.
We consider the collision term while accounting for the scattering of electrons with different optical phonon modes η in C η λλ , including both intraband (λ = λ ) and interband (λ = λ ) processes with the elastic scattering processes in C el λ (k). The carrier collision term C η λλ (k) for the interaction of the electron and optical phonons is expressed as where P η k λ kλ and P η kλk λ are the carrier scattering rates obtained by the optical phonon modes η between states (k , λ ) → (k, λ) and (k, λ) → (k , λ ), respectively. This is expressed by accounting for the phonon emission and absorption given by where D η kk is the electron-phonon coupling (EPC)matrix element defined by Ref. 46 (see APPENDIX B). ρ = 7.6 × 10 −7 kgm −2 is the area density of graphene, and ω η and n η are the angular frequency and occupation of the optical phonons, respectively. Density functional theory (DFT) calculations have demonstrated that only three optical phonon modes contribute significantly to the inelastic scattering of electrons in graphene 56 . The first two relevant modes are longitudinal optical (LO) and transversal optical (TO) phonons near the Γ point with energies of 196 meV 57 , which contribute to intravalley scattering. Moreover, zone-boundary phonons have ω K = 161 meV 57 close to the K point and are responsible for the intervalley scattering process. The carrier scattering rates obtained by the optical phonons in Eq. (7) account for phonon emission and absorption. The delta functions δ (ε λk − ε λ k ∓ ω η ) and δ k − k ∓ q in Eq. (7) arise from Fermis golden rule, thereby ensuring the conservation of energy and momentum, respectively. The EPC-matrix elements |D η kk | 2 for Γ-LO, Γ-TO, and K phonons are expressed by 46,56 Here, θ denotes the angle between k and k − k, θ denotes the angle between k and k − k, and θ denotes the angle between k and k . In the case of Γ-LO and K phonons, the plus sign refers to interband processes, and for Γ-TO phonons, it refers to intraband processes. The EPC coefficients were obtained via DFT calculations 58 . Their values are D 2 Γ F = 45.6 (eV/Å) 2 and D 2 K F = 92.1 (eV/Å) 2 . Note that the value of the EPC coefficient D 2 Γ F has been debated owing to the renormalization effect resulting from electron-electron interaction [59][60][61][62][63] . In the numerical calculations, we consider the DFT value for simplicity. The elastic term is given by where P s k k and P s k k are the scattering rates for the elastic scatterings. The index s refers to the different elastic or quasi-elastic scattering modes such as charged impurities and acoustic phonons. Substituting Eqs. (5) and (9) into Eq. (2), we get where E = |E | and k = |k | are the magnitude of the electric field and wavevector, respectively.
are the net in-and out-scattering rates for inelastic scattering. Furthermore, is the total relaxation rate, which is the summation of the relaxation rate, ν s caused by the different elasticscattering processes indicated by s using RTA. We consider the elastic carrier scattering via charged impurities as weak scatterers and acoustic phonons (for the formula of ν el , see Appendix B). Using the contraction mapping principle 64 , Eq. (10) is numerically solved using an iterative procedure for the given inelastic and elastic scattering rates. The (j + 1)th iteration of g j+1 (ε λk ) is considered to satisfy where we arbitrarily select g 0 (ε lk ) = g 0 (ε −lk ) = 0. In this case, S in λ = 0, and Eq. (14) provides the first solution as Equation (15) can be regarded as the solution that has a relaxation time of τ f = 1/ S out i + v el , which demonstrates that the first step of this iterative process can be regarded as the RTA momentum relaxation rate. However, the iterative process must continue until it converges to an appropriate accuracy. In addition to intrinsic optical phonon modes, remote scattering via surface polar optical phonons (SPOP) modes is known to be a limiting factor of electron mobility in graphene on polar substrates and 2D artificial structures such as silicon metal-oxide-semiconductor field-effect transistors. Although we consider only the intrinsic optical phonons of graphene for the inelastic term, a similar procedure can be applied for the SPOP modes using the quasiparticle scattering rate [65][66][67][68] . The effect of carrier-carrier scattering on the intraband conductivity appears from the deviation of the carrier distribution from the FD distribution and the asymmetry of the conduction and valence bands with a prominent contribution when the carrier distribution is far from the equilibrium 40 . In this study, we ignore carrier-carrier scattering because a low electric field causes a small disturbance for the carrier distribution in the 2D-MDF.

B. Iterative solutions of BTE under time-dependent electric field
Here, we extend Eq. (14) regarding the perturbed distribution g (ε λk ) in the steady state under a constant electric field to include arbitrarily time-dependent driving forces. In this case, ∂f λ (k , t)/∂t = 0. We explain the derivation of the iterative solution of the BTE by following the procedure in Ref. 50 . To include time dependence, a constant Ω s ≥ 0 is added to the scattering-out, and the term Ω s g j (ε λk ) is added to the numerator of Eq. (14), which does not affect the solution of g (ε λk ). If there exists a unique g ∞ (ε λk ) of then g ∞ (ε λk ) is independent of Ω s . Furthermore, it is equal to g ∞ (ε λk ) for Ω s = 0, which is the solution to Eq. (14). While the condition Ω s ≥ 0 lowers the convergence rate of the sequence, g j (ε λk ) because g j+1 (ε λk ) approaches g j (ε λk ) as Ω s approaches infinity. Further, we can relate lim Ωs→∞ to ∂g λ (ε λk )/∂t as follows: From Eq. (16), where the final equation follows from the fact that g j+1 (ε λk ) is indistinguishable from g j (ε λk ) when Ω s approaches infinity. Recalling the definitions of S in λ , S out λ , and ν el in Eqs. (11)-(13), we have the Boltzmann equation where f j λ (k ) = f 0 (ε λk ) + g j (ε λk ) cos α. The left-hand side of Eq. (18) is simply identified by ∂f j λ (k )/∂t at time t j = j/Ω s , where 1/Ω s is the time increment between successive iterations. Therefore, the sequence {g j (ε λk )} yields f j λ (k ) versus time when Ω s is sufficiently large compared with S out λ + ν e . Further, Eqs. (17) and (18) may be adopted with a slight modification to include the arbitrarily time-dependent electric field E (t). In this instance, E (t) becomes a function of time through the iteration index j such that The sequence {g j (ε λk )} is used to calculate the fieldinduced current J(t) in graphene as where v F is the Fermi velocity and N (ε λk ) = 2 |ε λk | / π 2 v 2 F is the density of states for 2D-MDF. Subsequently, the intraband conductivity can be obtained by where J(ω) and E(ω) are the Fourier transformations of J(t) and E(t), respectively. Because g j (ε λk ) is a function of electron energy ε λk , the computation cost in the proposed method is decreased significantly compared with the full calculation of f λ (k , t) when solving Eq. (1) in a 2D momentum space. The convergence of Eq. (14) and (19) is demonstrated by calculating the DC and AC conductivities in Appendix E. Note that Eqs. (14) and (19) are valid under a sufficiently weak electric field, which causes the small distribution changes g (ε λk ) from the corresponding equilibrium state.
C. Calculation of carrier temperature and optical phonon occupations in hot carrier cooling process Next, we explain the THz-conductivity calculation procedure for hot carriers in photoexcited graphene at quasiequilibrium. Because the iterative solution of Eq. (19) for the BTE is valid under near equilibrium, it cannot be used directly for calculating the hot-carrier distribution at quasi-equilibrium, which requires the thermalized electron and optical phonon distributions in the cooling process following photoexcitation. During the cooling process, hot carriers lose their energy by emitting strongly coupled optical phonons (Γ-LO, Γ-TO, and K) resulting in a change in the optical-phonon occupation. To consider the cooling process for the hot carriers, we employ a comprehensive model based on the rate equations that describe the temporal evolution of the electron temperature T e and the optical phonon occupations n η 42,43,45,51 (for details, see Appendix C).
In this case, I p represents the energy injected into the graphene sample during laser irradiation, which is assumed to be of the hyperbolic secant form I p (t) = (F ab /2τ exc ) sech 2 (t/τ exc ), where F ab is the absorbed pump fluence, and 2τ exc is the pump-pulse duration. Furthermore, C is the sum of the specific heat of the electrons in the conduction and valence bands. R Net η = R η −G η denotes the total balance between the optical phonon emission and absorption rate, and J sc denotes the energy loss rate for the SC carrier process 43 . R Net M,η = R M,η − G M,η denotes the total balance between the optical phonon emission and absorption rate per number of phonon modes. Moreover, n Γ0 and n K0 represent the phonon occupation near Γ and K points, respectively, in equilibrium at room temperature, and τ ph is the phenomenological phonon decay time via the phonon-phonon interaction caused by lattice anharmonicity. The acoustic phonon occupation is assumed to remain unchanged from the equilibrium state for the picosecond time range following photoexcitation 69 . By substituting the solution of T e and n η obtained by numerically solving the coupled Eqs. (22) to (25) into Eq. (19) during the iteration process, the sequence g j (ε λk ) for the disturbed distribution caused by the applied THz electric field E(t) that includes the temporal evolution of T e and n η in the cooling process can be obtained. The solution of Eq. (19) using T e and n η is not valid for the THz-conductivity calculation of hot carriers with a highly nonequilibrium distribution (non-FD distribution) or hot FD distribution with separate quasi-Fermi levels because Eq. (3) assumes f 0 (ε λk ) as the equilibrium or quasi-equilibrium FD distribution with a single chemical potential.

III. NUMERICAL RESULTS
In this section, we present the numerical simulations performed to investigate the intrinsic and extrinsic carrier scatterings on the intraband optical conductivity of hot carriers in monolayer graphene that has different carrier concentrations. We considered intrinsic carrier scattering mechanisms caused by intrinsic optical and acoustic phonons and the extrinsic scattering mechanisms caused by the charged impurities on the substrate and weak scatterers such as defects and neutral impurities 68 . Figure 1(a) presents a schematic of the OPTP experimental setup used in the simulation. An optical pump pulse with a temporal duration of 2τ exc = 35 fs irradiates a graphene sample at room temperature (T 0 = 295 K) on a substrate at a normal incidence angle. The dielectric constants of the substrate are listed in Table I. Owing to the presence of the substrate, 1.37% of the incident pump pulse is absorbed in the graphene layer if we neglect the saturable absorption of pump pulses in graphene 70,71 . The created nonequilibrium electron and hole pairs are rapidly recombined by the Auger recombination process and interband optical phonon scattering. The nonequilibrium distribution is assumed to change into the quasi-equilibrium hot-carrier state with a single chemical potential, at which point the present procedure can be applied. In doped graphene, the rapid recombination of photoexcited carriers within 150 fs is caused by the limited phase space of the impact ionization process, as reported in Ref 24 . Figure 1(b) presents the temporal waveforms of THz-probe pulse calculated by the second derivative of the Gaussian function exp(−t 2 /τ 2 p ) with pulse durations of 2τ p = 300 and 500 fs. Figure. 1(c) presents the corresponding normalized FFT power spectra that have center frequencies ω/2π of 2.2 and 1.2 THz. The THz-probe pulses are transmitted to the graphene sample at a normal incidence angle at time τ 1 after the photoexcitation, and the waveforms of the transmitted Optical phonon decay time τ ph (ps) 1.0 Pulse duration of pump pulse 2τexc (fs) 35 Pulse duration of THz probe 2τp (fs) 300, 500 THz pulses are assumed to be detected using a timeresolved detection scheme such as via electro-optic sampling by varying the delay time τ 2 between the THz probe and the trigger pulses. The parameters of the graphene sample and experimental setups used in the simulation are summarized in Table I and Figs. 12-16 in Appendices A-C.

A. Undoped graphene
First, we present the numerical results of undoped graphene with different charged impurity concentrations of n i = 0.1 and 1.0 × 10 12 cm −2 , as illustrated in Figs. 2 and 3, respectively. The charge inhomogeneity and disorder in graphene smear out the intrinsic behavior near the Dirac point, thereby resulting in unintended carrier doping and a finite ε F 72 . Thus, we set the finite p-type carrier concentration as ε F = −0.01 eV for the undoped graphene 30 . The corresponding hole concentration at T = 0 K is n c = 7.0 × 10 9 cm −2 . The calculated DC conductivity of the undoped graphene with n i = 0.1 and 1.0 × 10 12 cm −2 at equilibrium without pump fluence are approximately σ DC = 6G 0 and 0.8G 0 , respectively, where G 0 = 2e 2 /h is the quantum conductance. Figure 2(a) presents the temporal evolutions of the carrier temperature (T e ) and the optical phonon temperatures (T Γ , T K ) in the undoped graphene with the low-charged impurity concentration of n iL = 0.1 × 10 12 cm −2 under absorbed pump fluences of F ab = 0.04 and 0.13 µJcm −2 . The optical phonon temperatures were calculated by inverting the Bose-Einstein distribution function n η = 1/(e ωη/kBTη −1). T e increased to almost 1,000 and 1, 500 K. Then, it relaxed to the equilibrium via double exponential decay with fast decay times of τ T1 = 0.25 and 0.37 ps and slow decay times of τ T2 = 2.6 and 2.8 ps. The fast decay with τ T1 corresponds to hot-carrier relaxation for the increased optical phonon emission process. It is governed by the total balance of the energy exchange rate for optical phonon emission and absorption R Net Fig. 15). The slow decay with τ T2 is related to R Net η , the optical phonon decay time τ ph , and the energy-loss rate J sc via the SC carrier cooling process (see Appendix C, Fig. 14(b)). The complex conductivity of graphene σ (ω, τ 1 ) = σ 1 (ω, τ 1 ) + iσ 2 (ω, τ 1 ) in Fig. 2(c) measured at a pump-probe delay time τ 1 after photoexcitation is calculated using the temporal evolutions of the carrier and optical phonon temperatures illustrated in Fig. 2(a). Thereafter, the electric field of the THzprobe pulse E t (τ 2 , τ 1 ) transmitted through the photoexcited graphene at τ 1 is calculated as a function of the difference in time of arrival τ 2 between the THz-probe pulse and an trigger pulse using σ (ω, τ 1 ) (see Appendix D for details).
To understand the temporal evolutions of σ (ω, τ 1 ) in Fig. 2(c) and Fig. 3(c), we exploited the extended Drude model to determine the Drude weight D (ω, τ 1 ) and the momentum relaxation rate Γ (ω, τ 1 ), expressed as D (ω, τ 1 ) and Γ (ω, τ 1 ) are derived from σ (ω, τ 1 ), as fol- lows: Note that while a simple Drude model expressed by has constant values of D and Γ in the frequency space, both D (ω, τ 1 ) and Γ (ω, τ 1 ) in the extended Drude model is frequency dependent indicating the deviation from the simple Drude model. Figures 2(e) and (f) present the temporal evolutions of the frequency dependent D (ω, τ 1 ) and Γ (ω, τ 1 ) for undoped graphene with n iL , respectively. The D (ω, τ 1 ) for n iL at τ = 0.5 ps reaches four times that at equilibrium indicating substantial increases in carrier concentration caused by the high T e . The Γ (ω, τ 1 ) for n iL at τ = 0.5 ps increases to less than two times that at equilibrium indicating the enhancement of the optical phonon scattering (see Fig. 13 in Appendix B). These changes in D (ω, τ 1 ) and Γ (ω, τ 1 ) result in the approximately three times enhancement of σ (ω, τ 1 ) for n iL . Figures 3(e) and (f) present the temporal evolutions of the frequency dependent D (ω, τ 1 ) and Γ (ω, τ 1 ) for n iL , respectively. While the D (ω, τ 1 ) increases similar to that for n iL , Γ (ω, τ 1 ) for n iH decreases significantly as the temperatures increases, thereby resulting in the larger enhancement of σ (ω, τ 1 ) compared with that for n iL . The change of Γ (ω, τ 1 ) for n iH is attributed to the reduction of charged impurity scattering because it is strongly suppressed by the carrier-screening effect that increases with the Drude weight 55 (see also Fig. 13 in Appendix B).
In terms of frequency dependence, the deviation of σ (ω, τ 1 ) from the simple Drude model depends on the charged impurity concentration n i . σ (ω, τ 1 ) for n iL slightly deviates from the simple Drude model, which can be clearly seen in the frequency dependence of the D (ω, τ 1 ) and Γ (ω, τ 1 ) in Fig. 2(e) and (f), respectively. Although D (ω, τ 1 ) and Γ (ω, τ 1 ) for n iL at equilibrium have nearly flat frequency dependencies, both of them show a clear frequency dependence after photoexcitation, which increases with ω; further, their slope becomes larger after photoexcitation. The temporal evolution of D (ω, τ 1 ) can be roughly traced by the temperature dependent Drude weight of 2D-MDF given by D (T e ); however, the difference between D (T e ) and D (ω, τ 1 ) at ω/2π = 2.2 THz, which is the center frequency of the THz probe, becomes larger at a higher temperature as indicated by the colored star symbols in Fig. 2(e). Here, D (T e ) was calculated using T e at the corresponding τ 1 based on BTE considering the temperature-dependent chemical potential µ (T e ) of the 2D-MDF and the constant momentum relaxation rate in graphene (See, Appendix A), which gives the simple Drude-type optical conductivity. On the other hand, σ (ω, τ 1 ) for n iH exhibits a different frequency dependence from that for n iH , i.e., σ 1 (ω, τ 1 ) / σ 2 (ω, τ 1 ) decreases / increases with frequency for ω/2π < 1 THz, which is a clear signature indicating a non-Drude type conductivity. As a result, the similar non-Drude type frequency dependence appears in the photoconductivity ∆σ (ω, τ 1 ) in Fig. 3(d). Such non-Drude type photoconductivity was reported in a previous study 30 where an undoped graphene with ε F = −0.01 eV on a α-quartz substrate obtained by nitrogen doping compensation was used.
The origin of non-Drude type photoconductivity can be understood by considering the temporal evolutions of the hot-carrier distribution and momentum relaxation rate, as well as the broader carrier distribution and energy dependence of the momentum relaxation rate 29,35,36 . To demonstrate this, we calculated the temporal evolutions on the current densities J (τ 2 ) and J F (τ 2 ) induced by the THz-probe pulse with 2τ p = 300 fs as shown in Fig. 4(a) and (c) for n iL and n iH , respectively. Here, J (τ 2 ) and J F (τ 2 ) are calculated by varying the temperatures (T e , T Γ , andT K ) with τ 2 and by fixing the temperatures at τ 2 = 0 ps. For both n iL and n iH , J (τ 2 ) slightly deviate from J F (τ 2 ). However, the normalized difference ∆J/J max for n iH is clearly larger than that for n iL as seen in Figs. 4(b) and (d), where ∆J = J − J F originates from the temporally varying carrier distribution and momentum relaxation rate during the THz probing, and J max is the maximum amplitude of J. The corresponding optical conductivities σ (ω, τ 1 ) and σ F (ω, τ 1 ) for n iL and n iH are illustrated in Figs. 4(e) and (g), respectively. The normalized difference ∆σ F (ω, τ 1 ) /σ max where ∆σ F (ω, τ 1 ) = σ (ω, τ 1 ) − σ F (ω, τ 1 ) and σ max is σ (ω, τ 1 ) at ω/2π = 0.4 THz are plotted in Figs 4(f) and (h).
The frequency dependence of σ (ω, τ 1 ) for n iL is similar to σ F (ω, τ 1 ), and the small deviation arises from the temporally varying hot-carrier distribution and carrier scattering during THz probing. In addition, both σ (ω, τ 1 ) and σ F (ω, τ 1 ) clearly deviate from σ D (ω) obtained by the Drude model fitting for σ (ω, τ 1 ) with n iL . Since σ F (ω, τ 1 ) is calculated by fixing the temperatures, the deviation of σ F (ω, τ 1 ) from σ D (ω) is attributed to the broader carrier distribution and energy-dependent momentum relaxation rate. On the contrary, σ (ω, τ 1 ) for n iH exhibits a larger ∆σ F (ω, τ 1 ) /σ max than that for n iF as seen in Fig. 4(f) and (h), respectively. In this case, the origin of non-Drude-type σ (ω) for n iH is attributed to the charged impurity scattering. Because the carrier screening effect depends on the Drude weight, the charged impurity scattering changes greatly following the temporal evolution of T e . Moreover, we demonstrate that the non-Drude-type σ (ω) originated from the temporally varying hot-carrier distribution and carrier scattering during THz probing, which shows the dependence on the THz probe waveform. Figure 5 presents the results of using the THz probe with the broader pulse duration of 2τ p = 500 fs, which exhibits more pronounced non-Drude-type conductivity because a larger change occurs in the charged impurity scattering during the THz probing time. Compared with ∆J F /J max for 2τ p = 300 fs, that for 2τ p = 500 fs has an additional peak around τ 2 = −0.5 ps when T e increases immediately after the photoexcitation. This peak results in σ (ω, τ 1 ) with the oscillation and the considerably larger deviation ∆σ F (ω, τ 1 ) /σ max .
Unlike the undoped graphene, the frequency dependence of σ (ω, τ 1 ) and ∆σ (ω, τ 1 ) in the heavily doped graphene for both n iL and (h) n iH exhibit slight deviations from the simple Drude model, as indicated in Figs. 6(c) and (d), and 7(c) and (d), respectively. Such frequency dependence was experimentally observed in a large area gated graphene device structure, where ε F could be controlled by the gate voltage as reported in Ref. 26 . Figures 8(a) and (c) present the temporal waveforms of J (τ 2 ) and J F (τ 2 ) in the heavily doped graphene with n iL and n iH , respectively. Compared with undoped graphene, the ∆J F /J max values are significantly suppressed for both n iL and n iH as seen in Fig. 8(b) and (d) because of fewer temporal variations of D (ω, τ 1 ) and Γ (ω, τ 1 ). As a result, the frequency dependence of σ F (ω, τ 1 ) for both n iL and (h) n iH are similar to those for σ F (ω, τ 1 ); however, they slightly deviate from σ D (ω, τ 1 ), as can be seen in Fig. 8(e)-(h). As shown in Fig. 18 in Appendix E, the perturbed carrier distributions g (ε λk ) in the highly doped graphene at T e = 295 and 500 K are confined around Fermi energy ε F = −0.43 eV in the energy space because of the larger Fermi energy ε F compared with the thermal energy. This supports the frequency dependence of σ F (ω, τ 1 ), which resembles that of σ D (ω, τ 1 ) because a simple Drude model assumes a constant momentum relaxation rate in the energy space. However, the g (ε λk ) at higher T e becomes broader resulting in the deviation from σ D (ω, τ 1 ) observed below ω/2π = 1.5 THz, as seen in Fig. 8(e) and (g) because of the energy dependence of the momentum relaxation rate.

C. Lightly doped graphene
In this section, we present the numerical results of lightly p-type doped graphene with ε F = −0.15 eV (the corresponding hole concentration at T = 0 K is n c = 1.6 × 10 12 cm −2 ) with charged impurity concentrations set to n iL = 0.1 × 10 12 cm −2 and n iH = 1.0 × 10 12 cm −2 , as illustrated in Figs. 9 and 10. The corresponding DC conductivities at T 0 = 295 K are σ DC = 21.5G 0 and 7.8G 0 , respectively. Compared with the undoped and heavily doped graphene, the lightly doped graphene exhibits unique temporal evolutions of −∆E t (τ 1 ) /E 0 and ∆σ (ω, τ 1 ), and it strongly depends on n i and F ab . The similar temporal evolutions of −∆E t (τ 1 ) /E 0 presented in this study were observed experimentally 25,26 ; however, their mechanisms have not been fully understood yet. Figure 9(a) presents the temporal evolutions of T e and T η in the lightly doped graphene with n iL . The T e curves with F ab = 0.04 and 0.13 µJ/cm −2 exhibit double exponential decays with τ T1 = 0.55 and 0.21 ps, τ T2 = 2.1, and 2.1 ps, respectively. Moreover, T e increases to 1,000 and 1, 500 K, respectively, which is close to the values of the undoped graphene owing to the similar specific heat capacity C and R Net η (for details, see Appendix C). The −∆E t (τ 1 ) /E 0 curves for both F ab = 0.04 and 0.13 µJ/cm −2 in Fig. 9(b) exhibit negative changes indicating negative photoconductivity as observed in the highly doped graphene. However, the temporal evolution of −∆E t (τ 1 ) /E 0 depends on F ab , which is different from the heavily doped graphene. For a weak pump fluence of F ab = 0.04 µJ/cm −2 , −∆E t (τ 1 ) /E 0 exhibits a sharp negative peak immediately following photoexcitation at τ 1 = 0 ps, which subsequently decreases gradually, taking the minimum value at around τ 1 = 0.8 ps, and then recovering to that of the equilibrium. In con-trast, for F ab = 0.13 µJcm −2 , −∆E t (τ 1 ) /E 0 decreases in two steps. At the first step at τ 1 = 0 ps, it decreases slightly to −∆E t (τ 1 ) /E 0 = −0.0025 and maintains its value until τ 1 = 0.2 ps. At the second step, it starts to decrease gradually again and exhibits its minimum at around τ 1 = 1.8 ps. Finally, it recovers to that of the equilibrium. Figures 9(e) and (i) present σ (ω, τ 1 ) with F ab = 0.04 and 0.13 µJ/cm −2 , respectively, with the corresponding photoconductivity, ∆σ (ω, τ 1 ), as indicated in Figs. 9(f) and (j). It can be observed that the temporal evolution of ∆σ 1 (ω, τ 1 ) at ω/2π = 2.2 THz reproduces that of −∆E t (τ 1 ) /E 0 .
To understand the behaviors of −∆E t (τ 1 ) /E 0 , we calculated the temporal evolutions of ∆D (T e ) = D (T e ) − D (T e = 295 K) and the optical phonon occupation n η following photoexcitation, as illustrated in Figs. 9(c) and (d), respectively. ∆D (T e ) at F ab = 0.04 µJ/cm −2 exhibits two negative peaks. The first sharp negative peak of ∆D (T e ) around τ 1 = 0 ps corresponding to the negative peak of −∆E t (τ 1 ) /E 0 is caused by the sharp increase in T e immediately following photoexcitation because D (T e ) of the lightly doped graphene decreases with increasing T e and takes the minimum around T e = 750 K (see Fig. 12(b) of Appendix A). The second broad negative peak of −∆E t (τ 1 ) /E 0 appears around τ 1 = 0.75 ps. At that time, D (T e ) approaches the minimum around τ 1 = 1.0 ps and the carrier scattering caused by the optical phonons following n η is close to the maximum around τ 1 = 0.5 ps as seen in Fig. 10(c) and (d), respectively. These behaviors can be roughly observed in D (ω, τ 1 ) and Γ (ω, τ 1 ) at around ω/2π = 2.2 THz, as indicated in Figs. 9(g) and (h), respectively. For F ab = 0.13 µJcm −2 , ∆D (T e ) exhibit a sharp negative peak around τ 1 = 0 ps, followed by a large positive peak and a broad negative peak around τ 1 = 2.2 ps. The sharp negative and large positive peaks are caused by the sharp increase in T e up to 1, 500 K at τ 1 = 10 ps because ∆D (T e ) becomes positive for T e > 1, 000 K. While the positive ∆D (T e ) contributes to the positive photoconductivity, the optical phonon scattering following n η approaches its maximum around τ 1 = 0.73 ps and contributes to the negative photoconductivity. Thus, the total balance of the increased D (T e ) and optical phonon scattering lead to the observed two-step drops in −∆E t (τ 1 ) /E 0 . In addition, the broad negative peak of −∆E t (τ 1 ) /E 0 appears around τ 1 = 1.8 ps when the optical phonon scattering is close to the maximum and D (T e ) approaches the minimum.
In the case of the lightly doped graphene with n iH , −∆E t (τ 1 ) /E 0 also exhibits different behaviors depending on F ab because of the increased contribution of the charged impurity scattering. Figures 10(a) and (b) present the temporal evolutions of the temperatures (T e and T η ) and −∆E t (τ 1 ) /E 0 . The temporal evolutions of ∆D (T e ) and n η with n iH are the same with those with n iL in Fig. 9(c) and (d), respectively. After the photoexcitation, the −∆E t (τ 1 ) /E 0 value at F ab = 0.04 µJcm −2 decreases gradually, exhibiting its minimum around τ 1 = 0.8 ps. It then recovers to that of the equilibrium. This behavior is different from that with n iL , which has the sharp negative peak around τ 1 = 0 ps. Because both ∆D (T e ) and the charged impurity scattering contribute to −∆E t (τ 1 ) /E 0 in the same direction, the sharp negative peak of −∆E t (τ 1 ) /E 0 around τ 1 = 0 ps should be more distinguishable than that for n iL . The disappearance of the negative peak around τ 1 = 0 ps for n iH can be understood by considering the temporal evolution of the frequency-dependent Γ (ω, τ 1 ). Whereas Γ (ω, τ 1 ) at ω/2π = 2.2 THz around τ 1 = 0 ps for n iL increases significantly at higher temperatures, and that for n iH decreases slightly (not shown in Fig. 9(h) and 10(f)), thereby resulting in the disappearance of the negative peak. However, σ (ω, τ 1 ) around τ 1 = 0 ps cannot be calculated accurately by the proposed method, and a more detailed discussion on this point requires a full solution of BTE or semiconductor Bloch equation, because the highly nonequilibrium carrier distribution around τ 1 = 0 ps for n iL cannot be described by the perturbed solution of Eq. (19). In contrast, −∆E t (τ 1 ) /E 0 at F ab = 0.13 µJcm −2 in Fig. 10(b) exhibits a sharp negative peak around τ 1 = 0 ps, followed by a large positive peak, which decreases and reaches its minimum at τ 1 = 2.0 ps. It then recovers to that of the equilibrium. This behavior is similar to that of ∆D (T e ) at F ab = 0.13 µJcm −2 in Fig. 10 (c). In this case, the temporal evolution of the frequency-dependent Γ (ω, τ 1 ) also plays an important role. Because Γ (ω, τ 1 ) at ω/2π = 2.2 THz are nearly constant after the photoexcitation as seen in Fig. 10(j), it has only small contribution to ∆σ (ω, τ 1 ) around ω/2π = 2.2 THz. Thus, −∆E t (τ 1 ) /E 0 follows the temporal evolution of ∆D (T e ).
In terms of the frequency dependence, σ (ω, τ 1 ) with n iL for both F ab = 0.04 and 0.13 µJ/cm −2 in Fig. 9(e) and (i) exhibit Drude-like behavior as observed in those of undoped and heavily doped graphene with n iL with the slight deviation caused by a broader carrier distribution and the energy dependence of the momentum relaxation rate. σ (ω, τ 1 ) with n iH in Fig. 10(c) and (g) exhibit the non-Drude-type behaviors depending on F ab and T e . Further, this can be attributed to the temporal variation of τ −1 i of the charged impurity scattering. The non-Drude behaviors in σ (ω, τ 1 ) with n iH appear above T e = 800 K where τ −1 i starts to decrease significantly in lightly doped graphene.

IV. DISCUSSION
In Section III, we presented the numerical results of the intraband optical conductivity at quasi-equilibrium in photoexcited graphene considering the intrinsic and extrinsic carrier scattering mechanisms that exhibited positive and negative photoconductivity depending on the Fermi energy. In undoped graphene, the positive photoconductivity appears and is governed by the positive change in the Drude weight and the negative change of the charged impurity scattering caused by the significant increase in thermally excited carriers. In heavily doped graphene where the charged impurity scattering is greatly suppressed, the Drude weight decreases with T e and the optical phonon scattering becomes dominant as the temperatures (T e and T η ) increase. Both contribute to the negative photoconductivity. In lightly doped graphene, the photoconductivity shows unique temporal evolution depending on the total balance of the contributions from the Drude weight and momentum relaxation rates for the optical phonons and charged impurities.
The deviation of σ (ω, τ 1 ) from the Drude model is remarkable, particularly in the undoped and lightly doped graphene with n iH , which was explained by the temporal variation in the charged impurity scattering during THz probing. In the undoped and lightly doped graphene, the carrier screening effect is weak and changes depending on T e greatly, thereby resulting in the large temporal variation of τ i . It is expected that SPOP modes would also be crucial for non-Drude type photoconductivity of graphene on polar substrates because the temperature dependent carrier screening effect and the dielectric screening govern the coupling of carriers and SPOP modes [65][66][67][68] . In addition, it was found that the measured photoconductivity changed significantly depending on the waveform of the THz-probe pulse. Although the various frequency dependencies of ∆σ (ω, τ 1 ) with non-Drude behaviors resulting from the deviation of the σ (ω, τ 1 ) from the Drude model were reported in previous studies 25,26,29,30,35,36 , the mechanisms were mostly attributed to the broader carrier distribution and energydependent momentum relaxation rate in the hot carrier state, which has not been analyzed quantitatively. Our results provide the answer to the question why graphene exhibits a variety of non-Drude type photoconductivity.
In Ref. 25 , it was reported that the ∆σ (ω, τ 1 ) of doped graphene with ε F = −0.3 eV exhibited a strong dependence on different atmospheric environments with the presence of N 2 , air, or O 2 . The observed ∆σ 1 (ω, τ 1 ) measured at τ 1 = 2 ps after the photoexcitation were negative and exhibited the minimum value around ω/2π = 1.8-2.0 THz, and ∆σ 2 (ω, τ 1 ) changed from positive to negative as ω increased. The observed ∆σ (ω, τ 1 ) could be fitted by the Drude-Lorentz model with the resonant frequency ω 0 /2π; however, it did not depend on pump fluence. Therefore, its origin was attributed to the stimulated emission via population inversion and not plasmon excitation, where the energy gap was assumed to be opened by the adsorption of gas molecules. However, it was reported the population inversion decays within few hundreds of femtoseconds in the monolayer and bilayer graphene 24,73 and the population inversion and stimulated emission is not feasible beyond this time window. A simple Drude model can also reproduce the observed Drude-Lorentz type photoconductivity, ∆σ (ω, τ 1 ), in the frequency range of ω/2π = 0.5 − 3.0 THz if both ∆D and ∆Γ are negative, as demonstrated in Appendix F. However, this situation cannot really occur in photoexcited graphene. Because the negative ∆D in the doped graphene makes the carrier screening weaker, the momentum relaxation rates caused by charged impurities and by optical phonons increase at a high temperature leading to a positive ∆Γ. We performed numerical calculations of ∆σ (ω, τ 1 ) in the lightly doped graphene under various conditions for pump fluence and pump-probe delay time, thereby finding that this Drude-Lorentz type photoconductivity could be reproduced under particular conditions even without considering the stimulated emission as illustrated in Fig. 11.
The temporal evolutions of T e and T η in lightly doped graphene with n iH under strong pump condition at F ab = 0.77 µJcm −2 are presented in Fig. 11(a). The carrier temperature was elevated up to T e = 2600 K by the photoexcitation. It then decreased to 1200 K at τ 1 = 2.0 ps. At this temperature, the positive ∆D (T e ) in Fig. 11(b) and the reduced charged impurity scattering contribute to the positive photoconductivity. On the other hand, the increased phonon scattering following n η in Fig. 11(b) leads to negative photoconductivity. Figs. 11(c) and (d) compare σ (ω, τ 1 ) and ∆σ (ω, τ 1 ) at τ 1 = 2.0 ps under different pump fluences of F ab = 0.56, 0.77, and 1.02 µJcm −2 . The σ 1 (ω, τ 1 ) value at F ab = 0.56 µJcm −2 decreases slightly, compared with that at the equilibrium, thereby exhibiting negative photoconductivity that changes into positive photoconductivity as F ab increases. σ 2 (ω, τ 1 ) at F ab = 0.77 µJcm −2 increases/decreases compared with that of the equilibrium below/above ω/2π 2 THz. As a result, ∆σ 1 (ω, τ 1 ) at F ab = 0.77 µJcm −2 exhibits negative and the minimum value around ω/2π = 1.5 THz, and the corresponding ∆σ 2 (ω, τ 1 ) decreased with the frequency, thereby exhibiting the zero crossing around ω/2π = 2 THz, which effectively reproduced the observed Drude-Lorentz-type ∆σ (ω, τ 1 ) between ω/2π = 0.6 and 3.0 THz. A simple Drude-type photoconductivity with negative ∆D and positive ∆Γ cannot reproduce the Drude-Lorentz-type ∆σ (ω, τ 1 ) as demonstrated in Ap-pendix F. Therefore, the frequency dependent ∆D (ω, τ 1 ) and ∆Γ (ω, τ 1 ) for F ab = 0.77 µJcm −2 in Figs. 11(e) and (f), respectively, which stem from the temporal variation of the carrier distribution and momentum relaxation rate during THz probing, are necessary for the observed photoconductivity. This result suggests that the strong dependence of the Drude-Lorentz type ∆σ (ω, τ 1 ) on the different atmospheric environment would be caused by the change of charged impurity scattering through the gas adsorption that may affect the carrier and/or dielectric screening.
In summary, we studied the intraband optical conductivity of hot carriers at quasi-equilibrium in photoexcited graphene based on a combination of the iterative solutions for the BTE and comprehensive temperature model. The proposed method enables us to consider extrinsic effects such as charged impurity scattering and the intrinsic effect of optical phonon scattering on the hot carrier THz conductivity, with a reduced computational cost compared with the full BTE solution.
In the examples of photoexcited graphene with different Fermi energies, we demonstrated the temporal evolution and frequency dependence of the photoconductivity during the cooling process exhibit a strong dependence on the temporal variation of carrier distribution and momentum relaxation rate. Our method provides a quantitative analysis of hot-carrier THz conductivity with intrinsic and extrinsic microscopic parameters such as the electron-phonon coupling and charged impurity concentration, which are important for the development of future graphene optoelectronic devices. In n-type doped graphene, the carrier concentration n 0 at T e = 0 K in the conduction band is expressed as a function of the Fermi energy ε F : At a finite T e , the total carrier concentration n (T e ) in the conduction band is given by Here, µ (T e ) is the finite temperature chemical potential, and n (T e ) is expressed as the sum of n 0 and the thermally (A3) As a result, n 0 is expressed as a function of the temperature-dependent chemical potential µ (T e ) as Because n 0 is constant and given by Eq. (A1) for T e ≥ 0, the temperature dependence of the chemical potential can be obtained by numerically inverting Eq. (A4) 26 . Figure 12(a) presents the chemical potential of the charge carriers in graphene as a function of temperature.
In p-type doped graphene, the temperature dependence of the chemical potential can be obtained similarly. In the Boltzmann theory, assuming a constant carrier relaxation rate, the Drude weight of the 2D-MDF exhibits a unique temperature dependence as given by 26,[74][75][76] .
The temperature-dependence of D (T e ) for different Fermi energies are shown in Fig. 12(b).

Appendix B: Calculation of relaxation rate using RTA
We demonstrate the numerical results of the calculation of the inelastic scattering by an intrinsic non-polar optical phonon in monolayer graphene. Figures 13(a)-(c) present S out λ as a function of the carrier energy for graphene with |ε F | = 0.01, 0.15, 0.43 eV, respectively, calculated using Eqs. (7), (8), and (12). For comparison, we show the momentum relaxation rate calculated using RTA 47,52 . In the case of nondegenerate semiconductors, the Maxwell-Boltzmann distribution n MB = 1/Exp[(ε − µ)/k B T ] with µ = 0 eV is used instead of the FD distribution for the calculation of the momentum relaxation rate by non-polar optical phonons with mode η, which is given by where the top sign applies for phonon absorption, and the bottom sign applies for phonon emission. Furthermore, D η is the deformational potential of the optical phonon mode, η. We used D 2 η = D 2 η F for the RTA calculation. The total optical phonon scattering by RTA is then expressed by Figure 13(d) presents the energy dependence of τ −1 op , which is larger than S out λ . This is because the scatteringangle dependence is not included for the calculation of τ −1 op . Moreover, we demonstrate elastic or quasielastic carrier scattering by charged impurities 54,66,72,[77][78][79] , weak scatterers 68,[80][81][82][83][84][85][86] , and acoustic phonons 66,[77][78][79]83,84,[87][88][89][90][91] considered in the iterative calculation of the THz conductivity. For elastic scattering, the RTA under the low field condition is valid 47,52 , and the momentum relaxation rate in graphene is equal to the relaxation rate of the distribution function.
Charged impurity scattering plays an important role in the carrier transport and intraband conductivity in undoped or lightly doped graphene on a substrate. The carrier scattering by charged impurities at the graphenesubstrate interface limits the carrier mobility significantly. The relaxation rate τ −1 i (ε λk ) for the charged impurities under RTA is expressed by Here, n i is the charged impurity concentration, p = |k − k | is the scattering wave vector, and v i (p) = e 2 /2 ave p is the Fourier transform of the 2D Coulomb potential in an effective background lattice permittivity ave = (1 + s ) 0 /2, given by the average static dielectric constant of the vacuum and substrate. Here, 0 is the vacuum permittivity. Furthermore, (p, T e ) is the static electronic dielectric function of graphene calculated by the random-phase approximation (RPA), and it is responsible for the screening effect.
where Π (p, T e ) is the graphene irreducible finitetemperature polarizability function expressed by  Table I. in which g = 4 is the spin and valley degeneracy factor, A is the area of the system, and f λ (k ) is the carrier distribution function. Moreover, f λ (k) f 0 (ε λk ) is used in the calculation of Π (p, T e ) because we consider the small perturbed distribution g (ε λk ). The temperature dependence of the charged impurity scattering arises from the temperature-dependent carrier screening of the Coulomb disorder, which depends on the D (T e ) 55 . Figures 13(e)-(g) present the energy dependence of the momentum relaxation rate via charged impurities with n i = 1.0 × 10 12 cm −2 .
The possible physical origins of weak scatterers are ripples and point defects. The relaxation rate τ −1 s (ε λk ) for the weak scatterers is expressed by 68 where ρ s is the resistivity of the weak scatterers and ranges from 40-100 Ω 83-86 . The acoustic phonon scattering is treated as quasielastic, and therefore, the relaxation rate τ −1 ac (ε λk ) can be expressed as cited in 87 .
Here, D ac is the acoustic deformation potential, which ranges from 10-30 eV 66,[92][93][94][95] , and v ph = 2.0 × 10 4 ms −1 is the sound velocity in graphene. In this study, we used ρ s = 100 Ω and D ac = 30 eV in the calculation of the energy dependence of τ −1 s and τ −1 ac , respectively, as illustrated in Fig. 13(h). As a result, the elastic scattering term ν el in Eq. (13) reads In the RTA, the collision term in Eq. (1) under the low field is expressed in the form Here, τ −1 (ε λk ) is the relaxation rate of the distribution function. Using τ −1 (ε λk ), the intraband optical conductivity in graphene is given by 80 ∂ε λk .

(B10)
Appendix C: Calculation of temperature evolution of hot-carrier quasi-equilibrium state The temperature model was used in previous studies 15,37,45,51,69,96 ; however, we modified the model to include the carrier energy relaxation such as the SC carrier cooling process and the optical phonon emission and absorption process. This can be described using the coupled rate equations in Eqs. (22) to (25). Moreover, we used the specific heat capacity C considering the temperature dependence of µ(T e ) of 2D-MDF as shown in Fig.12(b). Here, C = C c + C v is the sum of the specific heat of the electrons in the conduction bands and valence bands given by 97 where U c and U v are the thermal kinetic energy of the electrons in the conduction band and in the valence band, respectively. Note that C (T e ) , ε λk , and ν (ε λk ) are functions of the Fermi velocity, v F , which is renormalized by electron-electron, electron-phonon, and electron-plasmon coupling and depends on the carrier concentration because of the carrier screening effect in graphene 61, . However, the screening constant of the electron-electron interaction has become a subject of considerable debate [123][124][125][126][127] . For simplicity, we used v F = 1.1 × 10 6 ms −1 in the numerical calculation. Figure 14  potential. Here, α = 16.3 × 10 −8 nJK −3 cm −2 is twice the value of α = 8.14 × 10 −8 nJK −3 cm −2 , which only considers the electron heat capacity in the conduction band, as reported in Ref. 96 . As can be observed, C (T e ) for graphene with |ε F | = 0.01, 0.15 and 0.43 eV exhibits different behaviors from the quadratic dependence below T e = 2, 000 K. These deviations are attributed not only to the finite Fermi energy, but also to the temperaturedependent chemical potential µ (T e ).
R Net η = R η − G η in Eq. (22) denotes the total balance between the optical phonon emission and absorption rate.
(C3) Figure 15 compares the T e dependence of the energy exchange rates R η ω η and G η ω η for optical phonon emission and absorption in graphene at different Fermi energy and phonon temperatures T η . For T η < T e , R η ω η is larger than G η ω η , as indicated in Figs. 15(a)-(c); the hot-carrier energy is transferred to the optical phonons. The total balance R Net η ω η = R η ω η − G η ω η becomes zero when T η is equal to T e , as illustrated in Figs. 15(d)-(f). However, the energy relaxation of the carriers and phonons is further driven by the optical phonon decay rate τ −1 ph caused by the anharmonic phonon-phonon interaction and SC process.
J sc in Eq. (22) denotes the energy loss rate for the SC process, which are disorder-mediated electron-acoustic phonon scatterings and take place via a three-body collision involving a carrier, a defect, and an acoustic phonon. Although energy relaxation by acoustic phonon scattering is essential for low-energy carriers, the small Fermi surface and momentum conservation severely constrain the phase space of the acoustic phonon-scattering process. As a result, acoustic phonon scattering becomes an inefficient cooling channel. However, additional momentum exchanges with disorders enable acoustic phonons to use a much wider phase space, thereby enabling a larger dissipation of energy from the hot carriers. Thus, SC becomes an efficient cooling pathway. According to Ref. 43 , J sc is given by where g ac = D ac /2ρv 2 ph is the electron-acoustic phonon coupling constant, N (ε F ) is the density of states at the Fermi energy per one spin or valley flavor, k F is the Fermi wave number, l is the mean free path of the short-range weak scatterers, and T ac is the acoustic phonon temperature, which is assumed to remain unchanged from the equilibrium state. Further, J sc can be expressed as a function of the carrier mobility µ m = σ DC /n c e for shortrange weak scatterers, as discussed in Ref. 45 : For simplicity, we used the carrier mobility of µ m = 4, 800 cm 2 V −1 s −1 for the SC carrier cooling throughout the calculation. The T e dependence of J sc is plotted in Fig. 14(b). In Eqs. (23) to (25), R Net M,n = R M,η − G M,η denotes the total balance between the optical phonon emission and absorption rate per number of phonon modes.
(C7) In the above equation, M − η (ε λk ) and M + η (ε λk ) are the numbers of η phonon modes (per unit area) that participate in the carrier-phonon scattering of the emission and absorption processes for carriers that have energy ε λk , respectively.
In this case, a Γ = 1 for Γ-LO and Γ-TO phonons, and a K = 2 for K phonons. The factor of a K = 2 for the K phonons represents the degenerate phonon valleys at the K and K points. Fig. 16 compares the T e dependence of the energy-exchange rates R M,η and G M,η for the optical phonon emission and absorption in graphene at different Fermi energies and optical phonon temperatures T η , that demonstrate similar behaviors to those of R η ω η and G η ω η Appendix D: Calculation of transmission change of THz probe pulses When applying the standard thin-film approximation 128 , the THz-amplitude transmission coefficient of monolayer graphene with complex conductivity σ(ω) on a substrate with a thickness of d at normal incidence is expressed by the ratio of the incident wave E i (ω), and the transmitted wave E t (ω).

(D4)
The normalized negative transmission change −∆E t (τ 2 , τ 1 ) /E t (τ 2 ) as a function of the probe-trigger delay τ 2 at τ 1 is expressed by Here, E t (τ 2 ) is the transmitted THz field through the graphene sample without pump fluence. Figure 17 plots ∆E t (τ 1 ) /E 0 ≡ ∆E t (0, τ 1 ) /E t (0) of the undoped graphene with ε F = −0.01 eV and n iH calculated using the THz probe with pulse durations of 2τ p = 300 and 500 fs. The corresponding central frequencies of which are ω/2π = 2.2 and 1.2 THz, respectively. In this case, −∆E t (τ 1 ) /E 0 for 2τ p = 500 fs exhibits a faster decay time than that of 2τ p = 300 fs, which indicates the difference of the decay time of ∆σ (ω, τ 1 ) between 1.2 and 2.2 THz. Appendix E: Convergence of electrical and optical conductivity by iterative calculation The convergence of the iterative calculation including inelastic scattering is presented. Figure 18(a) shows the dependence of the DC conductivity σ DC in the highly doped graphene under a constant electric field 100 Vcm −1 on the number j of iterations using Eq. (14) considering only optical phonon scattering at different temperatures. At T e = 295 K, σ DC converged very rapidly, and even at j = 1, the error σ c DC − σ j DC /σ c DC < 10 −2 , where σ c DC and σ j DC are the converged DC conductivity and that at the jth iteration, respectively. However, the convergence rate became slower and the error of the first iteration increased as the temperature increased, which could be attributed to the broader carrier distribution and increased scattering rate by optical phonons at higher temperatures. Figure 18(b) presents the converged g (ε λk ) of the heavily doped graphene under an electric field of 100 Vcm −1 calculated by Eq. (14). The g (ε λk ) value with a clear peak around ε F = −0.43 eV at T e = 295 K spread and decreased with an increasing temperature that resulted in a reduction in σ c DC at high temperature. Moreover, we determined the convergence of the iterative solution for the time-dependent process given by Eq. (19) in which we considered the inelastic and elastic scatterings by optical phonons, acoustic phonons, charged impurities, and weak scatterers for the highly doped graphene. The momentum relaxation rates and S out λ used in the calculation at T e = 295 and 1, 500 K are illustrated in Fig. 19. Fig. 20(a) depicts the temporal evolution of the current density J(t) when an electric field with a step function shape was applied. The rising time of J(t) decreased with an increasing Ω s in Eq. (19) and converged for Ω s > 250 THz, which was more than 10 times the momentum relaxation rate (approximately 25 THz) at |ε λk | = |ε F | = 0.43 eV. The J(t) value converged to 0.22 Acm −1 , yielding the DC conductivity σ c DC = 28.5G 0 . Further, we calculated the temporal variation in the current density induced by applying the THz-probe pulse with 2τ p = 300 fs as indicated in Fig. 20(b). The temporal waveforms of J(t) and the corresponding σ(ω) are plotted in Figs. 20(c) and (d), respectively. Whereas the temporal waveform of J(t) appears to be converged around Ω s = 250 THz, the convergence of σ(ω) is strongly dependent on the frequency. Below ω/2π = 0.2 THz, the convergence is achieved even at Ω s = 10 THz, and the convergence value of the DC conductivity of σ(ω) = 28.5G 0 is equal to σ c DC as obtained from the calculation by the step function-type electric field presented in Fig. 20(a). However, σ(ω) at a higher frequency requires a larger Ω s to be converged. In the numerical calculations explained in the main manuscript, we set Ω s = 200 and 1000 THz for ∆E t (τ 1 ) /E 0 and ∆σ (ω, τ 1 ), respectively.