# Stoner ferromagnetism in a momentum-confined interacting 2D electron gas

Ohad Antebi  Department of Condensed Matter Physics, Weizmann Institute of
Science, Rehovot 76100, Israel  Ady Stern  Department of Condensed Matter
Physics, Weizmann Institute of Science, Rehovot 76100, Israel  Erez Berg
Department of Condensed Matter Physics, Weizmann Institute of Science, Rehovot
76100, Israel

###### Abstract

In this work we investigate the ground state of a momentum-confined
interacting 2D electron gas, a momentum-space analog of an infinite quantum
well. The study is performed by combining analytical results with a numerical
exact diagonalization procedure. We find a ferromagnetic ground state near a
particular electron density and for a range of effective electron (or hole)
masses. We argue that this observation may be relevant to the generalized
Stoner ferromagnetism recently observed in multilayer graphene systems. The
collective magnon excitations exhibit a linear dispersion, which originates
from a diverging spin stiffness.

_Introduction.–_ The ground state of a 2D electron gas is determined by a
competition between kinetic and interaction energies. In the absence of
external fields or ionic potentials, numerical calculations[1, 2, 3, 4, 5]
suggest the ground state is either a paramagnetic Fermi liquid or a Wigner
crystal[6], depending on the density, while a Stoner ferromagnetic[7] Fermi
liquid state is a close competitor. By applying an out-of-plane magnetic
field, the energy spectrum of the electron gas forms Landau levels. In the
absence of Zeeman coupling the ground state is known to be ferromagnetic for
densities near one electron per flux quantum[8, 9].

Recent experiments in Bernal bilayer and rhombohedrally stacked multilayer
graphene show spin and valley ferromagnetic phases for certain regimes of
electron density and out-of-plane displacement field[10, 11, 12, 13, 14, 15].
The band dispersion of such systems is approximately flat up to some momentum
scale, at which point the kinetic energy increases rapidly with momentum. This
corresponds to a high density of states in a bounded region of momentum space,
and is an example of a partially flat band. Here, we argue that such
dispersion is favorable for Stoner ferromagnetism when the flat region is
nearly fully occupied. We summarize these results in Fig. 1. Examples of
various correlated phenomena that had been previously studied on models with
partially flat bands can be found in Refs.[16, 17, 18, 19, 20, 21].

Motivated by multilayer rhombohedrally stacked graphene, we present a toy
model for an interacting 2D electron gas in which the kinetic energy of the
electrons diverges beyond a limited region in momentum space. We will refer to
this as a momentum-confined gas. We find the ground state of this gas to be
spin-polarized. We find a linearly dispersing magnon branch of excitations, in
contrast with the commonly predicted quadratic dispersion[22]. We trace this
apparent anomaly to a divergence of the spin stiffness associated with the
infinite slope of the kinetic energy dispersion. We show that the spin
polarization is robust to an addition of weak dispersion to the kinetic energy
within the allowed momentum region.

|  Low-density semiconductor |  Rhombohedral multilayer graphene |  Landau levels
---|---|---|---
Single-particle dispersion |  |  |
Stoner Magnetism |  No* |  Yes** |  Yes
Figure 1: Comparison of dispersion and Stoner ferromagnetism between different
2D systems. The dashed horizontal line represents the chemical potential.
*Based on state of the art quantum Monte Carlo studies [5]. **Experimentally
known for bilayer[11, 12, 13, 14] and trilayer[10] graphene, both with a
strong perpendicular displacement field, and analytically shown in this work
for infinitely many layers.

_Hamiltonian.–_ Our goal is to construct a simple model that captures non-
trivial interaction effects of electrons confined in momentum space, with a
minimal set of parameters. We consider a 2D electron gas with an artificially
constructed kinetic dispersion, and turn on the conventional Coulomb
repulsion. We imagine the kinetic dispersion to be such that the electrons are
only allowed to occupy momenta in a disk of radius $k_{0}$ in momentum space,
whose area is small compared to the Brillouin zone. Idealizing the dispersion,
we set the kinetic energy to zero within the disk and to infinity outside it:

$E_{\text{K}}(\mathbf{k})=\begin{cases}0,&|\mathbf{k}|\leq k_{0},\\\
\infty,&|\mathbf{k}|>k_{0}.\end{cases}$ (1)

This dispersion leads to the strict confinement of the electrons to a disk in
momentum space, forming a momentum space analog of a 2D circular infinite
quantum well in real space.

The confinement in Eq. (1) can be thought of as mimicking the low energy band
dispersion of $N_{l}$ layers of rhombohedral graphite subject to a
displacement field $D$. For this system, there is a band gap proportional to
$D$, the dispersion is very flat up to a momentum scale $k_{0}$, and rises as
$E(k)\sim(k/k_{0})^{N_{l}}$ for $k>k_{0}$. In the limit of
$N_{l}\rightarrow\infty$, the scale $k_{0}$ is set by the ratio of the
interlayer tunneling $t_{\perp}$ and the monolayer Dirac velocity $v_{D}$ (see
S1 in the Supplemental Material - SM). In this limit, the dispersion resembles
the idealized momentum confined model of Eq. (1).

Next, we examine the effects of electron-electron interactions in such a
momentum-confined setup. We consider our Hamiltonian to be the normal-ordered
2D screened Coulomb interaction strictly confined in momentum space:

$\displaystyle\mathcal{H}_{\text{int}}$
$\displaystyle=\frac{1}{2A}\sum_{\mathbf{q}}V_{\mathbf{q}}:\overline{\rho_{\mathbf{q}}}^{\dagger}\overline{\rho_{\mathbf{q}}}:,$
(2) $\displaystyle\overline{\rho_{\mathbf{q}}}$
$\displaystyle=\sum_{\sigma=\uparrow,\downarrow}\sum_{\mathbf{k}\in
D_{\mathbf{q}}}c_{\sigma,\mathbf{k}+\mathbf{q}}^{\dagger}c_{\sigma,\mathbf{k}}.$
(3)

where $c_{\sigma,\mathbf{k}}$ annihilates an electron with spin
$\sigma\in\\{\uparrow,\downarrow\\}$ and momentum $\mathbf{k}$, $A$ is the
system area, $V_{\mathbf{q}}=2\pi e^{2}\tanh{(|\mathbf{q}|d)}/|\mathbf{q}|$ is
the Fourier transform of the Coulomb potential, and $d$ is the distance to the
screening gate. The sum over momenta in the confined particle density
operator, $\overline{\rho_{\mathbf{q}}}$, is limited to the domain
$D_{\mathbf{q}}=\\{\mathbf{k}\enspace|\enspace(|\mathbf{k}|\leq
k_{0})\cap(|\mathbf{k}+\mathbf{q}|\leq k_{0})\\}$ such that the fermion
operators are within the disk of radius $k_{0}$ (see Fig. 2). While this
Hamiltonian has only two (spin) flavors, some of the conclusions we draw below
are applicable also in the presence of multiple valleys. We remark that any
physical realization of the kinetic dispersion would be accompanied by form
factors in Eq. (3) that originate from the Bloch wavefunctions. For
simplicity, in this work we set the form factors to one.

Figure 2: The kinetic dispersion confines electrons to the disk $D_{0}$ of
radius $k_{0}$ in momentum space. The region $D_{\mathbf{q}}$ is the
intersection of two disks mutually shifted by $\mathbf{q}$. The strict
momentum confinement forbids electrons with momenta in $D_{0}-D_{\mathbf{q}}$
from scattering with a momentum transfer $\mathbf{q}$.

The Hamiltonian in Eq. (2) is the one we will explore for the rest of the
paper. It has a continuous rotational symmetry around $\mathbf{k}=0$, and a
global $SU(2)$ symmetry for spin rotations. The most interesting property
emerges in the limit of unscreened Coulomb interaction,
$\left.k_{0}d\rightarrow\infty\right.$. In this limit, the Hamiltonian has
only one length scale, given by $l\sim k_{0}^{-1}$. Consequently, we find a
single energy scale $E\sim e^{2}k_{0}$, and an electron density scale
$n_{0}=k_{0}^{2}/(4\pi)$, which corresponds to completely filling the disk
with a single spin flavor. Defining the filling factor by $\nu=n/n_{0}$, we
will focus our discussion in this paper on $\nu=1$ and on
$\nu=1\pm\varepsilon$ for $0<\varepsilon\ll 1$.

_The role of normal-ordering.–_ The notion of the single-particle kinetic
dispersion in the presence of many-body interactions is not without
subtleties. In our model, the flat nature of the dispersion is tied to our
choice to consider the normal-ordered form of the interaction operator.
Physically, the normal-ordering prevents an electron from interacting with
itself. By undoing the normal-ordering of the confined Coulomb interaction, we
find that Eq. (2) can be separated into a positive semi-definite density-
density operator and a single-particle operator:

$\mathcal{H}_{\text{int}}=\underbrace{\frac{1}{2A}\sum_{\mathbf{q}}V_{\mathbf{q}}\overline{\rho_{\mathbf{q}}}^{\dagger}\overline{\rho_{\mathbf{q}}}}_{=\mathcal{H}_{\rho-\rho}}-\sum_{\sigma=\uparrow,\downarrow}\sum_{|\mathbf{k}|\leq
k_{0}}E_{0}(\mathbf{k})c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}},$
(4)

where we introduced the dispersion:

$E_{0}(\mathbf{k})=\frac{1}{2A}\sum_{|\mathbf{k}^{\prime}|\leq
k_{0}}V_{\mathbf{k}-\mathbf{k}^{\prime}}.$ (5)

This dispersion is, up to a sign, the exchange contribution to the self-energy
of a particle at momentum $\mathbf{k}$ in the presence of a completely filled
disk of radius $k_{0}$. In contrast with the standard unconfined interaction,
we find that the momentum cutoff introduces a non-trivial dispersion term
associated with the normal-ordering.

Figure 3: The dispersion $E_{0}(|\mathbf{k}|)$ for different screening lengths
of the Coulomb interaction. In the case of unscreened Coulomb interaction
($k_{0}d\rightarrow\infty$) the dispersion is singular at $k=k_{0}$. This
singularity is smoothed out for screened Coulomb interaction.

The dispersion in Eq. (5) plays an important role in several properties of our
model and thus warrants an explicit discussion. In the limit of unscreened
interaction, $k_{0}d\rightarrow\infty$, this dispersion is given by:

$E_{0}(\mathbf{k})=e^{2}k_{0}\frac{1}{\pi}\tilde{E}\left(\frac{\pi}{2},\frac{|\mathbf{k}|}{k_{0}}\right)$
(6)

where $\tilde{E}\left(\frac{\pi}{2},x\right)$ is the complete elliptic
integral of the second kind whose argument $x$ is the elliptic modulus. The
derivative of Eq. (6) diverges logarithmically at $|\mathbf{k}|=k_{0}$. This
divergence is cut off for finite $k_{0}d$. Expanding Eq. (5) in
$(k_{0}d)^{-1}$, we find that the series does not uniformly converge on the
disk, but nevertheless, the correction is small everywhere for $k_{0}d\gg 1$.
The effect of screening length on the dispersion is plotted in Fig. 3. The
details are described in the SM S2.

_Exact ground state.–_ The ground state of Eq. (2) is not analytically
solvable. We can, however, solve for the exact ground state upon introducing a
particular dispersion to the otherwise flat disk, which would eliminate the
contribution of Eq. (5). With this dispersion, we are left with the density-
density Hamiltonian $\mathcal{H}_{\rho-\rho}$ as defined in Eq. (4). The
Hamiltonian $\mathcal{H}_{\rho-\rho}$ is positive semi-definite, which
provides a direct route to finding its exact ground state. For $\nu=1$,
consider completely filling the spin $\sigma=\downarrow$ flavor, i.e. a
maximally spin-polarized state. We denote this state by
$\ket{\Psi_{\text{SP}}}$. By applying Eq. (3) one easily finds that:

$\forall\mathbf{q}:\overline{\rho_{\mathbf{q}}}\ket{\Psi_{\text{SP}}}=0\implies\mathcal{H}_{\rho-\rho}\ket{\Psi_{\text{SP}}}=0,$
(7)

which proves that $\ket{\Psi_{\text{SP}}}$ is a ground state of
$\mathcal{H}_{\rho-\rho}$. We emphasize that this ground state is degenerate
with $SU(2)$ spin rotation symmetry. For a model with multiple valleys, the
above statement holds for any integer filling $\nu$ and any completely spin
and valley polarized state, provided that a large momentum separation between
different valleys allows for a neglect of inter-valley scattering in Coulomb
processes.

The spin polarized state $\ket{\Psi_{\text{SP}}}$ is also an eigenstate of
$\mathcal{H}_{\text{int}}$. We hypothesize that $\ket{\Psi_{\text{SP}}}$ is a
ground state of $\mathcal{H}_{\text{int}}$. We provide numerical evidence for
this claim below. Assuming this hypothesis is correct, we proceed to calculate
the excitations with respect to this state.

_Single-particle excitations.–_ The single electron and hole excitations
energies relative to the fully spin $\sigma=\downarrow$ polarized state,
$E_{e}(\mathbf{k})$ and $E_{h}(\mathbf{k})$ respectively, are defined by:

$\displaystyle\begin{split}[\mathcal{H}_{\text{int}},c_{\uparrow,\mathbf{k}}^{\dagger}]\ket{\Psi_{\text{SP}}}&=E_{e}(\mathbf{k})c_{\uparrow,\mathbf{k}}^{\dagger}\ket{\Psi_{\text{SP}}},\\\
[\mathcal{H}_{\text{int}},c_{\downarrow,\mathbf{k}}]\ket{\Psi_{\text{SP}}}&=E_{h}(\mathbf{k})c_{\downarrow,\mathbf{k}}\ket{\Psi_{\text{SP}}},\end{split}$
(8)

and by direct calculation are found to be:

$\displaystyle\begin{split}E_{e}(\mathbf{k})&=\frac{1}{A}\sum_{|\mathbf{k}^{\prime}|\leq
k_{0}}V_{\mathbf{0}},\\\
E_{h}(\mathbf{k})&=\frac{1}{A}\sum_{|\mathbf{k}^{\prime}|\leq
k_{0}}\left(V_{\mathbf{k}-\mathbf{k}^{\prime}}-V_{\mathbf{0}}\right).\end{split}$
(9)

The identity
$E_{0}(\mathbf{k})=\left(E_{e}(\mathbf{k})+E_{h}(\mathbf{k})\right)/2$ holds
also if one adds a spin-independent single-particle dispersion to the
Hamiltonian $\mathcal{H}_{\text{int}}$. The electron excitation energy is the
charging energy of the system’s geometric capacitance, and thus has a flat
dispersion. The hole excitation energy has the same contribution (with an
opposite sign), along with the exchange interaction of the missing electron.
The dispersion of the hole excitation energy is therefore shown in Fig. 3 up
to a factor of two. We find that the hole dispersion is such that the lowest
energy is obtained by removing a hole from the edge of the disk, thereby
reducing the Fermi sea radius by an infinitesimal amount. This is identical to
the Fermi liquid behavior.

_Collective excitations.–_ The collective particle-hole excitations of
momentum $\mathbf{Q}$ and spin $\hbar$ with respect to
$\ket{\Psi_{\text{SP}}}$ are eigenstates of $\mathcal{H}_{\text{int}}$ spanned
by wavefunctions of the form:

$\ket{\Psi_{\text{ph}}^{\mathbf{k},\mathbf{Q}}}=c_{\uparrow,\mathbf{k}+\mathbf{Q}}^{\dagger}c_{\downarrow,\mathbf{k}}\ket{\Psi_{\text{SP}}}.$
(10)

The solution is given by diagonalizing the following matrix:

$[H(\mathbf{Q})]_{\mathbf{k},\mathbf{k}^{\prime}}=\bra{\Psi_{\text{ph}}^{\mathbf{k}^{\prime},\mathbf{Q}}}\mathcal{H}_{\text{int}}\ket{\Psi_{\text{ph}}^{\mathbf{k},\mathbf{Q}}},\quad\mathbf{k},\mathbf{k}^{\prime}\in
D_{\mathbf{Q}},$ (11)

This restriction of both indices to the domain
$D_{\mathbf{Q}}=\\{\mathbf{k}\enspace|\enspace(|\mathbf{k}|\leq
k_{0})\cap(|\mathbf{k}+\mathbf{Q}|\leq k_{0})\\}$ is due to our infinite
kinetic energy dispersion for states with momentum outside the disk. This
sharp cutoff leads to a striking result - the lowest-lying particle-hole
excitation is massless. For small $|\mathbf{Q}|/k_{0}$ we find:

$E_{ph}(\mathbf{Q})\approx\frac{4E_{0}(k_{0})}{\pi k_{0}}|\mathbf{Q}|.$ (12)

This is surprising, as for a magnon we normally expect to find a finite spin
stiffness $\rho_{s}$ such that
$E_{ph}(\mathbf{Q})\approx\rho_{s}|\mathbf{Q}|^{2}$. The above result holds in
the presence of any rotationally symmetric single-particle dispersion, see SM
S3.

To resolve this apparent discrepancy, we replace our strictly confining
kinetic dispersion with the following:

$\mathcal{H}_{\text{K}}=U_{\text{K}}\sum_{\sigma=\uparrow,\downarrow}\sum_{\mathbf{k}}\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{N_{D}}c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}},$
(13)

where $U_{\text{K}}$ is chosen such that $0<U_{\text{K}}\ll e^{2}k_{0}$. Under
this softened momentum confinement the momentum summations are unconstrained
since particles are allowed to be excited beyond the disk. We emphasize that
in the limit of $N_{D}\rightarrow\infty$ our strictly confining kinetic
dispersion is restored. For $N_{D}\gg 1$ we find, using second order
perturbation theory, that the spin stiffness is finite, and the magnon
dispersion is quadratic:

$E_{ph}(\mathbf{Q})\approx\frac{1}{2}\frac{U_{\text{K}}}{k_{0}^{2}}N_{D}|\mathbf{Q}|^{2}.$
(14)

However, in the limit $N_{D}\rightarrow\infty$, the spin stiffness
$\left.\rho_{s}\sim N_{D}\right.$ diverges. This leads to a breakdown of
perturbation theory as $|\mathbf{Q}|/k_{0}$ can no longer be used as a small
parameter. The details are described in SM S3.

_Phase diagram near $\nu=1$.–_ We have seen above that when removing electrons
(or adding holes) to the spin-polarized ground state $\ket{\Psi_{\text{SP}}}$,
it is energetically favorable to start with the outermost states, thereby
reducing the radius of the occupied disk. Consider removing a small fraction
of the electrons from the system, i.e. setting the system at filling fraction
$\nu=1-\varepsilon$ for $0<\varepsilon\ll 1$. We now argue that the system has
a spin-polarized Fermi liquid ground state for small enough $\varepsilon>0$.
In the limit of $\varepsilon\rightarrow 0$, the spin polarization is
undiminished, and the system can be thought of as having a single species of
particles (in this case, holes), with a single particle dispersion
$E_{h}(\mathbf{k})=2E_{0}(\mathbf{k})$ and subject to the 2D screened Coulomb
interaction. Partially filling this system leads to a circular Fermi surface.
The finite interaction strength (due to screening) and the constraint to 2D
imply a Fermi liquid ground state. Presumably, continuing to decrease the
filling will eventually destroy the spin polarization, and we expect a phase
transition to a spin-depolarized state.

When the dispersion is exactly flat ($N_{D}\rightarrow\infty$) there is
difficulty making a similar argument for filling $\nu=1+\varepsilon$, as the
single-electron excitation spectrum is completely flat. However, as previously
discussed, the average of single-electron and single-hole excitation spectra
is constrained to be $E_{0}(\mathbf{k})$. Therefore, upon introduction of some
weak dispersion to the disk, for example - by setting $N_{D}$ to a finite
value, we expect both electron and hole sides to display a Fermi liquid phase
by the exact same reasoning. Generally, we note that the system has no
particle-hole symmetry with respect to $\nu=1$.

_Numerical Analysis.–_ In order to support our hypothesis of a spin-polarized
ground state, we have performed an exact diagonalization study of Eq. (2) for
finite systems with periodic boundary conditions at filling factor $\nu=1$ for
$N=7$ and $N=19$ particles. The diagonalization is done using the matrix-free
implicitly restarted Lanczos method. We find the ground state at $\nu=1$ to be
the fully spin-polarized state for both system sizes considered. Further
technical details are given in the SM S4.

Figure 4: The ground state spin polarization of the Hamiltonian in Eq. (15) at
$\nu=1$ as a function of inverse effective mass $m_{\text{eff}}^{-1}$,
computed for two system sizes of $N=7$ and $N=19$ particles. The horizontal
axis increments are $m_{0}^{-1}=V_{k_{0}}N/(2Ak_{0}^{2})$. The unpolarized
phase has a sub-extensive polarization explained by Hund’s rule.

To further solidify our claim, we have examined the robustness of the spin-
polarized ground state to additional dispersion. We repeat the exact
diagonalization at $\nu=1$ for the confined Coulomb Hamiltonian
$\mathcal{H}_{\text{int}}$ in Eq. (2) with an additional quadratic dispersion
($\hbar=1$):

$\mathcal{H}_{\text{int+mass}}=\mathcal{H}_{\text{int}}+\sum_{\sigma=\uparrow,\downarrow}\sum_{|\mathbf{k}|\leq
k_{0}}\frac{|\mathbf{k}|^{2}}{2m_{\text{eff}}}c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}}.$
(15)

In Fig. 4 we plot the total spin $S$ of the ground state of Eq. (15) at
$\nu=1$ vs $m_{\text{eff}}^{-1}$, as the latter is swept from negative to
positive values in increments of $m_{0}^{-1}=V_{k_{0}}N/(2Ak_{0}^{2})$. The
flat kinetic dispersion in Eq. (1) corresponds to $m_{\text{eff}}^{-1}=0$. For
the range of finite systems we considered, we find that the ground state is
fully spin-polarized, i.e. $S=N/2$, for $-3m_{0}^{-1}\leq
m_{\text{eff}}^{-1}\leq 2m_{0}^{-1}$. Outside this region we find the spin
polarization drops to $S=5/2$, precisely the value expected by Hund’s rule
applied to partially filling the degenerate finite-size Fermi surface (see SM
S4). We associate this result with a spin unpolarized phase, and expect that
in the thermodynamic limit we will obtain $S/N\rightarrow 0$ in this regime of
parameters. The asymmetry of the spin-polarized regime boundary with respect
to the sign of $m_{\text{eff}}^{-1}$ is due to the lower exchange energy of
the unpolarized disk-shaped Fermi sea compared to that of an annular-shaped
Fermi sea, the two different shapes that correspond to positive and negative
$m_{\text{eff}}^{-1}$, respectively. In multilayer graphene systems with
typical values of $k_{0}=0.2\textrm{nm}^{-1}$ and dielectric coefficient
$\epsilon_{d}=6$, the mass scale is given by
$m_{0}\sim\epsilon_{d}k_{0}/e^{2}\sim 0.06m_{e}$ where $m_{e}$ is the free
electron mass.

_Conclusions.–_ We considered a model of interacting electrons confined in
momentum space, and showed that when nearly half of the available states are
filled, its ground state is spin polarized. We then numerically demonstrated
that this spin polarization is stable to small variations of the model, and is
not a fine-tuned consequence of the kinetic dispersion. We find this ground
state to have a non-trivial excitation spectrum, including a diverging spin
stiffness which leads to massless particle-hole excitations at low momenta. We
argued in favor of a ferromagnetic Fermi liquid phase at
$\nu=1\pm\varepsilon$, thereby demonstrating Stoner ferromagnetism in this
model. The prevalence of spin and valley Stoner ferromagnets in multilayer
rhombohedral graphene in a displacement field, where for a range of densities
electrons are confined to a flat region of momentum space, may be understood
within the framework of our model.

###### Acknowledgements.

OA would like to thank Johannes S. Hofmann, Daniel Kaplan, Tobias Holder and
Shahar Barkai for the useful discussions. We acknowledge support from the
Israeli Science Foundation Quantum Science and Technology grant no. 2074/19,
and from the CRC 183 of the Deutsche Forschungsgemeinschaft (Project C02).
This work has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme [grant
agreements No. 788715 and 817799, Projects LEGOTOP (AS) and HQMAT (EB)].

## References

  * Ceperley and Alder [1980] D. M. Ceperley and B. J. Alder, Ground state of the electron gas by a stochastic method, Phys. Rev. Lett. 45, 566 (1980).
  * Tanatar and Ceperley [1989] B. Tanatar and D. M. Ceperley, Ground state of the two-dimensional electron gas, Phys. Rev. B Condens. Matter 39, 5005 (1989).
  * Rapisarda and Senatore [1996] F. Rapisarda and G. Senatore, Diffusion Monte Carlo study of electrons in two-dimensional layers, Aust. J. Phys. 49, 161 (1996).
  * Attaccalite  _et al._ [2002] C. Attaccalite, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet, Correlation energy and spin polarization in the 2D electron gas, Phys. Rev. Lett. 88, 256601 (2002).
  * Drummond and Needs [2009] N. D. Drummond and R. J. Needs, Phase diagram of the low-density two-dimensional homogeneous electron gas, Phys. Rev. Lett. 102, 126402 (2009).
  * Wigner [1934] E. Wigner, On the interaction of electrons in metals, Phys. Rev. 46, 1002 (1934).
  * Stoner [1938] E. C. Stoner, Collective electron ferronmagnetism, Proc. R. Soc. Lond. 165, 372 (1938).
  * Nomura and MacDonald [2006] K. Nomura and A. H. MacDonald, Quantum Hall ferromagnetism in graphene, Phys. Rev. Lett. 96, 256602 (2006).
  * Girvin [1999] S. M. Girvin, The Quantum Hall Effect: Novel excitations and broken symmetries, in  _Aspects topologiques de la physique en basse dimension. Topological aspects of low dimensional systems_ (Springer Berlin Heidelberg, 1999) pp. 53–175.
  * Zhou  _et al._ [2021] 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 (2021).
  * 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).
  * Zhang  _et al._ [2023] Y. Zhang, R. Polski, A. Thomson, É. Lantagne-Hurtubise, C. Lewandowski, H. Zhou, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Enhanced superconductivity in spin–orbit proximitized bilayer graphene, Nature 613, 268 (2023).
  * de la Barrera  _et al._ [2022] S. C. de la Barrera, S. Aronson, Z. Zheng, K. Watanabe, T. Taniguchi, Q. Ma, P. Jarillo-Herrero, and R. Ashoori, Cascade of isospin phase transitions in Bernal-stacked bilayer graphene at zero magnetic field, Nature Physics 18, 771 (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).
  * Shi  _et al._ [2020] Y. Shi, S. Xu, Y. Yang, S. Slizovskiy, S. V. Morozov, S.-K. Son, S. Ozdemir, C. Mullan, J. Barrier, J. Yin, _et al._ , Electronic phase separation in multilayer rhombohedral graphite, Nature 584, 210 (2020).
  * Mielke [1999] A. Mielke, Ferromagnetism in single-band Hubbard models with a partially flat band, Phys. Rev. Lett. 82, 4312 (1999).
  * Tang and Fu [2014] E. Tang and L. Fu, Strain-induced partially flat band, helical snake states and interface superconductivity in topological crystalline insulators, Nat. Phys. 10, 964 (2014).
  * Huang  _et al._ [2019] E. W. Huang, M.-S. Vaezi, Z. Nussinov, and A. Vaezi, Enhanced correlations and superconductivity in weakly interacting partially flat-band systems: A determinantal quantum Monte Carlo study, Phys. Rev. B Condens. Matter 99, 235128 (2019).
  * Sayyad  _et al._ [2020] S. Sayyad, E. W. Huang, M. Kitatani, M.-S. Vaezi, Z. Nussinov, A. Vaezi, and H. Aoki, Pairing and non-Fermi liquid behavior in partially flat-band systems: Beyond nesting physics, Phys. Rev. B Condens. Matter 101, 014501 (2020).
  * Aoki [2020] H. Aoki, Theoretical possibilities for flat band superconductivity, J. Supercond. Novel Magn. 33, 2341 (2020).
  * Sayyad  _et al._ [2023] S. Sayyad, M. Kitatani, A. Vaezi, and H. Aoki, Nematicity-enhanced superconductivity in systems with a non-Fermi liquid behavior, Journal of Physics: Condensed Matter 35, 245605 (2023).
  * Watanabe [2020] H. Watanabe, Counting rules of Nambu–Goldstone modes, Annu. Rev. Condens. Matter Phys. 11, 169 (2020).
  * Min and MacDonald [2008] H. Min and A. H. MacDonald, Chiral decomposition in the electronic structure of graphene multilayers, Phys. Rev. B 77, 1 (2008).
  * Zhang  _et al._ [2010] F. Zhang, B. Sahu, H. Min, and A. H. MacDonald, Band structure of ABC-stacked graphene trilayers, Phys. Rev. B 82, 1 (2010).
  * Gershgorin [1931] S. A. Gershgorin, Uber die Abgrenzung der Eigenwerte einer Matrix, lzv. Akad. Nauk. USSR. Otd. Fiz-Mat. Nauk 7, 749 (1931).
  * Varga [2011] R. S. Varga,  _Geršgorin and His Circles_ (Springer Science & Business Media, 2011).
  * Lehoucq  _et al._ [1998] R. B. Lehoucq, D. C. Sorensen, and C. Yang,  _ARPACK users’ guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods_ (Society for Industrial and Applied Mathematics, 1998) pp. i–xv.
  * [28] opencollab, ARPACK-NG: A new generation of implicitly restarted Arnoldi library, https://github.com/opencollab/arpack-ng, accessed: NA-NA-NA.

Supplemental Material

###### Contents

  1. S1 Relation to rhombohedral graphite
  2. S2 The role of normal-ordering
  3. S3 Collective particle-hole excitations
    1. S3.1 The strictly confining kinetic dispersion
      1. S3.1.1 Constructing the collective excitation matrix
      2. S3.1.2 Excitation spectrum bounds
      3. S3.1.3 Dispersion of the lowest excitation energy
    2. S3.2 The softened momentum confinement
  4. S4 Numerical Analysis
    1. S4.1 Discretization of momentum space
    2. S4.2 Exact diagonalization
      1. S4.2.1 Procedure to obtain the ground state’s energy and spin polarization
      2. S4.2.2 Simulation results

##  S1 Relation to rhombohedral graphite

In this section we show how the confining kinetic dispersion in Eq. (1) can be
thought of as a limit of infinitely many layers of rhombohedral graphite in
the presence of a perpendicular electric field. Consider a graphene multilayer
whose sublattices at layer $n$ are denoted by $A_{n}$ and $B_{n}$. We say that
the graphene is stacked rhombohedrally, or chirally[23], if the sublattice
$B_{n+1}$ is directly above the sublattice $A_{n}$ for every $n$. For two
layers this stacking configuration is known as Bernal stacking, while for
three layers it is sometimes referred to as ABC stacking[24]. In the
derivation below we follow Ref. [23], and assume the only interlayer tunneling
is between $A_{n}$ and $B_{n+1}$. For $N_{l}$ layers, the effective low energy
theory is hosted on sublattices $A_{1}$ and $B_{N_{l}}$, and is given (near
one valley) by:

$H_{\text{eff}}(\mathbf{k})=-t_{\perp}\begin{pmatrix}0&\left(\frac{v_{D}}{t_{\perp}}(k_{x}-ik_{y})\right)^{N_{l}}\\\
\left(\frac{v_{D}}{t_{\perp}}(k_{x}+ik_{y})\right)^{N_{l}}&0\end{pmatrix}.$
(S1)

Here, $v_{D}$ is the monolayer Dirac velocity, and $t_{\perp}$ is the
interlayer tunneling amplitude. This Hamiltonian displays a natural momentum
scale given by $k_{0}=t_{\perp}/v_{D}$ which will simplify the notation below.

Next, consider applying a constant out-of-plane displacement field $D$. The
two sites $A_{1}$ and $B_{N_{l}}$ are vertically separated by
$(N_{l}-1)a_{\perp}$, where $a_{\perp}$ is the interlayer distance. The
effective Hamiltonian in the presence of this displacement field can be
written as:

$H_{\text{eff}}(\mathbf{k})=\begin{pmatrix}\frac{eDa_{\perp}(N_{l}-1)}{2}&-t_{\perp}\left(\frac{k_{x}-ik_{y}}{k_{0}}\right)^{N_{l}}\\\
-t_{\perp}\left(\frac{k_{x}+ik_{y}}{k_{0}}\right)^{N_{l}}&-\frac{eDa_{\perp}(N_{l}-1)}{2}\end{pmatrix},$
(S2)

which yields two bands, separated in energy by a gap that is roughly
proportional to the displacement field and the number of layers:

$\epsilon_{\pm}(\mathbf{k})=\pm\sqrt{\left(\frac{eDa_{\perp}(N_{l}-1)}{2}\right)^{2}+t_{\perp}^{2}\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{2N_{l}}}.$
(S3)

Considering many layers of graphene, i.e. $N_{l}\gg 1$, we find the following
approximate energy dispersion:

$\epsilon_{\pm}(\mathbf{k})\approx\begin{cases}\pm\frac{e\left|D\right|a_{\perp}N_{l}}{2}+\mathcal{O}\left(\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{2N_{l}}\right),&|\mathbf{k}|<k_{0}\\\
\pm\left|t_{\perp}\right|\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{N_{l}}+\mathcal{O}\left(N_{l}^{2}\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{-N_{l}}\right),&|\mathbf{k}|>k_{0}\end{cases}.$
(S4)

In the limit of $N_{l}\rightarrow\infty$ this kinetic dispersion has a simple
form,

$\epsilon_{+}(\mathbf{k})-\epsilon_{+}(0)=\begin{cases}0,&|\mathbf{k}|<k_{0}\\\
\infty,&|\mathbf{k}|>k_{0}\end{cases},$ (S5)

which is the confining kinetic dispersion in Eq. (1).

We note that, in practice, it is difficult to open a gap in rhombohedral
graphite for too large $N_{l}$, since the perpendicular applied field is
screened by the bulk electrons.

##  S2 The role of normal-ordering

In this section we derive the dispersion in Eq. (5) and calculate its explicit
form. The momentum-confined Hamiltonian, as given by Eq. (2), is rewritten
here as:

$\displaystyle\begin{split}\mathcal{H}_{\text{int}}=\frac{1}{2A}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}\sum_{\mathbf{q}}\sum_{\mathbf{k}\in
D_{\mathbf{q}}}\sum_{\mathbf{k}^{\prime}\in
D_{-\mathbf{q}}}V_{\mathbf{q}}c_{\sigma^{\prime},\mathbf{k}^{\prime}-\mathbf{q}}^{\dagger}c_{\sigma,\mathbf{k}+\mathbf{q}}^{\dagger}c_{\sigma,\mathbf{k}}c_{\sigma^{\prime},\mathbf{k}^{\prime}},\quad
V_{\mathbf{q}}=2\pi
e^{2}\frac{\tanh{(|\mathbf{q}|d)}}{|\mathbf{q}|},\end{split}$ (S6)

and we reiterate our definition of
$D_{\mathbf{q}}=\\{\mathbf{k}\enspace|\enspace\left(|\mathbf{k}|\leq
k_{0}\right)\cap\left(|\mathbf{k}+\mathbf{q}|\leq k_{0}\right)\\}$, shown in
Fiq. 2. By undoing the normal ordering, we find that:

$\displaystyle\begin{split}\mathcal{H}_{\text{int}}&=\frac{1}{2A}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}\sum_{\mathbf{q}}\sum_{\mathbf{k}\in
D_{\mathbf{q}}}\sum_{\mathbf{k}^{\prime}\in
D_{-\mathbf{q}}}V_{\mathbf{q}}c_{\sigma^{\prime},\mathbf{k}^{\prime}-\mathbf{q}}^{\dagger}c_{\sigma^{\prime},\mathbf{k}^{\prime}}c_{\sigma,\mathbf{k}+\mathbf{q}}^{\dagger}c_{\sigma,\mathbf{k}}-\frac{1}{2A}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}\sum_{\mathbf{q}}\sum_{\mathbf{k}\in
D_{\mathbf{q}}}\sum_{\mathbf{k}^{\prime}\in
D_{-\mathbf{q}}}V_{\mathbf{q}}c_{\sigma^{\prime},\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}-\mathbf{k}}\delta_{\sigma^{\prime},\sigma}.\end{split}$
(S7)

The first term is simply a density-density Hamiltonian for the projected
particle density operator $\overline{\rho}_{\mathbf{q}}$ defined in Eq. (3),
and it is denoted in the main text by $\mathcal{H}_{\rho-\rho}$. The second
term is the single particle dispersion we wish to calculate here. The
dispersion contains a summation over momenta which is restricted such that
$\left(|\mathbf{k}|\leq k_{0}\right)\cap\left(|\mathbf{k}+\mathbf{q}|\leq
k_{0}\right)$ and $\left(|\mathbf{k}^{\prime}|\leq
k_{0}\right)\cap\left(|\mathbf{k}^{\prime}-\mathbf{q}|\leq k_{0}\right)$. The
delta function $\delta_{\mathbf{q},\mathbf{k}^{\prime}-\mathbf{k}}$ allows us
to make the following change in the summation order:

$-\frac{1}{2A}\sum_{\sigma=\uparrow,\downarrow}\sum_{\mathbf{q}}\sum_{\mathbf{k}\in
D_{\mathbf{q}}}\sum_{\mathbf{k}^{\prime}\in
D_{-\mathbf{q}}}V_{\mathbf{q}}c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}-\mathbf{k}}=-\frac{1}{2A}\sum_{\sigma=\uparrow,\downarrow}\sum_{\mathbf{k}\in
D_{0}}\sum_{\mathbf{k}^{\prime}\in
D_{0}}\sum_{\mathbf{q}}V_{\mathbf{q}}c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}-\mathbf{k}}.$
(S8)

Performing the summation over $\mathbf{q}$, we find the dispersion:

$-\frac{1}{2A}\sum_{\sigma=\uparrow,\downarrow}\sum_{\mathbf{k}\in
D_{0}}\sum_{\mathbf{k}^{\prime}\in
D_{0}}\sum_{\mathbf{q}}V_{\mathbf{q}}c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}-\mathbf{k}}=-\sum_{\sigma=\uparrow,\downarrow}\sum_{\mathbf{k}\in
D_{0}}E_{0}(\mathbf{k})c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}},\quad
E_{0}(\mathbf{k})=\frac{1}{2A}\sum_{|\mathbf{k}^{\prime}|\leq
k_{0}}V_{\mathbf{k}-\mathbf{k}^{\prime}},$ (S9)

which is the result presented in Eq. (5).

We now proceed to explicitly calculate the dispersion $E_{0}(\mathbf{k})$ in
the thermodynamic limit. Rewriting the sum as an integral, we find:

$\displaystyle\begin{split}E_{0}(\mathbf{k})=\frac{1}{2A}\sum_{|\mathbf{k}^{\prime}|\leq
k_{0}}V_{\mathbf{k}-\mathbf{k}^{\prime}}&=\frac{1}{2}\int_{|\mathbf{k}^{\prime}|\leq
k_{0}}\frac{d^{2}k^{\prime}}{(2\pi)^{2}}V_{\mathbf{k}-\mathbf{k}^{\prime}}\\\
&=\frac{1}{2}\int_{|\mathbf{p}-\mathbf{k}|\leq
k_{0}}\frac{d^{2}p}{(2\pi)^{2}}V_{\mathbf{p}}\\\ &=\pi
e^{2}\int_{|\mathbf{p}-\mathbf{k}|\leq
k_{0}}\frac{d^{2}p}{(2\pi)^{2}}\frac{\tanh{(|\mathbf{p}|d)}}{|\mathbf{p}|}\\\
&=\pi
e^{2}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\int_{0}^{\sqrt{k_{0}^{2}-|\mathbf{k}|^{2}\sin^{2}{\theta}}+|\mathbf{k}|\cos{\theta}}\frac{pdp}{2\pi}\frac{\tanh(pd)}{p},\end{split}$
(S10)

and by performing the radial integration we obtain the result:

$E_{0}(k)=\frac{e^{2}k_{0}}{2k_{0}d}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\ln\cosh\left[\left(\sqrt{k_{0}^{2}-k^{2}\sin^{2}{\theta}}+k\cos{\theta}\right)d\right],$
(S11)

where we denoted $|\mathbf{k}|=k$, as the dispersion has rotational symmetry.

The result in Eq. (S11) is plotted in Fig. 3 in the main text. This integral
has a familiar form in the limit of unscreened interactions
$k_{0}d\rightarrow\infty$. We now present it along with an expansion in powers
of $(k_{0}d)^{-1}$. Before continuing, it is useful to switch to dimensionless
variables. Denote $f(x)=E_{0}(xk_{0})/(e^{2}k_{0})$ and
$g(x,\theta)=\sqrt{1-x^{2}\sin^{2}{\theta}}+x\cos{\theta}$ such that:

$f(x)=\frac{1}{2k_{0}d}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\ln\cosh\left[g(x,\theta)k_{0}d\right].$
(S12)

Evaluation of the logarithm reveals the leading order terms in
$(k_{0}d)^{-1}$:

$\displaystyle\begin{split}f(x)&=\frac{1}{2k_{0}d}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\left(-\ln{2}+g(x,\theta)k_{0}d+\ln{\left(1+e^{-2g(x,\theta)k_{0}d}\right)}\right)\\\
&=\frac{1}{2}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\sqrt{1-x^{2}\sin^{2}{\theta}}+\frac{1}{2k_{0}d}\left(-\ln{2}+\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\ln{\left(1+e^{-2g(x,\theta)k_{0}d}\right)}\right)\\\
&=\frac{1}{\pi}\tilde{E}\left(\frac{\pi}{2},x\right)+\frac{1}{2k_{0}d}\left(-\ln{2}+\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\ln{\left(1+e^{-2g(x,\theta)k_{0}d}\right)}\right),\\\
\end{split}$ (S13)

where $\tilde{E}\left(\frac{\pi}{2},x\right)$ is the complete elliptic
integral of the second kind.

The first term in Eq. (S13) is a closed form solution for $f(x)$ to zeroth
order in $(k_{0}d)^{-1}$. For the next order in $(k_{0}d)^{-1}$, we need to
evaluate the integral in the third term and determine its behavior in the
limit of $k_{0}d\rightarrow\infty$. To do so, we separate our discussion into
two regimes: $x<1$ and $x=1$. For $x<1$ we have $g(x,\theta)>0$ for all
$\theta$, while for $x=1$, however, $g(x,\theta)>0$ only for
$\frac{\pi}{2}<|\theta|<\pi$, otherwise it vanishes. This implies the
following limits:

$\lim_{k_{0}d\rightarrow\infty}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\ln{\left(1+e^{-2g(x,\theta)k_{0}d}\right)}=\begin{cases}0,&x<1,\\\
\frac{\ln{2}}{2},&x=1.\end{cases}$ (S14)

We conclude that the power series for $f(x)$ does not converge uniformly:

$f(x)\approx\begin{cases}\frac{1}{\pi}\tilde{E}\left(\frac{\pi}{2},x\right)-\frac{\ln{2}}{2k_{0}d}+\mathcal{O}\left((k_{0}d)^{-2}\right),&x<1,\\\
\frac{1}{\pi}\tilde{E}\left(\frac{\pi}{2},x\right)-\frac{\ln{2}}{4k_{0}d}+\mathcal{O}\left((k_{0}d)^{-2}\right),&x=1.\end{cases}$
(S15)

Re-introducing the physical units of $e^{2}$ and $k_{0}$, we find the
dispersion to leading order in $(k_{0}d)^{-1}$ to be:

$E_{0}(\mathbf{k})=\begin{cases}e^{2}k_{0}\left(\frac{1}{\pi}\tilde{E}\left(\frac{\pi}{2},\frac{|\mathbf{k}|}{k_{0}}\right)-\frac{1}{2}\frac{\ln{2}}{k_{0}d}\right),&|\mathbf{k}|<k_{0},\\\
\\\
e^{2}k_{0}\left(\frac{1}{\pi}\tilde{E}\left(\frac{\pi}{2},1\right)-\frac{1}{4}\frac{\ln{2}}{k_{0}d}\right),&|\mathbf{k}|=k_{0},\end{cases}$
(S16)

and in the limit of $k_{0}d\rightarrow\infty$ we obtain the simple form
presented in Eq. (6). The seeming discontinuity in Eq. (S16) comes from the
divergence of the derivative of
$\tilde{E}\left(\frac{\pi}{2},\frac{|\mathbf{k}|}{k_{0}}\right)$ at
$|\mathbf{k}|=k_{0}$.

##  S3 Collective particle-hole excitations

In this section we examine collective particle-hole excitations with momentum
$\mathbf{Q}$ above the spin-polarized state, which by definition are magnons.
We consider both the strictly confining kinetic dispersion and a softened
version of the momentum confinement. For both cases we are interested in the
dispersion relation of the lowest-lying magnon excitation. We show how the
former results in a linear dispersion (and calculate its velocity), while for
the latter we find the usual quadratic magnon dispersion and calculate its
spin stiffness.

###  S3.1 The strictly confining kinetic dispersion

####  S3.1.1 Constructing the collective excitation matrix

The collective particle-hole excitations are eigenstates of
$\mathcal{H}_{\text{int}}$ which are spanned by single particle-hole excited
states denoted by $\ket{\Psi_{\text{ph}}^{\mathbf{k},\mathbf{Q}}}$ and defined
in Eq. (10). In what follows we derive the matrix in Eq. (11) whose
eigenstates are the collective excitations. We therefore consider the
following matrix element:

$\displaystyle\begin{split}[H(\mathbf{Q})]_{\mathbf{k},\mathbf{k}^{\prime}}&=\bra{\Psi_{\text{ph}}^{\mathbf{k}^{\prime},\mathbf{Q}}}\mathcal{H}_{\text{int}}\ket{\Psi_{\text{ph}}^{\mathbf{k},\mathbf{Q}}}\\\
&=\bra{\Psi_{\text{SP}}}c_{\downarrow,\mathbf{k}^{\prime}}^{\dagger}c_{\uparrow,\mathbf{k}^{\prime}+\mathbf{Q}}\mathcal{H}_{\text{int}}c_{\uparrow,\mathbf{k}+\mathbf{Q}}^{\dagger}c_{\downarrow,\mathbf{k}}\ket{\Psi_{\text{SP}}},\end{split}$
(S17)

for $\mathbf{k},\mathbf{k}^{\prime}\in D_{\mathbf{Q}}$, and
$D_{\mathbf{Q}}=\\{\mathbf{k}\enspace|\enspace(|\mathbf{k}|\leq
k_{0})\cap(|\mathbf{k}+\mathbf{Q}|\leq k_{0})\\}$. The wavefunction
$\ket{\Psi_{\text{SP}}}=\prod_{\mathbf{k}\in
D_{0}}c_{\downarrow,\mathbf{k}}^{\dagger}\ket{0}$ is the one described in the
main text, and is obtained by completely filling the disk $D_{0}$ with spin
down electrons. The Hamiltonian $\mathcal{H}_{\text{int}}$ is defined in Eq.
(2), and is written more explicitly in Eq. (S6). From the explicit form of
$\mathcal{H}_{\text{int}}$ we see that the matrix element
$[H(\mathbf{Q})]_{\mathbf{k},\mathbf{k}^{\prime}}$ in Eq. (S17) is a sum over
contractions of the form:

$\displaystyle\begin{split}&\bra{\Psi_{\text{SP}}}c_{\downarrow,\mathbf{k}^{\prime}}^{\dagger}c_{\uparrow,\mathbf{k}^{\prime}+\mathbf{Q}}\left(c_{\sigma,\mathbf{p}+\mathbf{q}}^{\dagger}c_{\sigma^{\prime},\mathbf{p}^{\prime}-\mathbf{q}}^{\dagger}c_{\sigma^{\prime},\mathbf{p}^{\prime}}c_{\sigma,\mathbf{p}}\right)c_{\uparrow,\mathbf{k}+\mathbf{Q}}^{\dagger}c_{\downarrow,\mathbf{k}}\ket{\Psi_{\text{SP}}}\\\
&=1_{\mathbf{k}\in
D_{0}}\left[\delta_{\sigma,\downarrow}\delta_{\sigma^{\prime},\downarrow}\delta_{\mathbf{k},\mathbf{k}^{\prime}}1_{\mathbf{p}\in
D_{0}}1_{\mathbf{p}^{\prime}\in
D_{0}}(\delta_{\mathbf{q},0}-\delta_{\mathbf{q},\mathbf{p}^{\prime}-\mathbf{p}})\right.\\\
&+\delta_{\mathbf{k},\mathbf{k}^{\prime}}\delta_{\mathbf{q},0}\left(\delta_{\sigma,\downarrow}\left(\delta_{\sigma^{\prime},\uparrow}\delta_{\mathbf{p}^{\prime},\mathbf{k}+\mathbf{Q}}-\delta_{\sigma^{\prime},\downarrow}\delta_{\mathbf{p}^{\prime},\mathbf{k}}\right)1_{\mathbf{p}\in
D_{0}}+\delta_{\sigma^{\prime},\downarrow}\left(\delta_{\sigma,\uparrow}\delta_{\mathbf{p},\mathbf{k}+\mathbf{Q}}-\delta_{\sigma,\downarrow}\delta_{\mathbf{k},\mathbf{p}}\right)1_{\mathbf{p}^{\prime}\in
D_{0}}\right)\\\
&+\delta_{\mathbf{k},\mathbf{k}^{\prime}}\left(\delta_{\sigma,\downarrow}\delta_{\sigma^{\prime},\downarrow}\delta_{\mathbf{q},\mathbf{p}^{\prime}-\mathbf{k}}\delta_{\mathbf{k},\mathbf{p}}+\delta_{\sigma,\downarrow}\delta_{\sigma^{\prime},\downarrow}\delta_{\mathbf{q},\mathbf{k}-\mathbf{p}}\delta_{\mathbf{k},\mathbf{p}^{\prime}}\right)1_{\mathbf{p}\in
D_{0}}1_{\mathbf{p}^{\prime}\in D_{0}}\\\
&\left.-\delta_{\sigma,\uparrow}\delta_{\sigma^{\prime},\downarrow}\delta_{\mathbf{p},\mathbf{k}+\mathbf{Q}}\delta_{\mathbf{q},\mathbf{k}^{\prime}-\mathbf{k}}\delta_{\mathbf{k}^{\prime},\mathbf{p}^{\prime}}1_{\mathbf{p}^{\prime}\in
D_{0}}-\delta_{\sigma,\downarrow}\delta_{\sigma^{\prime},\uparrow}\delta_{\mathbf{p}^{\prime},\mathbf{k}+\mathbf{Q}}\delta_{\mathbf{q},\mathbf{k}-\mathbf{k}^{\prime}}\delta_{\mathbf{k}^{\prime},\mathbf{p}}1_{\mathbf{p}\in
D_{0}}\right].\end{split}$ (S18)

The matrix element is calculated by carefully performing the summation over
Eq. (S18) as constrained by Eq. (S6) for $\mathbf{k},\mathbf{k}^{\prime}\in
D_{\mathbf{Q}}$:

$\bra{\Psi_{\text{ph}}^{\mathbf{k}^{\prime},\mathbf{Q}}}\mathcal{H}_{\text{int}}\ket{\Psi_{\text{ph}}^{\mathbf{k},\mathbf{Q}}}=\delta_{\mathbf{k},\mathbf{k}^{\prime}}\left(E_{\text{SP}}+\frac{1}{A}\sum_{\mathbf{p}\in
D_{0}}V_{\mathbf{k}-\mathbf{p}}\right)-\frac{1}{A}V_{\mathbf{k}-\mathbf{k}^{\prime}},$
(S19)

where we identified the spin polarized energy
$E_{\text{SP}}=\braket{\Psi_{\text{SP}}}{\mathcal{H}_{\text{int}}}{\Psi_{\text{SP}}}=\frac{1}{2A}\sum_{\mathbf{p},\mathbf{p}^{\prime}\in
D_{0}}\left(V_{0}-V_{\mathbf{p}-\mathbf{p}^{\prime}}\right)$.

The collective excitation matrix is confined to
$\mathbf{k},\mathbf{k}^{\prime}\in D_{\mathbf{Q}}$ and is given by:

$[H(\mathbf{Q})]_{\mathbf{k},\mathbf{k}^{\prime}}=\delta_{\mathbf{k},\mathbf{k}^{\prime}}\left(E_{\text{SP}}+\frac{1}{A}\sum_{\mathbf{p}\in
D_{0}}V_{\mathbf{k}-\mathbf{p}}\right)-\frac{1}{A}V_{\mathbf{k}-\mathbf{k}^{\prime}},\quad\mathbf{k},\mathbf{k}^{\prime}\in
D_{\mathbf{Q}}.$ (S20)

####  S3.1.2 Excitation spectrum bounds

Diagonalizing the matrix in Eq. (S20) can be done numerically for any
$\mathbf{Q}$, and the results are given in Fig. S1. Analytically, we can
derive bounds on parts of the excitation spectrum which will be useful for the
analytical structure of the lowest-lying eigenvalue in the vicinity of
$\mathbf{Q}=0$. We demonstrate the existence of a gapless particle-hole
excitation mode, which is expected, and prove that there is precisely one such
mode.

We are interested in the excitation energies of Eq. (S20) with respect to the
ground state energy. We therefore shift our reference energy by defining:

$\tilde{H}(\mathbf{Q})=H(\mathbf{Q})-E_{SP},$ (S21)

and denote its corresponding eigensystem ($n=1,\dots,|D_{\mathbf{Q}}|$) by:

$\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(n)}(\mathbf{Q})}=E_{\text{ph},n}(\mathbf{Q})\ket{\Psi_{ph}^{(n)}(\mathbf{Q})},\quad
E_{\text{ph},1}(\mathbf{Q})\leq E_{\text{ph},2}(\mathbf{Q})\leq\dots\leq
E_{\text{ph},|D_{\mathbf{Q}}|}(\mathbf{Q}).$ (S22)

In what follows we will find upper and lower bounds on the excitation energies
$E_{\text{ph},n}(\mathbf{Q})$, and show that the only excitation mode that has
zero energy is at $n=1$ and $\mathbf{Q}=0$.

The matrix $\tilde{H}(\mathbf{Q})$ in Eq. (S21) is a Hermitian diagonally
dominant matrix with positive diagonal entries, as can be easily verified from
Eq. (S20). Using the Gershgorin circle theorem[25, 26], it is easy to show
that the eigenvalues of $\tilde{H}(\mathbf{Q})$,
$E_{\text{ph},n}(\mathbf{Q})$, are bounded by:

$\min_{\mathbf{k}\in D_{\mathbf{Q}}}\\{\frac{1}{A}\sum_{\mathbf{k}^{\prime}\in
D_{0}}V_{\mathbf{k}-\mathbf{k}^{\prime}}-\frac{1}{A}\sum_{\begin{subarray}{c}\mathbf{k}^{\prime}\in
D_{\mathbf{Q}}\\\
\mathbf{k}^{\prime}\neq\mathbf{k}\end{subarray}}V_{\mathbf{k}-\mathbf{k}^{\prime}}\\}\leq
E_{\text{ph},n}(\mathbf{Q})\leq\max_{\mathbf{k}\in
D_{\mathbf{Q}}}\\{\frac{1}{A}\sum_{\mathbf{k}^{\prime}\in
D_{0}}V_{\mathbf{k}-\mathbf{k}^{\prime}}+\frac{1}{A}\sum_{\begin{subarray}{c}\mathbf{k}^{\prime}\in
D_{\mathbf{Q}}\\\
\mathbf{k}^{\prime}\neq\mathbf{k}\end{subarray}}V_{\mathbf{k}-\mathbf{k}^{\prime}}\\}.$
(S23)

For $\mathbf{Q}\neq 0$, this implies that $E_{\text{ph},n}(\mathbf{Q})>0$ are
strictly positive. For $\mathbf{Q}=0$, however, we find the following bounds
(in the thermodynamic limit):

$0\leq E_{\text{ph},n}(\mathbf{Q}=0)\leq
2\max_{\mathbf{k}}\\{\frac{1}{A}\sum_{\mathbf{k}^{\prime}\in
D_{0}}V_{\mathbf{k}-\mathbf{k}^{\prime}}\\}=\left.4E_{0}(\mathbf{k})\right|_{\mathbf{k}=0}\leq
2e^{2}k_{0},$ (S24)

where $E_{0}(\mathbf{k})$ is defined in Eq. (5).

The upper bound in Eq. (S24) will suffice for our purposes. The lower bound,
however, will be improved momentarily. Before we do so, we make a detour to
simplify notation. We denote the basis vectors
$\ket{\Psi_{\text{ph}}^{\mathbf{k},\mathbf{Q}}}$ by $\ket{\mathbf{k}}$ where
the momentum $\mathbf{Q}$ is implied by context, such that:

$\tilde{H}(\mathbf{Q})=\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
D_{\mathbf{Q}}}[\tilde{H}(\mathbf{Q})]_{\mathbf{k},\mathbf{k}^{\prime}}\ket{\mathbf{k}}\bra{\mathbf{k}^{\prime}}.$
(S25)

For $\mathbf{Q}=0$, there are $N=|D_{0}|$ such basis vectors. Returning to the
effort to tighten the lower bound, consider a normalized eigenvector
$\ket{\Psi_{ph}^{(n)}(0)}$ of $\tilde{H}(\mathbf{Q}=0)$. Its eigenvalue
$E_{\text{ph},n}(0)$ is given by:

$\displaystyle\begin{split}E_{\text{ph},n}(0)&=\braket{\Psi_{ph}^{(n)}(0)}{\tilde{H}(\mathbf{Q}=0)}{\Psi_{ph}^{(n)}(0)}\\\
&=\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
D_{0}}[\tilde{H}(\mathbf{Q}=0)]_{\mathbf{k},\mathbf{k}^{\prime}}\braket{\Psi_{ph}^{(n)}(0)}{\mathbf{k}}\braket{\mathbf{k}^{\prime}}{\Psi_{ph}^{(n)}(0)}\\\
&=\sum_{\mathbf{k}\in
D_{0}}[\tilde{H}(\mathbf{Q}=0)]_{\mathbf{k},\mathbf{k}}\braket{\Psi_{ph}^{(n)}(0)}{\mathbf{k}}\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}+\sum_{\begin{subarray}{c}\mathbf{k},\mathbf{k}^{\prime}\in
D_{0}\\\
\mathbf{k}\neq\mathbf{k}^{\prime}\end{subarray}}[\tilde{H}(\mathbf{Q}=0)]_{\mathbf{k},\mathbf{k}^{\prime}}\braket{\Psi_{ph}^{(n)}(0)}{\mathbf{k}}\braket{\mathbf{k}^{\prime}}{\Psi_{ph}^{(n)}(0)}\\\
&=\sum_{\begin{subarray}{c}\mathbf{k},\mathbf{k}^{\prime}\in D_{0}\\\
\mathbf{k}\neq\mathbf{k}^{\prime}\end{subarray}}[\tilde{H}(\mathbf{Q}=0)]_{\mathbf{k},\mathbf{k}^{\prime}}\left(\braket{\Psi_{ph}^{(n)}(0)}{\mathbf{k}}\braket{\mathbf{k}^{\prime}}{\Psi_{ph}^{(n)}(0)}-\braket{\Psi_{ph}^{(n)}(0)}{\mathbf{k}}\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}\right)\\\
&=\frac{1}{2}\sum_{\begin{subarray}{c}\mathbf{k},\mathbf{k}^{\prime}\in
D_{0}\\\
\mathbf{k}\neq\mathbf{k}^{\prime}\end{subarray}}\left(-[\tilde{H}(\mathbf{Q}=0)]_{\mathbf{k},\mathbf{k}^{\prime}}\right)\left|\braket{\mathbf{k}^{\prime}}{\Psi_{ph}^{(n)}(0)}-\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}\right|^{2}\\\
&\geq\frac{1}{2}\min_{\begin{subarray}{c}\mathbf{k},\mathbf{k}^{\prime}\in
D_{0}\\\
\mathbf{k}\neq\mathbf{k}^{\prime}\end{subarray}}\\{-[\tilde{H}(\mathbf{Q}=0)]_{\mathbf{k},\mathbf{k}^{\prime}}\\}\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
D_{0}}\left|\braket{\mathbf{k}^{\prime}}{\Psi_{ph}^{(n)}(0)}-\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}\right|^{2},\end{split}$
(S26)

where we used the fact that the diagonal is given by the sum of the rest of
the row, the symmetry of the matrix, and that its off-diagonal elements are
negative. The minimal off-diagonal element is given by:

$\min_{\begin{subarray}{c}\mathbf{k},\mathbf{k}^{\prime}\in D_{0}\\\
\mathbf{k}\neq\mathbf{k}^{\prime}\end{subarray}}\\{-[\tilde{H}(\mathbf{Q}=0)]_{\mathbf{k},\mathbf{k}^{\prime}}\\}=\frac{V_{2k_{0}}}{A},$
(S27)

which implies:

$E_{\text{ph},n}(0)\geq\frac{V_{2k_{0}}}{2A}\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
D_{0}}\left|\braket{\mathbf{k}^{\prime}}{\Psi_{ph}^{(n)}(0)}-\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}\right|^{2}.$
(S28)

The lower bound in Eq. (S28) can be computed exactly and therefore compared to
the lower bound found in Eq. (S24). For $n=1$ the two lower bounds coincide,
as they are satisfied by the uniform vector
$\ket{\Psi_{ph}^{(1)}(\mathbf{Q}=0)}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}\in
D_{0}}\ket{\mathbf{k}}$, i.e.
$\tilde{H}(\mathbf{Q}=0)\ket{\Psi_{ph}^{(1)}(0)}=E_{\text{ph},1}(0)\ket{\Psi_{ph}^{(1)}(\mathbf{Q}=0)}=0$.
For $n>1$, the bound in Eq. (S28) is stricter than the bound in Eq. (S24). To
show that, we make use of the orthogonality of eigenvectors of
$\tilde{H}(\mathbf{Q}=0)$:

$n>1:\quad\braket{\Psi_{ph}^{(1)}(0)}{\Psi_{ph}^{(n)}(0)}=0\implies\sum_{\mathbf{k}\in
D_{0}}\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}=0,$ (S29)

and use the above result to calculate the sum in Eq. (S28) exactly for $n>1$:

$\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
D_{0}}\left|\braket{\mathbf{k}^{\prime}}{\Psi_{ph}^{(n)}(0)}-\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}\right|^{2}=2\left(\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
D_{0}}\left|\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}\right|^{2}-\sum_{\mathbf{k}^{\prime}\in
D_{0}}\braket{\Psi_{ph}^{(n)}(0)}{\mathbf{k}^{\prime}}\sum_{\mathbf{k}\in
D_{0}}\braket{\mathbf{k}}{\Psi_{ph}^{(n)}(0)}\right)=2N,$ (S30)

where we used the normalization of $\ket{\Psi_{ph}^{(n)}(0)}$, $|D_{0}|=N$,
and Eq. (S29). Inserting this result back into Eq. (S28), we find the lower
bound on the gap to the next excitation energy, $E_{\text{ph},n}(0)\geq
V_{2k_{0}}(N/A)$. Taking the thermodynamic limit, at filling $\nu=1$ we have
$N/A=k_{0}^{2}/(4\pi)$, which implies $E_{\text{ph},n}(0)\geq e^{2}k_{0}/4$.
Combining this result with the upper bound in Eq. (S24), we find:

$\begin{cases}E_{\text{ph},1}(\mathbf{Q}=0)=0,&n=1,\\\
\frac{1}{4}e^{2}k_{0}\leq E_{\text{ph},n}(\mathbf{Q}=0)\leq
2e^{2}k_{0},&n>1.\end{cases}$ (S31)

We conclude there is precisely one excitation energy that starts at zero
energy, and the rest are gapped excitations.

Figure S1: Particle-hole excitation spectrum for $k_{0}d\rightarrow\infty$. In
figure (a) we plot the entire spectrum for the strict momentum confinement.
The lowest branch begins with a linear slope. The dashed line is linear whose
slope is the one theoretically predicted (computed for finite system size
$N$). Note that the number of finite excitations decreases with increasing $Q$
as the domain of our Hamiltonian becomes smaller. In figure (b) we compare the
lowest branch of the strict momentum confinement (and its linear slope) with
that of the softened confinement of different powers $N_{D}$ and the same
energy scale $U_{\text{K}}=0.1V_{k_{0}}N/A$. The crossover to linear slope
increases with $N_{D}$, and in the limit $N_{D}\rightarrow\infty$ it reaches
the origin.

####  S3.1.3 Dispersion of the lowest excitation energy

We are interested in calculating the evolution of the lowest excitation energy
with $\mathbf{Q}$. We will show here that it is linear in $|\mathbf{Q}|$, and
calculate its velocity as shown in Eq. (12). We will do so by bounding the
energy from both sides by two bounds with the same linear behavior at small
$|\mathbf{Q}|$.

Our starting point is the collective excitation matrix as defined in Eq.
(S21):

$\displaystyle\begin{split}\tilde{H}(\mathbf{Q})&=\sum_{\mathbf{k}\in
D_{\mathbf{Q}}}2E_{0}(\mathbf{k})\ket{\mathbf{k}}\bra{\mathbf{k}}-\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
D_{\mathbf{Q}}}\frac{1}{A}V_{\mathbf{k}-\mathbf{k}^{\prime}}\ket{\mathbf{k}}\bra{\mathbf{k}^{\prime}}.\end{split}$
(S32)

We decompose it as
$\tilde{H}(\mathbf{Q})=\tilde{H}_{0}(\mathbf{Q})+\tilde{H}_{1}(\mathbf{Q})$
such that:

$\displaystyle\tilde{H}_{0}(\mathbf{Q})$ $\displaystyle=\sum_{\mathbf{k}\in
D_{\mathbf{Q}}}\left(\frac{1}{A}\sum_{\mathbf{p}\in
D_{\mathbf{Q}}}V_{\mathbf{k}-\mathbf{p}}\right)\ket{\mathbf{k}}\bra{\mathbf{k}}-\sum_{\mathbf{k},\mathbf{k}^{\prime}\in
D_{\mathbf{Q}}}\frac{1}{A}V_{\mathbf{k}-\mathbf{k}^{\prime}}\ket{\mathbf{k}}\bra{\mathbf{k}^{\prime}},$
(S33) $\displaystyle\tilde{H}_{1}(\mathbf{Q})$
$\displaystyle=\sum_{\mathbf{k}\in
D_{\mathbf{Q}}}\left(\frac{1}{A}\sum_{\mathbf{p}\in
D_{0}-D_{\mathbf{Q}}}V_{\mathbf{k}-\mathbf{p}}\right)\ket{\mathbf{k}}\bra{\mathbf{k}}.$
(S34)

The purpose of this decomposition will be clear momentarily. The diagonal
elements of $\tilde{H}_{1}(\mathbf{Q})$ can be bounded as follows:

$\displaystyle\begin{split}\frac{1}{A}\sum_{\mathbf{p}\in
D_{0}-D_{\mathbf{Q}}}V_{\mathbf{k}-\mathbf{p}}&=\int_{\mathbf{p}\in
D_{0}-D_{\mathbf{Q}}}\frac{d^{2}p}{(2\pi)^{2}}V_{\mathbf{k}-\mathbf{p}}\\\
&\leq\int_{-\frac{|\mathbf{Q}|}{2}}^{\frac{|\mathbf{Q}|}{2}}\frac{dq_{x}}{2\pi}\int_{-k_{0}}^{k_{0}}\frac{dq_{y}}{2\pi}V_{\mathbf{q}}\\\
&\leq\frac{2}{\pi}e^{2}\int_{0}^{\frac{|\mathbf{Q}|}{2}}dq_{x}\int_{0}^{k_{0}}\frac{dq_{y}}{\sqrt{q_{x}^{2}+q_{y}^{2}}}\\\
&=\frac{2}{\pi}e^{2}\int_{0}^{\frac{|\mathbf{Q}|}{2}}dq_{x}\sinh^{-1}\left(\frac{k_{0}}{q_{x}}\right)\\\
&=\frac{1}{\pi}e^{2}\frac{|\mathbf{Q}|}{k_{0}}\log{\left(\frac{k_{0}}{|\mathbf{Q}|}\right)}+\mathcal{O}\left(\frac{|\mathbf{Q}|}{k_{0}}\right).\end{split}$
(S35)

In the bound above, the first inequality is a result of slicing the region
$D_{0}-D_{\mathbf{Q}}$ into slabs of length $|\mathbf{Q}|$ along $q_{x}$ (i.e.
$\hat{Q}$) and width $dq_{y}$, and shifting them along $q_{x}$ to maximize the
integrand. The second inequality is simply bounding the (possibly) screened
Coulomb potential by the unscreened potential with $k_{0}d\rightarrow\infty$.

Next, consider a uniform vector
$\ket{\phi_{0}(\mathbf{Q})}=\frac{1}{\sqrt{|D_{\mathbf{Q}}|}}\sum_{\mathbf{k}\in
D_{\mathbf{Q}}}\ket{\mathbf{k}}$. One can easily verify that
$\tilde{H}_{0}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}=0$. The expectation value
for the energy in this state is given by:

$\displaystyle\begin{split}\bra{\phi_{0}(\mathbf{Q})}\tilde{H}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}&=\bra{\phi_{0}(\mathbf{Q})}\tilde{H}_{1}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}\\\
&=\frac{2}{|D_{\mathbf{Q}}|}\sum_{\mathbf{k}\in
D_{\mathbf{Q}}}\left(\frac{1}{2A}\sum_{\mathbf{k}^{\prime}\in
D_{0}-D_{\mathbf{Q}}}V_{\mathbf{k}-\mathbf{k}^{\prime}}\right)\\\
&=\frac{2}{|D_{\mathbf{Q}}|}\sum_{\mathbf{k}^{\prime}\in
D_{0}-D_{\mathbf{Q}}}\left(E_{0}(\mathbf{k}^{\prime})-\frac{1}{2A}\sum_{\mathbf{k}\in
D_{0}-D_{\mathbf{Q}}}V_{\mathbf{k}-\mathbf{k}^{\prime}}\right),\end{split}$
(S36)

Keeping only the linear order in $|\mathbf{Q}|/k_{0}$, and making use of the
bound in Eq. (S35), we find:

$\displaystyle\begin{split}\bra{\phi_{0}(\mathbf{Q})}\tilde{H}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}&\approx
2\frac{|D_{0}-D_{\mathbf{Q}}|}{|D_{0}|}E_{0}(k_{0})+\mathcal{O}\left(\frac{|\mathbf{Q}|^{2}}{k_{0}^{2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right)\\\
&\approx 2\frac{2k_{0}|\mathbf{Q}|}{\pi
k_{0}^{2}}E_{0}(k_{0})+\mathcal{O}\left(\frac{|\mathbf{Q}|^{2}}{k_{0}^{2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right).\end{split}$
(S37)

By definition, the lowest eigenvalue is bounded from above by the expectation
value of $\tilde{H}(\mathbf{Q})$ with respect to any normalized wavefunction.
Therefore, we find that:

$E_{\text{ph},1}(\mathbf{Q})\leq\frac{4E_{0}(k_{0})}{\pi
k_{0}}|\mathbf{Q}|+\mathcal{O}\left(\frac{|\mathbf{Q}|^{2}}{k_{0}^{2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right)$
(S38)

We now show that the lower bound is equal to the upper bound in Eq. (S38).
Consider $\ket{\Psi_{ph}^{(1)}(\mathbf{Q})}$, the exact ground state of
$\tilde{H}(\mathbf{Q})$. By projecting $\ket{\Psi_{ph}^{(1)}(\mathbf{Q})}$
onto $\ket{\phi_{0}(\mathbf{Q})}$, we can generally decompose it as:

$\ket{\Psi_{ph}^{(1)}(\mathbf{Q})}=a_{\mathbf{Q}}\ket{\phi_{0}(\mathbf{Q})}+b_{\mathbf{Q}}\ket{\chi(\mathbf{Q})},$
(S39)

where the unknown wavefunction $\ket{\chi(\mathbf{Q})}$ is orthogonal to
$\ket{\phi_{0}(\mathbf{Q})}$, i.e.
$\braket{\phi_{0}(\mathbf{Q})}{\chi(\mathbf{Q})}=0$, and
$|a_{\mathbf{Q}}|^{2}+|b_{\mathbf{Q}}|^{2}=1$. Now, let us examine the ground
state energy:

$\displaystyle\begin{split}E_{\text{ph},1}(\mathbf{Q})&=\bra{\Psi_{ph}^{(1)}(\mathbf{Q})}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(\mathbf{Q})}\\\
&=\left[a_{\mathbf{Q}}^{*}\bra{\phi_{0}(\mathbf{Q})}+b_{\mathbf{Q}}^{*}\bra{\chi(\mathbf{Q})}\right]\left(\tilde{H}_{0}(\mathbf{Q})+\tilde{H}_{1}(\mathbf{Q})\right)\left[a_{\mathbf{Q}}\ket{\phi_{0}(\mathbf{Q})}+b_{\mathbf{Q}}\ket{\chi(\mathbf{Q})}\right]\\\
&=|a_{\mathbf{Q}}|^{2}\bra{\phi_{0}(\mathbf{Q})}\tilde{H}_{1}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}+|b_{\mathbf{Q}}|^{2}\bra{\chi(\mathbf{Q})}\left(\tilde{H}_{0}(\mathbf{Q})+\tilde{H}_{1}(\mathbf{Q})\right)\ket{\chi(\mathbf{Q})}+2\textrm{Re}\\{b_{\mathbf{Q}}^{*}a_{\mathbf{Q}}\bra{\chi(\mathbf{Q})}\tilde{H}_{1}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}\\}.\end{split}$
(S40)

Taking the limit $|\mathbf{Q}|/k_{0}\rightarrow 0$, and using (S35), we know
that all matrix elements of $\tilde{H}_{1}(\mathbf{Q})$ vanish. The only term
that does not necessarily vanish is:

$\lim_{|\mathbf{Q}|/k_{0}\rightarrow
0}E_{\text{ph},1}(\mathbf{Q})=\lim_{|\mathbf{Q}|/k_{0}\rightarrow
0}|b_{\mathbf{Q}}|^{2}\bra{\chi(\mathbf{Q})}\tilde{H}_{0}(\mathbf{Q})\ket{\chi(\mathbf{Q})}.$
(S41)

Using the same method as in Eqs. (S26) – (S30), we can write the lower bound
for the expectation value in Eq. (S41) as:

$\bra{\chi(\mathbf{Q})}\tilde{H}_{0}(\mathbf{Q})\ket{\chi(\mathbf{Q})}\geq|D_{\mathbf{Q}}|\min_{\begin{subarray}{c}\mathbf{k},\mathbf{k}^{\prime}\in
D_{\mathbf{Q}}\\\
\mathbf{k}\neq\mathbf{k}^{\prime}\end{subarray}}\\{-[\tilde{H}(\mathbf{Q}=0)]_{\mathbf{k},\mathbf{k}^{\prime}}\\}.$
(S42)

Inserting this bound back into Eq. (S41), we find:

$\lim_{|\mathbf{Q}|/k_{0}\rightarrow
0}E_{\text{ph},1}(\mathbf{Q})\geq\frac{1}{4}e^{2}k_{0}\lim_{|\mathbf{Q}|/k_{0}\rightarrow
0}|b_{\mathbf{Q}}|^{2}.$ (S43)

Combined with the upper bound on the lowest excitation energy found in Eq.
(S38), we find that:

$\lim_{|\mathbf{Q}|/k_{0}\rightarrow
0}\frac{E_{\text{ph},1}(\mathbf{Q})}{|\mathbf{Q}|/k_{0}}\leq\frac{4}{\pi}E_{0}(k_{0})\implies\lim_{|\mathbf{Q}|/k_{0}\rightarrow
0}\frac{|b_{\mathbf{Q}}|^{2}}{|\mathbf{Q}|/k_{0}}\leq\frac{16}{\pi}\frac{E_{0}(k_{0})}{e^{2}k_{0}}\implies|b_{\mathbf{Q}}|\sim\sqrt{\frac{|\mathbf{Q}|}{k_{0}}}.$
(S44)

We can now derive the lower bound. Continuing from Eq. (S40) and using the
positive-definiteness of $\tilde{H}(\mathbf{Q})$:

$\displaystyle\begin{split}E_{\text{ph},1}(\mathbf{Q})&\geq|a_{\mathbf{Q}}|^{2}\bra{\phi_{0}(\mathbf{Q})}\tilde{H}_{1}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}+2\textrm{Re}\\{b_{\mathbf{Q}}^{*}a_{\mathbf{Q}}\bra{\chi(\mathbf{Q})}\tilde{H}_{1}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}\\}\\\
&\geq\bra{\phi_{0}(\mathbf{Q})}\tilde{H}_{1}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}\left(1+\mathcal{O}\left(\frac{|\mathbf{Q}|}{k_{0}}\right)\right)+\mathcal{O}\left(\frac{|\mathbf{Q}|^{3/2}}{k_{0}^{3/2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right)\\\
&=\frac{4E_{0}(k_{0})}{\pi
k_{0}}|\mathbf{Q}|+\mathcal{O}\left(\frac{|\mathbf{Q}|^{3/2}}{k_{0}^{3/2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right),\end{split}$
(S45)

where in the last equality we inserted the expectation value as calculated in
Eq. (S37). Combining this lower bound with the upper bound in Eq. (S38), we
find that:

$\frac{4E_{0}(k_{0})}{\pi
k_{0}}|\mathbf{Q}|+\mathcal{O}\left(\frac{|\mathbf{Q}|^{3/2}}{k_{0}^{3/2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right)\leq
E_{\text{ph},1}(\mathbf{Q})\leq\frac{4E_{0}(k_{0})}{\pi
k_{0}}|\mathbf{Q}|+\mathcal{O}\left(\frac{|\mathbf{Q}|^{2}}{k_{0}^{2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right),$
(S46)

which implies that:

$E_{\text{ph},1}(\mathbf{Q})\approx\frac{4E_{0}(k_{0})}{\pi
k_{0}}|\mathbf{Q}|.$ (S47)

For $k_{0}d\rightarrow\infty$, the velocity of these excitations is given by
$v=\left(\frac{2}{\pi}\right)^{2}e^{2}=\left(\frac{2}{\pi}\right)^{2}\alpha
c$, where $\alpha\approx\frac{1}{137}$ is the fine-structure constant and $c$
is the speed of light. This linear slope is added to Fig. S1 by evaluating
$E_{0}(k_{0})$ for a finite system, and is found to be in agreement with the
numerically computed excitation spectrum.

The isotropic result above is derived for the idealized dispersion in Eq. (1).
One may naturally wonder what would happen upon the introduction of single-
particle dispersion to the disk. Would we still have a linearly dispersive
magnon? If so, what would be the velocity? Would it remain isotropic?

Let us assume our Hamiltonian is now of the form:

$\mathcal{H}_{\text{int}}\rightarrow\mathcal{H}_{\text{int}}+\sum_{\sigma=\uparrow,\downarrow}\sum_{|\mathbf{k}|\leq
k_{0}}E_{\text{disp}}(\mathbf{k})c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}},$
(S48)

where $E_{\text{disp}}(\mathbf{k})$ is a dispersion that is differentiable on
the disk with a bounded gradient. Note that $E_{0}(\mathbf{k})$ in Eq. (5) for
the case on an unscreened Coulomb interaction ($k_{0}d\rightarrow\infty$) does
not possess this property, as its derivative diverges as
$|\mathbf{k}|\rightarrow k_{0}$. We will address this dispersion separately.
The Hamiltonian $\tilde{H}(\mathbf{Q})$ defined in Eq. (S21) acquires a new
correction,

$\displaystyle\begin{split}\tilde{H}_{2}(\mathbf{Q})&=\sum_{\mathbf{k}\in
D_{\mathbf{Q}}}\left(E_{\text{disp}}(\mathbf{k}+\mathbf{Q})-E_{\text{disp}}(\mathbf{k})\right)\ket{\mathbf{k}}\bra{\mathbf{k}}\\\
&\approx\sum_{\mathbf{k}\in
D_{\mathbf{Q}}}\mathbf{v}_{\mathbf{k}}\cdot\mathbf{Q}\ket{\mathbf{k}}\bra{\mathbf{k}}+\mathcal{O}\left(\frac{|\mathbf{Q}|^{2}}{k_{0}^{2}}\right),\end{split}$
(S49)

such that
$\tilde{H}(\mathbf{Q})=\tilde{H}_{0}(\mathbf{Q})+\tilde{H}_{1}(\mathbf{Q})+\tilde{H}_{2}(\mathbf{Q})$,
and we defined $\mathbf{v}_{\mathbf{k}}=\frac{\partial
E_{\text{disp}}(\mathbf{k})}{\partial{\mathbf{k}}}$. Following the previous
arguments, the upper bound on the energy is given by:

$\displaystyle\begin{split}E_{\text{ph},1}(\mathbf{Q})&\leq\bra{\phi_{0}(\mathbf{Q})}\tilde{H}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}\\\
&=\bra{\phi_{0}(\mathbf{Q})}(\tilde{H}_{1}(\mathbf{Q})+\tilde{H}_{2}(\mathbf{Q}))\ket{\phi_{0}(\mathbf{Q})}\\\
&\approx\frac{4E_{0}(k_{0})}{\pi
k_{0}}|\mathbf{Q}|+\frac{1}{|D_{0}|}\sum_{\mathbf{k}\in
D_{0}}\mathbf{v}_{\mathbf{k}}\cdot\mathbf{Q}+\mathcal{O}\left(\frac{|\mathbf{Q}|^{2}}{k_{0}^{2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right).\end{split}$
(S50)

The lower bound is derived similarly,

$\displaystyle\begin{split}E_{\text{ph},1}(\mathbf{Q})&\geq\frac{4E_{0}(k_{0})}{\pi
k_{0}}|\mathbf{Q}|+\frac{1}{|D_{0}|}\sum_{\mathbf{k}\in
D_{0}}\mathbf{v}_{\mathbf{k}}\cdot\mathbf{Q}+\mathcal{O}\left(\frac{|\mathbf{Q}|^{3/2}}{k_{0}^{3/2}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}\right),\end{split}$
(S51)

such that the particle-hole dispersion is given by:

$E_{\text{ph},1}(\mathbf{Q})\approx\frac{4E_{0}(k_{0})}{\pi
k_{0}}|\mathbf{Q}|+\frac{1}{|D_{0}|}\sum_{\mathbf{k}\in
D_{0}}\mathbf{v}_{\mathbf{k}}\cdot\mathbf{Q}.$ (S52)

We find that the smooth single-particle dispersion adds its average velocity
on the disk to the isotropic velocity given by the strict confinement of the
Coulomb interaction. Therefore, our result in Eq. (12) holds for any smooth
dispersion whose average velocity on the disk is zero.

Finally, consider adding the single-particle dispersion $E_{0}(\mathbf{k})$ in
Eq. (5), such that:

$\displaystyle\begin{split}\tilde{H}_{2}(\mathbf{Q})&=\sum_{\mathbf{k}\in
D_{\mathbf{Q}}}\left(E_{0}(\mathbf{k}+\mathbf{Q})-E_{0}(\mathbf{k})\right)\ket{\mathbf{k}}\bra{\mathbf{k}}.\end{split}$
(S53)

The difficulty with blindly repeating the previous arguments occurs for the
unscreened Coulomb interaction, where logarithmically singular terms of the
form $\frac{|\mathbf{Q}|}{k_{0}}\log{\left(\frac{|\mathbf{Q}|}{k_{0}}\right)}$
appear in the small $|\mathbf{Q}|/k_{0}$ expansion of this Hamiltonian. Using
the rotational symmetry of $E_{0}(\mathbf{k})$ it is easy to show that, in
this case,
$\bra{\phi_{0}(\mathbf{Q})}\tilde{H}_{2}(\mathbf{Q})\ket{\phi_{0}(\mathbf{Q})}=0$.
Consequently, the upper bound on $E_{\text{ph},1}(\mathbf{Q})$ is unchanged by
this dispersion. The lower bound is also unchanged, as the logarithmically
singular terms in $\tilde{H}_{2}(\mathbf{Q})$ play a similar role to such
terms that are already present in $\tilde{H}_{1}(\mathbf{Q})$. We therefore
conclude that that Eq. (12) holds even when one adds the single-particle
dispersion $E_{0}(\mathbf{k})$ to Eq. (1), i.e. for $\mathcal{H}_{\rho-\rho}$
as defined in Eq. (4).

###  S3.2 The softened momentum confinement

So far we have only considered the strict momentum confinement of Eq. (1).
Now, we would like to consider a softened momentum confinement,
$\mathcal{H}_{\text{K}}$, as defined in Eq. (13), and calculate the dispersion
of the particle-hole excitations above the same spin-polarized state
$\ket{\Psi_{\text{SP}}}$. In Fig. S1(b) we plot the lowest lying excitation
for various powers of $N_{D}$ (the strength of the soft confinement), and
compare the result to the limit of $N_{D}\rightarrow\infty$ obtained above. We
emphasize that the state $\ket{\Psi_{\text{SP}}}$ is not the ground state of
the Hamiltonian $\mathcal{H}_{\text{K}}$, but for large enough $N_{D}$ we
expect it to be a good approximation.

The softened momentum confinement means we allow particle excitations beyond
the disk. Consequently, the Coulomb interaction is no longer restricted to the
disk, thus we denote it by $\mathcal{V}_{\text{int}}$ to distinguish it from
$\mathcal{H}_{\text{int}}$:

$\mathcal{V}_{\text{int}}=\frac{1}{2A}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}V_{\mathbf{q}}c_{\sigma^{\prime},\mathbf{k}^{\prime}-\mathbf{q}}^{\dagger}c_{\sigma,\mathbf{k}+\mathbf{q}}^{\dagger}c_{\sigma,\mathbf{k}}c_{\sigma^{\prime},\mathbf{k}^{\prime}}.$
(S54)

The first step is to compute the matrix elements of
$\mathcal{V}_{\text{int}}+\mathcal{H}_{\text{K}}$ between different
wavefunctions $\ket{\Psi_{\text{ph}}^{\mathbf{k},\mathbf{Q}}}$. For the
softened confinement, the momentum $\mathbf{k}$ is the entire disk for any
$\mathbf{Q}$, and we find the following matrix for
$\mathbf{k},\mathbf{k}^{\prime}\in D_{0}$:

$\displaystyle\begin{split}[H(\mathbf{Q})]_{\mathbf{k},\mathbf{k}^{\prime}}&=\bra{\Psi_{\text{ph}}^{\mathbf{k}^{\prime},\mathbf{Q}}}(\mathcal{V}_{\text{int}}+\mathcal{H}_{\text{K}})\ket{\Psi_{\text{ph}}^{\mathbf{k},\mathbf{Q}}}\\\
&=\delta_{\mathbf{k},\mathbf{k}^{\prime}}\left[E_{\text{SP}}+\frac{1}{A}\sum_{|\mathbf{p}|\leq
k_{0}}V_{\mathbf{k}-\mathbf{p}}+U_{\text{K}}\left(\left(\frac{|\mathbf{k}+\mathbf{Q}|}{k_{0}}\right)^{N_{D}}-\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{N_{D}}\right)\right]-\frac{1}{A}V_{\mathbf{k}-\mathbf{k}^{\prime}}\end{split}$
(S55)

We retain the notation of the Hamiltonian, the eigensystem and the basis (Eqs.
(S21),(S22) and Eq. (S25), respectively), as the two problems can, for the
most part, be considered aspects of the same one with a different choice of a
parameter (finite vs infinite $N_{D}$). The only relevant difference is the
number of finite eigenvalues, which is now unconstrained by $\mathbf{Q}$
(i.e., $n=1,\dots,N$ where $N=|D_{0}|$):

$\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(n)}(\mathbf{Q})}=E_{\text{ph},n}(\mathbf{Q})\ket{\Psi_{ph}^{(n)}(\mathbf{Q})},\quad
E_{\text{ph},1}(\mathbf{Q})\leq E_{\text{ph},2}(\mathbf{Q})\leq\dots\leq
E_{\text{ph},N}(\mathbf{Q}).$ (S56)

We are interested in evaluating $E_{\text{ph},1}(\mathbf{Q})$ to leading order
in $|\mathbf{Q}|/k_{0}$. To second order in perturbation theory, we find:

$\displaystyle\begin{split}E_{\text{ph},1}(\mathbf{Q})&\approx
E_{\text{ph},1}(0)+\bra{\Psi_{ph}^{(1)}(0)}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}+\sum_{n=2}^{N}\frac{\left|\bra{\Psi_{ph}^{(n)}(0)}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}\right|^{2}}{E_{\text{ph},1}(0)-E_{\text{ph},n}(0)},\end{split}$
(S57)

where $\tilde{H}(\mathbf{Q})$ should be expanded up to second order in
$|\mathbf{Q}|/k_{0}$. The exact ground state of at $\mathbf{Q}=0$ is known
exactly:

$\ket{\Psi_{ph}^{(1)}(0)}=\frac{1}{\sqrt{N}}\sum_{|\mathbf{k}|\leq
k_{0}}\ket{\mathbf{k}},$ (S58)

and we have $E_{\text{ph},1}(0)=\tilde{H}(0)\ket{\Psi_{ph}^{(1)}(0)}=0$.
Therefore, the zeroth order in $|\mathbf{Q}|/k_{0}$ vanishes, as expected. For
the first and second order in $|\mathbf{Q}|/k_{0}$, we expand the Hamiltonian
as:

$\displaystyle\begin{split}[\tilde{H}(\mathbf{Q})]_{\mathbf{k},\mathbf{k}^{\prime}}&\approx[\tilde{H}(0)]_{\mathbf{k},\mathbf{k}^{\prime}}+\delta_{\mathbf{k},\mathbf{k}^{\prime}}U_{\text{K}}\frac{N_{D}}{2}\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{N_{D}-4}\frac{2|\mathbf{k}|^{2}\mathbf{k}\cdot\mathbf{Q}+|\mathbf{k}|^{2}|\mathbf{Q}|^{2}+(N_{D}-1)\left(\mathbf{k}\cdot\mathbf{Q}\right)^{2}}{k_{0}^{4}}+\mathcal{O}\left(\frac{|\mathbf{Q}|^{3}}{k_{0}^{3}}\right).\end{split}$
(S59)

Applying this expansion to the exact ground state at $\mathbf{Q}=0$, we find:

$\displaystyle\begin{split}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}&\approx\frac{1}{\sqrt{N}}U_{\text{K}}\frac{N_{D}}{2}\sum_{|\mathbf{k}|\leq
k_{0}}\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{N_{D}-4}\frac{2|\mathbf{k}|^{2}\mathbf{k}\cdot\mathbf{Q}+|\mathbf{k}|^{2}|\mathbf{Q}|^{2}+(N_{D}-1)\left(\mathbf{k}\cdot\mathbf{Q}\right)^{2}}{k_{0}^{4}}\ket{\mathbf{k}}.\end{split}$
(S60)

Equipped with the result above, we can evaluate the perturbative corrections
to the ground state energy.

The first order correction term, which has also $|\mathbf{Q}|^{2}/k_{0}^{2}$
terms, is given by:

$\displaystyle\begin{split}\bra{\Psi_{ph}^{(1)}(0)}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}&\approx\frac{1}{N}U_{\text{K}}\frac{N_{D}}{2}\sum_{|\mathbf{k}|\leq
k_{0}}\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{N_{D}-4}\frac{2|\mathbf{k}|^{2}\mathbf{k}\cdot\mathbf{Q}+|\mathbf{k}|^{2}|\mathbf{Q}|^{2}+(N_{D}-1)\left(\mathbf{k}\cdot\mathbf{Q}\right)^{2}}{k_{0}^{4}}\\\
&\approx\frac{4\pi}{k_{0}^{2}}U_{\text{K}}\frac{N_{D}}{2}\int_{0}^{k_{0}}\frac{kdk}{2\pi}\int_{0}^{2\pi}\frac{d\theta}{2\pi}\left(\frac{k}{k_{0}}\right)^{N_{D}-4}\frac{k^{2}Q^{2}+(N_{D}-1)k^{2}Q^{2}\cos^{2}{\theta}}{k_{0}^{4}}\\\
&=\frac{1}{2}U_{\text{K}}\left(N_{D}+1\right)\frac{Q^{2}}{k_{0}^{2}},\end{split}$
(S61)

where we used $\frac{N}{A}=\frac{k_{0}^{2}}{4\pi}$ as the density of electrons
in the system.

The second order correction is strictly negative. We will now bound its
magnitude. For $\mathbf{Q}=0$, the finite and infinite $N_{D}$ cases coincide,
as we consider the same ground state. Therefore, the energy spectrum is
bounded between the same two values as in Eq. (S31):

$\displaystyle\frac{1}{2e^{2}k_{0}}\leq\frac{1}{|E_{\text{ph},1}(0)-E_{\text{ph},n}(0)|}\leq\frac{4}{e^{2}k_{0}}$
(S62)

Therefore, up to a constant $\frac{1}{2}\leq C\leq 4$, we find the second
order correction to be:

$\displaystyle\begin{split}\sum_{n=2}^{N}\frac{\left|\bra{\Psi_{ph}^{(n)}(0)}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}\right|^{2}}{|E_{\text{ph},1}(0)-E_{\text{ph},n}(0)|}&=\frac{C}{e^{2}k_{0}}\sum_{n=2}^{N}\left|\bra{\Psi_{ph}^{(n)}(0)}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}\right|^{2}\\\
&=\frac{C}{e^{2}k_{0}}\left(\left|\left|\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}\right|\right|^{2}-\left|\bra{\Psi_{ph}^{(1)}(0)}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}\right|^{2}\right)\end{split}$
(S63)

Including terms up to second order in $|\mathbf{Q}|/k_{0}$, only the first
term has a $Q^{2}$ component:

$\displaystyle\begin{split}\sum_{n=2}^{N}\frac{\left|\bra{\Psi_{ph}^{(n)}(0)}\tilde{H}(\mathbf{Q})\ket{\Psi_{ph}^{(1)}(0)}\right|^{2}}{|E_{\text{ph},1}(0)-E_{\text{ph},n}(0)|}&\approx\frac{C}{N}\frac{U_{\text{K}}^{2}}{e^{2}k_{0}}N_{D}^{2}\sum_{|\mathbf{k}|\leq
k_{0}}\left(\frac{|\mathbf{k}|}{k_{0}}\right)^{2N_{D}-4}\frac{(\mathbf{k}\cdot\mathbf{Q})^{2}}{k_{0}^{4}}\\\
&\approx\frac{C}{2}\frac{U_{\text{K}}^{2}}{e^{2}k_{0}}N_{D}\frac{Q^{2}}{k_{0}^{2}}\end{split}$
(S64)

where $\frac{1}{2}\leq C\leq 4$.

Taking $N_{D}\gg 1$, we find that for small enough $Q/k_{0}$:

$E_{\text{ph},1}(\mathbf{Q})\approx\frac{1}{2}U_{\text{K}}\left(1-C\frac{U_{\text{K}}}{e^{2}k_{0}}\right)N_{D}\frac{Q^{2}}{k_{0}^{2}}$
(S65)

Therefore, if $U_{\text{K}}$ is small compared to the next particle-hole
excitation energy, i.e. $U_{\text{K}}\ll e^{2}k_{0}$, the second term in the
parenthesis becomes unimportant, and we obtain the positive and divergent spin
stiffness in Eq. (14). The divergence can be seen numerically in Fig. S1 where
$E_{\text{ph},1}(\mathbf{Q})$ is computed numerically for several values of
$N_{D}$. We remind the reader that while the particle-hole excitation here is
calculated with respect to a state that is not the ground state of the system,
we expect that the statement above will be qualitatively correct for large
enough $N_{D}$.

##  S4 Numerical Analysis

###  S4.1 Discretization of momentum space

We begin by describing our finite-size model. Consider a finite 2D system of
area $A$ with periodic boundary conditions. In an effort to utilize as many
symmetries as possible, we choose the simulation box to be invariant under
$\frac{\pi}{3}$ rotations. In real space, this box is a rhombus with an angle
of $\frac{\pi}{3}$ and a side of length $L=\sqrt{\frac{2}{\sqrt{3}}A}$.
Correspondingly, the momentum space is discretized as a triangular lattice
with periodicity $\Delta k=\frac{4\pi}{\sqrt{3}L}$. We are interested in
considering only momenta up to a cutoff $k_{0}$. By adjusting the area $A$, we
set $k_{0}=R\cdot\Delta k$ for an integer $R$ which determines the number of
grid points up to the cutoff. The system area is therefore set to be
$A=\frac{2}{\sqrt{3}}\frac{(2\pi)^{2}}{k_{0}^{2}}R^{2}$. For $1\leq R\leq 6$,
the number of grid points is $N=3R^{2}+3R+1$, such that at $R=1$ and $R=2$ we
obtain $N=7$ and $N=19$ grid points, respectively. In the thermodynamic limit
(taking $A\rightarrow\infty$, which implies $R\rightarrow\infty$), we find
$N\approx\frac{2\pi}{\sqrt{3}}R^{2}$. For $R\geq 3$ the number of grid points
is too large for exact diagonalization, as the Hilbert space dimension is
exponentially large. The particle-hole excitation calculation has a Hilbert
space dimension that is linear in $N$, and there we use $R=6$. At filling
factor $\nu=1$ we place one electron at each grid point, and the electron
density $n=\frac{N}{A}$ is given by:

$\frac{N}{A}=\begin{cases}\frac{\sqrt{3}}{2}\frac{3R^{2}+3R+1}{R^{2}}\frac{k_{0}^{2}}{(2\pi)^{2}},&1\leq
R\leq 6\\\ \pi\frac{k_{0}^{2}}{(2\pi)^{2}},&R\rightarrow\infty.\end{cases}$
(S66)

We use the density factor obtained for small $R$ when presenting numerical
results in Fig. S3, and the one for infinite $R$ when discussing analytical
results.

###  S4.2 Exact diagonalization

The numerical calculations were performed on the WEXAC cluster at the Weizmann
Institute of Science. The exact diagonalization is implemented using the
matrix-free implicitly restarted Lanczos method, converging $7$ lowest
eigenvalues to an accuracy of $10^{-5}$. The algorithm is implemented using
the ARPACK [27] library as maintained by opencollab ARPACK-NG [28].

####  S4.2.1 Procedure to obtain the ground state’s energy and spin
polarization

Next, we describe the procedure used to exactly diagonalize the Hamiltonians
in Eq. (2) and Eq. (15) for the finite system of $N$ grid points defined
above. We are interested in finding the ground state energy at a filling
factor of $\nu=1$, determining its total spin, and repeating the calculation
for different effective masses. The first step is to separate the Hilbert
space into small sectors that are invariant subspaces of the Hamiltonian. This
reduction in size will render the exact diagonalization tractable. The second
step will be to extract the information about the spin polarization.

The entire Hilbert space $\mathcal{H}$ associated with the Hamiltonians of
interest consists of two orbitals at each grid point, one for each spin
projection, such that $\dim\mathcal{H}=2^{2N}$. Focusing on a fixed number of
particles equal to the number of grid points $N$, i.e. for filling factor
$\nu=1$, the relevant Hilbert space is reduced to
$\dim\mathcal{H}={{2N}\choose{N}}$. We can further decrease our Hilbert space
dimension by considering sectors of total momentum,
$\mathbf{k}_{\text{tot}}=\sum_{\sigma}\sum_{|\mathbf{k}|\leq
k_{0}}\mathbf{k}c_{\sigma,\mathbf{k}}^{\dagger}c_{\sigma,\mathbf{k}}$ , and of
total spin projection, $S_{z}=\frac{1}{2}\sum_{|\mathbf{k}|\leq
k_{0}}(c_{\uparrow,\mathbf{k}}^{\dagger}c_{\uparrow,\mathbf{k}}-c_{\downarrow,\mathbf{k}}^{\dagger}c_{\downarrow,\mathbf{k}})$.
The largest sector is the
$\left(\mathbf{k}_{\text{tot}},S_{z}\right)=\left(0,\frac{1}{2}\right)$ sector
which for $N=19$ has $\dim\mathcal{H}\approx 9\cdot 10^{7}$. While still very
large, this sector is small enough to allow for its exact diagonalization.

There are more symmetries we wish to exploit - the point group symmetries of
the hexagonal lattice (rotation and mirror), and the $SU(2)$ symmetry
associated with rotations of the total spin. These symmetries have important
yet very different consequences on our numerical calculations. The point group
symmetries imply that there is complete redundancy in information from
different sectors of $\mathbf{k}_{\text{tot}}$ that are related by symmetry,
thereby the number of distinct $\mathbf{k}_{\text{tot}}$ sectors we need to
consider is significantly reduced. In a similar manner, the spin rotation
symmetry allows us to consider only spin sectors with $S_{z}\geq 0$ (actually,
$S_{z}\geq\frac{1}{2}$, as we have an odd number of particles), as the
negative values of $S_{z}$ contain the same information as their positive
counterparts.

However, unlike the point group symmetries, the spin rotational degree of
freedom allows us to explore the total spin of the eigenstates. To see that,
consider a state with total spin $S$. Due to $SU(2)$ symmetry, this state will
be $2S+1$ degenerate and will appear in sectors with $S_{z}=-S,\dots,S$. Since
total spin is a conserved number, any eigenstate with an odd number of
particles necessarily is degenerate with a representative in the
$S_{z}=\frac{1}{2}$ sector. Therefore, the $S_{z}=\frac{1}{2}$ sector contains
the entire eigenvalue spectrum for an odd number of particles. Sweeping
different $S_{z}$ sectors allows us to identify the total spin $S$ of every
state as the maximal $S_{z}$ sector in which it appears.

We summarize the argument above with the following recipe to find the ground
state and its spin polarization:

  1. 1.

List all sectors of $S_{z}=\frac{1}{2}$ and $\mathbf{k}_{\text{tot}}$ that are
unrelated by the point group symmetries.

  2. 2.

Diagonalize each sector to find its lowest eigenvalue.

  3. 3.

Identify the momentum sector hosting the lowest eigenvalue across all the
sectors. This momentum is the ground state momentum (up to point group
symmetries), and the eigenvalue is the ground state energy.

  4. 4.

Repeat the exact diagonalization for $S_{z}=S_{z}+1$ only in the
$\mathbf{k}_{\text{tot}}$ sector that hosts the ground state.

  5. 5.

The spin polarization $S$ of the ground state is the maximal value of $S_{z}$
that displays the ground state energy found at $S_{z}=\frac{1}{2}$.

Figure S2: The lowest energy per particle of Eq. (15) at filling factor of
$\nu=1$, as obtained by exact diagonalization, as a function of total momentum
sector $|\mathbf{k}_{\text{tot}}|$, and for a range of inverse effective
masses $m_{\text{eff}}^{-1}$. Figures (a) and (b) correspond to $N=7$ and
$N=19$ grid points, respectively.

####  S4.2.2 Simulation results

The results of the first and second steps of the procedure described above are
shown in Fig. S2, where the lowest energy eigenvalue of the Hamiltonian in Eq.
(15) as obtained by exact diagonalization is plotted as a function of the
momentum sector $\mathbf{k}_{\text{tot}}$ (and fixed $S_{z}=\frac{1}{2}$).
Since the Hamiltonian in Eq. (15) has an effective mass as a parameter, we
repeat this simulation for a range of inverse effective masses
$m_{\text{eff}}^{-1}$, increased in steps of size
$m_{0}^{-1}=\frac{V_{k_{0}}N}{2Ak_{0}^{2}}$. For the special value of
$m_{\text{eff}}^{-1}=0$ where the Hamiltonian in Eq. (15) is reduced to Eq.
(2) we find the relation between energy and momentum is approximately linear.
This linearity is due to the single velocity scale ($E/k_{0}\sim e^{2}$) in
this Hamiltonian.

Proceeding with steps (3)-(5), we identify the momentum sector hosting the
lowest energy (for each $m_{\text{eff}}^{-1}$) and repeat the diagonalizing in
this sector for different values of $S_{z}$ to identify its total spin. The
results for the ground state momentum, energy and total spin are shown in Fig
S3. We find the ground state for $-3m_{0}^{-1}\leq m_{\text{eff}}^{-1}\leq
2m_{0}^{-1}$ is fully spin polarized (which implies it has
$\mathbf{k}_{\text{tot}}=0$) for both system sizes. Beyond this regime, i.e.
for $m_{\text{eff}}^{-1}<-3m_{0}^{-1}$ or $2m_{0}^{-1}<m_{\text{eff}}^{-1}$,
we find the ground state has a total spin of $S=\frac{5}{2}$ for both systems
sizes, with a total momentum of $|\mathbf{k}_{\text{tot}}|=k_{0}$ and
$|\mathbf{k}_{\text{tot}}|=\frac{\sqrt{3}}{2}k_{0}$ for $N=7$ and $N=19$,
respectively.

At first sight, in the range of $m_{\text{eff}}^{-1}<-3m_{0}^{-1}$ or
$2m_{0}^{-1}<m_{\text{eff}}^{-1}$, the ground state seems to be both partially
spin-polarized and to have a non-zero momentum. However, we now argue that
this is a finite-size artifact. In the limit of
$m_{\text{eff}}^{-1}\rightarrow\pm\infty$, we expect a completely depolarized
state, i.e. $S_{z}=\frac{1}{2}$. This is because the Coulomb interaction is
completely negligible compared to the kinetic dispersion. In this limit,
Hund’s rule dictates that the favorable configuration among states with the
same kinetic energy is the one that maximizes the total spin. Consider placing
$N=19$ electrons on our momentum grid one at a time while minimizing kinetic
energy. For $m_{\text{eff}}^{-1}\rightarrow\infty$
($m_{\text{eff}}^{-1}\rightarrow-\infty$) we find that the last $5$ electrons
($7$ electrons) have to be placed on $6$ spin-degenerate sites with momenta
$k=\frac{\sqrt{3}}{2}k_{0}$ (see Fig. S4). This leads to a ${12}\choose{5}$
degenerate ground state for either sign of $m_{\text{eff}}^{-1}$. In this
subspace of ground states, the states with maximal total spin have
$S=\frac{5}{2}$, and a $12$-fold degeneracy due to $C_{6}$ rotations and
$S_{z}=\pm\frac{5}{2}$. All the maximal spin configurations occur with
$|\mathbf{k}_{\text{tot}}|=\frac{\sqrt{3}}{2}k_{0}$. By fixing
$\mathbf{k}_{\text{tot}}$ and $S_{z}>0$, it follows from Hund’s rule that the
ground state has this total spin and momentum. A similar argument for the
$N=7$ system leads to $S=\frac{5}{2}$ with $|\mathbf{k}_{\text{tot}}|=k_{0}$,
in perfect agreement with our numerical results. We expect that increasing the
system size will retain this residual spin polarization which gradually
becomes insignificant: $S/N\rightarrow 0$. In the thermodynamic limit, this
artifact would not occur. This is hinted in Fig. 4 by considering the total
spin per particle for different system sizes. We note that if we were to focus
our attention on the $\mathbf{k}_{\text{tot}}=0$ sector where the
thermodynamic limit’s ground state is expected, we find (numerically) its
lowest energy state has $S=\frac{3}{2}$, as expected from the exact same
Hund’s rule argument above. We conclude that the phase outside
$-3m_{0}^{-1}\leq m_{\text{eff}}^{-1}\leq 2m_{0}^{-1}$ for
$m_{0}^{-1}=\frac{V_{k_{0}}N}{2Ak_{0}^{2}}$ is consistent with a spin
unpolarized phase up to finite-size effects.

Figure S3: The ground state energy per particle of Eq. (15) at filling factor
of $\nu=1$, as obtained by exact diagonalization, as a function of inverse
effective mass. Figures (a) and (b) correspond to $N=7$ and $N=19$ grid
points, respectively. The dashed yellow line corresponds to the energy of the
spin-polarized state, which exactly matches the ground state for data points
between the vertical dashed lines. The overlayed text specifies the total spin
and momentum of the obtained ground state. Figure S4: The ground state spin
configuration of the Hamiltonian in Eq. (15) at filling factor $\nu=1$ for
$N=19$ in the limit of (a) $m_{\text{eff}}^{-1}\rightarrow\infty$ and (b)
$m_{\text{eff}}^{-1}\rightarrow-\infty$, as dictated by Hund’s rule. The small
solid circles denote the momentum grid, while the dashed concentric circles
are shells that connect sites with degenerate kinetic energy. The partially
filled shell of momentum $k=\frac{\sqrt{3}}{2}k_{0}$ is occupied in the
configuration that maximizes the total spin. This configuration is found to
have a total spin of $S=\frac{5}{2}$ and a total momentum (indicated by the
black arrow) of $|\mathbf{k}_{\text{tot}}|=\frac{\sqrt{3}}{2}k_{0}$.