High-Harmonic Generation with a twist: all-optical characterization of magic-angle twisted bilayer graphene

If we stack up two layers of graphene while changing their respective orientation by some twisting angle, we end up with a system that has striking differences when compared to single-layer graphene. For a very specific value of this twist angle, known as magic angle, twisted bilayer graphene displays a unique phase diagram that cannot be found in other systems. Recently, high harmonic generation spectroscopy has been successfully applied to elucidate the electronic properties of quantum materials. The purpose of the present work is to exploit the nonlinear optical response of magic-angle twisted bilayer graphene to unveil its electronic properties. We show that the band structure of magic-angle twisted bilayer graphene is imprinted onto its high-harmonic spectrum. Specifically, we observe a drastic decrease of harmonic signal as we approach the magic angle. Our results show that high harmonic generation can be used as a spectroscopy tool for measuring the twist angle and also the electronic properties of twisted bilayer graphene, paving the way for an all-optical characterization of moir\'e materials.


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 the 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][8][9], ferromagnetism [10] and also strange metal behavior [11,12].Such 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 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 alignement among others.This has for example allowed observing 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][16][17] at the so-called magic angle θ mag ≈ 1.1 • .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 bandstructure 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 of new spectroscopic tools as the large moiré unit cell and the two dimensionality of MATBG compli-Figure 1.Schematic figure of the moiré superlattice alongside with the laser directions.Note that the Γ − M direction connects AA to AA zones while Γ − K direction connects AA to AB/BA zones.We note that the moiré Brillouin zone is not shown at scale.cate 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 non-equilibrium 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 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,27], topological transitions [28,29], and light-driven band structure [30].
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 [31], higher-order ones are caused by the latter contributions [32].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 [33].
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 [34][35][36][37][38][39], its ultrafast dynamics have been hardly explored, with very few exceptions [40][41][42][43].However and to the knowledge of the authors, no previous work has dealt with highly nonlinear optical response in MATBG.Here we show that the emitted spectrum reveals the formation of the flat bands at the magic angle.

RESULTS
As it has been discussed earlier, the HHG of solids is extremely dependent on the electronic band structure of the system [44,45], up to the point that one can reconstruct the band structure from its high-harmonic spectrum [46].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 θ 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 Methods for more details on the model and the time propagation scheme).To test the effect of the flat bands on the emission spectra, we will scan a set of twist angles, near θ mag , ranging from 1.00 • to 1.35 • .In Fig. 2, we show the band structure for 3 different angles: the magic angle itself (θ = 1.12 • in our model), one angle above (θ = 1.22 • ) and one below it (θ = 1.05 • ).One can notice how the bands lose their flatness when we move away from θ mag .
Figure 3(a) shows the harmonic spectrum of the TBG for a laser in the Γ−M direction (see Fig. 1) as a function of the twist angle.A main feature emerge: a strong suppression of the intensity of the odd non-linear harmonics for twist angles near θ 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 • versus 1.12 • .
Even though the process here described is rather complex, the physical intuition behind it, is not.At a semiclassical level, the intraband current can be expressed as [32] where ε n (k(t)) is the energy of the n-th band and the crystal momentum is given by k(t) = k + eA(t).Therefore, the most relevant contribution to the emission spectra is the curvature of the bands.Hence, when the band starts to flatten, the harmonic intensity will be reduced.While this semiclassical explanation cannot account for all the subtle, yet interesting, physics of the system, it provides an understandable picture of this process.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.Figure 3(c,d) shows the same quantities as in 2(a,b) but for a laser in the Γ − K direction, see Fig. 1.Independently of the laser direction, the same physics appears; when the twist angle gets close to θ 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. Intuitively, one would expect that when upper bands get populated, the current will have relevant intraband contributions and hence, the decrease on the emission will be smoothed out, as the harmonic signal will not be exclusively originated from the flat bands.
To check this hypothesis, we performed a wide temperature scan for a set of angles near θ mag .Fig. 4 shows the yield for the 3rd 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 which signify the minimum of the 3rd harmonic, the minimum is located at the magic angle (1.12 • ).However, this depression near θ 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.
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 ω = 1 meV), spurious contributions from higher conductions bands can be avoided.Such small laser frequencies are crucial for our study and to detect experimentally the drop in the harmonic efficency 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 [47], making the experiments proposed experimentally feasible.

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 that 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.

METHODS
To model TBG we start from the continuum model [15,17,20] taking a graphene Fermi velocity of v F = 2.34a 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 [48].With these parameters the magic angle is found at θ ≈ 1.12 • .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 the model satisfies the symmetries: C 3 , C 2 T (the combined operation involving time reversal symmetry T and C 2 ), and a valley preserving mirror symmetry M .Here C 3 and C 2 , involve respectively 120 • and 180 • rotations with respect to an axis perpendicular to the TBG and M is a layer exchanging two fold rotation around an axis parallel to the TBG.We then adapt the Wannier function model of TBG with 8 orbitals per valley [49,50] 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 r, and we restrain ourselves to the diagonal approximation [45,51].
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 [52] that are well suited for the model describe above, following the formalism presented in [45].During this work, we have considered the interaction of the twisted bilayer graphene with a relatively weak (E 0 = 1 • 10 5 V/m) pulse in the microwave regime (ω = 1 meV) and a gaussian envelope with a full width at half-maximum of 23.5 ps.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.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 × 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 [32,45] was set to T 2 = 1000 fs ≈ Ω/4, where Ω is the laser period.However, we checked that the same results were obtained if the value of T 2 was varied.

Figure 2 .Figure 3 .
Figure 2. Band structure of TBG for 3 distinct twist angles alongside its Fermi-Dirac distribution.a corresponds to θ = 1.05 • , b to θ = 1.12 • and c to θ = 1.22• , showing the appearance of flat bands at the so-called magic angle θ = 1.12 • .Red color indicates the band structure associated with the K valley, while blue lines depicts its K counterpart (see Methods for a discussion on the model).

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