# Itinerant ferromagnetism in transition metal dichalcogenides moiré
superlattices

Pawel Potasz1 Nicolás Morales-Durán2 Nai Chao Hu2 Allan H. MacDonald2
1Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus
Copernicus University, Grudziadzka 5, 87-100 Toruń, Poland  2Department of
Physics, University of Texas at Austin, Austin, Texas, 78712, USA

###### Abstract

Moiré materials are artificial crystals formed at van der Waals
heterojunctions that have emerged as a highly tunable platform that is able to
realize much of the rich quantum physics of electrons in atomic scale solids,
and in several cases even new quantum phases of matter. Here we use finite-
size exact diagonalization methods to explore the physics of single-band
itinerant electron ferromagnetism in semiconductor moiré materials. We predict
where ferromagnetism is likely to occur in triangular-lattice moiré systems,
and where it is likely to yield the highest Curie temperatures.

##  I Introduction

Moiré materials have already been established as hosts of Mott [1, 2, 3] and
topological insulators [4], a rich variety of magnetic states [5, 6, 7, 8],
and recently even fractional Chern insulators [9, 10]. Moiré materials also
provide an alternative platform for studies of itinerant electron
ferromagnetism. Ferromagnets are many-electron ground states that break time-
reversal but not translational symmetry, have finite macroscopic
magnetization, and are more common in metals than in insulators. Ferromagnetic
metals exhibit a rich variety of interesting hysteretic magneto-resistive
effects that lie at the heart of spintronics [11] and are valuable for
technology. Theoretical studies of metallic ferromagnetism in the context of
simple one-band Hubbard models [12, 13, 14, 15, 16, 17, 18], although rarely
physically realistic, have nevertheless helped provide an understanding of the
necessary conditions to stabilize such ground states in crystalline materials.
The moiré material case, in which isolated bands are common, offers the
opportunity to compare theories of single-band itinerant electron
ferromagnetism directly with experiment.

In this article we use exact diagonalization methods (ED) to explore metallic
ferromagnetism in the single-band triangular lattice moiré materials realized
in transition-metal dichalcogenide (TMD) heterobilayers [19, 20, 21, 22] such
as WSe2/MoSe2 and WSe2/WS2. We predict where ferromagnetism is most likely to
occur and where ferromagnetic transition temperatures are maximized. The
restriction of our study to the case in which a single band is partially
occupied and well separated from other bands [23] is motivated by a technical
consideration, namely the need to restrict the dimensions of the many-electron
Hilbert spaces studied to manageable sizes 111Hilbert space dimensions grow
exponentially with the number of retained bands.. Metallic ferromagnetism is
interesting in both single-band and multi-band systems In the multi-band case
local moments from one subset of bands that supply local Hunds magnetism can
combine with large spin-stiffnesses supplied by another set of band that
validate simple mean-field descriptions - using using density functional
theory for true atomic scale materials. In contrast, single-band systems are
often more difficult to understand, requiring non-perturbative approaches as
the one we take here. Although it seems likely that the highest ferromagnetic
transition temperatures that can be realized in moiré systems are in multi-
band systems 222 The highest critical temperature metallic ferromagnets in
atomic scale crystals are after all multi-band systems. we nevertheless
anticipate that scientific progress can be achieved by comparisons between
theory and experiment across a broad range of band filling factors and band
widths in the single isolated-band regime.

Our paper is organized as follows. In Section II we specify the model that we
study \- a triangular lattice moiré material model with the Hilbert space
truncated to the lowest energy moiré band and interaction matrix elements
calculated exactly. In Section III we present our numerical results. We
examine three different ferromagnetism indicators that are available from
finite-size ED calculations; i) ground state spin quantum numbers, ii) magnon
energy estimates from the total-momentum dependence of the low-energy many-
body excitation spectrum and iii) Lanczos spin-susceptibility calculations.
All are consistent with the notion that ferromagnetism occurs when the band
filling factor of the lowest energy hole miniband is around $\nu\sim 3/4$. We
estimate that Curie temperatures that can reach $T\sim 10$ K. Finally in
Section IV we summarize and discuss our findings, estimating conditions for
which the single-band model is realistic. We conclude that the single-band
approximation is not applicable at $\nu\sim 3/4$ in the TMD moiré materials
studied experimentally to date, but that it can be realized by choosing
systems with the strongest possible moiré potentials and maximizing background
screening of the Coulomb interaction.

##  II Finite Size Moiré Material Model

In this paper we will focus on transition metal dichalcogenide heterobilayer
moiré materials [19] in which the topmost valence miniband is energetically
isolated, so that holes only populate this band upon doping. Because we are
interested mainly in understanding where ferromagnetism has a substantial
ordering temperature, we focus on the range of twist angles for which the
topmost band is relatively dispersive. The single-particle part of the
continuum model Hamiltonian describing these systems is [19]

$\displaystyle H_{0}$
$\displaystyle=-\frac{\hbar^{2}}{2m^{*}}{\bm{k}}^{2}+\Delta({\bm{r}}),$ (1)
$\displaystyle\Delta({\bm{r}})=$ $\displaystyle
2V_{m}\sum_{j=1,3,5}\cos({\bm{b}}_{j}\cdot{\bm{r}}+\psi).$ (2)

where the ${\bm{b}}_{j}$ are members of the first shell of moiré reciprocal
lattice vectors and $m^{*}$, $V_{{\rm m}}$ and $\psi$ are heterojunction
specific parameters. The specific calculations we report on below take
effective mass $m^{*}=0.35\,m_{0}$, where $m_{0}$ is the rest mass of the
electron, moiré modulation strength $V_{{\rm m}}=25$ meV, and moiré potential
shape parameter [20] $\psi=-94^{\circ}$. These numerical values correspond to
WSe2/MoSe2 heterobilayer moirés [19]. It is known [26, 27] that strain
relaxation of the moiré pattern strengthens the moiré modulation potential, an
effect that can be incorporated approximately simply by increasing the value
of $V_{{\rm m}}$. For this reason we take a slightly larger value for the
moiré modulation than the one reported for the unstrained bilayer [19].
(Approximate scaling relations relating our results to those at larger values
of $V_{{\rm m}}$ are explained in the discussion section.)

Figure 1: (a) Particle-hole transformed (hole-picture) bandstruture of moiré
TMD heterobilayers at twist angle $\theta=3.0$, moiré modulation strength
$V_{{\rm m}}=25$ meV and shape parameter $\psi=-94^{\circ}$. Note that the
lowest energy hole miniband is partially occupied and isolated from the remote
bands. (b) Density of states (DOS) of the lowest energy hole miniband
$\it{vs.}$ band filling $\nu$ for twist angles $\theta=2.5$, $\theta=3.0$,
$\theta=3.5$. The inset indicates the discrete momenta of a $M=16$ unit cell
finite-size system within a color scale band contour plot for the $\theta=3.0$
case. These bands have a van Hove singularity at energy $E_{\rm VH}\approx 15$
meV and band filling factor $\nu_{\rm VH}\approx 0.75$ in the thermodynamic
limit, $M\to\infty$. (c) and (d): The discrete energies of the (c) $M=16$ and
(d) $M=36$ finite size systems discussed in the text.

Figs. 1(a) and (b) illustrate the implied moiré band structures and densities-
of-states. The density-of-states maximum occurs at the energy of a saddle-
point van Hove singularity (VHS) at band filling $\nu\approx 3/4$, where
$\nu=\frac{N}{2M}$ with $N$ the number of valence band holes in the system (we
call them particles from now on), and $M$ the number of moiré unit cells. We
will find that ferromagnetism occurs when the van Hove singularity is close to
the Fermi level of the competing paramagnetic state. The position of the van
Hove singularity (VHS) shifts slightly to larger band filling factors $\nu$
with increasing twist angle. The VHS is manifested in finite size calculations
with $M$ unit cells by a bunching of discrete states in a small energy
interval. In Fig. 1(c) we show the discrete single-particle spectra of (c)
$M=16$ and (d) $M=36$ meshes. When momentum space is discrete, the
thermodynamic limit VHS results in a set of closely spaced discrete energies
slightly below $E=15$ meV. When these states are occupied only by majority
spins and all other states are doubly occupied the filling factor is
$\nu=0.72$ for $M=16$ and $\nu=0.74$ for $M=36$ system sizes, respectively.
Note that single particle states at general momenta in the Brillouin-zone
interior are six fold degenerate simply due to triangular lattice rotational
symmetries; this property is responsible for the bunching near $E=5.0$ meV for
$M=16$ and near $E=2.5$, $E=7.0$, $E=9.0$ meV for $M=36$. ($\gamma$ point
(${\bf k}=0$) states are non-degenerate and Brillouin zone corner states are
doubly-degenerate - the degeneracy between $K$ and $K^{\prime}$ points.) As is
commonly recognized, the bunching of single particles energy levels has an
impact on finite-size many-body results, and limits the types of conclusions
that can be reached. We will consider a variety of different finite size
geometries, each with a corresponding discretization of the moiré Brillouin
zone. In order to correctly capture the VHS physics we seek meshes that
neither underrepresent nor overrepresent the associated high density of states
close to $\nu=0.75$. In the appendix we discuss how we choose finite-size
geometries for the calculations discussed in the main text.

The full Hamiltonian is obtained by projecting the two-particle Coulomb
interaction term to the topmost valence band shown in Fig. 1(a):

$\displaystyle H$
$\displaystyle=\sum_{{\bm{k}}}\epsilon_{{\bm{k}}}c^{\dagger}_{{\bm{k}}\sigma}c_{{\bm{k}}\sigma}+\frac{1}{2}\sum_{\begin{subarray}{c}i,j,k,l\\\
\sigma,\sigma^{\prime}\end{subarray}}\;V_{i,j,k,l}^{\sigma,\sigma^{\prime}}\;c^{\dagger}_{{\bm{k}}_{i}\sigma}c^{\dagger}_{{\bm{k}}_{j}\sigma^{\prime}}c_{{\bm{k}}_{l}\sigma^{\prime}}c_{{\bm{k}}_{k}\sigma},$
(3)

where $c^{\dagger}_{{\bm{k}}\sigma}$ ($c_{{\bm{k}}\sigma}$) creates
(annihilates) a particle with momentum ${\bm{k}}$ and spin $\sigma$,
$\epsilon_{{\bm{k}}}$ are band energies, and the Coulomb matrix elements are
given by

$\displaystyle
V_{i,j,k,l}^{\sigma,\sigma^{\prime}}=\frac{1}{A}\sum^{{}^{\prime}}_{\begin{subarray}{c}{\bm{G}}_{i},{\bm{G}}_{j}\\\
{\bm{G}}_{k},{\bm{G}}_{l}\end{subarray}}\left(z^{*}_{{\bm{k}}_{i},{\bm{G}}_{i}}z^{*}_{{\bm{k}}_{j},{\bm{G}}_{j}}z_{{\bm{k}}_{k},{\bm{G}}_{k}}z_{{\bm{k}}_{l},{\bm{G}}_{l}}\right)\frac{2\pi
e^{2}}{\epsilon\,q},$ (4)

with $z_{{\bm{k}},{\bm{G}}}$ eigenstate coefficients obtained from
diagonalization of Hamiltonian $H_{\rm 0}$ given by Eq. (1) in a basis of
plane waves ${\bm{G}}$. In Eq. (4) $A$ is moiré unit cell area, momentum
conservation implies that matrix elements are non-zero only if
${\bm{k}}_{i}+{\bm{k}}_{j}={\bm{k}}_{k}+{\bm{k}}_{l}$ modulo a moiré
reciprocal lattice vector, the prime on the sum over the ${\bf G}$’s implies
that ${\bm{k}}_{i}+{\bf
G}_{i}+{\bm{k}}_{j}+{\bm{G}}_{j}={\bm{k}}_{k}+{\bm{G}}_{k}+{\bm{k}}_{l}+{\bm{G}}_{l}$,
and $q=|{\bm{q}}|=|{\bm{k}}_{i}+{\bf G}_{i}-{\bm{k}}_{k}-{\bm{G}}_{k}|$ is the
momentum transfer. As we have shown previously [28], by working in a Wannier
representation the matrix elements can be reexpressed in terms of a single
large parameter, the on-site Coulomb interaction $U_{\rm 0}$, and a series of
smaller parameters including non-local exchange, interaction assisted hopping,
and longer range local interactions. The strength of interactions depends on
the value used for the effective dielectric constant $\epsilon$, which
represents screening by the three-dimensional dielectric environment of the
moiré system. We return to this issue in the discussion section.

The physics of ferromagnetism is often viewed qualitatively as a competition
between band energies, which favor states with minimal spin-polarization and
interaction energies, which favor spin-polarized states because many-electron
wavefunctions must vanish when electrons with parallel spins approach each
other, thereby avoiding strong repulsive interactions. The gain in interaction
energy per unit cell is often referred to as the Stoner energy $I$. In Fig. 2
we compare finite size kinetic energies for single-Slater-determinant (SD)
states with maximal and minimal spin-polarization in triangular lattice moiré
materials, $\Delta E_{{\rm kin}}=E^{{\rm min}}_{{\rm kin}}(S^{z}_{\rm
max})-E^{{\rm min}}_{{\rm kin}}(S^{z}_{\rm min})$, where the superscripts
‘min’ emphasize that the occupation numbers are chosen to minimize the kinetic
energy subject to the spin-polarization constraint. The energy difference per
moiré period reaches its maximum when the band is half-filled because this is
the filling factor with the maximum possible spin-polarization per moiré cell.
The kinetic energy cost increases with twist angle $\theta$ because of
increasing band widths, see the appendix. Note that the kinetic energy cost of
spin-alignment is, for the most part, reasonably well approximated at
relatively small system sizes, and that the kinetic energy cost is very small
for large band filling factors because of the VHS near the top of the first
hole miniband. This is the filling factor regime where itinerant
ferromagnetism might be expected.

Figure 2: Kinetic energy difference between maximal and minimal spin polarized
states per moiré cell $\Delta E_{{\rm kin}}/M=(E^{{\rm min}}_{{\rm
kin}}(S^{z}_{\rm max})-E^{{\rm min}}_{{\rm kin}}(S^{z}_{\rm min}))/M$ as a
function of band filling factor $\nu=N/2M$ for $\theta=3.0$.

##  III Exact Diagonalization results

We will discuss three different indicators for ferromagnetism that are
available from finite-size calculations. First of all we consider the total
spin quantum number of the finite-size many-electron ground state. The absence
of spin-orbit coupling in our model allows a ferromagnet to be defined as a
system in which the ground state total spin quantum number $S$ is extensive.
We find that maximal spin-polarization is common in finite-size systems at
band filling factors larger than about $3/4$, and conclude that ferromagnetism
will occur through much of this filling factor range. In the following
subsections we estimate the temperature to which ferromagnetism survives in
two different ways: i) by extracting magnon-energies from the momentum
dependence of the many-body excitation spectrum and ii) by extracting finite
temperature Stoner energies $I$ from the temperature-dependent spin-
susceptibilities calculated using finite-temperature Lanczos methods.

Figure 3: The ground state total spin $S$ as a function of filling factor
$\nu$ from exact diagonalization calculations for a system with $M=N_{\rm
1}\times N_{\rm 2}=16$ unit cells ($N_{\rm 1}$ and $N_{\rm 2}$ are defined in
the Appendix). (a) Spin polarization map: Total spin as a function of filling
factor $\nu$ and dielectric constant $\epsilon$ for twist angle $\theta=2.5$
and moiré potential strength $V_{{\rm m}}=11$ meV. A horizontal blue line
labels the metal-insulator transition at half-filling [29]. (b) Comparison of
the ground state spin polarization of the moiré continuum Hamiltonian and the
corresponding on-site Hubbard model for dielectric constant
$\epsilon^{-1}=0.04$, twist angle $\theta=3.0$, moiré strength $V_{{\rm
m}}=25$ meV and moiré shape $\psi=-94^{\circ}$. A dashed line indicates the
position of the van Hove singularity for finite size mesh.

###  III.1 Ground State Spin

We first assess the tendency toward ferromagnetism by comparing ground state
energies in different total spin $S$ sectors. Typical results are summarized
in Fig. 3(a), where we plot ground state spin quantum numbers vs $\nu$ and the
interaction strength parameter $\epsilon^{-1}$. Large ground state spins
appear in several different regimes in this plot. First of all they appear at
small band filling factors and weak interactions. We view ferromagnetism in
this regime as an artifact of the symmetry-related momentum-space shell
degeneracy of the finite-size mesh used to produce these results, which we
have illustrated in Fig. 1. Secondly, ferromagnetism is seen near half-filling
of the band at large interaction strengths. The ground state at $\nu=1/2$ for
this range of interaction parameters is an interaction-induced insulator (the
blue line in Fig. 3(a) labels a metal-insulator transition), but the ground
state is ferromagnetic rather than antiferromagnetic because spatially
indirect exchange interactions ($\propto\epsilon^{-1}$) exceed
antiferromagnetic superexchange interactions ($\propto\epsilon$). The property
that Mott insulators are sometimes ferromagnetic in moiré materials has been
discussed previously [28]. Our main interest here is in the very robust
ferromagnetic states that appear near band filling $\nu=3/4$, where the ground
state is metallic. In Fig. 3(b) we plot the ground state spin vs. $\nu$ in the
moderate interaction strength regime, where non-local exchange is unimportant,
demonstrating that its value is unchanged when the interaction model is
truncated to include only the on-site Hubbard-like Coulomb interaction term.
In the appendix we show that the magnetic competition in the insulating state
at $\nu=1/2$ is shifted in favor of antiferromagnetism with increasing twist
angle, as expected since larger band widths imply stronger superexchange
interactions.

In Fig. 4 we analyze the competition between ferromagnetism and paramagnetism
by partitioning the total energy into four different contributions: kinetic
energy $E_{{\rm kin}}$, Hartree energy $E_{{\rm H}}$, Fock (exchange) energy
$E_{{\rm exch}}$, and correlation energy $E_{{\rm corr}}$. Convergence to the
thermodynamic limit is easily obtained for the first three terms, whereas the
fourth part, the correlation energy, must be estimated from finite-size
calculations and extrapolated to the thermodynamic limit. For the purposes of
the qualitative point that we wish to make in this paragraph, we define the
sum of the first three terms as the expectation value of the full Hamiltonian
in the single Slater determinant (SD) state constructed by occupying the
lowest energy single-particle states for a given spin-polarization. We define
the mean-field interaction energy difference $\Delta E_{{\rm HF}}=\Delta
E_{\rm H}+\Delta E_{\rm exch}$ between maximally and minimally spin polarized
SD states by subtracting the kinetic energy contribution to the energy
difference:

$\displaystyle\Delta E_{{\rm HF}}=\Delta E_{{\rm SD}}-\Delta E_{{\rm kin}}.$
(5)

Note that $\Delta E_{{\rm HF}}$ accounts for the fact that the shape of the
charge distribution within the unit cell is different in spin-polarized and
unpolarized state, an effect that is absent in the Hubbard model. Because of
this effect the lowest energy SD state is not always the one constructed from
the lowest energy single-particle states. In the Appendix we show results for
$\Delta E_{{\rm exch}}$ obtained from multi-band self-consistent Hartree-Fock
calculations. These energies have larger negative values because of the
additional band-mixing degrees-of-freedom that are optimized.

The correlation energy is defined as the difference between the ED ground
state energy and the lowest energy SD ground state energy in a given spin
sector with subtracted kinetic energies contributions

$\displaystyle E_{{\rm corr}}(S_{\rm max})=$ $\displaystyle E_{{\rm
tot}}(S_{\rm max})-E^{kin}_{{\rm tot}}(S_{\rm max})$ (6) $\displaystyle-$
$\displaystyle(E_{{\rm SD}}(S_{\rm max})-E^{kin}_{{\rm SD}}(S_{\rm max})).$

The correlation energy difference is

$\displaystyle\Delta E_{{\rm corr}}=E_{{\rm corr}}(S_{\rm max})-E_{{\rm
corr}}(S_{\rm min}).$ (7)

With the above definitions, the total energy difference is

$\displaystyle\Delta E_{{\rm tot}}=\Delta E_{{\rm SD}}+\Delta E_{{\rm corr}}.$
(8)

In Fig. 4 we see that mean-field interaction energies $\Delta E_{{\rm HF}}$
strongly favor spin-polarized states, and that the degree to which
interactions favor spin-polarized states is strongly reduced when correlations
are included. For the parameters of this calculation, increasing the strength
of interactions actually does not substantially increase the degree to which
interactions favor spin-polarization. This is precisely the problem in
estimating where ferromagnetism occurs; once correlations are strong,
electrons avoid each other well even if they have the same spin, and even in
metallic states. Ferromagnetism is most likely when one subset of states has a
high density-of-states so that it is easily polarized, and the remaining
states are strongly dispersive so that correlations are suppressed. Conditions
favorable for ferromagnetism are regularly achieved in multi-band systems,
like the paradigmatic late $3d$ transition metals. In single-band systems
somewhat less favorable conditions can be achieved by having a sharp maximum
in the density-of-states. For 2D materials, maxima always appear at saddle
points in the band structure. It follows that single-band ferromagnetism in 2D
moiré materials is most likely when the Fermi level of the paramagnetic state
is close to a saddle point in the band structure.

A typical result for the competition in total energy between fully spin
polarized and depolarized states is summarized in Fig. 5 where we see that
ferromagnetism is most likely near $\nu=3/4$ as expected. The Hartree-Fock
theory results for the weaker of the two interaction strengths considered tell
a cautionary tale about finite-size effects since they predict ferromagnetism
for $M=16$ finite-size systems and paramagnetism for $M=441$ finite-size
systems; the $M=16$ mesh overstates the van Hove singularity, see Fig. 2. In a
vicinity of half-filling ferromagnetism is predicted for $M=441$ but the
energy of SD state with $S=S_{\rm min}$ is not the lowest one here; instead a
state with broken translation symmetry, the three sublattice Néel state, is
expected to have lower energy and competes with FM; both of these two states
have been indeed observed in experiment [1, 6]. For stronger interactions,
ferromagnetism is predicted in a vinicty of $\nu=0.75$ for both meshes. In the
following sections, we focus on estimates of transition temperatures for
ferromagnetism around this particular filling, indicated by a black dashed
line in Fig. 3.

Figure 4: Exchange energy and correlation energy difference between maximally
$S_{\rm max}$ and minimally $S_{\rm min}$ spin polarized states normalized per
moiré unit cell. These plots are based on finite-size ED calculations for
$M=16$ and on non-self-consistent Hartree-Fock, single Slater determinant SD,
for $M=441$. A dashed line indicates the position of the van Hove singularity.
(a) and (c) for interaction strength $\epsilon^{-1}=0.04$ and (b) and (d) for
interaction strength $\epsilon^{-1}=0.1$. These plots are for twist angle
$\theta=3.0$, moiré modulation strength $V_{{\rm m}}=25$ meV, and potential
shape $\psi=-94^{\circ}$. Figure 5: The total energy difference between
maximally $S_{\rm max}$ and minimally $S_{\rm min}$ spin polarized states
$\Delta E_{{\rm tot}}=E_{{\rm tot}}(S_{\rm max})-E_{{\rm tot}}(S_{\rm min}))$
per moiré unit cell for (a) $\epsilon^{-1}=0.04$ and for (b)
$\epsilon^{-1}=0.1$. $\Delta E_{{\rm tot}}$ for $M=16$ is obtained from ED
calculations and for $M=441$ from exchange energy and extrapolated correlation
energy from ED. These results were obtained with model parameters
$\theta=3.0$, $V_{{\rm m}}=25$ meV and $\psi=-94^{\circ}$.

###  III.2 Magnon Energies

In metallic ferromagnets with large splitting between majority spin and
minority spin quasiparticle energies, the ordering temperature is typically
limited by collective thermal fluctuations. The Curie temperature then scales
with the energies of the magnon modes, just as it does in insulating magnets.
In Fig. 6(a) we show the spin-flip excitation spectrum of a typical maximally
spin-polarized state near $\nu=3/4$. We associate the 16 lowest energy
excitations (one for each momentum) with magnon modes and the higher-energy
excitations with unbound quasiparticle spin-flip excitations. We see that the
magnon energies are several times smaller than the quasiparticle spin-
splitting energy. In Fig. 6(b) we show the twist angle dependence of the
highest magnon energy, which grows with the band width, suggesting that spin-
stiffness is supplied mainly by band dispersion.

Since we neglect spin-orbit interactions, our two-dimensional model is spin-
rotationally invariant and its critical temperature therefore vanishes. We
defer to a separate study the issue of engineering strong spin-orbit
interactions in TMD triangular lattice moiré materials in order to suppress
long-wavelength thermal fluctuations. Fig. 6(b) suggests that ferromagnetic
critical temperatures approaching $100$ K could be achievable at large twist
angles for sufficiently strong spin-orbit interactions. However, it is
important to realize that the single-band approximation could fail at large
twist angles. We return to this point again in the discussion section.

Figure 6: Spin-flip excitation spectrum of a $M=N_{\rm x}\times N_{\rm y}=16$
fully polarized ground state. (a) Energy spectrum for total spin $S=S_{{\rm
max}}-1$ for the system with $N_{{\rm h}}=9$ holes ($N=23$, $\nu=0.73$),
corresponding to the filling factor $\nu$ indicated by a dashed line in Fig.
3, $\epsilon^{-1}=0.04$, $\psi=-94^{\circ}$. The 16 lowest energy excitations
can be associated with magnon collective modes, and the higher energy
excitations with unbound spin-flip particle-hole excitations. $\Delta E$
indicates the width of the magnon spectrum, which scales with the transition
temperature. $N_{\rm x}$ ($N_{\rm y}$) is the number of unit cells along two
directions determined by real space lattice vectors ${\bf a}_{1}$ and ${\bf
a}_{2}$ on a triangular lattice, $K_{\rm x}$($K_{\rm y}$) are total momenta
along two directions determined by reciprocal space lattice vectors ${\bf
b}_{1}$ and ${\bf b}_{2}$. (b) The width of magnon spectrum $\Delta E$ as a
function of a twist angle $\theta$ for the moiré superlattice and its
corresponding Hubbard model.

###  III.3 Finite Temperature Lanczos Method

One of the interesting aspects of moiré materials physics from a fundamental
point of view is that the regime in which the temperature is comparable to or
larger than the band width is experimentally accessible. In the following
paragraphs we address the temperature dependence of magnetic properties over
this wide energy interval.

For the evaluation of thermodynamic properties in the canonical ensemble, we
need to calculate thermal expectation values of relevant operators $A$:

$\displaystyle\langle A\rangle=\frac{\sum_{n=1}^{N_{st}}\langle n|e^{-\beta
H}A|n\rangle}{\sum_{n=1}^{N_{st}}\langle n|e^{-\beta H}|n\rangle},$ (9)

where $\beta=1/k_{B}T$ with $k_{B}$ the Boltzman constant, the partition
function $Z=\sum_{n=1}^{N_{st}}\langle n|e^{-\beta H}|n\rangle$, and
$|n\rangle$ is summed over orthonormal basis states. The exponential increase
of $N_{st}$ with system size places severe limits on the direct application of
these fundamental formulas.

The problem can be avoided if an appropriate statistical average of the full
Hilbert space is generated. In the finite temperature Lanczos method (FTLM)
[30] one starts with the high temperature expansion:

$\displaystyle\langle A\rangle_{\beta\rightarrow
0}=Z^{-1}\sum_{n=1}^{N_{st}}\sum_{k=0}^{\infty}\frac{(-\beta)^{k}}{k!}\langle
n|H^{k}A|n\rangle,$ (10)

where

$\displaystyle
Z=\sum_{n=1}^{N_{st}}\sum_{k=0}^{\infty}\frac{(-\beta)^{k}}{k!}\langle
n|H^{k}|n\rangle.$ (11)

The Lanczos algorithm is an iterative method for finding extreme eigenvalue of
a large matrix in which expectations of high powers of the Hamiltonian
naturally appear. During Lanczos iteration steps, a set of orthogonal basis
vectors is generated (a Krylov space), spanning a finite-size space that
contains approximations to eigenvectors corresponding to extreme eigenvalues
of a full Hilbert space with accuracy controlled by the number of iteration
steps. In the Lanczos method the Hamiltonian is diagonalized in this Krylov
space obtaining Lanczos eigenvectors $|l\rangle$ and the associated Lanczos
energy eigenvalues $\epsilon_{l}$. When the number of Lanczos steps $N_{l}\geq
k$ one can write

$\displaystyle\langle n|H^{k}A|n\rangle\approx\sum_{l=0}^{N_{L}}\langle
n|H^{k}|l(n)\rangle\langle l(n)|A|n\rangle=$
$\displaystyle\sum_{l=0}^{N_{L}}(\epsilon_{l(n)})^{k}\langle
n|l(n)\rangle\langle l(n)|A|n\rangle$ (12)

and

$\displaystyle\langle
n|H^{k}|n\rangle\approx\sum_{l=0}^{N_{L}}(\epsilon_{l(n)})^{k}|\langle
l(n)|n\rangle|^{2}.$ (13)

$N_{L}$ is a parameter of the approximation that needs to be large enough to
reach accurate extremal energy eigenvalues; for the calculations we present
below we take $N_{L}=150$. Inserting Eq. (12) and Eq. (13) into Eq. (10) and
Eq. (11) and replacing the sum over all orthonormal basis states by a much
smaller sum over $R$ random Lanczos seed states, in analogy to Monte Carlo
methods, yields

$\displaystyle\langle A\rangle\approx Z^{-1}\frac{N_{\rm st}}{N_{\rm
R}}\sum_{\nu\in N_{{\rm R}}}\sum_{l}^{N_{\rm
L}}e^{-\beta\epsilon_{l(\nu)}}\langle
l(\nu)|A|\nu\rangle\langle\nu|l(\nu)\rangle,$ (14)

where the partition function is

$\displaystyle Z\approx\frac{N_{st}}{N_{{\rm R}}}\sum_{\nu}^{N_{{\rm
R}}}\sum_{l}^{N_{L}}e^{-\beta\epsilon_{l(\nu)}}|\langle
l(\nu)|\nu\rangle|^{2}.$ (15)

The exponential-size Hilbert space of the Hamiltonian is thereby approximated
by its spectral representation in a Krylov space spanned by the $N_{\rm L}$
Lanczos vectors starting from each random vector. The chosen random vectors
$|\nu\rangle$ should ideally be mutually orthogonal, but for practical
purposes this is not really necessary since two vectors with random components
in a large dimensional space are always nearly orthogonal.

Figure 7: Comparison of FTLM and exact susceptibility calculations for $M=16$
moiré unit cells, $N=4$ electrons, and different number of random vectors
$N_{{\rm R}}$. (a) $\chi/\beta$ as a function of inverse temperature $\beta$.
The inset shows low energy many-body spectrum with total spin indicated by
color. (b) susceptibility $\chi(T)$ and (c) inverse susceptibility
$\chi^{-1}(T)$ as a function of temperature $T$. The blue line in (c) is a
linear fit to estimate a transition temperature $T_{C}$. The number of Lanczos
steps is taken to be $N_{\rm L}=150$. The size of the Hilbert space for
$S_{z}=0$ and fixed total momentum ${\bf K}$ is weakly momentum dependent and
around 900. The parameters for this illustration are interaction strength
$\epsilon^{-1}=0.04$, twist angle $\theta=3.0$, $V_{{\rm m}}=25$ meV, and
$\psi=-94^{\circ}$. $N_{{\rm R}}=20,10,5$ means that for $S_{\rm z}=0$ we take
$N_{{\rm R}}=20$, $S_{\rm z}=\pm 1$ we take $N_{{\rm R}}=10$ and so on. If one
numbers is given, for all subspaces we take the same $N_{{\rm R}}$

.

In general calculations using this approach are less sensitive to finite size
effects as temperature increases, and most sensitive to finite size at $T=0$.
This property is related to the fact that at $T=0$ both static and dynamical
quantities are calculated from one eigenstate only, and the selection of this
state can be dependent on the size and on the shape of the finite-size system.
$T>0$ introduces thermodynamic averaging over a larger number of eigenstates
and this directly reduces finite-size effects for static quantities.
Calculational efficiency can be improved by taking symmetries into account, so
that $N_{st}$ corresponds to the number of states with a given symmetry.

In our view, the finite temperature Lanczos method (FTLM) is ideally suited to
exploring the high-temperature physics that is observable for the first time
in moiré materials. In this work we focus on calculations of the spin magnetic
susceptibilty $\chi=\beta\langle S_{z}^{2}\rangle$ where

$\langle
S_{z}^{2}\rangle=\frac{\sum_{n}\exp{(-\beta\epsilon_{n})}S_{z}(n)^{2}}{\sum_{n}\exp{(-\beta\epsilon_{n})}}.$
(16)

Beacuse $[H,S_{z}]=0$, the Lanczos method can be applied to each $S_{z}$
sector separately. The FTLM formula for the susceptibilty is

$\displaystyle\chi=Z^{-1}\sum_{s}\frac{N_{\rm st}(s)}{N_{{\rm
R}}(s)}\sum_{\nu=1}^{N_{{\rm
R}}(s)}\sum_{l=1}^{N_{L}}e^{-\beta\epsilon_{l(\nu)}}|\langle
l(\nu)|\nu\rangle|^{2}s^{2},$ (17)

where $s$ is the $S_{z}$ value for the subspace. We find that the most
accurate results are obtained for $N_{{\rm R}}(s)$ chosen such that the ratio
between the Hilbert subspace size and the number of vectors is kept constant.

Figure 8: FTLM susceptibility calculations for $M=16$ moiré unit cells, $N=23$
electrons (9 holes), and different number of random vectors $N_{{\rm R}}$. (a)
$\chi/\beta$ normalized per moiré cell as a function of inverse temperature
$\beta$. Analytically calculated values are given for the limit of infinite
temperature (left) and zero temperature (right). (b) susceptibility $\chi(T)$
and (c) inverse susceptibility $\chi^{-1}(T)$ as a function of temperature
$T$. The red line is a linear fit to estimate the transition temperature
$T_{C}$. The inset shows the Stoner parameter
$I(T)\equiv\chi^{-1}(T)/M-\chi_{\rm NI}^{-1}(T)/M$, where ${\rm NI}$ means
non-interacting. The number of Lanczos steps was $N_{\rm L}=150$. The size of
Hilbert space for $S_{z}=0$ at fixed total momentum ${\bf K}$ is around
$5\times 10^{5}$. $N_{{\rm R}}=20,10,5,5$ means that for $S_{\rm z}=\pm 0.5$
we take $N_{{\rm R}}=20$, $S_{\rm z}=\pm 1.5$ we take $N_{\rm R}=10$ and for
all other subspaces $N_{{\rm R}}=5$. If one numbers is given, for all
subspaces we take the same $N_{{\rm R}}$. The parameters used for this
calculation are: dielectric constant $\epsilon^{-1}=0.04$, twist angle
$\theta=3.0$, $V_{{\rm m}}=25$ meV, and $\psi=-94^{\circ}$.

The accuracy of FTLM finite-size calculations is assessed in Fig. 7 by
comparing $\chi$ as calculated by performing the full sum over all states with
the FTLM sum. The three plots in Fig. 7 ((a) $\chi/\beta$ as a function of
inverse temperature $\beta$, (b) the susceptibility $\chi(T)$ and the (c)
inverse susceptibility $\chi^{-1}(T)$ as a function of temperature $T$)
emphasize different aspects of the temperature dependence of $\chi$. The
$\beta\to 0$ ($T\to\infty$) and $\beta\to\infty$ ($T\to 0$) limits of
$\chi/\beta$ can be calculated analytically by averaging $S_{z}^{2}$ over the
full Hilbert space ($\chi^{-1}/\beta\to M\nu(1-\nu)/2$ for $\beta\to 0$) and
over the ground state spin multiplet ($\chi^{-1}/\beta\to S(S+1)/3$ for
$\beta\to\infty$) respectively. For the test case ($M=16$ and $N=4$)
illustrated in Fig. 7, the susceptibility can be calculated exactly from the
full many-body spectrum because the Hilbert space dimension for a given
$S_{z}$ subspace does not exceed 1000. The exact result is indicated by a red
line in Fig. 7, and compared with FTLM estimates based on different numbers of
random vectors $N_{{\rm R}}$. All lines overlap for temperatures $T>20$ K,
demonstrating the high accuracy of the method in the high temperature limit.
The ground state of the system in this case has $S=0$ (see the inset), which
leads to vanishing susceptibility in (b) and divergence of the inverse
susceptibility in (c) as $T\to 0$. The susceptibility reaches a maximum at
around $T=4$ K. The blue line in (c) is a high-temperature linear fit that
extrapolates to a finite value for $T=0$, consistent with a paramagnetic
state.

The FTLM estimates have the advantage that they can be drawn from larger
Hilbert spaces. In Fig. 8 we show a typical result obtained for $N=23$
particles ($N_{\rm h}=9$ holes) in the $M=16$ case, in the regime of filling
factors where ferromagnetism is expected on the basis of the many-body ground
state calculations. In this case the many-body ground state has non-zero total
spin $S=9/2$. The exact value of $\chi/M\beta$ normalized per moiré unit cell
in the $\beta\rightarrow\infty$ limit is therefore $0.515$, as indicated by a
black arrow in Fig. 8(a). (The $\beta\to 0$ limit $0.10433$, which is
independent of interactions, is also indicated by a black arrow.) We see that
the FTLM method gives accurate results in both limits, irrespective of
$N_{{\rm R}}$ at $T\rightarrow\infty$ and for $N_{{\rm R}}\geq 10$ at $T=0$.
Generally speaking, the ratio between $N_{{\rm R}}$ and the dimension of a
given Hilbert subspace is a good accuracy indicator. The increase in $\chi$ at
intermediate temperatures relative to the high-temperature limit shows that on
average interactions lower the energies of states with larger $S_{z}$ relative
to those with smaller $S_{z}$. The linear fit to the inverse susceptibility
shown in 8(c) estimates the Curie temperature $T_{\rm C}\approx=9$ K for this
case, and the estimate is not strongly affected by $N_{{\rm R}}$ in reasonable
ranges. The inset shows the finite size Stoner parameter $I$ that has the
expected linear-in-$T$ dependence up to around $T\approx 12$ K.

Having established the efficacy of the FTLM, we now employ it to study trends
in ferromagnetism in triangular lattice moiré materials. In SM in Fig. 10(b)
and (d) we compare inverse susceptibility results for two other twist angles
for the same moiré modulation potential. Extrapolating from high temperatures
where finite-size effects are less severe, we see that the susceptibility at
higher temperatures decreases with twist angle. We attribute this decrease to
an increase in bandwidth, which decreases the Pauli susceptibility of non-
interacting electrons. At the same time, the high-temperature estimate of the
Curie temperature at which the susceptibility diverges (the inverse
susceptibility vanishes) increases with twist angle. We attribute this
increase also to increasing bandwidth, which increases magnon energies by
increasing the kinetic energy cost of spatial modulation of the magnetization.

##  IV Discussion

We have used three different indicators available from finite-size exact-
diagonalization calculations to address the physics of itinerant
ferromagnetism in single-band triangular lattice moiré materials: i) ground
state spin quantum numbers, ii) magnon excitation energies, and iii)
temperature dependent spin-susceptibilities. All indicate that ferromagnetism
is common at hole band filling factors near $\nu=3/4$ at temperatures up to
$\sim 10$K. Our calculations were performed for particular values of the moiré
modulation strength and shape parameters. These are however expected to be
strongly dependent on the specific heterojunction at which the moiré pattern
is formed, and in particular on strain relaxations at those heterojunctions
which will tend to increase modulation strengths [31, 27]. When
$V_{m}\to\lambda V_{m}$, twist angle $\theta\to\sqrt{\lambda}\theta$, and
dielectric screening parameter $\epsilon\to\sqrt{\lambda}\epsilon$, the three
terms in the continuum model Hamiltonian (interaction, moiré potential, and
kinetic energy) all increase by a factor of $\lambda$. Since the properties of
interest here are relatively insensitive to the interaction strength parameter
within reasonable ranges, it follows that the properties of systems with
stronger moiré potentials can be read off from our results by increasing
temperature scales and twist angles. In particular the larger energy scales
increase the temperatures at which ferromagnetism can occur.

It is interesting to compare TMD triangular lattice moiré materials, with
graphene multilayer moiré materials that also support ferromagnetic states. In
the latter case, it is known that because of topological obstructions
inherited from the individual layer Dirac cones [32, 33, 34], a faithful
representation of the flat moiré minibands requires multi-band tight-binding
[35] models, for which the exact diagonalization approach is not practical. In
the TMD moiré material case, however, the lowest energy moiré bands have
Wannier functions that are similar to harmonic oscillator ground states
centered on moiré potential extrema [19]. Although we do not approximate the
interaction matrix elements in our one-band model, we have verified that all
properties related to ferromagnetism are similar to those of simple triangular
lattice Hubbard models.

It is also interesting to compare TMD triangular lattice materials with
rhomohedral graphene multilayers [36, 37, 38, 39, 40, 41, 42, 43, 44], a class
of two-dimensional materials in which metallic ferromagnetism has been
discovered recently. These graphene multilayer systems are like TMD moiré
materials in that they have peaks in their densities of states, related in
that case to Liftshitz transitions of distorted Dirac cones, but they do not
have minibands and are not approximated by Hubbard models. The magnetism that
appears in these systems is consistent with the notion that the key to
ferromagnetism is a sharp density-of-states peak in a low-density-of-states
background.

The exact diagonalization method we have employed is most suitable when the
many-electron Hilbert space can be truncated to a single moiré miniband. The
small parameter which controls the applicability of this approximation is the
ratio of the largest interaction scale, the on-site Hubbard interaction
$U_{0}$, to the sub-band separation. As explained in Ref. 19 these can be
estimated by making a harmonic approximation for the moiré potential. We find
that $U_{0}\sim\rm{Ry}^{3/4}(zV_{m})^{1/4}(a_{B}/a_{M})^{1/2}$, where $z=6$ is
the triangular lattice coordination number and $\rm{Ry}$ and $a_{B}$ are the
host 2D semiconductor Rydberg energy scale $\sim 0.3$ eV and Bohr radius
length scale $\sim 1$ nm. Similarly the subband separation
$\hbar\omega\sim\rm{Ry}^{1/2}(zV_{m})^{1/2}a_{B}/a_{M}$. It follows that

$\frac{U_{0}}{\hbar\omega}\sim(\rm{Ry}/zV_{m})^{1/4}(a_{M}/a_{B})^{1/2}.$ (18)

Truncation to the lowest moiré band is justified at all band filling factors
$\nu\in(0,1)$ when the right hand side of Eq. 18 is smaller than $\sim 1$.
Most systems [45, 46] that have been studied to date do not satisfy this
criterion. Since continuum model approximations are valid only for
$a_{M}\gtrsim a_{B}$, it follows that single-band ferromagnetism will occur
only when the first factor on the right side of Eq. 18 is made small, for
example by increasing the dielectric screening environment of the moiré system
to decrease ${\rm Ry}$, or by choosing a system with a particularly large
value of $V_{m}$. From exponentially localized Wannier functions obtained for
the topmost valence band used in our calculations, we get, for $\theta=3.0$,
$U_{\rm 0}\epsilon\approx 1121$ meV [28], $\hbar\omega\approx 58.5$ meV. For
$\epsilon^{-1}=0.1$, $\frac{U_{0}(\epsilon^{-1}=0.1)}{\hbar\omega}>1$, while
for $\frac{U_{0}(\epsilon^{-1}=0.04)}{\hbar\omega}<1$. Thus for the limit of
weaker interaction strength the single band approximation is justified. This
suggest that our predictions are relevant for systems with sufficiently close
nearby gates.

We note that Coulomb repulsion will increase the energy of the lowest energy
hole miniband, as it is filled, by more than it increases the energies of
states in higher energy moiré minibands. For this reason the regime of
parameter space in which occupation of higher energy minibands can be
neglected decreases as band filling factor increases. When correlations are
included, the ground state at hole filling factor $\nu=1/2$ is often an
insulator. When its lowest energy hole-charged excitation is dominantly in a
higher hole miniband, the insulator is referred to as a charge transfer
insulator [45, 46]. Since single-band ferromagnetism is most likely near band-
filling factor $\nu=3/4$, the present single-band study is never relevant when
the ground state of the half-filled band is a charge transfer insulator, which
already involves higher energy subbands in an essential way. If systems could
be realized in which the sign of $V_{m}$ is reversed (or equivalently
$\psi\to\psi+180^{\circ}$), ferromagnetism would be expected for minibands
that are less than half-filled. For the standard sign of $V_{m}$ however, any
ferromagnetism that occurs when the interaction parameter that is the subject
of Eq. 18 is large must be of multi-band character. We leave the analysis of
this situation for a future study, for it requires a different approach.

The authors acknowledge helpful interactions with L. Fu and Y. Zhang. This
work was supported by the U.S. Department of Energy, Office of Science, Basic
Energy Sciences, under Award $\\#$ DE-SC0022106. PP acknowledges support from
the Polish National Science Centre based on Decision No. 2021/41/B/ST3/03322.
We acknowledge the Texas Advanced Computing Center (TACC) at The University of
Texas at Austin for providing high-performance computer resources.

## References

  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Simulation of hubbard model physics in wse2/ws2 moiré superlattices, Nature 579, 353 (2020).
  * Li  _et al._ [2021a] T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, J. Shan, and K. F. Mak, Continuous mott transition in semiconductor moiré superlattices, Nature 597, 350 (2021a).
  * Ghiotto  _et al._ [2021] A. Ghiotto, E.-M. Shih, G. S. S. G. Pereira, D. A. Rhodes, B. Kim, J. Zang, A. J. Millis, K. Watanabe, T. Taniguchi, J. C. Hone, L. Wang, C. R. Dean, and A. N. Pasupathy, Quantum criticality in twisted transition metal dichalcogenides, Nature 597, 345 (2021).
  * Li  _et al._ [2021b] 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 (2021b).
  * Tang  _et al._ [2023] Y. Tang, K. Su, L. Li, Y. Xu, S. Liu, K. Watanabe, T. Taniguchi, J. Hone, C.-M. Jian, C. Xu, K. F. Mak, and J. Shan, Evidence of frustrated magnetic interactions in a wigner–mott insulator, Nature Nanotechnology 18, 233 (2023).
  * Anderson  _et al._ [2023] E. Anderson, F.-R. Fan, J. Cai, W. Holtzmann, T. Taniguchi, K. Watanabe, D. Xiao, W. Yao, and X. Xu, Programming correlated magnetic states with gate-controlled moiré geometry, Science 381, 325 (2023).
  * Ciorciaro  _et al._ [2023] L. Ciorciaro, T. Smolenski, I. Morera, N. Kiper, S. Hiestand, M. Kroner, Y. Zhang, K. Watanabe, T. Taniguchi, E. Demler, and A. Imamoglu, Kinetic magnetism in triangular moiré materials (2023), arXiv:2305.02150 [cond-mat.str-el] .
  * Tao  _et al._ [2023] Z. Tao, W. Zhao, B. Shen, P. Knüppel, K. Watanabe, T. Taniguchi, J. Shan, and K. F. Mak, Observation of spin polarons in a frustrated moiré hubbard system (2023), arXiv:2307.12205 [cond-mat.str-el] .
  * Cai  _et al._ [2023] J. Cai, E. Anderson, C. Wang, X. Zhang, X. Liu, W. Holtzmann, Y. Zhang, F. Fan, T. Taniguchi, K. Watanabe, Y. Ran, T. Cao, L. Fu, D. Xiao, W. Yao, and X. Xu, Signatures of fractional quantum anomalous hall states in twisted MoTe2, Nature 10.1038/s41586-023-06289-w (2023).
  * Zeng  _et al._ [2023] Y. Zeng, Z. Xia, K. Kang, J. Zhu, P. Knüppel, C. Vaswani, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, Thermodynamic evidence of fractional chern insulator in moiré MoTe2, Nature 10.1038/s41586-023-06452-3 (2023).
  * Bader and Parkin [2010] S. D. Bader and S. Parkin, Spintronics, Annu. Rev. Condens. Matter Phys. 1, 71 (2010).
  * Tasaki [1998] H. Tasaki, From Nagaoka’s Ferromagnetism to Flat-Band Ferromagnetism and Beyond: An Introduction to Ferromagnetism in the Hubbard Model, Progress of Theoretical Physics 99, 489 (1998), https://academic.oup.com/ptp/article-pdf/99/4/489/5209512/99-4-489.pdf .
  * Hanisch  _et al._ [1997] T. Hanisch, G. S. Uhrig, and E. Müller-Hartmann, Lattice dependence of saturated ferromagnetism in the hubbard model, Phys. Rev. B 56, 13960 (1997).
  * Held and Vollhardt [1998] K. Held and D. Vollhardt, Microscopic conditions favoring itinerant ferromagnetism: Hund’s rule coupling and orbital degeneracy, The European Physical Journal B - Condensed Matter and Complex Systems 5, 473 (1998).
  * Obermeier  _et al._ [1997] T. Obermeier, T. Pruschke, and J. Keller, Ferromagnetism in the large-$u$ hubbard model, Phys. Rev. B 56, R8479 (1997).
  * Kamogawa  _et al._ [2019] Y. Kamogawa, J. Nasu, and A. Koga, Ferromagnetic instability for the single-band hubbard model in the strong-coupling regime, Phys. Rev. B 99, 235107 (2019).
  * Lebrat  _et al._ [2023] M. Lebrat, M. Xu, L. H. Kendrick, A. Kale, Y. Gang, P. Seetharaman, I. Morera, E. Khatami, E. Demler, and M. Greiner, Observation of nagaoka polarons in a fermi-hubbard quantum simulator (2023), arXiv:2308.12269 [cond-mat.quant-gas] .
  * Prichard  _et al._ [2023] M. L. Prichard, B. M. Spar, I. Morera, E. Demler, Z. Z. Yan, and W. S. Bakr, Directly imaging spin polarons in a kinetically frustrated hubbard system (2023), arXiv:2308.12951 [cond-mat.quant-gas] .
  * Wu  _et al._ [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard model physics in transition metal dichalcogenide moiré bands, Phys. Rev. Lett. 121, 026402 (2018).
  * Pan  _et al._ [2020] H. Pan, F. Wu, and S. Das Sarma, Band topology, hubbard model, heisenberg model, and dzyaloshinskii-moriya interaction in twisted bilayer ${\mathrm{wse}}_{2}$, Phys. Rev. Research 2, 033087 (2020).
  * 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).
  * Kennes  _et al._ [2021] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio, Moiré heterostructures as a condensed-matter quantum simulator, Nature Physics 17, 155 (2021).
  * Morales-Durán  _et al._ [2023] N. Morales-Durán, P. Potasz, and A. H. MacDonald, Magnetism and quantum melting in moiré-material wigner crystals, Phys. Rev. B 107, 235131 (2023).
  * Note [1] Hilbert space dimensions grow exponentially with the number of retained bands.
  * Note [2] The highest critical temperature metallic ferromagnets in atomic scale crystals are after all multi-band systems.
  * Shabani  _et al._ [2021a] S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy, Deep moiré potentials in twisted transition metal dichalcogenide bilayers, Nature Physics 17, 720 (2021a).
  * Wang  _et al._ [2023] C. Wang, X.-W. Zhang, X. Liu, Y. He, X. Xu, Y. Ran, T. Cao, and D. Xiao, Fractional chern insulator in twisted bilayer mote $\\_2$, arXiv preprint arXiv:2304.11864 (2023).
  * Morales-Durán  _et al._ [2022] N. Morales-Durán, N. C. Hu, P. Potasz, and A. H. MacDonald, Nonlocal interactions in moiré hubbard systems, Phys. Rev. Lett. 128, 217202 (2022).
  * Morales-Durán  _et al._ [2021] N. Morales-Durán, A. H. MacDonald, and P. Potasz, Metal-insulator transition in transition metal dichalcogenide heterobilayer moiré superlattices, Phys. Rev. B 103, L241110 (2021).
  * Jaklič and Prelovšek [1994] J. Jaklič and P. Prelovšek, Lanczos method for the calculation of finite-temperature quantities in correlated systems, Physical Review B 49, 5065 (1994).
  * Shabani  _et al._ [2021b] S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy, Deep moiré potentials in twisted transition metal dichalcogenide bilayers, Nature Physics 17, 720 (2021b).
  * Zou  _et al._ [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions, Phys. Rev. B 98, 085435 (2018).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, All magic angles in twisted bilayer graphene are topological, Phys. Rev. Lett. 123, 036401 (2019).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Faithful tight-binding models and fragile topology of magic-angle bilayer graphene, Phys. Rev. B 99, 195455 (2019).
  * 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).
  * Zhou  _et al._ [2022] H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and A. F. Young, Isospin magnetism and spin-polarized superconductivity in bernal bilayer graphene, Science 375, 774 (2022), https://www.science.org/doi/pdf/10.1126/science.abm8386 .
  * Zhou  _et al._ [2021a] H. Zhou, T. Xie, A. Ghazaryan, T. Holder, J. R. Ehrets, E. M. Spanton, T. Taniguchi, K. Watanabe, E. Berg, M. Serbyn, and A. F. Young, Half- and quarter-metals in rhombohedral trilayer graphene, Nature 598, 429 (2021a).
  * Zhou  _et al._ [2021b] H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Superconductivity in rhombohedral trilayer graphene, Nature 598, 434 (2021b).
  * Huang  _et al._ [2023] C. Huang, T. M. R. Wolf, W. Qin, N. Wei, I. V. Blinov, and A. H. MacDonald, Spin and orbital metallic magnetism in rhombohedral trilayer graphene, Phys. Rev. B 107, L121405 (2023).
  * Lee  _et al._ [2022] Y. Lee, S. Che, J. Velasco Jr., X. Gao, Y. Shi, D. Tran, J. Baima, F. Mauri, M. Calandra, M. Bockrath, and C. N. Lau, Gate-tunable magnetism and giant magnetoresistance in suspended rhombohedral-stacked few-layer graphene, Nano Letters 22, 5094 (2022).
  * Seiler  _et al._ [2022] A. M. Seiler, F. R. Geisenhof, F. Winterer, K. Watanabe, T. Taniguchi, T. Xu, F. Zhang, and R. T. Weitz, Quantum cascade of correlated phases in trigonally warped bilayer graphene, Nature 608, 298 (2022).
  * Han  _et al._ [2023] T. Han, Z. Lu, G. Scuri, J. Sung, J. Wang, T. Han, K. Watanabe, T. Taniguchi, H. Park, and L. Ju, Correlated insulator and chern insulators in pentalayer rhombohedral stacked graphene (2023), arXiv:2305.03151 [cond-mat.mes-hall] .
  * Jang  _et al._ [2023] Y. Jang, Y. Park, J. Jung, and H. Min, Chirality and correlations in the spontaneous spin-valley polarization of rhombohedral multilayer graphene, Phys. Rev. B 108, L041101 (2023).
  * Das and Huang [2023] M. Das and C. Huang, Unconventional metallic magnetism: Non-analyticity and sign-changing behavior of orbital magnetization in abc trilayer graphene (2023), arXiv:2308.01996 [cond-mat.mes-hall] .
  * Zhang  _et al._ [2020] Y. Zhang, N. F. Q. Yuan, and L. Fu, Moiré quantum chemistry: Charge transfer in transition metal dichalcogenide superlattices, Phys. Rev. B 102, 201115 (2020).
  * Zhang  _et al._ [2021] Y. Zhang, T. Liu, and L. Fu, Electronic structures, charge transfer, and charge order in twisted transition metal dichalcogenide bilayers, Phys. Rev. B 103, 155142 (2021).
  * Oitmaa and Betts [1978] J. Oitmaa and D. D. Betts, The ground state of two quantum models of magnetism, Canadian Journal of Physics 56, 897 (1978), https://doi.org/10.1139/p78-120 .
  * Avella and Mancini [2013] A. Avella and F. Mancini, Strongly correlated systems, Springer-Verlag https://doi.org/10.1007/978-3-642-35106-8 (2013).
  * Repellin  _et al._ [2014] C. Repellin, B. A. Bernevig, and N. Regnault, ${\mathbb{z}}_{2}$ fractional topological insulators in two dimensions, Phys. Rev. B 90, 245401 (2014).

## Supplemental materials

###  IV.1 Exact diagonalization Method

The model we study has orbital and spin degrees of freedom, discrete
triangular lattice translational symmetry, $SU(2)$ spin-rotational invariance,
and no spin-orbit coupling. The Hilbert space can be divided into smaller
subspaces with total momentum ${\bf K}$ with discrete values determined by the
number of Moire unit cells $M$, total spin ${\bf S}$, and azimuthal spin
$S^{z}$. The basis is constructed in an occupation number representation,
distributing particles among single-particle states labeled by spin $z$
azimuthal quantum number and quasi-particle $(k_{x},k_{y})$ momentum. The
total number of possible configurations $N_{\rm st}$ for particles distributed
on $M$ single particle states with a given spin $N_{\downarrow}$ or
$N_{\uparrow}$ is determined by a product of binomial coefficients, $N_{\rm
st}=\binom{M}{N_{\uparrow}}\cdot\binom{M}{N_{\uparrow}}$. The many-body
Hamiltonian projected to a given total momentum ${\bf K}$ is diagonalized in
$S^{z}$ subspaces. We do not rotate the Hamiltonian matrix to a ${\bf S}$
basis as this is an additional computational cost, and instead determine the
ground state from calculations of expectation value of total total spin $S$
for each energy eigenstate.

A high symmetry of a triangular lattice induces a six fold degeneracy of
single particle states. This affects finite-size ED results leading to
spurious results at particular fillings; we refer to spin polarization around
$\nu=0.25$ for $M=16$ as one of these effects. A typical finite-size sample is
obtained by taking a given number of Moire unit cells along two real space
primitive unit vectors and the corresponding plaquette in real space is
spanned by the vectors $N_{1}\,{\bf a_{1}}$ and $N_{1}\,{\bf a_{1}}$, where
${\bf a_{i}}$ are the real space lattice vectors. The points forming the
momentum space mesh are given by ${\bf k}=n_{1}\,{\bf b_{1}}/N_{1}+n_{2}\,{\bf
b_{2}}/N_{2}$, where ${\bf b_{i}}$ are the reciprocal lattice vectors and
$n_{i}$ are integers between $0$ and $N_{i}-1$. Such meshes sustain the
symmetry of the lattice

An alternative way of getting finite-size plaquettes is using tilted meshes
[47, 48] (see Appendix A in Ref. 49 for a nice introduction to derivation of
the tilted periodic boundary conditions constrains). In this case, momentum
space points from the interior of the Brillouin zone do not have to occur in
multiplies of six. An example of a tilted momentum space mesh with $M=24$, we
call it T24, is shown in the inset in Fig. 9(a) with single particle energies
in Fig. 9(a). Single particle states are non-degenerate, double, or four fold
degenerate. Van Hove singularity is represented by a set of states in an
energy range between $E=10$ meV and $E=13$ meV. We populate these states by up
to $N_{h}=9$ holes; we are still limited by a maximum size of a dimension of
the Hilbert space subspace. The corresponding total spin ground states for
three different twist angles are shown in Fig. 9(b). For fillings $\nu<0.9$,
finite spin polarization occurs for the smallest considered twist angle
$\theta=2.5$ for the system with odd number of particles, $N_{h}=5$,
$N_{h}=7$, and $N_{h}=9$. This is again a finite size effect, as one can see
that below the two highest energetically single particles states, there are
five closely lying states, one single and four degenerate states. In order to
minimize the total energy, one can populate these two states by two or four
holes, and singly occupied five states by particles with the same spin in
order to maximize exchange interaction. For $N_{h}=7$ and $N_{h}=9$, the
ground state total spin is $S=5/2$, which corresponds to five particles with
the same spin on these five closely lying degenerate states. Above these
states, there are two sets of closely lying four degenerate states which
allows scatterings from/to these states. Thus, here the correlation effects
are taken into account, while they were significantly reduced in the case of
$M=16$ and $N_{h}=9$, where the set of nine states were separated from other
single particles states by a large energy gap. While for T24, correlation
effects are more realistic than for symmetric meshes, at the same time one can
suspect that due to lack of the full symmetry of the lattice system, exchange
interactions are reduced. Thus, we consider T24 results as the other extreme
limit, in opposite to $N_{1}=N_{2}$ symmetric meshes, where factors
responsible for magnetization are too strongly reduced. Indeed, finite spin
polarization for the system with $N_{h}=7$ is quite unstable; for larger twist
angles it disappears. Notice that we consider here a limit of stronger
interactions, $\epsilon^{-1}=0.1$. We confirm high sensitivity of the
magnetization at this filling by looking at inverse susceptibility as a
function of temperature, what is shown in Fig. 9(c). Inverse susceptibility at
high temperature limit extrapolates to $T\approx 0$ for $\theta=2.5$ and
$\theta=3.0$, and to negative values for $\theta=3.5$; the magnetization at
finite temperature is not predicted from these calculations.

Figure 9: Tilted momentum space mesh $T24$ ED and FTLM results. (a) Energies
of finite size momentum space meshes for $M=24$ for $\theta=3.0$. Insets show
momentum space grid in the first Brillouin zone. (b) Total spin as a function
of filling factor $\nu$ for twist angles $\theta=2.5$, $\theta=3.0$, and
$\theta=3.5$. (c) Inverse susceptibility $\chi^{-1}(T)$ as a function of
temperature for the system with $N=7$ holes indicated by a dashed line in (b).
Parameters for moiré superlattice are: dielectric constant
$\epsilon^{-1}=0.1$, $V_{{\rm m}}=25$ meV. Parameters for FTLM: $N_{\rm
L}=150$, $N_{{\rm R}}=4$.

###  IV.2 Twist angle dependence

To complement our results for twist angle $\theta=3.0$, we analyze twist angle
dependence of conditions potentially leading to ferromagnetism. We investigate
the properties of the systems for $\theta=2.5$ and $\theta=3.5$. Fig. 10(a)
and (c) show kinetic energy difference between maximal and minimal spin
polarized states per moiré unit cell. For smaller (larger) twist angle, the
energetic cost for spin polarization is reduced (enlarged) because of decrease
(increase) of the bandwidth. Similarity between shape of lines of energetic
cost of spin polarization as a function of the filling factor in these two
figures, and the one in the main article, suggest that the shape of the band
does not change significantly when a twist angle is increased; it is mainly
nearest neighbor hopping $t$ that increases.

Figure 10: Kinetic energy difference between maximal and minimal total spin
states per moiré cell $(\Delta E_{\rm kin}=E^{kin}_{\rm SD}(S_{\rm
max})-E^{kin}_{\rm SD}(S_{\rm min}))/M$ as a function of band filling factor
$\nu=N/2M$ for different system sizes and for (a) $\theta=2.5$ and (c)
$\theta=3.5$. Inverse susceptibility $\chi^{-1}(T)$ as a function of
temperature for (b) $\theta=2.5$ and (d) $\theta=3.5$. Blue lines are linear
fits to estimate a transition temperature $T_{C}$ giving $T_{\rm C}\approx 4$
K for $\theta=2.5$ and $T_{\rm C}\approx 16$ K for $\theta=3.5$. The moiré
superlattice parameters used for this calculation: interaction strength
$\epsilon^{-1}=0.04$, modulation strength $V_{{\rm m}}=25$ meV and shape
$\psi=-94^{\circ}$. Parameters for FTLM: $N_{\rm L}=150$, $N_{{\rm
R}}=20,10,5$.

Changes in the bandwidth, translates into changes in transition temperatures
read from $\chi^{-1}$, shown in Fig. 10(b) and (d). However, one needs to be
caution with the band-mixing that might be significant for larger twist
angles, see the Discussion section. In Fig. 11 we show total spin as a
function of the filling factor for all three considered twist angles. In a
vicinity of van Hove singularity, indicated by a dashed line, finite spin
polarization is expected for all three systems. Increase of the twist angle to
$\theta>3.0$ leads to a transition between maximially spin polarized state at
half-filling, to minimal spin polarization; antiferromagnetic superexchange
dominates when the bandwidth increases (hopping integral increases) over
direct exchange interaction.

Figure 11: The ground state total spin $S$ as a function of filling factor
$\nu$ from exact diagonalization calculations for a system with $M=N_{\rm
1}\times N_{\rm 2}=16$ unit cells ($N_{\rm 1}$ and $N_{\rm 2}$. Comparison of
the spin polarization of moiré superlattice for different twist angles for
dielectric constant $\epsilon^{-1}=0.1$ and $V_{{\rm m}}=25$ meV,
$\psi=-94^{\circ}$. A dashed line indicates the position of the van Hove
singularity for finite size mesh.

###  IV.3 Hartree-Fock

We employ the self-consistent Hartree-Fock approximation to compute the
exchange energy for finite-size systems exactly. Hartree-Fock theory
approximates interaction through an additional self-energy term $\Sigma^{HF}$,
which can be expressed using a self-consistently determined single-particle
density matrix $\rho=\sum_{n}\ket{\psi_{n}}\bra{\psi_{n}}$. The sum here is
over the occupied moiré Bloch states. To evaluate the mean-field electronic
structure of moiré superlattices, a plane-wave representation is employed,
where the Hartree-Fock self-energy $\Sigma^{HF}$ at each $\mathbf{k}$ in the
Brillouin-zone is a matrix in moiré reciprocal lattice vectors $\mathbf{G}$:

$\displaystyle\Sigma^{HF}_{\alpha,\mathbf{G};\beta,\mathbf{G}^{\prime}}(\mathbf{k})=$
$\displaystyle\Sigma^{H}_{\alpha,\mathbf{G};\beta,\mathbf{G}^{\prime}}(\mathbf{k})+\Sigma^{F}_{\alpha,\mathbf{G};\beta,\mathbf{G}^{\prime}}(\mathbf{k})$
$\displaystyle=$
$\displaystyle\frac{\delta_{\alpha,\beta}}{A}\sum_{\alpha^{\prime}}V(\mathbf{G}^{\prime}-\mathbf{G})\sum_{\mathbf{k}^{\prime},\mathbf{G}^{\prime\prime}}\rho_{\alpha^{\prime},\mathbf{G}+\mathbf{G}^{\prime\prime};\alpha^{\prime},\mathbf{G}^{\prime}+\mathbf{G}^{\prime\prime}}(\mathbf{k}^{\prime})-\frac{1}{A}\sum_{\mathbf{G}^{\prime\prime},\mathbf{k}^{\prime}}V(\mathbf{G}^{\prime\prime}+\mathbf{k}^{\prime}-\mathbf{k})\rho_{\alpha,\mathbf{G}+\mathbf{G}^{\prime\prime};\beta,\mathbf{G}^{\prime}+\mathbf{G}^{\prime\prime}}(\mathbf{k}^{\prime}).$
(19)

In Eq. 19, Greek letters label spin, and
$\rho_{\alpha,\mathbf{G};\beta,\mathbf{G}^{\prime}}$ is the momentum-space
density matrix with a chosen symmetry. We do not project the interaction to
the lowest energy band in Hartree-Fock, all the bands are taken into account,
and the Coulomb interaction takes the standard form: $V(q)=2\pi e^{2}/\epsilon
q$. Starting with a physically plausible density matrix $\rho_{0}$, we
minimize the energy by performing self-consistent iterations. We consider a
state to be converged when the total Hartree-Fock energy difference between
iterations is less than $10^{-3}$ meV. The exchange energy is evaluated after
convergence by
$E_{\text{exch}}^{HF}=\operatorname{tr}(\Sigma^{F}(\mathbf{k})\rho(\mathbf{k}))/2$.
Exchange energy difference between maximally $S_{\rm max}$ and minimally
$S_{\rm min}$ spin polarized states normalized per moiré cell from Hartree-
Fock self-consistent calculations for are shown in Fig. 12(a)
for$\epsilon^{-1}=0.04$ and 12(b) for$\epsilon^{-1}=0.1$ interaction strength.
Clearly, mixing with other bands enhance this difference in comparison to
single Slater determinant calculations shown in Fig. 4(a) and 4(b).

To compare with ED calculations that are projected to the lowest-energy band,
the two relevant competing states for the present work are minimally spin-
polarized ($\nu_{\uparrow}=\nu_{\downarrow}=\nu/2$) and maximally spin-
polarized ferromagnetic ($\nu_{\uparrow}=\nu$ for $\nu\leq 1/2$;
$\nu_{\uparrow}=1/2,\nu_{\downarrow}=\nu-1/2$ for $\nu>1/2$) states. Our
calculations are carried out such that the particle number for each spin is
conserved. We use the same scheme of discretization for the Brillouin zone and
the same continuum model parameters as ED in Fig. 12.

Figure 12: Exchange energy difference between maximally $S_{\rm max}$ and
minimally $S_{\rm min}$ spin polarized states normalized per moiré cell from
Hartree-Fock self-consistent calculations for (a) $\epsilon^{-1}=0.04$ and (b)
$\epsilon^{-1}=0.1$ interaction strength. The results for twist angle
$\theta=3.0$, moiré modulation strength $V_{{\rm m}}=25$ meV, and potential
shape $\psi=-94^{\circ}$.

###  IV.4 Role of Long Range Interactions

Moire superlattice many-body model inherits many features from on-site Hubbard
model but in general long-range and non-local interactions affects the ground
state properties. This can be seen as ferromagnetism at half-filling for
sufficiently strong interactions, due to direct exchange, and from formation
of generalized Wigner crystal formation. These extra terms affects also charge
distributions of spin polarized and depolarized states because population of
single particles states is interaction dependent. In Fig. 13 we analyze how
long range direct interaction, together with direct exchange and assisted
hoppings, affects magnon spectrum at van Hove singularity for $M=16$ for
$\epsilon^{-1}=0.04$. The width of magnon spectrum from onsite and extended
Hubbard models are $\Delta E_{\rm Hubb}=2.8$ meV. This energy is increased by
around 1 meV, when non-local interaction terms are included. This also
translates to stronger tendency for ferromagnetism around van Hove singularity
for moiré systems in comparison to onsite Hubbard model, what can be read from
lower values of inverse susceptibility at low temperatures in Fig. 14(b).

Figure 13: $\theta=3.0$ analysis. Energy spectrum for total spin $S=S_{{\rm
max}}-1$ for the system with $N_{{\rm h}}=9$ holes, corresponding to filling
factor $\nu$ indicated by a dashed line in Fig. 3, $\epsilon^{-1}=0.04$, for
different models. $\Delta E$ indicates the width of magnon spectrum, which
determines transition temperature to a ferromagnetic phase. $N_{\rm
x}$($N_{\rm y}$) is a number of unit cells along two directions determined by
real space lattice vectors ${\bf a}_{1}$ and ${\bf a}_{2}$ on a triangular
lattice, $K_{\rm x}$($K_{\rm y}$) are total momenta along two direction
determined by reciprocal space lattice vectors ${\bf b}_{1}$ and ${\bf
b}_{2}$. Parameters for moiré superlattice are: dielectric constant
$\epsilon^{-1}=0.04$, twist angle $\theta=3.0$, $V_{{\rm m}}=25$ meV.  Figure
14: Comparison between moiré system and Hubbard model results. Inverse
susceptibility $\chi^{-1}(T)$ as a function of temperature. Calculations for
$M=16$ with (a) $N=16$ electrons $\nu=0.5$ and (b) $N=23$ particles $\nu=0.73$
($N_{{\rm h}}=9$ holes). The number of random vectors $N_{{\rm R}}=4$.
Parameters for moiré superlattice are: dielectric constant
$\epsilon^{-1}=0.04$, twist angle $\theta=3.0$, $V_{{\rm m}}=25$ meV.