Emergence and stability of spin-valley entangled quantum liquids in moir\'e heterostructures

Twisting moir\'e heterostructures to the flatband regime allows for the formation of strongly correlated quantum states, since the dramatic reduction of the bandwidth can cause the residual electronic interactions to set the principal energy scale. An effective description for such correlated moir\'e heterostructures, derived in the strong-coupling limit at integer filling, generically leads to spin-valley Heisenberg models. Here we explore the emergence and stability of spin liquid behavior in an SU(2)$^{\mathrm{spin}}\otimes$SU(2)$^{\mathrm{valley}}$ Heisenberg model upon inclusion of Hund's-induced and longer-ranged exchange couplings, employing a pseudofermion functional renormalization group approach. We consider two lattice geometries, triangular and honeycomb (relevant to different moir\'e heterostructures), and find, for both cases, an extended parameter regime surrounding the SU(4) symmetric point where no long-range order occurs, indicating a stable realm of quantum spin liquid behavior. For large Hund's coupling, we identify the adjacent magnetic orders, with both antiferromagnetic and ferromagnetic ground states emerging in the separate spin and valley degrees of freedom. For both lattice geometries the inclusion of longer-ranged exchange couplings is found to have both stabilizing and destabilizing effects on the spin liquid regime depending on the sign of the additional couplings.

Several model constructions for correlated moiré materials have been put forward in terms of effective tightbinding descriptions on the moiré superlattice, augmented by various interaction terms [36][37][38]. Whereas details of the models may differ, they feature a series of universal traits: (1) an emergent hexagonal superlattice, (2) a multi-orbital structure, and (3) Hund's and extended Hubbard interactions. More specifically, while TBG is preferably described using a honeycomb superlattice [36][37][38], related structures such as twisted doublebilayer graphene (TDBG) or trilayer graphene/hexagonal boron nitride heterostructures (TLG/h-BN) are better captured by a triangular superlattice [6,38,39]. The orbital degrees of freedom are inherited from the valleys in the original bands, e.g., the two Dirac valleys in the Brillouin zone of graphene.
These universal aspects can be combined into a minimal model, with a two-orbital extended Hubbard model [16,17] serving as a paradigmatic starting point. Its kinetic term H t = −t ij 4 α=1 (c † iα c jα + h.c.) for the electrons combines the spin projection s ∈ {↑, ↓} and valley quantum number l ∈ {+, −} in a flavor index α ∈ {(↑,+), (↑,−), (↓,+), (↓,−)}, reflecting an effective SU(4) symmetry. On the triangular lattice this results in a set of four degenerate bands, which can describe, e.g., the set of minibands above charge neutrality in TDBG or TLG/h-BN. On the honeycomb lattice, with its additional sublattice degree of freedom, this results in eight bands, two sets of four degenerate bands, which describe the minibands above and below charge neutrality in TBG [40].
The simplest conceivable interaction term, which also retains the SU(4) symmetry, is a Hubbard interaction H int = U i ( 4 α=1 n iα ) 2 , which can arise in the limit of large lattice periods where the interaction depends primarily on the total charge on a site and becomes the dominant interaction scale. In this strongcoupling limit, the kinetic term can then be treated perturbatively [16,17,38]. With an integer number of electrons per site this leads to an effective spin-valley Heisenberg Hamiltonian with SU(4) symmetric superexchange coupling J H ∝ t 2 /U . Additional symmetry-breaking interactions are, however, to be expected, in particular in the form of Hund's-type couplings in either the spin or valley degrees of freedom [16,17]. Moreover, Wannier state constructions suggest that further-neighbor interactions can become sizable [36] and should augment any minimal model.
In this work, the above considerations naturally lead us to explore a nearest-neighbor spin-valley Heisenberg model with SU (2) spin ⊗ SU ( (4) spins. This is in contrast to the four-dimensional fundamental representation of SU(4) relevant to, e.g., the case of quarter-filling.
For both lattice geometries, we find extended parameter regimes surrounding the SU(4) symmetric point where no long-range symmetry-breaking order occurs, indicating a stable realm for a spin-valley entangled quantum liquid. Moving further away from the SU(4) symmetric point, we find magnetic order in the spin and valley degrees of freedom that can be either antiferromagnetic or ferromagnetic, providing a possible explanation for the clear signatures of spin-polarization observed in TDBG [9][10][11]. To explore the effect of longerrange interactions, we augment our model by a nextnearest neighbor coupling and determine its role in stabilizing quantum spin-valley liquid (QSVL) behavior versus long-range order for different signs of the coupling and the two lattice geometries. Our work complements earlier work for the case of quarter-filling, where it was argued that a QSVL state with neutral gapless fermionic excitations forms on the honeycomb lattice [41], while on the triangular lattice extended parameter regimes without any net magnetization have been identified in DMRG simulations [42].
Spin-valley model. The starting point of our study is an SU(4) spin-valley Heisenberg model [16,38], where J H is the antiferromagnetic exchange coupling between nearest neighbors on either the triangular or honeycomb lattice, and T i denote SU(4) spins. The µ = 1, . . . , 15 components of the spin operators can be represented on a fermionic Hilbert space via the parton constructionT µ i = f † iα T µ αβ f iβ , where the index α enumerates four different fermion flavors and the matrices T µ are the SU(4) generators [43]. At half-filling of the underlying Hubbard model, the local spin-valley Hilbert space is sixdimensional (4 choose 2), which leads to a local filling constraint of two partons per lattice site α f † iα f iα = 2. Upon inclusion of Hund's couplings, the SU(4) symmetry of the model is explicitly broken [16]. Omitting other sources of SU(4) breaking, a residual separate spinvalley SU (2) s ⊗ SU(2) v symmetry remains which is reflected by the extended Hamiltonian Instead of enumerating the four fermion types by a single index, we have exposed the spin quantum number s ∈ {↑, ↓} and the valley quantum number l ∈ {+, −} explicitly; Pauli matrices are denoted by θ a , a ∈ {1, 2, 3}. At the high-symmetry point J = J s = J v the full SU(4) symmetry is restored. We assume that the Hund's interactions are weak enough such that all exchange couplings are antiferromagnetic [42], i.e. J, J v , J s > 0.
Pseudofermion functional RG. Parton-decomposed quartic Hamiltonians of the general type defined in Eq. (1) can readily be analyzed by the pseudofermion functional renormalization group (pf-FRG) [44,45]. For SU(N ) spins, the approach is already naturally formulated with a local constraint of N/2 fermions per site. It combines aspects of an expansion in spin length S [46] (which naturally favors magnetic order) and in the SU(N ) spin symmetry [47,48] (which typically favors quantum spin liquid states), and it becomes exact on a mean-field level in the separate limits of large S and large N . It is thus suited to resolve the competition between ordered ground states and QSVL phases in the spinvalley model at hand. We extend the standard implementation of pf-FRG to incorporate the SU(2) s ⊗ SU(2) symmetry, thereby obtaining flow equations for the oneparticle irreducible vertices as a function of an RG frequency cutoff scale Λ [40]. Numerically solving the set of O(10 6 ) flow equations at up to 84 Matsubara frequencies and using a real-space vertex truncation of L = 7 lattice bonds in each spatial direction, spontaneous symmetry breaking, e.g., the onset of long-range magnetic or valence bond order, is indicated by an instability of the RG flow [44,49] which occurs at some critical scale Λ c .
In the case of long-range order, to identify the precise nature of the ordered state we can separately gain access to the elastic component (ω = 0) of the correlation functions in the spin sector and in the valley sector, Sharp features emerging in the respective structure factors ij allow us to deduce the type of long-range order in either the spin or the valley degrees of freedom, cf. Figs. 1 and 2.
Emergent spin-valley liquid behavior. We begin our analysis with the SU(4) symmetric point, J s /J = J v /J = 1. For both the triangular and honeycomb lattice, no instabilities are detected in the pf-FRG flow, indicating a fully symmetric ground state. In addition, upon varying the vertex range L we observe no finite-size dependence of the flows, consistent with a ground state without symmetry-breaking long-range order (see [40] for details on the finite-size analysis). This rules out not just magnetically ordered states, but also valence bond or dimer crystals [50], an ordering which spins in the self-conjugate representation are often prone to [51,52].
For SU(4) spins in the self-conjugate representation we can further use the Lieb-Schultz-Mattis-Hastings [53][54][55] theorem to rule out a featureless Mott insulator as the ground state in the case of the triangular lattice, whereas such a state is in principle still a possibility on the honeycomb lattice. We note that the spin/valley structure factors have features resembling 120 • /Néel order, albeit significantly broadened, see Figs. 1(c) and 2(c).
Stability of spin-valley liquid and adjacent magnetism. Moving towards parameter regimes with broken SU(4) symmetry, J s /J, J v /J = 1, we find that an extended paramagnetic region emanates from the SU(4) symmetric point, see the white wedges in Figs. 1 and 2. Importantly, this finding supports the stability of the emergent spin-valley liquid behavior even in the presence of SU(4) breaking perturbations such as the Hund's coupling. Comparing the two lattice geometries, the triangular lattice gives rise to a parametrically larger QSVL phase than the bipartite honeycomb lattice, which can likely be traced back to the geometric frustration of the former. Along the diagonal line of equal coupling J v = J s , the QSVL region eventually collapses and disappears, being replaced by long-range antiferromagnetic order. Moving along the dotted diagonal line in the respective phase diagrams we observe a strongly suppressed breakdown scale Λ c , relative to the surrounding parameter space, indicating that quantum fluctuations are strongest when J v = J s .
For sufficiently strong dominance of either spin or valley coupling, different ordered phases occur for both lattice geometries. The transition towards an ordered state  is indicated by a leading instability in the RG flow, either in the spin or valley sector. To explore the subleading instabilities in the remaining sector, we employ a heuristic mean-field-like approach to estimate the effective spin or valley couplings between nearest-neighbor sites i and j, Note that for 120 • or Néel order in one of the SU(2) sectors the corresponding nearest-neighbor correlation becomes negative. Therefore, the effective couplings J eff v and J eff s may, too, turn negative and drive a ferromagnetic instability in the other sector, despite the antiferromagnetic nature of all couplings in the microscopic spinvalley model [42]. This kind of mechanism may be at the origin of the spin polarization observed at half-filling in TDBG [9,10], as first pointed out in Ref. [42] for quarter-filling. Extracting the sign of the effective coupling according to Eq. (3) at the transition scale of the leading sector allows us to distinguish two regimes with either ferro-or antiferromagnetic correlations in the subleading sector [56]. In Figs. 1 and 2 the so-determined order in the subleading regimes is indicated by triangle (ferromagnetic) or square (antiferromagnetic) symbols.
Longer-range interactions. In the ongoing search for an effective microscopic description for moiré het-erostructures it has been pointed out that longer-ranged Coulomb interactions should not be neglected [36], which in the effective spin model will give rise to exchange couplings beyond nearest-neighbor. To probe the stability of the QSVL regime in our model we here consider the effect of a next-nearest neighbor coupling J 2 .
Let us first recapitulate the effect a next-nearestneighbor coupling J 2 for the spin-1/2 SU(2) case on the triangular and the honeycomb lattices. Here the bare nearest neighbor coupling leads to magnetic ordering and only an antiferromagnetic J 2 of intermediate coupling strength facilitates the formation of a narrow quantum spin liquid (QSL) regime [57,58], as indicated by the grey boxes in Fig. 3. Notably, the induced QSL regime is somewhat larger for the honeycomb lattice where the next-nearest neighbor interaction introduces geometric frustration.
For the model at hand, we first concentrate on the SU(4) symmetric point and explore the effect of J 2 /J 1 ∈ [−1, 1]. As shown in Fig. 3, the QSVL region for the SU(4) model is significantly expanded for both lattice geometries in comparison to the SU(2) QSL case. The impact of J 2 on the the full spin-valley (J s , J v ) phase diagrams of Figs. 1 and 2 is illustrated in Fig. 4 for both ferromagnetic and antiferromagnetic J 2 . While an antiferromagnetic J 2 is found to further widen the wedgeshaped QSVL region, the converse occurs for ferromagnetic J 2 , which drives the system closer to the ordered states. This means that, depending on the sign of J 2 , longer-range interactions can actually stabilize and even expand the region of QSVL behavior.
Conclusions. In this work, we studied SU(2) s ⊗ SU (2) v -symmetric spin-valley Heisenberg models in the self-conjugate representation for both the triangular and honeycomb lattice. Seen as the effective Hamiltonians generated in the strong-coupling limit of an underlying Hubbard model, such models are relevant as minimal models in the exploration of the correlated insulating states of recently synthesized moiré heterostructures, including TBG (honeycomb) or TDBG and TLG/h-BN (triangular). Depending on which set of minibands the Hubbard model is designed to describe, the half-filling case studied here can potentially describe different candidate correlated insulators [40], e.g. the insulator at halffilling n = +n s /2 in the triangular system TDBG or the honeycomb system TBG at charge neutrality n = 0.
In particular, we focused on the study of Hund'sinduced as well as longer-ranged exchange couplings and their impact on the spin-valley liquid which has been found to emerge in the limit of SU(4) symmetry in both lattice geometries. We find extended parameter regimes where this phase is stabilized, with no signatures of long-range order, providing evidence for a stable realm of spin-valley liquid behavior. Experimentally, such a phase would be consistent with a correlated insulator lacking spin and valley polarisation. However, the precise nature of the phase and potential experimental fingerprints are left for future study, though we note that a recent projective-symmetry-group classification of fermionic partons on the half-filled triangular lattice suggests the possibility of a U(1) spin liquid with four Fermi surfaces [59], which would be consistent with our analysis. Our findings hint at the possibility of spin-valley entangled quantum liquids lurking within the correlated insulating regimes of moiré heterostructures. A. Laeuchli, and A. Paramekanti for discussions. C.H. thanks the Aspen Center for Theoretical Physics, where parts of this paper were written, for hospitality and support under Grant No. NSF 1066293. We acknowledge partial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Projektnummer 277146847 -SFB 1238 (project C03) and Projektnummer 277101999 -TRR 183 (projects B01 and A02). The numerical simulations were performed on the CHEOPS cluster at RRZK Cologne and the JURECA Booster at the Forschungszentrum Juelich.

I. Hexagonal moiré structures
As noted in the main text, the minimal model that covers the necessary universal aspects of the various moiré heterostructures is a two-orbital extended Hubbard model. With four flavors of fermions per site, two spin and two valley degrees of freedom, this leads to a four band model on the triangular lattice and an eight band model on the honeycomb lattice (where the doubling is simply due to the doubling of the unit cell). Which of these lattices is appropriate to use depends on the particular moiré heterostructure one is interested in.
For TBG, TLG/h-BN and TDBG there are a total of eight minibands near charge neutrality, four above and four below, that are separated from the rest of the spectrum by trivial band gaps. Filling of these minibands is thus typically denoted as ranging from n = −n s to n = +n s , as indicated in Fig. 5 (where, for convenience, we plot n/(n s /4)). In the case of TBG, the bands above/below charge neutrality are connected via Dirac points, meaning that any effective Hubbard model must describe all eight bands. This naturally motivates the use of the honeycomb lattice Hubbard model. Half-filling, i.e. the scenario focused on in the main text, thus corresponds to charge neutrality n = 0. On the other hand, in the case of TDBG and TLG/h-BN the bands above/below charge neutrality are disconnected from one another, meaning that an effective Hubbard model description need only focus on one or the other set of four bands. This naturally leads to a triangular lattice description, with half-filling now corresponding to n = ±n s /2.

II. Pseudofermion functional RG approach
The pseudofermion functional renormalization group (pf-FRG) has recently been established as a versatile tool for the investigation of ground state phase diagrams for a wide class of spin models [44,45]. In doing so, the free fermion propagator G 0 = (iω) −1 of a pseudofermion decomposed quartic Hamiltonian, e.g., Eq. (1), is modified by a step-like regularization function Θ(|ω| − Λ) with frequency cutoff scale Λ, i.e. G 0 → G Λ 0 = G 0 Θ Λ . The artificial scale dependence of this theory results in a hierarchy of coupled one-loop RG flow equations for the one-particleirreducible (1PI) interaction vertices. We employ a standard approximation scheme, where the hierarchy is truncated to exclusively account for the frequency-dependent self-energy Σ Λ and two-particle interaction vertex Γ Λ , see, e.g., Ref. [47] for more details and technicalities.
Here, we describe the aspects of the pf-FRG which are particular to the present spin-valley model, i.e. the vertex parametrization for the SU(2) ⊗ SU(2) symmetry and the implementation of the filling constraint.

II. A. Vertex parametrization for SU(2)⊗SU(2) symmetry
The pseudofermion decomposition of the spin-valley operatorsσ a i ,τ b i ,σ a i ⊗τ b i in the Hamiltonian, Eq. (1), exhibits a local U(1) symmetry. Consequentially, the self-energy Σ Λ can be efficiently parametrized as being local and the two-particle interaction vertex Γ Λ as being bilocal. Translation invariance in imaginary time additionally reduces the number of independent frequency arguments by one. The spin/valley dependence of the 1PI irreducible vertices can be II. C. Finite-size analysis of the RG flow An instability in the vertex function during the RG flow indicates a spontaneous breaking of symmetries that have been implemented in the initial conditions [47]. Most prominently, magnetic instabilities appear as pronounced kinks or cusps in the flow of the momentum resolved two-spin correlations. Alternatively, one may check the behavior of an on-site correlation function, i.e., χ Λ ii , for different values of the vertex range L. Formally, L does not determine the system size (which is in fact infinite in pf-FRG) but rather sets the scale on which spins can be correlated. It is then natural to expect sensitivity to this parameter near the critical scale since the physics is governed by the collective behavior of all spins. On the other hand, if the system does not develop an instability down to the smallest energy scales, i.e. the pf-FRG flow stays regular, real space correlations should be robust with respect to variations of L.
Indeed as shown in Figs. 6 and 7, flows of the spin correlation in the dominant interaction channel for different L are aligned within the paramagnetic regions of the spin-valley phase diagrams, but deviate from each other around the critical scale in the ordered phases. We find, however, that this effect is more subtle for the triangular than the honeycomb lattice, which we attribute to the inherent geometric frustration of the former. The spin-valley entangled liquid ground states of the nearest-neighbor SU(4) Heisenberg models (on both the triangular and honeycomb lattice) remain stable upon inclusion of moderate longer-ranged exchange interactions as illustrated in Fig. 3 of the main text.
Here, we provide further information about the evolution of the structure factors upon varying J 2 /J 1 . First, we recall that for J 2 /J 1 = 0, local correlations are reminiscent of 120 • (Néel) order for the triangular (honeycomb) model. Going to large antiferromagnetic J 2 > 0, stripe (spiral) order emerges with the evolution of the structure factor being plotted in Figs. 8 and 9 at the onset of these orders. Around J 2 /J 1 ≈ 0.2 for the triangular and at about J 2 /J 1 ≈ 0.3 for the honeycomb lattice, the topology of the momentum resolved correlation functions changes visibly, indicating a Lifshitz transition.