Abstract
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 \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(2) spins is capable of treating diagonal and off-diagonal couplings of generic spin–valley exchange Hamiltonians in the self-conjugate representation of the \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(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)\(_{\text {spin}}\) \(\otimes \) U(1)\(_{\text {valley}}\) model on the triangular lattice which exhibits a rich phase diagram of spin and valley ordered phases.
Graphic Abstract
Similar content being viewed by others
1 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.
a Moiré pattern emerging in two stacked layers of graphene with a relative twist angle \(\theta \). 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
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 models can then be mapped to \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(4)Footnote 1 spin–valley Hamiltonians that may be used as a starting point to investigate the nature of the correlated insulating states. The so-derived \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(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 \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(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 \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(4) spin–valley (or spin–orbital) 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, self-conjugate representation of \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(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)\(_{\text {spin}}\) \(\otimes \) U(1)\(_{\text {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 \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(4) models as well as the implementation of model specific symmetries (Sect. 4). Finally, numerical results for the phase diagram of a SU(2)\(_{\text {spin}}\) \(\otimes \) U(1)\(_{\text {valley}}\) model on the triangular lattice are presented and examined in Sect. 5.
2 Model
Microscopically, the SU(4) models of interest in this manuscript can be cast in terms of a general Hamiltonian
which couples two elementary \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(2) degrees of freedom, captured by the operators \(\sigma \) and \(\tau \), which might denote a spin and valley (or oribtal) degree of freedom. 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 \(\varvec{\sigma }_i \varvec{\sigma }_j\) on sites i and j (with \(\varvec{\sigma }_i = (\sigma _j^x, \sigma _j^y, \sigma _j^z)^T\)) and a similar interaction of the valley degrees of freedom \(\varvec{\tau }_i \varvec{\tau }_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 \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(4), where the spin–valley operators can be represented in terms of fermionic creation and annihilation operators as
with a local half-filling constraint
subject to every lattice site, where summation over repeated spin indices s and valley indices l is implied. Here, \(\theta ^{\mu }\) denotes a Pauli matrix with \(\mu \in \{0, 1, 2, 3\}\) and \(\theta ^{0} = {\mathbbm {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}^{\mu \nu } \otimes J_{v, ij}^{\kappa \lambda }\) is understood as the Kronecker product of the spin and valley exchange matrices and summation over repeating \(\mu ,\nu ,\kappa \) or \(\lambda \) is again implied. Here, \({\hat{n}}_i\) is the density operator \({\hat{n}}_i \equiv \sigma _i^0 \tau _i^0 = f^\dagger _{isl}f^{{\dagger }}_{isl}\), and the term proportional to the coupling \(I_{ij}\) is needed to potentially cancel the density term \(\sim \sigma _i^0 \tau _i^0 J^{00}_{s, ij}J^{00}_{v, ij} \sigma _j^0 \tau _j^0\), which does not appear in pure \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \({\mathbb {Z}}_2 \times {\mathbb {Z}}_2 \times {\mathbb {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-neighbour (\(\sim J_1\)) and next-nearest-neighbour (\(\sim J_2\)) interactions, contains both diagonal \(\sim J^1_{p, ij}\) and off-diagonal \(\sim J^2_{p,ij}\) valley exchange that breaks the SU(4) symmetry down to an SU(2)\(_\mathrm {spin}\) \(\otimes \) U(1)\(_\mathrm {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).
3 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 off-diagonal couplings as defined by Eq. (5). To set the stage, we will first revisit the flow equations of the conventional pf-FRG approach for \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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.
3.1 Pseudo-fermion functional renormalization group
Let us set the stage by revisiting some of the conceptual steps of the pseudo-fermion FRG, which has originally been formulated for bilinear \({{\mathfrak {s}}}{{\mathfrak {u}}}(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
with \(G_0(\omega ) = (i\omega )^{-1}\). The multi-index \(1=(i_1, s_1, l_1, \omega _1)\) consists of a lattice site index \(i_1\), a spin index \(s_1\), a fermionic Matsubara frequency \(\omega _1\) and, for spin–valley models, the additional valley index \(l_1\). To implement the RG scale, or cut-off, \(\varLambda \), 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 hereFootnote 2. In the Katanin truncation, only the self-energy \(\varSigma ^\varLambda \) and the two-particle vertex \(\varGamma ^\varLambda \) 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 \(S^{\varLambda } \equiv -\partial _\varLambda G^{\varLambda }|_{\varSigma ^{\varLambda } = \text {const.}}\). Note that the flow equations are formulated in the \(T \rightarrow 0\) limit and the sums should therefore be understood as \(\sum _{1} \equiv \sum _{i_1 s_1 l_1} \int d\omega _1\).
Flow of the spin and valley structure factor in a magnetically ordered phase (a) and a paramaganetic phase (b) for different values of the vertex truncation length L. 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 \(\varLambda /J = 0.02\), which is the smallest scale for which our calculations are numerically reliable
The fermionic representation of the spin–valley operators, as presented in the previous section, necessitates the enforcement of a local half-filling constraint \((f^\dagger _{isl}f^{{\dagger }}_{isl} = 2)\) to determine the dimensionality of the local Hilbert space. To this end, we employ the same technique used for \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \(\langle f^\dagger _{isl}f^{{\dagger }}_{isl} \rangle \) 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) spin models, the physically relevant filling for spin–valley models depends on the considered material and may also be at quarter-filling \((f^\dagger _{isl}f^{{\dagger }}_{isl} = 1)\) [11, 35], corresponding to the fundamental representation of \({{\mathfrak {s}}}{{\mathfrak {u}}}(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_{\tau }\) is the time-ordering operator. From this general definition, we can then read off the form of spin–spin correlations
as well as valley–valley correlations
A thermal phase transition to long-range, symmetry-breaking 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 \(\varLambda _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 \(\chi ^{\varLambda s/v}_{ij}(\omega = 0)\), then indicates the type of symmetry-breaking. Since the solution of the flow equation below the breakdown scale \(\varLambda _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, mean-field 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 \(\chi ^{\varLambda 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 \(\chi ^{\varLambda 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 (\(\varLambda = 0.02 J\)), 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).
3.2 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 parametrized by a single function \(\varSigma (\omega )\) 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) spin models—apart from an increased number of components due to the valley sector \(\sim \theta _{l_{1'}l_1}^\kappa \theta _{l_{2'}l_2}^\lambda \) with the corresponding indices \(\kappa \) and \(\lambda \). If we assume the Hamiltonian to be diagonal in the spin sector, we will only need to consider components diagonal in the spin \(\sim \theta _{s_{1'}s_1}^\mu \theta _{s_{2'}s_2}^\mu \), with the corresponding index \(\mu \) (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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 two-particle irreducibility. This results in a particle–particle (pp), direct particle–hole (dph), and crossed particle–hole (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 \(\varGamma ^{\varLambda \rightarrow \infty }\) is the bare two-particle vertex at infinite cut-off. Each channel \(g_c(\omega _c, v_c, v^\prime _c)\) is parametrized by one bosonic transfer frequency \(\omega _c\) and two additional fermionic frequencies \(v^{{\prime }}_c, v^\prime _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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 dependenceFootnote 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 \(\varLambda \rightarrow \infty \), which are given by
with the couplings \(J^\mu _{s,i_1i_2}\) and \(J^{\kappa \lambda }_{i_1i_2}\) defined in Eq. (5).
3.3 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 \(\varGamma ^\varLambda _{i_0 j}\) (considering only the lattice site dependence), we define a finite length scale L and truncate the vertex \(\varGamma ^\varLambda _{i_0, 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 \sim 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 linear-logarithmic meshes, which are updated using a scanning routine between each step of the ordinary differential 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_\varOmega \) and (positive) fermionic frequencies \(N_\nu \) as \({\mathcal {O}}(N_\varOmega \cdot N_\nu ^2)\). A typical setup for which the two-particle vertex is sufficiently well resolved is \(N_\varOmega = 40\) and \(N_\nu =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_\varSigma = 250\) frequencies. In the \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 particle–particle and crossed particle–hole channel (discussed in Sect. 3.2), we have to also consider negative frequencies for either \(\nu _c\) or \(\nu '_c\). This results in an additional factor of two in computation time compared to \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) case, more complex spin–valley models usually require a much larger numerical effort, as the extra valley index greatly increases the number of independent two-particle vertex components \(\varGamma ^{\varLambda , \mu \kappa \lambda }_{i_1i_2}\), in which the computation time scales linearly. With the coupling matrices given in Eq. (5), there would be \(N_\varGamma = 4 \cdot 4^2 = 64\) independent vertex components (only considering the spin–valley dependence). In comparison, the parametrization for generic \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) models only has \(N_\varGamma = 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_\varGamma = 2 \cdot 6 = 12\). For these models, the numerical effort is similar to \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) models with off-diagonal interactions and even allows for computations of relatively large-phase diagrams as will be presented in Sect. 5.
4 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 particle–hole symmetry that is present in the \({{\mathfrak {s}}}{{\mathfrak {u}}}(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.
4.1 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 \(\varphi _i \in [0, 2\pi )\) to the fermionic operators as
which clearly leaves all spin–valley operators invariant. Interpreting the spin–valley Hamiltonian as a fermionic representation of an \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 pseudo-fermion representation. For our functional renormalization group approach, we are interested in the implication of the symmetry on the functional formFootnote 4 of the one-particle correlation function
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\).
4.2 Global particle–hole symmetry
In the pf-FRG for \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) spin models, spin operators \(S_i^a\) are represented using fermions with one spin index \(\alpha = \pm 1\) as
with \(a \in \{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 \({\bar{\alpha }}\equiv -\alpha \). It leaves the fermionic representation of the \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 = \pm 1\) and valley index \(l = \pm 1\) is the transformation
under which the spin–valley operator transforms as
which can be shown straightforwardly using the anti-commutation relations of the fermionic operators and the identity
Spin–valley operators with either the spin index \(\mu \) or the valley index \(\kappa \) set to zero—which correspond to the individual spin and valley operators as defined in Eq. (2)—are invariant under this transformation. General spin–valley operators, on the other hand, may change their sign. The Hamiltonian is, therefore, not invariant under the local particle–hole symmetry that acts on the Hilbert space of just one lattice site. Spin–valley operators, however, only appear in pairs in the spin–valley Hamiltonian. Performing the local particle–hole symmetry transformation on all lattice sites, such a pair of spin–valley operators transform as
If an odd number of spin and valley indices is set to zero, this again implies a sign change. Recalling the definition of the spin–valley Hamiltonian in Eq. (4) and the following definition of the exchange matrices in Eq. (5), such terms are not included in our definition of the Hamiltonian. All terms that do appear in the Hamiltonian are indeed invariant. The main difference to the \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) pseudo-fermion Hamiltonian is, therefore, that the spin–valley Hamiltonian is invariant only under the global transformation, while the former was invariant under the local transformation. For the local single-particle correlation function, the global particle–hole 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 particle–hole channel with each other. As already discussed in Sect. 3, this constraint is, consequently, missing for spin–valley models.
4.3 Generalized time-reversal symmetry
For \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \rightarrow - 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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\pi (\alpha - \alpha ')/2} = \alpha \alpha '\) 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 two-particle correlation function, it implies
Apart from extra valley indices, this is exactly the same as in the \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) case.
4.4 Hermitian symmetry
Just as the \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) case.
4.5 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) spin models or spin–valley 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 one-particle 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.
4.6 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \(\theta ^\mu \theta ^\kappa \) (\(\mu ,\kappa = 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 list of symmetry constraints is very similar to the \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 one-particle 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 \(\sim \theta ^a_{ss'}\theta ^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 \({\mathbb {Z}}_2 \times {\mathbb {Z}}_2 \times {\mathbb {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 \(\mu \) individually. This simply reverses the signs of all \(\sigma ^\mu _i \otimes \tau ^\kappa _i\) with \(\mu > 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 \(\sim \theta ^\mu \) (with \(\mu >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(\omega )\). 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 \(\sim \theta ^\mu \theta ^\mu \) (no sum over \(\mu \)) are allowed, resulting in the constraint
Imposing these additional constraints, all factors of \(\xi (\mu )\xi (\nu )\) 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 \(\xi (\kappa , \lambda )\), 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 one-particle irreducible correlation functions, i.e., the self-energy and the vertex, is the same as given for the \({{\mathfrak {s}}}{{\mathfrak {u}}}(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]
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 two-particle correlation function carry over to the connected correlation function can be proven by their definition via generating functionals [30].
4.7 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 \(\varLambda \), 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 \(\frac{d}{d\varLambda }\varSigma \) and \(\frac{d}{d\varLambda }\varGamma \) 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 -\omega _1 s_1 l_1)\) and \({\bar{1}} = (i_1 \omega _1 {\bar{s}}_1 {\bar{l}}_1)\). 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 cut-off \(\varLambda \).
5 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)\(_\mathrm {spin}\) \(\otimes \) U(1)\(_\mathrm {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 \(\langle ij \rangle \).
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)\(_\mathrm {spin}\) \(\otimes \) SU(2)\(_\mathrm {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].
5.1 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 out-of-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 \(\varLambda _c\). Due to the SU(2)\(_\mathrm {spin}\) \(\otimes \) U(1)\(_\mathrm {valley}\) symmetry of the Hamiltonian all non-vanishing components of the spin–spin correlation are equivalent and we calculate only \(\chi ^{\varLambda s}_{ij} \equiv \chi ^{\varLambda s, xx}_{ij} = \chi ^{\varLambda s, yy}_{ij} = \chi ^{\varLambda s, zz}_{ij}\). For the valley–valley correlation, we can distinguish between in-plane and out-of-plane order by calculating the in-plane valley–valley correlation \(\chi ^{\varLambda v, x}_{ij} \equiv \chi ^{\varLambda v, xx}_{ij} = \chi ^{\varLambda v, yy}_{ij}\) and out-of-plane valley-valley correlation \(\chi ^{\varLambda v, z}_{ij} \equiv \chi ^{\varLambda v, zz}_{ij}\).
Starting at the SU(4) point with \(J_x/J = J_z/J = 0\), where all spin–spin and valley–valley correlations are equivalent, we observe no flow breakdown of the structure factors, as depicted by the grey line in Fig. 3. This indicates that no magnetic order is present in both the spin and the valley sector even for very-low-energy scales and indicates a putative quantum spin–valley 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 \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) XXZ model on the triangular lattice, where by varying the out-of-plane coupling a phase transition from in-plane \(120^\circ \) 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.
Flow of the structure factor at points of higher symmetry. All structure factors are shown at the momentum where they are maximal. The grey line shows the structure factor at the SU(4) point, where the considered spin–valley model corresponds to the SU(4) symmetric Heisenberg model. Here, all structure factor components are identical. The flow is smooth and convex down to the lowest numerically reliable cut-off and no flow breakdown occurs, indicating a putative quantum spin–valley liquid (QSVL) ground state. The purple and green lines show the spin and valley structure factor for \(J_x/J = J_z/J = -1\), where all terms containing valley operators cancel and the spin–valley model resembles an SU(2) symmetric Heisenberg model. In this case, the valley structure factors are strongly suppressed and the spin structure factor shows a sharp peak at the K and \(K'\) points, indicating \(120^\circ \) order in the spin 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^\circ \) 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 \(\varvec{\sigma }_{i}\varvec{\sigma }_{j} \approx \langle \varvec{\sigma }_{i}\varvec{\sigma }_{j}\rangle \), with \(\langle \varvec{\sigma }_{i}\varvec{\sigma }_{j}\rangle = 1\) for ferromagnetic (FM) and \(\langle \varvec{\sigma }_{i}\varvec{\sigma }_{j}\rangle = \cos (2\pi /3) = -0.5\) for \(120^\circ \) order. The resulting mean-field Hamiltonian is then given, up to a constant, by
where the spin expectation value only appears in the positive factor \(E_0^s \equiv J(1+\langle \varvec{\sigma }_{i}\varvec{\sigma }_{j}\rangle )\) and, therefore, has no influence on the type of valley order. Approximating the valley operators by classical vectors with \({\left| \varvec{\tau }_i\right| }=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 \({\left| 1+J_x/J\right| }\) (\({\left| 1+J_z/J\right| }\)), which is either FM for positive, or \(120^\circ \) like for negative values. The precise phase boundaries along with the ground-state energies \(E_g\) are depicted in Fig. 4.
Classical phase diagram in valley space for fixed spin ordering obtained from a Luttinger–Tisza analysis. The grey lines depict the phase boundaries and the color illustrates the (normalized) ground-state energy, where blue denotes out-of-plane and orange denotes in-plane ordering. At \(J_x/J = J_z/J = - 1\), where the phase boundaries meat, the classical mean-field Hamiltonian vanishes. Away from this point the Luttinger–Tisza analysis predicts the following types of valley order: (II) in-plane ferromagnetic (FM), (III) out-of-plane FM, (IV) in-plane \(120^\circ \), and (VI) out-of-plane \(120^\circ \). The so-obtained valley order is independent from the fixed nearest-neighbour spin order
Phase diagram of the SU(2) \(_\mathrm {spin}\) \(\otimes \) U(1)\(_\mathrm {valley}\) symmetric spin–valley model. The color indicates which structure factor is dominant at the breakdown scale \(\varLambda _c\), where purple implies dominant spin order, colors between orange and yellow indicate dominant order in the valley sector, and the opacity determines the magnitude of the breakdown scale \(\varLambda _c/J\). In the case where we observe a flow breakdown in both the in-plane (\(\chi ^{\varLambda v, x}\)) and out-of-plane (\(\chi ^{\varLambda v, x}\)) valley structure factor the color determines the angle \(\phi \) illustrated in the cones on the right of the figure. (I–VI) show the structure factors at \(\varLambda _c\) for the different types of order we observe: (I) 120\(^\circ \) spin order, (II) out-of-plane FM valley order, (III) in-plane FM valley order, (IV-VI) 120\(^\circ \) valley order shifting from an out-of-plane (IV) to an in-plane (VI) orientation, with competing order (V) in between. For \(J_x/J = J_z/J = 0\), indicated by the star, the model is equivalent to the SU(4) symmetric Heisenberg model for which no flow breakdown is observed
Flow of the structure factors for different types of order as described under Fig. 5. The dashed lines show the breakdown scale \(\varLambda _c\). (I) shows dominant spin order (purple) with the valley structure factors strongly suppressed. (IV–VI) shows dominant valley order which shifts from an in-plane (orange) to an out-of-plane (blue) orientation
Special attention needs to be paid to the point at \(J_x/J = J_z/J = -1\) where the phase boundaries meet. Exactly at this point, the couplings in front of the valley operators are equal to zero and the mean-field Hamiltonian vanishes. Going back to the full quantum Hamiltonian, it reduces to only the term \(\sum _{ij} J (1 + \varvec{\sigma }_i \varvec{\sigma }_{j})\), which resembles an SU(2) symmetric Heisenberg model 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^\circ \) spin order, which again agrees with the results for the conventional \({{\mathfrak {s}}}{{\mathfrak {u}}}(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 \(\varGamma \) 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^\circ \) 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 \(\chi ^{\varLambda v, x}\) to a dominant \(\chi ^{\varLambda 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)
Horizontal cuts through the phase diagram of Fig. 5 at \(J_z/J = 1\) (top), \(J_z/J = 0\) (middle) and \(J_z/J = -1\) (bottom). The color-coding and the labels (I-VI) denote different types of order and are explained in Fig. 5. The transition between these phases is always accompanied by a dip or kink in the breakdown scale. In the transitions between dominant spin and valley order, there are regions where both the spin and valley structure factor show flow breakdowns at a similar \(\varLambda _c\) and with similar magnitudes. These regions are shaded and colored both orange and purple, as, e.g., visible in the transition between II and I along the \(J_z/J =-1.0\) axis. At the SU(4) point (colored white), no flow breakdown occurs
Vertical cuts through the phase diagram of Fig. 5 at \(J_x/J = 1\) (top), \(J_x/J = 0\) (middle) and \(J_x/J = -1\) (bottom). The color-coding and the labels (I–VI) denote different types of order and are explained in Fig. 5. The transition between these phases is always accompanied by a dip or kink in the breakdown scale. In the transitions between dominant spin and valley order, there are regions where both the spin and valley structure factor show flow breakdowns at a similar \(\varLambda _c\) and with similar magnitudes. These regions are shaded and colored both orange and blue, as e.g. visible in the transition between III and I along the \(J_x/J =0\) axis. The white regime close to the SU(4) point marks points for which no flow breakdown is observed
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-of-plane 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 \approx -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 \(\varLambda _c\).
6 Summary
In this manuscript, we have presented a generalization of the established pf-FRG approach to generic spin–valley Hamiltonians in the self-conjugate representation of \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(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)\(_{\text {spin}}\) \(\otimes \) U(1)\(_{\text {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 symmetry group from SU(2) to SU(4), quantum fluctuations are boosted, ultimately resulting in a smooth RG flow down to the lowest energy scales, indicative of a spin–valley 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 spin–valley 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, \(\alpha \)-\(\hbox {ZrCl}_3\) – a 4d sister compound of the isostructural Kitaev material \(\hbox {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.
Data availability statement
This manuscript has no associated data or the data will not be deposited. [Authors’ comment: The numerical data generated and analyzed in this paper is available from the corresponding author on reasonable request.]
Notes
With \({{\mathfrak {s}}}{{\mathfrak {u}}}\)(4), we refer to the Lie algebra of the Lie group SU(4).
More recently, an alternative multi-loop truncation has been introduced in the context of electronic FRG calculations [32], which was subsequently also adapted in the context of pf-FRG [23, 24]. Such a multi-loop approach can also be applied in the context of spin–valley pf-FRG calculations, but will be left to future exploration. Additionally, one may further improve the efficiency of the method by adapting the cluster FRG scheme [33].
Fortunately, as we will explain in Sect. 3.3, this only results in an increase of numerical complexity by a factor of two, making numerical calculations only slightly more costly.
Note that our definition deviates from normal ordering to be in line with the conventional definition of retarded Greens functions.
References
Y. Cao, V. Fatemi, A. Demir, S. Fang, S.L. Tomarken, J.Y. Luo, J.D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R.C. Ashoori, P. Jarillo-Herrero, Nature 556(7699), 80 (2018). https://doi.org/10.1038/nature26154
C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe et al., Nat. Phys. 16(5), 520–525 (2020). https://doi.org/10.1038/s41567-020-0825-9
X. Liu, Z. Hao, E. Khalaf, J.Y. Lee, Y. Ronen, H. Yoo, D. Haei Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, et al., Nature 583(7815), 221–225 (2020). https://doi.org/10.1038/s41586-020-2458-7
Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J.M. Park, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, Nature 583(7815), 215–220 (2020). https://doi.org/10.1038/s41586-020-2260-6
Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, P. Jarillo-Herrero, Nature 556(7699), 43 (2018). https://doi.org/10.1038/nature26160
X. Lu, P. Stepanov, W. Yang, M. Xie, M.A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A.H. MacDonald, D.K. Efetov, Nature 574(7780), 653 (2019). https://doi.org/10.1038/s41586-019-1695-0
G. Chen, A.L. Sharpe, P. Gallagher, I.T. Rosen, E.J. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, J. Jung, Z. Shi, D. Goldhaber-Gordon, Y. Zhang, F. Wang, Nature 572(7768), 215 (2019). https://doi.org/10.1038/s41586-019-1393-y
A.L. Sharpe, E.J. Fox, A.W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M.A. Kastner, D. Goldhaber-Gordon, Science 365(6453), 605 (2019). https://doi.org/10.1126/science.aaw3780
M. Koshino, N.F.Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, L. Fu, Phys. Rev. X 8, 031087 (2018). https://doi.org/10.1103/PhysRevX.8.031087
L. Xian, D.M. Kennes, N. Tancogne-Dejean, M. Altarelli, A. Rubio, Nano Lett. 19(8), 4934–4940 (2019). https://doi.org/10.1021/acs.nanolett.9b00986
Y.H. Zhang, T. Senthil, Phys. Rev. B 99, 205150 (2019). https://doi.org/10.1103/PhysRevB.99.205150
C. Xu, L. Balents, Phys. Rev. Lett. 121(8), 087001 (2018). https://doi.org/10.1103/PhysRevLett.121.087001
N.F.Q. Yuan, L. Fu, Phys. Rev. B 98, 045103 (2018). https://doi.org/10.1103/PhysRevB.98.045103
H.C. Po, L. Zou, A. Vishwanath, T. Senthil, Phys. Rev. X 8, 031089 (2018). https://doi.org/10.1103/PhysRevX.8.031089
L. Classen, C. Honerkamp, M.M. Scherer, Phys. Rev. B 99(19), 195120 (2019). https://doi.org/10.1103/PhysRevB.99.195120
K.I. Kugel, D.I. Khomskiĭ, Sov. Phys. Usp. 25(4), 231 (1982). https://doi.org/10.1070/pu1982v025n04abeh004537
L.F. Feiner, A.M. Oleś, J. Zaanen, Phys. Rev. Lett. 78, 2799 (1997). https://doi.org/10.1103/PhysRevLett.78.2799
P. Corboz, M. Lajkó, A.M. Läuchli, K. Penc, F. Mila, Phys. Rev. X 2, 041013 (2012). https://doi.org/10.1103/PhysRevX.2.041013
W.M.H. Natori, E.C. Andrade, R.G. Pereira, Phys. Rev. B 98(19), 195113 (2018). https://doi.org/10.1103/PhysRevB.98.195113
W.M.H. Natori, R. Nutakki, R.G. Pereira, E.C. Andrade, Phys. Rev. B 100(20) (2019). https://doi.org/10.1103/PhysRevB.100.205131
D. Kiese, F.L. Buessen, C. Hickey, S. Trebst, M.M. Scherer, Phys. Rev. Res. 2(1), 013370 (2020). https://doi.org/10.1103/PhysRevResearch.2.013370
J. Reuther, P. Wölfle, Phys. Rev. B 81(14), 144410. https://doi.org/10.1103/PhysRevB.81.144410
J. Thoenniss, M.K. Ritter, F.B. Kugler, J. von Delft, M. Punk, arXiv:2011.01268
D. Kiese, T. Mueller, Y. Iqbal, R. Thomale, S. Trebst, arXiv:2011.01269
R. Bistritzer, A.H. MacDonald, PNAS 108(30), 12233 (2011). https://doi.org/10.1073/pnas.1108174108
D.I. Khomskii, Transition Metal Compounds (Cambridge University Press, 2014)
F.L. Buessen, V. Noculak, S. Trebst, J. Reuther, Phys. Rev. B 100(12), 125164 (2019). https://doi.org/10.1103/PhysRevB.100.125164
F.L. Buessen, D. Roscher, S. Diehl, S. Trebst, Phys. Rev. B 97(6), 064415 (2018). https://doi.org/10.1103/PhysRevB.97.064415
C. Wetterich, Phys. Lett. B 301, 90 (1993). https://doi.org/10.1016/0370-2693(93)90726-X
P. Kopietz, L. Bartosch, F. Schütz, Introduction to the Functional Renormalization Group (Springer, Berlin, 2010)
A.A. Katanin, Phys. Rev. B 70(11), 115109 (2004). https://doi.org/10.1103/PhysRevB.70.115109
F.B. Kugler, J. von Delft, Phys. Rev. B 97(3), 035162 (2018). https://doi.org/10.1103/PhysRevB.97.035162
J. Reuther, R. Thomale, Phys. Rev. B 89, 024412 (2014). https://doi.org/10.1103/PhysRevB.89.024412
F.L. Buessen, M. Hering, J. Reuther, S. Trebst, Phys. Rev. Lett. 120, 057201 (2018). https://doi.org/10.1103/PhysRevLett.120.057201
J.W.F. Venderbos, R.M. Fernandes, Phys. Rev. B 98(24), 245103 (2018). https://doi.org/10.1103/PhysRevB.98.245103
A. Keselman, L. Savary, L. Balents, SciPost Phys. 8, 76 (2020). https://doi.org/10.21468/SciPostPhys.8.5.076
N. Wentzell, G. Li, A. Tagliavini, C. Taranto, G. Rohringer, K. Held, A. Toschi, S. Andergassen, Phys. Rev. B102(8), 085106 (2020). https://doi.org/10.1103/PhysRevB.102.085106
P. Bogacki, L.F. Shampine, Appl. Math. Lett. 2(4), 321 (1989). https://doi.org/10.1016/0893-9659(89)90079-7
J. Reuther, PhD thesis, Karlsruhe (2011). https://doi.org/10.5445/IR/1000023236
D. Sellmann, X.F. Zhang, S. Eggert, Phys. Rev. B 91(8), 081104 (2015). https://doi.org/10.1103/PhysRevB.91.081104
J.M. Luttinger, L. Tisza, Phys. Rev. 70(11–12), 954 (1946). https://doi.org/10.1103/PhysRev.70.954
J.M. Luttinger, Phys. Rev. 81(6), 1015 (1951). https://doi.org/10.1103/PhysRev.81.1015
J. Reuther, R. Thomale, Phys. Rev. B 83, 024402 (2011). https://doi.org/10.1103/PhysRevB.83.024402
S. Trebst, C. Hickey, Phys. Rep. 950, 1 (2022). https://doi.org/10.1016/j.physrep.2021.11.003
M. Yamada, M. Oshikawa, G. Jackeli, Phys. Rev. Lett. 121, 097201 (2018). https://doi.org/10.1103/PhysRevLett.121.097201
Acknowledgements
We thank F. L. Buessen, Y. Iqbal, and M. Scherer for discussions. Partial support from the Deutsche Forschungsgemeinschaft (DFG)—Projektnummer 277146847—SFB 1238 (project C02) is gratefully acknowledged. The numerical simulations were performed on the JURECA Booster and JUWELS cluster at the Forschungszentrum Juelich as well as the CHEOPS cluster at RRZK Cologne.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Contributions
The main part of the PFFRG code was written by Dominik Kiese (with the Julia package PFFRGSolver.jl). Lasse Gresista (with help of Dominik Kiese) performed the symmetry analysis for the considered models, adapted the PFFRG code accordingly, and performed the calculations and data analysis. Simon Trebst supervised the study. All authors were involved in the preparation of the manuscript.
Corresponding author
Appendices
Appendix A: Symmetry constraints in the asymptotic frequency parametrization
In Eq. (19), we stated the symmetry constraints of the two-particle vertex in the frequency parametrization using the three transfer frequencies s, t, and u. As was discussed in Sect. 3.2, in our implementation of the pf-FRG, we use a refined frequency parametrization [23, 24, 37], where the vertex is split into three channels \(g_c(\omega _c, v_c, v'_c)\) as defined in Eq. (22), with our choice of frequencies given in Eq. (23). We can obtain symmetry constraints for the different channels by employing the same parametrization in the spin, valley, and site indices as for the full vertex
and utilizing that the frequencies \(\omega _c, v_c, v'_c\) can be written as linear combinations of the transfer frequencies. Combining one or more symmetry constraints of the two-particle vertex, this results in symmetry constraints for the particle–particle channel
the direct particle–hole channel
and the crossed particle–hole channel
where for the fermionic frequencies \(v_c, v'_c\) we used the subscripts s, t, u instead of pp, dph, cph for brevity. Using these symmetry relations, we only have to explicitly calculate the two-particle vertex for positive values of \(\omega _c\) and \(v_c\), but have to also consider negative values for \(v'_c\). Additionally, we only have to calculate components with \({\left| v'_c\right| } < {\left| v_c\right| }\). We note again that, compared to the \({{\mathfrak {s}}}{{\mathfrak {u}}}2\) case, no constraints relating the particle–particle and crossed particle–hole channel are present.
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 \({{\mathfrak {s}}}{{\mathfrak {u}}}2\) 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 self-energy and two-particle vertex themselves already do. That the constraints are fulfilled in the initial conditions is easy to see, as for \(\varLambda \rightarrow \infty \) 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
differs from the \({{\mathfrak {s}}}{{\mathfrak {u}}}(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'}s_1l_{1'}l_1 s_{2'}s_2l_{2'}l_2 \varGamma ^\varLambda ({\bar{1}}, {\bar{2}}, {\bar{1}}', {\bar{2}}')\) and then manipulating the RHS
In step I, we applied Eq. (B.5) and transformed the sum indices \({\bar{3}}\), \({\bar{4}}\) to 3 and 4 using that the propagator is odd in frequency space. In step II, we exchanged the indices \(3 \leftrightarrow 4\) and applied the particle exchange symmetry [Eq. (57)] to the last term. This concludes the proof for the two-particle vertex. Due to the additional vertex components compared to the \({{\mathfrak {s}}}{{\mathfrak {u}}}(2)\) case, we have to repeat the proof for the self-energy, although we will closely follow Ref. [27]. To this end, we first rewrite the relations in Eqs. (57, 60) for the parametrized vertex, but using natural frequencies
which directly implies
Using these relations, we can simplify the self-energy flow equation
In step I, we simply wrote out the self-energy flow equation using the vertex parametrization. In step II, we performed the sum over \(s_2\) using \((\theta ^\mu )^2 = {\mathbbm {1}}\) and \(\mathrm {Tr}\, \theta ^\mu = 2 \delta _{\mu ,0}\) and applied Eq. (B.7) to rearrange the sum over \(\kappa \) and \(\lambda \). In step III, we perform the sum over \(l_2\) again using \(\mathrm {Tr}\, \theta ^\lambda = 2 \delta _{\lambda ,0}\) and \((\theta ^\kappa )^2 = {\mathbbm {1}}\) and the (anti)commutation relations \([\theta ^\kappa , \theta ^0] = 0\) and \(\{\theta ^\kappa , \theta ^\lambda \} = 2 \delta _{\kappa \lambda }\) (for \(\kappa , \lambda >0\)). In this form of the flow equation, it is clear that the self-energy indeed stays diagonal during the flow and, as all two-particle vertex components appearing in the last line are real and the single-scale propagator is imaginary, the self-energy is completely imaginary. That the self-energy is also odd in frequency space can easily be seen using Eq. (58) and the particle-exchange symmetry [Eq. (57)]. This concludes the proof for the self-energy.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits 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://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Gresista, L., Kiese, D. & Trebst, S. Moments and multiplets in moiré materials. Eur. Phys. J. B 95, 119 (2022). https://doi.org/10.1140/epjb/s10051-022-00373-2
Received:
Accepted:
Published:
DOI: https://doi.org/10.1140/epjb/s10051-022-00373-2