Angle-tunable intersubband photoabsorption and enhanced photobleaching in twisted bilayer graphene

Van der Waals heterostructures obtained by artificially stacking two-dimensional crystals represent the frontier of material engineering, demonstrating properties superior to those of the starting materials. Fine control of the interlayer twist angle has opened new possibilities for tailoring the optoelectronic properties of these heterostructures. Twisted bilayer graphene with a strong interlayer coupling is a prototype of twisted heterostructure inheriting the intriguing electronic properties of graphene. Understanding the effects of the twist angle on its out-of-equilibrium optical properties is crucial for devising optoelectronic applications. With this aim, we here combine excitation-resolved hot photoluminescence with femtosecond transient absorption microscopy. The hot charge carrier distribution induced by photo-excitation results in peaked absorption bleaching and photo-induced absorption bands, both with pronounced twist angle dependence. Theoretical simulations of the electronic band structure and of the joint density of states enable to assign these bands to the blocking of interband transitions at the van Hove singularities and to photo-activated intersubband transitions. The tens of picoseconds relaxation dynamics of the observed bands is attributed to the angle-dependence of electron and phonon heat capacities of twisted bilayer graphene.

Van der Waals heterostructures composed of stacked single-layer crystals represent a promising platform for novel electronic, photonic and spintronic devices [1][2][3][4][5][6][7][8]. Wide tunability of the optoelectronic properties is achieved by both selection of the stacked materials and control of the relative rotation between the crystal registries [9,5,[10][11][12], i.e. the twist angle. In the case of two layers of graphene, strong interlayer coupling drives the hybridization of the electronic states and can lead to the formation of moiré super-lattices for finite twist angle.
The rotational misalignment of twisted bilayer graphene (tBLG) enables to explore also new properties as topological valley transport [13], superconductivity [14] at the magic twist angle θ ~ 1.1°, or the topologically protected 1D chiral states appearing at the moiré domain walls [15]. The Dirac cones associated to the two layers in tBLG are no longer aligned in the energy-momentum landscape and two saddle points appear at their lowenergy intersections. The saddle points contribute to van Hove singularities (vHs) in the density of states, which occur at energies EvHs that depend on the twist angle [16]. The vHs are detected as a peak in the optical conductivity of tBLG [17][18][19] and they have been related to the increase of hot photoluminescence (PL) for resonant excitation [20]. These vHs appear at much lower energies than in single-layer graphene (SLG) [16,[21][22][23][24] and allow angle controlled visible and near-infrared (NIR) absorption enhancement. This property has been exploited in a tBLG based photo-detector to selectively increase the photocurrent generation at the vHs [25].
The changes in crystal symmetry, electronic band structure and optical conductivity come along with modifications of the relaxation pathways available to charge carriers driven out of equilibrium. The carrier dynamics in SLG is well-known as it has been intensively studied by means of ultrafast transient absorption [26][27][28][29][30][31] (TA), time-resolved photocurrent [32,33] and photoemission spectroscopies [34][35][36]. Excitation with femtosecond light pulses produces out-of-equilibrium distributions of electrons and holes in the conduction and valence bands, which thermalize in ~ 50 fs via Coulomb two-body scattering to hot Fermi-Dirac (FD) distributions [30,31]. The relaxation of the hot carriers proceeds through scattering with strongly coupled optical phonons (SCOP) till a common temperature is reached in about 200 fs. The cooling of the thermalized electron-phonon bath is finally driven by the lattice and achieved in less than 10 picoseconds via anharmonic decay of the hot phonons. In defected SLG, hot carriers can release their excess energy also by defect-mediated emission of acoustic phonons (supercollision mechanism), which occurs on a ps-time scale dependent on substrate and defect densities [31,33,[37][38][39][40][41][42][43]. The hot carrier distribution inhibits interband absorption over a broad energy range due to Pauli blocking, so that in the TA spectra of SLG a decreased absorption, or photo-bleaching (PB), signal is detected.
In Bernal stacked bilayer graphene, together with the PB band, Limmer et al. [44] observed a photoinduced absorption (PA) with a single decay time of 5 ps, about twice as long as the PB signal. The PA has been attributed to transitions between subbands formed within the conduction and the valence bands, named intersubband (ISB) transitions, which are absent in SLG since it lacks a subband structure. Also time-and angle-resolved photoemission spectroscopy (TARPES) has demonstrated long-living (up to 10 ps) population inversion in bilayer graphene [45,46] attributed to the presence of a small gap. Previous studies in tBLG by TA microscopy [47,48] and TARPES [49] have explored the twist-angle dependence and shown long decay times attributed to bound excitons and to states reached by two-photon absorption [48,50].
Here we study the out-of-equilibrium optical properties of tBLG as a function of the twist angle by femtosecond TA microscopy with high sensitivity, which allows to explore low excitation regimes in which angle-dependent contribution from ISB transitions emerges. Specifically, we investigate heterostructures with the vHs located in the NIR and twist angle θ in the range of ~ 6°-8°. The tBLG samples are prepared by chemical vapor deposition as described in the Methods and they have tens of μm linear dimension. A lattice rotation in the real space by an angle θ corresponds to a rotation of θ in the phase-space of the hexagonal Brillouin zones (BZs), as sketched in Fig. 1(a-b). Accordingly, the Dirac points of the two layers (see Figure 1(c)), KA and KB, no longer coincide and are displaced by a momentum ΔK = 2K sin(θ/2), where K = 4π/3a and a ≈ 0.246 nm is the SLG lattice parameter [51]. The interlayer interactions perturb the Dirac cones of each layer causing the formation of minigaps with the same energy in the conduction and valence bands at the vHs, as is schematically shown in Fig. 1(c).
Figure 1(d) shows the optical microscope image of one of the investigated samples, named stack 1, measured in reflection. It consists of a large single-crystal SLG, with an inner brighter region, where a second graphene layer is stacked to form a bilayer single crystal. A first determination of the mutual orientation of the lattices in the bilayer domain is obtained by extrapolation of the hexagonal habits of the outer/inner growth regions, corresponding respectively to the 1 st /2 nd layer [52]. In stack 1, the hexagons are rotated by an angle θ = 7°± 1°.
The angle determination from the optical images is affected by a large uncertainty related to the definition of the single-crystals shape. Twist angle is usually measured by field high-resolution transmission electron microscopy [53,54], scanning tunneling microscopy [16] and lattice resolution atomic force microscopy [55]. Here we adopt an all-optical, far-field approach, based on measuring the hot-photoluminescence excitation (hot-PLE) spectrum under femtosecond laser excitation. As previously observed [20], due to the enhancement of pump absorption, a peak is expected in the bilayer emission when the laser is tuned at the vHs. In Fig. 1 We obtain the value θ = 7.31° ± 0.04° for stack 1, in good agreement with the evaluation from the optical image inspection. The hot PL map in Fig. 1(e) for laser excitation tuned at 1.4 eV, close to EvHs, reports indeed a ≈ 2 times higher emission intensity at the tBLG compared to the SLG, uniform across the bilayer domain.
Extra peaks have been observed in visible, IR and UV Raman spectra of tBLG compared to SLG. Accordingly, we characterize the heterostructure by micro-Raman spectroscopy using a visible laser tuned at 633 nm. Figure   1(f) shows the Raman spectra at the SLG (red line) and tBLG (blue line) regions. The G-band peak at 1584 cm -1 , the 2D peak at 2648 cm -1 and an additional peak at 1621 cm -1 are observed for tBLG, while for SLG these band are observed at 1588 cm -1 (G) and 2645 cm -1 (2D).
In the bilayer, the G-band peak has a two times higher intensity, as shown in the inset of Fig. 1(f). The enhancement of the G-band and the extra peak, La, testify the formation of a new band structure by interlayer coupling. It is thus an indication of the high crystalline quality of the samples and of a clean interface between the two monolayers forming the tBLG. The La peak above the G-band position for the bilayer is the mode activated by intralayer electron-phonon scattering process from the unfolded longitudinal optical phonon (LO) [56] branch.
The measured La frequency of 1621 cm -1 is consistent with the theoretical predictions [56] for a tBLG with θ ~ 7°.
The LO dispersion is expected to have a non-monotonic dependence on the twist angle, with a maximum frequency at θ ≈ 10° but only slowly varying in the range of interest, thus making it unpractical to precisely define the twist angle only from the Raman data.
We apply ultrafast TA microscopy to study the photophysics of tBLG samples as function of the twist angle. A NIR pump pulse with photon energy Epump is used to photo-excite the sample off-resonance with the vHs (Epump = 0.8 eV stack 1 and Epump = 1.19 eV for stacks 2-3), while a second tunable pulse, the probe, monitors the relaxation of the induced hot carriers' distributions as function of pump-probe time delay t. Pump and probe beams are collinearly coupled to a confocal microscope to achieve a ≈ 1 μm spatial resolution. The differential transmission (ΔT/T) is measured as a function of the probe photon energy Eprobe, using two distinct high sensitivity TA microscopes described in the Methods, guaranteeing a spectral coverage from 0.85 eV to 1.65 eV. Selected ΔT/T maps at fixed Eprobe and time delay t = 150 fs are reported in Figure 2, for three tBLG stacks with different θ. As expected, the SLG regions exhibit a PB signal at all the investigated Eprobe, detected as a positive ΔT/T signal. The PB signal only shows a smooth monotone increase towards lower energy, following the momentumspreading of the hot-carriers FD distribution.
This behavior contrasts to that of tBLG, where a complex dependence on Eprobe emerges showing also a PA signal. Specifically, under low photo-excitation fluences, we identify two characteristic components in the transient response: a clear peak in the PB and the appearance of PA at lower energy. The PB peak corresponds to a factor 4-6 enhancement with respect to the SLG signal and it is centered in proximity to EvHs (see first column of maps of Figure 2). The PA is detected for Eprobe < EvHs. At lower probe energies, about 300 meV below EvHs, the ΔT/T signal of tBLG becomes barely distinguishable from that of SLG. While a similar behavior is observed in all the investigated samples, the energy position of PB and PA bands changes with θ. The larger is the twist angle, the higher is Eprobe at which similar tBLG/SLG contrast is observed. In stacks 2 and 3, bilayer regions with different twist angle are present, so that the signal intensity evolves differently according to their orientation.
For all samples, the precise value of the twist angle is determined by recording the hot-PLE spectrum. The experimental ΔT/T spectra of monolayer and bilayer regions at t = 150 fs for several tBLG samples are reported in Figure 3(a), together with the PLE spectra used to extract EvHs and θ (blue symbols).
As θ decreases, the overall TA spectrum experiences a red-shift. The PLE peak almost overlaps with the PB peak, with the larger mismatch appearing in the stack with the most intense PA band. The signal near the vHs, indeed, results from the interplay of PA and PB contributions whose relative intensity varies with the twist angle and excitation fluence. Remarkably, the PA emerges at low photo-excited carrier densities of the order of n = 1 × 10 13 cm -2 , as estimated from an incident fluence of 50 μJ cm -2 considering 4% absorption. At an incident fluence of 350 μJ cm -2 , corresponding to n = 7 × 10 13 cm -2 , the PA band disappears and a PB is observed over all the spectral range (see Fig. S-2 and S-3 of SM). The lack of previous evidences [47,48] of contribution from the ISB to the TA of tBLG could be traced back to the higher absorbed photon flux attained by pumping resonantly to the vHs. This dominance of the PB over the PA for high laser fluence was also observed for the intra-and inter-band transitions in suspended single layer graphene [31].
The dependence of the experimental TA spectrum on crystal orientation is traced back to the modifications of the phonon modes, electronic structure and of the JDOS, which we theoretically compute as a function of the twist angle and the hot-electron temperature Te (see details in Methods section) [57]. Figure 3 Figure 3(c) shows a comparison of both IB and ISB JDOS for carrier temperature of 1500 K and room temperature (300 K); the blue arrows indicate that the ISB JDOS intensity increases with the temperature while the IB JDOS intensity decreases. The intensity of the JDOS is related to the dynamical conductivity (or absorption). In this way, the difference of the JDOS at these two temperatures is proportional to the ΔT/T spectrum and qualitatively explains the opposite sign of the contributions from ISB and IB transitions to the TA signal.
The increase of the JDOS associated to ISB transitions due to photo-excitation corresponds to the PA signal observed at energies just below E vHs in the experimental data of Figure 3  The ΔT/T(t) dynamics of tBLG, on the contrary, strongly depends on the transition energy even though, away from the PA and the vHs, it is very similar to that of SLG. At the vHs, the PB not only has a much higher initial intensity than for SLG, but it decays on a longer time scale. The dynamics, after deconvolution with the IRF, is best fit by a biexponential decay with τ1 tBLG = 140 ± 20 fs and τ2 SLG = 2500 ± 100 fs. The PA band is best described by a three-exponential decay with time constants τ1 tBLG = 120 ± 30 fs, τ2 tBLG = 1960 ± 40 fs and τ3 tBLG = 70 ± 2 ps with the additional decay component responsible for the tail extended to 100 ps in Figure 4(b). Such long lifetimes can be explained by the peculiar temperature dependence of the JDOS that emerges at finite twist angles. We

Sample preparation and characterization
The twisted bilayer graphene single crystals were prepared by ambient pressure chemical vapor deposition (AP-CVD) synthesis on polycrystalline Cu foils. The nucleation conditions on Cu foils are controlled to allow both constituent layers to form as hexagonal single crystals, and each tBLG grain has a specific twist angle. The samples were transferred to a transparent amorphous quartz substrate for the optical images and spectroscopy measurements. More details about the sample preparation are presented in Ref. [62].
To select the samples for TA studies we first measured the sample twist angle by an analysis of the flake shapes on optical microscopy white light images. The optical images were obtained using an Olympus microscope (DP72) in reflection mode with a 100X objective lens. To obtain the twist angle we drew a hexagon for each layer and measured the angle between the edges of the two hexagons [52], as can be seen for the sample in Figure 1(d).
We also characterized the sample twist angle by Raman spectroscopy using an alpha300R WITec confocal Raman spectrometer at excitation energy of 1.96 eV (633 nm). The Raman spectra were measured in backscattering geometry on a triple monochromator spectrometer equipped with a liquid-N2 CCD. All measurements were performed at room temperature.
For the PLE measurements a PL image is acquired for each excitation laser energy, that is tuned to cover the range of the vHs transition for each sample. The PL intensity is averaged at the bilayer region and divided by the PL at the monolayer region to obtain each point in the PLE spectra. The PL signal is measured in back scattering, separated from the excitation laser by a dichroic mirror (Semrock FF665-Di02). A band pass filter (75 nm centered at 610 nm) is used in front of the detector to further block any scattered laser light (Chroma 610/75M). The cross-correlation between pump and probe pulses at the sample is measured to be ~ 0.3 ps, and the time delay interval between adjacent images is set to Δt = 133 fs. A single t scan of 20 ps takes about 1.5 min, and a complete spectral scan takes about 30 min.

Theoretical model
The electronic structures were obtained by folding the SLG calculation results. For SLG, we follow the procedure given in Ref. [60,57], where the electronic structure calculations are based on a fifth neighbors tight-binding approach, with one orthonormalized pz orbital per site, in which the hopping parameters are fitted to reproduce density functional theory (DFT) calculations with many-body corrections.
The electronic transitions between bands in tBLG are captured considering two SLG BZs twisted by a certain angle θ. In this way, we considered only four electronic bands for the tBLG, one valence and one conduction band for each layer. For low temperatures, only the electronic IB transitions are possible. The magenta arrow in Fig. 1(c) shows the IB that coincides with the transition between the vHs, that gives rise to a peak in the JDOS at an energy EvHs which is a function of θ. When the temperature increases, ISB transitions between valence, or between conduction bands, are also possible (red arrows in Fig. 1(c)).
The allowed IB and ISB transitions do not depend on the size of the minigap, which appears due to the hybridization of the eigenstates near to the band intersections. Our model cannot reproduce the interaction between layers and the atomic relaxation is not taken into account. However, Havener et al. [19] have shown that, although the interlayer interactions perturb the bands from each layer producing minigaps, which splits the valence and the conduction bands in tBLG, the difference in energy between the minigaps is generally too small to be seen directly in the PL and absorption spectra, contributing only to the peak broadening. Also, since the allowed electronic transitions do not depend on the size of the minigap and the PL peak is usually symmetric, which means that both transitions from minigaps to the conduction and the valence bands contribute equally, the energy position of the JDOS peak obtained by the folding approach is supposed to give a good estimation of the twisted angle. Therefore, the folding approach used in our simulations, which does not consider interaction between layers, is supposed to give good results [19,56].
In our model, it is possible to analyze the JDOS for IB and ISB transitions separately, as a function of the electronic temperature. We compute the JDOS as where the sum is performed on a uniform grid of Nk = 2400 × 2400 k-points over the two dimensional BZ and for the different pairs of bands i; j while γ =1 × 10 -2 eV is a small energy broadening parameter related to the inverse of the electron-hole lifetime. The function � , , � = [1 + exp ( , / )] −1 is the Fermi function for a specific band energy i of the state k and temperature Te (the chemical potential μ is close to the Fermi level and both are set to zero).
Using only IB transitions in the JDOS calculation we find the energy peaks corresponding to the EvHs for a specific twist angle θ. For the same values of EvHs obtained by the PLE peaks, we compute the normalized JDOS for IB and ISB transitions, as a function of the probe energy for Te = 1500 K, which is the estimated temperature for the hot electrons excited by the pump laser (see Fig. S-4 and S-5 in the SM). We see that the peaks corresponding to ISB transitions are always at lower energies than the IB transitions between the vHs, as shown in Fig. 3(c).
An estimation of the ΔT/T dynamics can be obtained assuming that electron-radiation matrix element is constant for all transitions. In this way, the dynamical conductivity σ (or absorption α) becomes proportional to the JDOS which we compute including both the IB and the ISB as function of the Te. Details of the Te time-evolution calculation are described in the SM.     The inset shows for the vHs energy in our experimental region of twist angles together with the best fit function (red solid line).
To calculate the uncertainty in the angle estimation we derive the expression: considering the as the PLE peak width. From the PLE results, the fitted Gaussian widths give ≤ 0.01 eV for all the stacks. For the bilayer stack with EvHs = 1.388 ± 0.007 eV, for example, we find that θ= 7.312° ± 0.04°.

Fluence dependence of transient absorption spectra
The intersubband (ISB) transitions energetically overlap with interband optical transitions. The hot-electron distribution causes simultaneously the bleaching of interband absorption due to Pauli blocking and the activation of ISB absorption with contributions to the transient transmissivity of opposite signs. The interplay of these contributions depends on the fluence and, at high density of photo-excited carriers, we observe that the bleaching is dominating over the photo-induced absorption over the entire probe photon energy range.
Accordingly, while the TA of SLG scales linearly with fluence, the spectra of tBLG change shape with fluence, as While the photobleaching signal increases with excitation power, due to the increase in initial electron temperature and to the superlinear dependence of ΔT/T on Te, the instantaneous PA signal decreases. In Fig. S-3 we report the pump power dependence of relaxation dynamics of the stack with θ = 7.31° at probe energy Eprobe = 1.240 eV corresponding to the ISB, and Eprobe= 0.956 eV in the low-energy PB band.

Calculating the electron temperature
In order to obtain the electron and phonon temperatures as a function of time we adopt a three temperature model [2]. This model considers that after a transition to an out-of-equilibrium state the electrons rapidly exchange energy with phonons via electron-phonon interaction achieving a quasi-equilibrium Fermi distribution. Considering the occupation number of electrons and phonons as Fermi-Dirac and Bose-Einstein distributions respectively, we can associate an effective temperature for the electrons, Te, and for phonons, Tν which are separated into strongly coupled (ν = sc) and weakly coupled (ν = wc) phonons. Thus, a system of coupled equations for the three temperatures (Te, Tsc, Twc) is used to simulate the dynamics of the electrons and phonons after the laser pulse excitation.
For a Gaussian pulse I(t) with fluence F and pulse duration τp, the temporal dependence of the effective temperatures is given by: where β is a dimensionless parameter which determines the fraction of the pulse energy absorbed by the electrons and τo = 2.5 ps is the optical phonon lifetime [3][4]. The respective specific heats for the electrons (ce), strongly coupled phonons (csc) and weakly coupled phonons (cwc) are obtained by the following integrals: Here, N( ) and Fν(ω) are the density of states for electrons, strongly coupled (ν = sc) and weakly coupled (ν = wc) phonons obtained, for both SLG and tBLG, from tight-binding [5,6] and force constant [7] methods parameterized by first principles calculations, and ( ) is the Fermi-Dirac (Bose-Einstein) distribution. The difference between weakly and strongly coupled phonons in graphene is addressed considering that the optical A'1 (TO) and E 2 g (LO) phonons in the vicinity of K and Γ points of the Brillouin zone are those with high electron-phonon coupling. In this way, in our model for SLG, only phonons with frequency close to (1350 ± 2) cm -1 and (1580 ± 2) cm -1 are considered in the calculation of the strongly coupled quantities, while the others modes are used for the weakly coupled ones.
For tBLG, the additional low frequency and high density phonon modes close to ZA2 (88 ± 2) cm -1 at the Γ point are also considered as strongly coupled phonons. This mode is not present in the SLG because its polarization represents the out-of-plane breathing between two adjacent layers due to the weak van der Waals interaction.

Relaxation dynamics of tBLG stacks
Similar relaxation dynamics is measured in all the probed tBLG stacks presented in Figure 2