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

TL;DR

This paper investigates the conditions under which spin-valley quantum liquids emerge and remain stable in moiré heterostructures, using a theoretical model and renormalization group analysis to identify regimes of spin liquid behavior.

Contribution

It introduces a comprehensive analysis of spin liquid stability in SU(2)$^{ ext{spin}} imes$SU(2)$^{ ext{valley}}$ Heisenberg models with extended couplings in moiré heterostructures, highlighting the role of Hund's coupling and lattice geometry.

Findings

Existence of an extended parameter regime with no long-range order indicating stable spin liquids.

Identification of magnetic orders for large Hund's coupling, including antiferromagnetic and ferromagnetic states.

Longer-ranged exchange couplings can both stabilize and destabilize the spin liquid regime depending on their sign.

Abstract

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.

Equations

$$\mathcal{H}$$

$${\chi }_{ij}^{s\Lambda }=⟨{\hat{\sigma }}_i^a{\hat{\sigma }}_j^a{⟩}^{\Lambda },and/or{\chi }_{ij}^{v\Lambda }=⟨{\hat{\tau }}_i^b{\hat{\tau }}_j^b{⟩}^{\Lambda }.$$

$$J_v^{\text{eff}}=J_v+J\cdot {\chi }_{ij}^{s{\Lambda }_c},\mathrm{and}{J}_s^{\text{eff}}=J_s+J\cdot {\chi }_{ij}^{v{\Lambda }_c}.$$

$${\Sigma }^{\Lambda }(1^{\prime };1)$$

$${\Gamma }^{\Lambda }(1^{\prime },2^{\prime };1,2)$$

$${\Sigma }^{\Lambda }(1^{\prime };1)$$

$${\Gamma }^{\Lambda }(1^{\prime },2^{\prime };1,2)$$

$$\displaystyle\phantom{\bigg{[}}+\Gamma^{\Lambda ds}_{i_{1}i_{2}}(w_{1^{\prime}}w_{2^{\prime}};w_{1}w_{2})\delta_{s_{1^{\prime}}s_{1}}\delta_{s_{2^{\prime}}s_{2}}\theta^{b}_{l_{1^{\prime}}l_{1}}\theta^{b}_{l_{2^{\prime}}l_{2}}+\Gamma^{\Lambda dd}_{i_{1}i_{2}}(w_{1^{\prime}}w_{2^{\prime}};w_{1}w_{2})\delta_{s_{1^{\prime}}s_{1}}\delta_{s_{2^{\prime}}s_{2}}\delta_{l_{1^{\prime}}l_{1}}\delta_{l_{2^{\prime}}l_{2}}\bigg{]}$$

$$\displaystyle\phantom{\bigg{[}}\times\delta_{i_{1^{\prime}}i_{1}}\delta_{i_{2^{\prime}}i_{2}}\delta(w_{1^{\prime}}+w_{2^{\prime}}-w_{1}-w_{2})-(1\leftrightarrow 2)\,,$$

$${\Gamma }_{i_1{i}_2}^{\infty ss}(w_{1^{\prime }}{w}_{2^{\prime }};w_1{w}_2)=J_{},{\Gamma }_{i_1{i}_2}^{\infty sd}(w_{1^{\prime }}{w}_{2^{\prime }};w_1{w}_2)=J_s,$$

$${\Gamma }_{i_1{i}_2}^{\infty ds}(w_{1^{\prime }}{w}_{2^{\prime }};w_1{w}_2)=J_v,{\Gamma }_{i_1{i}_2}^{\infty dd}(w_{1^{\prime }}{w}_{2^{\prime }};w_1{w}_2)=0,$$

Full text

Spurred by the discovery of a plethora of insulating and superconducting states in twisted bilayer graphene (TBG) Cao et al. (2018a, b), a growing stream of experimental evidence points to the generic emergence of correlated electronic behavior in various moiré heterostructures Yankowitz et al. (2019); Sharpe et al. ; Lu et al. ; Chen et al. (2019, ); Codecido et al. ; Liu et al. ; Cao et al. ; Shen et al. ; Serlin et al. . The basic mechanism that gives rise to strongly-enhanced correlation effects in these materials is the formation of long-wavelength moiré patterns with (almost) flat low-energy bands whose narrow bandwidth becomes comparable to the otherwise negligible energy scale of the electronic interactions Trambly de Laissardière et al. (2010); Bistritzer and MacDonald (2011); Lopes dos Santos et al. (2012). Due to a high degree of control, e.g., in the regulation of the twist angle, tunable bandwidths or fillings, and a low level of disorder, such systems are discussed as ideal platforms for detailed studies of quantum many-body states. Despite a vast amount of concomitant theoretical activity Xu and Balents (2018); Yuan and Fu (2018); Dodaro et al. (2018); Kennes et al. (2018); Tang et al. (2019); Isobe et al. (2018); Lin and Nandkishore (2018, ); Rademaker and Mellado (2018); Huang et al. (2019); Fidrysiak et al. (2018); Guo et al. (2018); Roy and Juričić (2019); Liu et al. (2018); Da Liao et al. ; You and Vishwanath (2019); Sherkunov and Betouras (2018); Kang and Vafek (2019); Classen et al. (2019); Seo et al. (2019), the precise nature of the observed insulators and superconductors, however, remains to be explored and settled through the construction of faithful models and application of appropriate quantum many-body approaches.

Several model constructions for correlated moiré materials have been put forward in terms of effective tight-binding descriptions on the moiré superlattice, augmented by various interaction terms Koshino et al. (2018); Kang and Vafek (2018); Po et al. (2018). 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 Koshino et al. (2018); Kang and Vafek (2018); Po et al. (2018), related structures such as twisted double-bilayer graphene (TDBG) or trilayer graphene/hexagonal boron nitride heterostructures (TLG/h-BN) are better captured by a triangular superlattice Chen et al. (2019); Po et al. (2018); Xian et al. . 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 Xu and Balents (2018); Yuan and Fu (2018) serving as a paradigmatic starting point. Its kinetic term $H_{t}\!=\!-t\sum_{\langle ij\rangle}\sum_{\alpha=1}^{4}(c_{i\alpha}^{\dagger}c_{j\alpha}^{\phantom{\dagger}}+\mathrm{h.c.})$ for the electrons combines the spin projection $s\in\{\uparrow,\downarrow\}$ and valley quantum number $l\in\{+,-\}$ in a flavor index $\alpha\in\{(\uparrow,\!+),(\uparrow,\!-),(\downarrow,\!+),(\downarrow,\!-)\}$, 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 SM .

The simplest conceivable interaction term, which also retains the SU(4) symmetry, is a Hubbard interaction $H_{\mathrm{int}}=U\sum_{i}(\sum_{\alpha=1}^{4}n_{i\alpha})^{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 strong-coupling limit, the kinetic term can then be treated perturbatively Xu and Balents (2018); Yuan and Fu (2018); Po et al. (2018). 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}\propto 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 Xu and Balents (2018); Yuan and Fu (2018). Moreover, Wannier state constructions suggest that further-neighbor interactions can become sizable Koshino et al. (2018) 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 $\mathrm{SU(2)}^{\mathrm{spin}}\otimes\mathrm{SU(2)}^{\mathrm{valley}}$ symmetry for both triangular and honeycomb lattice geometries, which we later supplement with further-neighbor interaction terms. Our focus is on the case of half-filling of the underlying Hubbard model, i.e. two electrons per site. For the effective Heisenberg model at strong coupling, this implies that we are working with the six-dimensional self-conjugate representation of 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 Liu et al. ; Cao et al. ; Shen et al. . To explore the effect of longer-range interactions, we augment our model by a next-nearest 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 Corboz et al. (2012), while on the triangular lattice extended parameter regimes without any net magnetization have been identified in DMRG simulations Wu et al. .

Spin-valley model.

The starting point of our study is an SU(4) spin-valley Heisenberg model Xu and Balents (2018); Po et al. (2018), $\mathcal{H}_{\mathrm{SU(4)}}=J_{H}\sum_{\langle ij\rangle}\hat{T}^{\mu}_{i}\hat{T}^{\mu}_{j}$, where $J_{H}$ is the antiferromagnetic exchange coupling between nearest neighbors on either the triangular or honeycomb lattice, and $\hat{T}_{i}$ denote SU(4) spins. The $\mu\!=\!1,\dots,15$ components of the spin operators can be represented on a fermionic Hilbert space via the parton construction $\hat{T}^{\mu}_{i}=f^{\dagger}_{i\alpha}T^{\mu}_{\alpha\beta}f^{\phantom{\dagger}}_{i\beta}$, where the index $\alpha$ enumerates four different fermion flavors and the matrices $T^{\mu}$ are the SU(4) generators Arovas and Auerbach (1988). At half-filling of the underlying Hubbard model, the local spin-valley Hilbert space is six-dimensional (4 choose 2), which leads to a local filling constraint of two partons per lattice site $\sum_{\alpha}f^{\dagger}_{i\alpha}f^{\phantom{\dagger}}_{i\alpha}\!=2$.

Upon inclusion of Hund’s couplings, the SU(4) symmetry of the model is explicitly broken Xu and Balents (2018). Omitting other sources of SU(4) breaking, a residual separate spin-valley $\mathrm{SU(2)}^{\rm s}\otimes\mathrm{SU(2)}^{\rm v}$ symmetry remains which is reflected by the extended Hamiltonian

[TABLE]

where the spin-valley operators read $\hat{\sigma}^{a}_{i}=f^{\dagger}_{is^{\prime}l^{\prime}}\theta^{a}_{s^{\prime}s}\delta_{l^{\prime}l}f^{\phantom{\dagger}}_{isl}$, $\hat{\tau}^{b}_{i}=f^{\dagger}_{is^{\prime}l^{\prime}}\delta_{s^{\prime}s}\theta^{b}_{l^{\prime}l}f^{\phantom{\dagger}}_{isl}$, and $\hat{\sigma}^{a}_{i}\otimes\hat{\tau}^{b}_{i}=$ $f^{\dagger}_{is^{\prime}l^{\prime}}\theta^{a}_{s^{\prime}s}\theta^{b}_{l^{\prime}l}f^{\phantom{\dagger}}_{isl}$. Instead of enumerating the four fermion types by a single index, we have exposed the spin quantum number $s\in\{\uparrow,\downarrow\}$ and the valley quantum number $l\in\{+,-\}$ explicitly; Pauli matrices are denoted by $\theta^{a}$, $a\in\{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 Wu et al. , 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) Reuther and Wölfle (2010); Wetterich (1993). 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$ Baez and Reuther (2017) (which naturally favors magnetic order) and in the SU($N$) spin symmetry Buessen et al. (2018); Roscher et al. (2018) (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 spin-valley model at hand. We extend the standard implementation of pf-FRG to incorporate the $\mathrm{SU(2)}^{\rm s}\otimes\mathrm{SU(2)}^{\rm v}$ symmetry, thereby obtaining flow equations for the one-particle irreducible vertices as a function of an RG frequency cutoff scale $\Lambda$ SM . Numerically solving the set of $\mathcal{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 Reuther and Wölfle (2010); Roscher et al. which occurs at some critical scale $\Lambda_{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 ($\omega=0$) of the correlation functions in the spin sector and in the valley sector,

[TABLE]

Sharp features emerging in the respective structure factors $\chi^{s/v}(\vec{q})\propto\sum_{ij}e^{i\vec{q}\cdot(\vec{r_{i}}-\vec{r_{j}})}\chi_{ij}^{s/v}$ 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 SM for details on the finite-size analysis). This rules out not just magnetically ordered states, but also valence bond or dimer crystals 111Generally, an instability can be observed for any (quantum) order which is captured by a fermionic bilinear operator, see Ref. Buessen et al. (2018). More complicated order whose order parameter involves higher order operators does not necessarily manifest in a flow breakdown, an ordering which spins in the self-conjugate representation are often prone to Marston and Affleck (1989); Rokhsar (1990). For SU(4) spins in the self-conjugate representation we can further use the Lieb-Schultz-Mattis-Hastings Lieb et al. (1961); Hastings (2004); Yao et al. 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^{\circ}$/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\neq 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 $\Lambda_{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$,

[TABLE]

Note that for $120^{\circ}$ or Néel order in one of the SU(2) sectors the corresponding nearest-neighbor correlation becomes negative. Therefore, the effective couplings $J_{v}^{\text{eff}}$ and $J_{s}^{\text{eff}}$ may, too, turn negative and drive a ferromagnetic instability in the other sector, despite the antiferromagnetic nature of all couplings in the microscopic spin-valley model Wu et al. . This kind of mechanism may be at the origin of the spin polarization observed at half-filling in TDBG Liu et al. ; Cao et al. , as first pointed out in Ref. Wu et al. 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 222 We finally note that, in principle, the FRG setup can also be extended to allow for a continuation of the flow into the ordered regime Salmhofer et al. (2004); Roscher et al. (2018, ), however, such modifications are beyond the scope of this work. . 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é heterostructures it has been pointed out that longer-ranged Coulomb interactions should not be neglected Koshino et al. (2018), 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-nearest-neighbor 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 Zhu and White (2015); Hu et al. (2015), 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}\in[-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 wedge-shaped 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 $\mathrm{SU(2)}^{\rm s}\otimes\mathrm{SU(2)}^{\rm 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 SM , e.g. the insulator at half-filling $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’s-induced 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 Zhang and Mao , 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.

Acknowledgments.

We thank L. Classen, P. Corboz, 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.

Appendix A 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=\pm n_{s}/2$.

Appendix B 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 Reuther and Wölfle (2010); Wetterich (1993). In doing so, the free fermion propagator $G_{0}=(i\omega)^{-1}$ of a pseudofermion decomposed quartic Hamiltonian, e.g., Eq. (1), is modified by a step-like regularization function $\Theta(|\omega|-\Lambda)$ with frequency cutoff scale $\Lambda$, i.e. $G_{0}\to G_{0}^{\Lambda}=G_{0}\Theta^{\Lambda}$. The artificial scale dependence of this theory results in a hierarchy of coupled one-loop RG flow equations for the one-particle-irreducible (1PI) interaction vertices. We employ a standard approximation scheme, where the hierarchy is truncated to exclusively account for the frequency-dependent self-energy $\Sigma^{\Lambda}$ and two-particle interaction vertex $\Gamma^{\Lambda}$, see, e.g., Ref. Buessen et al. (2018) 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 $\mathrm{SU(2)}\otimes\mathrm{SU(2)}$ symmetry and the implementation of the filling constraint.

B.0.1 II. A. Vertex parametrization for SU(2)$\otimes$SU(2) symmetry

The pseudofermion decomposition of the spin-valley operators $\hat{\sigma}^{a}_{i},\hat{\tau}^{b}_{i},\hat{\sigma}^{a}_{i}\otimes\hat{\tau}^{b}_{i}$ in the Hamiltonian, Eq. (1), exhibits a local U(1) symmetry. Consequentially, the self-energy $\Sigma^{\Lambda}$ can be efficiently parametrized as being local and the two-particle interaction vertex $\Gamma^{\Lambda}$ 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 expanded in terms of an $\mathfrak{su}(2)$ basis. This scheme is augmented by symmetry-allowed SU(2)-invariant density terms with the most general parametrization reading

[TABLE]

where $1=\{i_{1},s_{1},l_{1},w_{1}\}$ and $\alpha,\beta\in\{0,1,2,3\}$ with $\theta^{0}=\mathbf{1}$. Exploiting SU(2) symmetry in both spin and valley indices we are left with pure density contributions for the self-energy, while the two-particle vertex may also contain off-diagonal terms albeit with equal spin directions, i.e.

[TABLE]

where $a,b\in\{1,2,3\}$.

The initial conditions at the UV scale then read $\Sigma^{\infty}_{i_{1}}(w)=0$ for the self-energy and

[TABLE]

for the two-particle interaction vertices.

B.0.2 II. B. Particle-hole symmetry and the half-filling constraint

The pf-FRG approach considered here explicitly implements the particle-hole symmetry of the spin-valley model on the level of the 1PI vertices. As a consequence, the half-filling constraint is naturally fulfilled on average upon setting the chemical potential to zero and local fluctuations around the average are energetically penalized Reuther and Wölfle (2010). Therefore, contributions from unphysical sectors of the Hilbert space violating the local half-filling constraint are strongly suppressed at zero temperature. Consequently, the local half-filling constraint is expected to be well-enforced within the pf-FRG approach as long as the flow does not develop an instability.

Specifically, in the model studied here, the local Hilbert space for fermionic flavors $\alpha\in\{(\uparrow\!+),(\uparrow\!-),(\downarrow\!+),(\downarrow\!-)\}$ is equipped with the particle-number basis $B=\{|n_{1},...,n_{4}\rangle\}$. We define the linear unitary operator $\mathcal{P}$ acting on the basis by exchanging each occupied state with an empty state $\mathcal{P}|n_{1},...,n_{4}\rangle=|1-n_{1},...,1-n_{4}\rangle$. By computing the corresponding matrix elements, one finds that $\mathcal{P}$ transforms creation and annihilation operators into each other, i.e $\mathcal{P}^{\dagger}f^{\dagger}_{\alpha}\mathcal{P}^{\phantom{\dagger}}\!=\!f^{\phantom{\dagger}}_{\alpha}\,,\ \mathcal{P}^{\dagger}f^{\phantom{\dagger}}_{\alpha}\mathcal{P}^{\phantom{\dagger}}\!=\!f^{\dagger}_{\alpha}$. This transformation leaves the spin-valley Hamiltonian invariant. On the level of vertex functions, we obtain the identities $\Sigma^{\Lambda}(1^{\prime};1)=-\Sigma^{\Lambda}(-1;-1^{\prime})$ and $\Gamma^{\Lambda}(1^{\prime},2^{\prime};1,2)=\Gamma^{\Lambda}(-1,-2;-1^{\prime},-2^{\prime})$, where the minus sign applies only to Matsubara frequencies. The vertex components therefore obey $\Sigma^{\Lambda}_{i_{1}}(w)=-\Sigma^{\Lambda}_{i_{1}}(-w)$ and $\Gamma^{\Lambda\zeta}_{i_{1}i_{2}}(s,t,u)=\Gamma^{\Lambda\zeta}_{i_{1}i_{2}}(-s,-t,-u)$ where $\zeta\in\{ss,sd,ds,dd\}$.

These symmetries are explicitly implemented in our pf-FRG approach. However, only one local subspace, namely the one with two occupied states, is mapped to itself by $\mathcal{P}$, i.e. by enforcing the symmetries of that respective subspace we can be sure that the RG flow starts from half-filling at each lattice site. Furthermore, since the particle number per site must be conserved as a consequence of local U(1) symmetry, hopping processes that alter the filling would trigger a measurable non-magnetic instability of the flow, which we do not observe here.

B.1 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 Buessen et al. (2018). 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., $\chi^{\Lambda}_{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.

B.2 III. Structure factor evolution in the spin-valley liquid of the $J_{1}$-$J_{2}$ model

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^{\circ}$ (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}\approx 0.2$ for the triangular and at about $J_{2}/J_{1}\approx 0.3$ for the honeycomb lattice, the topology of the momentum resolved correlation functions changes visibly, indicating a Lifshitz transition.