Moments and multiplets in moiré materials

The observation of strongly correlated states in moiré systems has renewed the conceptual interest in magnetic systems with higher SU(4) spin symmetry, e.g., to describe Mott insulators where the local moments are coupled spin–valley degrees of freedom. Here, we discuss a numerical renormalization group scheme to explore the formation of spin–valley ordered and unconventional spin–valley liquid states at zero temperature based on a pseudo-fermion representation. Our generalization of the conventional pseudo-fermion functional renormalization group approach for su\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathfrak {s}}}{{\mathfrak {u}}}$$\end{document}(2) spins is capable of treating diagonal and off-diagonal couplings of generic spin–valley exchange Hamiltonians in the self-conjugate representation of the su\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathfrak {s}}}{{\mathfrak {u}}}$$\end{document}(4) algebra. To achieve proper numerical efficiency, we derive a number of symmetry constraints on the flow equations that significantly limit the number of ordinary differential equations to be solved. As an example system, we investigate a diagonal SU(2)spin\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\text {spin}}$$\end{document}⊗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\otimes $$\end{document} U(1)valley\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\text {valley}}$$\end{document} model on the triangular lattice which exhibits a rich phase diagram of spin and valley ordered phases.


Introduction
Moiré materials that exhibit flat bands such as twisted bilayer graphene (tBG) or certain van der Waals heterostructures such as trilayer graphene on hexagonal boron nitride (TLG/h-BN) have recently been established as novel, highly tunable platforms for the study of strongly correlated electrons. Relative to an almost vanishing bandwidth, residual interactions in these materials can induce a plethora of different many-body phenomena ranging from the formation of correlated insulators [1][2][3][4] and superconductors [5][6][7] to anomalous quantum Hall effects [8]. However, a microsopic description of these phenomena is a formidable challenge as the number of of low-energy degrees of freedom is often increased [9][10][11] in comparison to conventional Mott insulators.
More specifically, it has been argued [12,13] that multi-orbital Hubbard models can describe the flat band physics in, e.g., TLG/h-BN within the topologically trivial regime, where fully symmetric Wannier states may be constructed [14]. The proposed interaction terms for the corresponding Hamiltonians usually include a generalized Hubbard U [12,13,15] as well as Hund's type couplings. Performing a strong coupling expansion where one treats the interactions as the dominant energy scale, these extended Hubbard mod- 1 With su(4), we refer to the Lie algebra of the Lie group SU (4). a e-mail: gresista@thp.uni-koeln.de (corresponding author) els can then be mapped to su(4) 1 spin-valley Hamiltonians that may be used as a starting point to investigate the nature of the correlated insulating states. The so-derived su(4) models bear a close resemblance to Kugel-Khomskii models [16] that have a long history in the study of transition metal oxides, where they are used to capture the Jahn-Teller physics of intertwined spin and orbital degrees of freedom. Increasing the number of relevant microscopic degrees of freedom (in comparison to conventional quantum magnets) has been particularly appreciated to boost quantum fluctuations independent of, e.g., lattice geometries [17], which has made Kugel-Khomskii models a recurring target in the search for unusual many-body states such as quantum spin-orbital liquids [18][19][20][21]. As such, one might expect the su(4) spin-valley physics relevant to the correlated insulating states of moiré materials to hold similar promise for the observation of spin-valley liquid states with macroscopic entanglement and potentially long-range, topological order.
In this manuscript, we present a powerful numerical scheme to analyze such su(4) spin-valley (or spinorbital) models based on a functional renormalization group (FRG) technique. Our approach is based on the pseudo-fermion FRG (pf-FRG) [22], approximating the elementary spin operators of the six-dimensional, selfconjugate representation of su(4) by auxiliary complex fermions combined with an on-average constraint on the number of particles per site. Our approach allows to go beyond mean-field level by treating competing instabilities in different interaction channels on equal footing, and is able to capture both, long-ranged spin and/or valley ordered states as well as spin-valley liquid phases. In expanding previous work (by some of us) [21], we extend the range of applicability of this approach to models with off-diagonal interactions in either spin or valley space by formulating an efficient vertex parametrization derived from a meticulous symmetry analysis. We demonstrate the feasibility of this method by studying a spin-valley Hamiltonian with SU(2) spin ⊗ U(1) valley symmetry where we identify a plethora of spin and valley orderded phases from a state-of-the-art numerical implementation of pf-FRG [23,24].
The remainder of this manuscript is structured as follows. To begin with, we introduce the spin-valley Hamiltonian of interest on a general level and discuss its specific form for TLG/h-BN as a concrete example in Sect. 2. We will continue by reviewing the pf-FRG approach (Sect. 3), its generalization for su(4) models as well as the implementation of model specific symmetries (Sect. 4). Finally, numerical results for the phase diagram of a SU(2) spin ⊗ U(1) valley model on the triangular lattice are presented and examined in Sect. 5.

Model
Microscopically, the SU(4) models of interest in this manuscript can be cast in terms of a general Hamiltonian which couples two elementary su(2) degrees of freedom, captured by the operators σ and τ , which might denote a spin and valley (or oribtal) degree of free-dom. The overall SU(4) symmetry of the Hamiltonian arises from the balanced couplings of equal strength in both degrees of freedom, i.e., J is identical for the Heisenberg-like coupling of spins σ i σ j on sites i and j (with σ i = (σ x j , σ y j , σ z j ) T ) and a similar interaction of the valley degrees of freedom τ i τ j . Such valley degrees of freedom arise, in the context of tBG and related moiré materials, from the Dirac cones in the original graphene bands, which hybridize between the two layers upon twisting and thereby add an extra index [25] to the moiré bands, as illustrated in Fig. 1. Before drawing broad attention in the context of moiré materials, the spin-orbital variant of this model has been widely studied as Kugel-Khomskii model [16], often in connection with Jahn-Teller physics in transition metal oxides where spin and orbital ordering are intertwined [26]. We note that while we will frame our discussion of the SU(4) model (1) in the language of spin-valley physics relevant to moiré materials, the presented pf-FRG approach is equally applicable in the study of such spin-orbital models. We will return to this point in the discussion section at the end.
In what we will discuss in the following, we will put a focus on the self-conjugate representation of su (4), where the spin-valley operators can be represented in terms of fermionic creation and annihilation operators as with a local half-filling constraint . 1 a Moiré pattern emerging in two stacked layers of graphene with a relative twist angle θ. Clearly visible are the different regions with AA, AB, and BA stacking leading to a triangular super-lattice structure. b Construction of the two degenerate mini Brillouin zones from the difference of the K (or K') points of the two layers of graphene. In addition to the spin degree of freedom, indicated by the grey arrows, the electrons obtain a valley degree of freedom due to the possibility of being in either one of the mini Brillouin zones at the two valleys (at the K and K' points) of the graphene band structure subject to every lattice site, where summation over repeated spin indices s and valley indices l is implied.
Here, θ μ denotes a Pauli matrix with μ ∈ {0, 1, 2, 3} and θ 0 = 1. Allowing also for more generic, i.e., SU(4) breaking, interactions, any bilinear spin-valley Hamiltonian can be written as where J μν s,ij ⊗ J κλ v,ij is understood as the Kronecker product of the spin and valley exchange matrices and summation over repeating μ, ν, κ or λ is again implied. Here,n i is the density operatorn i ≡ σ 0 , and the term proportional to the coupling I ij is needed to potentially cancel the density term ∼ σ 0 i τ 0 i J 00 s,ij J 00 v,ij σ 0 j τ 0 j , which does not appear in pure su(4) spin models as, e.g., the SU(4) symmetric Hamiltonian in Eq. (1).
To keep the numerical effort for employing our pf-FRG approach at a manageable level, we assume a specific form of the exchange matrices, namely, that both, the spin and the valley exchange only couple bilinears of spin/valley or density operators and that the spin exchange is Z 2 × Z 2 × Z 2 symmetric, and thus This form, although it spoils the generality of Eq. (4) it is nevertheless relevant to certain practical applications. For instance, the effective Hamiltonian for TLG/h-BN [11] can be recast to this form. Originally, the former is often given as which, in addition to SU(4) symmetric nearest-neigh bour (∼ J 1 ) and next-nearest-neighbour (∼ J 2 ) interactions, contains both diagonal ∼ J 1 p,ij and off-diagonal ∼ J 2 p,ij valley exchange that breaks the SU(4) symmetry down to an SU(2) spin ⊗ U(1) valley symmetry. Comparing this model to the form of the general spin-valley Hamiltonian defined in Eq. (4), the nearest-neighbour exchange matrices can be written as and the next-nearest-neighbour exchange is fully SU(4) symmetric, showing that they are indeed captured by the exchange matrices defined in Eq. (5).

pf-FRG for spin-valley models: an overview
We now proceed to the core methodological advancement of this manuscript, which will be laid out in this section-the extension of the conventional pf-FRG to spin-valley models described by Hamiltonians of the form given in Eq. (4), with general, diagonal, and offdiagonal couplings as defined by Eq. (5). To set the stage, we will first revisit the flow equations of the conventional pf-FRG approach for su(2) spins and explain how the numerical solution of the flow equations can be used to determine whether and what type of magnetic order forms for a particular spin Hamiltonian at zero temperatures. We then proceed to the adapted pf-FRG approach for spin-valley models, for which we derive an efficient parametrization of the self-energy and two-particle vertex in what is a direct extension of the parametrization for su(2) spin models with generic two-spin interactions [27]. Our particular focus is on constraints that symmetries of the spin-valley Hamiltonian pose on the parametrized vertex functions-very similar to the su(2) case but with slight differences which we especially highlight. To put these equations into numerical practice, we discuss our implementation of the spin-valley pf-FRG approach and its algorithmic scaling. This section is intended as an overview stating the main results of our study important for the implementation of the pf-FRG for spin-valley models. Readers looking for a more detailed discussion of how the symmetries of the Hamiltonian lead to the parametrization and symmetry constraints are referred to Sect. 4.

Pseudo-fermion functional renormalization group
Let us set the stage by revisiting some of the conceptual steps of the pseudo-fermion FRG, which has orig-inally been formulated for bilinear su(2) spin models [22] with generic (diagonal and off-diagonal) interactions [27] and later generalized to SU(N) Heisenberg models [28], in the context of the spin-valley models at hand. By going to a pseudo-fermion representation of the original degrees of freedom, one arrives at a fermionic representation of the original model (with an additional half-filling constraint) as outlined in the previous section. One can then proceed to apply the well-established methods of the fermionic FRG [29,30]. An important distinction to electronic systems is that the pseudo-fermion Hamiltonian exhibits only a quartic interaction term and no quadratic kinetic terms. This readily implies that the free propagator is diagonal in all its arguments and takes the simple form consists of a lattice site index i 1 , a spin index s 1 , a fermionic Matsubara frequency ω 1 and, for spin-valley models, the additional valley index l 1 . To implement the RG scale, or cut-off, Λ, we multiply a regulator to the free propagator where we choose a smooth regulator for improved numerical stability. The pf-FRG flow equations are then given as a special case of the general fermionic FRG equations by assuming that the flowing self-energy is, just as the free propagator, diagonal in all its arguments. This assumption is true for arbitrary spin-model bilinear in su(2) spin operators [27]. For spin-valley Hamiltonians, however, we will show in Sect. 4 that this is only the case if the couplings are diagonal in either the spin or valley sector. That is why, in this work, we always consider couplings diagonal in the spin sector as stated in Eq. (5). In the context of moiré materials, most physically relevant spin-valley models are indeed of this form. This additional assumption, therefore, leaves our method still generally applicable to most models of interest.
In the original implementation of the pf-FRG [22] and most works since, then the flow equations are truncated using the Katanin truncation scheme [31], which we also adapt here 2 . In the Katanin truncation, only the selfenergy Σ Λ and the two-particle vertex Γ Λ are considered, while higher order vertex functions are neglected.
The flow equations are then given by for the self-energy and for the two-particle vertex. Here, the single-scale propagator is defined as . Note that the flow equations are formulated in the T → 0 limit and the sums should therefore be understood as The fermionic representation of the spin-valley operators, as presented in the previous section, necessitates the enforcement of a local half-filling constraint (f † isl f † isl = 2) to determine the dimensionality of the local Hilbert space. To this end, we employ the same technique used for su(2) models, where the constraint is fulfilled only on average by explicitly implementing particle-hole symmetry on the level of the vertex functions [22][23][24], as will be discussed in detail in Sect. 4. Numerically, computing the expectation value f † isl f † isl from the self-energy and two-particle vertex, we confirm that the average constraint is indeed fulfilled during the entire pf-FRG flow. Although particle-number fluctuations violating the exact constraint have been shown to be sizeable, recent studies suggest that they leave the physical results obtained from the pf-FRG qualitatively unaffected [21,23,34]. We note that in contrast to su(2) spin models, the physically relevant filling for spinvalley models depends on the considered material and may also be at quarter-filling (f † isl f † isl = 1) [11,35], corresponding to the fundamental representation of su(4). At quarter-filling, however, the spin-valley Hamiltonian is no longer particle-hole symmetric and the constraint cannot be enforced in the same efficient manner. In this work, we therefore limit ourselves to half-filling.
To identify the ground state of a model of interest, we numerically solve the flow equations (as discussed in more detail below in Sect. 3.3) and thereby calculate the flow of various correlation functions from the flow of the vertex functions. In its most general we define a spin-valley-spin-valley correlation function where T τ is the time-ordering operator. From this general definition, we can then read off the form of spinspin correlations as well as valley-valley correlations A thermal phase transition to long-range, symmetrybreaking order in the spin or valley sector at some finite temperature can formally be detected by a divergence in the RG flow of the corresponding correlation at some breakdown scale Λ c [28], as shown in Fig. 2a.
Due to finite numerical resolution, however, they often manifest as a kink or a peak in the susceptibility. The momentum space profile of the dominant structure factor close to the breakdown scale, i.e., the Fourier transform of the static correlation χ Λs/v ij (ω = 0), then indicates the type of symmetry-breaking. Since the solution of the flow equation below the breakdown scale Λ c is no longer physical, this only allows us to detect the phase transition that occurs at the largest breakdown scale if there are multiple subsequent transitions. This might be the case when spin and valley degrees of freedom exhibit different ordering transitions at two distinct energy scales. If, in this scenario, the spin sector orders at the larger of the two energy scales, we cannot directly determine the ground-state order of the valley sector from the flow of the valley-valley correlations. Instead, we need to fall back to, for instance, meanfield arguments as proposed in [21] to determine the most likely valley order. If, on the other hand, the correlations show no flow breakdown, both spin and valley degrees of freedom do not order, indicative of a ground state that remains paramagnetic or exhibits spin-valley liquid behavior.
These two scenarios are illustrated in Fig. 2a, b. Both panels show the flow of the structure factor at the dominant momentum for a magnetically ordered phase with dominant valley order (a) and the paramagnetic state at the SU(4) point [36] (b) where the spin-valley Hamiltonian corresponds to Eq. (1). In the magnetically ordered phase of panel (a), we see a clear flow breakdown in the valley structure factor χ Λv , which manifests as a peak or divergence, depending on the vertex truncation length L (further discussed in Sect. 3.3). The spin structure factor χ Λs shown by the purple lines is strongly suppressed. At the SU(4) point, on the other hand, the flow of the structure factor is smooth and convex down to the lowest energy scale we can reliable calculate (Λ = 0.02J), indicating a paramagnetic ground state. Here, spin and valley correlations are identical due to the global SU(4) symmetry of the Hamiltonian (and indistinguishable in our plot).

Vertex parametrization and symmetry constraints
To make the solution of the flow equations numerically feasible, one needs to keep the overall number of differential equations needed to capture the flow equations as small as possible. Practically, this can be achieved by eliminating redundant calculations through implementing the symmetry constraints which the Hamiltonian poses on the self-energy and the two-particle vertex. A comprehensive symmetry analysis of this sort has been carried out for generic su(2) spin models [27], which here will be generalized to the spin-valley Hamiltonians of interest. Details of this symmetry analysis will be discussed in Sect. 4, while we will report its main findings in the following.
The first important finding is that symmetries dictate that the self-energy is completely diagonal and can be All structure factors are shown at the momentum at which they are maximal. The insets zoom into the flow at small cut-offs. In the magnetically ordered phase, we clearly see a breakdown of the flow in the valley sector, which manifests as a peak for small L and a more clear divergence when increasing L. In the paramagnetic phase, the flow is smooth and convex down to about Λ/J = 0.02, which is the smallest scale for which our calculations are numerically reliable parametrized by a single function Σ(ω) as We emphasize again that this is only the case if the interactions remain diagonal in either the spin or valley sector. For Hamiltonians with off-diagonal interactions in both sectors, the self-energy will not be diagonal in the spin and valley indices, greatly increasing the numerical cost for the solution of the flow equations. The two-particle vertex can be parametrized as with the three bosonic transfer frequencies This parametrization is of the same form as for su (2) spin models-apart from an increased number of components due to the valley sector ∼ θ κ l 1 l1 θ λ l 2 l2 with the corresponding indices κ and λ. If we assume the Hamiltonian to be diagonal in the spin sector, we will only need to consider components diagonal in the spin ∼ θ μ s 1 s1 θ μ s 2 s2 , with the corresponding index μ (and vice versa for a system with a diagonal valley Hamiltonian). The basis functions of the parametrization are constrained by the symmetries of the Hamiltonian as where we defined the sign function These are the same relations as for the su(2) case, apart from a missing constraint relating the s und u frequencies in the two-particle vertex (c.f. Eq. (14) in Ref. [27]). This is a consequence of the Hamiltonian only being invariant under a global particle-hole symmetry instead of the local particle-hole symmetry under which the su(2) Hamiltonian is invariant. We discuss this in more detail in Sect. 4. The missing relation, however, does not change the key implications of the constraints, namely that the basis functions are either completely real or imaginary, and that values of the vertex functions at negative transfer frequencies can be inferred from the positive frequency axes. The parametrization of the two-particle vertex using the three transfer frequencies in Eq. (17) is convenient for deriving the flow equations and symmetry constraints. However, to better capture the asymptotic frequency dependence of the two-particle vertex one can further refine the frequency parametrization [23,24,37]. The first step is to group the contributions in the flow equation of the two-particle vertex given in Eq. (11) into three channels according to their twoparticle irreducibility. This results in a particle-particle (pp), direct particle-hole (dph), and crossed particlehole (cph) channel, which correspond to the three contributions on the right-hand side (RHS) of Eq. (11), in the respective ordering. In these terms, the flow equation for the two-particle vertex can be written as and the vertex is parametrized (stating only the frequency dependence) as where Γ Λ→∞ is the bare two-particle vertex at infinite cut-off. Each channel g c (ω c , v c , v c ) is parametrized by one bosonic transfer frequency ω c and two additional fermionic frequencies v c , v c . The precise definition of the frequencies can be chosen in numerous ways. It is, however, advantageous to choose them, so that the symmetry constraints of the two-particle vertex given in Eq. (19) result in equally simple relations for each channel in the new parametrization. Here, we adapt the choice of Ref. [24] and give the resulting symmetry constraints for the channels in A. Compared to su(2) spin models, no constraint relating the particle-particle and crossed particle-hole channel with each other is present, which can be traced back to the missing symmetry constraint relating the s and u frequency dependence 3 .
To complete the discussion, we still need to state the initial conditions of the flow equations corresponding to the self-energy and two-particle vertex in the limit Λ → ∞, which are given by with the couplings J μ s,i1i2 and J κλ i1i2 defined in Eq. (5).

Numerical implementation
The numerical solution of the pf-FRG flow equations poses several challenges and necessitates further approximations to be made. To overcome these challenges, we employ the state-of-the-art numerical implementation of Refs. [23,24], where additional details of the implementation are discussed. Here, we only give a short overview and discuss some slight technical differences in the implementation for spin-valley models.
First, one has to truncate the infinite lattice geometry by a finite lattice graph. Employing the symmetries of the lattice geometry for which the spin-valley model is formulated and the local U(1) symmetry present in all pseudo-fermion Hamiltonians, the spatial dependence of the two-particle vertex can be reduced to just one site index j and one arbitrary fixed reference site i 0 , as will be derived in Sect. 4. To obtain a finite number of vertex components Γ Λ i0j (considering only the lattice site dependence), we define a finite length scale L and truncate the vertex Γ Λ i0,j for bond distances d(i 0 , j) > L, effectively enforcing a maximal correlation length. The finite-size effect of this truncation can be observed in Fig. 2, where several calculations with increasing values of L were performed for a magnetically ordered and a paramagnetic phase. In the ordered phase, the flow breakdown sharpens from a relatively broad peak for low values of L to a clear divergence for larger values of L, which is a typical observation. The paramagnetic phase is, in contrast, not affected by the increase of L (at least qualitatively). From an algorithmic point of view, the asymptotic scaling of the computation time is quadratic in the number of lattice points N L ∼ L d , where d is the number of spatial dimensions. This is due to the fact that the number of vertex components as well as the sum over all lattice sites included in the flow equations scale linearly with N L . In this work, we typically perform calculations at L = 9, above which the breakdown scale does not significantly change anymore and the numerical effort is still reasonable.
Since the pf-FRG approach is formulated at zero temperature, another point we need to address is how to discretize the continuous Matsubara frequencies. To accurately resolve all features of the two-particle vertex, it turns out that particular care needs to be taken in the choice of frequency meshes [23,24]. To this end, the frequencies are discretized on adaptive, hybrid linearlogarithmic meshes, which are updated using a scanning routine between each step of the ordinary dif-ferential equation (ODE) solver. In addition to continuous Matsubara frequencies, the flow equations at T = 0 include frequency integrals which have to be performed numerically. To calculate these integrals, we employ an adaptive quadrature which takes both the relevant features around the origin and the algebraic decay for large frequencies into account. Values of the vertex for frequencies not lying on the discrete frequency meshes are obtained by multi-linear interpolation. The computation time asymptotically scales with the number of (positive) bosonic frequencies N Ω and (positive) fermionic frequencies N ν as O(N Ω · N 2 ν ). A typical setup for which the two-particle vertex is sufficiently well resolved is N Ω = 40 and N ν = 30, which we use for all calculations in this work. The computational effort to compute the self-energy is, compared to the vertex, negligible, as it only depends on one frequency. Here, we choose a frequency mesh with N Σ = 250 frequencies. In the su(2) case, only positive frequencies were required, as the symmetry constraints map all negative frequency components to some positive counterpart. For spin-valley models, however, due to the missing symmetry constraint relating the particleparticle and crossed particle-hole channel (discussed in Sect. 3.2), we have to also consider negative frequencies for either ν c or ν c . This results in an additional factor of two in computation time compared to su(2) spin models.
The adaptive frequency meshes and integration routine allow for an efficient evaluation of the RHS of the flow equations. For the solution of the ODEs themselves, we choose the Bogacki-Shampine method [38], which is a third-order Runge-Kutta method with adaptive step size control. We find that this method is a good compromise between computational cost and numerical precision.
Although the asymptotic scaling of the computation time with the number of lattice points and frequencies is the same as for the su(2) case, more complex spinvalley models usually require a much larger numerical effort, as the extra valley index greatly increases the number of independent two-particle vertex components Γ Λ,μκλ i1i2 , in which the computation time scales linearly. With the coupling matrices given in Eq. (5), there would be N Γ = 4 · 4 2 = 64 independent vertex components (only considering the spin-valley dependence). In comparison, the parametrization for generic su(2) models only has N Γ = 4 2 = 16 components. Fortunately, in almost all physical models, extra symmetries in the spin and valley space will greatly reduce the number of independent components. Considering, e.g., an SU(2) symmetry in the spin space and a U(1) symmetry in valley space, which is present in several models for moiré materials [11,35], the number is already reduced to N Γ = 2 · 6 = 12. For these models, the numerical effort is similar to su(2) models with off-diagonal interactions and even allows for computations of relatively large-phase diagrams as will be presented in Sect. 5.

Symmetry classification
To proof the validity of the parametrization and the symmetry constraints presented in the previous section, we repeat the symmetry analysis of Ref. [27], where the pseudo-fermion Hamitonian for su(2) spin models with generic diagonal and off-diagonal interactions is considered, but for the spin-valley Hamiltonian given in Eq. (4). We show that most of the symmetries of the su(2) pseudo-fermion Hamiltonian are either also present in the spin-valley Hamiltonian, or can be generalized in a straightforward fashion. There are, however, some differences that we will highlight in the following. Most notably, we show that, even at the SU(4) point, the spin-valley model does not posses a local particlehole symmetry that is present in the su(2) case, but only the corresponding global symmetry. Consequently, it is also not present in generalizations of the SU(2) Heisenberg model to SU(N), which might not have been clearly stated before. This is the reason for the missing symmetry constraint for the two-particle vertex as presented in the previous section.

Local U(1) symmetry
The first symmetry transformation we consider, a local U(1) transformation, directly follows from the form of the spin-valley operator given by Eq. (2). It acts on the fermionic Hilbert space at site i by multiplying a local phase ϕ i ∈ [0, 2π) to the fermionic operators as which clearly leaves all spin-valley operators invariant.
Interpreting the spin-valley Hamiltonian as a fermionic representation of an su(4) spin model, it is simply a consequence of the choice for the fermionic representation of the spin operators. It is therefore also present in all conventional pf-FRG implementations using the standard pseudo-fermion representation. In that sense, it is sometimes also referred to as a gauge redundancy instead of a symmetry, as it is not a symmetry of the original spin Hamiltonian, but only of the pseudofermion representation. For our functional renormalization group approach, we are interested in the implication of the symmetry on the functional form 4 of the one-particle correlation function (26) 4 Note that our definition deviates from normal ordering to be in line with the conventional definition of retarded Greens functions. and the two-particle correlation function where we suppress the time-ordering operator as it becomes trivial in the path integral framework that the function renormalization group is formulated in. Acting with the local U(1) transformation given in Eq. (25) on the definition of the correlation functions and demanding their invariance leads to the corresponding symmetry constraint. It directly implies that we can restrict ourselves to a local one-particle correlation function which only depends on one lattice site i 1 , and a bi-local two-particle correlation function which only depends on the two lattices sites i 1 and i 2 .

Global particle-hole symmetry
In the pf-FRG for su(2) spin models, spin operators S a i are represented using fermions with one spin index α = ±1 as with a ∈ {1, 2, 3}. In addition to the U(1) gauge redundancy, there exists another redundancy in this representation that can be formulated as a local particle-hole symmetry [27]. It acts on the fermionic Hilbert space as withᾱ ≡ −α. It leaves the fermionic representation of the su(2) spin operators invariant and is therefore a symmetry of the pseudo-fermion Hamiltonian. We note that this symmetry is not anti-unitary and therefore does not correspond to the usual physical particle-hole symmetry [27]. Instead, it is again a consequence of the representation of the spin operators. The natural extension for spin-valley models with spin index s = ±1 and valley index l = ±1 is the transformation under which the spin-valley operator transforms as which can be shown straightforwardly using the anticommutation relations of the fermionic operators and the identityᾱᾱ  For the local single-particle correlation function, the global particlehole symmetry implies and for the bi-local two-particle correlator, it implies These relations are, apart form the extra factors of valley indices, the same as for the su(2) case when considering the global transformation. The invariance under the local transformation would yield additional constraints on the two-particle correlator acting only on multi-indices with the same lattice site (i 1 or i 2 ). For the parametrized two-particle vertex, these result in a constraint relating the s and u dependence or, in the asymptotic frequency parametrization defined in Eqs. (22, 23), the particle-particle and crossed particlehole channel with each other. As already discussed in Sect. 3, this constraint is, consequently, missing for spin-valley models.

Generalized time-reversal symmetry
For su(2) spin models, a genuinely physical symmetry is the invariance under time-reversal. In this setting, time-reversal reverses the sign of all spin operators S a → −S a and, as it is an anti-unitary symmetry, additionally applies complex conjugation to all complex numbers. Hamiltonians with real couplings in which spin operators only appear in pairs are therefore always invariant under time-reversal. On the Hilbert space of the su(2) pseudo-fermions, it can be represented as We again consider a straightforward generalization of the transformation to spin-valley operators, which we define as the anti-unitary transformation Using the relation e iπ(α−α )/2 = αα and Eq. (34), it is straightforward to show that the spin-valley operator transforms as which, up to a minus sign, is the same transformation behavior as for the particle-hole symmetry in Eq. (33). As only pairs of spin-valley operators appear in the spin-valley Hamiltonian, for which the minus sign is irrelevant, the arguments for the invariance of Hamiltonian given there, consequently, also apply here. Applying this generalized version of time-reversal to the local one-particle correlator implies where the complex conjugation stems from the fact that the transformation is anti-unitary. For the bi-local twoparticle correlation function, it implies Apart from extra valley indices, this is exactly the same as in the su(2) case.

Hermitian symmetry
Just as the su(2) spin operator, the spin-valley operator is Hermitian. The spin-valley Hamiltonian only consists of pairs of spin-valley operators and we have restricted ourselves to real couplings, making it Hermitian aswell. Complex transposition, therefore, leaves the Boltzman factor in the thermal expectation value invariant. Applying complex transposition on both sides of Eqs. (26,27) and explicitly evaluating the RHS by "pulling" the complex transpose into the thermal expectation value, we obtain the constraint for the local one-particle correlator and for the two-particle correlator. These constraints are again of the same form as for the su(2) case.

Lattice symmetries
The spin models we consider are all formulated on lattices that can be specified in terms of an underlying Bravais lattice and a possibly multi-atomic basis. Therefore, lattice symmetries exist necessarily for any spin-valley model and are very important to efficiently implement the pf-FRG. Their implementation is the same whether one considers su(2) spin models or spinvalley models. We can therefore use the same approach as for the conventional pf-FRG as, e.g., explained in Ref. [27]. There, all sites are assumed to be identical, in the sense that one can map any site to any other site via a lattice automorphism T that leaves the lattice itself invariant. On the fermionic operators, such a transformation acts as In the case of bond-directional couplings, the transformation would additionally have to be combined with transformations in spin and valley space. For the oneparticle correlation function, this implies The locality constraint in Eq. (28), resulting from the local U(1) symmetry, already reduces the spatial dependence of the self-energy to only one site index i 1 .
Using lattice automorphisms, we can map all sites to an arbitrary reference site i 0 and therefore completely remove the spatial dependence of the one-particle correlation function. Similarly, for the two-particle correlation function, it implies Combining this with the bi-locality constraint in Eq. (29), and again mapping the first index i 1 to an arbitrary reference site i 0 , the spatial dependence of the two-particle correlator can be reduced to just one lattice site.

Parametrization of correlation functions
To make use of the symmetry constraints on the correlation functions, it is advantageous to parametrize them, so that the symmetry constraints manifest in a more practical form. To this end, we can extent the parametrization for the correlation functions for generic su(2) spin models introduced in [27] also to spin-valley models. This ultimately leads to the parametrization of the self-energy and two-particle vertex in Eqs. (15,16) and the symmetry constraints in Eqs. (18,19). Starting with the one-particle correlation function, we argued that due to the local U(1) symmetry and lattice symmetries, it is independent of the lattice site. Additionally, due to Matsubara frequency conservation, which is a consequence of translational invariance in imaginary time, it is diagonal in the frequency arguments. Expanding the spin and valley dependence in Pauli matrices θ μ θ κ (μ, κ = 0, 1, 2, 3), the one-particle correlation function can be parametrized as Similarly, the two-particle correlation function depends only on two lattice sites and three frequencies, for which we choose the three transfer frequencies defined in Eq. (17). The parametrization then reads Plugging this parametrization into the symmetry constraints derived in Sects. 4.1-4.5, we obtain the symmetry constraints for the basis functions of the parametrization listed in Table 1. In the derivation of these constraints, we make heavy use of Eq. (34) and the particle-exchange symmetry which is present in all purely fermionic models.
The labels specify which symmetries of the Hamiltonian were used in their derivation, where H stands for Hermitian, TR for generalized time-reversal, PH for global particle-hole and X for particle-exchange symmetry. The most notable implications are that all correlation functions will always be either only real or imaginary and all expression with negative frequencies can be related to those with positive frequencies The list of symmetry constraints is very similar to the su(2) case derived in Ref. [27], but has two significant differences. First, as already discussed in Sects. 3.2 and 4.2, the symmetry constraint relating s and u frequencies, or the particle-particle and crossed particle-hole channel, is missing because the spin-valley Hamiltonian is not invariant under a local particle-hole transformation but only under the global version. Second, the symmetry constraints do not imply that the oneparticle correlation function is completely diagonal in all spin and valley indices. In the parametrization, this would manifest in G 00 being the only non-vanishing basis function. Instead, for a general spin-valley Hamiltonian, also the terms G ab with a, b > 0, which come with the factor ∼ θ a ss θ b ll , are allowed. This would increase the number of flow equations and therefore also the numerical complexity significantly. Additionally, we could not use the conventional pf-FRG flow equations given in Eqs. (10,11), where a diagonal one particle correlator (and self-energy) was assumed. Fortunately, in the context of moiré materials, many Hamiltonians of physical relevance posses additional symmetries in the spin and valley space [11,35] that further constrain the spin and valley dependence of the self-energy and two-particle vertex. It turns out that the minimal symmetry needed in order for the one-particle correlator to be diagonal is a Z 2 × Z 2 × Z 2 symmetry in either the spin or valley sector. On the level of spin-valley operators, this means that the Hamiltonian is invariant under the transformation (for the case of the spin sector) for each μ individually. This simply reverses the signs of all σ μ i ⊗ τ κ i with μ > 0. Assuming a completely diagonal spin exchange matrix as in Eq. (5), the spin-valley Hamiltonian is indeed invariant under this transformation. This directly implies that all terms proportional to a single ∼ θ μ (with μ > 0) in the correlation functions have to vanish. More precisely, it imposes the constraint which in combination with the first equation in Table 1 implies resulting in a completely diagonal one-particle correlation function parametrized by a single basis function G(ω). For the coupling matrices stated in Eq. (5), we can therefore use the standard pf-FRG approach also for spin-valley models. Assuming this additional symmetry, in the two-particle correlator, only diagonal components in the spin sector ∼ θ μ θ μ (no sum over μ) are allowed, resulting in the constraint Imposing these additional constraints, all factors of ξ(μ)ξ(ν) in Table 1 are equal to one and the relations reduce exactly to the constraints given in Eqs. (18,19) with the self-energy and two-particle vertex replaced by the one-and two-particle correlation functions. We can therefore still consider a completely imaginary one-particle correlator that is odd in frequency space and completely diagonal. The two-particle correlator is either completely real or imaginary, depending on the sign of ξ(κ, λ), and all negative frequency components can be mapped to a positive counterpart. The argument why these constraints on the disconnected correlation functions carry over to the oneparticle irreducible correlation functions, i.e., the selfenergy and the vertex, is the same as given for the su(2) case in [27]. For the self-energy, it simply follows from the Dyson equation [30]: from which it is easy to see that all constraints carry over to the self-energy. For the two-particle vertex, the tree expansion (neglecting the three-particle vertex) relates it to the connected two-particle correlation function G (c) as [30] G (c) (1 , 2 , 1, 2 As the one-particle correlation function is diagonal in all indices, it is clear that all constraints carry over from the connected correlation function to the two-particle vertex. That the constraints from the disconnected twoparticle correlation function carry over to the connected correlation function can be proven by their definition via generating functionals [30].

Symmetries of the flow equations
To verify that the parametrization and the symmetry constraints derived in the previous sections are indeed preserved also for the flowing self-energy and two-particle vertex for any value of Λ, they can additionally be proven using the pf-FRG flow equations given in Eqs. (10,11) . That the parametrization for the self-energy in Eq. (15) and for the two-particle vertex in Eq. (16) is indeed complete can be seen by inserting them into the RHS of the flow equations and confirming that no additional terms are generated. For the additional symmetry constraints, the proof can be performed via induction, as already explained in Refs. [27,39]. This essentially amounts to verifying the fulfillment of the constraints in the initial conditions and then showing that the derivatives d dΛ Σ and d dΛ Γ given by the RHS of the flow equations also fulfill them, assuming the self-energy and two-particle vertex themselves already do. The proof that the self-energy is odd, imaginary, and completely diagonal has to be repeated for spin-valley models due to slight differences in the flow equations. This is quite lengthy and, therefore, done in B. For the two-particle vertex, the proof of the symmetry constraints is much easier on the level of the unparametrized vertex, as there the flow equations still have a much simpler form. We therefore postulate the relations where we defined −1 = (i 1 −ω 1 s 1 l 1 ) and1 = (i 1 ω 1s1l1 ). When translated to the parametrized two-particle vertex and then combined, these relations yield exactly the symmetry constraints given in Eq. (19). Proving the relations for the unparametrized vertex, therefore, directly proves the symmetry constraints of the parametrized vertex. As Eq. (57) simply amounts to a simple particle exchange, no further proof is required. Eq. (58) is proven in [27] and Eq. (59) in [39] using the general pf-FRG flow equations. The only remaining relation still left to prove is Eq. (60), which we also show in B. This proves that the parametrization and the symmetry constraints are indeed valid also for the flowing self-energy and vertex, at any value of the cutoff Λ.

Results
To give an explicit example for the application of the pseudo-fermion functional renormalization group approach introduced in the manuscript and its efficient implementation in terms of the aforementioned symmetries, we apply it to elucidate the phase diagram of an SU(2) spin ⊗ U(1) valley symmetric spin-valley Hamiltonian on the triangular lattice. The explicit Hamiltonian we consider is with a SU(4) symmetric term proportional to the coupling J and an in-plane J x and out-of-plane J z coupling that when non-zero break the SU(4) symmetry down to an SU(2) symmetry in the spin sector and a U(1) symmetry in the valley sector. We only include interactions between nearest neighbours ij . Such a model can be motivated, e.g., from including the effect of Hund's type couplings in a two-orbital extended Hubbard model and performing a strong coupling expansion [20]. It can therefore be regarded as a natural extension to previously studied models with either full SU(4) or reduced SU(2) spin ⊗ SU(2) valley symmetry [18][19][20][21] by adding an XXZ type perturbation to the orbital sector and likewise provides an intermediate, but important step towards the more complicated spin-valley Hamiltonians proposed for various moiré systems [10,11,20].

Phase diagram
To obtain the quantum phase diagram, we fix the coupling J in front of the SU(4) symmetric part of the Hamiltonian in Eq. (61) to a positive value and then vary the values of the in-plane coupling J x and outof-plane coupling J z which break the SU(4) symmetry. As described in Sect. 3, to determine the magnetic order for a particular pair of couplings (J x , J z ) we calculate the flow of the spin-spin and valley-valley correlations (and associated structure factors) as defined in Eqs. (13,14), check whether or not a flow breakdown occurs and if so, which type of order is visible in the structure factor at the breakdown scale Λ c .  Fig. 3. This indicates that no magnetic order is present in both the spin and the valley sector even for very-lowenergy scales and indicates a putative quantum spinvalley liquid (QSVL) state [21]. Going away from the SU(4) point, however, we almost immediately observe a flow breakdown in either the spin or valley sector, indicating that the putative QSVL state is highly unstable in the presence of XXZ like perturbations. This is in line with results for the su(2) XXZ model on the triangular lattice, where by varying the out-of-plane coupling a phase transition from in-plane 120 • order to an "umbrella" order is observed at the SU(2) symmetric point [40]. Similarly, we observe a rich ensemble off different spin and valley ordered phases with both in-and out-of-plane ordering in the valley sector.
Before we present the full quantum phase diagram, however, let us first consider a classical mean-field approach to better understand the origin of the different phases. To this end, we note that the spin sector by itself will order either ferromagnetically (FM) or in a 120 • order, depending on the sign of the exchange coupling. Assuming one of these states is realized, we decouple the spin and valley sector by approximating the pair of spin operators by its expectation value where the spin expectation value only appears in the positive factor E s 0 ≡ J(1 + σ i σ j ) and, therefore, has no influence on the type of valley order. Approximating the valley operators by classical vectors with |τ i | = 1, we perform a Luttinger-Tisza analysis [41,42] on the mean-field Hamiltonian. This analysis predicts in-plane (out-of-plane) valley order for large values of |1 + J x /J| (|1 + J z /J|), which is either FM for positive, or 120 • like for negative values. The precise phase boundaries along with the ground-state energies E g are depicted in with antiferromagnetic coupling J. Here, the flow of the spin structure factor shows a sharp peak, while the valley structure factors are strongly suppressed, as is depicted in Fig. 3. The same behavior occurs in a larger region around J x /J = J z /J = −1, which is shown in Fig. 5 along with the corresponding momentum resolved structure factor (annotated with the numeral I). The spin structure factor (in purple) shows strong peaks at the K and K points, while the in-plane (orange) and out-of-plane (blue) valley structure factors show no distinct features when shown on the same color scale. This indicates 120 • spin order, which again agrees with the results for the conventional su(2) Heisenberg model [40,43]. In all other regions of the quantum phase diagram, the valley structure factors are clearly dominant and the spin structure factor shows only weak features. We enumerate the different types of valley order we find by the numerals II-VI, as shown in Fig. 5. The valley order at large negative couplings (II and III) agrees with the classically predicted results, as either the in-or out-of-plane structure factors show strong peaks at the Γ point, indicating FM order. At larger positive values for either the in-plane or out-of-plane coupling (VI-IV), the valley structure factors show peaks at the K and K points indicating 120 • like order. In contrast to the sharp phase boundary in the classical case, however, the valley order seems to gradually shift from mostly in-plane (IV), over competing in-and out-of-plane (V) to out-of-plane (VI) order. This is well visualized by the flow of the structure factors in Fig. 6. The valley structure factors both show flow breakdowns at approximately the same breakdown scale, but the magnitude at the breakdown scale shifts from a dominant χ Λv,x to a dominant χ Λv,z when going from IV to VI. To quantify this transition, we define the angle illustrated in Fig. 5 by the cones on the right and by the color scale ranging from from blue (in-plane) over green (competing in-plane and out-of-plane) to orange (out-of-plane) To better illustrate the transitions between the different types of order, Figs. 7 and 8 show horizontal and vertical cuts through the phase diagram, respectively. The transitions between the phases are always accompanied by a dip or kink in the breakdown scale, indicating the positions of the phase boundaries. This is especially relevant for the transitions between mixed in-and out-of-plane valley order (V) to dominant in-or out-of-plane valley order (IV and VI), where the phase transition would not be as easily recognizable by just considering the evolution of the structure factors. The same is true for the transition from dominant valley to dominant spin order, as, e.g., depicted in the J x /J = 0 cut in Fig. 8. Here, the at first very dominant out-ofplane valley order (III) gradually transitions to dominant spin order (I), with a region in between where the spin and valley structure factors are of similar magnitude. The kink in the breakdown scale appears at the largest J z /J where the valley structure factor still shows a clear flow breakdown (J z /J ≈ −1.6), even though the spin structure factor is already dominant for smaller J z /J. This is similar at all boundaries of phase I, which also becomes evident in the phase diagram of Fig. 5 by noting that the minima of the breakdown scale are positioned slightly inwards in the region of dominant spin order (colored in purple).
Of special interest are the cuts across the SU(4) point (J x /J = 0 and J z /J = 0), which show that even for very small values of the in-and out-of-plane couplings, the flow develops a breakdown at a finite Λ c .

Summary
In this manuscript, we have presented a generalization of the established pf-FRG approach to generic spinvalley Hamiltonians in the self-conjugate representation of su(4), with either diagonal spin or valley interactions. We performed a careful symmetry analysis and derived a set of constraints on the vertex functions, which drastically lower the computational cost of tracking the flow of running couplings. Using a highly accurate solver for the functional flow equations, we subsequently applied this method to map out the quantum phase diagram of an SU(2) spin ⊗ U(1) valley model on the triangular lattice, which presents a simplified variant of the more general Hamiltonian proposed for TLG/h-BN, but already hosts a rich variety of spin and valley ordered ground states. In addition, we were able to demonstrate, that, by promoting the spin symme-  (2) to SU(4), quantum fluctuations are boosted, ultimately resulting in a smooth RG flow down to the lowest energy scales, indicative of a spinvalley liquid state. However, this QSVL state appears to be very sensitive even to weak XXZ anisotropies in the valley sector, and we almost immediately detect the emergence of long-range order, when perturbing it. While our focus in this manuscript has been on spinvalley Hamiltonians, we note that very similar models have been discussed for spin-orbit coupled systems that go beyond the celebrated Kugel-Khomskii model. The microscopic ingredients of such spin-orbital models are surprisingly similar to those of "Kitaev materials" [44]-a partially filled 4d or 5d orbital, the formation of a spin-orbital entangled local moment, and an edge-sharing octahedral crystalline environment. Specifically, a d 1 configuration can lead to local j = 3/2 moments subject to bond-directional exchanges that break the original SU(4) symmetry of the j = 3/2 moments. As a concrete material candidate exhibiting this microscopic mechanism, α-ZrCl 3 -a 4d sister compound of the isostructural Kitaev material RuCl 3 -has been put forward [45]. To study the phase diagram of spin-orbital ground states in such a setting with varying diagonal and off-diagonal couplings, one can again rely on the pseudo-fermion FRG approach put forward in this manuscript. use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix B: Proof of symmetry constraints via flow equations
In Sect. 4.7, we claim that the completeness of the parametrization given in Eqs. (15,16) and the symmetry constraints given in Eqs. (18,19) can also be proven by induction using the flow equations, as was already done for the su2 case [27,39]. The proof amounts to checking that the constraints are fulfilled in the initial conditions and then showing that the RHS of the pf-FRG flow equations in Eqs. (10,11) also fulfill the constraints, assuming the selfenergy and two-particle vertex themselves already do. That the constraints are fulfilled in the initial conditions is easy to see, as for Λ → ∞ the two-particle vertex is frequency independent and the self-energy vanishes. We will, therefore, only perform the induction step here. Starting with the two-particle vertex, it is straightforward to see that the parametrization is complete by plugging it in the pf-FRG flow equations and showing that no additional terms are generated. To proof the symmetry constrains, we postulated equivalent constraints for the unparametrized two-particle vertex in Eqs. (57-60), which, when combined, lead to the symmetry constraints of the parametrized vertex. Fortunately, only the relation Γ Λ (1 , 2 , 1, 2) = s 1 s1l 1 l1s 2 s2l 2 l2Γ Λ (1,2,1 ,2 ) (B.5) differs from the su(2) case and all other relations have already been proven [27,39]. The induction step for this relation is performed by writing down the flow equations for s 1 s1l 1 l1s 2 s2l 2 l2Γ Λ (1,2,1 ,2 ) and then manipulating the RHS