Correlated states in magic angle twisted bilayer graphene under the optical conductivity scrutiny

Moiré systems displaying flat bands have emerged as novel platforms to study correlated electron phenomena. Insulating and superconducting states appear upon doping magic angle twisted bilayer graphene (TBG), and there is evidence of correlation induced effects at the charge neutrality point (CNP) which could originate from spontaneous symmetry breaking. Our theoretical calculations show how optical conductivity measurements can distinguish different symmetry breaking states, and reveal the nature of the correlated states. In the specific case of nematic order, which breaks the discrete rotational symmetry of the lattice, we find that the Dirac cones are displaced, not only in momentum space but also in energy, inducing finite Drude weight at the CNP. We also show that the sign of the Drude weight anisotropy induced by a nematic order depends on the degree of lattice relaxation, the doping and the nature of the symmetry breaking.


INTRODUCTION
Instabilities in correlated materials arise when the interaction energy overcomes the kinetic energy gain. Hence, materials whose bands have a very small bandwidth (flat bands) are prone to show correlation effects. The band structure of some moiré systems can be tuned to display very narrow bands. A moiré pattern appears when there is a small mismatch between two overlapping grids, as it happens in twisted bilayer graphene (TBG) with one layer rotated by a relative angle θ with respect to the other. Very weakly dispersing bands are found in the so-called magic angle TBG with θ~1.1°around charge neutrality 1 . A plethora of insulating and superconducting states appear as the system is doped [2][3][4][5][6][7] and band widening has been observed in STM experiments when the chemical potential μ lies within the flat bands, in particular, in undoped systems (namely, at the charge neutrality point (CNP)) [8][9][10][11] . Some observations are sample dependent and are affected by whether the TBG is aligned to the underlying Boron Nitride (h-BN) as in-plane twofold rotation symmetry may then be broken 12 . Correlations have to be included to explain both the insulating states which appear at integer fillings and the doping dependence of the band widening [13][14][15][16] . The nature of these states is highly debated. An important issue, yet also unknown, is whether the correlated states affect only the flat bands or also involve higher energy bands, split by a few tens of meV gaps from the flat bands.
In the absence of symmetry breaking, the upper and lower flat bands touch at Dirac points and are fourfold degenerate, stemming from the spin, and the graphene valley degrees of freedom. For each valley the bands show M 2y , C 3 , and C 2 T symmetry 17 . Here M 2y is a mirror symmetry which flips the y coordinate producing a layer exchanging twofold rotation in 3D, C 3 , and C 2 rotations are referred to an axis perpendicular to the TBG and T stands for time reversal. Signatures of nematicity (C 3 symmetry breaking 18 ) have been found in measurements of the superconducting upper critical field in doped samples and observed in STM experiments when the chemical potential lies in the flat bands 9,10,19 . The latter happens at the same dopings at which STM measurements detect the band widening, namely an increased separation of the density of states (DOS) peaks [8][9][10][11] . These effects disappear if the flat bands are completely filled or empty and therefore they cannot be a simple consequence of strain. On the other hand, activated behavior at the CNP appears in some transport experiments 5,6 , reflecting a gap opening at the Dirac points. Whether this latter effect has an electronic origin or is due to coupling to the substrate has not been clarified. Nevertheless, no alignment with the h-BN has been detected in the experiments.
Here we study the effect of symmetry breaking states on the optical response of two tight binding models proposed for TBG with different degrees of relaxation. We show that the optical conductivity corresponding to C 2 T and C 3 symmetry breaking states display specific signatures that help identify these correlated states. In particular, optical conductivity is sensitive to the gapping of the Dirac points and to the breaking of discrete rotational symmetries. Unexpectedly, we find that a nematic order can induce Fermi pockets, and hence a finite Drude weight at the CNP. The dependence of the optical conductivity on doping gives information about whether the rotational symmetry breaking arises from the lattice or from the electronic degree of freedom and can clarify if the higher energy bands participate in the correlated state. Our calculations show a strong sensitivity of the nematic response to the still not well understood lattice relaxation, as well as to the nature of the nematic order. This sensitivity can be traced back to the flat bands structure and impacts dramatically the sign of the Drude weight and the dc conductivity anisotropies of doped TBG.

RESULTS
Optical conductivity experiments allow to probe the excitation spectrum of materials 20 . Specifically, they probe the q = 0 interband transitions between occupied and empty electronic states, giving accurate measurements of the direct excitation gap in insulators. If a material undergoes a phase transition, its band structure is modified and the inter-band transitions are sensitive to these changes, including a possible gap opening at low energies. In particular, the changes in the band structure responsible for the modification of the DOS observed in STM when the flat bands are partially filled, with respect to the one in fully empty or filled flat bands, can be addressed by optical conductivity measurements in the far-infrared regime. In the following, we present results for the optical conductivity calculated in linear response as given by Kubo formula, following the derivation for multiorbital systems in ref. 21 , see "Methods" and Supplementary Eq. 20.
Optical conductivity in non-correlated states Figure 1a, b shows the optical conductivity of two model systems for TBG, in the absence of correlations, for different values of the electronic filling. We focus on σ 0 λλ ðωÞ, the real part of the longitudinal electrical conductivity σ λλ (ω) which relates the electric current J λ generated in linear response in direction λ and the homogeneous electric field E λ applied in the same direction J α = σ λλ (ω)E λ . In metals, σ 0 λλ ðωÞ shows a peak at zero frequency, the Drude peak, related to the d.c. conductivity. At higher energies there is contribution from the inter-band transitions with energy ℏω and zero momentum q = 0, weighted by matrix elements that depend on the electronic momentum k and on the orbital composition of the bands. Following Pauli exclusion principle, the allowed transitions depend on doping. The arrows in Fig. 1c indicate the optical transitions at CNP (black), with empty flat bands (red), and with full flat bands (blue). The optical transitions are labelled γ m with m ranging from 1 to 6.
The noninteracting band structures of the two model systems, a fully relaxed (FR) TBG with θ FR~0 . 9°and a partially relaxed (PR) TBG with θ PR~1 .05°, are displayed, for a single valley, in Fig. 1d, e and in the Supplementary Information 22,23 . In both models the chemical potential of the undoped system μ CNP touches the Dirac points, located at the K and K′ points of the Brillouin zone. Besides the different bandwidths, a clear contrast between both models is the reduced particle-hole symmetry of the FR model.
The optical spectrum of the undoped compounds (black curves in Fig. 1d, e) is characterized by the absence of a Drude peak at zero frequency, but with inter-band transitions contributing from zero energy as expected from the presence of Dirac points at μ CNP [24][25][26] . At small frequencies the optical conductivity is controlled by the transitions γ 1 between the two flat bands and is finite up to a frequency of the order of the flat bands bandwidth. The inter-band transitions are weighted by matrix elements and the maxima in σ 0 λλ ðωÞ are not simply related to those in the DOS. For example, in the PR model the conductivity peaks around ω~2 meV while the DOS sharp peaks at ω~0.3 meV, associated with the presence of van Hove singularities at M and M′ (see for instance the black line in Fig. 2b), only produce a small feature in σ 0 λλ ðωÞ at ω~0.6 meV. Beyond the frequencies dominated by the γ 1 transitions, there is a gap up to the frequencies of the γ 2 and γ 4 transitions, in Fig. 1c. For the small twist angles considered here the threshold for these transitions is determined by a region in momentum space close to Γ. Interestingly, at this threshold there is a sharp peak for the PR model while a small shoulder is found in the FR model. This distinction arises, again, from differences in the matrix elements. A smaller peak at higher frequencies signals the onset of γ 3 transitions.
In doped TBG new transitions between the flat and higher energy bands (γ 5 and γ 6 in Fig. 1c) are allowed. As μ is moved away from the Dirac points, a gap opens in the optical spectrum corresponding to transitions γ 1 . Nevertheless, if the chemical potential still lies in the flat bands, the system is metallic and there is a Drude peak at zero frequency (see Supplementary Information). For fully empty or occupied bands, the optical conductivity is zero up to a frequency given by the γ 5 or γ 6 transitions, respectively.
Results in symmetry breaking states Figure 2 shows how the optical conductivity at CNP responds to different symmetry breaking states.
In a correlated order which breaks the C 2 T symmetry, named α here, the flat bands at Γ are barely affected but the Dirac points at K and K′ are gapped 9,13,27 , see Fig. 3a, b, producing activated behavior in transport measurements, as the one observed in some experiments. This gap shows up in the optical spectrum of the undoped system, see Fig. 2a, inducing a reorganization of the spectral weight which resembles the one in the DOS shown in Fig. 2b.
In contrast, a nematic state does not open a gap in the optical conductivity of undoped TBG 9,13,27 . Optical conductivity is a directional probe which reflects the discrete rotational symmetry of a solid. Hence nematicity can be detected by studying the spectrum along different directions. In the absence of correlations or external symmetry breaking the optical conductivity spectra along the directions related by the symmetry rotations of the crystalline lattice are equal. In TBG, the lattice has C 6 symmetry when both valleys are included (C 3 in each valley) and the optical conductivities along any two axes rotated by ±π/3 are equal. Moreover, for any system with C 3 symmetry, the equality σ 0 Xi Xi ¼ σ 0 YiYi is also fulfilled. A nematic state lowers the rotational symmetry and the optical conductivity spectra along the affected directions become different.
In STM experiments, when the flat bands are partially filled the X 1 direction becomes different to directions X 2 and X 3 and the DOS is modified [9][10][11] . Figure 2c shows the optical conductivity for the nematic order η proposed to explain the observed changes in the DOS 9 . As illustrated in Fig. 2c the nematicity is revealed in the spectrum via the inequalities σ 0 YiYi . The differences show up at frequencies corresponding to the γ 1 , γ 2 , γ 4 , and even γ 3 transitions, revealing that not only the flat bands but also the higher energy bands are affected by this nematic order. In this state the reorganization of the optical spectrum in Fig. 2d is not easily related to the enhanced peak separation in the DOS shown in Fig. 2d. The low frequency spectrum depends non-monotonically on the amplitude of the order parameter η C3 and, as the nematic order enters as an effective hopping term, the spectral weight at low frequencies is not conserved.
Unexpectedly, while the noninteracting state at charge neutrality is a semimetal, a Drude peak may arise for finite values of η C3 evidencing that a semimetal to metal Lifshitz transition takes place. This Lifshitz transition is illustrated in Fig. 3c-f, which show the evolution with η C3 of the band structure close to μ CNP . In the non-correlated state with η C3 ¼ 0 the Dirac points are placed at K and K′ as required by C 3 symmetry (Fig. 3c). For small η C3 the Dirac points move in k-space from K and K′ 13 but still lie at μ CNP (Fig. 3d). Further modifications of the flat bands with increasing η C3 shift the Dirac points away from μ CNP and generate small hole and electron Fermi pockets (Fig. 3e). The value of the order parameter at which the pockets emerge produces changes in the DOS, red curve in Fig. 2f, compatible with the band widening observed in STM experiments 9 . Additional Fermi pockets appear for larger values of η C3 , also leading to new band crossings between the lower and upper flat bands (Fig. 3f). These crossings produce pairs of Dirac points with opposite vorticity which drift away in energy and momentum as the nematicity increases.
The Lifshitz transition is not specific to the nematic state η or to the PR model. As shown in Fig. 2f, a finite Drude weight D X1 also appears in the optical spectrum of the FR model in a different nematic state β and originates in the Fermi pockets in Fig. 3h. The presence of several maxima in the lower flat band of the FR model (see Figs. 1b and 3h) makes it more prone to display Fermi pockets at the CNP even for small order parameters. See Supplementary Information for more details.
The anisotropy of the Drude weight, more specifically the sign of the A i B j -anisotropy ratio ðD Ai =D Bj Þ À 1 between the Drude weights D Ai , D Bj of two otherwise equivalent directions A i and B j , has been used to characterise the nematic state of other compounds, such as iron superconductors [28][29][30][31][32][33] . In Fig. 4 we plot the X 1 X 2 and X 1 Y 1 anisotropy ratios for the PR and FR models in different nematic states as a function of the electronic filling and the nematic order parameter. Black lines denote the values with zero anisotropy and delimitate the regions with different signs. The anisotropy vanishes if the order parameter η C3 or β C3 is zero. The vertical line at the CNP marks a region with no Drude weight. The end of this line indicates the value of the order parameter at which the Lifshitz transition takes place. The sign of the anisotropy, originating in our approach in the morphology and topology of the Fermi pockets, is similar for X 1 X 2 and X 1 Y 1 for a given TBG model and nematic order but it is notably distinct for different TBG models with the same order or for the same TBG model with different nematic orders: for example, the sign of the X 1 X 2 anisotropy of the FR model in state η is given by the sign of the order parameter η C3 while the PR model anisotropy in the same state or the FR model anisotropy in state β display several sign changes.

DISCUSSION
Our results show that the optical conductivity can distinguish different symmetry states and reveal their correlated nature. In particular, both the reduction of the rotational symmetry, detected in STM experiments, and the presence of a gap at the CNP, observed in transport measurements, can be directly seen in the optical conductivity spectrum. We have shown that the lack of symmetry may be observable not only at low frequencies but also at higher ones. We have restricted our discussion of the transitions to higher energy bands around Γ near the threshold for γ 2 or γ 4 transitions as they are expected to be more easily identifiable. At these frequencies the nematic order could be detected. Nevertheless, at higher frequencies the transitions around K or M involving flat and higher energy bands would also be modified by both C 3 and C 2 T symmetry breakings. If these transitions are identified 34 , they could be used to study the correlated state. Some of the signatures produced by strain or alignment with the substrate would be similar to the ones discussed here, even though the lattice symmetry breaking effects should enter in the model differently. However, the lattice effects on the spectrum should not vary significantly with doping or temperature, allowing the distinction between a lattice or electronic origin. In the absence of correlations, the chemical potential shifts rigidly and the peaks in the spectrum associated to these transitions do not shift in frequency. However if, as observed in STM experiments, the signatures of the correlated state, namely the breaking of the C 3 symmetry or other modifications of the band structure, disappear when the flat bands are empty or full, a frequency shift is expected in the peaks of the spectrum which originate in transitions between the bands affected by the correlations. Doping away from the charge neutrality point, the transitions γ 1 between the flat bands become progressively forbidden but the information can still be obtained from the evolution with doping of transitions γ 2 -γ 4 . These transitions may be easier to detect experimentally than the γ 1 , even in undoped systems, as they appear at significantly higher frequencies.
Detection of a shift with doping of the frequency of the transitions γ 3 , which do not involve the flat bands, would reveal the involvement of the higher energy bands in the correlated state, at present under debate. In general, the study of the different transitions could allow a momentum selective analysis. For example, the transitions γ 2 -γ 4 around Γ dominate the spectrum close to the threshold while those around K or M appear at higher frequencies. On the other hand, the frequency of the γ 1 transitions, whose contribution is very sensitive to the electronic filling, is largest at Γ.
The nematic order can be also detected by looking at the Drude weight anisotropy in optical or transport experiments. However, unlike for other materials, like iron superconductors [28][29][30][31][32][33] , the anisotropies' sign displays a complicated dependence on both the doping and the order parameter. Such a sensitivity to details suggests that a strong sample variability could be found in the sign of the anisotropy of TBG.
Our finding of Lifshitz transitions induced by nematicity impacts not only the optical conductivity but many other experimental measurements, including transport, quantum oscillations, or STM. The details of the Lifshitz transition would vary in different samples as the degree of relaxation of the lattice significantly affects the shape of the flat bands and hence how they are affected by the nematic order, making the relaxed lattices more sensitive to it. In general, Fermi pockets emerge already for magnitudes of the order parameter whose effect in the flat band structure is small or compatible with the reorganization of the DOS observed experimentally.
The gap observed at charge neutrality in transport experiments could be due to the breaking of C 2 T symmetry 35 , discussed here, but also to other electronic orders such as intervalley coherent states or valley/spin polarized orders [36][37][38][39] . However, note that different theoretical predictions do not agree on whether the polarized states at CNP are insulating or metallic [37][38][39] . Within the quasiparticle description adopted here, the contribution of different spins and valleys to the spectrum add. As the spin and valley are not coupled by the hopping, if the gap at CNP is due to valley or spin polarized order, the γ 1 transitions should be Fig. 3 Three-dimensional band structure in correlated and non-correlated states. Flat bands of the partially relaxed (PR) model in the a non-correlated state, and b the C 2 T symmetry breaking state. The non-correlated state has Dirac points protected by C 2 T symmetry at the K and K′ points as required by the C 3 symmetry. When the C 2 T symmetry is broken, a gap is open at the Dirac points while, in the particular state under study, the region around Γ is weakly affected. c-f Zoom of the flat band structure for the PR model and different values of η C3 . In c, η C3 ¼ 0 and the Dirac points are placed at K and K′. If η C3 is small, as in d, the Dirac points move in momentum but stay at μ CNP . When its magnitude increases, as in e, f, the bands at small energies are modified displacing the Dirac points away from the chemical potential of the undoped system and creating small Fermi pockets at CNP. Hence, nematicity produces a Lifshitz transition turning a semimetal into a metal. g, h Same as c, e for the FR model in the nematic state β. In all the panels a plane is plotted at μ CNP for reference.
suppressed for the spin/valley flavor whose bands become empty (filled). Concomitantly, the γ 5 (γ 6 ) peaks will emerge. Therefore, the temperature dependent spectral weight reorganization could help distinguish whether the gap originates in spin/valley polarization versus the breaking of C 2 T symmetry or intervalley coherence.
Finally, very recent experiments have shown that for certain fillings the population of the flat bands occurs through a sequence of phase transitions at which a single spin/valley flavor takes all the carriers and becomes completely empty/filled 40,41 . Such phase transitions should also be detectable via optical conductivity as they would lead to a redistribution of the total spectral weight as a function of doping and temperature.

Tight binding models
We use a 10-band tight binding model for each valley and spin previously proposed to satisfy all the symmetries and the fragile topology of the continuum model of TBG 22 , and later confirmed by a projection from an ab-initio k ⋅ p perturbation model which includes the lattice relaxation 23 . Different choices for the parameters reproduce the band structure corresponding to different angles and degrees of lattice relaxation. The model is based on effective moiré orbitals, not on carbon atomic orbitals. The ten effective moiré orbitals are located at the three different lattices formed by the symmetry points of the TBG (see Fig. 1c): 3p orbitals (p zT , p +T , and p −T ) at the triangular lattice determined by the AA regions, 4p orbitals (p A þH , p A ÀH , p B þH , and p B ÀH ) at the hexagonal lattice of the AB and BA regions with a p + and a p − at each of the two A and B lattice sites, and three s orbitals, s K1 , s K2 , and s K3 , at the kagome lattice formed by the SP points, with a single orbital at each of the three sites of the kagome lattice. The flat bands, namely the two bands close to zero energy in Fig. 1, are mostly composed by the p +T and p −T orbitals having, nevertheless, p zT or kagome s κ1 , s κ2 , and s κ3 character at the Γ point.
The FR model with θ FR~0 . 9°is derived from an ab-initio k ⋅ p approximation which includes the effect of in-plane and out of plane relaxation 23 . The tight binding for the θ PR~1 .05°PR TBG is obtained from the continuum model and only vertical corrugation is introduced through a 15% suppression of the interlayer hopping between carbon atoms in the same sublattice 22,42 .

Symmetry breaking states
The three orders α, η, and β are introduced phenomenologically at the mean field level, but the orders η and α were obtained self-consistently from a microscopic model in ref. 9 . To modify the energy of the flat bands the three orders involve the p +T and p −T orbitals. The C 2 T symmetry breaking α order breaks the degeneracy between p +T and p −T with an onsite diagonal term. This makes the two atomic sublattices in each layer inequivalent. The nematic states η and β arise from bond orders: η breaks the C 3 symmetry of the interorbital hopping between p +T and p −T while in β the hopping amplitude between the p +T and p −T orbitals and s κ1 differs from the corresponding hopping to s κ2 and s κ3 . In other words, in the non-nematic states the hopping amplitudes in the directions joining the nearest neighbor hexagon centers (Y 1 , Y 2 , and Y 3 directions, see Fig. 1f) are equal. In these nematic states, the hopping amplitude in the Y 1 direction is different from the other two, making this bond inequivalent to the other two 9 . See Supplementary Information for further details.

Optical conductivity
In the optical conductivity calculation, we assume a band picture and neglect possible excitonic effects. In a strongly correlated metallic state with both quasiparticle and Hubbard bands, this kind of description will approximate the behavior of the quasiparticle spectrum while Hubbard bands would give further contributions 20 . In the derivation of the explicit expression, Supplementary Eq. 20, for multiorbital systems 21 the coupling to the electromagnetic field is introduced via a Peierls substitution. Within the used approximation the matrix elements for the optical conductivity are written as derivatives of the orbital-dependent tight binding terms. The 10-band model used here has several atoms per unit cell and applicability of the expressions in ref. 21 requires writing the Hamiltonian such that the specific distances between the atoms are accounted for, as the phase acquired in the Peierls substitution is sensitive to these distances and not to unit cell labels 43 , see Supplementary Information. Provided the valley and spin degeneracies are not broken, in the used approximation, the optical conductivities of both valleys and spins are equal. We therefore focus on the optical conductivity of a single valley and spin. A small broadening 0.04 meV is used in the calculation of σ λλ (ω) except in Fig. 2e, f where it is 0.06 meV. The Drude weight is calculated directly from the hamiltonian and the current matrix elements without any broadening, see Supplementary Information. Fig. 4 Anisotropy of the Drude weight. Colour plots for the a X 1 X 2 and b X 1 Y 1 Drude anisotropy for the PR model in the nematic state η as a function of doping and order parameter η C3 . The A i B j anisotropy is defined as the ratio ðD Ai =D Bj Þ À 1 between the Drude weights D Ai and D Bj along two directions A i and B j . Doping is defined per spin and valley with respect to the CNP: doping is zero at the CNP and it is 1 and −1, respectively, for completely filled and empty bands. c, d Same as a but for the FR model respectively in η and β nematic states.