Twisted bilayer graphene at charge neutrality: competing orders of SU(4) Dirac fermions

We study possible patterns for spontaneous symmetry breaking in a Dirac fermion model, which is applicable to twisted bilayer graphene at charge neutrality. We show how a chiral SU(4) symmetry emerges and construct the corresponding low-energy model that includes a Fierz-complete set of symmetry-allowed four-fermion interactions. We employ an unbiased renormalization group treatment to identify the critical points that describe transitions into different ordered phases. The resulting phase diagram depends on the number of fermion flavours and we show that the coupling between ordering channels prevents many of the possible mean-field orders from being accessible at relevant, small flavour numbers. We argue that, as a consequence, twisted bilayer graphene is governed by a quantum Hall state or an SU(4) manifold of insulating spin-valley orders with emergent Lorentz symmetry that contains inter-valley coherent, spin Hall, and valley Hall states. We study how SU(4)-breaking perturbations affect the accessibility and can additionally stabilize symmetry-broken (semi-)metallic states.


I. INTRODUCTION
Spontaneous symmetry breaking of Dirac fermions plays an important role in many systems, ranging from the chiral phase transition in quantum chromodynamics to quantum critical points in semimetals [1][2][3][4] .Interactions must be sufficiently strong for such phase transitions to occur because the density of states of massless Dirac fermions vanishes at charge neutrality.Furthermore, if the phase transition takes place at zero temperature and is continuous, fermionic quantum critical behavior is expected which does not possess any classical analogue 3 .
As a well-known example, spontaneous symmetry breaking in graphene was intensely studied.However, interactions in graphene are estimated to be slightly too small to induce a phase transition 5 .
The recent discovery of strongly correlated moiré materials 6,7 provides a new opportunity for the investigation of quantum phase transitions in two-dimensional Dirac systems, where the relative interaction can be tuned via a twist angle [8][9][10] .This includes, for example, symmetrically twisted few-layer graphene systems 11 , or twisted Γ-valley transition metal dichalcogenides 12,13 .In particular, the most prominent moiré material -twisted bilayer graphene (TBG) -also hosts Dirac fermions at charge neutrality.Although the Dirac velocity vanishes in TBG at a so-called magic angle, strong interactions are expected in its vicinity where Dirac fermions remain present 14,15 .
In this work, we present a complementary analysis that investigates competing orders as instabilities of itinerant Dirac electrons.Focusing on the Dirac spectrum and employing an unbiased renormalization group approach allows us to draw universal conclusions and assess relevant orders.Our rationale comes from the observed signatures of a Dirac dispersion in TBG, and the relative reduction of interaction effects due to the vanishing density of states of Dirac electrons.
As a point of departure, we use an effective low-energy model valid around charge neutrality that respects the symmetries and topology of the narrow bands of TBG in the vicinity of the magic angle 48 .We argue that a chiral SU(4) symmetry emerges in this low-energy Dirac fermion model, which we generalise to an arbitrary number of flavours (Dirac points) N f .We study possible phase transitions based on all symmetry-allowed local interactions via the perturbative renormalisation group (RG).To do so, we analyse the fixed point structure of the RG equations for a Fierz-complete basis of the fermionic interactions.Our analysis is valid beyond the concrete application to TBG for other two-dimensional Dirac systems and thereby complements, e.g., previous fermionic RG investigations of phase transitions in graphene [49][50][51][52][53][54][55][56] .We classify possible quantum phase transitions based on flavour and Lorentz symmetry of their order parameter manifolds, and if they dynamically generate a mass that gaps out the Dirac spectrum or if they distort the spectrum but maintain (semi-)metallicity.
We find that, although there are many accessible transitions on the single-channel mean-field level, which is equivalent to large flavour numbers, the corresponding fixed points can disappear for small flavour numbers due to the coupling between different ordering channels.As a result, only a few transitions remain accessible for small flavour numbers, which is relevant to TBG.They include transitions with emergent Lorentz symmetry to a quantum anomalous Hall state and as well as a state with an SU(4) order parameter manifold that contains spin Hall, valley (spin) Hall and intervalley-coherent states with a gapped spectrum.We also determine how relaxing SU(4) symmetry affects the possible phase transitions.We show that the instabilities descending from the order parameter manifold of insulating states remain accessible separately.In addition, we argue that perturbing SU(4) stabilizes instabilities towards spin-(valley-)-polarized, spin nematic, and metallic intervalley-coherent states.

A. Dirac fermion model of twisted bilayer graphene
Our starting point is a low-energy effective model for Dirac fermions around charge neutrality in TBG 48 with the action at zero temperature where ψ is a 16-component spinor due to spin, sublattice, valley, and mini-valley degrees of freedom.We denote unity (γ = 0) and Pauli (γ = x, y, z) matrices in sublattice space with ρ γ and in valley space with τ γ .Analogously, we use σ γ and µ γ for spin and mini-valley space below.The sublattice and valley degrees of freedom form the Clifford algebra of the Dirac spinors (see below for a Lorentz-invariant formulation), while the kinetic part of the action S 0 is diagonal in spin and mini-valley space.The mini-valley degree of freedom originates from the Dirac points in the same valley of the two different graphene layers, i.e. the two mini-valleys per valley possess the same chirality.
We generalize the mini-valley degrees of freedom to an arbitrary number of flavours for the Dirac fermions ψ → ψ α with α = 1, 2, . . .N f .In the case of TBG N T BG f = 2, i.e., N f counts the number of eightcomponent spinors 108 .The Dirac fermion model (1)  implements the emergent C 2 symmetry of the continuum model of TBG via R C2 = ρ x τ x µ x , q → −q. and time-reversal via Θ = σ y τ x µ x K, q → −q.The threefold rotation C 3 is promoted to a full U(1) rotational symmetry R rot = exp[−iφρ z τ z ] in the low-energy limit around charge neutrality.Furthermore, since the hybridization between layers can be treated as a perturbation for states around the Dirac cones 15 , the spin-valley SU(2)×SU( 2)×U(1) symmetry of TBG is enhanced to an emergent SU(4) in the low-energy limit, similar to the chiral limit 38,57 .It is generated by the set of matrices with c = {x, y, z} and γ = {0, x, y, z}.The generalization of the flavour index leads to a unitary flavour symmetry acting via a transformation On the level of interactions, we include all symmetryallowed local four-fermion couplings according to the symmetries above.They will be generated by fluctuations, even if zero in the microscopic Hamiltonian so that it is important to include them for an unbiased analysis.
In particular, this means that we must include couplings that break the U(4)×U(4) symmetry of the interactions in TBG in the chiral limit because it is partially broken by the dispersion 30,[36][37][38] .We discuss effects that break chiral SU(4) symmetry further in Sec.III F. We also retain the emergent flavour symmetry in the interactions, i.e., different symmetry breaking patterns within mini valleys (which would translate to layer polarization or translational symmetry breaking on the scale of the moiré lattice) are degenerate.We then find that there are six distinct couplings allowed by these symmetries so that the interaction Lagrangian is given by where We assume the long-range Coulomb tail to be screened by the surrounding gates.The 64 matrices describing the couplings in Eq. ( 3) form a complete basis set for 8×8 matrices.Furthermore, we only included flavourdiagonal terms in L int because Fierz identities allow us to appropriately rewrite any flavor-symmetry-allowed terms as a linear combination of the ones considered in Eq. (3) (see App. A).Note that this includes terms with explicit mini-valley dependence (ψ † M µψ) 2 , where M ∈ {T i , ρ z τ z , ρ z τ z T i , ν, ω i }.Thus, the six couplings form a Fierz-complete basis for interacting SU(4)×SU(N f )symmetric Dirac fermions with N f > 1.The inclusion of all terms allowed by the Fierz identities is important for the correct identification of critical points 4,58 .For N f = 1, there are additional Fierz identities which reduce the number of couplings to three (see App. A).
In the following we study the full effective action We focus on particle-hole instabilities motivated by the experimental observations of semimetallic or insulating ground states at charge neutrality in TBG.Pairing instabilities are, in principle, also contained in our setup and can be analysed via the use of other Fierz identities that translate to the pairing channel.

B. Ordered states
Taken on their own, each of the couplings can induce an instability towards a state with spontaneously broken symmetry if strong enough.The ordered states are characterized by the condensation of the corresponding bilinears ϕ M = ⟨ψ α † M ψ α ⟩ ≠ 0 with M ∈ {T i , ρ z τ z , ρ z τ z T i , ν, ω i } (see also Tab.I).Note that n = ⟨ψ α † ψ α ⟩ is fixed by the density and does not break any symmetry.The condensation of ϕ ρzτz corresponds to the spontaneous formation of a quantum anomalous Hall state (QAH) and generates a gap in the spectrum.A finite ϕ ρzτzTi also gaps out the Dirac fermions.Its SU(4) order parameter manifold contains quantum spin Hall (QSH) ∼ ρ z τ z σ c , valley Hall (VH), and valley spin Hall (VSH) ∼ ρ z σ γ states, as well as inter-valley coherent order 29 ∼ ρ x τ x,y σ γ (IVC-1) that anticommutes with the Dirac Hamiltonian.In contrast, the SU(4) order parameter manifold described by ϕ Ti splits the degeneracy of the Dirac spectrum, resulting in a metallic state (except the splitting becomes on the order of the bandwidth).It contains spin-polarized ∼ σ c , valley-polarized ∼ τ z , and spin-valley-polarized ∼ τ z σ c orders, as well as inter-valley coherent order ∼ ρ y τ x,y σ γ (IVC-2) that commutes with the Dirac Hamiltonian.IVC orders can be classified into time-reversal symmetric T-IVC and nonsymmetric Kramers K-IVC 30 states.This classification generally depends on the projection in µ space.With their diagonal mini-valley configuration, IVC-1 is a T-IVC, and IVC-2 a K-IVC state.Instabilities towards other IVC states with triplet mini-valley configuration µ x,y,z can be probed via making use of the Fierz identities (see Sec. III C).Finally, condensation of ϕ ν and ϕ ωi spontaneously breaks rotation symmetry.This shifts the position of the Dirac points away from the corners of the Brillouin zone and preserves the semi-metallic state.A finite order parameter ϕ ν couples to the fermions in the same form like a vector potential ∼ (ρ x τ z , ρ y ) which leads to the integer quantum Hall (IQH) effect 59 .The SU(4) order parameter manifold ϕ ωi contains a spinful variant of this ∼ (ρ x τ z σ c , ρ y σ c ) (S-IQH), in addition to nematic ∼ (ρ x , ρ y τ z ) (NEM) and spin-nematic ∼ (ρ x σ c , ρ y τ z σ c ) (S-NEM) orders, as well as nematic and spin-nematic intervalley coherent states (N-IVC and S-N-IVC).The latter possess matrix order parameters built of two vectors under rotation (ρ z τ y σ γ , τ x σ γ ), (ρ z τ x σ γ , τ y σ γ ) which are related by valley U(1) v 60,61 .

A. RG flow and beta functions
To study the possible ordering tendencies, we employ a perturbative renormalization group scheme.Within the RG, we successively integrate out modes above a cut-off scale k and express the scale evolution of the couplings λ ∈ {g 1 , g 2 , g 4 , v 1 , v 2 , v 4 } via the differential equations β λ = k∂ k λ.This defines an effective action Γ k at scale k, and for k → 0, we recover the full effective action that includes all quantum corrections.We show in App.B in terms of a functional RG formulation [62][63][64][65][66][67][68] that the one-loop flow equations are independent of the cut-off scheme.For example, a Wilsonian scheme can be used that integrates over modes within a momentum shell.We obtain the flow equations where we rescaled the couplings k d−2 l f λ → λ with spacetime dimension d = 2+1 and loop integral l f (see App. B).For Wilson's momentum-shell cut-off, is the area of the unit sphere in d = 2 dimensions.For the Lorentz-invariant model, we again impose g

B. Fixed points and stability
We are interested in fixed points of the RG equations because they are connected to the possible quantum phase transitions that can be induced by strong couplings.Thus, we look for solutions λ * = (g * 1 , . . ., v * 4 ) where all the beta functions vanish To identify which solutions correspond to critical points, we consider the linearized flow of the beta functions around a fixed point and evaluate the stability matrix which describes how the scale evolution of the couplings is attracted to or repelled from the fixed point.The eigenvalues of the stability matrix determine the critical exponents of the corresponding second-order phase transition.They are universal quantities which do not depend on the microscopic details of the model.We are interested in stable fixed points that can be accessed by tuning only one parameter because they are associated to critical points.This means the spectrum of the stability matrix must have all negative eigenvalues except one, which defines a relevant repulsive direction in coupling space.This largest critical exponent determines the correlation-length exponent.The second largest exponent describes the corrections to scaling and decides over the stability of the fixed point.Fixed points with more than one relevant direction are considered multicritical or unstable.The trivial Gaussian fixed point (λ = 0), which defines a non-interacting theory in the IR, has a fully negative spectrum, reflecting the need of strong couplings to induce a phase transition.Any nontrivial solution describes an interacting fixed point with finite values of the couplings.If the bare couplings lie beyond the threshold set by the interacting fixed points, their flow diverges along the relevant direction towards the infrared indicating the formation of an ordered state with spontaneously broken symmetry (see Fig. 1).

C. Susceptibility analysis
The characterization of the possible ordered phases based on the stable fixed points is non-trivial for finite flavour numbers.Since the fixed points' coordinates generally have non-zero values, the identification of the ordering tendencies is complicated by multiple divergent couplings.In order to gain further insight within this fermionic description, we calculate the flow of susceptibilities.This is done by adding an infinitesimal test field to the effective action that explicitly breaks the symmetry along a specific ordering channel [69][70][71] where These terms represent the coupling of fermions to an external field, which is set to zero at the end.We define the corresponding susceptibilities as In the vicinity of the fixed points, the fields and susceptibilities scale according to with γ i = 2β i + 1 [72][73][74] (see App. C).If γ i < 0, the corresponding susceptibility diverges.In d = 2 + 1 dimensions, this means that β < − 1 2 .We associate the fixed point with the ordering channel whose susceptibility shows the strongest singularity (Fig. 4).
We extract the exponents from the flow equations for the infinitesimal fields at the fixed point solution.We can see that in the linear response regime the introduction of a symmetry breaking term along a certain channel renormalizes only the respective field with no feedback to the others, i.e the flow equations are decoupled so that we can evaluate them separately.Additionally, we note that these equations are Lorentz invariant if we again impose g = g 1 = −g 2 and v = v 1 = −v 2 .Analogously, to probe the susceptibilities for flavour-symmetry breaking, we introduce infinitesimal fields h t i ψ † M i µ j ψ and calculate the susceptibility exponents.In this case, the flow equations are expressed via the triplet couplings λ t i defined by the interaction terms We obtain the fixedpoint values of λ t i through a Fierz transformation of the fixed points in the flavour-diagonal basis (see App. A).

D. RG flow of the Lorentz invariant system
We first analyze the RG flow of the Lorentz invariant system, i.e. we consider the case where g = g 1 = −g 2 and v = v 1 = −v 2 .Our results for the system without Lorentz invariance and relaxed SU(4) symmetry are presented in the next sections.We perform the fixed point analysis described above for several values of the fermionic flavor number N f .In the large-N f limit, the RG analysis becomes equivalent to a single-channel, mean-field treatment.But for finite N f , the feedback between ordering channels is important and can qualitatively change possible ordering tendencies as we show below.We start with the large-N f limit, where the beta functions decouple and allow for an unambiguous characterization of the possible ordered phases.We identify four stable fixed points that in this limit are characterized by only one coupling of the scalar (S) or vector (V) channels being non-zero Upon lowering N f , the coupling between the different ordering channels starts to affect the fixed points and their properties in several ways.With the exception of FP S , other couplings become non-zero at the fixed points, necessitating the susceptibility analysis for an unambiguous identification of the corresponding instability.The scalar fixed point FP S is located at for general N f , and the vector fixed point FP V at with We give explicit analytical expressions in App.D. For FP S−SU(4) and FP V−SU(4) , the fixed-point values of all couplings become non-zero for general N f , see Fig. 2. We find λ * i,V−SU(4) ≥ 0 for all λ ∈ {g, v, g 4 , v 4 }, while g * S−SU(4) , v * 4S−SU(4) ≥ 0, and v * S−SU(4) , g * 4S−SU(4) ≤ 0. As we described above, the location of the fixed points determines the regime of strong coupling where the flow becomes unstable signalling an instability towards spontaneous symmetry breaking.The flow to strong coupling governed by FP S and FP S−SU(4) is shown in Fig. 1.
Importantly, we observe that not all of the fixed points are accessible for any N f .Specifically, we find critical values of the flavor number N c f at which they either disappear or at which they become multi-critical.These changes occur via a "collision" with other (multi-critical) fixed points.In the first case, the fixed point ceases to exist in the space of real-valued couplings at the critical N f , and instead a pair of complex conjugate solutions moves into the complex plane.We find that the fixed-point solution FP S−SU(4) exhibits this behavior and disappears at N c S−SU(4) = 1.89.In the second case, the second largest eigenvalue of the stability matrix θ 2 changes its sign at the critical value of the flavor number so that the fixed point becomes unstable.The largest eigenvalue of the stability matrix θ 1 = 1 independent of the critical point and N f as expected for critical four-fermion models 58 .The second largest critical exponent approaches θ 2 → −1 for large N f , but generally varies as function of N f .We show the exponent θ 2 for all four critical points in Fig. 3.The quantum critical points labeled as FP S and FP V−SU(4) remain accessible, i.e. θ 2 < 0, for all values 1 < N f < ∞, as does FP S−SU (4) in the range where it exists N f > 1.89.In contrast, FP V becomes multi-critical at N c V = 3.0 (see App. D for analytical expression).Furthermore, we underline the importance of including the analysis of the susceptibilities (Sec.III C) in the proper characterization of a fixed-point solution.
We find divergent susceptibilities at small values of N f ≈ 2 for the same instabilities as at large N f governed by critical points FP S , and FP S−SU (4) .These both correspond to states with a gap in the symmetry-broken regime.In the case of TBG, they describe a QAH state, or a state with SU(4) order parameter manifold connecting QSH, VH, VSH, and IVC-1 phases.In the case of FP V−SU (4) , no divergent singlet channel ∼ µ 0 exists.Instead we find an instability in the triplet channel, i.e., towards flavour symmetry breaking, ∼ ⟨ψ † ρ z τ z µψ⟩, which also opens a gap.In the case of TBG, this describes a QAH state with additional modulation on the moiré scale, which includes degenerate layer polarised states ∼ µ z and moiré density waves ∼ µ x,y due to the flavour symmetry.
E. RG flow of the Dirac fermion model for TBG Since a physical system like TBG is a priori not Lorentz-invariant, we are interested in extending the above discussion to the full Lagrangian spanned by all six couplings in Eq.( 3) and compare to what is captured in the Lorentz-invariant case.We thus look to find the fixed point solutions that describe quantum phase transitions for the group of beta functions in Eqs.(8-13).We find in total seven fixed points that are stable for some range of N f , but generally display a varied behavior as a function of the fermionic flavor number N f .In the large-N f limit, the RG equations again decouple and we can identify six of the seven stable fixed points corresponding to the separate ordering channels where all other couplings are of order O(1/N f ), respectively.We recover the fixed points FP S and FP S−SU(4) of the scalar channels from the Lorentz-invariant case.The two remaining Lorentz-invariant fixed points are unstable for N f → ∞.Instead, we find the four fixed points FP n , FP SM , FP M−SU(4) , and FP SM−SU(4) to be stable at large N f .FP n describes the density channel, FIG.4: The exponents β that determine the divergence of susceptibilities Eq. ( 18) for the vector SU(4) (FP V−SU( 4) ) and scalar (FP S ) fixed points in the Lorentz invariant case.If β < − 1 2 , the corresponding susceptibility diverges signaling a phase transition as marked by the orange dashed line.For the case of FP V−SU(4) the leading instability is in the triplet QAH channel.FIG.5: Schematic overview which summarizes the varied behavior of the different fixed point solutions as function of flavor number N f , where N f = 2 in the Dirac fermion model of TBG.Density (FPn) and semi-metallic (FP SM ) fixed points become multi-critical for small N f , the Lorentz-invariant vector fixed point (FP V ) is multi-critical for any N f , and the SU(4) vector fixed point (FP V−SU(4) ) at large N f (red dashed line).Metallic (FP M−SU(4) ), semi-metallic (FP SM−SU(4) ), and insulating (FP S−SU(4) ) SU(4) fixed points lie in the complex plane below a critical flavor number (blue dashed line).Only the scalar fixed point is stable for any N f .It is associated with a QAH phase transition.
which does not correspond to any symmetry breaking.A divergence of g 1 signals the singular response of the chemical potential to a density change at the Dirac point where the inverse compressibility diverges.The semi-metallic (SM) instabilities of FP SM and FP SM−SU(4) correspond to the breaking of rotation symmetry with or without an SU(4) order parameter manifold, and FP M−SU(4) describes polarizing instabilities which yield a metallic (M) state for small order parameters (see Sec. II B).
In a similar manner to the Lorentz invariant case, the existence and the stability of the fixed points can change when N f becomes small, and a single-channel mean-field description breaks down see App.D for the evolution of fixed-point couplings and critical exponents with N f ).Exceptions are the fixed point related to a QAH instability FP S given by Eq. ( 30), and the fixed point related to the density channel FP n given by g * 1 = N f /4(4N f − 1), which both possess just one non-zero coupling.A general overview of the of the fixed points' behavior as a function of N f is provided in Fig. 5.
We find that FP S and FP n remain accessible for any 1 < N f < ∞.As in the Lorentz-invariant case, FP S−SU(4) exists and is stable for 1.89 < N f < ∞.This means that Lorentz symmetry emerges in the vicinity of the associated phase transition towards the insulating SU(4) order.In contrast, the Lorentz-invariant fixed point FP V is always multi-critical if Lorentz invariance is not enforced.Interestingly, another fixed point with emergent Lorentz symmetry becomes critical when lowering N f , which is not part of the stable large-N f fixed points.The Lorentz invariant solution FP V−SU(4) collides with the semi-metallic SU(4) solution FP SM−SU(4) at N c f = 13.6 leaving FP SM−SU(4) unstable and FP V−SU(4) stable at smaller N f .At slightly smaller N f ≈ 12.89, FP SM−SU(4) and FP M−SU(4) collide and vanish into the complex plane.Finally, FP SM becomes multi-critical at N f ≈ 2.55.We also calculate the exponents characterizing the divergence of susceptibilities for the seven fixed points (see App.C).We find that all the divergent channels can be connected to the large-N f limit, with the exception of FP V−SU(4) , which is not stable at large N f .Furthermore, slightly before FP SM−SU(4) disappears, the related leading instability changes to the triplet channel (see Fig. 7 in App.C).
Therefore, a picture similar to the Lorentz-invariant case arises in the Dirac fermion model of TBG at small N f ≈ 2. We obtain instabilities related to FP S FP V−SU (4)  and FP S−SU (4) .These correspond to symmetry-broken phases with a gap due to QAH states with (∼ µ 0 ) or without (∼ µ x,y,z ) mini-valley symmetry, or SU(4)-related VH, VSH, QSH and IVC-1 states.Interestingly, there is neither an instability towards polarizing orders (SP, VP, IVC-2) associated with FP M−SU(4) nor towards nematic orders associated with FP SM or FP SM−SU(4) in the strict SU(4)-symmetric model.
If however, this symmetry is reduced, additional solutions become accessible also in the small-N f limit because the critical flavour number of the separate instabilities is lowered.This is further elaborated on in the next section.

F. Relaxing SU(4)
As described in Sec.II A, the Dirac fermion model of TBG possesses an emergent SU(4) symmetry, leading to a fourfold degeneracy of the Dirac fermions at charge neutrality.Since three out of the six quantum critical points are related to orders containing an SU(4) symmetric manifold, we are interested in studying whether these fixed-point solutions are robust under a possible breaking of this symmetry.The breaking can occur on the level of the interactions due to high-energy corrections above the UV cutoff of the Dirac fermion model, which lower the symmetry from SU(4) to SU(2)×SU( 2)×U(1), i.e., independent spin rotations in both valleys.
We thus split the initial manifold into five new channels that are invariant by the lowered symmetry T (1) = {σ c } T (2) = {τ z σ c } T (3) = τ z T (4) = {ρ y τ x , ρ y τ y } T (5) = {ρ y τ x σ c , ρ y τ y σ c } and study the Lagrangian that accounts for the new interaction terms.The above channel separation means that five new couplings are to be introduced for each term ∝ v i that contained an SU(4) manifold, yielding in total eighteen interaction terms in the lower symmetry Lagrangian where ω We derive the flow equations for the new set of couplings (see App. E) and investigate how the fixed-point structure changes.
To study the robustness of the SU(4) symmetry, we evaluate the eigenvalues of the stability matrix at the coupling coordinates of the SU(4) symmetric quantum critical points.We find that the fixed points related to FP S , FP n and FP SM retain the same behavior as in the previous cases since they are not related to any SU(4) manifold.Furthermore, we consider the stability matrix for N f = 2, 14, 20 for the SU(4) fixed points FP M−SU(4) , FP SM−SU(4) and FP S−SU (4) , which all become multi-critical.The number of relevant directions is increased to four for FP S−SU(4) and FP SM−SU(4) and to five for FP M−SU (4) .They are primarily along couplings of the channel they descend from, i.e along v (i) 1 for FP M−SU(4) , and along v (i) 4 for FP S−SU (4) .As such, we observe that the emergent SU(4) symmetry of the Dirac fermion model is not robust against perturbations.Mechanisms which induce a breaking of the SU(4) symmetry select a specific ground state out of the SU(4) order parameter manifolds.
To determine possible selections, we follow the N fevolution of the eighteen single-channel solutions that occur as quantum critical points at large N f in this case.The result is summarized in Tab.I.In total, out of the eighteen large-N f stable fixed points we find that seven remain critical down to small N f ≈ 2.Besides FP S and FP n , fixed-point solutions characterizing the transition to QSH, VH, VSH, and spinless IVC-1 orders are stable.These originate from the FP S−SU(4) fixed point solution in the SU(4) symmetric case.Lorentz invariance is again emergent at these solutions.Furthermore, a fixed point that originates from inaccessible FP M−SU( 4) is now stable for N f = 2 in the SU(2)×SU(2)×U(1) case.It is related to the transition to a spin-valley polarized state.Interestingly, the critical N c f for several other fixed points that describe transitions to orders previously related by SU(4) symmetry is also considerably lowered.The effect is the strongest for spin-polarization, IVC-2, and N-IVC order.Their critical N f as part of the SU(4) order parameter manifold and  N c f ≈ 5.2, respectively.In this sense, these ordering tendencies are stabilized by perturbations that break SU(4).

IV. CONCLUSION
In the present work we studied competing orders in the Dirac fermion model of TBG at charge neutrality in a universal, unbiased way.As such, our considerations generally apply to Dirac fermions with approximate SU(4) symmetry in 2 + 1 space-time dimensions.
We determined a Fierz-complete set of all symmetrypreserving, local 4-Fermi interactions, which we could classify in eighteen, six, or four different channels in the case of SU(2)×SU( 2)×U(1), chiral SU(4), or Lorentz in-TABLE I: Summary of possible orders introduced in Sec.II B and their tensor structure specified by Pauli matrices in sublattice ρ, valley τ and spin σ space (c = x, y, z).We also list the corresponding single-channel coupling in the limit of large flavour numbers N f and the critical values N c f in the SU(2) × SU(2) × U(1) case.The dash indicates stability for the whole interval of N f studied.Double lines group orders related by SU(4) symmetry Eq. ( 2).

Order
Tensor structure Large- 3.66

S-NEM
19.9 variant SU(4) symmetry, respectively.We calculated the perturbative RG flow of the couplings and investigated their fixed-point structure with the aim of identifying quantum critical solutions that are associated with instabilities towards different ordered states.We diagnosed the instabilities by a divergence in the corresponding susceptibilities.
The model contains a control-parameter in the form of the fermion flavor number.This allowed us to investigate the interplay of multiple interaction channels, which becomes important in the small-N f regime relevant to TBG.We connected the results to the large-N f limit, where the RG equations decouple and reduce to a single-channel mean-field approach.
We found a rich landscape of critical fixed points, which display a varied behavior as a function of N f .Importantly, we showed that the inter-channel feedback makes many of the single-channel mean-field solutions unstable at small N f , either because the solutions disappear or because they become multi-critical.We note that, while the critical values of N f where this happens cannot be determined quantitatively on the oneloop level, the qualitative findings which solutions become unstable are usually robust.
We showed that the instabilities that gap out the Dirac spectrum, and thus can gain much condensation energy, are particularly stable at any N f and symmetry setups.They correspond to quantum anomalous Hall, quantum spin Hall, valley (spin) Hall, and time-reversal symmetric inter-valley coherent states in the singlet sector, or a QAH together with moiré density waves or mini-valley polarisation in the triplet sector.These solutions possess emergent Lorentz symmetry and thus can already be described via a Lorentz-invariant version of the Dirac fermion model of TBG.The fixed point associated with QAH order remains stable through the entire range of 1 < N f < ∞, as well as upon breaking SU(4) symmetry.The fixed point related to a 'triplet QAH' instability is only stable below N f < 13.6 with SU(4) symmetry.Order parameters of QSH, VH, VSH, and IVC-1 form a degenerate manifold in the SU(4)-symmetric case since they are related by symmetry.
The fixed point of the SU(4) manifold FP S−SU(4) disappears at N f ≈ 1.89 and we found it to be unstable when SU(4) is broken.However, the separate solutions associated with QSH, VSH, VH, and gapped IVC that descend from this channel in the SU(2)×SU( 2)×U(1) case remain stable in the entire range 1 < N f < ∞.The only exception is the fixed point that describes a phase transition to a spinful, gapped IVC state, which disappears at N f ≈ 2.75.
In addition, we found that relaxing SU(4) symmetry stabilizes several orders which generically do not gap out the Dirac spectrum.In particular, the critical flavour number of spin-(valley-)polarized, metallic intervalley coherent, and spin nematic ordering tendencies is strongly reduced.Among them, a spin-valley-polarized state is stable at N f = 2.
With respect to selecting the ground state in TBG, we argued that it is crucial to account for the competition between ordering tendencies because several symmetrybroken states lie close in energy.In this context and in light of recent experiments 28,100 , it is interesting to note that the leading instability is not necessarily the one expected based on symmetries in the mean-field picture [30][31][32][41][42][43][44] due to inter-channel renormalizations of the interaction. It wil be enlightening to extend our unbiased treatment by including the breaking of additional symmetries either spontaneously or externally in future studies.In particular, it was argued that strain plays an important role and selects the so-called incommensurate Kekulé spiral (a time-reversal symmetric IVC phase with incommensurate wave vector q ≠ 0) as the ground state away from charge neutrality 28,42,43,100,101 .Furthermore, the near-degeneracy of symmetry broken states also suggests an investigation of phases of coexistence.[102][103][104][105][106] Appendix A: Fierz identities In principle, the generalized flavour symmetry permits terms of the form ∑ 2 in the interacting Lagrangian, where Ψ = (ψ 1 , . . ., ψ N f ) T is a collective 8N f -component spinor, M is an 8 × 8 matrix acting in spin-valley-sublattice space, and κ i are generators of SU(N f ), which replace the Pauli matrices µ c of mini-valley space for general N f .To make the flavour index explicit, we can use the completeness relation ∑ i κ αβ i κ γδ i = N f δ αδ δ βγ − δ αβ δ γδ and rewrite with 8-component spinors ψ α of flavour α.Furthermore, we can connect the terms non-diagonal in flavour (ψ α † M ψ β )(ψ β † M ψ α ) to the flavour-diagonal ones via Fierz identities using that the 64 matrices M form a basis for 8 × 8 matrices.The general form of a Fierz identity is (A2) The above equation effectively means that we can write any fermionic bilinear as a linear combination of others provided we know the coefficients F XY .To see how diagonal and non-diagonal flavour terms are related, we need to set a = b and c = d.Thus, the results we provide in the text can be translated to a basis which contains off-diagonal flavour terms.
The Fierz identities can be condensed in the form of a matrix whose entries are exactly the coefficients F XY .By labelling the 4-Fermi terms based on the six possible matrix channels M , we can define and Eq.(A2) can be cast in matrix form Note that F −1 = F .Using the completeness relation and the Fierz matrix, we can also rewrite X s using only triplet terms In the special case where N f = 1, the Fierz matrix relates flavour-diagonal terms and the number of independent couplings can be reduced to three, which is the degeneracy of the eigenvalue 1 of the Fierz matrix F .
For the Lorentz-invariant case, where the number of independent couplings is reduced to four, the Fierz identity is given by and Our starting point is the Wetterich equation 107 where ∂ t = k d dk and Γ (2) k is the matrix of functional derivatives of the effective action with respect to the fields, defined as with Ψ = (ψ, ψ † T ) T .The essence of Eq.(B1) is that it describes the evolution of an effective action as a function of a scale variable k ∈ [0, Λ] with Λ being a UV-cutoff.This evolution is encoded in the regulator function R t , which defines the way fluctuations in the interval [k, Λ] are integrated out.At k = 0, all fluctuations are integrated out, which yields the full quantum effective action and the complete solution to the problem.
Based on the discussion of Sec.II A we make an ansatz for the scale dependent action by introducing a scale dependence on all interactions and the regulator function.We neglect the wavefunction renormalization coefficient as diagrams that contribute to the anomalous dimension in a purely fermionic 1-loop approximation vanish in the regime of point-like interactions.To extract the beta functions for the scale-dependent couplings, we first redefine the scale derivative to only act on the regulator and rewrite Eq. (B1) as We split the inverse full propagator into a fieldindependent part Γ (2) 0 and ∆Γ k which incorporates the effects of fluctuations.We expand the logarithm around this point We then insert our ansatz and compare right and left hand side of Eq. (B4) to identify the beta functions for the running couplings.This is analogous to a 1-loop, Wilsonian RG scheme for a specific choice of the regulator function.However our analysis will turn out to be independent of the choice of the regulator.Following the procedure detailed above the right hand side of Eq (B4) evaluates to where p labels the internal momentum integration variable, M i,j are any of the 64 matrices considered in the interacting Lagrangian and G k is the scale dependent fermionic propagator defined as where r(k) is a dimensionless function encoding the regulator scheme R k = (iω + q x ρ x τ z + q y ρ y )r(q 2 /k 2 ).Evaluating the matrix products in Eq. (B5) yields the beta functions for the six couplings.
The quantity l f defined as is the threshold function.By rewriting our couplings as λi = λ i /l f , the threshold function is absorbed in the rescaling and thus the beta functions become regulatorindependent.

Appendix C: Susceptibilities
As mentioned in the main text, to gain insight into the ordered phases related to the quantum critical points, we analyze the divergence of the test-vertex susceptibilities.We start from an effective action The quantities r and χ will flow according to: where the flow of Γ k is given by Eq. (B1).We perform an expansion of the effective action, while noting that the field-independent part now gets a contribution from the inclusion of the linear vertex term.Keeping only terms that contribute to the flow of the above quantities we get where C ij is a matrix containing all constant prefactors for each term that appears in the beta functions and l f is given by Eq. ( 6).Eq.C2 is explicitly given by: We rescale the couplings λ = k d−2 λl f , with d being the space-time dimension of the model.At the fixed point solutions, λi = λ * i , we can solve the differential equations to get the explicit dependence as a function of k: We define β i ≡ C ij λi .Then, the dependence of the susceptibility can be written as We thus relate the two exponents via γ i = 2β i + d − 2. Setting d = 3, we get the condition for a divergent susceptibility used in the text β i < − 1 2 .To determine the tendency for breaking of flavour symmetry, we repeat this analysis in the triplet basis, i.e., we introduce test vertices r t i h t i ψ † M i ⊗κ j ψ and use the Fierz identities to express the fixed point couplings via combinations of triplet couplings 2 .Replacing M i → M i ⊗ κ j in Eq. (C3), we find that ∂ t r t i can be obtained from ∂ t r i by replacing singlet via triplet couplings λ i → λ t i and a relative sign between the two terms on the right hand side in Eq. (C3).Explicitly we obtain The treatment described above is condensed in Fig. 4 in the main text and Figs.6 and 7 in this appendix for the Lorentz invariant model (Eq.6) and full SU(4) model (Eq.3) respectively.We reiterate the importance of this analysis to identify which of the stable fixed points can describe a quantum phase transition.Specifically, as can be seen in Fig. 7 none of the singlet susceptibilities of FP V−SU(4) diverges in the regime of N f where it is stable.However, we identify a divergent susceptibility for FP V−SU(4) that is associated with a triplet QAH state ⟨ψ † ρ z τ z µψ⟩, i.e. a QAH state with additional modulation (µ x,y ) or polarisation (µ z ) on the moiré scale (see Fig. 7).
FIG. 8: An example of a complex conjugate collision of the scalar SU(4) fixed point FP S−SU(4) with another multicritical fixed point leading to the creation of a pair of real valued fixed point solutions.The points correspond to the coupling values for different N f .FIG. 9: The couplings of the scalar (FP S ) and vector (FP V ) SU(4) fixed points for several values of the fermion flavor number N f .As mentioned in the main text, for general N f , several fixed-point coupling values are non-zero so that different ordering channels are coupled.For large values of N f only one coupling is non-zero and a single-channel, mean-field description is possible.Out of the four stable, Lorentz invariant fixed point solutions, FP V does not emerge as a quantum critical point in the non Lorentz-symmetric model.

FIG. 1 :
FIG.1: Flow of couplings according to the beta functions βg 4 , βv 4 and βv 1 along a plane in coupling space illustrating the behavior around the critical fixed points FP S , FP S−SU(4) and FP M−SU(4) (red) and the trivial non-interacting Gaussian fixed point (black) for the Lorentz noninvariant case.If the bare values of the couplings are larger in magnitude than the values set by the critical fixed points, they flow to strong coupling in the infrared signaling a potential instability towards symmetry breaking.

FIG. 2 :
FIG. 2:The couplings of the scalar (FP S−SU(4) ) and vector (FP V−SU(4) ) SU(4) fixed points for several values of the fermion flavor number N f .For general N f , all fixed-point couplings are non-zero demonstrating that different ordering channels are coupled.For large values of N f only one coupling is non-zero and a single-channel, mean-field description is possible.