1. INTRODUCTION

The discovery of magic-angle twisted bilayer graphene (MATBG) [1,2] has sparked a huge increase of attraction into moiré quantum materials [3] in recent years. One of the key reasons for that increase is that MATBG shows a unique plethora of exotic phenomena [4], which includes correlated insulators [1,5], unconventional superconductivity [2,5,6], interacting topological phases [7–9], ferromagnetism [10], and also strange metal behavior [11,12]. Such a phase diagram has some resemblance with the phenomenology present in cuprates [13] and iron pnictides [14] superconductors. For instance, both systems share the presence of superconducting domes surrounded by correlated insulating states as doping is varied. Moiré materials, see Fig. 1, have opened completely new avenues in research thanks to the unprecedented in situ tunability of their electronic properties through electronic gates, strain, magnetic field, pressure, or substrate alignment, among others. This has, for example, allowed the observation of many different correlated insulating states in a single MATBG sample, including integer and fractional Chern insulators as well as trivial insulators with or without charge density modulations [9].

At the heart of the correlated phases of MATBG it is the existence of two flat bands [15–17] at the so-called magic angle ${\theta _{{\rm mag}}} \approx {1.1^ \circ}$. As the bandwidth is reduced, the role played by the interactions increases [18,19], and rich and complex physics arises. The nature of these flat bands is far from trivial, as they are formed by the hybridization of the Dirac cones of each layer of graphene [1,20]. A proper characterization of the low energy band structure of MATBG is key to unveil the origin of the correlated states. On the other hand, studying the electronic properties of these flat bands requires new spectroscopic tools as the large moiré unit cell and the two dimensionality of MATBG complicate the applicability of many experimental techniques of common use in other strongly correlated materials.

In parallel, the study of the interaction between ultrashort laser pulses and condensed matter systems has also experienced a huge growth in relevance [21]. Such short pulses can drive the system’s electrons into a nonequilibrium excited state within times below their own cycle. As a consequence of this excitation, the electrons emit coherent radiation in multiples (up to the hundreds) of the laser frequency, providing us with a tool to unveil their dynamics. Even though this high-harmonic generation (HHG) had its origins in atomic and molecular physics [22], it has been recently, and successfully, applied to condensed matter systems as a spectroscopy tool [23,24]. For instance, it has allowed for the observation of ultrafast electron-hole dynamics [25], insulator-to-metal transitions [26–29], topological transitions [30–32], and light-driven band structure [33]. In particular, in 2D materials, experimental efforts are being carried toward developing HHG as a novel probe of moiré patterns, with the latest demonstration of HHG from stacked layers [34].

Why HHG can be effectively used as a spectroscopy tool can be understood thanks to its underlying physical mechanisms. The HHG spectrum in solids comes from two distinct contributions, namely intraband and interband currents. Intraband harmonics are generated due to the accelerated motion of the electron (hole) in the conduction (valence) band after a multiphoton/tunneling excitation of the electron from the valence band. Afterwards, the recombination of such pair gives rise to the interband harmonics. While intraband contributions are generally associated with low order harmonics [35], higher-order ones are caused by the latter contributions [36]. In the presence of a strong low frequency laser field, the charge carriers are drifted through the whole Brillouin zone, exploring sections, and hence physics, that are not commonly seen by other optical tests, such as second harmonic generation [37].

In the present work, we aim to bridge the gap between these two distinct fields by studying HHG in magic-angle twisted bilayer graphene. Remarkably, in spite of extensive work in twisted bilayer graphene [38–43], its ultrafast dynamics have been hardly explored, with very few exceptions [44–47]. While the previous references show unique and interesting physics, we want emphasize that to the best of our knowledge, there has been no previous work that deals with the ultrafast dynamics of magic-angle twisted bilayer graphene. Furthermore, one of the main objectives of the present work is to show that we can identify the magic angle in an all-optical way, which also has not been achieved before.

2. RESULTS

As has been discussed earlier, the HHG of solids is extremely dependent on the electronic band structure of the system [48,49], up to the point that one can reconstruct the band structure from its high-harmonic spectrum [50]. Therefore, the harmonic spectra of twisted bilayer graphene must inevitably reflect the existence (or absence) of the flat bands.

The purpose of this work is to show that one can characterize the magic angle in an all-optical way, offering an alternative to other characterization techniques, such as scanning tunneling microscopy [16]. We show that when the twist angle of the TBG gets near ${\theta _{{\rm mag}}}$, the intensity of the high harmonics (from the third upwards) is reduced by several orders of magnitude. Moreover, it is shown that the effect is robust against important variations of the laser field and to temperature changes, providing a more robust probe of magic angle over the existing ones. To check this hypothesis, we performed numerical simulations. (See Section 4 for more details on the model and the time propagation scheme.) In order to optically isolate the flat bands, i.e., avoid one-/two-photon transitions between the lower valence bands and the upper conduction bands, we selected a laser field in the terahertz regime ($\omega = 1\,\,{\rm meV}$) and a relatively weak peak field of 100 kV/m. To test the effect of the flat bands on the emission spectra, we scanned a set of twist angles, near ${\theta _{{\rm mag}}}$, ranging from ${1.00^ \circ}$ to ${1.35^ \circ}$. In Fig. 2, we show the band structure for three different angles: the magic angle itself ($\theta ={ 1.12^ \circ}$ in our model), one angle above ($\theta ={ 1.22^ \circ}$), and one below it ($\theta ={ 1.05^ \circ}$). One can notice how the bands lose their flatness when we move away from ${\theta _{{\rm mag}}}$.

Figure 3(a) shows the harmonic spectrum of the TBG for a laser in the $\Gamma {-} M$ direction (see Fig. 1) as a function of the twist angle. A main feature emerges: a strong suppression of the intensity of the odd nonlinear harmonics for twist angles near ${\theta _{{\rm mag}}}$. In order to obtain a clearer grasp of the magnitude of this decrease, Fig. 3(b) presents the harmonic yield from the third up to the seventh harmonic. In spite of the small angle differences, the intensities of each of the three harmonics differ by several orders of magnitude. This disparity would allow for all-optical characterization of the magic-angle twisted bilayer graphene, even for twist angles close to each other; see, for instance, the difference between ${1.06^ \circ}$ and ${1.12^ \circ}$.

Even though the process here described is rather complex, the physical intuition behind it is not. The high-harmonic spectrum of TBG is primarily dominated by the intraband contribution. This is because the flat bands are extremely close to each other and, thus, the energy of the recombination processes that make up the interband term will be equally small due to the energy conservation. We have checked numerically (see Supplement 1) that the interband contribution is generally two orders of magnitude smaller than the intraband one. At a semiclassical level, the intraband current can be expressed as [36]

 

Fig. 1. (a) Schematic figure of the moiré superlattice alongside (b) laser directions. Note that the $\Gamma {-} M$ direction connects $AA$ to $AA$ zones, while the $\Gamma {-} K$ direction connects $AA$ to $\textit{AB}{/}\textit{BA}$ zones.

Download Full Size | PDF

 

Fig. 2. Band structure of TBG for three distinct twist angles alongside its Fermi–Dirac distribution. (a) corresponds to $\theta={ 1.05^ \circ}$, (b) to $\theta={ 1.12^ \circ}$, and (c) to $\theta={ 1.22^ \circ}$, showing the appearance of flat bands at the so-called magic angle $\theta={ 1.12^ \circ}$. Red color indicates the band structure associated with the graphene ${K_g}$ valley, while blue lines depict its counterpart at the graphene ${K^\prime _g}$ valley. Along $\Gamma {-\!\!-}K{-\!\!-}M{-\!\!-}\Gamma$ and $({-\!\!-}M){-\!\!-}({-\!\!-}K){-\!\!-}\Gamma$ directions both bands overlap and only the blue ones are seen. (See Section 4 for a discussion of the model.)

Download Full Size | PDF

Consequently, this phenomenon must be robust under important variations of the laser because it is only based on the existence of the flat bands. This intuition is confirmed by our numerical simulations. Figures 3(c) and 3(d) show the same quantities as in Figs. 2(a) and 2(b) but for a laser in the $\Gamma {-} K$ direction; see Fig. 1. Independently of the laser direction, the same physics appears; when the twist angle gets close to ${\theta _{{\rm mag}}}$, a depression of several orders of magnitude appears in the emission spectrum.

As seen above, the suppression of the harmonic intensity close to the magic angle crucially depends on the narrowing of the flat bands, but what happens when other bands are involved on the dynamics? For instance, when the upper conduction bands have some population due to an increase in the temperature (see Fig. 2), one would expect the current to have relevant interband and intraband contributions from these upper band and the decrease on the emission would be smoothed out, as the harmonic signal would not be exclusively originated from the flat bands.

 

Fig. 3. (a) and (c) High-harmonic spectrum for a set of twist angles. (b) and (d) Yield for the third (H3), fifth (H5), and seventh (H7) harmonics. (a) and (c) results are obtained for a laser in the $\Gamma {-} M$ direction, while (c) and (d) are obtained for the same laser but in the $\Gamma {-} K$ direction. All calculations have been done for a temperature of 80 K. One can see the clear depression of the harmonic signal, for all subfigures, near the magic angle.

Download Full Size | PDF

To check this hypothesis, we performed a wide temperature scan for a set of angles near ${\theta _{{\rm mag}}}$. Figure 4 shows the yield for the third harmonic in terms of the twist angle and the temperature. For temperatures below 100 K, the depression on the intensity is heavily located around the magic angle. As seen from the location of black circles that signify the minimum of the third harmonic, the minimum is located at the magic angle (${1.12^ \circ}$). However, this depression near ${\theta _{{\rm mag}}}$ disappears when the temperature is enough to populate the upper bands, and the minimum moves away from the magic angle. This effect is still visible for temperatures that are well above the ones at which the correlated insulators and superconductivity have been observed in MATBG.

 

Fig. 4. Yield of the third harmonic in terms of temperature. Black circles depict the minimum of the yield for each twist angle.

Download Full Size | PDF

For low frequency laser pulses, the most important contributions to the current, and therefore to the harmonic signal, come from bands near the Fermi energy, in the case of TBG, the flat bands. In particular, by choosing laser frequencies that are smaller than the typical bandwidth of the flat bands in TBG (here $\omega = 1\,{\rm meV}$), spurious contributions from higher conduction bands can be avoided. Such small laser frequencies are crucial for our study and to detect experimentally the drop in the harmonic efficiency due to the flattening of the bands. Laser frequencies and intensities similar to the ones considered here have been used in high-harmonic generation measurements on a graphene monolayer [51], making the proposed measurements experimentally feasible. In the latter work, the experiment was conducted using a superconducting accelerator-based source of THz radiation. Nevertheless, the same parameter set can be achieved with a table-top, laser-driven THz source [52]. This alternative approach would greatly simplify the experimental setup.

3. DISCUSSION

Our results have shown that the flatness of the bands of MATBG has a strong impact on its high-harmonic spectrum. Their narrowing close to the magic angle produces a decrease of several orders of magnitude in the emission intensity. Such discrepancy in the optical response can be used in a fully optical study of MATBG. This effect is crucially based only on the existence of the flat bands, and it holds for a wide range of laser parameters, even showing no anisotropy in the laser direction. In addition, we have confirmed that the suppression is maintained within a relatively wide range of temperatures up until 100 K, far away from the ultralow temperature regime used in most MATBG experiments. Nevertheless, we hope our work can not only be useful as a spectroscopy tool but also shed light onto the mechanisms behind electronic ultrafast dynamics in twisted materials.

4. METHODS

To model TBG we start from the continuum model [15,17,20] taking a graphene Fermi velocity of ${v_F} = 2.34a\, {\rm eV}$, where $a$ is the lattice constant of graphene, and a ratio between the interlayer tunneling at AA and AB regions of ${w_0}/{w_1} = 0.78$ [53]. With these parameters the magic angle is found at $\theta \approx {1.12^ \circ}$. In the continuum model the two valleys of graphene are assumed to be uncoupled. The spectrum of TBG equals the sum of the contributions of each valley. In each of the valleys, ${K_g}$ and ${K^\prime _g}$, the model satisfies the symmetries: ${{\cal C}_3}$, ${{\cal C}_2}{\cal T}$ (the combined operation involving time reversal symmetry ${\cal T}$ and ${{\cal C}_2}$), and a valley preserving mirror symmetry $M$. Here ${{\cal C}_3}$ and ${{\cal C}_2}$, involve, respectively, ${120^ \circ}$ and ${180^ \circ}$ rotations with respect to an axis perpendicular to the TBG and $M$ is a layer exchanging twofold rotation around an axis parallel to the TBG. We then adapt the Wannier function model of TBG with eight orbitals per valley [54–56] to approximate the bandstructure model obtained with the continuum model. The resulting tight binding model includes hopping up to a radius of 10 moiré lattice constant. This Wannier model satisfies the symmetries of the continuum model in each valley. The TBG is assumed to be undoped for all the calculations by setting the chemical potential for each temperature. Since we are interested in the optical response of the bilayer, we need to assume some form for the position operator ${\textbf r}$, and we restrain ourselves to the diagonal approximation [49,57].

Inversion and time reversal symmetry are satisfied only when the two valleys are included. The spectrum of the even harmonics is finite for each valley, but it vanishes when adding the contribution of both valleys.

The time-dependent calculations were performed by solving the semiconductor Bloch equations (SBE), using maximally localized Wannier functions [58] that are well suited for the model describe above, following the formalism presented in [49]. During this work, we have considered the interaction of the twisted bilayer graphene with a relatively weak $({E_0} = 1 \cdot {10^5}\,{\rm V/m})$ pulse in the microwave regime $(\omega = 1\,{\rm meV})$ and a Gaussian envelope with a full width at half-maximum of 23.5 ps. Even though the intensity may seem not strong enough, we have checked that it is enough to be in the nonpertubative regime for the third and fifth harmonic (see Supplement 1). This optical regime was chosen so to properly isolate the flat bands, i.e., to avoid transitions between them and the upper conduction bands. Otherwise, such processes would give relevant interband contributions to the spectrum. Nevertheless, the physical behavior is robust under variations of both the amplitude and the frequency of the field (see Supplement 1). We note that, even though the field amplitude may seem weak, for the case of MATBG it will be relatively strong due to the reduced bandwidth of the system. The reduced density matrix was propagated in the Wannier gauge using a $100 \times 100$ Monkhorst–Pack grid and a time step of 2.41 fs. We checked that the dynamics were properly converged for all numerical parameters. The dephasing term [36,49] was set to ${T_2} = 1000\,\,{\rm fs} \approx \Omega /4$, where $\Omega$ is the laser period. However, we have checked that the effects observed here remain unaffected as we vary ${T_2}$ (see Supplement 1).