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

Moir\'e 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 dc conductivity anisotropy induced by a nematic order depends on the degree of lattice relaxation, the doping and the nature of the symmetry breaking.


I. 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 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 twisted bilayer graphene (TBG) with θ ∼ 1.1 o 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 four fold degenerate, stemming from the spin and the graphene * calderon@icmm.csic.es † leni.bascones@icmm.csic.es 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 two-fold 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 sensi- tivity 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 inter-band 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 Equation 20.
Optical conductivity in non-correlated states for different values of the electronic filling. We focus on σ λλ (ω), 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, σ λλ (ω) 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. 1 (c) 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 non-interacting band structures of the two model  [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. 1(d) and (e)) is characterized by the absence of a Drude peak at zero frequency, but with interband 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 σ λλ (ω) 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. 2(b)), only produce a small feature in σ λλ (ω) at ω ∼ 0.6 meV.
Beyond the frequencies dominated by the γ 1 transi-tions, there is a gap up to the frequencies of the γ 2 and γ 4 transitions, in Fig. 1(c). 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. 1(c)) 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. In a correlated order which breaks the C 2 T symmetry, 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 σ XiXi = σ YiYi is also fulfilled. A nematic state lowers the rotational symmetry and the optical conductivity spectra along the affected directions become different.

Results in symmetry breaking states
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]. Fig. 2 (c) shows the optical conductivity for the nematic order η proposed to explain the observed changes in the DOS [9]. As illustrated in Fig. 2 (c) the nematicity is revealed in the spectrum via the inequalities σ X1X1 = σ X2X2 = σ X3X3 and σ XiXi = σ 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. 2 (d) is not easily related to the enhanced peak separation in the DOS shown in Fig. 2 (d). 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 non-interacting state at charge neutrality is a semi-metal, a Drude peak may arise for finite values of η C3 evidencing that a semi-metal to metal Lifshitz transition takes place. This Lifshitz transition is illustrated in Fig. 3 (c) to (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. 3 (c)). For small η C3 the Dirac points move in kspace from K and K' [13] but still lie at µ CNP (Fig. 3 (d)). 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. 3 (e)). The value of the order parameter at which the pockets emerge produces changes in the DOS, red curve in Fig. 2 (f), compatible with the band widening observed in STM experiments [9]. Additional Fermi pockets appear for larger values of η C3 , also leading to new band cross-  ings between the lower and upper flat bands ( Fig. 3 (f)). These crossings produce pairs of Dirac points with opposite vorticity which drift away in energy and momentum as the nematicity increases. The Lifshiftz transition is not specific to the nematic state η or to the PR model. As shown in Fig. 2 (f), 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. 3 (h). The presence of several maxima in the lower flat band of the FR model (see Fig. 1 (b) and Fig. 3 (h)) 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 Lifshiftz 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 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. 1(c)): 3 p orbitals (p zT , p +T and p −T ) at the triangular lattice determined by the AA regions, 4 p orbitals (p A +H , p A −H , p B +H and p B −H ) at the hexagonal lattice of the 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 fully relaxed (FR) model with θ FR ∼ 0.9 o 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 o partially relaxed (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 to 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. 1 (f)) 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  Equation 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 val-ley and spin. A small broadening 0.04 meV is used in the calculation of σ λλ (ω) except in Fig. 2(e) and (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.

DATA AVAILABILITY
All relevant data are available from the authors upon reasonable request. E.B. conceived the project. Both M.J.C. and E.B. designed and performed the calculations, analysed and discussed the results, and wrote the manuscript.

Tight-binding model
We describe the non-interacting band-structure with a 10-band tight-binding model per valley based on effective moiré orbitals. The ten effective moiré orbitals are located at the three different lattices formed by the symmetry points of the TBG (see Fig. 1(f)): 3 p orbitals (p zT , p +T and p −T ) at the triangular lattice determined by the AA regions, 4 p orbitals (p A +H , p A −H , p B +H and p B −H ) at the hexagonal lattice formed by the AB and BA regions with a p + and a p − at each of the two inequivalent 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.
This tight-binding was introduced in Ref. [22] where parameters were chosen to fit the bands obtained from the continuum theory for TBG with a twist angle of θ = 1.05 o [42]. These parameters correspond to what we call PR (partially relaxed) model in the main text. The same tight-binding was used in Ref. [23] where it was fitted to the band-structure resulting from an ab-initio k · p continuum model which account for the lattice relaxation for a twist angle of θ = 0.9 o . We use these parameters for our FR (fully relaxed) model. The parameters for the PR and FR models are listed in Table S1 and the corresponding band structure is plotted in Fig. S1. The in-plane relaxation strongly reduces the approximate particle-hole symmetry of the continuum model and results in a more complicated flat bands dispersion. In both cases, the considered angles are below the corresponding magic angles (1.08 o for PR and 1 o for FR).
Here we write down the 10 × 10 hamiltonian in the symmetric formulation required to calculate the optical conductivity. It is useful to define φ lm = e −ik·(la1+ma2) (S1) with a 1 = aŷ and a 2 = a √ 3 2x − 1 2ŷ the lattice vectors, l, m fractional numbers and a the moiré lattice constant. φ lm condenses the information about the hopping directions. The tight binding hamiltonian is written For convenience, in order to write h 10×10 (k) we differentiate the six orbitals at the triangular and kagome lattices from the 4 orbitals at the hexagonal lattice. In the following we do not write explicitly the dependence on k. Here H pz describes the p zT intraorbital term H ± gives the intra and interorbital hopping between the p ±T orbitals and C p±pz the interorbital hopping between p zT and p ±T : H κ includes the hopping between the orbitals in the kagome lattice and C κp± between s κ (1,2,3) and p ±T terms are: The onsite energies are written in terms of the hopping parameters: µ pz = −6t pz + δ pz , µ p± = 3t p± + δ p± , and µ κ = −4(t κ + t κ ) + δ κ . The hopping between these six orbitals and the four orbitals in the hexagonal lattice is given by H 6×4 Here h (A) and h (B) with ζ = e i2π/6 and Ω = ζ 2 . Finally, the terms within the hexagonal lattice are (S16)

Symmetry breaking states
To study the optical response to symmetry breaking correlated states we introduce phenomenologically three different terms which break any symmetry of the Hamiltonian of a given valley: S1. Values of the parameters used for the two models considered in the main text. We use the notation in Ref. [22].
Parameter PR model [22] (meV) FR model [23] H α breaks the C 2 T and the M 2y symmetries while H η and H β break the C 3 symmetry [9]. These terms can be interpreted as the interaction part of a mean-field hamiltonian and are introduced to explore some generic features of the optical conductivity and the modification of the band structure under a symmetry breaking phase transition. Nevertheless H α and H η were obtained in a self-consistent calculation in Ref. [9] using a model based on the PR tight-binding with a Hubbard interaction restricted to the three p-orbitals in the triangular lattice. The state η, namely with finite η C3 , was proposed to explain the enhanced peak separation in the density of states observed when the chemical potential lies in the flat bands, as compared to empty or full bands, observed in this work [9]. Fig. S2 shows the modification of the bands in the ordered states η and β along several directions in k-space, with Fermi pockets at the CNP, as in Fig. 3. These pockets may emerge even for a quite moderate modification of the band structure: an order parameter β C3 three times smaller than the one used in Fig. S2(b) moves the Dirac points below the CNP in Fig. 3(h). When α C2T , η C3 and β C3 vanish we refer to the system as non-correlated.
With these terms we aim to represent the symmetry breaking effects associated to an electronic transition. Nevertheless some of the features present in the optical conductivity are generic for a given symmetry breaking independently of whether the symmetry breaking is due to the electronic interactions or a lattice effect, even if the latter should be modelled with different terms in the hamiltonian. The alignment with h-BN has two effects on twisted bilayer graphene: (i) it induces a weak moiré potential due to the different lattice constant of h-BN and graphene, and (ii) it breaks the sublattice symmetry (and consequently the C 2 T symmetry in each valley) due to the different potential exerted by the boron and nitrogen atoms on the two carbon sublattices [44]. The second effect opens a gap at the Dirac points, as also found in the α order. Uniaxial strain in the device would break the C 3 symmetry making one of the directions X i different to the other two and producing an anisotropic response in the optical conductivity while shifting the Dirac points from K and K as in η and β states. If the strain is different in the two layers (heterostrain), the Dirac points would not remain degenerate [45,46].

Optical Conductivity
Here we reproduce the expressions for the conductivity used in the calculations and derived in Ref. [21] for a Hamiltonian H mf bilinear in fermionic operators using Kubo formula and introducing the coupling between the electrons and the electromagnetic field via a Peierls substitution: ×θ( n (k))θ(− n (k))δ(ω − n (k) + n (k)), n (k) and n (k) are the band energies, θ( n (k)) the Heaviside step function and the Drude weight D λ can be written D λ = −π kn t λλ nn (k)θ(− n (k)) − 2π V kn =n j λ n n (k) 2 n (k) − n (k) θ( n (k))θ(− n (k)) , with t λλ nn (k) = µν ∂ 2 µν (k) ∂k 2 λ a * µn (k)a νn (k) , with µν the elements of the hamiltonian written in the orbital basis H mf = Σ k,µ,ν µν (k)c † µ,k c ν,k = Σ k,n n (k)d † n,k d n,k and a µn (k) the rotation matrix between the orbital and the band basis c † kµσ = n a * µn (k)d † knσ . In the present case c † µ,k and c ν,k are the creation and annihilation operators for the ten orbitals p zT , p +T , p −T , s κ1 , s κ2 , s κ3 , p A +H , p A −H , p B +H and p B −H . H mf includes both the tight-binding h 10×10 (k) and the mean-field terms H η or H β . H η and H β can be written as hopping terms and couple to the electromagnetic field, contributing to the vertices t λλ nn (k) and j λ n n (k), while H α is a local onsite term and its contribution to the vertices vanishes. As shown in Fig. S3, in the absence of correlations and with the chemical potential within the flat bands, the low energy interband transitions γ 1 become progressively suppressed and a Drude peak emerges as the system is doped from the CNP. The transitions γ 6 , whose threshold correspond to transitions around Γ, are easily detected only when the bands are completely filled. At the same doping, the transitions γ 4 become completely suppressed.