Nematic insulator at charge neutrality in twisted bilayer graphene

We investigate twisted bilayer graphene near charge neutrality using a generalized Bistritzer-MacDonald continuum model, accounting for corrugation effects. The Fermi velocity vanishes for particular twist angles properly reproducing the physics of the celebrated magic angles. Using group representation theory, we identify all contact interaction potentials compatible with the symmetries of the model. This enables us to identify two classes of quartic interactions leading to either the opening of a gap or to nematic ordering. We then implement a renormalization group analysis to study the competition between these interactions for a twist angle approaching the first magic value. This combined group theory-renormalization study reveals that the proximity to the first magic angle favors the occurrence of a layer-polarized, gapped state with a spatial modulation of interlayer correlations, which we call nematic insulator.


INTRODUCTION
Within band theory, a reasonable estimate for the relative strength of the quasiparticles kinetic energy is the ratio between the bandwidth of the conducting bands, and some interaction energy. In normal metals, where the density of states at the Fermi energy is nonzero, the excitations of the Fermi sea largely screen the Coulomb interactions, which renormalizes their strength downwards in a dramatic way. But when the bands near the Fermi energy disperse very little, even small interactions can lead to significant, qualitative consequences.
The prominent example of such systems is provided by Landau levels and the associated fractional quantum Hall effects [1,2]. More recently, a different class of materials with vanishing bandwidth was uncovered in twisted bilayer graphene (TBG). When two sheets of graphene are rotated with respect to one another by a small angle of about 1.1 • , a large moiré pattern forms with several thousands of atoms per unit cell. Remarkably, for some twist angles-the so-called magic angles-the Fermi velocities of the Dirac cones originating from each layer vanish exactly [3]. Other phenomena accompany this effect such as a minimal bandwidth and a sizeable band gap between the conducting and excited bands, under some conditions [3][4][5][6].
TBG gained considerable momentum after the experimental discovery of correlated insulators at various fillings, with neighboring regions of possibly-unconventional superconductivity [7][8][9][10]. In addition, scanning tunneling microscopy and transport data point toward "ferromagnetism" and an "anomalous Hall effect" in TBG, but also in thicker van der Waals heterostructures like trilayer graphene [11][12][13][14]. The tunability of TBG through a rich phase diagram by electronic gating also sparked numerous works in new directions. While the origin of the superconductivity remains unclear [15][16][17][18][19][20][21][22], evidence is mounting toward the strongest insulator emerging at charge neutrality-where band theory alone would predict a semimetal-with a charge gap of around 0.86 meV in the most angle-homogeneous devices [23]. Other scanning measurements find a three-fold rotation symmetry breaking near the first magic angle at charge neutrality [24]. The aim of this article is to unveil the nature of this rotation symmetry breaking insulator at charge neutrality close to the first magic angle and to provide methodology to analyze its occurrence.
The study of interacting phases in systems with vanishing bandwidth is notoriously difficult. Our strategy to tackle this challenge in TBG is based on the combination of an algebraic identification of interactions preserving the symmetries of the low energy description of TBG and a renormalization group approach to select the generically favored interaction as the twist angle approaches its first magic value. In the present case, an additional obstacle lies in the absence of a simple description of the moir pattern in TBG, thus inhibiting the use of standard field theories.
Our starting point for the non-interacting continuum model of two twisted layers of graphene accounts for several channels of interlayer hoppings. These interlayer hoppings renormalize the Fermi velocity, leading to its vanishing at magic twist angles. We develop a diagrammatic technique to compute the Green's function and thereby the magic angles to arbitrary order in the interlayer hopping strength α and non-perturbatively in the imbalance between different hopping channels β. This non-interacting model is then complemented with interactions. By formal group theory considerations, we identify all symmetry-allowed contact or short-ranged interactions. To determine the most favorable one, we develop a renormalization group (RG) technique. The vanishing of the bandwidth at the magic angle leads to a singular behavior of the RG: indeed, any interacting potential, while usually treated in perturbations, now corresponds to a dominant energy scale. Moreover the first magic angle is not determined by a specific value of a parameter of the free field theory, but through a systematic resummation of interlayer hopping terms. To overcome these difficulties, we study the scaling behavior of all interacting potentials as the twist angle is varied. When approaching the first magic angle, we monitor the relevance in the RG sense of all potentials, thereby identifying the dominant interacting instability. We find that as 1. In a continuum model, the relative twist of the top (green) and bottom (red) layers by an angle θ leads to one mini-Brillouin zone (mBZ) for each valley of monolayer graphene. The two Dirac cones of the same valley, Kt and K b , set the sides the mBZ of size 2K sin(θ/2), where K = |K t,b | is the Dirac momentum of monolayer graphene. Electron hoppings between the two layers involve a small momentum transfer qj, j = 1, 2, 3 between Kt and each of the three nearest K b nodes of the mBZ.
the twist angle approaches the first magic value, a state with both a gap opening and a periodic modulation of interlayer correlations is favored. We call this phase a nematic insulator.

II. FREE ELECTRON MODEL
Following the seminal work of Ref. [25], we treat TBG as a periodic moiré superlattice characterized by a twist angle θ. The top and bottom Dirac cones of the same valley, denoted K t and K b , delineate the mini-Brillouin zone (mBZ) of the superlattice (Fig. 1). Focusing on the low energy and long wavelength description of TGB, we restrict ourselves to small momentum transfers that are diagonal in valley, and thus occur within a single mBZ [26]. The characteristic kinetic energy scale of the model, set by the typical difference of kinetic energy of electrons in different layers, is E c = 2v 0 K sin(θ/2), where v 0 and K are respectively the Fermi velocity and the Dirac momentum of monolayer graphene. In addition to the kinetic energy in each layer, the single-particle Hamiltonian involves two different interlayer hopping amplitudes. First, the amplitude w 1 of interlayer hopping that is off-diagonal in graphene sublattice is typically of order w 1 ≈ 110 meV [25,27]. Its strength relative to the kinetic energy is measured by the dimensionless parameter α = w 1 /E c . Second, the amplitude w 2 = βw 1 of interlayer hopping that is diagonal in graphene sublattice is measured by the relative strength β ∈ [0, 1] in comparison to off-diagonal hopping. This relative strength is difficult to determine precisely in experiments, being affected by corrugation effects, with typical values evaluated as β ≈ 0.82 [28,29]. Here we keep β as a free parameter. Note that our model thus interpolates between the Bistritzer-MacDonald continuum (BMC) model for FIG. 2. Diagrammatic expansion of the electron self-energy to order 6 in the interlayer hopping amplitude α = w1/Ec relative to the kinetic energy Ec. The wavy line represents a pair of opposite hopping processes, summed over all channels with a transfer of momentum ±qj, j = 1, 2, 3. Diagram (a) is of order α 2 , diagram (b) of order α 4 , and diagrams (c)-(e) are of order α 6 . The expansion is non perturbative in the relative strength β between hoppings off-diagonal and diagonal in sublattices. β = 1 [3] and a chirally symmetric continuum (CSC) model for β = 0 [6]. Following Ref. [3], we use a rotated basis where the Dirac cones K t,b of the two layers have the same (k x , k y ) coordinates in the mBZ, and measure all energies in units of E c (see Appendix A for details). The effective Hamiltonian then reads where ∂ = (∂ x , ∂ y ) and the hopping matrices T + j are Here we introduced two sets of Pauli matrices, σ and τ , which describe respectively the sublattice and layer sectors, with σ z = ±1 = A/B and τ z = ±1 = top/bottom. The low-energy physics of this model is nontrivial even without interactions. Indeed, the interlayer couplings prohibit diagonalizing H 0 . This forbids the use of a simple effective theory valid for all twisting angles θ in Fermi velocity v(α, β) renormalized by interlayer hoppings at order α 6 , as a function of the relative strength α of the hopping amplitude with respect to the kinetic energy. As the twist angle increases, so does α, and the renormalized velocity vanishes at the first magic angle encoded in the first magic value α0(β) where β sets the asymmetry between diagonal and off-diagonal in sublattice hoppings. Inset: α0(β) depends weakly on corrugation effects, i.e. on the value of β.
the vicinity of the magic values. As a result, we resort to a free electron model in which interlayer hopping effects are accounted for by a self-energy which is calculated in a perturbative expansion in α. Denoting G 0 and Σ the translationally-invariant components of the propagator corrected by interlayer hoppings and the selfenergy respectively, we have where ∂ τ represents the partial derivative with respect to imaginary time, N ψ is a wavefunction normalization and v(α, β) the Fermi velocity renormalized by the hopping processes (see Appendix B). An expansion to order 6 in α but exact in β, diagrammatically represented in Fig. 2, leads to We call α 0 (β) the lowest value of α for which this Fermi velocity vanishes, which sets the first magic angle value to be approximately 1.1 • . As shown in Fig. 3, this first magic value depends weakly on the parameter β, and thus on corrugation, and ranges from α 0 (1) = 0.598 for the BMC model to α 0 (0) = 0.585 for the CSC model. These constitute our first results.

III. SYMMETRY-ALLOWED INTERACTIONS
We now identify all short-ranged interaction potentials allowed by the symmetries of the model. In order to do so we turn to a field theoretic formalism and consider the Euclidean action, S = S 0 + S int , written as a sum of the free electron term and an interaction term S int which includes generic local quartic couplings between the fermionic fields ψ † and ψ.
Using group theoretic methods, detailed in Appendix C, we identify all couplings allowed by the symmetries of the low energy model (1) [30]. This amounts to identiyfing scalar invariants built as direct products of irreducible representations of the corresponding symmetry group. We find that the allowed couplings are (i) 8 channels originating from one-dimensional (1d) corepresentations; (ii) 4 channels originating from 2d corepresentations: where the densities ρ (i) (r) = ψ † R (i) (r)ψ and currents J (j) (r) = ψ † M (j) (r)ψ involve coupling matrices R (i) (r) and M (j) (r). Following our choice of coordinates for can either preserve ( ) or break the combination of inversion and time reversal symmetries IT , the mirror symmetry C2, and the particle-hole antisymmetry P , while preserving the three-fold rotational symmetry C3. The ± exponents label the eigenvalue of the IT symmetry.
the fields in Eq. (1), the coupling matrices in the rotated basis, which enter Eq. (5), are obtained through where A(r) is the transformation matrix of the field (see Appendix A). The coupling matricesR (i) andM (i) are provided in Tab. I. The couplings g i and λ j are the amplitudes associated with the corresponding coupling potentials. The moiré pattern is invariant under four discrete "symmetries": (i) the 2π/3 rotation C 3 = e 2iπ/3σz τ 0 around the z axis orthogonal to the bilayer together with (ii) the π rotation C 2 = σ x τ x around the x axis of Fig. 1 generate the point group D 3 , (iii) the composition IT = σ x τ 0 K of inversion and time reversal is an antiunitary symmetry, where K denotes complex conjugation, and (iv) the unitary particle-hole antisymmetry [31] P = σ x τ z , which satisfies {P, H 0 } = 0 [32]. The group generated by D 3 and P comprises all unitary operations that leave the Hamiltonian invariant up to a sign. We refer to this ensemble as the unitary groupD 3 of the model. It can be decomposed into the semi-direct prod-uctD 3 = {e, P,ē,P } D 3 , where e = σ 0 τ 0 is the identity operation,ē = (P C 2 ) 2 = −e andP =ēP . The dichromatic magnetic group M generated byD 3 and IT can be written as the direct product M =D 3 × {e, IT }. Using the Schur-Frobenius criterion [33][34][35][36][37][38], we determine the corepresentations (corep.) of M from the irreducible representations ofD 3 , which can be constructed from that of D 3 by induction and basic properties of linear representation theory (see Appendix C).
We combine the resulting coupling matrices of Tab. I into three sets. The eight interactions originating from 1d coreps. correspond to the density-density couplings diagonal in sublattice while preserving C 3 . Out of these 4 (a) 4. Schematic dispersion relation of (a) a gapped layerpolarized correlated phase and (b) a density modulated phase. While the gapped phase is characterized by the amplitude of the gap ∆z opening at the Kt and K b Dirac points, the second phase is characterized by a shift of these Dirac points of amplitude G z . This shift leads to a modulation of the relative amplitude of wavefunctions between the two layers, revealed in the density |ψt + ψ b /2| 2 probed e.g. by a STM tip located on the top layer. The behavior of this density is shown (c) without any shift and (d) for a shift G z = 0.3ey. The spatial C3-breaking of this phase is clearly manifested by the appearance of stripes for this density, perpendicular to G z . eight couplings, (i) the four interactions associated with 1d coreps. which preserves IT are those symmetric on the A/B sublattices, of the form (ψ † σ 0 τ µ ψ)(ψ † σ 0 τ µ ψ) for µ = 0, x, y, z [39]. They are distinguished by their breaking of C 2 or P symmetries. (ii) The four interactions associated with 1d coreps. which break IT are those which are antisymmetric in the A/B sublattices, with couplings of the form (ψ † σ z τ µ ψ)(ψ † σ z τ µ ψ). Similarly to set (i), they can break C 2 or P . Finally, the (iii) four interactions originating from 2d coreps. are current-current couplings between the layers, off-diagonal in sublattices, of the form (ψ † στ µ ψ) · (ψ † στ µ ψ). They break all symmetries of the free model, and in particular the three-fold rotational symmetry C 3 .

IV. NATURE OF THE CORRELATED PHASES
Let us first discuss the nature of the phases induced by these couplings. The interactions of type (i), sym-metric in sublattices, neither open a gap at the Dirac point nor induce a density modulation. On the other hand, interactions of type (ii) generate phases with a gap ∆ µ ∝ g µ ψ † σ z τ µ ψ , reminiscent of the gap opening in Boron Nitride, see Fig. 4(a). These various gapped phases are distinguished by their layer correlations. Current-current interactions of type (iii) lead to radically different phases, in which the Dirac cones of the two layers are shifted with respect to each other by a momentum 2G µ ∝ λ µ ψ † στ µ ψ , as shown in Fig. 4(b). They generate gapless phases with C 3 -breaking density modulations. These spatial modulations are detected when probing the electronic density from one side of the bilayer, which amounts to coupling the local probe asymmetrically to the top and bottom wavefunctions, thus scanning some interlayer density of the form |ψ t + rψ b | 2 , where 0 < r < 1 is the asymmetry parameter. In Fig. 4 we compare the corresponding density for an asymmetry r = 1/2 in the absence of any instability in Fig. 4(c) with that in the presence of a C 3 -breaking instability in Fig. 4(d). Stripe-like modulations of the interlayercorrelated density are readily observed in this last case.
To gain further insight into the behavior of these phases close to the first magic angle, we now study their mean-field behavior. As we will show later using a renormalization group analysis, only four out of the twelve couplings are sensitive to the proximity of the magic angle. They correspond to the interaction potentials diagonal in layers, and originate from the a − 1 , a − 2 coreps. of set (ii) with respective amplitude g z , g 0 , and the E − 4 , E + 2 coreps. of set (iii) with amplitudes λ 0 , λ z . The corresponding order parameters satisfy the self-consistency equations [40] The correlators in Eq. (6) are the translationallyinvariant parts of statistical averages computed over the Bloch Hamiltonian density The corrections by interlayer hoppings of the correlators in Eq. (6) are obtained within a perturbation expansion in α. Incorporating the hoppings H α leads to an enhancement of the order parameters by factors N , as is the case for the renormalization of the Fermi velocity; they are calculated diagrammatically to sixth order in α in Appendix D.
The resulting dependence of each separate order parameter on the proximity to the magic angle and for various strengths of the couplings is depicted in Fig. 5. The insulating phases, characterized by a gap ∆ 0 or ∆ z , develop at a critical coupling which decreases as the parameter α approaches its magic value α 0 . Such gapped phases generically occur in some range of twist angles around the magic value, in agreement with the experimental findings of Ref. [41]. At the mean-field level and for a fixed α, we find that a ∆ 0 insulator occurs for weaker couplings g 0 than the couplings g z required for the appearance of the ∆ z insulator. We will see that this hierarchy is modified when fluctuations are accounted for, demonstrating the necessity to develop the RG approach.
The situation for the C 3 symmetry breaking phases with periodic modulation is different: while the phase which is antisymmetric in layers, characterized by a G z momentum, is also favored by the vanishing bandwidth close to α 0 , the finite critical strength for the analogous phase symmetric in layer, associated with G 0 , does not depend on α.
Having identified all interacting instabilities of TBG and established their strong enhancement close to the magic angle, we now study the competition between them by resorting to a renormalization group technique.

V. RENORMALIZATION GROUP PICTURE
As we have seen, the vanishing of the kinetic energy scale set by the renormalized velocity v entails that the four non-trivial interactions are relevant close to the first magic angle. In order to identify the leading instability, we study via a RG approach the competition between them as α approaches the first magic value. Our starting point is the field theory described by the action S = S 0 + S int where the quadratic action S 0 given in Eq. (4) involves both the kinetic energy and the interlayer hoppings. We sum over the latter to capture the vanishing energy scale. Hence, we expand the correlation functions of the interacting theory not only in the coupling constants but also in the amplitude of interlayer hoppings.
Motivated by the experimental observations of an insulating behavior at charge neutrality, we carry out a "particle-hole" Hubbard-Stratonovich transformation suitable to describe gapped phases (as opposed to superconductors). We introduce the scalar bosonic fields φ i , i = 1, ..., 8 for each 1d channel and the vector bosonic fields ϕ j , j = 1, ..., 4 for each 2d channel, so that the interaction part of the action (5) becomes: We expand around the lower critical dimension, setting the space-time dimension to d = 2 + , and renormalize the theory using the minimal subtraction scheme. We introduce the renormalized couplings g ∈ {g i , λ j } related to the bare ones byg = µ − N 2 Here µ is the momentum scale at which we renormalize the theory and N ψ is the wavefunction normalization factor generated by interlayer hoppings which was introduced previously. The renormalization constants Z g and Z φ for φ ∈ {φ i , ϕ j } absorb the poles of the three-point vertex and the bosonic self-energy, respectively. We expand these functions to first order in the interaction couplings using the Green's function corrected by interlayer hoppings G 0 as the fermionic propagator. This expansion is represented by the diagrams shown in Fig. 6, where interlayer hoppings are treated perturbatively to second order in α. We obtain the RG flow equations for the coupling constants as a function of the parameters α and β. Crucially, we find that only four interactions out of twelve have non-zero divergent corrections. The eight other couplings have a trivial flow, either because the correction has no pole in -these correspond to the four channels with σ 0 sublattice structure; or because the correction vanishes at low energy -these correspond to the four channels that are off-diagonal in layer space. We thus restrict our study to the four-dimensional subspace corresponding to the instabilities of Tab. II, which are all associated with a phase transition toward a correlated phase. These relevant couplings are all diagonal in layer.
We now briefly discuss the essential features of this four-dimensional flow. The gaussian fixed point (FP) at the origin is always stable in d = 3. Besides, we identify four critical points, one for each non-trivial coupling, listed in Tab. II. They control phase transitions toward the four correlated phases discussed in the meanfield analysis. As α approaches the magic value α 0 , all four critical FPs collapse towards the gaussian FP. Meanwhile, the (Dirac) semimetallic region, which corresponds to the basin of attraction of the gaussian FP, shrinks and disappears completely. As a result, these four couplings are always relevant close enough to the magic angle, regardless of the value of the bare interaction strength. This scenario provides a natural way of identifying the dominant instabilities near the magic angle: they correspond to the couplings whose critical FPs collapse the fastest towards the origin. Hence, as follows from Tab. II, we discard the couplings g 0 and λ 0 (the amplitudes of interactions which are symmetric in layers) and focus on the competition between the couplings which are antisymmetric in layers, of amplitude g z associated with the layer-polarized gapped phase and λ z associated with the C 3 -breaking densitymodulated phase. We note that, while the gapped phase is reminiscent of the dynamical mass generation in the Gross-Neveu model [42,43], the C 3 -breaking densitymodulated phase is specific to TBG. The competition between the two most relevant instabilities is dictated by the following coupled RG flows (for derivation see Ap- As mentioned above, all FPs collapse to the origin at the magic angle. Thus to explore the competition between the phases, we plot the renormalization flow for the couplings rescaled by the vanishing velocity. The effect of the proximity to the magic angle on this competition is shown in Fig. 7, where we compare the flow close to the first magic angle (b) with that for the case when interlayer hopping is suppressed (a) [44]. This comparison shows that the proximity to the magic angle favors the occurrence of density modulations. The large scale behavior is dominated by the fastest diverging coupling, whether g z or λ z . Within our perturbative RG analysis, a crossover line separates the corresponding regions, whose parametric equation reads Around the crossover line, both order parameters coexist over a large range of length scales, corresponding to the appearance of a gapped, periodically modulated state, asymmetric in layers and breaking the C 3 and IT symmetries. We call it a nematic insulator by analogy with phases discussed in Ref. [45]. This nematic insulating behavior is characterized by a runaway RG flow of both λ z , g z . It appears for a wide range of coupling parameters, as a consequence of the proximity to the magic angle.

VI. DISCUSSION AND OUTLOOK
Applying group theory supplemented by a renormalization group approach, we found that a gapped nematic state with C 3 breaking modulation of density is favored at charge neutrality in TBG when the twist angle approaches its first magic value. A gap was observed at charge neutrality in TBG both in scanning tunneling microscopy and spectroscopy studies [14,24,46,47] as well as in four-terminal transport measurements [23]. Threefold symmetry breaking and nematic ordering were also reported [22,24,47]. Both of these experimental observations strongly support the occurence of a nematic insulating state at charge neutrality in TBG, as we obtain within our RG scenario. We also find that such a  [28,29], (a) for a weak interlayer hopping amplitude, α = 0.1α0; and (b) close to the first magic angle, α α0. As α approaches α0, the blue semimetal region shrinks to the origin. To study the competition between the couplings we rescale them by the vanishing velocity v. The red critical FPs control the transitions toward the gapped (red region) or density-modulated (green region) phases. The black source FP gives rise to a crossover region (mixed state). It migrates away from the vertical axis as we increase α, thus expanding the density-modulated region.
state persists even when the strength of interactions is weakened by screening as was experimentally observed in Ref. [48]. Let us stress that our RG approach identifies a gapped nematic behavior in the perturbative scaling regime, but does not rule out that other types of correlations develop at larger length scales, including those of intervalley-coherent and generalized ferromagnetic insulating states recently discussed in Refs. [49][50][51][52]. It is worth noting that the energies of these different ground states seem to be very close to each other, suggesting a strong sensitivity to experimental conditions: indeed, h-BN encapsulation, which induces chirality breaking, favors the layer-polarized insulators such as the nematic insulator discussed in this paper as opposed to intervalley-coherent or generalized ferromagnetic insulating states [52,53]. Finally, let us note that the occurence of an analogous nematic insulating state close to quantum spin Hall phase transitions raises the questions of its relation with the topological nature of the underlying semimetal [32,54,55]. The Hamiltonian describing the low-energy physics near the two twisted Dirac cones at K t,b originating from a single valley of graphene can be written as [56] The momentum q 1 = K t −K b gives the relative displacement of the Dirac momentum K of each layer due to the twist, while v 0 is the Fermi velocity of graphene. Notice that there are three equivalent K points in monolayer graphene, each leading to one copy of the Hamiltonian H with a relative displacement q j , j = 1, 2, 3, where the momenta q 2 and q 3 are obtained through a rotation of q 1 by an angle of 2π/3 and 4π/3 respectively. The interlayer hopping matrixT (r) reads where t AA and t AB are the hopping amplitudes in the AA and AB/BA regions, respectively. Hamiltonian (A1) is simplified by rotating the basis [3] ψ(r, τ ) = A 1 (r)ψ(r, τ ), A j (r) = e −i(qj ·r/2)τz , (A3) which brings the Dirac cones to the same momentum: where T (r) = 3 j=1 e −iqj ·r T + j . Applying the same change of basis to quartic terms in the action, e.g. corresponding to a density-density interaction of the form we arrive at with the rotated interaction matrix Though bothR and R(r) describe contact interactions, whileR is space-independent, R(r) depends in general on the position as a consequence of Eq. (A7).

8
(i) IfR is diagonal in layer, i.e. proportionnal to τ 0/z , it commutes with A j (r) so that R(r) =R.
(ii) IfR is not diagonal in layer, i.e. proportionnal to τ x/y , it does not commute with A j (r) so that R(r) differs fromR. In that case, R(r) is modulated periodically over a distance of the order of the moire lattice constant. Indeed, we have with f 1 (r) = 1 3 3 j=1 cos(q j ·r), f 2 (r) = 1 3 3 j=1 sin(q j ·r). These results also apply to current-current quartic interactions, whereR is replaced by a vector of matricesM .

Appendix B: Diagrammatic technique for the non-interacting theory
To expand any observable in interlayer hoppings in the absence of interactions, it is not mandatory to resort to a field theoretical approach. We do so however, because it is useful for applying RG when interactions are included. To that end we need to introduce the Feynman rules specific to this unusual field theory. The free fermionic propagator associated to the action of the decoupled bilayer in Fourier space, where k is the momentum, Ω the Mastubara frequency, and we omitted the identity matrices σ 0 and τ 0 for simplicity. When drawing Feynman diagrams, we represent the free propagator (B1) with a solid line. We reserve q and ω for the internal momentum and Matsubara frequency and use k and Ω for external ones. Any correlation function can be written as an ensemble average ... 0 over S 0 . In particular for the time-ordered two-point function we have where ... 0 denotes the ensemble average over the quadratic action S 0 = S 0 + S α , which includes the hop-ping action S α = d 2 r dτ ψ † H α ψ, from which an expansion order by order in α can be carried out. The twopoint function (B2) is non-diagonal in momentum space, since S α reduces the continuous translational symmetry to the discrete translational symmetry over the reciprocal lattice R, which is the Z-module generated by the (linearly dependent) family of vectors {q j , j = 1, 2, 3}. For every vector b in R, we define the component G 0 (b, k, Ω) of the two-point function such that for all momenta k, q and frequency Ω. We focus on how interlayer hoppings renormalize the dispersion relation, so that we are mainly interested in the translational invariant part G 0 (k, Ω) = G 0 (0, k, Ω) of the fermionic propagator, represented in Fig. 8(f) as a double line. Successive interlayer hoppings that transfer the momenta (η 1 q j1 , ..., η m q jm ) in this precise order -where η 1 , ..., η m = ±, with the plus sign for a hopping to the top layer, and a minus sign to the bottom layer -give a nonzero contribution to G 0 (k, Ω) if the following conditions are met.
In particular, condition (ii) forbids odd numbers of insertions, so that all correlation functions can be expanded in α 2 , and entails that a hopping sequence is determined by the momenta and the sign of only the first hopping process η = η 1 . Joined with condition (i), it also yields that the transfer of a momentum at one point of the diagram must be followed by the transfer of the opposite momentum at another point. Thus we can join insertions of opposite momenta by a wavy line like in Figs. 8(a)-8(d).
We now introduce the translational part Σ α (k, Ω) of the self-energy as The contributions to Σ α (k, Ω) come from the connected two-point diagrams that conserve total momentum and that cannot be cut by one stroke into two subdiagrams that conserve themselves total momentum. Expanding Eq. (B2) to sixth order in α, we can decompose it as Σ α (k, Ω) = Σ 2 α (k, Ω) + Σ 4 α (k, Ω) + Σ 6,nes α (k, Ω) + Σ 6,row α (k, Ω) + Σ 6,cro α (k, Ω). In the following we use the shortcutη = −η with η = ±, and j, l, k = 1, 2, 3. The second order contribution ( Fig. 8(a)) reads The fourth order contribution ( Fig. 8(b)) reads Σ 4 α (k, Ω) = α 4 η,j =l Tη j G 0 (k + ηq j , Ω)T η l G 0 (k + ηq j − ηq l , Ω)Tη l G 0 (k + ηq j , Ω)T η j . order parameters are defined such that the mean-field Hamiltonian corrected by interlayer hoppings in momentum space takes the form The order parameters must satisfy self-consistent equations, which we treat as mutually independent, that is we put to zero all parameters other than ∆ 0 in the self-consistent equation satisfied by ∆ 0 , for example. Upon introducing an ultraviolet cut-off Λ, these equations read where to facilitate computation, the regularized integral Λ runs over a circle of radius Λ in Eq. (C15a) and (C15b), and over a square of side Λ in Eq. (C15c) and (C15d). The coefficients u 0 , u, w 0 and w arise from interlayer hoppings, and are the analogs of v for the terms σ z τ 0 , etc. As such we can compute them in a similar way. The results are where the wavefunction normalization N ψ is given in Sec. C 1. The integrands can be computed analytically but the self-consistent equations can also be solved numerically.

Hubbard-Stratonovitch decoupling
We aim at finding the most relevant insulating state near charge neutrality. It is therefore practical to decouple the interactions in the direct, particle-hole channel, to evince order parameters of the form ψ † M ψ , where the bracket ... denotes the ensemble average over the complete action S = S 0 + S int . Using Hubbard-Stratonovitch transformations, we introduce one auxiliary bosonic field for each interaction, whose ground state value in the correlated phase is a constant solution of the classical equation of motion. Thanks to this method, fewer diagrams need be computed, but the flows for each coupling constant artificially decouple, by discarding the ladder diagrams appearing in the corrections to the four-fermion vertex. The sixth order contribution splits into three terms. The first diagram hosts three nested hopping lines (Fig. 8(c)), The second diagram hosts two hopping lines in a row, embedded in a third one (Fig. 8(d)), The third diagram consists in three crossing hopping lines (Fig. 8(e)), Within a low-energy theory where k, Ω 1, we can further expand to order two in momentum k, and one in Matsubara frequency Ω, which results in where k = k x + ik y . If one keeps only the correction to the linear dispersion, the translational part of the fermionic propagator corrected by interlayer hoppings can be massaged into with the normalisation of the wave function N ψ = 1+3α 2 (1+β 2 )+2α 4 (1+7β 2 +4β 4 )+ 3 28 α 6 (8+16β 2 +376β 4 +187β 6 ), and the Fermi velocity dressed by interlayer hoppings as given in Eq. (3). We remind that v is expressed in units of the Fermi velocity v 0 of monolayer graphene. nian of the decoupled bilayer H 0 = iσ · ∂τ 0 , is invariant under the layer pseudospin rotational group U(1), and the Poincaré group R 1+2 O (1, 2), where indicates a semi-direct product. The hopping Hamiltonian H α = α 3 j=1 e −iqj ·r T + j + h.c, breaks Lorentz invariance, continuous space translations and layer plus pseudospin rotational symmetry. The symmetry group is thus reduced to the symmorphic space-time group R × (a 1 Z + a 2 Z) D 3 , composed of time translation R, discrete translations on the moiré lattice a 1 Z + a 2 Z, where a 1 and a 2 are the superlattice vectors, and the point group D 3 generated by the rotation C 3 around the z axis and the rotation C 2 around the x axis, Henceforth we disregard the translational symmetries and focus on the magnetic group generated by the point group D 3 and the two special "symmetries" IT and P . These operations act by conjugation on the singleparticle Hamiltonian H 0 in a four-dimensional representation (4d rep.), denoted as Γ, whose unitary matrix representation we now define.
The operation C 3 rotates the bilayer by an angle 2π/3 around the z axis perpendicular to the bilayer. Only the sublattice pseudospin is rotated, while the layer pseudospin is unaffected, so that Γ(C 3 ) = e 2iπ/3σz τ 0 . (C2) The operation C 2 rotates the bilayer by an angle π around the x axis of Fig. 1 at mid-distance between the layers. Both sublattices and layers are flipped, so that The composition of inversion I and time reversal T , denoted as IT , is antiunitary and represented by where K denotes the complex conjugation of the matrix elements. The operations R defined in Eq. (C2) to (C4) are pure symmetries of the Hamiltonian, which means that Γ(R) −1 H 0 (Rr, Rt) Γ(R) = H 0 (r, t) for R = C 3 , C 2 , and Γ(R) −1 H 0 (Rr, Rt) * Γ(R) = H 0 (r, t) for the antiunitary element R = IT . Finally, the unitary particle-hole operation P reverses the energy, and acts in real space as a reflection x → −x. Its matrix representation reads Following Ref. 32 we define P as a unitary operation in order to have a single antiunitary generator (IT ), unlike the convention of Ref. 57. This operation is an antisymmetry of the Hamiltonian, which means that Γ(P ) −1 H 0 (P r, P t) Γ(P ) = −H 0 (r, t). This antisymmetry is lost when the angular dependence of the kinetic terms σ ±θ/2 · k is kept, when terms quadratic in momentum are included in the single-particle Hamiltonian, or when intervalley scattering is permitted [57]. Since Γ is faithful rep., we can infer the multiplication table of the magnetic group from that of the matrix representation of the four generators.

Unitary group
The unitary groupD 3 generated by D 3 and P can be cast into the semi-direct product whereē = (P C 2 ) 2 and a barred operator represents the product of this operator byē. The classes of conjugation ofD 3 are given in Tab. III. The double group D D 3 = D 3 × {e,ē} is a normal subgroup ofD 3 , while D 3 is not. To find the irreducible representatinos (irrep.) ofD 3listed in Tab. IV-we can either find them from scratch using the composite operator method [33], or construct them by induction and other tricks from that of D 3 . Let us illustrate the second method.
The two 1d irrep. of the quotient groupD 3 /D D 3 = {e, P } generate the irreps A 1 and a 1 . The irrep. A 1 and A 2 of D 3 induce the reps A 1 ↑D 3 ∼ A 1 ⊕ a 1 ⊕ E 3 and A 2 ↑D 3 ∼ A 2 ⊕ a 2 ⊕ E 3 respectively, where ⊕ denotes a direct sum and ∼ the equivalence of rep. The commutator subgroup [D 3 ,D 3 ] is isomorphic to D 3 , whose index |D 3 /D 3 | = 4 gives the number of 1d irrep. Hence we have found all 1d irrep. The cardinal of the group being |D 3 | = 24, the remaining irrep. are five 2d irrep., including E 3 . We can decompose the 4d rep. Γ defined by Eq. (C2) to (C4) and (C5) Each column corresponds to a class of conjugation, and each line to an irrep. We use the symbols A and E prescribed by Mulliken's notation; the symbol a denotes a 1d irrep. whose character differs than one on antisymmetric operators.
where the remaining irrep. E 2 and E 4 can be found by orthonormality of the characters.

Magnetic group
The group generated byD 3 and IT is the dichromatic magnetic group To find the corepresentations (corep.) of M, we apply the Schur-Frobenius criterion [34][35][36][37][38] to each irrep. of D 3 . The Schur-Frobenius criterion states the following. For any irrep ρ ofD 3 , let us define the rep. ρ :D 3 → M n (C), R → ρ(IT·R·IT −1 ) * . We are necessarily in one of the three following scenarios. (i,ii) Either ρ is equivalent to its primed counterpart, in which case there exists an invertible matrix U such that ρ = U ρU −1 . (i) if U U * = ρ(IT 2 ), there is no Kramer degeneracy: the corep. issued from ρ has the same dimension as ρ and satisfies ρ(IT ) = ±U . (ii) If U U * = −ρ(IT 2 ), the corep. issued from ρ has twice the dimension. (iii) Or ρ is not equivalent to its primed counterpart, in which case ρ is necessarily equivalent to another irrep. ofD 3 , and the corep. has again twice the dimension of ρ, and coincides with ρ ⊕ ρ onD 3 . It turns out that all irrep. ofD 3 pertain to case (i), except E 1 and E 5 = E 1 , which fall into case (iii). In the former case, each irrep. leads to two corep., with ρ(IT ) = ±1 for the 1d irrep. or ρ(IT ) = ±σ x for the 2d irrep., where σ x represents here a generic Pauli matrix, but has nothing to do with the pseudospin. In the latter case, the one corep. formed by E 1 and E 5 is equivalent to the 4d rep. Γ. In the following, we write each corep. with an exponent ± to indicate whether the eigenvalues of ρ(IT ) are +1 or −1.
(α, β), which are the counterparts of the renormalized Fermi velocity for the matrix structures corresponding to those order parameters. They are calculated diagrammatically to sixth order in α and satisfy where the wavefunction normalization N ψ is given in Sec. B. In the main text, we plotted the order parameters corrected by interlayer hoppings, i.e. the quantities ∆ 0/z = N (∆) 0/z ∆ 0/z and G 0/z = N (G) 0/z G 0/z , which satisfy the self-consistency equations 13 where for simplicity we assume the shift momenta to be aligned along a crystallographic axis of the moire pattern, here along the y axis. The dimensionless functions F and F 0/z read where y ± = 2 + x(x ± 2) and z ± = 1 + 2x(x ± 1).
Appendix E: Renormalization

Hubbard-Stratonovitch decoupling
We aim at finding the most relevant insulating state near charge neutrality. It is therefore practical to decouple the interactions in the direct, particle-hole channel, to evince order parameters of the form ψ † M ψ for M ∈ {R i , M j }, where the bracket ... denotes the ensemble average over the complete action S = S 0 + S int .
Using Hubbard-Stratonovitch transformations, we introduce one auxiliary bosonic field for each interaction, whose ground state value in the correlated phase is a constant solution of the classical equation of motion. We must distinguish between the 1d corep., for which a scalar field φ i for i = 1, ..., 8, is sufficient, and the 2d corep., for which a two-component field ϕ j = {ϕ j,1 , ϕ j,2 } must be introduced, for j = 9, ..., 12. Such transformation enables to recast the action for quartic fermion interactions (5) into where the sum runs over both 1d and 2d irreps.. For simplicity we dropped the spatial and time dependences of the fields in Eq. (E1). The action for quartic fermion interactions thus splits into a bosonic quadratic action (the first term in the parenthesis), and a three-point vertex which takes the form of a Yukawa coupling (the second term in the parenthesis).

Renormalization procedure
The field theory described by the sum of action S 0 and action S Hub (E1) has critical dimension d c = 2, which entails that in d = 3 space-time dimensions, we expect that all interactions lead to quantum critical points that are perturbative in the small parameter = d − 2. From now onwards, we work in Fourier space and express all fields and integrals in terms of the momentum q and the Matsubara frequency ω. A general Fourier-transformed field is written φ q,ω for the bosonic case, or {ψ q,ω , ψ † q,ω } for the fermionic case, where ψ † q,ω = (ψ q,ω ) † denotes the conjugate of the Fourier-transformed field ψ q,ω . To renormalize the field theory, we assume that the complete action S is actually expressed in terms of bare fields {φ,ψ} for φ ∈ {φ i , ϕ j } and couplingsg for g ∈ {g i , λ j } , which are ill-defined in the interacting theory. The physical parameters and fields-written without the˚symbolare connected to their bare counterparts through the socalled Z constants.
We define define the Z constants for the fields such thatφ To regularize the theory, we work in an isotropic spacetime of dimension d = 2 + , and introduce a mass scale µ to make the regularized couplings dimensionless. A renormalized coupling g is linked to its bare valueg bẙ We included the normalization of the wavefunction N ψ in the redefinition of the couplings in order to compensate at all loop orders those arising from the corrected fermionic propagator G 0 , given in Eq. (B11). Owing to dimensional regularization, we must promote the Pauli matrices in S 0 to a Clifford algebra in arbitrary dimension d, satisfying the anticommutation rules {σ i , σ j } = 2δ ij for i, j = 1, ..., d. Using Eqs. (E2) and (E3), we find the renormalized action S R = S R,0 + S R,α + S R,φ + S R,int , where the quadratic, decoupled action reads The quadratic hopping action reads 14 The renormalized bosonic part of the action is given by Finally, the renormalized interaction between fermionic and bosonic fields is In Eqs. (E4) -(E7), we have ommitted dependence of the fields on frequency and used the shorthand To fix values of the Z constants we use the minimal subtraction (MS) scheme, i.e. we absorb in them only the divergent parts of the diagrams. Inspecting Eqs. (E6) and (E7), we see that the constant Z φ can be found from the divergences of the polarization, i.e. the bosonic selfenergy, while the constant Z g can be determined by absorbing the divergences of the three-point vertices. Before computing explicitely the diagrams for the polarisation and vertices let us outline the general strategy.

Preliminary mathematical remarks
Taking into account the interlayer hopping is done into two steps. We first draw the diagrams without insertion the hopping matrices, and then replace all solid lines by double lines. This corresponds to replacing free propagators by the propagators dressed by interlayer hoppings, like in the polarization shown in Fig. 9(a). In second step we include wavy lines, i.e. hopping matrices, connecting different propagators (see Fig. 9(b)). This splitting allows us to explicitly extract factors of v −1 , where v is the Fermi velocity corrected by interlayer hoppings (B12) and which vanishes at the first magic angle. We restrict our computation to order α 2 , which is the first non-trivial order.
When expanding product of matrices and integrating the trace, useful relations can be found in Refs. [58][59][60]. The Feynman trick, valid for any expressions A and B, enable to linearize products of denominators. For the four relevant interactions we consider here, the space-time integrals are isotropic and can be computed in arbitrary dimension d using for any reals a and b, and where Q = (q, ω) is the relativistic d-momentum. The dummy mass m → 0 plays the role of an infrared regulator and Γ denotes Euler's Gamma function, which satisfies for all real x and integer n; this relation is usually used with x = . In Eq. (E11), Ψ = (ln Γ) is Euler's Digamma function, which does not intervene at one loop, since we discard all finite quantities in the MS scheme.

Polarization
The one-loop polarisation Π i is the self-energy of the auxiliary field φ i (for i = 1, ..., 12). If ∆ i (k, Ω) denotes the corrected propagator of the bosonic field, we have Ω). The one-loop diagrams contributing to the polarization at zero external momentum k = 0 and fixed frequency Ω are drawn in Fig. 9. The polarization at order α 0 reads Notice that for the sake of generality, we will write all interaction matrices as the vectors M i , which can either denote a single matrix R i for i = 1, ..., 8, or a twocomponent vector M i for i = 8, ..., 12. The polarization at order α 2 reads The pole of the integral in Eq. (E12) per number of fermion flavors n (equal to four in our case), is given In Eq. (E13) we can use again the separation of energy scales : the theory is meaningful only at low energy, i.e. for q, ω 1, so that G 0 (q + ηq j , ω) can be replaced by The one-loop polarisation Π i is the self-energy of the auxiliary field φ i (for i = 1, ..., 12). If ∆ i (k, Ω) denotes the corrected propagator of the bosonic field, we have ∆ −1 i (k, Ω) = σ 0 τ 0 − Π i (k, Ω), where the identity matrix σ 0 τ 0 is FIG. 9. One-particle irreducible diagrams at one loop, up to order two in interlayer hoppings. The double line stands for the fermionic propagator corrected by interlayer hoppings of Fig. 8(f), the dashed line for the bosonic propagator, and the wavy line for the sum of interlayer hoppings of opposite momenta ±ηqj, for η = ± and j = 1, 2, 3. Polarisation Πi at zero external momentum and fixed Matsubara frequency, for the field φ i at order (a) α 0 and (b) α 2 . Three-point vertex V il at order (c) α 0 and (d-g) α 2 , whose hopping line is (d) internal, (e) external, (f) isolated and (g) crossed.
G 0 (ηq j , 0). This results in Π 1 i = −3α 2 χ i h i (β)Π 0 i where χ i equals either +1 for the interaction matrices σ z τ 0 and στ z or −1 for the interaction matrices σ z τ z and στ 0 ; and the corrugation-dependent function h i (β) equals either 1 − β 2 or 1 if the interaction matrix matches σ 0 or σ z in the pseudospin sector, respectively. This fixes the renormalization constant to (E15)

Vertices
We denote the one-loop contribution to the three-point vertex of interaction i renormalised by interaction l by V il . The one-loop vertices at zero external momentum k = 0 and fixed frequency Ω are drawn in Fig. 9(c) to 9(g) and computed in Eq. (E16) to (E20). For the vertices correcting an interaction i associated to a 2d channel, we write the vertex for only one component of the matrix M i , simply denoted as M i . The three-point vertex at order α 0 , given by diagram shown in Fig. 9(c), reads The three-point vertex at order α 2 (mixed diagram with two interlayer hopping) have either multiplicity one, or two. Those with multiplicity one nest either an internal hopping line, as in Fig. 9(d), η,j q,ω M l G 0 (q, ω)Tη j G 0 (q + ηq j , ω)M i G 0 (q + ηq j , ω)T η j G 0 (q, ω) · M l , or an external hopping line, as in Fig. 9(e), Tη j G 0 (ηq j , ω)M l G 0 (q, ω)M i G 0 (q, ω) · M l G 0 (ηq j , ω)T η j .
The mixed diagrams with multiplicity two nest either an isolated hopping line, as in Fig. 9(f), η,j q,ω Tη j G 0 (ηq j , ω)M l G 0 (q + ηq j , ω)T η j G 0 (q, ω)M i G 0 (q, ω) · M l , or a hopping line that crosses the interaction line, shown in Fig. 9(g), η,j q,ω M l G 0 (ηq j , ω)Tη j G 0 (q, ω)M i G 0 (q, ω) · M l G 0 (ηq j , ω)T η j . (E20) List of the functions f il (α, β) appearing in the RG flows of the four non-trivial channels a − 2 , a − 1 E + 2 , and E − 4 . The interaction matrices associated to each of these channels are indicated in the first line of the table.
The integrals J il are numerical constants, dependent of neither the number of fermion flavors n nor the corrugation parameter β, while K il (β) depends on the corrugation parameter. Using commutation relations between interaction and hopping matrices, we can express all vertices (E16) -(E20) in terms of J il and K il only. We then find the vertex renormalization constant to be (E23)
and v is the Fermi velocity (B12). We compute the RG flow equations by deriving Eq. (E3) with respect to µ at constant bare couplings. This yields