Hot carrier dynamics and electron-optical phonon coupling in photoexcited graphene via time-resolved ultrabroadband terahertz spectroscopy

Electron-electron (e-e) interaction is known as a source of logarithmic renormalizations for Dirac fermions in quantum field theory. The renormalization of electron--optical phonon coupling (EPC) by e-e interaction, which plays a pivotal role in hot carrier and phonon dynamics, has been discussed after the discovery of graphene. We investigate the hot carrier dynamics and the EPC strength using time-resolved ultrabroadband terahertz (THz) spectroscopy combined with numerical simulation based on the Boltzmann transport equation and comprehensive temperature model. The large negative photoconductivity and the non-Drude behavior of THz conductivity spectra appear under high pump fluence and can be attributed to the temporal variation of the hot carrier distribution and scattering rate. We successfully estimate the dimensionless EPC matrix element of the $A_1^{\prime}$ optical phonon mode near the $\mathbf{K}$ point as $\lambda_{\mathbf{K}} \approx$0.09 from the fitting of THz conductivity spectra and temporal evolution of transient THz reflectivity, which is slightly larger than the prediction of the renormalization group.


I. INTRODUCTION
Hot carrier effects are regarded as insightful in studying many-body interactions in condensed matter, and play a crucial role in the operation of electronics and optoelectronic devices. For this reason, they have been investigated extensively in both metals and semiconductors 1,2 . The rise of graphene had offered new opportunities for this research field because the carriers thereof are 2D massless Dirac fermions (MDFs) with a linear energy dispersion. This fact has promoted graphene as an attractive platform for hot carrier physics and various applications  . Electron or hole relaxation mainly involves non-radiative electron-electron (e-e) and electron-phonon scatterings, depending on the excitation energy. Electron-electron interaction is dominant at high energy, redistributes the electrical or optical power within the electron gas, and builds up a hot carrier population. Electron-phonon interaction operates on a longer time scale to equilibrate the electron and phonon temperatures, and to cool the hot carriers 24 .
Hot carrier effects play a significant role in the optoelectronic properties of photoexcited graphene, in which the photocarriers are excited at high energies. The subsequent relaxation drives the working efficiency of optoelectronic devices. In this respect, spectroscopic investigations such as pump probe spectroscopy 25 and angleresolved photo-electron spectroscopy 26,27 of hot carriers complement transport studies. Optical pump terahertz (THz) probe spectroscopy (OPTP) is a powerful tool for investigating the hot carrier dynamics of graphene because it probes the intraband optical conductivity dominated not only by the hot carrier distribution, but also the carrier scattering process in contrast to optical pump optical probe spectroscopy. Extensive studies using OPTP [28][29][30][31][32][33][34][35][36][37][38][39][40][41] have revealed the unusual behaviors of graphene hot carriers, which undergo positive and negative changes in the intraband optical conductivity with non-Drude type frequency dependence. The negative change observed in heavily doped graphene is an indicative of enhanced carrier scattering and reduced Drude weight in quasi-equilibrium hot carrier state with a single chemical potential owing to ultrafast recombination of photoexcited carriers. However, most of these works were performed by THz probe with the relatively narrow band (1-3 THz) which was not sufficient for capturing the whole spectrum of non-Drude type conductivity and their results have been interpreted using the framework of the phenomenological model 31,33,36,37,40 . Such a phenomenological analysis for the narrow band spectra is not sufficient to understand the hot carrier and phonon dynamics quantitatively and to derive the microscopic parameters. Theoretical studies have been conducted by incorporating the microscopic theory based on the density matrix formalism 41 or Boltzmann transport equation (BTE) 42,43 .
The electron-optical phonon coupling (EPC) strength is a crucial factor that makes it difficult to understand the hot carrier and phonon relaxation process by numerical studies. The density functional theory (DFT) calculations demonstrated that only three strongly coupled optical phonon (SCOP) modes contribute signifi-cantly to the inelastic carrier scattering in graphene 44,45 . The first two relevant modes are associated with the G peak of the Raman spectrum and the highest optical branches at Γ (the E 2g mode) with the energy of ω Γ = 196.0 meV, which split into an upper longitudinal optical (LO) branch and a lower transverse optical (TO) branch near Γ. Owing to their long wavelengths, these phonons scatter electrons within one valley. Moreover, it is essential to take into account the highest optical branch of the zone boundary phonon ω K = 161 meV at the K point (the A 1 mode). This mode is responsible for intervalley processes and associated with the D and 2D peaks of the Raman spectrum. In Refs. 44,[46][47][48] , the EPCs D 2 η F for dominant optical phonon modes η (Γ LO , Γ TO , K) were defined as the average on the Fermi surface of the matrix element D η λkλ k of the Kohn-Sham potential, differentiated with respect to the phonon displacement. The EPC for LO and TO modes at the Γ point had D 2 Γ F = 45.6 (eVÅ −1 ) 2 , which was in good agreement with experimental results 49 . However, the EPC value at the K point has been debated 47,49-52 because it is renormalized by the e-e interaction and is affected by the presence of the substrate owing to the dielectric screening effect 53 . The amount calculated by DFT with generalized gradient approximation was D 2 K F = 92.0 (eVÅ −1 ) 244 . However, a GW calculation, which considers the e-e interaction by approximating the self-energy using the product of the Green function and screened Coulomb potential, but neglects the vertex corrections, yielded D 2 K F = 193 (eVÅ −1 ) 247, 54 .
In this work, we investigate the hot carrier dynamics in photoexcited heavily doped graphene on a polyethylene terephthalate (PET) substrate using an OPTP and estimate the EPC strength at the K point via a numerical simulation based on the combination of BTE and comprehensive temperature model 43 . Owing to the small change in the Drude weight of heavily doped graphene and negligible contribution of charged impurity and surface optical phonon (SOP) of PET substrate, the rise and relaxation dynamics of the SCOP are effectively captured by the transient THz reflectivity change measured by ultrabroadband THz probe. Using the renormalization group analysis, the obtained dimensionless EPC at K point is discussed and compared with those determined by other techniques.

II. SIMULATION METHOD AND RESULTS
In this section, we present a numerical simulation of the THz conductivity and the transient THz reflectivity measured by the OPTP experiment according to the following procedures. After photoexcitation, photoexcited carriers are quickly recombined and their energy is redistributed within electron gas forming hot carrier state in quasi-equilibrium with single chemical potential. A number of cooling pathways for hot carriers by inelastic scattering have been proposed such as SCOPs 41,55,56 , acoustic phonon 57-60 , SOP of substrate 61 . As we will explain later, the contribution of SOP and its coupled mode with plamons can be neglected by selecting the substrate with low polarizability and small phonon energy ω  . Effect of acoustic phonon on hot carrier cooling is considered by the supercollision process and the acoustic phonon occupation is assumed to remain unchanged from the equilibrium state in the picosecond time scale after photoexcitation 27 . Therefore, we use comprehensive temperature model to calculate the temporal evolutions of the temperature for hot carriers in quasi-equilibrium and the occupations for three dominant SCOP modes. Thereafter, the iterative solution of BTE 43 is used to calculate the intraband complex conductivity of the hot carriers in THz region. Because interband transition is forbidden at a THz probe energy of ω THz < 2ε F , the THz conductivity of doped graphene is dominated by the intraband transition. This scheme enables us to reduce the computational cost substantially compared to the calculation of the full solutions of coupled graphene Bloch equation and BTEs for hot carriers and hot phonon modes in 2D momentum space.

A. THz conductivity calculation
The iterative solution of the BTE for obtaining the steady-state and dynamical conductivity of semiconductors was introduced in Refs. 65,66 and was subsequently modified for 2D MDF in graphene 43 . The temporal evolution of the carrier distribution is described by the BTE under a time-dependent electric field, which is expressed as 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, and E(t) is the electric field of the THz probe pulse. ∂f λ (k, t)/ ∂t| c is the collision term that describes the change in the distribution function via carrier scattering. We consider the intrinsic carrier scattering mechanism by the optical and acoustic phonon modes 67-76 and the extrinsic mechanism by the charged impurities 67,68,[77][78][79] , and weak scatterers 69,70,[80][81][82][83][84][85] . For spherical bands under a low field E, the general solution of Eq. (1) is approximately provided by the first two terms of the zone spherical expansion.
where f 0 (ε λk ) = 1/ [exp {(ε λk − µ (T e )) /k B T e } + 1] is the Fermi-Dirac distribution for the corresponding equilibrium electron distribution at the electron temperature T e . ε λk = ± v F |k| (ε 1k ≥ 0 and ε −1k ≤ 0 for the conduction and valence bands, respectively) is the electron energy within the Dirac approximation of the graphene energy-band structure 86 , and v F is the Fermi velocity. In this expression, µ (T e ) is the temperature-dependent chemical potential of the 2D MDF 36,78,87 and is illustrated in Fig. 3(a). g (ε λk , t) is the perturbation part of the distribution, and α k is the angle between E and k. In Eq. (1), the collision term is given by while accounting for the scattering of the electrons with dominant optical phonon modes η, in C η λλ , including both the intraband (λ = λ ) and interband (λ = λ ) processes with 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 transition rate by the optical phonon modes, η, between states (k , λ ) → (k, λ) and (k, λ) → (k , λ ), respectively. P η kλk λ is expressed by which accounts for the phonon emission and absorption, given by Here, D η λkλ k is the EPC matrix element defined in Ref. 44 , k = k ± q, and q is the wave vector of the optical phonons. ρ = 7.6 × 10 −7 kgm −2 is the area density of graphene and A is the area of the graphene sample, whereas ω η and n η are the angular frequency and occupation of the optical phonons, respectively.
The carrier-scattering rates that are obtained by the optical phonons in Eq. (6) account for the phonon emission and absorption. For small q and k , the EPC matrix elements |D η λλ k | 2 for the Γ LO , Γ TO , and K phonons are expressed by 44,88 Here, θ k,q denotes the angle between k and q, θ k ,q denotes the angle between k and q, and θ k,k denotes the angle between k and k . In the case of Γ LO and K phonons, the plus sign refers to the interband processes, and for Γ TO phonons, it refers to the intraband processes. In Eq. (3), the elastic term C el λ (k) is calculated using the elastic scattering rate P s λkλk 43 . The index, s, refers to the different elastic scattering modes by weak scatterers, and charged impurities, which are characterized by resistivity of the weak scatterers ρ s , and charged impurity concentration n i , respectively. The reported ρ s ranges from 40-100 Ω 69,70,84,85 . Interactions with acoustic phonons are treated in a quasi-elastic and included in C el λ (k). Different electron-acoustic phonon coupling models have been proposed to extract the effective coupling constant J a from experimental data for graphene which ranges 10-30 eV 70,[73][74][75]81,[89][90][91][92][93][94][95][96][97][98][99] . A firstprinciple study suggests that the gauge-field contribution is more important than the screened deformation potential 100,101 .
The iterative solution of g j (ε λk ) = g (ε λk , t j ) is provided by Here, E j = |E (t j )| and k = |k | are the magnitudes of the electric field and wavevector, respectively. Ω s is known as the self-scattering rate, and 1/Ω s is the time increment between successive iterations, and S in λ and S out λ are the net in-and out-scattering rates for inelastic scattering, respectively. Furthermore, ν el is the total relaxation rate by the elastic scattering mechanisms. The sequence {g j (ε λk )} yields f λ (k , t j ) versus time when Ω s is sufficiently large compared to S out λ + ν e .

B. Temperature model of hot carriers
The hot carrier intraband optical conductivity σ(ω, τ 1 ) in the cooling process can be calculated from f λ (k , t j ), which is obtained by substituting the hot carrier and three dominant optical phonon temperatures (T e (t j ), T η (t j )) into Eq. (8) in the iteration process. Here, τ 1 is the pump probe delay. We employ the coupled rate equations for a comprehensive temperature model that describe the temperature evolutions of the electron temperature T e and optical phonon occupations n η by photoexcitation: In this case, I ab represents the pump intensity absorbed in graphene sample during laser irradiation, considering the multiple reflections inside the substrate with dielectric constant (ω pump )=2.4 for the pump wavelength and saturable absorption (SA) effect in graphene. 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 indicates the energy loss rate for the supercollision carrier-cooling process 58,102 . R Net M,η = R M,η − G M,η denotes the total balance between the optical phonon emission and absorption rate per number of phonon modes that participate the carrier scattering. In calculations of R Net η and R Net M,η , we include the scattering angle dependence of the |D η λk,λ k | 2 in Eq.(7) which have not been considered in the temperature model used in the previous study 43,[102][103][104] . Moreover, n η0 represents the phonon occupation near the Γ and K points, respectively, in equilibrium at room temperature, whereas τ ph is the phenomenological optical phonon decay time to other phonon modes via the phonon-phonon interaction caused by lattice anharmonicity 105 . The effective optical phonon temperatures are calculated by inverting the Bose-Einstein distribution function, n η = 1/(e ωη/kBTη − 1). The formula and temperature dependence of C, R Net η , J sc and R Net M,η can be found in Ref. 43 and Section SIII in Supplemental Material (SM).
The optical pump pulse is absorbed by interband transition and the absorption coefficient for free-standing graphene at the normal incidence is α inter = πα = 0.23 % under sufficiently weak pump condition, where α is the fine structure constant. However, the SA effect in graphene under the intense pump fluence 108-111 should be considered. The SA is a nonperturbative, nonlinear optical phenomenon that depends on the pump power as well as the temperature and Fermi energy. Based on the theory by Marini et al. 111 , we derived the formula of I ab considering the SA and multiple reflections inside the substrate at the oblique angle of incidence for the temperature calculation in the experimental condition (see Section IV in the SM): where I 0 (t) is the envelope function of the incident pump pulse, which is assumed to have hyperbolic secant form, I 0 (t) = (F 0 /2τ pump ) sech 2 (t/τ pump ). In this case, F 0 is the incident fluence and 2τ pump is the pump pulse duration. I n (t + n∆T ) = (F n /2τ pump ) sech 2 ((t + n∆T )/τ pump ) represents the pump pulse by the n-th multiple reflection of the incident pump pulse inside the substrate, where F n is the fluence and n∆T is the round-trip time for the n-th reflection pump pulse in the substrate. A s * ij (F 0 /2τ pump ) is the absorption coefficient including the carrier temperature dependence of the SA effect at the interface of layer i/graphene/layer j when the pump pulse excites the graphene from layer i (see Fig. S1 of Section I in the SM). In this model, the SA is characterized by the inelastic carrier relaxation time τ ie . The pump intensity dependence of the interband absorption coefficient α inter for the free standing graphene and A s * ij (F 0 /2τ pump ) for the graphene on substrate can be seen in Figs.S2 and S3 of Section IV in the SM.

C. Simualtion for graphene on PET substrate
In the simulation, the carrier scattering by SOPs of substrate are not included, while the SOPs play crucial roles for the carrier dynamics in graphene on polar substrate 90,[112][113][114][115][116] . The square of EPC matrix element between SOP and carries is proportional to Here, g SO = βe 2 ω SO /2 0 , ω SO is the angular frequency of the SOP, ε 0 is permittivity of vacuum and d is the equilibrium distance of the graphene sheet from the substrate surface. q is the angular wavenumber of the surface phonon, q s is the Thomas-Fermi screening constant of the 2D carriers and where s and ∞ are the low and high frequency dielectric constant, respectively. β is a measure of the polarizability of the dielectric interface. For example, in crystalline SiO 2 ( s = 3.9, ∞ = 2.5) has two SOP modes at ω s1 = 60.0 meV, ω s2 = 146.5 meV, with β 1 = 0.025 β 2 = 0.062, respectively. These values correspond to g SO1 = 0.14 (eV 2Å−1 ) and g SO2 = 0.82 (eV 2Å−1 ) and are enhanced by roughly 50 % in conventional SiO 2 glass with ∞ = 2.1. As a result, the temperature dependence of carrier transport is dominated by SOP scattering in graphene on polar substrate such as SiO 2 and HfO 2 90,114 . The energy loss rate of hot carrier by SOP modes is given as R NET SO ∝ g SO ω SO so that the large ω SO also affect the hot carrier dynamics significantly. The dispersion relation of SOP modes can be altered by the coupling of plasmon and SOP in doped graphene. These effects change significantly the hot carrier dynamics and makes the simulation more complex leading to hindering the estimation of EPC at K point. Therefore, in this study, we select graphene sample on a PET substrate which has the low polarizability ( s = 3.0, ∞ = 2.54) owing to the polar low frequency vibrational modes around 10 meV 106 . The g SO = 0.029 (eV 2Å−1 ) of PET is small and decreases significantly in doped graphene by carrier screening effect. The R NET SO between carriers and SOP of PET is expected to be smaller by 3 orders of magnitude than SiO 2 and makes the negligible contribution on hot carrier cooling and THz conductivity. Furthermore, the small static dielectric constant s = 3.0 of a PET substrate provides weak dielectric screening with an expected larger renormalization effect on the EPC by e-e interaction 53 .
We investigated the effect of the EPC on the hot carrier dynamics of photoexcited graphene on the PET substrate for different Fermi energies ε F . The parameters used in the simulation are summarized in Table I. Eq. (7). A comparison between Figs.1 (a) and (b) reveals that the rise and relaxation dynamics of the hot carrier and optical phonon temperatures depend significantly on D 2 K F . At D 2 K F = 703 (eVÅ −1 ) 2 , T K followed T e more rapidly and increases up to 1800 K much higher than T Γ LO/TO , indicating that substantially more hot carrier energy is mainly transferred into the K phonon owing to the stronger EPC. As a result, the maximum T e for D 2 K F = 703 (eVÅ −1 ) 2 becomes lower than that for D 2 K F = 193 (eVÅ −1 ) 2 . Figure 1(c) presents the D 2 K F dependence of the transient reflection change ∆E r (τ 1 )/E 0 calculated from the σ(ω, τ 1 ) using the THz probe pulse with 2τ p = 300 fs. The sign of ∆E r (τ 1 )/E 0 remains negative indicating the negative photoconductivity as varying the D 2 K F . The peak value of |∆E r (τ 1 )/E 0 | increases monotonically as D 2 K F increases and effectively reflects the enhancement of T K . Figure 2 depicts the simulation results on the lightly doped graphene with |ε F | = 0.15 eV. Although the same phonon decay time τ ph = 1 ps is used, the relaxation time of T e of the lightly doped graphene is longer than that of the heavily doped graphene owing to the weaker R Net η originated from the small density of state at the Fermi energy ε F . The sign of ∆E r (τ 1 )/E 0 indicated in Fig. 2  ∆E r (τ 1 )/E 0 exhibits positive photoconductivity, which is transformed into negative photoconductivity as D 2 The different behaviors in ∆E r (τ 1 )/E 0 between the heavily and lightly doped graphene can be understood by considering the temperature dependence of the Drude weight D (T e ) of the graphene 2D MDF, which is the oscillator strength of free carrier absorption and plays a crucial role in carrier screening. As can be observed in Fig. 3(a), the chemical potential µ(T e ) of graphene 2D MDF decreases with T e , leading to the unique temperature dependence of D (T e ) according to ε F 36,43,117-119 . In the case of a constant carrier relaxation rate, D (T e ) is expressed as The D(T e ) of the undoped graphene with |ε F | = 0.01 eV in Fig. 3(b) increases linearly with T e , yielding positive photoconductivity. However, D(T e ) of the heavily doped graphene with |ε F | = 0.43 eV decreases slightly as T e increases and exhibits the minimum at around T e = 2000 K, contributing to the negative photoconductivity below T e = 3000 K. At temperatures below 3000 K, the maximum change in D(T e ) is only 13% and the temperature dependence of THz conductivity change is mainly dominated by of the carrier scattering with the SCOPs. In the lightly doped graphene, D(T e ) increases significantly above T e = 1000 K and the contributions of D(T e ) and the carrier scattering with SCOPs to the photoconductivity compete with one another resulting in the positive and negative photoconductivity depending on T e and D 2 K F . We also investigated the effect of the charged impurity on the hot carrier dynamics in the heavily and lightly doped graphene because the charged impurity is one of the dominant scattering mechanism in graphene on substrate 67,68,77,90 . Figure 1(d) shows the ∆E r (τ 1 )/E 0 of the heavily doped graphene is almost unaffected by charged impurity scattering owing to the strong carrier screening effect. Here, the effective coupling constant J a of acoustic phonon is selected so that the DC conductivity is almost equal as shown in Table I. However, the ∆E r (τ 1 )/E 0 of the lightly doped graphene in Fig. 2(d) changes significantly by the presence of the low charged impurity concentration n i = 0.17 × 10 12 cm −2 , indicating a crossover from the negative ∆E r (τ 1 )/E 0 to the positive one and the reduction of the carrier scattering due to the enhanced carrier screening effect. Therefore, the information of the accurate charged impurity concentration is required to derive the D 2 K F from ∆E r (τ 1 )/E 0 of lightly doped graphene. These findings indicate that heavily doped graphene is suitable for the determination of D 2 K F from ∆E r (τ 1 )/E 0 .

III. EXPERIMENTAL RESULTS
The graphene sample (Graphene Platform Corporation) that was examined in this study was prepared using chemical vapor deposition. The single-layer graphene (area: 10 mm×10 mm) was transferred to a PET substrate. Raman scattering measurements confirmed the single-layer thickness of the sample and their low defect density. The equilibrium THz conductivity of the sample at room temperature (T 0 = 295 K) was characterized by ultrabroadband THz time domain spectroscopic ellipsometry (THz-TDSE) (see Section I in the SM for details), which enabled the broad Drude peak to be captured directly by measuring the ratio of the reflection coefficient r p (ω)/r s (ω) in the frequency range between 1.0 and 20 THz 120 , as illustrated in Fig. 4. The fitting of the THz conductivity spectrum obtained from r p (ω)/r s (ω) by the Drude model allows us to determine the Drude weight D 0 and carrier relaxation rate Γ 0 for the equilibrium state at room temperature T 0 = 295 K accurately. We estimated D 0 = 1.36 × 10 3 G 0 and Γ 0 = 21.4 meV, respectively. Here, G 0 = 2e 2 /h is the quantum conductance. The corresponding Fermi energy is |ε F | = 0.43 eV, indicating that the sample is heavily doped and suitable for estimating the EPC strength. The carrier concentration n c at T e = 0 K and the DC conductivity at T 0 were estimated as n c = 1.1 × 10 13 cm −2 and σ DC = 20G 0 , respectively, where we used v F = 1.1 × 10 6 m s −1 considering the carrier and dielectric screening effect in heavily doped graphene on PET substrate 121 . Figure 5(a) presents the optical setup of the reflectiontype OPTP used in the experiment. Amplified femtosecond laser pulses (1kHz repetition rate, 785 nm center wavelength) are used to generate ultrabroadband THz probe pulses from laser-excited air plasma 122 . Spolarized pump pulses with a pulse duration of 220 fs are loosely focused and excited the graphene sample at an incident angle of θ = 60 • and the created hot carrier state was probed by s-polarized THz pulses with a pump probe time delay τ 1 . The temporal waveforms of the reflected THz probe pulses are measured by air breakdown coherent detection, which detects the second harmonic generation of the trigger pulse induced by the THz electric field 123 . Figure 5(b) depicts the temporal waveforms of the THz probe pulse reflected from the photoexcited graphene. When the pump fluence is increased, the peak amplitude of THz probe decreases slightly, indicating negative photoconductivity. The ratio of the reflection coefficient r s (ω, τ 1 )/r s (ω) of graphene with and without pump fluence F 0 = 200 µJ cm −2 calculated by Fourier transformation of the THz waveforms at different τ 1 values, as plotted in Fig. 5(c), decreases and then recovers to the equilibrium reflecting the rise and subsequent relaxation process of the hot carrier dynamics, and this was used for the calculation of σ(ω, τ 1 ) (see Section II in the SM for details). Figure 5(d) presents the fluence dependence of ∆E r (τ 1 )/E 0 , which exhibits multiple negative peaks around τ 1 = 0.2, 1.4, 2.3 ps owing to the multiple reflections inside the PET substrate. As F 0 increases, the peak height ∆E r (τ 1 )/E 0 increases but it exhibits saturation behavior with an increased relaxation time.    less than half of that at the equilibrium (gray curve), indicating a significant increase in the carrier scattering by SCOPs at high temperatures. It is found that σ(ω, τ 1 ) for D 2 Γ F by the DFT (black curve) and GW (blue curve) calculations can not reproduce the observed negative photoconductivity, even if the SA effect is not considered. On the other hand, σ(ω, τ 1 ) for D 2 Γ F = 703 and 946 (eVÅ −1 ) 2 show the larger deviation than that for D 2 Γ F = 450 (eVÅ −1 ) 2 . Figures 7(a)-(c) depict the comparison of ∆E r (τ 1 )/E 0 between the experiment and calculations, which is significantly dependent on D 2 K F and the pump fluence F 0 . For D 2 K F by the DFT and GW calculations, the peak height and temporal evolution of ∆E r (τ 1 )/E 0 differ significantly from the experimental values and the higher values D 2 K F = 450-946 (eVÅ −1 ) 2 are required to reproduce the ∆E r (τ 1 )/E 0 . By comparing σ(ω, τ 1 ) and ∆E r (τ 1 )/E 0 with the calculation in Figs. 6 and 7, we estimated D 2 K F ≈ 450 (eVÅ −1 ) 2 and τ ie = 116 fs, at which the calculation (blue curves) best fits the experimental results. In this case, the obtained τ ie = 116 fs corresponds to the saturated pump intensity I s = 1.0 and 1.7 × 10 8 W cm −2 for α inter and A s * 12 respectively, which is slightly smaller than the reported value in Ref. 111,124 . Figure 8 presents the temporal evolution of T e and T η calculated for D 2 K F = 92.0 and 450 (eVÅ −1 ) 2 under the pump fluence F 0 = 200 µJ cm −2 indicating that hot carrier and phonon dynamics are significantly dependent on the EPC. For D 2 K F = 92.0 (eVÅ −1 ) 2 as shown in Fig. 8(a), the hot carrier temperature increases beyond T e =3000 K, and T K followed T e slowly owing to  Temperature (K) the weak EPC and reaches up to T K ≈ 1500 K. In this high temperature range, the carrier scattering by optical phonons is dominant and the Drude weight D(T e ) makes the positive contribution to σ(ω, τ 1 ) in contrast to the carrier scattering. The competition of these factors leads to broader peaks of ∆E r (τ 1 )/E 0 for DFT (black line) in Fig. 7(c) than those of T η in Fig. 8(a). For D 2 K F = 450 (eVÅ −1 ) 2 as seen in Fig. 8(b), the hot carrier temperature increases up to only T e ≈2000 K and T K follows T e rapidly and reaches up to T K ≈ 1400 K owing to the SA effect and strong EPC. In this case, D(T e ) makes the same contribution to σ(ω, τ 1 ) as the optical phonon scattering, resulting in sharper peaks of ∆E r (τ 1 )/E 0 and a successful reproduction of the experimental results. Furthermore, the frequency dependence of σ(ω, τ 1 ) at τ 1 = 0.1 ps in Fig. 6 deviates from the simple Drude model as F 0 increases. This originates from the rapid temporal variation in the carrier temperature and scattering rate during the THz probing time following the photoexcitation, and the calculation with D 2 K F = 450 (eVÅ −1 ) 2 effectively reproduces the observed large negative photoconductivity with non-Drude behavior. This indicates that most photoexcited carriers are recombined and the quasi-equilibrium hot carrier state is almost established at τ 1 = 0.1 ps owing to the strong Auger recombination in the heavily doped graphene, as reported in Ref. 26 . The parameters used in the calculation are displayed in Table II.

IV. DISCUSSION
Based on the fitting of ∆E r (τ 1 )/E 0 by the calculation considering the EPC, we estimated the phenomenological phonon decay time due to lattice anharmonicity as τ ph = 0.3, 0.45 and 0.57 ps for F 0 = 50, 100 and 200 µJ cm −2 , respectively. Refs. 125,126 reported longer τ ph = 0.8-1.5 ps for graphene on a SiO 2 substrate. However, these values were determined from the simple fitting of transient absorption or anti-stokes Raman intensity by exponential function and do not consider the EPC. The simple fitting of ∆E r (τ 1 )/E 0 with exponential curve results in τ ph = 1.15-1.5 ps which are comparable to the reported values. The theoretical study reported the phonon decay time τ ph ≈ 3.5 and 4.5 ps for Γ and K phonon by only considering the anharmonicity of lattice in graphene without substrate 105 . Therefore, the obtained τ ph indicates the dominant contribution of substrate for the op-   tical phonon decay channel.
The dimensionless coupling constants λ Γ and λ K for the optical phonons near the Γ and K points, respectively, are useful for comparing the EPC strengths determined from various experiments and calculations, which are defined as 49 In the above, M ≈ 2.00 × 10 −26 kg is the mass of the carbon atom and A u.c. ≈ 5.24Å 2 is the unit-cell area. F 2 Γ and F 2 K have the dimensionality of a force and are the proportionality coefficients between the change in the effective Hamiltonian and lattice displacement along the corresponding phonon mode. Subsequently, the matching rules are expresses as F 2 Γ = 4 D 2 Γ F and F 2 K = 2 D 2 K F . Note that λ K is subject to Coulomb renormalization, which implies that λ K is dependent on the electronic energy scale, such as the electron energy, Fermi energy, or temperature T, whichever is larger: .09 using Eqs. (13) and (14). Figure 9 presents the flow of λ Γ and λ K for different background static dielectric constants av = (1 + s )/2 = 1, 2, and 5 calculated by solving the renormalization group equation in Ref. 53 , which sum up the leading logarithmic corrections and go beyond the Hartree-Fock approximation. The bare values of the dimensionless EPCs λ Γ = 0.031 and λ K = 0.038 were selected to satisfy the relation λ Γ /λ K = ω K /ω Γ and to reproduce the experimental value λ Γ = 0.031 127 . The renormalization group analysis demonstrated that, although λ Γ was almost constant, λ K was strongly dependent on the energy scale as well as av . The obtained λ K (ε F ) = 0.09 slightly larger than the calculated value of λ K (ε F ) = 0.073. According to the ratio of the λ(ω K )/λ K (ε F ) = 1.21 for av = 2 in Fig. 9, we obtained λ(ω K ) = 0.11, which is a factor of 3.2 larger than the DFT value λ K (ω K ) = 0.034. Raman studies 127,129-132 using a field effect transistor based on the polymer electrolyte ( av = 5) reported λ Γ =0.028 and 0.031 from the ratio of the area between G and the 2D peak, which were comparable to λ Γ = 0.028 by the DFT calculation of g 2 Γ F using Eq. (14). However, λ K (E L /2) ranged between 0.05 and 0.15 as seen in Fig. 9, where E L is the laser excitation energy (for a typical Raman measurement E L /2 ∼ 1eV). The corresponding λ K (ω K ) are estimated as 0.063 and 0.19. The lower limit value is comparable to the calculated λ K (ω K ) for av = 5. Although Raman spectroscopy is a powerful tool for the determination of λ K (ω K ) as well as λ Γ (ω K ), it requires the accurate estimation of the gate capacitance of FET device which are not required in OPTP experiments.

V. CONCLUSION
In conclusion, we investigated the EPC of the optical phonons near the K point of heavily doped graphene on PET substrate and the hot carrier dynamics using a combination of the time-resolved THz spectroscopy and numerical simulations. The hot carrier dynamics in heavily doped graphene on PET substrate is less sensitive to the extrinsic charged impurity and surface polar phonons of the substrate and is dominated by the electron-optical phonon interactions. According to the quantitative analysis based on the BTE and comprehensive temperature model considering the SA effect on pump fluence, the ∆E(τ 1 )/E 0 value can be used for the determination of the EPC in graphene. The estimated D 2 K F ≈ 450 (eVÅ −1 ) 2 indicates the strong renormalization by e-e interaction and the corresponding dimensionless coupling constant λ K (E F ) ≈ 0.09 slightly larger than the calculation by the renormalization group theory. The extension of the simulation model for the undoped or lightly doped graphene on various substrate requires the accurate estimation of charged impurities and surface polar phonons of the substrate is a future issue that will be important to the development of graphene optoelectronic devices.  In this section, we explain the calculation procedure of the THz conductivity σ(ω THz ) of graphene on the substrate from the ratio of the complex reflection coefficient (r p (ω THz )/r s (ω THz )) for the p-and s-polarized THz waves measured by THz time domain spectroscopic ellipsometry (THz-TDSE) 1 . According to the standard thin-film approximation, the reflection coefficients of graphene on a substrate for p-and s-polarized THz wave are given by 2 (SI.1b) In the above, Z 0 = 376.7 (Ω) is the vacuum impedance and θ 1 = 60 • is the incidence angle of the THz wave. Furthermore, i (ω THz ) is the dielectric constant of layer i, as indicated in Fig.S1. From Eq. (SI.1), σ(ω) is expressed as By substituting (r p /r s ) measured by THz-TDSE into Eq. (SI.2), σ(ω THz ) can be determined.  In this section, we present the calculation procedure of the hot carrier THz conductivity σ(ω THz , τ 1 ) of photoexcited graphene at the pump probe delay τ 1 from the reflection-type OPTP measurement. The reflection-type OPTP measures the ratio of the complex reflection coefficient X s (ω THz , τ 1 ) = r s (ω THz , τ 1 )/r s (ω THz ) of graphene with and without photoexcitation. The reflection coefficient for the s-polarization of graphene with complex conductivity σ(ω THz ) at an incident angle of θ 1 is expressed by Eq. (SI.1b). Similarly, the THz-amplitude reflection coefficient for the s-polarization of graphene with hot carrier complex conductivity σ(ω THz , τ 1 ) on the substrate at an incident angle of θ 1 for the pump probe delay τ 1 is expressed by (SII.1) Using Eqs. (SI.1b) and (SII.1), we obtain where B and B are provided by Eqs. (SI.3c) and (SI.3d), respectively, and r s (ω THz ) is calculated using the equilibrium σ(ω THz ) obtained by THz-TDSE. We can obtain the σ(ω THz , τ 1 ) by substituting X s (ω THz , τ 1 ) into Eq. (SII.2).

III. RATE EQUATIONS FOR TEMPERATURE MODEL
In this section, we present the derivation of the hot carrier recombination and generation rate by optical phonon emission and absorption process, respectively, used in the temperature model. The Hamiltonian of electron-phonon interaction H cp is Here, V ep is the potential of the electron-phonon interaction, c † k (c k ) is the creation (annihilation ) operator with carrier wave vector k, b † q (b q ) is the creation (annihilation ) operator with phonon wave vector q. From Fermi's golden rule, the carrier transition rate from k to k by the emission and absorption of the Γ LO phonon or Γ LO phonon with the energy of ω Γ are given by Here, D Γ λkλ k 2 is the square of the EPC matrix element. For small q and k, the EPC matrix elements are D Γ . ρ is the mass density, A is the area of graphene sample, ε λk = λ v F |k| is the energy of 2D MDF and λ = ±1 is the band index. The upper and lower sign corresponds to the optical phonon emission and absorption process, respectively. The corresponding hot carrier recombination and generation rate per unit area including both intra-and inter-band transitions are written as (SIII.3b) Here, N (ε λk ) = 2|ε λk |/π( v F ) 2 is the density of state of 2D MDF. Furthermore, the electron distribution function f λ (k) can be replaced by Fermi-Dirac type distribution f 0 (ε λk , T e ) for hot carriers in quasi-equilibrium. Similarly, the hot carrier recombination and generation rate by K-phonon with the energy of ω K are given by Using Eqs. (SIII.3)-(SIII.4), the total balance between the optical phonon emission and absorption rate is given by R Net η = R η − G η . In Eq. (9), R N et M,η = R η − G η denotes the total balance between the optical phonon emission and absorption rate per number of phonon modes.
Here, M − η (λk) and M + η (λk) are the number of η-phonon modes (q) per unit area that participate the phonon emission and absorption processes for carries state (λ, k), respectively.
In this case, a Γ = 1 for Γ-LO and Γ-TO phonons, and a K = 2 for K phonon. The factor of a K = 2 represents the degenerate phonon valleys at the K and K points. Using Eqs. (SIII.5), the total balance between the optical phonon emission and absorption rate per number of phonon modes is given by R Net

IV. PUMP POWER INJECTED INTO GRAPHENE SAMPLE CONSIDERING SATURABLE ABSORP-TION
In this section, we present the derivation of the pump intensity F ab injected into the graphene sample, considering the multiple reflections inside the substrate and the saturable absorption (SA) effect. The SA is an extreme nonlinear phenomenon that consists of the quenching of the optical absorption under high-intensity illumination. Following Marini et al. 3 , we introduce the derivation of saturable absorption coefficient α inter in graphene. Thereafter, we explain the derivation of the absorbed pump intensity F ab by graphene on the substrate at an oblique incidence angle using α inter .
We study the response of a single electron in graphene under an in-plane x-direction applied field E(t) = E 0 e −iωtx . The extended Bloch equations describing the temporal variation in the interband coherence ρ k and population difference n k in photoexcited graphene are as follows: where ξ = (eτ ie E 0 / k) sin φ and ω ± = ω ± 2ω 0 . The macroscopic interband current density depending on the light intensity I 0 = (c/2π)|E 0 | 2 at the electronic temperature T e is determined by , and . (SIV.6b) Subsequently, by expressing the integral over the reciprocal space in polar coordinates, the following is obtained: Using the interband current, the interband absorption coefficient is determined as the ratio of the time-averaged absorbed power over an optical cycle to the incident intensity I 0 : Although the above results were obtained under the CW illumination conditions, these are also applicable to commonly used optical pulses that have a large duration compared to the optical period. Taking into account the SA for the interband transition by the pump irradiation, the transmission and reflection coefficients of the s-polarized pump pulse incident on the system of layer i/graphene/layer j from layer i, as illustrated in Fig. S1, are calculated by (SIV.9b) In this case, the pump pulse irradiates the graphene from layer i with the incidence angle of θ i and transmits it to layer j with the angle θ j . Moreover, γ ij is the correction factor. Although α inter (ω, I 0 , T e ) is appropriate for the case in which the optical pump pulse excites the suspended graphene at the normal incidence angle, the saturation behavior will change when graphene on a substrate is excited by a pump pulse at an oblique incidence angle, where the injected pump power becomes smaller by a factor of γ ij . The corresponding transmittance and reflectance are determined by Using Eq. (SIV.10), the absorption of the pump pulse by the graphene layer is provided by The correction factor γ ij is calculated by the ratio of the absorption coefficient , (SIV. 12) and can be determined self-consistently. Using the converged γ * ij , the transmittance, reflectance, and absorption coefficients in the experimental condition are obtained by The envelope function of the pump pulse considering the n th multiple reflections inside the substrate is given by In this case, I 0 (t) represents the incident pump pulse, which is assumed to have hyperbolic secant form I 0 (t) = (F 0 /2τ pump ) sech 2 (t/τ pump ), where F 0 is the fluence and 2τ pump is the pulse duration. I n (t) = (F n /2τ pump ) sech 2 (t/τ punmp ) represents the n th reflection of the incident pump pulse and F n is the fluence of the n th reflection pulse. ∆T is the time delay owing to one round trip in the substrate. Using Eq. (SIV.13) and I 0 = F 0 /2τ pump , F n for n ≥ 1 is obtained by In the above, R s 23 is the reflectance of the pump pulse incident at the substrate (layer 2) /N 2 purged (layer 3) interface from the substrate (α inter (I 0 ) = 0 in Eq. (SIV.13b)). Using Eqs. (SIV.13), (SIV.14), and (SIV.15), the absorbed pump intensity F ab (t) is determined by  Figure S3 (a) and (b) present the saturated pump intensities I s for α inter and A s * 12 , respectively, where I s is defined as α inter (I s ) = (1/2)α inter (0) and A s * 12 (I s ) = (1/2)A s * 12 (0). Figure S4 shows the absorbed pump fluence in graphene with |ε F | =0.15 and 0.43 eV, calculated using Eq. (SIV.16).       Pump-probe delay t 1 (ps)