Bilayer graphene (BLG) is a kind of graphene composed of two coupled honeycomb-like carbon layers^[1]. The strong interaction of the Dirac–Fermions in BLG has driven new insight into the strong correlation phenomena such as superconductor and correlated insulator phases in the magic-angle twisted BLG systems^[2–6]. During the phase diagram investigating the Bardeen–Cooper–Schrieffer (BCS) state and Bose–Einstein condensate (BEC) crossover in the magic-angle graphene system, the large critical temperature (T_c) and Fermi temperature ($T_F$) make the weak electron-phonon (e-p) coupling system in BCS theory not applicable to describe the emergent phenomena in the twist BLG (TBG)^[2]. Therefore, it is crucial to investigate the effect of the layer-layer interaction on the electron-electron (e-e) and e-p coupling in BLG in the ultrafast time scale. Kar et al. studied the response of the photoexcited hot carriers to terahertz (THz) fields in low and moderately doped BLG deposited on a quartz substrate, which shows that the relaxation dynamics is dominated by carrier disorder interaction along with carrier acoustic phonon interactions, as compared to charge-impurity Coulomb interaction or carrier surface phonon interaction^[7]. Recently, by means of high-pressure treatment, the ultrafast hot electron cooling process in BLG was interpreted as a strong interlayer shear phonon coupling, which could accelerate the relaxation dynamics^[8]. However, the impact of strong interlayer coupling on the THz response, which is more sensitive to the electrons near the Fermi energy, had been seldom discussed^[9].

In the case of doped monolayer graphene (MLG), the ultrafast THz photoconductivity is reduced by the intraband carrier-carrier thermalization near the Fermi energy, which is different from the interband transition origin compared with the optical pump optical probe and angle-resolved photoelectron spectroscopy (ARPES)^[8,10]. Similarly, THz is sensitive to the strong coupling induced flat band and the possible gap opening process near the Dirac cone, which is associated with the interlayer e-e and e-p interactions^[6,11,12]. Therefore, to disentangle the electron dynamics near the Fermi surface of graphene, the visible-infrared pump THz probe spectroscopy (VIPTP) is introduced to investigate the ultrafast thermalization and hot carrier cooling process in both MLG and BLG. The strong interlayer interaction is proved to be able to accelerate the ultrafast cooling dynamics in BLG. However, the subpicosecond thermalization efficiency is found to be reduced in BLG. Furthermore, the carrier scattering time obtained by the Drude fitting indicates that the strong layer-layer coupling in the BLG is probably correlated with the band renormalization process. Therefore, the interlayer strong coupling plays different roles in the thermalization and hot carrier cooling processes in BLG. Specifically, for thermalization, the interlayer interaction-induced less efficient thermal process could be attributed to the carrier-carrier scattering enhancement rather than the interlayer e-p coupling cooling mechanism, which is possible due to the A-A stacking (Turbostratic stacking)^[13]. In contrast, in the case of cooling process, the interlayer phonon scattering, which is related to the A-B stacking (Bernal stacking), the interlayer interaction could help to boost the cooling process^[8]. Our investigation would not only propose the different physical origin between the interlayer interaction and the intralayer e-p coupling, but also give new insight into the ultrafast investigation of the strong coupling in the BLG system.

The VIPTP is driven by a 1 kHz Ti:sapphire regenerative amplifier with 800 nm central wavelength and 35 fs pulse duration and then through an optical parameter amplifier to change the center wavelength to 190–2600 nm with the beam size of $0.25{\mathrm{cm}}^2$. The THz pulses are generated by optical rectification and detected by electro-optic sampling in a pair of 1 mm thick, (110)-oriented ZnTe crystals. The MLG and BLG were grown on a ${\mathrm{SiO}}_2$ substrate by chemical vapor deposition (CVD), respectively. As shown in Fig. S1 of Supplementary Material, the static THz transmission spectrum is shown to be 6% absorption for MLG and 8% absorption for BLG. Figure 1(c) shows the typical Raman spectra of MLG and BLG excited by a 532 nm laser, respectively. The 2D peak and G peak are related to the energy band structure of graphene^[14]. The G peak position of both samples indicates that MLG and BLG have approximately the same Fermi level^[13]. Compared with the MLG, the electronic band structure splits, both the conduction band and valence band are composed of two parabolas, and there are four possible double resonance scattering processes in BLG of the 2D peak, which have the $\mathrm{FWHM}=66{\mathrm{cm}}^{-1}$ in BLG and $\mathrm{FWHM}=45{\mathrm{cm}}^{-1}$ in MLG^[14–16]. As shown in Fig. 1(c), the Raman spectrum of MLG features a much higher 2D peak intensity than that of the G peak ($I_{2\mathrm{D}}/I_G\approx 2$)^[17]. Moreover, the peak width of the 2D peak increases, and the intensity ratio $I_{2\mathrm{D}}/I_G$ decreases in the case of BLG, which reconfirms the bilayer crystal structure.

The pump-induced photoconductivity ($\mathrm{\Delta }\sigma$) is probed by a THz pulse by monitoring the change in the transmitted electrical field ($\mathrm{\Delta }E=E_{\mathrm{pump}}-E_0$) following the optical pump as a function of pump-probe delay. The measurement is based on the principle that $\mathrm{\Delta }E$ is proportional to $\mathrm{\Delta }\sigma =\frac{1+n}{Z_0}(\frac{1}{1+\mathrm{\Delta }T/T_0}-1)$^[9], where, $n=1.95$ is the refractive index of the fused silica substrate, and $Z_0=377$ is the free space impedance. The time-resolved signals $\mathrm{\Delta }\sigma$ of MLG and BLG films dependent on the different pump fluence under pump photon energy of 0.78 eV, 1.13 eV, and 1.91 eV are shown in Fig. S2 (in Supplementary Material). We use the deconvolution equation of the e-e interaction in BLG with the deconvolution fitting, $$|\mathrm{\Delta }\sigma |=A\times \text{Conv}(e^{-\frac{t}{{\tau }_{\text{decay}}}},{\tau }_{\text{pulse}})+B\times \text{Conv}(e^{-\frac{t}{{\tau }_{\text{decay}}}},{\tau }_{\text{slowrise}}),$$(1)where $\text{Conv}(e^{-\frac{t}{{\tau }_{\text{decay}}}},{\tau }_{\text{pulse}})$ is the convolution of $e^{-\frac{t}{{\tau }_{\text{decay}}}}$ with a Gaussian pulse of pulse width $\tau$, ${\tau }_{\text{pulse}}$ is the pulse width used, and the calculated effective conductivity drop time can be obtained with ${\tau }_{\text{rise}}=(A+B)/[(\frac{A}{{\tau }_{\text{pulse}}})+(\frac{B}{{\tau }_{\text{slowrise}}})]$; this is how long it takes for the initial light excitation to produce a peak negative THz conductivity. As shown in Fig. 2(a), with the Gaussian deconvolution fitting of the thermalization process in BLG and MLG, the e-e thermalization time in BLG is about 0.75 ps; in contrast, the e-e process takes about 0.35 ps, which indicates a higher electron thermalization in BLG with the help of interlayer interaction. For the cooling process, BLG is faster than MLG. The photoexcited THz conductivity increases with increasing pump fluence and photon energy due to the presence of the carrier multiplication process^[18]. In contrast to the monolayer one, the carrier-carrier scattering-induced transient conductivity of BLG is reduced, which could be modulated either by the carrier mobility or the e-p coupling efficiency.

In order to disentangle the e-p and optical phonon-acoustic phonon interaction effect between the graphene layers, the mono-exponential fitting is presented in Fig. 3^[19]: $$-\mathrm{\Delta }\sigma (\mathrm{\Delta }t)=A·e^{-\frac{t}{{\tau }_{\text{decay}}}}+C,$$(2)where ${\tau }_{\text{decay}}$ is the relaxation dynamics of the lattice temperature. $C$ represents the longitude thermal conductivity diffusion process, which could be regarded as a plateau in our condition. As shown in Figs. 3(a)–3(c), the relaxation lifetime of MLG increases with the pump fluence increasing from $4\mathrm{\mu J}/{\mathrm{cm}}^2$ to $160\mathrm{\mu J}/{\mathrm{cm}}^2$, while in BLG the relaxation time of the hot carrier cooling process increases with the pump fluence at first and then reaches saturation after $40\mathrm{\mu J}/{\mathrm{cm}}^2$. For MLG, cooling is slower because electrons first cool via optical phonons and then there can be a hot-phonon bottleneck due to optical-to-acoustic phonon cooling^[19]. However, the bottleneck is less severe because there is more phonon coupling in BLG due to the interlayer interaction^[20]. As shown in Fig. 3(d), the cooling process of BLG is significantly more efficient than in MLG^[11], which indicates that an extra layer-layer interaction-induced cooling channel is possible to affect the cooling process in BLG^[21]. The interlayer phonon scattering, which is caused by the interlayer interaction, could help boost the cooling process. Considering that the Fermi energy and the degree of the disorder do not change obviously between these two samples, as represented by the Raman spectra shown in Fig. 1(c)^[15], the intrinsic origin such as the carrier-carrier strong interaction, the phonon-phonon interaction assisted interband transition, and the interlayer shear mode acoustic phonon coupling should be taken into account^{[2,8,10,11,22]}.

For the primary concern, the negative THz onductivity arises from the carrier-carrier thermalization process, which reduced the conductivity of the carriers near the Fermi surface; therefore, the measurement of the thermalization efficiency could represent the carrier-carrier Coulomb interaction strength in MLG and BLG^[18]. Therefore, we introduce a parameter to describe the degree of sub-linearity. Figures 4(a)–4(c) show the THz conductivity peak value $|\mathrm{\Delta }\sigma |$ of MLG (black) and BLG (red) at 0.78 eV, 1.13 eV, and 1.91 eV, respectively. The peak THz conductivities of MLG and BLG are fitted by the power exponent function. The degree of deviation of the thermalizing sub-linearity $\eta$ can be fitted from the following formula^[23]: $$\mathrm{\Delta }\sigma (\text{peak})=A{F}^{\eta },$$(3)where $\mathrm{\Delta }\sigma (\text{peak})$ is the peak value of the photo-induced THz conductivity at different pump wavelengths (photo energy), $F$ is the pump fluence, $\eta$ describes the degree of sub-linearity of the thermalization, and $\eta$ increases as the thermalization efficiency decreases. The degrees of deviation of the thermalizing sub-linearity are summarized in Table 1.

The degree of deviation of the thermalizing sub-linearity $\eta$ of the hot electron in ideal MLG is predicted to be 0.5 due to the electron capacity C being in direct proportion to the electron temperature T^[24]. In the case of CVD graphene with high mobility, the thermal efficiency is found to be dependent on the doping level and thus the Fermi–Dirac distribution near the Fermi surface^[25]. As shown in Figs. 4(a)–4(c) and Table 1, the sub-linearity index $\eta (\mathrm{BLG})$ is larger than $\eta (\mathrm{MLG})$ and monotonically increases with the pump photon energy, which indicates less photo-induced thermal efficiency in BLG. As described in Figs. 4(d) and 4(e), for BLG, the strong interaction between layers may affect the band structure near the Dirac point, resulting in the decrease of the interband transition-induced intraband thermalization efficiency. Therefore, it is reasonable to conclude that the layer-layer interaction would reduce the thermal efficiency and thus the carrier-carrier interaction in the initial subpicosecond thermal process^[26]. Moreover, the sublinear photon energy dependence of the sub-linearity index $\eta$ in graphene demonstrates that the thermal efficiency is still related to the electron-optical phonon coupling, while the thermal efficiency decrease indicates that the origins for the thermal and cooling processes in both graphene samples are different^[7,27–29].

Despite the e-p coupling, the strong layer interaction affected layer electron momentum scattering time should be investigated. We measured the 2D photoconductivity of two samples at the pump-probe delay time of 0 and 2 ps, respectively, as shown in Figs. 5(a)–5(d). To understand the conductivity dispersion relation for different samples, we fit the experimental conductivity by the Drude model, which is often employed for the transport of charge carriers in graphene^[9], $$\mathrm{\Delta }\tilde{\sigma }(\omega )=\frac{i\mathrm{\Delta }D}{\pi (\omega +i\mathrm{\Gamma })},$$(4)where $\mathrm{\Gamma }$ is the average scattering rate for momentum changing collisions of charge carriers. The Drude weight $\mathrm{\Delta }D$ characterizes the magnitude of the reduced conductivity. The carrier scattering times at different pump-probe delays are $0.023\pm 0.059\mathrm{ps}$ and $0.094\pm 0.057\mathrm{ps}$ in MLG, $0.477\pm 0.069\mathrm{ps}$ and $0.835\pm 0.133\mathrm{ps}$ in BLG, respectively. Due to the relationship below, $\sigma =n|e|{\mu }_e$, ${\mu }_e=\frac{e\tau }{m^{*}}$, one can find that the Drude conductivity is in direct proportion to the scattering time and the carrier density. As the Fermi energy and the pump fluence are fixed, the time dependent conductivity change should come from the carrier momentum scattering time $\tau$, which is related to the electronic band structure near the Dirac cone^[10,26]. Therefore, the thermalization of the Drude-like Dirac–Fermion should be modulated by the A-A stacking of the interlayer tunneling potential and thus the band structure^[3,13]. In addition, the e-e scattering time increase could possibly be induced by the asymmetric band structure of BLG^[30,31]. The momentum scattering time difference modulated by the interlayer interaction further confirms the different origins of the strong interlayer coupling modulation of the thermalization and the hot carrier cooling process in BLG.

In conclusion, the time-resolved THz conductivity dynamics of MLG and BLG have been studied by optical pump THz probe spectroscopy with different photon energies and fluences. By analyzing the thermal and cooling dynamics of the optical-induced negative THz conductivity, we find that the interlayer strong coupling process plays different roles in the ultrafast thermal and relaxation process of THz conductivity. In addition, the interlayer interaction of BLG could accelerate the hot carrier cooling process by means of the phonon-phonon coupling process due to the A-B stacking clusters in the illumination zone of BLG, and the carrier-carrier scattering time of the Dirac–Fermion reduction could be assigned to the band structure modulation induced by the A-A stacking interlayer tunneling process. Finally, the different carrier scattering times changed with the pump delay in both the MLG and BLG are obtained by utilizing the Drude model, and the momentum scattering time in BLG is much larger than that in MLG, indicating that the band structure was modulated by the layer-layer interaction. Our work will provide new insights for the application of ultrafast nonequilibrium heating and cooling pathways, and the improvement of multi-layer graphene-based nano optoelectronic devices^[32].