# Network of chiral one-dimensional channels and localized states emerging in
a moiré system

Jeyong Park  Institute for Theoretical Physics, University of Cologne, 50937
Cologne, Germany  Lasse Gresista  Institute for Theoretical Physics,
University of Cologne, 50937 Cologne, Germany  Simon Trebst  Institute for
Theoretical Physics, University of Cologne, 50937 Cologne, Germany  Achim
Rosch  Institute for Theoretical Physics, University of Cologne, 50937
Cologne, Germany  Jinhong Park  Institute for Theoretical Physics, University
of Cologne, 50937 Cologne, Germany

(August 30, 2024)

###### Abstract

Moiré systems provide a highly tunable platform for engineering band
structures and exotic correlated phases. Here, we theoretically study a model
for a single layer of graphene subject to a smooth moiré electrostatic
potential, induced by an insulating substrate layer. For sufficiently large
moiré unit cells, we find that ultra-flat bands coexist with a triangular
network of chiral one-dimensional (1D) channels. These channels mediate an
effective interaction between localized modes with spin-, orbital- and valley
degrees of freedom emerging from the flat bands. The form of the interaction
reflects the chiralilty and 1D nature of the network. We study this
interacting model within an $SU(4)$ mean-field theory, semi-classical Monte-
Carlo simulations, and an $SU(4)$ spin-wave theory, focusing on commensurate
order stabilized by local two-site and chiral three-site interactions. By
tuning a gate voltage, one can trigger a non-coplanar phase characterized by a
peculiar coexistence of three different types of order: ferromagnetic spin
order in one valley, non-coplanar chiral spin order in the other valley, and
120∘ order in the remaining spin and valley-mixed degrees of freedom. Quantum
and classical fluctuations have qualitatively different effects on the
observed phases and can, for example, create a finite spin-chirality purely
via fluctuation effects.

##  I Introduction

Stacking a two-dimensional van der Waals material on top of other van der
Waals materials (with or without a relative twist) defines a class of quantum
material known as moiré materials [1, 2]. Due to their highly tunable
experimental knobs for engineering band structures, thereby facilitating the
emergence of correlated phases [3, 4, 5, 6, 7], such moiré materials have
recently met with tremendous interest. A prototypical example is twisted
bilayer graphene (TBG) [8, 9, 10], where two sheets of graphene are stacked
with a relative twist. At twist angles $\sim 1.1^{\circ}$, the so-called
’magic angle’, flat bands emerge near the charge neutrality point [8, 9, 10,
11], which amplifies the effect of interaction to exhibit various correlated
phases [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28].
Besides TBG, a wealth of different types of exotic bands and interaction
effects have been discovered in multilayer moiré systems [1, 2, 29, 30, 31,
32, 33, 34, 35].

In this manuscript, we address one of the simplest models of a moiré system: a
single layer of graphene subject to a moiré potential induced by a substrate
layer. Despite its simplicity, it shows – even without fine tuning –
remarkably rich physics. For sufficiently large moiré unit cells two kinds of
moiré bands emerge: one-dimensional chiral channels (1DCCs) and ultra-flat
bands. Along lines where the gap arising from the moiré potential changes
sign, a network of topologically protected 1DCCs is developed, as depicted in
Fig. 1. At the same time, an extra set of localized modes emerges at the
junction where six 1DCCs join (red dots in Fig. 1). These modes only hybridize
weakly with the 1DCCs and with the neighboring localized modes giving rise to
ultra-flat bands. This coexistence of localized modes and propagating 1DCCs
and the resulting peculiar interaction physics are the main results of this
paper.

The emergence of a network of 1D chiral channels in moiré systems has been
previously discussed [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. In an early
study, San-Jose and Prada [36] pointed out that a network of topologically
protected 1D helical channels forms in TBG subject to an out-of-plane electric
field, see also Refs. [37, 45]. Experimentally, signatures of these 1D
channels have been observed in transport [39, 41, 42] and scanning tunneling
spectroscopy [40]. In contrast to our model, such systems do not exhibit the
coexistence of flat bands and 1D channels. Moreover, a coexistence of
propagating two-dimensional Dirac dispersing bands and flat bands has been
reported in mirror symmetric twisted trilayer graphene [30, 31]. In this
setting, Ramires and Lado discussed heavy fermion physics, emerging from the
interaction of localized and propagating modes [34]. From a more general point
of view, the emergence of localized and propagating bands in moiré systems has
been investigated in Ref. [47] using concepts of quantum chaos. Generic bands
tend not to be flat due to localization in momentum space, but these arguments
cannot be applied to the bands discussed in our paper arising from the
specific real-space structure of the moiré potential.

Figure 1:  Network model consisting of one-dimensional chiral channels (blue
arrows for the $K$ valley, black for the $K^{\prime}$ valley) and localized
states (red circles). Upper (Lower) right: the diagram for the two
(three)-spin interaction.  Figure 2: Moiré potential model. (a) Moiré pattern
from a layer of the hexagonal lattice (blue dots) stacked on top of another
layer of an hexagonal lattice (red). Due to a relative rotation of the layers,
different local stacking patterns (AA, AB, BA) occur (circles). The vectors
${\bm{r}}_{AB}=-{\bm{r}}_{BA}=(\frac{\sqrt{3}L}{6},\frac{L}{2})$ connect
regions of AA and AB stacking. (b) The uniform and staggered potential,
$V_{0}({\bm{r}})$ and $V_{s}({\bm{r}})$, in real space. (c) Band structure for
the effective Hamiltonian Eq. (1) in the moiré Brillouin zone. Flat bands are
classified by different representations (red: $\rho_{3}$, green: $\rho_{4}$,
orange: $\rho_{6}$) of the symmetry group $\textrm{Dic}_{3}$. (d) Density plot
of the energy for one of the 1D dispersing bands and the Fermi surfaces at
three different chemical potentials (red 0meV, blue 20meV, green 40meV). The
arrows indicate the three different propagating directions. (e) Bloch wave
functions $|\Psi_{n,{\bm{k}}}({\bm{r}})|^{2}$ of the green flat band at the
$\Gamma$ point. (f) Bloch wave functions of 1D propagating bands at momenta
indicated by the arrows in (d). Parameters: $u_{AA}=u_{BB}=0$ and
$u_{AB}=u_{BA}=\frac{8\pi v}{\sqrt{3}L}$.

##  II Model

We consider a single layer of graphene on top of some insulating substrate
which shares the hexagonal structure with graphene but has either a slightly
different lattice constant or is rotated by a small twist angle. As the
substrate is gapped, it mainly affects graphene via electrostatic potential
terms. Thus, at low-energies, the spinless single-particle Hamiltonian is
approximated by

$\displaystyle
H_{\textrm{eff}}=-iv(\partial_{x}s^{x}\tau^{z}+\partial_{y}s^{y}\tau^{0})+V_{0}({\bm{r}})\mathds{1}+V_{s}({\bm{r}})s^{z}.$
(1)

Here the Pauli matrices $\tau^{i}$ and $s^{i}$ act on the valley and
sublattice space, respectively, and $v$ is the graphene Fermi velocity. The
staggered term $V_{s}=(V_{A}-V_{B})/2$ describes the potential difference
between the $A$ and $B$ sublattice, $V_{A}$ and $V_{B}$, and a constant
$V_{s}$ opens a mass gap in the Dirac spectrum. The magnitude of $V_{s}$ has a
maximum in regions of the moiré lattice where the atoms of different
sublattices stack on the top of each other, i.e., AB or BA stacking as shown
in Fig. 2(a). The uniform potential $V_{0}$ is given by
$V_{0}=(V_{A}+V_{B})/2$. Due to the smoothness of moiré structures, we can
focus on the lowest Fourier components of the potentials. Denoting the six
smallest reciprocal lattice vectors of the moiré structure by ${\bm{G}}_{i}$,
$i=1,\dots,6$, with $|{\bm{G}}_{i}|\equiv G=\frac{4\pi}{\sqrt{3}L}$, we obtain

$\displaystyle V_{\beta}({\bm{r}})$
$\displaystyle=\sum\limits_{\beta^{\prime}=A,B}{\sum\limits_{i=1}^{6}{u_{\beta\beta^{\prime}}e^{i{\bm{G}}_{i}\cdot({\bm{r}}-{\bm{r}}_{\beta\beta^{\prime}})}}}\,\,$
(2)

with sublattice $\beta=A,B$, the size of moiré unit cell $L$, and
${\bm{r}}_{AA}={\bm{r}}_{BB}=0$, ${\bm{r}}_{AB}=-{\bm{r}}_{BA}$, see Fig.
2(a). We consider a hexagonal substrate with equivalent $A$ and $B$
sublattices such that $u_{AA}=u_{BB}$ and $u_{AB}=u_{BA}$. In this case, the
amplitudes of $V_{s}$ and $V_{0}$ are given by $u_{s}=u_{AB}$ and
$u_{0}=-2u_{AA}+u_{AB}$.

As shown in Fig. 2(b), $V_{s}$ vanishes along straight lines and thus changes
its sign across those lines. At the same time, $V_{0}$ has minima at high-
symmetry points (red dots) in the center of the moiré unit cell where the
lines cross. These two regions lead to two very different types of bands: 1D
dispersing bands and ultra-flat bands, see Fig. 2(c), computed by
diagonalizing Eq. (1) in momentum space.

###  II.1 One-dimensional chiral channels

The sign change of the mass term $V_{s}({\bm{r}})$ induces a 1DCC, propagating
along the straight lines in Fig. 2(b) with the full speed of the graphene
Fermi velocity. 1DCCs emerging from the $K$ and $K^{\prime}$ valley move in
opposite directions, as depicted by black and blue arrows in Fig. 1. Fig. 2(f)
shows the Bloch wave function of the propagating bands which perfectly tracks
the straight lines in Fig. 2(b). Surprisingly, the wave functions show almost
no modulation at their crossing points.

###  II.2 Localized states

The flat bands in Fig. 2(c) have their origins in states localized close to
the red dots in Fig. 2(b), where $V_{0}$ has a minimum whereas $V_{s}$ is
highly suppressed. A sufficiently strong moiré potential, $u_{0}\gg v/L$,
renders the states localized in real space as shown in Fig. 2(e). These
localized modes hybridize only weakly with the 1DCCs and neighboring localized
modes, leading to ultra-flat bands. The localized states with fixed valley
index can be classified by the dicyclic symmetry group $\text{Dic}_{3}$. From
$\text{Dic}_{3}$, one obtains three different types of localized states,
labeled by two one-dimensional representations $\rho_{3}$, $\rho_{4}$ and a
two-dimensional irreducible representation $\rho_{6}$, see supplement [48].

###  II.3 Network model

Combining localized and propagating states, we obtain the network model
depicted in Fig. 1. The kinetic Hamiltonian for the 1DCCs is given by

$\displaystyle H_{\rm kin}=-iv\sum_{\alpha,n,i,\sigma}\alpha\int
d\rho\Psi^{\dagger}_{n,i,\alpha,\sigma}(\rho)\partial_{\rho}\Psi_{n,i,\alpha,\sigma}(\rho).$
(3)

The operator $\Psi^{\dagger}_{n,i,\alpha,\sigma}(\rho)$ creates an electron
with spin $\sigma=\uparrow/\downarrow$ and valley $\alpha=\pm=K/K^{\prime}$ in
a 1DCC propagating along the lattice vector ${\bm{a}}_{i}$ with $i=1,2,3$; the
center of the corresponding wave packet is located at
${\bm{R}}_{n,i,\rho}=n{\bm{a}}_{i+1}+\rho{\bm{a}}_{i}$ with integer $n$. We
denote the location of crossing points of 1DCCs by ${\bm{R}}_{m}$ and define
$\rho_{n,i,m}$ as the solution of ${\bm{R}}_{m}={\bm{R}}_{n,i,\rho_{n,i,m}}$.
In these notations, the inter-channel tunneling $H_{w}$ and the coupling of
1DCCs to localized states, $H_{\lambda}$, are given by

$\displaystyle H_{w}=$ $\displaystyle\sum_{\text{crossing at
}{\bm{R}}_{m}}\hat{w}_{ii^{\prime}}\Psi^{\dagger}_{n,i,\alpha,\sigma}(\rho_{n,i,m})\Psi_{n^{\prime},i^{\prime},\alpha,\sigma}(\rho_{n^{\prime},i^{\prime},m}),$
$\displaystyle H_{\lambda}=$ $\displaystyle\sum_{\text{crossing at
}{\bm{R}}_{m}}\hat{\lambda}_{ij}\,d^{\dagger}_{m,\alpha,\sigma,j}\Psi_{n,i,\alpha,\sigma}(\rho_{n,i,m})+h.c.\
.$ (4)

We sum over all channels which cross at ${\bm{R}}_{m}$.
$d^{\dagger}_{m,\alpha,\sigma,j}$ creates localized electronic states where
$j$ denotes an extra orbital index if the localized states belong to the
$\rho_{6}$ representation. The form of the matrices $\hat{w}_{ii^{\prime}}$
and $\hat{\lambda}_{ij}$ is entirely determined by the symmetries of the
system and the representation of Dic3 of the localized states. Eqs. (3) and
(4) describe the bandstructure with high precision after fitting the amplitude
of $\hat{w}$ and $\hat{\lambda}$ and the energy of the localized states, see
supplement [48].

###  II.4 Local interaction

Since the flat bands are highly localized, there will be a Coulomb blockade
for adding electrons to the localized sites, described by

$\displaystyle
H_{U}=U\sum_{{\bm{R}}_{m}}\sum_{\xi\neq\xi^{\prime}}n_{m,\xi}n_{m,\xi^{\prime}},$
(5)

where $\xi=\\{\alpha,\sigma\\}$ includes all local quantum numbers, i.e.,
valley $\alpha$ and spin $\sigma$ (and an extra orbital quantum number for the
$\rho_{6}$ representation). Since for large moiré unit cells the dominant
contribution comes from the long-ranged part of the Coulomb interaction which
is only sensitive to charge, $H_{U}$ is approximately $SU(4)$ (or $SU(8)$ for
$\rho_{6}$) invariant. Using $U\sim\frac{e^{2}}{4\pi\epsilon_{0}d_{\rm loc}}$
with $d_{\rm loc}\approx 6.3$nm, we estimate $U\approx 230$meV for the
parameters of Fig. 2(e) which is more than an order of magnitude larger than
the hybridization of impurity levels, $\lambda$.

As $U\gg\lambda$, the system maps to a (generalized) Kondo lattice model,
where local degrees of freedom couple only via the network of 1DCCs. For a
localized state in the $\rho_{4}$ representation, one obtains an effective
SU(4) symmetric coupling

$\displaystyle H_{J}\approx JL\\!\\!\\!\sum_{\textrm{crossing at
}{\bm{R}}_{m}}\\!\\!\\!\Gamma^{\ell}({\bm{R}}_{m})\cdot\tilde{\Psi}^{\dagger}_{\xi}({\bm{R}}_{m})\gamma^{\ell}_{\xi\xi^{\prime}}\tilde{\Psi}_{\xi^{\prime}}({\bm{R}}_{m}).$
(6)

Here $\gamma^{\ell}$, $\ell=1,\dots,15$, are the 4$\times$4 generators of
$SU(4)$ acting on a linear combination of the three 1DCCs resulting from the
hybridization matrix $\hat{\lambda}$ [48],
$\tilde{\Psi}_{\xi}({\bm{R}}_{m})=\frac{1}{\sqrt{3}}\sum_{i=1,2,3}(-1)^{i}\Psi_{n_{i},i,\xi}(\rho_{n_{i},i,m})$.
$\Gamma^{\ell}({\bm{R}}_{m})=\sum_{\xi,\xi^{\prime}}d^{\dagger}_{m,\xi}\gamma^{\ell}_{\xi\xi^{\prime}}d_{m,\xi^{\prime}}$
describes the local $SU(4)$ degree of freedom and $J\sim\lambda^{2}/U$ is the
Kondo coupling.

Figure 3: Mean-field phase diagrams with the order parameter,
$\chi^{+}-\chi^{-}$, for small $J_{3},J_{2}^{\prime}$, and $\varphi$, defined
modulo $2\pi/3$ (cf. Eq. (III)). All phases shown in the figure show a
$120^{\circ}$ order in the ${\bm{\mu}}$ valley degrees of freedom. They differ
by their spin order which are distinguished by
$\chi^{\pm}\equiv\langle{\bm{\sigma}}P^{\pm}\rangle_{A}\cdot\langle{\bm{\sigma}}P^{\pm}\rangle_{B}{\bm{\times}}\langle{\bm{\sigma}}P^{\pm}\rangle_{C}$,
the valley-projected chirality, defined on a triangle (cf. the inset). Three
different types of spin-order occur: spin-ferromagnic order (green regions), a
chiral spin-order (red or blue regions) in one of the two valleys, and an in-
plane 120∘ spin order in one of the two valleys (white dotted line). $J_{2}$
is set to 1. Panel (a): $(\varphi,J_{3})$ space with $J_{2}^{\prime}=0$, panel
(b): $(\varphi,J_{2}^{\prime})$ with $J_{3}=0$, panel (c):
$(J_{3},J_{2}^{\prime})$ with $\varphi=0$. The yellow line in panel (b) shows
schematically how parameters change as a function of a gate voltage for
$k_{F}L$ close to $-\frac{\pi}{2}$. While the solid white lines represent
second order transitions, the dashed white line a first order transition where
the type of $120^{\circ}$ order changes from left- to right circulating on a
triangle $\triangleright$ (opening angle $\theta=0,\pi$, respectively, see
[48]).

The interaction between the localized states is mediated by the network of
1DCCs. The resulting RKKY interaction is obtained from a perturbation theory
both in $J$ and the inter-channel tunneling $w$. From a standard RKKY diagram
(see Fig. 1) to order $J^{2}w^{0}$, we obtain the two-spin interaction term

$\displaystyle H_{\textrm{2s}}$ $\displaystyle=-\frac{J^{2}L^{2}}{12\pi
v}\sum_{(m_{1}\to
m_{2})_{c}}\bigg{(}\frac{e^{2i\bm{k_{F}}^{i}\cdot{\bm{\rho}}_{m_{1}m_{2}}}}{|{\bm{\rho}}_{m_{1}m_{2}}|}$
$\displaystyle\times\left(1+{\bm{\sigma}}_{m_{1}}\cdot{\bm{\sigma}}_{m_{2}}\right)\tau^{-}_{m_{1}}\tau^{+}_{m_{2}}+h.c.\bigg{)}.$
(7)

${\bm{\sigma}}_{m}\equiv{\bm{\sigma}}({\bm{R}}_{m})$ and
${\bm{\tau}}_{m}\equiv{\bm{\tau}}({\bm{R}}_{m})$ are Pauli matrices acting on
spin and valley at ${\bm{R}}_{m}$, and $\tau^{\pm}\equiv\tau^{x}\pm
i\tau^{y}$. The summation $(m_{1}\to m_{2})_{c}$ runs only over localized
states connected by the same 1DCC, separated by
${\bm{\rho}}_{m_{1}m_{2}}\equiv{\bm{R}}_{m_{2}}-{\bm{R}}_{m_{1}}$ parallel to
the Fermi velocity of the $+$ valley channels. Importantly, the RKKY term
necessarily requires two valley flip processes,
$\tau^{-}_{m_{1}}\tau^{+}_{m_{2}}$: both a valley $+$ and $-$ channel running
in opposite directions are needed to form a closed loop connecting two sites,
see Fig. 1. This process breaks the $SU(4)$ symmetry.

Closed loops can also be formed by triangles in Fig. 1, remarkably, inducing a
chiral interaction to order $J^{3}$. From the diagram depicted in Fig. 1, we
obtain

$\displaystyle H_{3s}$ $\displaystyle=\frac{8J^{3}L^{3}}{27\sqrt{3}\pi
v^{2}}\sum_{p=\triangleright/\triangleleft,(m_{1},m_{2},m_{3})_{p}}\frac{\textrm{cos}(3k_{F}|{\bm{\rho}}_{m_{1},m_{2}}|)}{|{\bm{\rho}}_{m_{1},m_{2}}|}$
$\displaystyle\times p\
{\bm{\sigma}}_{m_{1}}\cdot({\bm{\sigma}}_{m_{2}}{\bm{\times}}{\bm{\sigma}}_{m_{3}})\Big{(}\prod_{i=1}^{3}P_{m_{i}}^{+}-\prod_{i=1}^{3}P_{m_{i}}^{-}\Big{)}.$
(8)

The summation runs over the right- and left-oriented triangles,
$p=\triangleright/\triangleleft=\pm 1$, in Fig. 1, where
$(m_{1},m_{2},m_{3})_{p}$ denotes the three sites of each triangle (in anti-
clockwise order). $k_{F}$ is the Fermi momentum of 1DCCs. We defined the
projector $P_{m}^{\pm}\equiv\frac{(\tau^{0}_{m}\pm\tau^{z}_{m})}{2}$ on valley
$\pm$ at ${\bm{R}}_{m}$. The chiral spin-interaction,
${\bm{\sigma}}_{m_{1}}\cdot({\bm{\sigma}}_{m_{2}}{\bm{\times}}{\bm{\sigma}}_{m_{3}})$,
is induced by the chiral motion of the 1DCCs within each triangle (even in the
absence of spin-orbit interaction). The direction of the chiral currents
determines the sign of the chiral interaction which changes when moving from
$\triangleright$ to $\triangleleft$ or from valley $+$ to $-$.

There is also a non-chiral contribution from the same diagram and from a
similar diagram to order $J^{2}w$,

$\displaystyle H^{\prime}_{2s}=\frac{J^{2}L^{3}}{27\sqrt{3}\pi
v^{2}}\sum_{p=\triangleright/\triangleleft,(m_{1},m_{2},m_{3})_{p}}\frac{\sin(3k_{F}|{\bm{\rho}}_{m_{1},m_{2}}|)}{|{\bm{\rho}}_{m_{1},m_{2}}|}$
$\displaystyle\times(1+{\bm{\sigma}}_{m_{1}}\cdot{\bm{\sigma}}_{m_{2}})\Big{(}J\tau^{z}_{m_{3}}(\tau^{z}_{m_{1}}+\tau^{z}_{m_{2}})+w(1+\tau^{z}_{m_{1}}\tau^{z}_{m_{2}})\Big{)}$
$\displaystyle+\text{permutations}.$ (9)

We sum over the 6 permutations for renaming $m_{1}$, $m_{2}$ and $m_{3}$.

##  III Mean-field Phase diagram

To study the interplay of Eqs. (II.4)-(II.4), we consider a simplified
Hamiltonian which contains only nearest neighbor interactions,
$H_{\rm{sv}}=H_{2}+H_{3}+H_{2^{\prime}}$, with

$\displaystyle H_{2}=$ $\displaystyle J_{2}\sum_{\langle m_{1}\to
m_{2}\rangle_{c}}(1+{\bm{\sigma}}_{m_{1}}\cdot{\bm{\sigma}}_{m_{2}})(e^{i\varphi}\tau^{+}_{m_{1}}\tau^{-}_{m_{2}}+h.c.)$
$\displaystyle H_{3}=$ $\displaystyle
J_{3}\\!\\!\\!\sum_{p=\triangleright/\triangleleft,(m_{1},m_{2},m_{3})_{p}}\\!\\!\\!\\!p\Big{(}\prod_{i=1}^{3}P_{m_{i}}^{+}-\prod_{i=1}^{3}P_{m_{i}}^{-}\Big{)}$
$\displaystyle\hskip
113.81102pt\times{\bm{\sigma}}_{m_{1}}\cdot({\bm{\sigma}}_{m_{2}}{\bm{\times}}{\bm{\sigma}}_{m_{3}})$
$\displaystyle H_{2^{\prime}}=$ $\displaystyle J_{2}^{\prime}\sum_{\langle
m_{1}\to
m_{2}\rangle_{c}}(1+{\bm{\sigma}}_{m_{1}}\cdot{\bm{\sigma}}_{m_{2}})(1+\tau^{z}_{m_{1}}\tau^{z}_{m_{2}})\,\,.$
(10)

Here $J_{2}>0$ is the largest coupling constant with $\varphi=2k_{F}L+\pi$,
while $J_{2}^{\prime}\sim\sin(3k_{F}L)$ and $J_{3}\sim\cos(3k_{F}L)$. From
$H_{2s}^{\prime}$, Eq. (II.4), we take, for simplicity, only the term $\sim
J^{2}w$ into account (assuming $w>J$) but we checked that the $J^{3}$
contribution to $H_{2s}^{\prime}$ does not lead to qualitative changes. The
continuous symmetries of $H_{\rm{sv}}$ are $U(1)\times SU(2)\times SU(2)$
generated by $\tau^{z}$, $P^{+}{\bm{\sigma}}$ and $P^{-}{\bm{\sigma}}$.
Remarkably, one can rotate the spin-orientation of the two valleys
independently.

Assuming that the localized states are filled with one electron, the states on
the $SU(4)$ space are spanned by a 4-component complex vector. In this basis,
we solve the self-consistent mean-field equations at $T=0$ iteratively. We
find that either a one- or a three-sublattice solution has the lowest energy.
As $J_{2}$ is the largest term, we first analyze the case
$J_{3}=J^{\prime}_{2}=0$. The parameter $\varphi$ in Eq. (III) can be viewed
as an Aharonov-Bohm phase arising from a staggered magnetic flux. As
$3\varphi$ is the total phase along a triangular loop, one can always ‘gauge
away’ changes of $\varphi$ by $\frac{2\pi}{3}$ using $\tau_{z}$ rotations by
$0$, $2\pi/3$, $4\pi/3$ on the A, B, C sublattices. For $\varphi=0$, we obtain
a variant of the Kugel-Khomskii model [49]

$\displaystyle H_{0}$ $\displaystyle=2J_{2}\sum_{\langle m_{1}\to
m_{2}\rangle_{c}}(1+{\bm{\sigma}}_{m_{1}}\cdot{\bm{\sigma}}_{m_{2}})(\tau^{x}_{m_{1}}\tau^{x}_{m_{2}}+\tau^{y}_{m_{1}}\tau^{y}_{m_{2}})$
$\displaystyle=2J_{2}\sum_{\langle m_{1}\to
m_{2}\rangle_{c}}{\bm{\mu}}^{1}_{m_{1}}\cdot{\bm{\mu}}^{1}_{m_{2}}+{\bm{\mu}}^{2}_{m_{1}}\cdot{\bm{\mu}}^{2}_{m_{2}}$
(11)

with four component vectors given by
${\bm{\mu}}^{1}_{m}=(\tau^{x}_{m},\tau^{y}_{m}\sigma^{x}_{m},\tau^{y}_{m}\sigma^{y}_{m},\tau^{y}_{m}\sigma^{z}_{m})$
and
${\bm{\mu}}^{2}_{m}=(\tau^{y}_{m},\tau^{x}_{m}\sigma^{x}_{m},\tau^{x}_{m}\sigma^{y}_{m},\tau^{x}_{m}\sigma^{z}_{m})$.
The ground states have a three-site unit-cell where the vectors
$\langle{\bm{\mu}}^{n}_{m}\rangle$, $n=1,2$, have the norm $1$, and show 120∘
order such that
$\langle{\bm{\mu}}^{n}_{m_{1}}\rangle\cdot\langle{\bm{\mu}}^{n}_{m_{2}}\rangle=\cos(2\pi/3)=-\frac{1}{2}$
for neighboring sites. Note that this specific type of 120∘ order is realized
with 4-component vectors. Surprisingly, the above described 120∘ order has an
extra degree of freedom that is revealed by the magnetization vectors $\langle
P^{\pm}{\bm{\sigma}}\rangle$ in the two valleys. These vectors have length
$1/2$ in the ground-state manifold. In one of the two valleys, the
magnetization is always ferromagnetic, but in the other valley a non-coplanar
spin configuration is possible, leading to a finite staggered chirality
$\chi^{\pm}$ with
$\chi^{\pm}=\langle{\bm{\sigma}}_{m_{1}}P^{\pm}_{m_{1}}\cdot({\bm{\sigma}}_{m_{2}}P^{\pm}_{m_{2}}{\bm{\times}}{\bm{\sigma}}_{m_{3}}P^{\pm}_{m_{3}})\rangle$.
In the supplementary material [48], we describe how the mean-field solution
can be parameterized by a continuous angle $\theta$ and a discrete variable
$\pm$, describing the opening angle of non-coplanar valley-projected spins on
the three sublattices and also which of the valley sector exhibits
ferromagnetic order.

States with an arbitrary chirality,
$-\frac{1}{8}\leq\chi^{\pm}\leq\frac{1}{8}$, are degenerate (within mean-field
theory) if only $H_{0}$, Eq. (III), is considered, see supplement [48]. Thus,
$H_{0}$ defines a highly singular point in the phase diagram and even small
perturbations can select one of the states in the ground-state manifold of
$H_{0}$. For example, for an infinitesimal $J_{3}>0$ perturbation, states are
selected which have either the minimal value $\chi^{+}=-\frac{1}{8}$ with
$\chi^{-}=0$ or the maximal value of $\chi^{-}=\frac{1}{8}$ with $\chi^{+}=0$.
Such a staggered (or uniform) chiral order has, e.g., been extensively studied
in the spin-$1/2$ [50, 51] or the half-filled Hubbard model [52, 53, 54, 55]
on the triangular lattice.

In contrast, the perturbation by a finite $\varphi$ stabilizes a phase where
$\langle P^{\pm}\sigma\rangle$ orders ferromagnetically for both valleys,
$\chi^{+}=\chi^{-}=0$, while $\langle{\bm{\mu}}^{n}\rangle$ displays a
coplanar 120∘ ordered phase. The presence of both $\varphi$ and $J_{3}$ leads
to the phase diagram of Fig. 3(a). A finite $J_{2}^{\prime}>0$, however,
suppresses such ferromagnetic configuration, selecting a state where $\langle
P^{\pm}\sigma\rangle$ is non-collinear but coplanar, forming a 120∘ order in
either $\langle P^{+}\sigma\rangle$ or $\langle P^{-}\sigma\rangle$ on top of
the 120∘ order in $\langle{\bm{\mu}}^{n}\rangle$. The resulting phase diagrams
are shown in Fig. 3(b) and (c).

##  IV Classical fluctuations

The mean-field theory discussed above, ignores the effect of both quantum and
classical fluctuations. To capture fluctuation effects, we have (i) performed
an $SU(4)$ spin-wave calculation (or, more precisely, spin-valley-wave
calculation) both in the classical and quantum regime. Details of the $SU(4)$
spin-wave theory are given in the supplementary material [48]. Furthermore, we
have (ii) calculated finite temperature properties of the semi-classical
version of our $SU(4)$ model using Monte Carlo calculations.

A semi-classical variant of our $SU(4)$ model can formally be obtained by
making a product ansatz for the wavefunction,
$|\Psi\rangle=\prod_{m}|\Psi_{m}\rangle$, where $|\Psi_{m}\rangle$ is a
single-site 4-component normalized wave function with an arbitrary phase per
site. A semi-classical state for a system of size $N\times N$ is thus
parameterized by $(4\cdot 2-2)N^{2}$ real numbers. At $T=0$, this semi-
classical model reproduces the mean-field results discussed in Sec. III.
Thermal expectation values at a finite temperature $T=1/\beta$ can be
approximately calculated by sampling the space of product-state wavefunctions
according to the Boltzmann distribution
$\sim\exp(-\beta\langle\psi|H|\psi\rangle)$ [56, 57] using a standard Markov
chain Monte Carlo algorithm [58]. Employing local Metropolis updates, a
typical Monte Carlo run consists of $N_{m}=1\cdot 10^{6}$ thermalization
sweeps followed by $N_{m}=4\cdot 10^{6}$ measurement sweeps, or up to
$N_{m}=10^{7}$ sweeps close to the transition temperature. We use linear
lattice sizes of up to $N=72$ with periodic boundary conditions. Additional
details on the simulations are provided in the supplementary material [48].

Figure 4: Thermodyamics and thermal order-by-disorder transition. Shown are
Monte Carlo results for the $J_{2}$ model, Eq. (III), in the semi-classical
approximation for different linear system sizes $N=12,24,36$ and $72$.
Numerical errors are smaller than the size of the symbols. Panel (a): The
specific heat shows a pronounced peak indicating a thermal phase transition
which gets sharper upon increasing $N$. Panel (b): At finite temperature $T$
and inside the ordered phase a finite spin-chirality develops. The
$T$-dependence at low temperature is singular and approximately proportional
to $T\log 1/\sqrt{T}$, see inset. Inset: Fit to the analytical result (12) at
low $T$ (dashed orange line). Panel (c): By an order-by-disorder mechanism,
the system selects a state with ferromagnetic spin order at low $T$. At finite
$T$ the order parameter is suppressed by thermal fluctuations linear in $T$.
The prefactor of the linear correction increases with system size, reflecting
the suppression of long-ranged order by thermal fluctuations, consistent with
the Mermin-Wagner theorem [59].  Figure 5: Phase diagram from semi-classical
Monte Carlo calculations. Shown is Monte Carlo data for the specific heat (a),
spin-chirality (b), and valley-magnetization (c), obtained from simulations of
a model with $\varphi=J_{3}=0$ as function of temperature and $J_{2}^{\prime}$
for a fixed linear system size $N=36$. Panel (a): Specific heat. At high
temperatures one obtains a paramagnetic phase while at low $T$ one finds two
ordered phases. All phase transitions appear to be of first-order type. For
$J_{2}^{\prime}<0$ one obtains 120∘ valley order coexisting with spin-
ferromagnetic order, while for $J_{2}^{\prime}>0$ a state with 120∘ valley
order and coplanar 120∘ spin order in one of the two valleys is realized. The
transition temperature between the two phases follows the analytically
estimated slope $T_{c}\approx 8.21J_{2}^{\prime}$ (white dashed line in panel
(b)). Panel (b): The spin-ferromagnetic order supports a finite spin-chirality
$|\chi^{+}-\chi^{-}|$ at finite $T$. Panel (c): As the spins order coplanar in
one valley but remain ferromagnetic in the other, $\tau_{z}$ develops a finite
expectation value.

As the mean-field ground state of the $J_{2}$-only model $H_{0}$, Eq. (III),
is degenerate, we focus our discussions on fluctuation effects around this
state. The specific heat of the semi-classical model, Fig. 4(a), shows a sharp
peak indicating a finite-temperature phase transition. The numerical data is
both consistent with a weak first-order or a second-order transition, see
supplementary material [48] which also discusses energy distributions at
criticality. We analyze two types of order parameters, the spin-chirality,
Fig. 4(b), and the valley-projected ferromagnetic order, Fig. 4(c), which show
very different finite-size and temperature behavior as discussed below.

For $T\to 0$, the spin-chirality, Fig. 4(b), vanishes while the ferromagnetic
magnetization in both valley sectors takes the value $1/2$, Fig. 4(c). This
shows that thermal fluctuations select the spin-ferromagnetic states,
$\theta=0,\pi$ from the ground-state manifold. At the same time, the valley,
more precisely ${\bm{\mu}}^{1,2}$, exhibits 120∘ order (not shown). This
‘order-by-disorder’ selection [60] of the classical ground state in the limit
$T\to 0$ is also found within our $SU(4)$ spin wave calculation, see Fig. 6(b)
below and supplement [48]: a fluctuation correction to the free energy linear
in $T$ selects the ferromagnetic state.

At finite $T$, the ferromagnetic order parameter shown in Fig. 4(c) is
suppressed linearly in $T$. The prefactor of this suppression increases with
system size $N$. This is explained by an order-parameter suppression $\propto
T\ln N$, well known from the Mermin-Wagner theorem [59] in two spatial
dimensions. Thus, there is nominally no long-ranged spin-order in the
thermodynamic limit.

A remarkable result is that in the spin-ferromagnetic state the spin-chirality
becomes finite at finite $T$, Fig. 4(b), showing a highly singular $T$
dependence which is almost independent on system size $N$. In the
supplementary material [48], we use an $SU(4)$ spin-wave calculation to
compute $\langle\hat{\chi}^{+}-\hat{\chi}^{-}\rangle$. The $SU(4)$ spin wave
theory is formally derived using a $1/M$ expansion, where $M$ are the number
of local bosons, $\sum_{\xi=1}^{4}b^{\dagger}_{m,\xi}b_{m,\xi}=M$, used to
describe the local $SU(4)$ degree of freedom, see [48]. For $M\to\infty$ one
recovers mean-field and spin-waves are computed to leading order in $1/M$,
where $M$ is set to its physical value, $M=1$ at the end of the calculation,
corresponding to one localized electron per site. In the classical limit, at
low-$T$ deviations from mean-field are small, which allows to make
quantitative predictions based on spin wave theory.

The naive spin-wave calculation in the classical limit predicts a divergent
result reflecting the ground-state degeneracy of the $T=0$ state. This
degeneracy is lifted by the order-by-disorder mechanism discussed above which
provides a mass linear in $T$ to the chirality-mode. Taking this higher-order
(in $1/M$ [48]) effect into account we obtain

$\displaystyle\langle\hat{\chi}^{+}-\hat{\chi}^{-}\rangle\approx\pm
0.22\,\frac{T}{J_{2}}\ln\left[\sqrt{T_{0}/T}\right]\,,$ (12)

in perfect agreement with the numerical data, see inset of Fig. 4(b). The
prefactor is fixed by our analytical results, see supplement [48], and the
only fitting parameter is $T_{0}$. We expect that such a non-analytic $T$
dependence is generic for classical systems with a degenerate ground-state
manifold where a ground state of the manifold is selected by thermal
fluctuations. Thus some ‘pseudo Goldstone modes’ obtain masses linear in $T$,
leading to non-analytic $T\log 1/T$ corrections in spatial dimension $d=2$ or
a $c_{1}T+c_{2}T^{3/2}$ correction in spatial dimension $d=3$ for observables
coupling to the mode, see supplementary material [48].

The sign in Eq. (12) is related to the spontaneous breaking of the
$\mathbb{Z}_{2}$ symmetry, $e^{i\pi\tau_{x}/2}=i\tau_{x}$, which maps
$\chi^{+}$ to $\chi^{-}$. Therefore,
$\langle\hat{\chi}^{+}-\hat{\chi}^{-}\rangle$ can be used, at $T>0$, as an
Ising order parameter of this symmetry. The extremely sharp rise of
$\langle\hat{\chi}^{+}-\hat{\chi}^{-}\rangle$ at the phase transition, see
Fig. 4(b), is both consistent with an Ising phase transition,
$\langle\hat{\chi}^{+}-\hat{\chi}^{-}\rangle\sim(T_{c}-T)^{1/8}$, or a first-
order transition, see supplement [48].

In Fig. 5 we show the phase diagram of the $J_{2}-J_{2}^{\prime}$ model as a
function of temperature $T$ and coupling $J_{2}^{\prime}$. At $T=0$, this
simply reproduces the mean-field result. While for $J_{2}^{\prime}<0$ a
ferromagnetic spin-order coexists with a 120∘ valley order, one obtains a
coplanar spin order in one of the two valley sectors for $J_{2}^{\prime}>0$.
Thus the valley symmetry is spontaneously broken in this phase, leading to a
finite expectation value for $\tau_{z}$, see Fig. 5(c). Numerically, we find
that the phase transition into the spin-coplanar phase at $J_{2}^{\prime}>0$
both as a function of $T$ or $J_{2}^{\prime}$ is always of first order; an
analysis of the energy distribution is given in the supplement [48]. At low
$T$, the first-order phase transition separating the two ordered phases has a
linear slope, $T_{c}\propto J_{2}^{\prime}$. This arises because at
$J_{2}^{\prime}=0$ the spin-ferromagnetic state gains energy linear in $T$ due
to the order-by-disorder mechanism described above. This linear-in-$T$ energy
gain competes with a linear-in-$J_{2}^{\prime}$ energy gain of the spin-
coplanar phase, Fig. 6(b), which arises because $J_{2}^{\prime}$ selects at
$T=0$ one of the states from the ground-state manifold of $H_{0}$.
Analytically, we obtain from this argument $T_{c}\approx 8.21J_{2}^{\prime}$,
which quantitatively explains the numerically observed slope, as shown in Fig.
5(b). As discussed above, the finite-$T$ transition from the paramagnetic into
spin ferromagnetic phase is accompanied by a $\mathbb{Z}_{2}$ symmetry
breaking.

For all considered values of $J_{2}^{\prime}$ the specific heat shows a low
temperature saturation of $c_{v}(T\to 0)=3$, indicating that both ordered
states in Fig. 5 feature six harmonic modes [61] per site as expected for an
$SU(4)$ model locally described by 6 parameters as discussed above.

Figure 6: Order-by-disorder mechanisms. (a) Ground-state energy per site
computed from $SU(4)$ spin-wave theory [48] as function of $\theta$,
parameterizing the opening angle of valley-projected spins in the three
sublattices. By the quantum order-by-disorder mechanism, a state with the an
opening angle $\theta=\pi/2$ (i.e., the spin-coplanar 120 degree order in one
of the two valleys) is selected from the mean-field ground-state manifold. (b)
In the classical model, the free energy obtains at low $T$ a correction linear
in $T$ from thermal fluctuations. In contrast with the quantum order-by-
disorder mechanism, the thermal order-by-disorder leads to the selection of
the spin ferromagnetic order in both of the valley sectors ($\theta=0$ or
$\pi$).

##  V Quantum fluctuations

Above, we discussed the effect of thermal fluctuations and showed that at low
$T$ an $SU(4)$ spin-wave calculation in the classical regime reproduces the
main numerical findings qualitatively and quantitatively including the order-
by-disorder mechanism and non-analytic $T$ dependences arising from pseudo
Goldstone modes. While the $SU(4)$ spin-wave theory becomes exact in the
classical case for $T\to 0$, this is not the case in the quantum model, where
quantum fluctuations in the ground state can be large. $SU(4)$ spin wave
theory only becomes exact in a large $M$ limit, see supplement [48] but we
expect that qualitative features of ordered phases (in contrast to spin-liquid
phases) are well captured by this approach.

In Fig. 6(a), we show the corrections due to quantum fluctuations to the
ground-state energy of $H_{0}$, Eq. (III), as a function of the spin-opening
angle $\theta$. The state with $\theta=\pi/2$, i.e., a coplanar 120∘ order of
the spins in one of the valley sectors, is selected by quantum fluctuations.
In contrast, as discussed above, thermal fluctuations select spin-
ferromagnetic order ($\theta=0$ or $\pi$). Thus, our system is one of the rare
cases where quantum and classical fluctuations select very different types of
ground states. As we show in the supplement, this arises, technically, because
classically a state is selected where the geometric average of the excitation
energies $E_{{\bm{k}},n}$ is lowest, while quantum fluctuations select the
state with the lowest arithmetic average of all $E_{{\bm{k}},n}$. While in
most systems the two averages show the same qualitative behavior, this is not
the case in our system.

How will the quantum fluctuations modify the ground-state phase diagrams shown
in Fig. 3? The main effect of quantum fluctuations is that they break the
degeneracy of mean-field ground state of $H_{0}$. As shown in Fig. 6(a), the
ground-state energy obtains a $\theta$ dependence. An almost identical
$\theta$ dependence can be obtained in the purely classical model by adding a
$J_{2}^{\prime}$ to the Hamiltonian with $J_{2}^{\prime}\approx 0.45\,J_{2}$.
Thus, we speculate that the quantum fluctuations have a similar effect on the
phase diagram as increasing $J_{2}^{\prime}$ within mean-field theory. This
procedure is well-controlled in an $1/M$ expansion, see supplement [48]: for
large $M$, quantum corrections of order $1/M$ can be fully compensated by a
shift of $J_{2}^{\prime}$ by $-0.45\,J_{2}/M$ (up to corrections of order
$1/M^{2}$).

Thus, we expect that the main effect of quantum fluctuations will be that in
Fig. 3(b) and (c) the phase boundaries are shifted along the $y$ direction,
most likely accompanied by a rounding of the sharp kink where the phases meet.
This extra rounding would be a $1/M^{2}$ effect, which is more difficult to
calculate.

##  VI Conclusions

Our study reveals that one of the simplest feasible moiré systems, a single
layer of graphene on a substrate, can exhibit surprisingly rich physics.
Ultra-flat bands generating localized modes coexist with a network of chiral
one-dimensional channels where electrons move very fast with a speed set by
the Fermi velocity of graphene. A main advantage of such large-unit-cell
system is that one can tune the electron density by external gates.

Different types of localized modes with spin-, valley- and orbital degrees of
freedom can be realized depending on how many electrons are loaded into the
local level and the quantum numbers of the localized states, fixed by the
representation of the relevant dicyclic group. The chiral nature of the
channels connecting the localized modes gives rise to characteristic chiral-
and non-chiral interactions. We expect that a wealth of different phases with
commensurate and incommensurate spin-, valley- and orbital order can be
realized.

As an example, we studied one such model, focusing on commensurate order
stabilized by two- and three spin interactions. We use mean-field theory, an
$SU(4)$ spin-wave theory both in the classical and quantum regime, and Monte
Carlo simulations of a semi-classical model. By tuning gate voltages one can
control $k_{F}$ and thus the effective interactions. One can, for example,
tune parameters along the yellow line shown in Fig. 3b. This triggers a
transition from a coplanar phase with ferromagnetic spin and 120∘ valley order
into a non-coplanar phase characterized by a peculiar coexistence of three
different types of order: ferromagnetic spin order in one valley, non-coplanar
chiral spin order in the other valley, and 120∘ order in remaining spin and
valley-mixed degrees of freedom.

The peculiar form of the mean-field phase diagrams, where tiny perturbations
can profoundly change the ground state, is governed by the proximity to a
variant of the Kugel-Khomskii model, $H_{0}$, where the mean-field ground
state is highly degenerate. For example, the tiniest chiral interactions
arising from 3-spin interactions mediated by the chiral electronic channels,
induce a state with a huge spin-chirality in one of the valleys. The
degeneracy of the mean-field ground state of $H_{0}$ is, however, lifted by
quantum and thermal fluctuations. While in most systems, quantum and thermal
fluctuations stabilize the same type of order by such an order-by-disorder
mechanism, this is not the case in our model where quantum fluctuations prefer
coplanar spin-order, while classical fluctuations favor ferromagnetic spin
order on top of a 120∘ valley order. The ferromagnetic spin order is, however,
highly unconventional. Due to the coupling of spin- and valley degrees of
freedom, quantum or thermal fluctuations around the spin-ferromagnetic state
are _always_ chiral with a finite spin-chirality. In the classical limit, this
fluctuation effect is enhanced, Eq (12), due to the coupling to a pseudo
Goldstone mode characteristic for the classical order-by-disorder mechanism.

Our results on quantum fluctuations are based on a spin-wave calculation,
which formally becomes exact in a large $M$ limit. As $M=1$, this result
remains speculative. An alternative scenario is that for $M=1$ quantum
fluctuations around the highly degenerate mean-field state induce a spin-
valley entangled quantum liquid. It would be interesting to test these very
different scenarios in future numerical studies. Furthermore, we expect that
the system will host many more commensurate and incommensurate phases and,
potentially, quantum liquids, when other localized modes and different
electronic fillings are considered.

###### Acknowledgements.

We thank Martin Zirnbauer, Guo-Yi Zhu, Ciarán Hickey, Shahal Ilani, Nick
Bultinck, Johannes Hofmann, Peter Cha and Hongki Min for useful discussions.
Financial support of the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) within CRC1238 (project number 277146847, C02 and C04) and CRC183
(project number 277101999, A01 and A04) is acknowledged. Jeyong Park also
acknowledges BCGS (Bonn-Cologne Graduate School) and ML4Q (Matter and Light
for Quantum Computing) for support. The numerical simulations were performed
on the Noctua 2 cluster at the Paderborn Center for Parallel Computing (PC2)
and the CHEOPS cluster at RRZK Cologne.

Note added: Recently, a preprint by Wittig et al. was published on ArXiv [62]
that also considers localized states coupled to a network of chiral modes in
the twisted bilayer graphene subject to an interlayer bias.

## References

  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva,  _Van der Waals heterostructures_ , Nature 499, 419 (2013).
  * Andrei  _et al._ [2021] E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young,  _The marvels of moiré materials_ , Nature Reviews Materials 6, 201 (2021).
  * Chen  _et al._ [2019] 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, and F. Wang,  _Signatures of tunable superconductivity in a trilayer graphene moiré superlattice_ , Nature 572, 215 (2019).
  * Xie  _et al._ [2022] Y.-M. Xie, C.-P. Zhang, J.-X. Hu, K. F. Mak, and K. T. Law,  _Valley-Polarized Quantum Anomalous Hall State in Moiré ${\mathrm{MoTe}}_{2}/{\mathrm{WSe}}_{2}$ Heterobilayers_, Phys. Rev. Lett. 128, 026402 (2022).
  * Jin  _et al._ [2021] C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, L. Fu, J. Shan, and K. F. Mak,  _Stripe phases in WSe2/WS2 moiré superlattices_ , Nature Materials 20, 940 (2021).
  * Xu  _et al._ [2020] Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan,  _Correlated insulating states at fractional fillings of moiré superlattices_ , Nature 587, 214 (2020).
  * Li  _et al._ [2021] T. Li, S. Jiang, B. Shen, Y. Zhang, L. Li, Z. Tao, T. Devakul, K. Watanabe, T. Taniguchi, L. Fu, J. Shan, and K. F. Mak,  _Quantum anomalous Hall effect from intertwined moiré bands_ , Nature 600, 641 (2021).
  * Cao  _et al._ [2018a] 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, and P. Jarillo-Herrero,  _Correlated insulator behaviour at half-filling in magic-angle graphene superlattices_ , Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero,  _Unconventional superconductivity in magic-angle graphene superlattices_ , Nature 556, 43 (2018b).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald,  _Moiré bands in twisted double-layer graphene_ , Proceedings of the National Academy of Sciences of the United States of America 108, 12233 (2011).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath,  _Origin of Magic Angles in Twisted Bilayer Graphene_ , Phys. Rev. Lett. 122, 106405 (2019).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil,  _Origin of Mott Insulating Behavior and Superconductivity in Twisted Bilayer Graphene_ , Phys. Rev. X 8, 031089 (2018).
  * Kennes  _et al._ [2018] D. M. Kennes, J. Lischner, and C. Karrasch,  _Strong correlations and $d+\mathit{id}$ superconductivity in twisted bilayer graphene_, Phys. Rev. B 98, 241407 (2018).
  * Wu  _et al._ [2018] F. Wu, A. H. MacDonald, and I. Martin,  _Theory of Phonon-Mediated Superconductivity in Twisted Bilayer Graphene_ , Phys. Rev. Lett. 121, 257001 (2018).
  * Lian  _et al._ [2019] B. Lian, Z. Wang, and B. A. Bernevig,  _Twisted Bilayer Graphene: A Phonon-Driven Superconductor_ , Phys. Rev. Lett. 122, 257002 (2019).
  * Roy and Juričić [2019] B. Roy and V. Juričić,  _Unconventional superconductivity in nearly flat bands in twisted bilayer graphene_ , Phys. Rev. B 99, 121407 (2019).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean,  _Tuning superconductivity in twisted bilayer graphene_ , Science 363, 1059 (2019).
  * Isobe  _et al._ [2018] H. Isobe, N. F. Yuan, and L. Fu,  _Unconventional Superconductivity and Density Waves in Twisted Bilayer Graphene_ , Phys. Rev. X 8, 041041 (2018).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge,  _Electronic correlations in twisted bilayer graphene near the magic angle_ , Nature Physics 15, 1174 (2019).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy,  _Maximized electron interactions at the magic angle in twisted bilayer graphene_ , Nature 572, 95 (2019).
  * Xie  _et al._ [2021] Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby,  _Fractional Chern insulators in magic-angle twisted bilayer graphene_ , Nature 600, 439 (2021).
  * Ledwith  _et al._ [2020] P. J. Ledwith, G. Tarnopolsky, E. Khalaf, and A. Vishwanath,  _Fractional Chern insulator states in twisted bilayer graphene: An analytical approach_ , Phys. Rev. Res. 2, 023237 (2020).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig,  _Magic-Angle Twisted Bilayer Graphene as a Topological Heavy Fermion Problem_ , Phys. Rev. Lett. 129, 047601 (2022).
  * Thomson  _et al._ [2018] A. Thomson, S. Chatterjee, S. Sachdev, and M. S. Scheurer,  _Triangular antiferromagnetism on the honeycomb lattice of twisted bilayer graphene_ , Phys. Rev. B 98, 075109 (2018).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel,  _Ground State and Hidden Symmetry of Magic-Angle Graphene at even Integer Filling_ , Phys. Rev. X 10, 031034 (2020).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck,  _Kekulé Spiral Order at All Nonzero Integer Fillings in Twisted Bilayer Graphene_ , Phys. Rev. X 11, 041063 (2021).
  * Hofmann  _et al._ [2022] J. S. Hofmann, E. Khalaf, A. Vishwanath, E. Berg, and J. Y. Lee,  _Fermionic Monte Carlo Study of a Realistic Model of Twisted Bilayer Graphene_ , Phys. Rev. X 12, 011061 (2022).
  * Chou and Sarma [2022] Y.-Z. Chou and S. D. Sarma,  _Kondo lattice model in magic-angle twisted bilayer graphene_ ,  (2022), arXiv:2211.15682 .
  * Mora  _et al._ [2019] C. Mora, N. Regnault, and B. A. Bernevig,  _Flatbands and Perfect Metal in Trilayer Moiré Graphene_ , Phys. Rev. Lett. 123, 026402 (2019).
  * Khalaf  _et al._ [2019] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath,  _Magic angle hierarchy in twisted graphene multilayers_ , Phys. Rev. B 100, 085109 (2019).
  * [31] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero,  _Tunable Phase Boundaries and Ultra-Strong Coupling Superconductivity in Mirror Symmetric Magic-Angle Trilayer Graphene_ ,  arXiv:2012.01434 .
  * Chen  _et al._ [2020] S. Chen, M. He, Y. H. Zhang, V. Hsieh, Z. Fei, K. Watanabe, T. Taniguchi, D. H. Cobden, X. Xu, C. R. Dean, and M. Yankowitz,  _Electrically tunable correlated and topological states in twisted monolayer–bilayer graphene_ , Nature Physics 17, 374 (2020).
  * He  _et al._ [2020] M. He, Y. Li, J. Cai, Y. Liu, K. Watanabe, T. Taniguchi, X. Xu, and M. Yankowitz,  _Symmetry breaking in twisted double bilayer graphene_ , Nature Physics 17, 26 (2020).
  * Ramires and Lado [2021] A. Ramires and J. L. Lado,  _Emulating Heavy Fermions in Twisted Trilayer Graphene_ , Phys. Rev. Lett. 127, 26401 (2021).
  * Park  _et al._ [2022] J. M. Park, Y. Cao, L.-Q. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero,  _Robust superconductivity in magic-angle multilayer graphene family_ , Nature Materials 21, 877 (2022).
  * San-Jose and Prada [2013] P. San-Jose and E. Prada,  _Helical networks in twisted bilayer graphene under interlayer bias_ , Phys. Rev. B 88, 121408 (2013).
  * Efimkin and MacDonald [2018] D. K. Efimkin and A. H. MacDonald,  _Helical network model for twisted bilayer graphene_ , Phys. Rev. B 98, 35404 (2018).
  * Ramires and Lado [2018] A. Ramires and J. L. Lado,  _Electrically Tunable Gauge Fields in Tiny-Angle Twisted Bilayer Graphene_ , Phys. Rev. Lett. 121, 146801 (2018).
  * Rickhaus  _et al._ [2018] P. Rickhaus, J. Wallbank, S. Slizovskiy, R. Pisoni, H. Overweg, Y. Lee, M. Eich, M.-H. Liu, K. Watanabe, T. Taniguchi, T. Ihn, and K. Ensslin,  _Transport Through a Network of Topological Channels in Twisted Bilayer Graphene_ , Nano Letters 18, 6725 (2018).
  * Huang  _et al._ [2018] S. Huang, K. Kim, D. K. Efimkin, T. Lovorn, T. Taniguchi, K. Watanabe, A. H. MacDonald, E. Tutuc, and B. J. LeRoy,  _Topologically Protected Helical States in Minimally Twisted Bilayer Graphene_ , Phys. Rev. Lett. 121, 37702 (2018).
  * Yoo  _et al._ [2019] H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, P. Cazeaux, S. H. Sung, R. Hovden, A. W. Tsen, T. Taniguchi, K. Watanabe, G.-C. Yi, M. Kim, M. Luskin, E. B. Tadmor, E. Kaxiras, and P. Kim,  _Atomic and electronic reconstruction at the van der Waals interface in twisted bilayer graphene_ , Nature Materials 18, 448 (2019).
  * Xu  _et al._ [2019] S. G. Xu, A. I. Berdyugin, P. Kumaravadivel, F. Guinea, R. Krishna Kumar, D. A. Bandurin, S. V. Morozov, W. Kuang, B. Tsim, S. Liu, J. H. Edgar, I. V. Grigorieva, V. I. Fal’ko, M. Kim, and A. K. Geim,  _Giant oscillations in a triangular network of one-dimensional states in marginally twisted graphene_ , Nature Communications 10, 4008 (2019).
  * Chou  _et al._ [2020] Y.-Z. Chou, F. Wu, and S. Das Sarma,  _Hofstadter butterfly and Floquet topological insulators in minimally twisted bilayer graphene_ , Phys. Rev. Res. 2, 033271 (2020).
  * De Beule  _et al._ [2020] C. De Beule, F. Dominguez, and P. Recher,  _Aharonov-bohm oscillations in minimally twisted bilayer graphene_ , Phys. Rev. Lett. 125, 096402 (2020).
  * De Beule  _et al._ [2021] C. De Beule, F. Dominguez, and P. Recher,  _Network model and four-terminal transport in minimally twisted bilayer graphene_ , Phys. Rev. B 104, 195410 (2021).
  * Chou  _et al._ [2021] Y.-Z. Chou, F. Wu, and J. D. Sau,  _Charge density wave and finite-temperature transport in minimally twisted bilayer graphene_ , Phys. Rev. B 104, 045146 (2021).
  * Attig  _et al._ [2021] J. Attig, J. Park, M. M. Scherer, S. Trebst, A. Altland, and A. Rosch,  _Universal principles of moiré band structures_ , 2D Materials 8, 044007 (2021).
  * [48] See the Supplemental Material at DOI for details on (i) a symmetry analysis to obtain the form of the matrices $\hat{w}_{ii^{\prime}}$ and $\hat{\lambda}_{ij}$ in Eq. (4), (ii) a simple effective Hamiltonian to fully reproduce the band structure in Fig. 2c, (iii) the mean-field phase diagram, (iv) the $SU(4)$ spin wave theory, and (v) the semi-classical Monte Carlo simulations. It includes Refs. [63, 64, 65, 66, 67, 68].
  * Kugel’ and Khomskiĭ [1982] K. I. Kugel’ and D. I. Khomskiĭ,  _The Jahn-Teller effect and magnetism: transition metal compounds_ , Soviet Physics Uspekhi 25, 231 (1982).
  * Wietek and Läuchli [2017] A. Wietek and A. M. Läuchli,  _Chiral spin liquid and quantum criticality in extended S= 12 Heisenberg models on the triangular lattice_ , Phys. Rev. B 95, 035141 (2017).
  * Gong  _et al._ [2017] S. S. Gong, W. Zhu, J. X. Zhu, D. N. Sheng, and K. Yang,  _Global phase diagram and quantum spin liquids in a spin- 12 triangular antiferromagnet_ , Phys. Rev. B 96, 075116 (2017).
  * Szasz  _et al._ [2020] A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore,  _Chiral Spin Liquid Phase of the Triangular Lattice Hubbard Model: A Density Matrix Renormalization Group Study_ , Phys. Rev. X 10, 021042 (2020).
  * Chen  _et al._ [2022] B. B. Chen, Z. Chen, S. S. Gong, D. N. Sheng, W. Li, and A. Weichselbaum,  _Quantum spin liquid with emergent chiral order in the triangular-lattice Hubbard model_ , Phys. Rev. B 106, 094420 (2022).
  * Sur  _et al._ [2022] S. Sur, A. Udupa, and D. Sen,  _Driven Hubbard model on a triangular lattice: Tunable Heisenberg antiferromagnet with a chiral three-spin term_ , Phys. Rev. B 105, 054423 (2022).
  * [55] C. Kuhlenkamp, W. Kadow, A. Imamoglu, and M. Knap,  _Tunable topological order of pseudo spins in semiconductor heterostructures_ , arXiv:2209.05506 .
  * Stoudenmire  _et al._ [2009] E. M. Stoudenmire, S. Trebst, and L. Balents,  _Quadrupolar correlations and spin freezing in $S=1$ triangular lattice antiferromagnets_, Phys. Rev. B 79, 214436 (2009).
  * Hickey and Paramekanti [2014] C. Hickey and A. Paramekanti,  _Thermal Phase Transitions of Strongly Correlated Bosons with Spin-Orbit Coupling_ , Phys. Rev. Lett. 113, 265302 (2014).
  * Landau and Binder [2014] D. P. Landau and K. Binder,  _A Guide to Monte Carlo Simulations in Statistical Physics_ (Cambridge University Press, 2014).
  * Mermin and Wagner [1966] N. D. Mermin and H. Wagner,  _Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models_ , Phys. Rev. Lett. 17, 1133 (1966).
  * Villain  _et al._ [1980] J. Villain, R. Bidaux, J.-P. Carton, and R. Conte,  _Order as an effect of disorder_ , J. Phys. 41, 1263 (1980).
  * Chalker  _et al._ [1992] J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender,  _Hidden order in a frustrated system: Properties of the Heisenberg Kagomé antiferromagnet_ , Phys. Rev. Lett. 68, 855 (1992).
  * Wittig  _et al._ [2023] P. Wittig, F. Dominguez, C. D. Beule, and P. Recher,  _Localized states coupled to a network of chiral modes in minimally twisted bilayer graphene_ (2023), arXiv:2303.03901 [cond-mat.mes-hall] .
  * Wietek  _et al._ [2022] A. Wietek, J. Wang, J. Zang, J. Cano, A. Georges, and A. Millis,  _Tunable stripe order and weak superconductivity in the Moiré Hubbard model_ , Phys. Rev. Res. 4, 043048 (2022).
  * Colpa [1978] J. Colpa,  _Diagonalization of the quadratic boson hamiltonian_ , Physica A: Statistical Mechanics and its Applications 93, 327 (1978).
  * [65] L. Gresista, D. Kiese, S. Trebst, and M. M. Scherer,  _Spin-valley magnetism on the triangular moiré lattice with SU(4) breaking interactions_ , arXiv:2210.080256 .
  * Muller [1959] M. E. Muller,  _A Note on a Method for Generating Points Uniformly on N-Dimensional Spheres_ , Commun. ACM 2, 19–20 (1959).
  * Alzate-Cardona  _et al._ [2019] J. D. Alzate-Cardona, D. Sabogal-Suárez, R. F. L. Evans, and E. Restrepo-Parra,  _Optimal phase space sampling for Monte Carlo simulations of Heisenberg spin systems_ , Journal of Physics: Condensed Matter 31, 095802 (2019).
  * Bauer and Freyer [2020] C. Bauer and F. Freyer,  _BinningAnalysis.jl: Standard error estimation tools_ (2020).

##  Appendix A Symmetry analysis for the network model

In this Appendix, we discuss the symmetry of the moiré system considered in
the main text, and identify the form of the coupling $\hat{w}$ and
$\hat{\lambda}$ (Eq. (4)) using a symmetry analysis.

###  A.1 Symmetry

The underlying lattice or discrete symmetries of our moiré model are the moiré
translation, the 120∘ rotation ($C_{3}$), the mirror ($M_{y}$) with respect to
the $y$ axis, the inversion ($\mathcal{I}$), and the time reversal symmetry
($\mathcal{T}$) and combinations thereof. Since the moiré potential varies
smoothly so that the large-momentum transfer is highly suppressed, it is a
good approximation to consider the two valley sectors separately. Therefore,
we only consider the valley-conserving symmetries generated by $C_{3}$,
$M_{y}$, and $\mathcal{T}\mathcal{I}$, in the following.

The underlying lattice symmetry group is the dehedral group $D_{3}$ with
$C_{3}$ and $M_{y}$. More precisely, in order to deal with the spinor wave
function properly, one has to take into account a minus sign under $2\pi$
rotation, and therefore consider the dicyclic group $\text{Dic}_{3}$, which
extends $D_{3}$. The character table of $\text{Dic}_{3}$ is as follows:

$\mathbf{Dic}_{3}$ |  $\mathbf{A}$ |  $\mathbf{B}$ |  $\mathbf{C}$ |  $\mathbf{D}$ |  $\mathbf{E}$ |  $\mathbf{F}$
---|---|---|---|---|---|---
$\rho_{1}$ |  $1$ |  $1$ |  $1$ |  $1$ |  $1$ |  $1$
$\rho_{2}$ |  $1$ |  $1$ |  $1$ |  $-1$ |  $-1$ |  $1$
$\rho_{3}$ |  $1$ |  $-1$ |  $-1$ |  $i$ |  $-i$ |  $1$
$\rho_{4}$ |  $1$ |  $-1$ |  $-1$ |  $-i$ |  $i$ |  $1$
$\rho_{5}$ |  $2$ |  $2$ |  $-1$ |  $0$ |  $0$ |  $-1$
$\rho_{6}$ |  $2$ |  $-2$ |  $1$ |  $0$ |  $0$ |  $-1$

with the 6 irreducible representations $\rho_{j=1,\cdots,6}$ and the
equivalent classes
$\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D},\mathbf{E},\mathbf{F}$ which are
given by

$\mathbf{A}\equiv\\{1\\},\ \mathbf{B}\equiv\\{C^{3}_{3}\\},\
\mathbf{C}\equiv\\{C_{3},C^{5}_{3}\\},\
\mathbf{D}\equiv\\{M_{y}C_{3},M_{y}C^{3}_{3},M_{y}C^{5}_{3}\\},\
\mathbf{E}\equiv\\{M_{y},M_{y}C^{2}_{3},M_{y}C^{4}_{3}\\},\
\mathbf{F}\equiv\\{C^{2}_{3},C^{4}_{3}\\}.$

For states with one electron per mode, only three of the representations,
$\rho_{3}$, $\rho_{4}$, $\rho_{6}$, are relevant since the condition
$C^{3}_{3}=-1$ should be fulfilled for the spinor wave functions. While
$\rho_{3}$ and $\rho_{4}$ are one-dimensional representations, $\rho_{6}$ is a
two-dimensional representation.

As shown in Fig. 1, the moiré system is effectively described by the network
model where localized states form at the junction of three 1D channels (for
each valley) and are weakly coupled to the channels. Below, employing a
symmetry analysis, we shall find the form of the inter-channel coupling
$\hat{w}$ and the coupling of 1D channels to localized states $\hat{\lambda}$
(Eq. (4)).

###  A.2 Localized modes

Localized states can be labeled by representation $\rho_{j=3,4,6}$ of the
group $\text{Dic}_{3}$. From the character table, one can find relevant
matrices for the transformation $C_{3}$ and $M_{y}$ in each of the
representations:
$C_{3}^{\rho_{3}}=-1,C_{3}^{\rho_{4}}=-1,C_{3}^{\rho_{6}}=e^{i\pi\sigma_{z}/3}$
and $M_{y}^{\rho_{3}}=-i,\
M_{y}^{\rho_{4}}=i,M_{y}^{\rho_{6}}=e^{i\pi\sigma_{y}/2}$, respectively. A
simple way to identify the symmetry of a given localized state from the band-
structure calculation is to analyze the symmetry properties of the
eigenfunction of flat bands at the $\Gamma$ point. For each of the three
representations we find examples in our band structure calculations, see Fig.
2c of the main text.

###  A.3 One-dimensional channels

From scaling, one finds that the width of the 1D channels is given by
$d_{1D}\sim(vL/u_{s})^{1/2}$ with $d_{1D}\ll L$ for $u_{s}\gg vG$. $d_{1D}\ll
L$ renders that the coupling of 1D channels being far apart is highly
suppressed, and thus it is a good approximation to only consider the coupling
of the neighboring channels at the junction. The symmetry properties are
determined by the spinor structure of the eigenfunctions. The three channels
at the junction are related to each other by a $120^{\circ}$ rotation matrix,
$C_{3}^{\text{1D}}$, and the mirror transformation matrix,
$M_{y}^{\text{1D}}$, e.g., for the $K$ valley, written by

$\displaystyle C_{3}^{\text{1D}}=\left(\begin{array}[]{ccc}0&0&-1\\\ 1&0&0\\\
0&1&0\end{array}\right),\
M_{y}^{\text{1D}}=\left(\begin{array}[]{ccc}-i&0&0\\\ 0&0&i\\\
0&i&0\end{array}\right).$ (19)

For the $K^{\prime}$ valley, the mirror transformation matrix has an extra
overall minus sign as $K$ and $K^{\prime}$ are related by time-reversal and
thus by complex conjugation. The relative minus sign in one of the matrix
elements of $C_{3}$ and $M_{y}$ reflects the fact that the $2\pi$ rotation
preinor wave function gets a minus sign. Also
$(\mathcal{I}\mathcal{T})^{\text{1D}}=C\mathds{1}_{3}$, where $C$ is the
complex conjugation and $\mathds{1}_{3}$ is the 3$\times$3 unit matrix. The
inter-channel tunneling $H_{w}$ is written as

$\displaystyle H_{w}=$ $\displaystyle\sum_{\text{crossing at
}{\bm{R}}_{m}}\hat{w}_{ii^{\prime}}\Psi^{\dagger}_{n,i,\alpha,\sigma}(\rho_{n,i,m})\Psi_{n^{\prime},i^{\prime},\alpha,\sigma}(\rho_{n^{\prime},i^{\prime},m}).$
(20)

By imposing the symmetry constraints,
$\hat{w}=(C_{3}^{\text{1D}})^{\dagger}\hat{w}C_{3}^{\text{1D}}$,
$\hat{w}=(M_{y}^{\text{1D}})^{\dagger}\hat{w}M_{y}^{\text{1D}}$, and
$\hat{w}=((\mathcal{I}\mathcal{T})^{\text{1D}})^{\dagger}\hat{w}(\mathcal{I}\mathcal{T})^{\text{1D}}$,
one can obtain $\hat{w}$, parameterized by a single real parameter $w$ as

$\displaystyle\hat{w}=$ $\displaystyle wL\left(\begin{array}[]{ccc}0&1&-1\\\
1&0&1\\\ -1&1&0\end{array}\right).$ (24)

The size of moiré unit cell, $L$, is used such that $w$ has units of energy.

###  A.4 Coupling between localized and propagating modes

The coupling of 1D channels to localized modes, Eq. (4), depends on to which
representation of $\text{Dic}_{3}$ the localized mode belongs. Similarly to
the $\hat{w}$ matrix, this coupling matrix $\hat{\lambda}$ can be obtained by
the symmetry constraints,
$\hat{\lambda}=(C_{3}^{\rho_{j}})^{\dagger}\hat{\lambda}C_{3}^{\text{1D}}$,
$\hat{\lambda}=(M_{y}^{\rho_{j}})^{\dagger}\hat{\lambda}M_{y}^{\text{1D}}$,
and
$\hat{\lambda}=((\mathcal{I}\mathcal{T})^{\rho_{j}})^{\dagger}\hat{\lambda}(\mathcal{I}\mathcal{T})^{\text{1D}}$,
for each of the representation $j=3,4,6$. For the $K$ valley, it is given by

$\displaystyle H_{\lambda}=$ $\displaystyle\sum_{\text{crossing at
}{\bm{R}}_{m}}\hat{\lambda}_{ij}\,d^{\dagger}_{m,\alpha,\sigma,j}\Psi_{n,i,\alpha,\sigma}(\rho_{n,i,m})+h.c.,$
$\displaystyle\hat{\lambda}=$
$\displaystyle\lambda\frac{\sqrt{L}}{\sqrt{3}}\left\\{\begin{array}[]{ll}(1,-1,1)&\rho_{3}\
\text{representation}\\\\[2.84526pt] (0,0,0)&\rho_{4}\
\text{representation}\\\\[2.84526pt]
\left(\begin{array}[]{ccc}1&e^{i\pi/3}&e^{2i\pi/3}\\\
-i&-ie^{-i\pi/3}&-ie^{-2i\pi/3}\end{array}\right)&\rho_{6}\
\text{representation}\end{array}\right.$ (30)

Interestingly, the localized modes in the $\rho_{4}$ representation do not
couple to the closest 1D channels at all. This vanishing coupling can be
understood by the symmetry of the system. For example, let us consider the
coupling of a localized mode in the $\rho_{4}$ representation to the
neighboring channel propagating in the $y$ direction. While the localized mode
has the eigenvalue $i$ under the mirror transformation $M_{y}$, the channel
has the eigenvalue $-i$. It implies that the coupling matrix $\hat{\lambda}$
has to have an extra minus sign under $M_{y}$, and therefore has to vanish.
The same argument can be applied to the two other neighboring channels with
the symmetry transformation $M_{y}C_{3}$, and $M_{y}C_{3}^{2}$. Nevertheless,
tunneling to channels further away are still possible but exponentially
suppressed in the ratio of potential and $vG$ as $e^{-L/d_{1D}}\sim
e^{-(u_{s}L/v)^{1/2}}$.

The coupling matrices for the $K^{\prime}$ valley can be obtained by time
reversal, which implies that the time-reversed partner of $\rho_{3}$ and
$\rho_{4}$ have an identical coupling matrix $\hat{\lambda}$. Formally, under
time-reversal $\rho_{3}$ maps to $\rho_{4}$ but this effect is compensated
because also in Eq. (A.4) the entries for $\rho_{3}$ and $\rho_{4}$ are
exchanged when one switches from $K$ to $K^{\prime}$. The matrix
$\hat{\lambda}$ for $\rho_{6}$ in the $K^{\prime}$ valley is obtained by
complex conjugation of $\hat{\lambda}$.

##  Appendix B Effective Hamiltonian in momentum space

In this section, we construct an effective Hamiltonian for the network model
in momentum space and show that the band structure obtained from the
diagonalization of the moiré potential model, Eq. (1), is nicely fitted by the
band structure obtained form this effective model.

As shown in Fig. 2(c) for the band structure in the main text (also the blue
solid curves in Fig. 7), the inter-channel coupling and the coupling of 1D
channels to localized states are very weak. For example, a blowup of the band
structure near the $\Gamma$ point (indicated by the black curcle) in Fig. 7
shows that a level-repulsion of the three 1D channels is an order of hundreds
of $\mu$eV. Such weak couplings (i.e., $w,\lambda\ll vG$) allow us to
construct an effective Hamiltonian from Eqs. (3) and (4). Depending on the
number of localized modes, $N_{\rm{loc}}$, within the energy window of
$\epsilon\in[-vG/2,vG/2]$ near the chemical potential, this effective
Hamiltonian $H_{\rm{eff}}(k_{x},k_{y})$ for a given valley can be described by
a $(3+N_{\rm{loc}})\times(3+N_{\rm{loc}})$ matrix for each momentum. For
example, the band structure for the $K$ valley shown in Fig. 2c can be
reproduced with high precision by the $7\times 7$ matrix, given by

$\displaystyle
H_{\rm{eff}}(k_{x},k_{y})=\left(\begin{array}[]{ccccccc}-vk_{y}+\epsilon_{\rm{1D}}&w&-w&\frac{\lambda_{6}}{\sqrt{3}}&i\frac{\lambda_{6}}{\sqrt{3}}&\frac{\lambda_{3}}{\sqrt{3}}&0\\\
w&-\frac{v}{2}\left(\sqrt{3}k_{x}-k_{y}\right)+\epsilon_{\rm{1D}}&w&\frac{\lambda_{6}}{\sqrt{3}}e^{-i\frac{\pi}{3}}&i\frac{\lambda_{6}}{\sqrt{3}}e^{i\frac{\pi}{3}}&-\frac{\lambda_{3}}{\sqrt{3}}&0\\\
-w&w&\frac{v}{2}\left(\sqrt{3}k_{x}+k_{y}\right)+\epsilon_{\rm{1D}}&\frac{\lambda_{6}}{\sqrt{3}}e^{-2i\frac{\pi}{3}}&i\frac{\lambda_{6}}{\sqrt{3}}e^{2i\frac{\pi}{3}}&\frac{\lambda_{3}}{\sqrt{3}}&0\\\
\frac{\lambda_{6}}{\sqrt{3}}&\frac{\lambda_{6}}{\sqrt{3}}e^{i\frac{\pi}{3}}&\frac{\lambda_{6}}{\sqrt{3}}e^{2i\frac{\pi}{3}}&\epsilon_{6}&0&0&0\\\
-i\frac{\lambda_{6}}{\sqrt{3}}&-i\frac{\lambda_{6}}{\sqrt{3}}e^{i\frac{\pi}{3}}&-i\frac{\lambda_{6}}{\sqrt{3}}e^{2i\frac{\pi}{3}}&0&\epsilon_{6}&0&0\\\
\frac{\lambda_{3}}{\sqrt{3}}&-\frac{\lambda_{3}}{\sqrt{3}}&\frac{\lambda_{3}}{\sqrt{3}}&0&0&\epsilon_{3}&0\\\
0&0&0&0&0&0&\epsilon_{4}\end{array}\right).$ (38)

The first three rows and columns of this effective Hamiltonian describe the 1D
channels propagating along the three directions, the 4- and 5th rows and
columns describe degenerate localized states that belong to the $\rho_{6}$
representation with energy $\epsilon_{6}$, and the 6 and 7th row and column
correspond to localized states that belongs to the $\rho_{3}$ and $\rho_{4}$
representation with energy $\epsilon_{3}$ and $\epsilon_{4}$, respectively. As
shown in Eqs. (20) and (A.4), the form of the coupling terms is determined by
the symmetry of the system and the representation of the localized states.

Figure 7: Comparison of two band structures. Left: The band structure (blue
solid curves) obtained from diagonalizing the original moiré potential model,
Eq. (1), and the band structure (red dashed curves) using the effective
Hamiltonian, Eq. (38). The parameters of Eq. (38) can be obtained by fitting
the blue solid curves; $w=0.2$meV, $\lambda_{6}=4.3$meV,
$\lambda_{3}=11.3$meV, $\epsilon_{\rm{1D}}=63.4$meV, $\epsilon_{6}=114$meV,
$\epsilon_{4}=154$meV, and $\epsilon_{3}=162$meV. Those band structures match
remarkably, confirming the validity of our network model. For the blue curves,
we have used the parameters, $u_{AA}=u_{BB}=0$ and $u_{AB}=u_{BA}=\frac{8\pi
v}{\sqrt{3}L}$. Right: a blowup of the band structures near the $\Gamma$ point
(indicated by the black circle) where the three propagating bands meet. Weak
level repulsion between the three bands is clearly shown, not visible in the
left band structure.

Figure 7 compares band structures obtained from two different Hamiltonians.
The band structure plotted in blue dashed curves is numerically obtained from
the diagonalization of the moiré model, Eq. (1). On the other hand, the band
structure plotted in blue solid line is obtained from the $7\times 7$
effective Hamiltonian in Eq. (38) with the parameters of $w=0.2$meV,
$\lambda_{6}=4.3$meV, $\lambda_{3}=11.3$meV, $\epsilon_{\rm{1D}}=63.4$meV,
$\epsilon_{6}=114$meV, $\epsilon_{4}=154$meV, and $\epsilon_{3}=162$meV, which
can be found by fitting to the blue dashed curves. Those band structures
remarkably match well, including silent features. It confirms the validity of
our network model.

##  Appendix C Mean-field phase diagram

In this section, we discuss the mean-field phase diagram of the Hamiltonian,
Eq. (III), given by

$\displaystyle H_{\rm{sv}}$ $\displaystyle=J_{2}\sum_{\langle m_{1}\to
m_{2}\rangle_{c}}(1+{\bm{\sigma}}_{m_{1}}\cdot{\bm{\sigma}}_{m_{2}})(e^{i\varphi}\tau^{+}_{m_{1}}\tau^{-}_{m_{2}}+h.c.)+J_{3}\\!\\!\\!\\!\sum_{p=\triangleright/\triangleleft,(m_{1},m_{2},m_{3})_{p}}\\!\\!\\!\\!(-1)^{p}{\bm{\sigma}}_{m_{1}}\cdot({\bm{\sigma}}_{m_{2}}{\bm{\times}}{\bm{\sigma}}_{m_{3}})\Big{(}\prod_{i=1}^{3}P_{m_{i}}^{+}-\prod_{i=1}^{3}P_{m_{i}}^{-}\Big{)}$
$\displaystyle+J_{2}^{\prime}\sum_{\langle
m_{1},m_{2}\rangle_{c}}(1+{\bm{\sigma}}_{m_{1}}\cdot{\bm{\sigma}}_{m_{2}})(1+\tau^{z}_{m_{1}}\tau^{z}_{m_{2}})\,\,\,\,.$
(39)

Compared with the Hamiltonians (Eqs. (II.4)-(II.4)) obtained from a
perturbation theory in $J$ and $w$, the Hamiltonian contains only the nearest
neighbor interactions. The continuous symmetries of $H_{\rm{sv}}$ are
$U(1)\times SU(2)\times SU(2)$ generated by $\tau^{z}$, $P^{+}{\bm{\sigma}}$
and $P^{-}{\bm{\sigma}}$. Assuming that the filling of the localized states is
unity, the states on the $SU(4)$ space are spanned by a 4-component complex
vector. In this basis, we solve the self-consistent mean-field equations
iteratively with the fixed unit cell. We use 4 different unit cells with one-,
two-, three-, four-sublattices. We find that either a one-sublattice or a
three-sublattice solution has the lowest energy. Therefore we consider only
these cases in the following.

The model is parametrized by 3 dimensionless parameters $J_{3}/J_{2}$,
$J_{2}^{\prime}/J_{2}$ and $\varphi$. Experimentally, we expect
$|J_{3}/J_{2}|,|J_{2}^{\prime}/J_{2}|\ll 1$ while $\varphi$ can take
arbitarily large values. The parameter $\varphi$ can be viewed as an Aharonov-
Bohm phase acquired by a particle with charge $\tau^{z}$ moving along a side
of the triangular loop. As $3\varphi$ is the total phase accumulated along the
triangular loop, it is possible to use a transformation to change the total
phase by $2\pi$ which is equivalent to a change of $\varphi$ by $n2\pi/3$,
$n\in\mathbb{Z}$. This is achieved by doing $\tau^{z}$ rotations of the spins
on the A, B, C (see the inset of Fig. 3) by $0$, $2\pi/3$, $4\pi/3$

$\displaystyle\varphi\rightarrow\varphi+\frac{2\pi}{3},\,\,\,{\bm{\tau}}_{\rm
A}\rightarrow{\bm{\tau}}_{\rm A},\,\,\,{\bm{\tau}}_{\rm B}\rightarrow
e^{i\tau^{z}_{\rm B}\frac{2\pi}{3}}{\bm{\tau}}_{\rm B}e^{-i\tau^{z}_{\rm
B}\frac{2\pi}{3}},\,\,\,{\bm{\tau}}_{\rm C}\rightarrow e^{i\tau^{z}_{\rm
C}\frac{4\pi}{3}}{\bm{\tau}}_{\rm C}e^{-i\tau^{z}_{\rm C}\frac{4\pi}{3}}.$
(40)

If one uses the Aharonov-Bohm analogy described above, this would be a gauge
transformation. In our system, however, a rotation by $\tau^{z}$ changes the
physical state. The transformation maps, for example, a ferromagnetic state
obtained for $\varphi=\pi$ to a state with $120^{\circ}$ order for
$\varphi=\pm\pi/3$, see below.

Figure 8:  Mean-field phase diagram with the order parameter
$(\chi^{+}-\chi^{-})$ in $(\varphi,J_{3})$ parameter space. $J_{2}^{\prime}=0$
and $J_{2}=1$. This diagram extends Fig. 3a, which only covers small $\varphi$
(defined modulo $2\pi/3$) and weak $J_{3}$ near the highly singular points.
While three different types of coplanar phase (two $120^{\circ}$ orders and a
ferromagnetic order) are realized in weak $J_{3}$, the chiral order achieves
the maximal values of chirality $\chi^{+}-\chi^{-}=-1$ for sufficiently strong
$J_{3}>0$ (purple region).

As $J_{2}>0$ is the largest term, we first analyze the case
$J_{3}=J^{\prime}_{2}=0$. For $\varphi=0$, the Hamiltonian, Eq. (C), becomes a
variant of the Kugel-Khomskii model [49]

$\displaystyle H_{0}$ $\displaystyle=H_{\rm{sv}}=2J_{2}\sum_{\langle m_{1}\to
m_{2}\rangle_{c}}(1+{\bm{\sigma}}_{m_{1}}\cdot{\bm{\sigma}}_{m_{2}})(\tau^{x}_{m_{1}}\tau^{x}_{m_{2}}+\tau^{y}_{m_{1}}\tau^{y}_{m_{2}})$
$\displaystyle=2J_{2}\sum_{\langle m_{1}\to
m_{2}\rangle_{c}}{\bm{\mu}}^{1}_{m_{1}}\cdot{\bm{\mu}}^{1}_{m_{2}}+{\bm{\mu}}^{2}_{m_{1}}\cdot{\bm{\mu}}^{2}_{m_{2}},$
(41)

with the 4-component vectors given by
${\bm{\mu}}^{1}_{m}=(\tau^{x}_{m},\tau^{y}_{m}\sigma^{x}_{m},\tau^{y}_{m}\sigma^{y}_{m},\tau^{y}_{m}\sigma^{z}_{m})$
and
${\bm{\mu}}^{2}_{m}=(\tau^{y}_{m},\tau^{x}_{m}\sigma^{x}_{m},\tau^{x}_{m}\sigma^{y}_{m},\tau^{x}_{m}\sigma^{z}_{m})$.
Numerically, we find from our mean-field analysis that the ground state has
the properties that the vectors $\langle{\bm{\mu}}^{n}_{m}\rangle$, $n=1,2$,
have the norm $1$, and show 120∘ order, such that
$\langle{\bm{\mu}}^{n}_{m_{1}}\rangle\cdot\langle{\bm{\mu}}^{n}_{m_{2}}\rangle=\cos(2\pi/3)=-\frac{1}{2}$
for neighboring sites. Thus, this specific type of 120∘ order is realized with
4-component vectors. More precisely, two different types of the 120∘ order are
realized in the ground-state manifold of $H_{0}$: a right-handed and left-
handed 120∘ order. The states of the right (left)-handed 120∘ order rotate in
anti-clockwise (clockwise) order around the $\tau^{z}$ orientation along the
triangular loop ($A\to B\to C$) as

$\displaystyle\psi_{B}=e^{\mp
i\tau^{z}\frac{2\pi}{3}}\psi_{A},\,\,\,\,\,\,\psi_{C}=e^{\mp
i\tau^{z}\frac{4\pi}{3}}\psi_{A}.$ (42)

Those 120∘ orders are staggered such that the states in the neighboring
triangles circulate in the opposite direction.

Surprisingly, the above described 120∘ orders have an extra degree of freedom.
This degree of freedom can be revealed by analyzing the valley-projected
magnetization vectors, $\langle{\bm{\sigma}}P^{\pm}\rangle$, which are length
$1/2$ in the ground-state manifold. In one of the two valley, the
magnetization is always ferromagnetic, but in the other valley a non-coplanar
spin configuration is possible. This non-coplanar spin configuration can be
fully characterized by the opening angle $0\leq\theta\leq\pi$ of the
magnetization vectors of the three neighboring sites. While the opening angle
$0$ and $\pi$ correspond to two distinct spin ferromagnetic states, see below,
$\pi/2$ corresponds to the coplanar $120^{\circ}$ states.

It is possible to write down analytically the spinor wavefunctions in the A, B
and C sublattices within the ground-state manifold. Up to rotations using the
continuous symmetries $U(1)\times SU(2)\times SU(2)$, they take the form

$\displaystyle\psi^{+}_{A}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\ 1\\\ 0\\\
0\end{array}\right),\,\,\,\,\,\psi^{+}_{B}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
-\frac{1}{2}+i\frac{\sqrt{3}}{2}\cos\theta\\\ 0\\\
i\frac{\sqrt{3}}{2}\sin\theta\end{array}\right),\,\,\,\,\,\psi^{+}_{C}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\
-\frac{1}{2}-i\frac{\sqrt{3}}{2}\cos\theta\\\ 0\\\
-i\frac{\sqrt{3}}{2}\sin\theta\end{array}\right)$ (55)

if the spins in the $+$ valley order ferromagnetically. If the spins in the
$-$ valley order ferromagnetically, one finds instead (again up to
transformations by the continuous symmetries $U(1)\times SU(2)\times SU(2)$)

$\displaystyle\psi^{-}_{A}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\\ 1\\\ 0\\\
0\end{array}\right),\,\,\,\,\,\psi^{-}_{B}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}-\frac{1}{2}-i\frac{\sqrt{3}}{2}\cos\theta\\\
1\\\ -i\frac{\sqrt{3}}{2}\sin\theta\\\
0\end{array}\right),\,\,\,\,\,\psi^{-}_{C}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}-\frac{1}{2}+i\frac{\sqrt{3}}{2}\cos\theta\\\
1\\\ i\frac{\sqrt{3}}{2}\sin\theta\\\ 0\end{array}\right).$ (68)

In both cases, $\theta$ is the opening angle characterizing the chiral spin
order. An opening angle $\varphi\neq 0,\pi,\pi/2$ leads to a finite valley-
projected chirality
$\chi^{\pm}=\langle{\bm{\sigma}}P^{\pm}\rangle_{A}\cdot(\langle{\bm{\sigma}}P^{\pm}\rangle_{B}{\bm{\times}}\langle{\bm{\sigma}}P^{\pm}\rangle_{C})$.
Using the $\Psi^{+}$ solutions above, we obtain

$\displaystyle\chi^{+}=0,\qquad\chi^{-}=\frac{3}{16}\sqrt{3}\cos\theta\sin^{2}\theta,$
(69)

while for the $\Psi^{-}$ solution we find

$\displaystyle\chi^{+}=-\frac{3}{16}\sqrt{3}\cos\theta\sin^{2}\theta,\qquad\chi^{-}=0.$
(70)

Ground states have an arbitrary $\theta$, and thus an arbitrary
$-\frac{1}{8}\leq\chi^{-}\leq\frac{1}{8}$ for the $\psi^{+}$ states or
$-\frac{1}{8}\leq\chi^{+}\leq\frac{1}{8}$ for the $\psi^{-}$ states. For
$\theta=0$ or $\theta=\pi$, the $\psi^{+}$ and $\psi^{-}$ solutions coincide
(up to trivial phases). The states with $\theta=0$ or $\pi$ in Eq. (55) are
spin-ferromagnet in both of the valleys, but have the distinct 120∘ order with
the left or right handedness, respectively.

Expanding $H_{0}$, at a energy minimum, in terms of three 4-component complex
states,
$|\psi_{m}^{0}\rangle\rightarrow|\psi_{m}^{0}\rangle+|\delta\psi_{m}\rangle$
with $m=A,B,C$, up to the second order, we obtain seven zero-modes among
$3\times(8-2)=18$ degrees of freedom. Those seven zero-modes show explicitly
that the ground-state manifold is 7-dimensional. Six of this seven modes arise
from the spontaneously broken continuous symmetries, the 7th mode, in
contrast, is related to a change of $\theta$, thus links states which are not
related by symmetry. An exception of this counting argument are the
ferromagnetic states at $\theta=0$ or $\theta=\pi$, which have a higher
symmetry. They have nine zero modes, from which five arise from spontaneously
broken symmetries while two each describe changes of the opening angle
$\theta$ either in the valley $+$ or the $-$ sector. Note that there are two
such modes per sector as magnetization vectors can tilt in two different
directions starting from the ferromagnetic configuration.

In the main text, we discuss the mean-field phase diagram arising from small
perturbations around the $\varphi=0$ point, see Fig. 3. In Fig. 8 we show the
analog of Fig. 3a but for an extended parameter range where $\varphi$ varies
from $-\pi$ to $\pi$ and we also allow for large values of $J_{3}$. We find
three types of non-chiral phases (red) which show either ferromagnetic or a
$120^{\circ}$ order in the vectors $\langle{\bm{\mu}}_{n}\rangle$ with
$n=1,2$.

Importantly, the phase diagram shows singular points not only at $\varphi=0$
but also at $\varphi=\pm\frac{2\pi}{3}$. Those singular points can be
understood by the enhanced symmetry of the $\varphi=0$ state discussed above
and the transformation of Eq. (40) which can be used to map the states at
$\varphi=\pm\frac{2\pi}{3}$ to $\varphi=0$. Due to the $\tau^{z}$ rotation,
Eq. (40), as adding $\varphi$ by $2\pi/3$ successively, the phase changes from
a left-handed $120^{\circ}$ ordered phase $\rightarrow$ a right-handed
$120^{\circ}$ phase $\rightarrow$ a ferromagnetic phase (more precisely, in
the vectors $\langle{\bm{\mu}}_{n}\rangle$ with $n=1,2$), and back to the
left-handed $120^{\circ}$ phase again. Such a transition between different
types of $120^{\circ}$ order was recently studied in the moiré Hubbard model
[63].

Since the $\varphi=0$ and $\varphi=\pm 2\pi/3$ points with
$J_{3}=J_{2}^{\prime}=0$ have a large, degenerate ground-state manifold, even
small perturbations which lift this degeneracy can lead to a giant effect
close to all three points as shown in Fig. 8. Note that the state with finite
chiral order in either the $+$ or $-$ sector, also breaks the discrete valley
symmetry leading to a finite $\langle\tau^{z}\rangle\neq 0$. For sufficiently
large $J_{3}$ (cf. Fig. 8), the $\tau^{z}$ symmetry is maximally broken with
$|\langle\tau^{z}\rangle|=1$ and also the chirality takes its maximal value,
$\chi^{+}-\chi^{-}=-1$ (purple region).

##  Appendix D Spin-wave theory

In this section, we perform a spin wave calculation in
$J_{2}^{\prime}=J_{3}=0$ and $\varphi=2\pi n/3$, $n\in\mathbb{Z}$, where
ground states are highly degenerate. The motivation of this spin wave
calculation is to investigate the effect of thermal and quantum fluctuation in
such a highly degenerate ground-state manifold. We show that by the thermal
order-by-disorder mechanism, the system selects spin ferromagnetic states in
both of the valley sectors from the ground-state manifold. This is contrasted
with that the quantum order-by-disorder mechanism favors 120∘ spin order in
one valley and ferromagnetic order in the other valley. The thermal order-by-
disorder mechanism leads to a mass gap linear in temperature for the soft
modes related with the opening angle $\theta$. Expanding around the classical
ground states (i.e., spin-ferromagnet in both of the valleys), we show that a
finite chirality $\chi^{\pm}$ is induced at finite temperatures by the
$\theta$-related soft modes. The non-analytic temperature dependence of the
chirality arises from the linear mass gap of the soft modes in temperatures.
These results remarkably match with results of the classical Monte-Carlo
simulation.

We first start by developing a spin wave theory for $SU(4)$ operators. Spin
wave theories become exact in certain large $M$ limits, where $M$ parametrize
representations of the group. In the $SU(2)$ case one uses the size of the
spin $s$ (with $M=2s$) and performs a $1/s$ expansion. In the $SU(4)$ case, we
choose a totally symmetric representation of the $SU(4)$ operators by (i)
writing the operators with bosonic creation and annihilation operators,
$a^{\dagger}_{i},a_{i}$ ($i=1,\dots,4$ and we suppress an extra site index
here), and (ii) fix the number of bosons (per site) to be $M$ using

$\displaystyle\hat{\Gamma}^{\ell}$
$\displaystyle=\sum_{i,i^{\prime}=1,\dots,4}\hat{a}_{i}^{\dagger}\gamma_{ii^{\prime}}^{\ell}\hat{a}_{i^{\prime}},\qquad
M=\sum_{i=1}^{4}\hat{a}_{i}^{\dagger}\hat{a}_{i}.$ (71)

Furthermore, we add an extra factor $1/M$ in front of the Hamiltonian, $H\to
H/M$, to make the large $M$ limit well defined, see below. The $M=1$ case
corresponds to the fundamental representation of $SU(4)$, realized if a single
electron is localized on each site. While we are interested in this limit, the
large $M$ theory to useful to derive a spin-wave theory in a controlled way.

Within the functional integral formalism, the partition function of our system
is expressed as

$\displaystyle\mathcal{Z}=\int D[\bar{a},a]D\lambda e^{-\int
d\tau\left(\sum_{i,s,{\bm{n}}}\bar{a}_{i,s,{\bm{n}}}\partial_{\tau}a_{i,s,{\bm{n}}}+\frac{1}{M}H(\bar{a},a)+i\sum_{s,{\bm{n}}}\lambda_{s,{\bm{n}}}(M-\sum_{i=1}^{4}\bar{a}_{i,s,{\bm{n}}}a_{i,s,{\bm{n}}})\right)}.$
(72)

where $a_{i,s,{\bm{n}}}$ are complex fields and the Lagrange multipliers
$\lambda_{s,{\bm{n}}}$ are used to implement the constraint on each site.
Rescaling the boson fields as
$a_{i,s,{\bm{n}}}=\sqrt{M}\tilde{a}_{i,s,{\bm{n}}}$ and using that $H$ is
quartic in these operators, we arrive

$\displaystyle\mathcal{Z}=\int D[\bar{\tilde{a}},\tilde{a}]D\lambda
e^{-MS_{\textrm{eff}}},\qquad
S_{\textrm{eff}}=\int_{0}^{\beta}d\tau\left(\sum_{i,s,{\bm{n}}}\bar{\tilde{a}}_{i,s,{\bm{n}}}\partial_{\tau}\tilde{a}_{i,s,{\bm{n}}}+H(\bar{\tilde{a}},\tilde{a})-i\sum_{s,{\bm{n}}}\lambda_{s,{\bm{n}}}\left(\sum_{i=1}^{4}\bar{\tilde{a}}_{i,s,{\bm{n}}}\tilde{a}_{i,s,{\bm{n}}}-1\right)\right)$
(73)

Due to the factor $M$ in front of $S_{\textrm{eff}}$, the functional integral
in the large $M$ limit is dominated by its saddle point and fluctuations
around the saddle point, which are controlled by $1/M$. Saddle point solutions
can be obtained by solving $\frac{\partial
S_{\textrm{eff}}}{\partial\tilde{a}_{i,s,{\bm{n}}}}|_{\tilde{a}_{\textrm{sp}},\lambda_{\textrm{sp}}}=\frac{\partial
S_{\textrm{eff}}}{\partial\bar{\tilde{a}}_{i,s,{\bm{n}}}}|_{\tilde{a}_{\textrm{sp}},\lambda_{\textrm{sp}}}=\frac{\partial
S_{\textrm{eff}}}{\partial\lambda_{s,{\bm{n}}}}|_{\tilde{a}_{\textrm{sp}},\lambda_{\textrm{sp}}}=0$.
Static saddle point solutions exactly correspond to the zero-temperature mean-
field solutions with $\lambda_{s,{\bm{n}}}=\lambda_{\textrm{sp}}$ being the
mean-field energy per site. Expanding $S_{\textrm{eff}}$ up to second order
around the saddle point solutions, the resulting action captures physics of
spin wave excitations. Although this saddle point approximation becomes more
accurate with large $M$, it also provides a good approximation even to the
$M=1$ case. Hereafter $M$ is set to 1, unless otherwise stated.

To understand the effect of quantum or thermal fluctuation on the degenerate
ground-state manifold of $H_{0}$ (Eq. (C)), we use the approach explained
above for $J_{2}^{\prime}=J_{3}=0$ and $\varphi=2\pi n/3$, $n\in\mathbb{Z}$
where the mean-field solution is highly degenerate. As a reference state
around which the action (Eq. (73)) is expanded, we take the spinor
wavefunctions $\psi^{0}_{s}(\alpha=\pm,\theta)=\psi_{s}^{\alpha}(\theta)$
(Eqs. (55) and (68)) that depend on the opening angle $\theta$ and $\alpha$.
$\alpha$ represents the valley sector in which the spin has the ferromagnetic
order. Then, low-energy states associated with the spin wave excitation can be
generally written, up to normalization, as

$\displaystyle\psi_{s}(\alpha,\theta)\sim\psi_{s}^{0}(\alpha,\theta)+\frac{1}{\sqrt{M}}\sum_{i=1}^{3}\hat{a}_{s,i}\delta\psi_{s,i}(\alpha,\theta).$
(74)

Here $\delta\psi_{s,i}$ are three 4-component unit vectors perpendicular to
the reference state $\psi_{s}^{0}(\pm,\theta)$. At this point, we find it
useful to switch back from the functional integral formalism to the operator
formalism, where it is more easy to keep track of commutation relation.
Expanding to the second order in $\hat{a}_{s,i}$ and performing the fourier
transform to momentum space, the Hamiltonian has a Bogoliubov-de Gennes (BdG)
form

$\displaystyle
H_{\textrm{BdG}}=\frac{1}{2}\sum_{{\bm{k}}}\begin{pmatrix}\hat{a}_{{\bm{k}}}^{\dagger}&\hat{a}_{-{\bm{k}}}\end{pmatrix}\begin{pmatrix}A({\bm{k}})&B({\bm{k}})\\\
B^{*}(-{\bm{k}})&A^{*}(-{\bm{k}})\end{pmatrix}\begin{pmatrix}\hat{a}_{{\bm{k}}}\\\
\hat{a}_{-{\bm{k}}}^{\dagger}\end{pmatrix}-\frac{1}{2}\sum_{{\bm{k}}}\textrm{Tr}A({\bm{k}}).$
(75)

Here $A({\bm{k}})$ and $B({\bm{k}})$ are $9\times 9$ matrices in the basis of
sublattice $s=A,B,C$ and $i=1,2,3$ representing three directions perpendicular
to $\psi^{0}_{s}(\alpha,\theta)$, and fulfill the condition
$A({\bm{k}})=A^{\dagger}({\bm{k}})$ and $B({\bm{k}})=B^{T}(-{\bm{k}})$,
respectively. The diagonalization of the bosonic BdG Hamiltonian, Eq. (75),
should be taken with special care. To fulfill the bosonic commutation relation
for the eigenmodes, the transformation matrix $T({\bm{k}})$ for the
diagonalization has to satisfy the paraunitarity condition
$T^{\dagger}({\bm{k}})\Sigma^{z}T({\bm{k}})=T({\bm{k}})\Sigma^{z}T^{\dagger}({\bm{k}})=\Sigma^{z}$
with the third Pauli matrix $\Sigma^{z}$ acting on the Nambu space. As a
consequence, $T({\bm{k}})$ and the corresponding eigen energies are obtained
from diagonalizing the matrix $\Sigma^{z}H({\bm{k}})$ instead [64]. Using this
diagonalization scheme, one obtains

$\displaystyle
H_{\textrm{BdG}}=\frac{1}{2}\sum_{{\bm{k}}}\begin{pmatrix}\hat{\gamma}_{{\bm{k}}}^{\dagger}&\hat{\gamma}_{-{\bm{k}}}\end{pmatrix}\begin{pmatrix}E({\bm{k}})&0\\\
0&E(-{\bm{k}})\end{pmatrix}\begin{pmatrix}\hat{\gamma}_{{\bm{k}}}\\\
\hat{\gamma}_{-{\bm{k}}}^{\dagger}\end{pmatrix}-\frac{1}{2}\sum_{{\bm{k}}}\textrm{Tr}A({\bm{k}})=\sum_{{\bm{k}}}\sum_{n=1}^{9}E_{{\bm{k}},n}(\alpha,\theta)\left(\hat{\gamma}^{\dagger}_{{\bm{k}},n}\hat{\gamma}_{{\bm{k}},n}+\frac{1}{2}\right)-\frac{1}{2}\sum_{{\bm{k}}}\textrm{Tr}A({\bm{k}}).$
(76)

Here $E({\bm{k}})$ is a diagonal matrix with positive elements
$E_{{\bm{k}},n}>0$. In the ground state the occupation of the finite-energy
states vanishes and therefore the ground-state energy is given by

$\displaystyle E_{0}\approx
E_{\text{MF}}+\frac{1}{2}\sum_{{\bm{k}},n}E_{{\bm{k}},n}(\alpha,\theta)-\frac{1}{2}\sum_{{\bm{k}}}\textrm{Tr}A({\bm{k}})$
(77)

where $E_{\text{MF}}$ is the mean-field energy. It turns out that
$\textrm{Tr}A({\bm{k}})$ is independent on the expansion point
$(\alpha,\theta)$. $E_{0}$ also includes corrections arising from quantum
fluctuations, computed to leading order in $1/M$. Similarly, the free energy
at low temperatures can be approximated by

$\displaystyle F\approx
E_{\text{MF}}+\frac{1}{2}\sum_{{\bm{k}},n}E_{{\bm{k}},n}(\alpha,\theta)-\frac{1}{2}\sum_{{\bm{k}}}\textrm{Tr}A({\bm{k}})+\frac{1}{\beta}\sum_{{\bm{k}},n}\log(1-e^{-\beta
E_{{\bm{k}},n}(\alpha,\theta)}).$ (78)

Figure 9:  Spin-valley excitation spectra for an expansion around (a) a state
with ferromagnetic spin order ($\theta=0$), (b) a state with chiral spin
chirality ($\theta=\pi/4$) and (c) a coplanar spin-state ($\theta=\pi/2$). The
spectra are drawn along the red line ($M-\Gamma-M$) in the Brillouin zone
shown in the inset of panel (a). The excitation spectra differ in the number
of Goldstone modes. Panel (a): four quadratic modes, (more precisely, for each
valley, a Goldstone mode corresponding to the ferromagnetic order and a
$\theta$-related soft mode) and a linear mode, Panel (b): three quadratic
modes and a linear mode, Panel (c): two quadratic and three linear modes.

The spectrum of excitation $E_{{\bm{k}},n}(\alpha,\theta)$ depends on the
chosen mean-field state $\psi^{0}_{s}(\alpha,\theta)$ even in cases where the
mean-field energy is exactly the same. This is shown in Fig. 9 where the
excitation spectrum is shown for an expansion around (i) a state with
ferromagnetic spin-order ($\theta=0$), (ii) a state with chiral spin chirality
($\theta=\pi/4$) and (iii) a coplanar spin-state ($\theta=\pi/2$). The
excitation spectrum differs in the number of Goldstone modes and also in its
high-energy spectrum. Thus both $E_{0}$ and $F$ will depend on the chosen
ground state within the mean-field ground-state manifold. Nature will select
the state with the lowest (free-) energy. This is an example of the “order by
disorder” mechanism, where quantum or thermal fluctuations select one specific
ordered state out of a larger manifold.

To be able to compare with the result of classical Monte Carlo calculations,
it is useful to evaluate the free energy (78) in the classical limit,
$E_{{\bm{k}},n}\ll T$, where we obtain

$\displaystyle F^{\text{cl}}\approx E_{\text{MF}}+T\sum_{{\bm{k}},n}\log(\beta
E_{{\bm{k}},n}(\alpha,\theta)).$ (79)

Figure 10(a-b) shows the free energy with different ground states,
characterized by the opening angle $\theta$. As shown in Fig 10(a), the free
energy at zero temperature has a minimum at the state with the opening angle
$\theta=\pi/2$, i.e., the spin-coplanar 120 degree order in one valley
(keeping the ferromagnetic order in the other valley). The selection of the
states are achieved by the quantum order-by-disorder mechanism. In contrast,
the thermal fluctuations select distinct states from the quantum fluctuations
as shown from the classical free energy in 10(b). The thermal order-by-
disorder mechanism leads to a selection of the spin ferromagnetic order in
both of the valley sectors ($\theta=0$ or $\pi$). Thus, our system is one of
the rare cases where quantum and classical fluctuations select very different
types of ground states. Technically this arises, because the classical
fluctuations select the state where the geometric average (sum of logarithms)
of the energies $E_{{\bm{k}},n}$ is lowest, while quantum fluctuations select
the state with the lowest arithmetic average.

In Fig. 10(c), we show the free energy, Eq. (78), as function of temperature
for the states with $\theta=\pi/2$ and $\theta=0,\pi$. Formally, the
calculation predicts a first order transition from the spin-coplanar state to
the spin-ferromagnetic state upon increasing $T$. The transition temperature,
$T\approx 3.8\,J_{2}$, is, however, so high that the expansion around the
$T=0$ mean-field, which underlies Eq. (78), is not expected to be valid any
more. Our classical Monte-Carlo simulations (see main text) show that there is
no long-ranged order at this temperature.

Figure 10: Order-by-disorder mechanisms. (a) The ground-state energy per site
(Eq. (77)) as function of $\theta$, parameterizing the opening angle of
valley-projected spins in the three sublattices. By the quantum order-by-
disorder mechanism, the state with the opening angle $\theta=\pi/2$ (i.e., the
spin-coplanar 120 degree order) is selected from the ground-state manifold.
(b) In the classical model, the free energy $F^{\text{cl}}$, Eq. (79), obtains
at low $T$ a correction linear in $T$ from thermal fluctuations. In contrast
with the quantum order-by-disorder mechanism, the thermal order-by-disorder
leads to the selection of the spin ferromagnetic order in both of the valley
sectors ($\theta=0$ or $\pi$). (c) Free energy as function of temperature. As
the thermal order-by-disorder mechanism selects the distinct ground states
from the quantum order-by-diorder mechanism, a first-order phase transition
occurs at an intermediate temperature, $k_{B}T\sim 3.8J_{2}$.

The discussion of the free energy given above explains that the classical
Monte Carlo calculations reported in the main text obtain a spin-ferromagnetic
ground state ($\theta=0,\pi$) in the limit $T\to 0$. For the classical model,
the spin-wave theory should become exact for low $T$, as it captures Gaussian
fluctuations around the classical ground state. We can therefore use it to
explain two numerical results shown in Fig. 4(b) and (c). The spin-chirality,
$\chi^{+}-\chi^{-}$, obtains a finite expectation value at $T>0$, which rises
in a singular way as function of temperature. Furthermore, the ferromagnetic
order parameter obtains a correction linear in $T$ with a prefactor which
increases for increasing system size. Remarkably, such a system-size
dependence is largely absent for the spin-chirality and only visible at the
lowest temperatures, see inset of Fig. 4(b).

The suppression of the order parameter is a well-known consequence of the
Mermin-Wagner theorem: the thermal occupation of the Goldstone modes gives
rise to a correction of order $-T\log{1/N}$ to the order parameter, where $N$
is the linear system size. This effect is clearly visible in the numerics. In
an infinite system, $N\to\infty$, long-range order is expected to be absent at
any finite $T$ with a correlation length which is exponentially large in
$1/T$.

More surprising is the finite chirality and its unusual temperature
dependence. It is straighforward to expand the valley-projected chirality
operators

$\displaystyle\chi^{\pm}=\langle{\bm{\sigma}}_{m_{1}}P^{\pm}_{m_{1}}\cdot({\bm{\sigma}}_{m_{2}}P^{\pm}_{m_{2}}{\bm{\times}}{\bm{\sigma}}_{m_{3}}P^{\pm}_{m_{3}})\rangle\approx\chi^{\pm}_{\text{cl}}+\frac{1}{V}\sum_{{\bm{k}},n}\chi^{\pm}_{{\bm{k}},n}(\,\langle\hat{\gamma}^{\dagger}_{{\bm{k}},n}\hat{\gamma}_{{\bm{k}},n}\rangle+1/2)$
(80)

where $V=N^{2}$ is the number of sites in the system, $\chi^{\pm}_{\text{cl}}$
is the chirality of the mean-field ground state and $\chi^{\pm}_{{\bm{k}},n}$
is a numerically determined weight factor which encodes how much chirality an
excitation with quantum numbers $n$ and ${\bm{k}}$ carries.

To compare to our classical Monte Carlo simulation, we use this formula
expanding around the spin-ferromagnetic state, $\theta=0,\pi$, where
$\chi^{\pm}_{\text{cl}}=0$. In the classical limit, $T\gg E_{{\bm{k}},n}$, we
have to replace
$\langle\hat{\gamma}^{\dagger}_{{\bm{k}},n}\hat{\gamma}_{{\bm{k}},n}\rangle+\frac{1}{2}$
by $T/E_{{\bm{k}},n}$ and we obtain

$\displaystyle\langle\hat{\chi}^{\pm}\rangle=\frac{1}{V}\sum_{{\bm{k}},n}\chi^{\pm}_{{\bm{k}},n}\frac{T}{E_{{\bm{k}},n}(\alpha,\theta)}.$
(81)

Importantly, $\chi^{\pm}_{{\bm{k}},n}$ turns out to be finite for ${\bm{k}}\to
0$ for one of the modes, which we label by $n_{\pm}$, with
$E_{{\bm{k}},n_{\pm}}\approx c_{\theta}k^{2}$ for ${\bm{k}}\to 0$.
Numerically, we obtain $\chi^{\pm}_{0,n_{\pm}}=\mp 2.60$ (when expanding
around $\theta=0$, signs are opposite when expanding around $\theta=\pi$) and
$c_{\theta}\approx 3.1$. The $n_{\pm}$ modes describe fluctuations of $\theta$
which naturally give rise to a finite spin chirality (note that
$\chi^{\pm}\approx\mp\frac{3\sqrt{3}}{16}\theta^{2}$ according to Eqs. (69)
and (70)).

Thus, Eq. (84) predicts for the $\theta=0$ state a nominally divergent
contribution to the chirality of the form

$\displaystyle\langle\hat{\chi}^{\pm}\rangle\approx\frac{1}{V}\sum_{{\bm{k}}}\chi^{\pm}_{0,n_{\pm}}\frac{T}{J_{2}c_{\theta}k^{2}}\approx\mp
0.11\,\frac{T}{J_{2}}\ln\left[\frac{k_{0}}{k_{\text{min}}}\right].$ (82)

where $k_{0}$ denotes a UV cutoff to the $k$ sum and we introduced ad hoc an
minimal momentum $k_{\text{min}}$ as an infrared cutoff. In three dimensions,
the analog calculation would give $T(c_{1}-c_{2}k_{\text{min}})$.

Thus the question arises, what sets the value of the infrared cutoff
$k_{\text{min}}$. Importantly, it is not set by the system size but by the
fact that the $\theta$-modes $n_{\pm}$ are not true Goldstone modes. While
within mean-field changes of $\theta$ in the two valley-sectors do not cost
any energy. This is, however, not the correct result. As shown in Fig. 10(b),
the free energy near the classical minima ($\theta_{0}=0$, $\pi$) is
approximately described by

$\displaystyle F_{\alpha,\theta}^{\textrm{cl}}\approx
F_{\alpha,\theta_{0}}^{\textrm{cl}}+\frac{1}{2}c_{1}T(\theta-\theta_{0})^{2}.$
(83)

This term induces a finite mass $\sim T$ to the $n_{\pm}$ modes resulting in
an effective IR cutoff $k_{\text{min}}\sim\sqrt{T}$,

$\displaystyle\langle\hat{\chi}^{+}-\hat{\chi}^{-}\rangle\approx\pm
0.22\,\frac{T}{J_{2}}\ln\left[\sqrt{T_{0}/T}\right].$ (84)

where the sign depends on whether we expand around $\theta=0$ or $\theta=\pi$
and $T_{0}$ is some UV cutoff energy. Eq. (84) explains the singular
temperature dependence observed in the Monte Carlo numerics, Fig. 4 and also
the approximate absence of finite-size effects in this quantity as long as
$k_{\text{min}}N\gtrsim 1$. In three dimensions, the analog calculation would
give a correction of the form $c_{1}T-c_{2}T^{3/2}$.

To see how the mass gap enters the dispersion of the $n_{\pm}$ modes, we add a
Hamiltonian, given by

$\displaystyle
H_{\text{mass}}=J_{\text{m}}\sum_{m,m^{\prime}}(({\bm{\sigma}}P^{+})_{m}\cdot({\bm{\sigma}}P^{+})_{m^{\prime}}+({\bm{\sigma}}P^{-})_{m}\cdot({\bm{\sigma}}P^{-})_{m^{\prime}}).$
(85)

Note that this Hamiltonian stabilizes the spin-ferromagnetic state
($\theta=0,\pi$), but opens up a mass gap in the $n^{\pm}$ modes. Expanding
around the $\theta=0$ state, one of the classical minima with
$H_{0}+H_{\text{mass}}$, we obtain the excitation spectra shown in Fig. 11.
Each of the $n_{\pm}$ modes acquires a mass gap with dispersion
$E_{{\bm{k}},n_{\pm}}\sim{\bm{k}}^{2}+J_{\text{m}}^{2}$. Therefore, the mass
gap nicely provides an effective infrared cutoff $k_{\text{min}}\sim
J_{\text{m}}$, as discussed above.

Figure 11: Mass gap of $\theta$-modes, $n_{\pm}$ modes (see the text), induced
by the Hamiltonian Eq. (85) to stabilize the ferromagnetic spin order. The
excitation spectra are obtained from an expansion of a ferromagnetic state
($\theta=0$), which is an energy minimum of the mean-field theory. The two
$n^{\pm}$ modes from two valleys acquire the mass gap as shown in near
$\Gamma$ point (the red circle).

The discussion given above applies to the classical model. In the quantum
case, when we expand around a ferromagnetic solution (stabilized, e.g., by
$J_{2}^{\prime}<0$), $|\langle\hat{\chi}^{+}-\hat{\chi}^{-}\rangle|$ is finite
even for $T=0$. In the spin-planar phase, $\theta=\pi/2$ (stabilized by
quantum fluctuations in the pure $J_{2}$ model and also obtained for
$J_{2}^{\prime}>0$), in contrast, $\chi^{\pm}$ vanishes by symmetry, see
below.

Two symmetries are most important for the discussion of the spin chirality.
First, a 180∘ rotation of both spin and space around, e.g., the $\hat{y}$ axis
maps $\chi^{\pm}$ to $-\chi^{\pm}$. Second, the inversion symmetry maps
$\chi^{\pm}$ to $-\chi^{\mp}$.

The spin-ferromagnetic state (with a 120∘ valley order) breaks 180∘ rotation
symmetry but is inversion symmetric. Thus $\chi^{+}+\chi^{-}=0$ while
$\chi^{+}-\chi^{-}$ is finite. The state with coplanar 120∘ spin-order
(realized in the quantum $J_{2}$ model), in contrast, has both symmetries and
thus $\chi^{\pm}=0$. This is also reflected in Fig. 5 of the main text which
shows that $|\chi^{+}-\chi^{-}|$ is only finite in the spin-ferromagnetic
phase.

##  Appendix E Semi-classical Monte Carlo

In this section, we provide a more detailed description of the our semi-
classical Monte Carlo implementation. We note that a very similar description
(by some of us) for a filling of two instead of one electron per site can be
found in Ref. [65]. We then conclude this manuscript by presenting additional
numerical data elucidating the type of phase transitions separating the
disordered and the two ordered states found in the $J_{2}-J_{2}^{\prime}$
model (see Fig. 5).

###  E.1 Implementation

To calculate finite-temperature observables, we perform semi-classical Monte
Carlo calculations using the Metropolis algorithm [58] with local updates.
Instead of a classical spin configuration, however, we need to update the
product-state wavefunction $|\Psi\rangle=\prod_{m}|\Psi_{m}\rangle$, where
$|\Psi_{m}\rangle$ is a single-site, 4-component, normalized wave function. To
this end, we parameterize the single-site wave-function as

$|\psi_{m}\rangle=\sum_{j=1}^{4}b_{m}^{j}|\gamma^{j}\rangle\,,$ (86)

with normalized, 4-dimensional, complex-valued vectors $|\mathbf{b_{m}}|=1$.
The states $|\gamma^{j}\rangle$ constitute a basis of the local Hilbert space,
for which we simply choose

$|\gamma^{j}\rangle\in\left\\{|\uparrow+\rangle,|\downarrow+\rangle,|\uparrow-\rangle,|\downarrow-\rangle\right\\}$
(87)

where $\sigma=(\uparrow,\downarrow)$ is the spin and $\alpha=(+,-)$ the valley
quantum number, labeling the eigenvalues of $\sigma^{z}$ and $\tau^{z}$,
respectively. Subtracting the normalization and a local arbitrary phase, a
state can therefore be parametrized by $N^{2}\cdot(8-2)$ real numbers. In the
Monte Carlo calculation, however, it turns out beneficial to simply include
the redundancy of the phase, which does not affect any of our observables, and
work with all components of $\mathbf{b}_{m}$. To perform a local Metropolis
update, we consequently need to be able to uniformly sample the space of
normalized, complex valued, 4-dimensional vectors. Such vectors can be
understood to live on a 7-dimensional hypersphere (7-sphere), parameterized by
the real- and imaginary part of each component. To uniformly sample on a
7-sphere, one can simply draw $7+1$ normally distributed numbers and then
normalize the resulting vector [66]. Sampling on the full sphere, however,
leads to very low acceptance rates for low temperatures, which in turn results
in a slow convergence of the results. Instead, we adapt Ref. [67] and utilize
the _Gaussian trial move_ , which generates a new local state in the
‘vicinity’ of the original as

$\mathbf{b}_{m}^{\prime}=\frac{\mathbf{b}_{m}+\sigma_{g}\mathbf{\Gamma}}{|\mathbf{b}_{m}+\sigma_{g}\mathbf{\Gamma}|},$
(88)

where $\mathbf{\Gamma}$ is a 4-dimensional complex vector, with the real and
imaginary part of each component sampled from a normal distribution. The value
of $\sigma_{g}$ controls the ‘step-size’ of the update. Staring with a large
$\sigma_{g}=60$ and then adjusting $\sigma_{g}$ every ten Monte Carlo sweeps
according to

$\sigma_{g}\to\frac{0.5}{1-R}\sigma_{g},$ (89)

where $R$ is the acceptance rate during the last ten sweeps, this very quickly
tunes the overall acceptance rate to approximately 50 % leading to significant
speedup in convergence at lower temperatures.

We begin each Monte Carlo run with a thermalization phase, typically lasting
for $N_{t}=1\cdot 10^{6}$ sweeps, in which the temperature is continuously
lowered from a large initial value of $T_{I}=3|J_{2}|$ to the desired
temperature $T$. More precisely, for the first $3/4N_{t}$ sweeps the
temperature is lowered by a multiplication with the factor
$(T/T_{I})^{\frac{4}{3N_{t}}}$ after each sweep. For the remaining
$\frac{1}{4}N_{t}$ sweeps the temperature is kept constant. During the
thermalization phase $\sigma_{g}$ is adjusted using the procedure described
above. After thermalization, we start the measurement phase, typically for
$N_{m}=10\cdot 10^{6}$ sweeps, where we keep $T$ and $\sigma_{g}$ constant and
perform measurements every tenth sweep. The statistical evaluation of the
measurements is done using the BinningAnalysis Julia Package [68].

###  E.2 Phase transitions

Figure 12: Thermal phase transitions in the $\mathbf{J_{2}-J_{2}^{\prime}}$
model. Monte Carlo data showing the specific heat $c_{v}$, energy per site
$\epsilon$ and energy distribution $p(\epsilon)$ in the three principal phase
transitions we observe, obtained for linear system sizes $N=24,36,48,72$. At
$J_{2}^{\prime}/J_{2}=-0.4$ (left column) a transition from the disordered
state to the ordered state with $120^{\circ}$ order in $\mu^{1,2}$ and
ferromagnetic order in $\bm{\sigma}P^{\pm}$ occurs. The continuous decrease in
energy and the single peak structure of the energy distribution at the
transition temperature suggest a continuous or a weak first-order phase
transition. At $J_{2}^{\prime}/J_{2}=0.06$ (middle column) the peak in $c_{v}$
at higher temperature corresponds to the same transition as for
$J_{2}^{\prime}/J_{2}=-0.4$. The second peak, which is shown in the insets,
corresponds to the transition between the two ordered states, where the spin
changes from ferromagnetic to $120^{\circ}$ order in one valley. The kink in
the energy and the bimodal structure of the energy distribution, which become
more pronounced with increasing $N$, imply a first-order transition. For
$J_{2}^{\prime}/J_{2}=0.4$ (right column) the system directly transitions from
the disordered state into the ordered state with $120^{\circ}$ spin order,
also showing signatures of a first-order transition. Figure 13: Latent heat of
the first-order phase transitions. The latent heat $\mathbf{\Delta\epsilon}$
is obtained from the Monte Carlo data at lattice size $N=72$ by fitting double
Gaussians to the energy distribution at the critical scale and calculating the
distance between the two peaks. For $0<J_{2}^{\prime}/J_{2}<0.1$ a first-order
phase transition between the two ordered states of Fig. 5 occurs, at which in
one valley the spin transitions from ferromagnetic to coplanar $120^{\circ}$
order. At $J_{2}^{\prime}>0.1J_{2}$ the disordered state directly transitions
into the state with $120^{\circ}$ spin order. Both transitions show a sizable
latent heat, indicating a strong first-order transition.

The finite-temperature phase diagram of the $J_{2}-J_{2}^{\prime}$ model shown
in Fig. 5 features three distinct phases: A disordered phase at high
temperature, a state with $120^{\circ}$ order in $\bm{\mu}^{1,2}$ and
ferromagnetic order in $\bm{\sigma}$, and a similar state where the spin
instead shows $120^{\circ}$ order in one valley. Fig. 12 shows Monte Carlo
data for all three of the corresponding phase transitions separating the
different phases.

The transition from the disordered state into the state with ferromagnetic
spin order ($J_{2}^{\prime}/J_{2}<0.1)$ features a seemingly continuous energy
as a function of temperature and the energy distribution at the transitions
shows only one Gaussian peak, indicating a continuous phase transition, a
thermal crossover or a weak first-order transition. As discussed in the main
text and shown in Fig. 8, the sharp rise of the chirality
$|\langle\chi^{+}-\chi^{-}\rangle|$, accompanied by the breaking of a discrete
$\mathbb{Z}_{2}$ symmetry, is mostly independent from $N$, strongly suggesting
a phase transition instead of a crossover. Very close to $T_{c}$, however,
$|\langle\chi^{+}-\chi^{-}\rangle|$ strongly fluctuates between different
Monte Carlo runs, even when repeating runs at fixed $T$, leading to large
statistical errors which prohibit us from determining the precise nature of
the phase transition.

In contrast, the transition separating the two ordered states
($0<J_{2}^{\prime}/J_{2}<0.1$), as well as the transition between the
disordered phase and the phase with spin $120^{\circ}$ order
($J_{2}^{\prime}/J_{2}>0.1$) show a discontinuity in the energy and a bimodal
energy distribution at the critical temperature, both becoming more pronounced
for larger lattice sizes $N$. This suggests that both phase transitions are of
first-order. To measure the strength of the first-order transitions, we obtain
the associated latent heat by fitting double Gaussians to the energy
distribution at the critical temperature and calculating the distance between
the two peaks. The resulting latent heat for both first-order transitions is
shown in Fig. 13, which exhibit a sizable latent heat of up to
$\Delta\epsilon/J_{2}\approx 0.07$ and $\Delta\epsilon/J_{2}\approx 0.08$,
respectively, indicating strong first-order transitions. When approaching the
transition into the ferromagnetic spin order, where the thermal phase
transition appears continuous, the latent heat smoothly vanishes.