# Exact ground state of interacting electrons
in magic angle graphene

Simon Becker<EMAIL_ADDRESS>ETH Zurich, Institute for Mathematical
Research, Rämistrasse 101, 8092 Zurich, Switzerland  ,  Lin Lin
<EMAIL_ADDRESS>Department of Mathematics, University of California,
Berkeley, CA 94720, USA; Applied Mathematics and Computational Research
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA  and
Kevin D. Stubbs<EMAIL_ADDRESS>Department of Mathematics, University of
California, Berkeley, CA 94720, USA

###### Abstract.

One of the most remarkable theoretical findings in magic angle twisted bilayer
graphene (TBG) is the emergence of ferromagnetic Slater determinants as exact
ground states for the interacting Hamiltonian at the chiral limit. This
discovery provides an explanation for the correlated insulating phase which
has been experimentally observed at half filling. This work is the first
mathematical study of interacting models in magic angle graphene systems.
These include not only TBG but also TBG-like systems featuring four flat bands
per valley, and twisted trilayer graphene (TTG) systems with equal twist
angles. We identify symmetries of the Bistritzer-MacDonald Hamiltonian that
are responsible for characterizing the Hartree-Fock ground states as zero
energy many-body ground states. Furthermore, for a general class of
Hamiltonian, we establish criteria that the ferromagnetic Slater determinants
are the unique ground states within the class of uniformly half-filled,
translation invariant Slater determinants. We then demonstrate that these
criteria can be explicitly verified for TBG and TBG-like systems at the chiral
limit, using properties of Jacobi-$\theta$ and Weierstrass-$\wp$ functions.

##  1\. Introduction

Over the past few years, moiré structures, notably “magic angle” twisted
bilayer graphene (TBG), have attracted significant attention in the condensed
matter physics community. This surge of interest began after numerous
experiments [14, 31, 13, 22, 41] which demonstrate that TBG has nearly flat
energy bands and can host intricate phases, such as the correlated insulator
(CI) phases at integer fillings and superconducting (SC) phases at non-integer
fillings. Before these exciting experimental results however, the Bistritzer-
MacDonald (BM) model [8] successfully predicted these nearly flat energy bands
near the magic angle of roughly $1.1^{\circ}$. Yet, the BM model does not
directly include electron-electron repulsion (i.e., it is a non-interacting
model) and in fact, due to symmetry restrictions, the BM model predicts that
twisted bilayer graphene is always in a simple metallic phase. As a result,
the CI and SC phases must arise from electron-electron interactions beyond the
non-interacting BM model.

At the chiral limit, we can define a flat-band interacting (FBI) Hamiltonian
for TBG that only consists of terms describing electron-electron interactions
[10, 15, 32, 39, 38, 16, 7, 20, 30, 18, 29, 21, 17]. One of the most
remarkable theoretical results on this model is that, at half filling, the FBI
Hamiltonian is _frustration-free_ , i.e., the Hamiltonian can be written as a
sum of terms (in this case, these terms do not commute), and there exists a
ground state that minimizes the energy of each term. For the FBI Hamiltonian,
one possible ground state is a single Slater determinant, which is the
simplest type of fermionic many-body wavefunctions. Moreover, it can be shown
that for this ground state, in the thermodynamic limit, a strictly positive
amount of energy is needed to either add or remove an electron. This offers
_one_ possible explanation for the CI phase in TBG at half-filling.

To our knowledge, mathematical studies so far have focussed on understanding
the BM model of twisted bilayer and multilayer graphene at the chiral limit
[2, 1, 3, 5, 4, 6, 40, 37], and systematic derivations of the BM Hamiltonian
from more complex models such as the atomistic tight binding models [36, 11,
12]. In this paper, we provide the _first_ mathematical study of interacting
systems, and specifically the FBI Hamiltonian and its ground states. We note
that our analysis applies to the FBI Hamiltonian defined on a single valley of
the moiré Brillouin zone.

Our analysis is based on the symmetries of the BM Hamiltonian at the chiral
limit. These symmetries are preserved as the number of layers and flat bands
varies, which allows us to study TBG as well as TBG-like systems, including
twisted multilayer graphene structures. For example, the standard chiral model
for TBG has two layers and two flat bands; a system we refer to as TBG-2. By
changing the interlayer potential, it is possible for two layers of graphene
to host four flat bands while still preserving the relevant symmetries [5]; we
refer to this system as TBG-4. We can also consider a twisted trilayer
graphene (TTG) system, where the relative twist between the top and middle and
the relative twist between the middle and bottom are the same, and can host
four flat bands [28]. We refer to this system as the equal angle twisted
trilayer graphene (eTTG-4).

Although the FBI Hamiltonian is an idealized and simplified representation of
interacting electrons, it still contains many intricate details. Our paper
does not discuss the derivation of the FBI Hamiltonian from the BM model, or
the derivation of the BM model from the finer level atomistic models. We
assume readers are familiar with second quantization and standard treatments
of non-interacting periodic systems.

###  1.1. Notation

In this paper, operators and matrices acting on the Fock space are represented
using the hat notation, such as $\hat{f}^{{\dagger}},\hat{f},\hat{H}_{FBI}$.
Operators and matrices that operate in the single particle space (such as
$L^{2}({\mathbb{R}}^{2};{\mathbb{C}}^{2}\times{\mathbb{C}}^{2})$ for TBG) are
indicated without the hat notation, as seen in $H,D$. For a single particle
operaotr $H$, the Bloch-Floquet transformed Hamiltonian is denoted by
$H_{\mathbf{k}}$, where $\mathbf{k}$ is the Brillouin zone index. With some
slight abuse of the hat notation, vectors defined in real space are denoted
without the hat, for example, $u(\mathbf{r})$. Their corresponding Fourier
transforms are indicated with the hat notation, as in $\hat{u}(\mathbf{q})$.
For any given matrix $A$, the operations of entrywise complex conjugation,
transpose, and Hermitian conjugation are represented by
$\overline{A},A^{\top},$ and $A^{{\dagger}}$, respectively.

Due to the symmetries present in twisted graphene systems, the number of flat
bands is always an even number which we will denote by $2M$. We denote the set
of indices of the flat bands by $\mathcal{N}=\\{-M,\cdots,-1,1,\cdots,M\\}$.
After a proper discretization of a single valley of the moiré Brillouin zone
into a discrete set $\mathcal{K}$ with $N_{\mathbf{k}}$ points (see the
definition of $\mathcal{K}$ in (2.12)), we can consider a finite dimensional
Hilbert space $\mathcal{F}$ consisting of $2MN_{\mathbf{k}}$ fermionic modes.
This space is spanned by $2^{2MN_{\mathbf{k}}}$ basis vectors of the form

$\ket{s_{1},\ldots,s_{2MN_{\mathbf{k}}}}=\prod_{n\in\mathcal{N}}\prod_{\mathbf{k}\in\mathcal{K}}(\hat{f}^{{\dagger}}_{n\mathbf{k}})^{s}\ket{\mathrm{vac}},\quad
s_{j}\in\\{0,1\\},$ (1.1)

where $\hat{f}_{n\mathbf{k}}^{{\dagger}},\hat{f}_{n\mathbf{k}}$ are fermionic
creation and annihilation operators satisfying the canonical anticommutation
relation (CAR), and the vacuum state $\ket{\mathrm{vac}}=\ket{0,\ldots,0}$
satisfies $\hat{f}_{n\mathbf{k}}\ket{\mathrm{vac}}=0$ for each
$n\in\mathcal{N}$ and $\mathbf{k}\in\mathcal{K}$.

The FBI Hamiltonian takes the form

$\hat{H}=\sum_{\mathbf{q}^{\prime}}\hat{A}_{\mathbf{q}^{\prime}}^{{\dagger}}\hat{A}_{\mathbf{q}^{\prime}},$
(1.2)

where each $\hat{A}_{\mathbf{q}^{\prime}}$ takes the form

$\hat{A}_{\mathbf{q}^{\prime}}=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\alpha_{\mathbf{k},\mathbf{k}+\mathbf{q}^{\prime}}]_{m,n}\hat{f}^{{\dagger}}_{m\mathbf{k}}\hat{f}_{n(\mathbf{k}+\mathbf{q}^{\prime})}+[\beta_{\mathbf{k},\mathbf{k}+\mathbf{q}^{\prime}}]_{m,n}.$
(1.3)

Here, the sum over $\mathbf{q}^{\prime}\in\mathbb{R}^{2}$ is taken over a
discrete set which will be specified later; $\mathbf{q}^{\prime}$ represents a
momentum difference and so the terms in $\hat{A}_{\mathbf{q}^{\prime}}$
connect states at momentum $\mathbf{k}$ and $\mathbf{k}+\mathbf{q}^{\prime}$.
The coefficient matrices $\alpha_{\mathbf{k},\mathbf{k}+\mathbf{q}^{\prime}}$
and $\beta_{\mathbf{k},\mathbf{k}+\mathbf{q}^{\prime}}$ are matrices of size
$2M\times 2M$ (defined in Eq. 2.22) depend on the eigenstates of the non-
interacting BM Hamiltonian at momenta $\mathbf{k}$ and
$\mathbf{k}+\mathbf{q}^{\prime}$ and inherit a number of symmetries from the
BM Hamiltonian.

The total number operator
$\hat{N}=\sum_{n\in\mathcal{N}}\sum_{\mathbf{k}}\hat{f}^{{\dagger}}_{n\mathbf{k}}\hat{f}_{n\mathbf{k}}$
counts the number of electrons in a state, i.e., if $\hat{N}\ket{\psi}=\nu
N_{\mathbf{k}}\ket{\psi}$, then the number of electrons in $\ket{\psi}$ is
$\nu N_{\mathbf{k}}$. The number operator $\hat{N}$ commutes with $\hat{H}$ so
we can restrict to functions that are simultaneously eigenfunctions of
$\hat{H}$ and $\hat{N}$. The integer filling regime refers to the case when
$\nu$ is an integer ($0\leq\nu\leq 2M$). One particular integer filling is
$\nu=M(=\frac{1}{2}(2M))$, and the corresponding $\ket{\psi}$ is also called a
_half-filled_ state. If we further have
$\sum_{n\in\mathcal{N}}\hat{f}^{{\dagger}}_{n\mathbf{k}}\hat{f}_{n\mathbf{k}}\ket{\psi}=M\ket{\psi}$
for all $\mathbf{k}\in\mathcal{K}$, then $\ket{\psi}$ is a _uniformly half-
filled_ state. The one-body reduced density matrix (1-RDM) is defined as
$[P(\mathbf{k},\mathbf{k}^{\prime})]_{nm}=\braket{\psi}{\hat{f}^{{\dagger}}_{m\mathbf{k}^{\prime}}\hat{f}_{n\mathbf{k}}}{\psi}$.
If
$[P(\mathbf{k},\mathbf{k}^{\prime})]_{nm}=[P(\mathbf{k})]_{nm}\delta_{\mathbf{k},\mathbf{k}^{\prime}}$,
then $\ket{\psi}$ is a _translation invariant_ state.

This paper considers the uniformly half-filled, translation invariant states,
which forms a subspace of $\mathcal{F}$ denoted by $\mathcal{F}_{u}$. The
ground-state energy in $\mathcal{F}_{u}$ is the solution to the following
optimization problem

$E=\min_{\ket{\psi}\in\mathcal{F}_{u}}\frac{\braket{\psi}{\hat{H}}{\psi}}{\braket{\psi}{\psi}},$
(1.4)

and its minimizer (which may not be unique) is called a ground state. For any
$\ket{\psi}\in\mathcal{F}$, by construction
$\braket{\psi}{\hat{H}}{\psi}=\sum_{\mathbf{q}^{\prime}}\braket{\hat{A}_{\mathbf{q}^{\prime}}\psi}{\hat{A}_{\mathbf{q}^{\prime}}\psi}\geq
0$. Therefore $\hat{H}$ is a positive semidefinite (PSD) Hamiltonian, and
$\ket{\psi}\in\mathcal{F}_{u}$ satisfying $\hat{H}\ket{\psi}=0$ is a ground
state.

We further consider a subset of $\mathcal{F}_{u}$

$\mathcal{S}=\left\\{\prod_{i=1}^{M}\prod_{\mathbf{k}\in\mathcal{K}}\hat{b}_{i\mathbf{k}}^{{\dagger}}\ket{\mathrm{vac}}\Big{|}\leavevmode\nobreak\
\hat{b}^{{\dagger}}_{i\mathbf{k}}=\sum_{n\in\mathcal{N}}\hat{f}^{{\dagger}}_{n\mathbf{k}}[\Xi(\mathbf{k})]_{ni},\quad\sum_{n\in\mathcal{N}}[\overline{\Xi(\mathbf{k})}]_{ni}[\Xi(\mathbf{k})]_{nj}=\delta_{ij}\right\\}.$
(1.5)

Each element of $\mathcal{S}$ is called a uniformly half-filled, translation
invariant Slater determinant. The _Hartree-Fock_ theory solves a much simpler
optimization problem

$E_{\text{HF}}=\min_{\begin{subarray}{c}\ket{\psi}\in\mathcal{S}\end{subarray}}\frac{\braket{\psi}{\hat{H}}{\psi}}{\braket{\psi}{\psi}}.$
(1.6)

By definition we have $E\leq E_{\text{HF}}$ and in general $E<E_{\text{HF}}$.

###  1.2. Main results

In this section, we summarize the main findings of this paper, emphasizing the
algebraic structures, while reserving some of the more technical specifics for
later in the paper.

Our first result is that for the FBI Hamiltonian, there exist two states in
$\mathcal{S}$ which are exact ground states. This also implies that the
Hartree-Fock theory is exact in this case.

######  Result 1 (Informal version of Proposition 5.2).

The single Slater determinants

$\ket{\Psi_{+}}=\prod_{n>0,n\in\mathcal{N}}\prod_{\mathbf{k}\in\mathcal{K}}\hat{f}_{n\mathbf{k}}^{\dagger}\ket{\mathrm{vac}},\quad\ket{\Psi_{-}}=\prod_{n<0,n\in\mathcal{N}}\prod_{\mathbf{k}\in\mathcal{K}}\hat{f}_{n\mathbf{k}}^{\dagger},$
(1.7)

satisfy

$\hat{A}_{\mathbf{q}^{\prime}}\ket{\Psi_{\pm}}=0$ (1.8)

for all $\mathbf{q}^{\prime}$. Therefore $E=E_{\operatorname{HF}}=0$, and
$\ket{\Psi_{\pm}}$ are exact ground states of $\hat{H}$.

The proof of Result 1 is a generalization of the results in [7, 10] for TBG.
The mechanism of constructing these exact ground states is to verify that
$\ket{\Psi_{\pm}}$ are ground states of _each_
$\hat{H}_{\mathbf{q}^{\prime}}=\hat{A}^{{\dagger}}_{\mathbf{q}^{\prime}}\hat{A}_{\mathbf{q}^{\prime}}$,
even though the terms $\hat{H}_{\mathbf{q}^{\prime}}$ do not commute with each
other. Therefore the FBI Hamiltonian is an example of a frustration-free
Hamiltonian with nonlocal interactions.

However, there may be other ground states. For instance, any non-vanishing
linear combination $a_{+}\ket{\Psi_{+}}+a_{-}\ket{\Psi_{-}}$ is automatically
a ground state. It is in general a difficult task to identify the entire
ground state manifold. Our main result identifies additional assumptions under
which $\ket{\Psi_{\pm}}$ are the _unique_ states111With slight abuse of
language, we refer to the statement that the ferromagnetic Slater determinant
states are the only elements of $\mathcal{S}$ which are exact ground states of
$\hat{H}$ as “unique states”, even though technically there are two of them.
in $\mathcal{S}$ that are ground states of $\hat{H}$.

######  Result 2 (Main result. Informal version of Theorem 4).

Assume that the set of matrices
$\\{\alpha_{\mathbf{k},\mathbf{k}+\mathbf{q}^{\prime}},\beta_{\mathbf{k},\mathbf{k}+\mathbf{q}^{\prime}}\\}$
satisfy additional non-degeneracy assumptions (Theorem 4), then
$\ket{\Psi_{\pm}}$ in Eq. 1.7 are the unique states in $\mathcal{S}$ that are
exact ground states $\hat{H}$.

Both Result 1 and Result 2 are apply to general FBI Hamiltonians which can be
used to describe the interacting electrons in twisted $N$-layer graphene
systems. The non-degeneracy assumptions can be explicitly verified for
specific systems. We verify these conditions for TBG-2, TBG-4, and eTTG-4
systems using the analytic expression for
$\\{\alpha_{\mathbf{k},\mathbf{k}+\mathbf{q}^{\prime}},\beta_{\mathbf{k}+\mathbf{q}^{\prime}}\\}$.
This leads to the following results.

######  Result 3 (Informal version of Theorems 9, LABEL:, 10, LABEL: and 11).

For TBG-2, TBG-4, and eTTG-4, $\ket{\Psi_{\pm}}$ in Eq. 1.7 are the unique
states in $\mathcal{S}$ that are exact ground states of $\hat{H}$.

While Result 1 and Result 2 provide a characterization of the ground state,
these results alone do not explain why the uniformly half-filled TBG-2 system
(and similar systems) are insulators. In Section 5.4, we compute the ground
state energy of systems with $MN_{\mathbf{k}}+1$ and $MN_{\mathbf{k}}-1$
electrons. Both energies are finite and non-zero in the thermodynamic limit
($N_{\mathbf{k}}\to\infty$). In other words, it costs a finite amount of
energy to add or to remove an electron from the system. Even if the precise
nature of the state can be debated, this shows that the system cannot be in a
simple metallic state.

###  1.3. Discussion and open questions

The states $\ket{\Psi_{\pm}}$ are sometimes referred to as “ferromagnetic
Slater determinants”. This term is inherited from the physics context where
the modes created by $\hat{f}_{n\mathbf{k}}$ with $n=-M,\ldots,-1$ and with
$n=1,\ldots,M$ are vectors of the sublattice symmetry operator $\mathcal{Z}$
with eigenvalues $+1$ and $-1$, respectively. These modes can be interpreted
as “pseudo” up spins and down spins. With this analogy, we can identify
$\ket{\Psi_{+}}$ with $\ket{\uparrow\uparrow\uparrow\cdots}$ and
$\ket{\Psi_{-}}$ with $\ket{\downarrow\downarrow\downarrow\cdots}$, which are
ground states of a ferromagnet. The states $\ket{\Psi_{\pm}}$ also have
nonvanishing Chern numbers, and are also referred to as (integer) quantum hall
states [9].

The mechanism for the uniqueness of the ferromagnetic Slater determinant
ground states in twisted graphene is rather different from previous uniqueness
results in many-body systems (for example, on the Hubbard model on a planar
graph [25, 26, 24]). Unlike these results, the FBI Hamiltonian involves
interactions between all momenta and the ground state is not determined by the
connectivity of the underlying graph structure. Instead, the selection of
ferromagnetic Slater determinant as the ground states seems to be related to
the joint invariant subspaces of the collection of matrices
$\\{\alpha_{\mathbf{k},\mathbf{k}+\mathbf{q}^{\prime}},\beta_{\mathbf{k}+\mathbf{q}^{\prime}}\\}$;
any state whose one-body reduced density matrix does not lie in an invariant
subspace acquires an energy penalty.

These conditions in Result 2 are physically relevant. This paper focuses on
the spinless, valleyless (or in physical terms, spin and valley polarized) FBI
Hamiltonian. When valley degrees of freedom are taken into account, these
assumptions become invalid. This leads to additional ground states in
$\mathcal{S}$, such as the inter-valley coherent (IVC) states [10, 7]. The
generalization of our result to FBI Hamiltonians taking into account both
spins and valleys will be our future work.

On the other hand, the conditions in Result 2 are difficult to verify. They
can be relaxed and explicitly verified for $M=1$ or $2$, as shown in Result 3.
However, for $M>2$, a generally computationally verifiable method for
reformulating these conditions remains unknown.

Even when the conditions are satisfied, Result 2 only establishes uniqueness
among uniformly half-filled, translation-invariant Slater determinants. We
conjecture that at half filling, the uniqueness result should extend to all
single Slater determinants. This requires techniques beyond the scope of this
initial work and will be our future research.

###  1.4. Organization

Our article is structured as follows:

  * •

In Section 2 we define the flat-band interacting model for twisted graphene
systems.

  * •

In Section 3 we review the chiral model of twisted graphene sheets.

  * •

In Section 4, we outline the symmetries that are relevant for our analysis and
fix a gauge of Bloch functions.

  * •

In Section 5, we review basics on Hartree-Fock theory for FBI Hamiltonians and
prove some of their important many-body properties.

  * •

In Section 6, we characterize the ground states of the flat band interacting
model and state our main result, Theorem 4.

  * •

In Section 7, we give the proof of our main result, Theorem 4.

  * •

In Section 8, we verify the assumptions of Theorem 4 for TBG-2, TBG-4, and
eTTG-4.

  * •

Our article contains three technical appendices, Appendices A, LABEL:, B,
LABEL: and C where derive Eq. 7.14, state the proof of Lemma 7.1, and derive a
real space condition for the assumptions of Theorem 4, respectively.

### Acknowledgments

This work was supported by the Simons Targeted Grants in Mathematics and
Physical Sciences on Moiré Materials Magic (K.D.S.) as well as the SNF Grant
PZ00P2 216019 (S.B.). This material is based upon work supported by the U.S.
Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research and Office of Basic Energy Sciences, Scientific Discovery
through Advanced Computing (SciDAC) program under Award Number DE-SC0022198
(L.L.). L.L. is a Simons Investigator in Mathematics. We thank Dumitru
Calugaru, Eslam Khalaf, Patrick Ledwith, and Oskar Vafek for their insightful
discussions.

##  2\. The Flat-Band Interacting Hamiltonian for Twisted Graphene

###  2.1. Notational Setup for Twisted Graphene

Before defining the flat-band interacting model for twisted graphene, we
introduce a general notation for twisted $N$-layer graphene based on the
Bistritzer-MacDonald model. Twisted bilayer and trilayer graphene correspond
to choosing $N=2$ and $N=3$ respectively.

After a proper choice of units, the lattice vectors for the moiré unit cell in
the real space are

$\mathbf{v}_{1}:=-\begin{bmatrix}\frac{2\pi}{\sqrt{3}},\frac{2\pi}{3}\end{bmatrix}^{\top}\qquad\mathbf{v}_{2}:=\begin{bmatrix}\frac{2\pi}{\sqrt{3}},-\frac{2\pi}{3}\end{bmatrix}^{\top}.$
(2.1)

These vectors generate the real space moiré lattice $\Gamma$:

$\Gamma:=\mathbf{v}_{1}{\mathbb{Z}}+\mathbf{v}_{2}{\mathbb{Z}}=\\{a_{1}\mathbf{v}_{1}+a_{2}\mathbf{v}_{2}:a_{1},a_{2}\in{\mathbb{Z}}\\}.$
(2.2)

The moiré unit cell in the real space is denoted by
$\Omega:=\mathbb{R}^{2}/\Gamma$ and can be identified with

$\Omega:=\left\\{t_{1}\mathbf{v}_{1}+t_{2}\mathbf{v}_{2};t_{1},t_{2}\in\left[-\frac{1}{2},\frac{1}{2}\right)\right\\}.$
(2.3)

The dual lattice (or the reciprocal lattice) of $\Gamma$ is denoted by
$\Gamma^{*}$, and is generated by

$\mathbf{g}_{1}:=-\begin{bmatrix}\frac{\sqrt{3}}{2},\frac{3}{2}\end{bmatrix}^{\top}\qquad\mathbf{g}_{2}:=\begin{bmatrix}\frac{\sqrt{3}}{2},-\frac{3}{2}\end{bmatrix}^{\top}.$
(2.4)

Note that $\mathbf{v}_{i}\cdot\mathbf{g}_{j}=2\pi\delta_{ij}$. The unit cell
in the reciprocal space (or the Brillouin zone) is denoted by
$\Omega^{*}:=\mathbb{R}^{2}/\Gamma^{*}$ and can be identified with

$\Omega^{*}=\left\\{t_{1}\mathbf{g}_{1}+t_{2}\mathbf{g}_{2}:t_{1},t_{2}\in\left[-\frac{1}{2},\frac{1}{2}\right)\right\\}.$
(2.5)

Throughout the paper, a vector in the reciprocal lattice $\Gamma^{*}$ is often
denoted by $\mathbf{G}$, while a vector in the Brillouin zone $\Omega^{*}$ is
often denoted by $\mathbf{k}$ or $\mathbf{q}$. Note that a generic vector
$\mathbf{q}^{\prime}\in{\mathbb{R}}^{2}$ can always be uniquely decomposed as
$\mathbf{q}^{\prime}=\mathbf{q}+\mathbf{G}$ for some $\mathbf{q}\in\Omega^{*}$
and $\mathbf{G}\in\Gamma^{*}$.

In addition to the generating vectors $\mathbf{g}_{1}$ and $\mathbf{g}_{2}$,
we also identify three special momentum vectors
$\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}$ which are related to each other
by a $\frac{2\pi}{3}$-counterclockwise rotation $R_{3}$:

$R_{3}:=\frac{1}{2}\begin{bmatrix}-1&-\sqrt{3}\\\ \sqrt{3}&-1\end{bmatrix}.$

In particular,

$\mathbf{q}_{1}:=[0,1]^{\top}\quad\mathbf{q}_{2}:=R_{3}\mathbf{q}_{1}=[-\sqrt{3}/2,-1/2]^{\top}\quad\mathbf{q}_{3}:=R_{3}\mathbf{q}_{2}=[\sqrt{3}/2,-1/2]^{\top}$
(2.6)

Notice that $\mathbf{g}_{1}=\mathbf{q}_{2}-\mathbf{q}_{1}$ and
$\mathbf{g}_{2}=\mathbf{q}_{3}-\mathbf{q}_{1}$.

Each moiré unit cell consists of many atoms located on one of the two
sublattices, denoted by $\\{A=1,B=-1\\}$. In the BM model for twisted bilayer
and multilayer graphene, all atomistic details in each layer are smeared out
except the sublattice information. The BM Hamiltonian $H$ (called a single
particle Hamiltonian) acts on a dense subset of the function space
$L^{2}({\mathbb{R}}^{2};{\mathbb{C}}^{2}\times{\mathbb{C}}^{N})$, where $N$ is
the number of layers. A function in
$L^{2}({\mathbb{R}}^{2};{\mathbb{C}}^{2}\times{\mathbb{C}}^{N})$ can be
written as $\psi(\mathbf{r};\sigma,\ell)$, where
$\mathbf{r}=[x_{1},x_{2}]^{\top}\in{\mathbb{R}}^{2}$ is the real space
coordinate, $\sigma\in\\{\pm 1\\}$ is the sublattice index, and
$\ell\in\\{1,\ldots,N\\}$ is the layer index. The BM Hamiltonian satisfies
certain translation symmetry with respect to the moiré lattice $\Gamma$.
Therefore an eigenfunction of $H$ can be labeled as
$\psi_{n\mathbf{k}}(\mathbf{r};\sigma,\ell)$. Here $n$ is called the band
index, $\mathbf{k}=(k_{1},k_{2})\in\Omega^{*}$ is the Bloch vector. For
example, a wavefunction for the Hamiltonian of TBG can be written as

$\begin{bmatrix}\psi_{n\mathbf{k}}(\mathbf{r};A,1)&\psi_{n\mathbf{k}}(\mathbf{r};A,2)&\psi_{n\mathbf{k}}(\mathbf{r};B,1)&\psi_{n\mathbf{k}}(\mathbf{r};B,2)\end{bmatrix}^{\top}\\\
$ (2.7)

and for TTG:

$\begin{bmatrix}\psi_{n\mathbf{k}}(\mathbf{r};A,1)&\psi_{n\mathbf{k}}(\mathbf{r};A,2)&\psi_{n\mathbf{k}}(\mathbf{r};A,3)&\psi_{n\mathbf{k}}(\mathbf{r};B,1)&\psi_{n\mathbf{k}}(\mathbf{r};B,2)&\psi_{n\mathbf{k}}(\mathbf{r};B,3)\end{bmatrix}^{\top}.$
(2.8)

###  2.2. The Flat-Band Interacting (FBI) Hamiltonian

While we adopt first quantization in expressing the BM Hamiltonian (see
Section 3), the FBI Hamiltonian is most conveniently expressed in second
quantization. The FBI Hamiltonian is defined in terms of the flat-band
eigenfunctions of the single particle Hamiltonian $H$; therefore, our first
step will be to construct such eigenfunctions. We make the following mild
assumptions on $H$:

######  Assumption 1.

We assume that $H$ is periodic with respect to the lattice $\Gamma$ and has
Bloch eigenpairs
$\\{(\varepsilon_{n\mathbf{k}},\psi_{n\mathbf{k}}):n\in{\mathbb{Z}}\setminus\\{0\\},\mathbf{k}\in\Omega^{*}\\}$
so that the eigenvalues $\varepsilon_{n\mathbf{k}}$ are ordered in non-
decreasing order and

$\begin{split}H\psi_{n\mathbf{k}}(\mathbf{r})&=\varepsilon_{n\mathbf{k}}\psi_{n\mathbf{k}}(\mathbf{r}),\\\
\psi_{n\mathbf{k}}(\mathbf{r}+\mathbf{a})&=e^{i\mathbf{a}\cdot\mathbf{k}}\psi_{n\mathbf{k}}(\mathbf{r})\qquad\forall\mathbf{a}\in\Gamma,\\\
\cdots\leq\varepsilon_{-1\mathbf{k}}\leq&\leavevmode\nobreak\
0\leq\varepsilon_{1\mathbf{k}}\leq\cdots\qquad\forall\mathbf{k}\in\Omega^{*}.\end{split}$
(2.9)

We furthermore define the Bloch-Floquet transformed Hamiltonian
$H_{\mathbf{k}}:=e^{i\mathbf{k}\cdot\mathbf{r}}He^{-i\mathbf{k}\cdot\mathbf{r}}$
for all $\mathbf{k}\in\Omega^{*}$ and assume that

  1. (1)

The resolvent operator $(H_{\mathbf{k}}-i)^{-1}$ is compact.

  2. (2)

If we identify $\mathbf{k}\in\Omega^{*}$ with the complex number
$k_{x}+ik_{y}\in{\mathbb{C}}$, the operator valued complex map
$\mathbf{k}\mapsto(H_{\mathbf{k}}-i)^{-1}$ depends analytically on
$\mathbf{k}$.

  3. (3)

There exists a non-empty set of flat bands $\mathcal{N}$

$\mathcal{N}:=\\{n\in{\mathbb{Z}}:\forall\mathbf{k}\in\Omega^{*},\,\varepsilon_{n\mathbf{k}}=0\\}.$
(2.10)

Figure 1. Magic angles with multiplicities in the chiral limit of 7 equally
twisted layers of graphene with multiplicities and $U_{+}=U_{0}$, see (3.14).
Band configuration for largest magic angle $\alpha=0.6922$ showing the six
lowest bands with positive energy including the four flat ones showing three
flat band crossings at high symmetry points.

In some systems, such as twisted trilayer graphene (TTG), there can be energy
bands which intersect with the flat bands but are not flat themselves. This
can happen in all $N$-layer systems with $N\geq 3$ (see Figure 1 for an
example). Let us define such the set of momenta where such crossings occur as
follows:

$\mathcal{K}_{\rm
crossing}=\\{\mathbf{k}\in\Omega^{*}:\varepsilon_{n\mathbf{k}}=0,\text{ for
some }n\not\in\mathcal{N}\\}.$ (2.11)

Due to the analyticity assumption on $H_{\mathbf{k}}$, standard results on
analytic families of operators [19, VII.5, Theo 1.10] implies that the set of
such crossings must have an empty interior. Furthermore, Rellich’s theorem
implies that the spectral projectors of $H_{\mathbf{k}}$ may be chosen
analytically throughout the entire Brillioun zone [19, VII.5, Theo 3.9].
Therefore, we will assume without loss of generality that the flat-band
eigenfunctions
$\\{\psi_{n\mathbf{k}}:n\in\mathcal{N},\mathbf{k}\in\mathcal{K}\\}$ have been
chosen so that the spectral projection onto these eigenfunctions is analytic
in $\mathbf{k}$.

In order to write down a many-body wavefunction with a finite number of
particles, we discretize the Brillouin zone using a uniform grid $\mathcal{K}$
of size $(n_{k_{x}},n_{k_{y}})$

$\mathcal{K}:=\left\\{\frac{i}{n_{k_{x}}}\mathbf{g}_{1}+\frac{j}{n_{k_{y}}}\mathbf{g}_{2}:i\in\\{0,1,\cdots,n_{k_{x}}-1\\},j\in\\{0,1,\cdots,n_{k_{y}}-1\\}\right\\}\subseteq\Omega^{*}.$
(2.12)

Here $\mathbf{g}_{1}$, $\mathbf{g}_{2}$ are a pair of generating vectors for
the moiré dual lattice $\Gamma^{*}$, see (2.4), and we define
$N_{\mathbf{k}}:=\\#|\mathcal{K}|=n_{k_{x}}n_{k_{y}}$.

The set $\mathcal{K}$ is known as a Monkhorst-Pack grid, and the
_thermodynamic limit_ (TDL) can be reached by taking
$n_{k_{x}},n_{k_{y}}\to\infty$. We also make the following technical
assumption on the grid $\mathcal{K}$ in relation to the Hamiltonian $H$:

######  Assumption 2.

For each $\mathbf{k}\in\mathcal{K}$, let $\Pi(\mathbf{k})$ be the orthogonal
projector onto the flat-band eigenfunctions at $\mathbf{k}$. We assume that
the grid $\mathcal{K}$ has been chosen so that for all pairs of momenta
$\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}$ there exists a sequence of
momenta $\\{\mathbf{k}_{i}\\}_{i=1}^{L}$ so that $\mathbf{k}_{1}=\mathbf{k}$,
$\mathbf{k}_{L}=\mathbf{k}^{\prime}$ and
$\|\Pi(\mathbf{k}_{i+1})-\Pi(\mathbf{k}_{i})\|<1$ for all
$i\in\\{1,\ldots,L\\}$.

######  Remark 1.

Since $H_{\mathbf{k}}$ is assumed to be analytic (Assumption 1), the resolvent
mapping $\mathbf{k}\mapsto(H_{\mathbf{k}}-i)^{-1}$ is Lipschitz in
$\mathbf{k}$. Therefore, for Assumption 2 to be true, it is sufficient for
$n_{k_{x}}\gtrsim L$ and $n_{k_{y}}\gtrsim L$ where $L$ is the Lipschitz
constant for $(H_{\mathbf{k}}-i)^{-1}$.

Having fixed the grid $\mathcal{K}$ and the flat-band eigenfunctions, we may
now define the FBI Hamiltonian. For each flat-band
$\psi_{n\mathbf{k}}(\mathbf{r}),\mathbf{k}\in\mathcal{K}$, we first define the
band creation and annihilation operators, $\hat{f}_{n\mathbf{k}}^{\dagger}$
and $\hat{f}_{n\mathbf{k}}$, which create or annihilate a particle in state
$\psi_{n\mathbf{k}}(\mathbf{r})$ respectively, which satisfy the canonical
anti-commutation relation (CAR):

$\\{\hat{f}^{{\dagger}}_{n\mathbf{k}},\hat{f}_{n^{\prime}\mathbf{k}^{\prime}}\\}=\delta_{nn^{\prime}}\delta_{\mathbf{k}\mathbf{k}^{\prime}},\quad\\{\hat{f}^{{\dagger}}_{n\mathbf{k}},\hat{f}^{{\dagger}}_{n^{\prime}\mathbf{k}^{\prime}}\\}=\\{\hat{f}_{n\mathbf{k}},\hat{f}_{n^{\prime}\mathbf{k}^{\prime}}\\}=0.$
(2.13)

The definition of the creation and annihilation operators can be periodically
extended outside the Brillouin zone as

$\hat{f}^{{\dagger}}_{n(\mathbf{k}+\mathbf{G})}=\hat{f}^{{\dagger}}_{n\mathbf{k}},\quad\hat{f}_{n(\mathbf{k}+\mathbf{G})}=\hat{f}_{n\mathbf{k}},\quad\mathbf{G}\in\Gamma^{*}.$
(2.14)

Next, we define the flat-band periodic Bloch functions

$u_{n\mathbf{k}}(\mathbf{r}):=e^{-i\mathbf{k}\cdot\mathbf{r}}\psi_{n\mathbf{k}}(\mathbf{r})$
(2.15)

which are normalized in the unit cell, i.e.,
$\int_{\Omega}|u_{n\mathbf{k}}(\mathbf{r})|^{2}\,\mathrm{d}\mathbf{r}=1$. Let
$\hat{u}_{n\mathbf{k}}(\mathbf{G})$ denote the Fourier coefficients of
$u_{n\mathbf{k}}(\mathbf{r};\sigma,j)$

$\hat{u}_{n\mathbf{k}}(\mathbf{G};\sigma,j):=\int_{\Omega}e^{-i\mathbf{G}\cdot\mathbf{r}}u_{n\mathbf{k}}(\mathbf{r};\sigma,j)d\mathbf{r}.$
(2.16)

and define the form factor

$[\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})]_{mn}:=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}\in\Gamma^{*}}\sum_{\sigma,j}\overline{\hat{u}_{m\mathbf{k}}(\mathbf{G}^{\prime};\sigma,j)}\hat{u}_{n(\mathbf{k}+\mathbf{q}+\mathbf{G})}(\mathbf{G}^{\prime};\sigma,j).$
(2.17)

Given $\mathbf{k},\mathbf{q}\in\Omega^{*},\mathbf{G}\in\Gamma^{*}$,
$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ is a
$\\#|\mathcal{N}|\times\\#|\mathcal{N}|$ momentum dependent matrix. Note that
$\\{\mathbf{q}+\mathbf{G}\leavevmode\nobreak\ |\leavevmode\nobreak\
\mathbf{q}\in\Omega^{*},\mathbf{G}\in\Gamma^{*}\\}={\mathbb{R}}^{2}$; this
notation facilitates later discussions.

The following identities which will be useful at various points in our proof.

######  Lemma 2.1.

The form factor matrix $\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ satisfies
the following identities for all $\mathbf{k},\mathbf{q}\in\mathcal{K}$ and all
$\mathbf{G},\mathbf{G}^{\prime}\in\Gamma^{*}$

$\displaystyle\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})^{\dagger}$
$\displaystyle=\Lambda_{\mathbf{k}+\mathbf{q}}(-\mathbf{q}-\mathbf{G}),$
(2.18)
$\displaystyle\Lambda_{\mathbf{k}+\mathbf{G}^{\prime}}(\mathbf{q}+\mathbf{G})$
$\displaystyle=\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G}).$ (2.19)

###### Proof.

This follows from the following simple calculation

$\begin{split}\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})^{\dagger}&=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}\in\Gamma^{*}}\sum_{\sigma,j}\hat{u}_{n\mathbf{k}}(\mathbf{G}^{\prime};\sigma,j)\overline{\hat{u}_{m(\mathbf{k}+\mathbf{q}+\mathbf{G})}(\mathbf{G}^{\prime};\sigma,j)}\\\
&=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}\in\Gamma^{*}}\sum_{\sigma,j}\overline{\hat{u}_{m(\mathbf{k}+\mathbf{q}+\mathbf{G})}(\mathbf{G}^{\prime};\sigma,j)}\hat{u}_{n\mathbf{k}}(\mathbf{G}^{\prime};\sigma,j)\\\
&=\Lambda_{\mathbf{k}+\mathbf{q}}(-\mathbf{q}-\mathbf{G}).\end{split}$ (2.20)

As for the second identity, one easily checks that
$\hat{u}_{n(\mathbf{k}+\mathbf{G}^{\prime})}(\mathbf{G})=\hat{u}_{n\mathbf{k}}(\mathbf{G}+\mathbf{G}^{\prime})$
and so

$\begin{split}\Lambda_{\mathbf{k}+\mathbf{G}^{\prime}}(\mathbf{q}+\mathbf{G})&=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime\prime}\in\Gamma^{*}}\sum_{\sigma,j}\overline{\hat{u}_{m(\mathbf{k}+\mathbf{G}^{\prime})}(\mathbf{G}^{\prime\prime}+\mathbf{G}^{\prime};\sigma,j)}\hat{u}_{n(\mathbf{k}+\mathbf{q}+\mathbf{G}+\mathbf{G}^{\prime})}(\mathbf{G}^{\prime\prime};\sigma,j)\\\
&=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime\prime}\in\Gamma^{*}}\sum_{\sigma,j}\overline{\hat{u}_{m\mathbf{k}}(\mathbf{G}^{\prime\prime}+\mathbf{G}^{\prime};\sigma,j)}\hat{u}_{n(\mathbf{k}+\mathbf{q}+\mathbf{G})}(\mathbf{G}^{\prime\prime}+\mathbf{G}^{\prime};\sigma,j)\\\
&=\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G}).\end{split}$ (2.21)

∎

We can now define the flat-band interacting (FBI) Hamiltonian for the twisted
$N$-layer graphene [10, 7]. When we compare the definition of $\hat{H}_{FBI}$
with that in Eqs. 1.2 and 1.3, we need to show that
$\hat{\rho}^{{\dagger}}(\mathbf{q}^{\prime})=\hat{\rho}(-\mathbf{q}^{\prime})$,
and that $\hat{H}_{FBI}$ is a positive semi-definite Hamiltonian. These
relations will be verified in Section 5.2.

######  Definition 2.2.

Let $\hat{f}_{m\mathbf{k}}^{\dagger},\hat{f}_{m\mathbf{k}}$ denote the flat-
band creation and annihilation operators, and let
$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ be defined as in Eq. 2.17. The
flat-band interacting (FBI) Hamiltonian is:

$\begin{split}\hat{H}_{FBI}&:=\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\widehat{\rho}(\mathbf{q}^{\prime})\widehat{\rho}(-\mathbf{q}^{\prime})\\\
\widehat{\rho}(\mathbf{q}^{\prime})&:=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}\left(\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n(\mathbf{k}+\mathbf{q}^{\prime})}-\frac{1}{2}\delta_{mn}\sum_{\mathbf{G}}\delta_{\mathbf{q}^{\prime},\mathbf{G}}\right).\end{split}$
(2.22)

Here, $\hat{V}(\mathbf{q}^{\prime})=\hat{V}(|\mathbf{q}^{\prime}|)$ is the
Fourier transform of a radially symmetric electron-electron potential which
satisfies $\hat{V}(\mathbf{q}^{\prime})>0$ for all
$\mathbf{q}\in\mathbb{R}^{2}$.

For instance, for TBG, the double gate-screened Coulomb interaction in Fourier
space reads (Note that
$\lim_{\mathbf{q}^{\prime}\to\bm{0}}\hat{V}(\mathbf{q}^{\prime})=\frac{\pi
d}{\varepsilon}$ is well defined)

$\hat{V}(\mathbf{q}^{\prime})=\frac{2\pi}{\varepsilon}\frac{\tanh(\left\lvert\mathbf{q}^{\prime}\right\rvert
d/2)}{\left\lvert\mathbf{q}^{\prime}\right\rvert}.$ (2.23)

Here $\varepsilon,d>0$ parametrize the strength and length of the screened
Coulomb interaction, respectively (see e.g., [7, Appendix C]).

##  3\. A Review of Twisted Bilayer Graphene and equal Twist Angle Trilayer
Graphene

To make the assumptions of Result 2 more concrete, we review the properties of
the TBG and eTTG [34, 2, 28, 3, 36, 11] in the Bistritzer-MacDonald
Hamiltonian at the chiral limit. The single particle Hamiltonian for both
chiral TBG and chiral eTTG depend on a parameter $\alpha$, that is inversely
proportional to the twisting angle, and take the form of a matrix-valued
differential operator

$H(\alpha)=\begin{bmatrix}0&D(\alpha)^{\dagger}\\\ D(\alpha)&0\end{bmatrix},$
(3.1)

For TBG, $D(\alpha)$ acts on $H^{1}(\mathbb{R}^{2};{\mathbb{C}}^{2})$

$D_{\rm TBG}(\alpha)=\begin{bmatrix}D_{x_{1}}+iD_{x_{2}}&\alpha
U_{+}(\mathbf{r})\\\ \alpha
U_{-}(\mathbf{r})&D_{x_{1}}+iD_{x_{2}}\end{bmatrix},$ (3.2)

and for eTTG, $D(\alpha)$ acts on $H^{1}(\mathbb{R}^{2};{\mathbb{C}}^{3})$

$D_{\rm eTTG}(\alpha)=\begin{bmatrix}D_{x_{1}}+iD_{x_{2}}&\alpha
U_{+}(\mathbf{r})&0\\\ \alpha U_{-}(\mathbf{r})&D_{x_{1}}+iD_{x_{2}}&\alpha
U_{+}(\mathbf{r})\\\ 0&\alpha
U_{-}(\mathbf{r})&D_{x_{1}}+iD_{x_{2}}\end{bmatrix}.$ (3.3)

Similarly, single particle Hamiltonians for $N$-layer systems can be defined
in terms of operators $D(\alpha)$ acting on
$H^{1}(\mathbb{R}^{2};{\mathbb{C}}^{N})$ (see, for example, [40]).

We remark that the specific form of the single particle Hamiltonian will not
be used in our proofs of Results 1, LABEL:, 3, LABEL: and 2, however we will
make use of two important properties which hold for Bistritzer-MacDonald-type
Hamiltonians at the chiral limit: (1) the existence of a magic angle with
exactly flat bands, and (2) the single particle Hamiltonian commutes with a
number of symmetry operations (see Section 4).

For small twist angles, due to the lattice structure of graphene, the
tunneling potentials $U_{\pm}(\mathbf{r})$ satisfy the following symmetry
properties independent of the number of layers [8, 34]:

$\displaystyle U_{\pm}(\mathbf{r}+\mathbf{a})$
$\displaystyle=\overline{\omega}^{\pm(a_{1}+a_{2})}U_{\pm}(\mathbf{r}),$ (3.4)
$\displaystyle U_{\pm}(R_{3}\mathbf{r})$ $\displaystyle=\omega
U_{\pm}(\mathbf{r}),$ (3.5) $\displaystyle U_{\pm}(x_{1},-x_{2})$
$\displaystyle=\overline{U_{\pm}(x_{1},x_{2})},$ (3.6)

where $\mathbf{r}=[x_{1},x_{2}]^{\top}$,
$\mathbf{a}=[a_{1},a_{2}]^{\top}\in\Gamma$, $\omega:=e^{2\pi i/3}$, and
$R_{3}$ is the $\frac{2\pi}{3}$-counterclockwise rotation. Potentials
$U_{\pm}$ satisfying symmetries Eqs. 3.4, 3.5 and 3.6 are of the general form

$U_{\pm}(\mathbf{r})=\sum_{n,m\in{\mathbb{Z}}}c_{nm}e^{\pm
i(m\mathbf{g}_{1}-n\mathbf{q}_{2}+\mathbf{q}_{1})\cdot\mathbf{r}}$ (3.7)

where we recall the definitions of the dual lattice vectors
$\mathbf{g}_{1},\mathbf{g}_{2},$ (Eq. 2.4) and $\mathbf{q}_{1}$ (Eq. 2.6). The
coefficients $c_{nm}$ satisfy, see [3, Prop.2.1],

$\begin{split}c_{nm}&=\omega c_{(m-n-1)(-n)}=\omega^{2}c_{(-m)(n-m+1)}\text{
and }\\\ \overline{c_{nm}}&=c_{(-m)(-n)}=\omega
c_{(-n+m-1)m}=\omega^{2}c_{n(n-m+1)}\text{ for all
}n,m\in{\mathbb{Z}}^{2}.\end{split}$

It is important to note that due to Eq. 3.4, both TBG and eTTG as defined
above are not periodic with respect to $\Gamma$, however they can both be made
into periodic Hamiltonian by performing a unitary transformation (see Section
3.3). While this non-periodic formulation is more convenient for the analysis
of the magic angles; we will use the periodic formulation in our definition of
the FBI Hamiltonians (Eq. 2.22).

###  3.1. Magic Angles in TBG

The TBG Hamiltonian (Eq. 3.1 with Eq. 3.2) commutes with the following
symmetries

$\mathcal{T}_{\mathbf{a}}\mathbf{u}(\mathbf{r})=\begin{bmatrix}\omega^{a_{1}+a_{2}}&&&\\\
&1&&\\\ &&\omega^{a_{1}+a_{2}}&\\\
&&&1\end{bmatrix}\mathbf{u}(\mathbf{r}+\mathbf{a}),\qquad\mathcal{R}\mathbf{u}(\mathbf{r})=\begin{bmatrix}1&&&\\\
&1&&\\\ &&\overline{\omega}&\\\
&&&\overline{\omega}\end{bmatrix}\mathbf{u}(R_{3}\mathbf{r}).$ (3.8)

Therefore, we can define

$L^{2}_{\ell,p}:=\\{\mathbf{u}\in
L^{2}_{\text{loc}}({\mathbb{R}}^{2};{\mathbb{C}}^{4});\mathcal{T}_{\mathbf{a}}\mathbf{u}=\omega^{\ell(a_{1}+a_{2})}\mathbf{u},\text{
for }\mathbf{a}\in\Gamma\text{ and
}\mathcal{R}\mathbf{u}=\overline{\omega}^{p}\mathbf{u}\\}$ (3.9)

where $\ell,p\in{\mathbb{Z}}_{3}.$ We may also consider symmetries (3.8) just
acting on ${\mathbb{C}}^{2}$-valued spinors by considering merely the first
two components.

The $\mathcal{T}_{\mathbf{a}}$-periodicity of the TBG Hamiltonian in Eq. 3.1
allows us to define the Bloch-Floquet transformed TBG Hamiltonian

$H_{\mathbf{k}}(\alpha)=\begin{bmatrix}0&D_{\rm
TBG}(\alpha)^{\dagger}+(k_{1}-ik_{2})\operatorname{id}_{\mathbb{C}^{2}}\\\
D_{\rm
TBG}(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{\mathbb{C}^{2}}&0\end{bmatrix}$
(3.10)

acting on $L^{2}_{\ell}:=\bigoplus_{p\in{\mathbb{Z}}_{3}}L^{2}_{\ell,p}$.

The TBG Hamiltonian additionally has a layer symmetry $\mathcal{L}$ which acts
as follows

$\mathcal{L}\mathbf{v}(\mathbf{r}):=\begin{bmatrix}&&&-1\\\ &&1&\\\ &-1&&\\\
1&&&\\\ \end{bmatrix}\mathbf{v}(-\mathbf{r}).$ (3.11)

One furthermore can check that $\mathcal{L}:L^{2}_{\ell,p}\to
L^{2}_{-\ell+1,p}$. In particular, for $\ell=2$ the map $\mathcal{L}$ leaves
spaces $L^{2}_{2,p}$ invariant. Thus, we shall without loss of generality,
consider $H_{\mathbf{k}}(\alpha)$ on $L^{2}_{\ell=2}$ in the sequel, since the
family of $H_{\mathbf{k}}(\alpha)$ on $L^{2}_{\ell}$ are all equivalent.

One then defines a magic angle of the chiral Hamiltonian as a parameter
$\alpha\in{\mathbb{C}}$ such that

$0\in\bigcap_{\mathbf{k}\in\mathbb{R}^{2}}\operatorname{Spec}(H_{\mathbf{k}}(\alpha)),$
(3.12)

where $H_{\mathbf{k}}(\alpha)$ is defined as in (3.10).

For the Hamiltonian (3.10) the first $2$ components correspond to lattice
sites $A$ and the remaining ones to lattice sites $B$. Since
$D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{\mathbb{C}^{2}}$ is a Fredholm
operator of index $0$, at a magic angle $\alpha$, the null space of the
Hamiltonian decomposes into

$\ker(H_{\mathbf{k}}(\alpha))=\ker(D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{\mathbb{C}^{2}})\oplus\ker(D(\alpha)^{*}+(k_{1}-ik_{2})\operatorname{id}_{\mathbb{C}^{2}})$
(3.13)

where $\ker(D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{\mathbb{C}^{2}})$ and
$\ker(D(\alpha)^{*}+(k_{1}-ik_{2})\operatorname{id}_{\mathbb{C}^{2}})$ are of
equal dimension. Elements of
$\ker(D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{\mathbb{C}^{2}})$ therefore
correspond to states at zero energy that are $A$-lattice polarized. This leads
us to make the following definition

######  Definition 3.1 (Multiplicity).

A magic angle $\alpha$ is $n$-fold degenerate if the number of zero energy
flat bands of the Hamiltonian is $2n$-fold degenerate. A $1$-fold degenerate
magic angle is also called _simple_.

In the case of TBG and eTTG we indicate the number $2M$ of flat bands at a
magic angle by appending a number $2M$, e.g. TBG$-2$ for simple magic angles
and TBG$-4$ for $2$-fold degenerates ones.

The simplest one-parameter family of potentials satisfying (3.7) is, for
$\varphi\in\mathbb{R}/\mathbb{Z}$, given by

$U_{\varphi}(\mathbf{r})=\cos(2\pi\varphi)\sum_{i=0}^{2}\omega^{i}e^{-i\mathbf{q}_{i}\cdot\mathbf{r}}+\sin(2\pi\varphi)\sum_{i=0}^{2}\omega^{i}e^{2i\mathbf{q}_{i}\cdot\mathbf{r}}.$

In twisted bilayer graphene, higher Fourier modes in the tunnelling potential
are necessary to describe lattice relaxation effects, see e.g. [36, Remark
$2.4$] and references therein. We then consider special cases

$U_{0}(\mathbf{r})=\sum_{i=0}^{2}\omega^{i}e^{-i\mathbf{q}_{i}\cdot\mathbf{r}}\text{
and
}U_{7/8}(\mathbf{r})=\frac{1}{\sqrt{2}}\Bigg{(}U_{0}(\mathbf{r})-\sum_{i=0}^{2}\omega^{i}e^{2i\mathbf{q}_{i}\cdot\mathbf{r}}\Bigg{)}.$
(3.14)

While both $U_{0}$ and $U_{7/8}$ satisfy the above symmetries (3.4), (3.5),
and (3.6), they give rise to a different number of flat bands for _magic_
$\alpha\in\mathbb{R}$. Numerical experiments suggest, see Figure 2, that in
case of $U=U_{0}$ and $\alpha\in{\mathbb{R}}$ magic, precisely two bands of
the chiral Hamiltonian become flat, in the case of $U=U_{7/8}$ and
$\alpha\in\mathbb{R}$ the Hamiltonian exhibits four flat bands at _magic
angles_. Thus, the chiral model with the two potentials $U=U_{0}$ and
$U=U_{7/8}$ serve as models for TBG$-2$ and TBG$-4$. While simple and two-fold
degenerate magic angles are the only types of magic angles that appear for
generic choices of tunnelling potentials in chiral limit TBG, see [5, Theo.
$3$], they exhibit an almost equidistant spacing for the potentials $U_{0}$
and $U_{7/8}$, as shown in Table 1.

Figure 2.  Magic angles $\alpha$ derived from potentials $U=U_{0}$ (left) and
$U=U_{7/8}$ (right). The dimension of
$\operatorname{ker}(D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{\mathbb{C}^{2}})$
is indicated by the numbers in the figure. No number indicates a one-
dimensional subspace. The number of flat bands of the Hamiltonian is twice the
dimension of
$\operatorname{ker}(D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{\mathbb{C}^{2}}).$

$k$ |  | $\alpha_{k}$ |  | $\alpha_{k}-\alpha_{k-1}$
---|---|---|---|---
1 |  | 0.58566 |  |
2 |  | 2.22118 |  | 1.6355
3 |  | 3.75140 |  | 1.5302
4 |  | 5.27649 |  | 1.5251
5 |  | 6.79478 |  | 1.5183
6 |  | 8.31299 |  | 1.5182
7 |  | 9.82906 |  | 1.5161
8 |  | 11.34534 |  | 1.5163
9 |  | 12.86061 |  | 1.5153
10 |  | 14.37607 |  | 1.5155
11 |  | 15.89096 |  | 1.5149

$k$ |  | $\alpha_{k}$ |  | $\alpha_{k}-\alpha_{k-1}$
---|---|---|---|---
1 |  | 0.853799 |  |
2 |  | 2.691433 |  | 1.8376
3 |  | 4.507960 |  | 1.8165
4 |  | 6.332311 |  | 1.8244
5 |  | 8.157130 |  | 1.8248
6 |  | 9.983510 |  | 1.8264
7 |  | 11.809376 |  | 1.8259
8 |  | 13.635446 |  | 1.8261
9 |  | 15.460894 |  | 1.8255
10 |  | 17.286231 |  | 1.8253
11 |  | 19.111041 |  | 1.8248

Table 1. First $11$ real magic angles, rounded to 6 digits, for $U=U_{0}$
(left) and $U=U_{7/8}$ (right).

###  3.2. Magic Angles in eTTG

When studying twisted trilayer graphene (TTG), the case of equal twisting
angles (eTTG) is of particular interest.

In a recent paper by Popov and Tarnopolsky [28], the authors demonstrated that
the collection of magic parameters $\mathcal{A}_{\mathrm{eTTG}}$ exhibits a
particularly interesting relation, which connects the set of magic angles of
chiral limit TBG to the magic angles of equal twisting angle TTG (eTTG)

$\sqrt{2}\mathcal{A}_{\mathrm{TBG}-M}=\mathcal{A}_{\mathrm{eTTG}-2M},$ (3.15)

see also the top right figure in Figure 3. Here
$\mathcal{A}_{\mathrm{eTTG}-2M}$ is the set of magic parameters for eTTG-2M.
Moreover, the multiplicity on the right is at least twice the one on the left.
The argument provided in [28] constructs protected flat bands for eTTG-2M from
flat bands of TBG-M. This construction works for eTTG-2M, however it does not
extend to more than three layers nor does not seem to extend beyond the equal
twisting angle case, as well.

At the center of this connection, is the map with obvious indices
corresponding to $D(\alpha)=D_{\text{TBG}}(\alpha)$ as in (3.2) and
$D(\alpha)=D_{\text{eTTG}}(\alpha)$ as in (3.3)

$\begin{split}\times:\ker_{L^{2}_{k,j}}(D_{\text{TBG}}(\alpha_{0}))&\oplus\ker_{L^{2}_{k^{\prime},j^{\prime}}}(D_{\text{TBG}}(\alpha_{0}))\to\ker_{L^{2}_{k+k^{\prime},j+j^{\prime}}}(D_{\text{eTTG}}(\sqrt{2}\alpha_{0}))\\\
(v,w)&\mapsto(v_{1}w_{1},2^{-1/2}(v_{1}w_{2}+v_{2}w_{1}),v_{2}w_{2}).\end{split}$
(3.16)

This construction allows us to construct from any two elements of the
nullspaces of $D_{\text{TBG}}$, for some magic $\alpha_{0}$, with respect to
different representations, an element of the nullspace of $D_{\text{eTTG}},$
but for a rescaled magic parameter $\sqrt{2}\alpha_{0}.$ We observe that the
number of zeros doubles by applying the map above. This means that the number
of flat bands for eTTG, with $\sqrt{2}\alpha_{0}$ is (at least) twice the
number of flat bands for TBG with $\alpha_{0}.$

Figure 3. Magic angles of TBG multiplied by $\sqrt{2}$,
$\sqrt{2}\mathcal{A}_{\text{TBG}-2}$ (left) and eTTG
$\mathcal{A}_{\text{eTTG}-4}$ (right) with $U_{+}=U_{0}$, see (3.14). The
multiplicity of the flat bands of the Hamiltonian is twice the multiplicity
indicated in the figures.

###  3.3. Periodicity of TBG and eTTG

As discussed previously, the Hamiltonians for TBG and eTTG as defined in the
beginning of this section are not periodic with respect to $\Gamma$. However,
by using the general form of the potentials $U_{\pm}(\mathbf{r})$ (Eq. 3.7),
we see that we can turn TBG into a $\Gamma$-periodic Hamiltonian by
conjugating it by
$V(\mathbf{r}):=\operatorname{diag}(1,e^{i\mathbf{q}_{1}\cdot\mathbf{r}},1,e^{i\mathbf{q}_{1}\cdot\mathbf{r}})$
to obtain the equivalent periodic Hamiltonian

$\mathscr{H}(\alpha)=\begin{bmatrix}0&\mathscr{D}(\alpha)^{\dagger}\\\
\mathscr{D}(\alpha)&0\end{bmatrix}\text{ with
}\mathscr{D}(\alpha)=\begin{bmatrix}D_{x_{1}}+iD_{x_{2}}&\alpha
e^{iq_{0}\cdot\mathbf{r}}U_{+}(\mathbf{r})\\\ \alpha
e^{-iq_{0}\cdot\mathbf{r}}U_{-}(\mathbf{r})&D_{x_{1}}+iD_{x_{2}}+i\end{bmatrix}.$

Similarly, eTTG can be made periodic by conjugating by

$V(\mathbf{r})=\operatorname{diag}(1,e^{i\mathbf{q}_{1}\cdot\mathbf{r}},e^{-i\mathbf{q}_{1}\cdot\mathbf{r}},1,e^{i\mathbf{q}_{1}\cdot\mathbf{r}},e^{-i\mathbf{q}_{1}\cdot\mathbf{r}})$

and using the fact that $3\mathbf{q}_{1}=(0,3)^{\top}\in\Gamma^{*}$.

##  4\. Symmetries and Gauge Fixing

We now discuss the symmetry assumptions required for our main theorem which
mirror the symmetries present in TBG and eTTG. Due to the form of the chiral
Hamiltonian Eq. 3.1 and the symmetries of the $U_{\pm}(\mathbf{r})$ (Eqs. 3.4,
3.5 and 3.6), the Hamiltonians for TBG and eTTG satisfy a number of symmetries
in addition to the translation $\mathcal{T}_{\mathbf{a}}$ and rotation
$\mathcal{R}$. In particular, they satisfy a sublattice symmetry
$\mathcal{Z}$, a layer symmetry $\mathcal{L}$, and a composite symmetry
$\mathcal{Q}$, whose actions are defined as follows:

$\displaystyle\mathcal{Z}\psi(\mathbf{r};\sigma,j)$
$\displaystyle=\sigma\psi(\mathbf{r};\sigma,j)$ (4.1)
$\displaystyle\mathcal{L}\psi(\mathbf{r};\sigma,j)$
$\displaystyle=(-1)^{j}\psi(-\mathbf{r};-\sigma,N-j).$ (4.2)
$\displaystyle\mathcal{Q}\psi(\mathbf{r};\sigma,j)$
$\displaystyle=\overline{\psi(-\mathbf{r};-\sigma,j)}$ (4.3)

In words, the sublattice symmetry multiplies the $A$ sublattice by $+1$ and
states supported on the $B$ sublattice by $-1$. The layer symmetry
$\mathcal{L}$, reverses the order of the layers, swaps the $A$ and $B$
sublattices, multiplies by an alternating minus sign, and maps
$\mathbf{r}\mapsto-\mathbf{r}$. The symmetry $\mathcal{Q}$ is a composition of
a $C_{2z}$ rotation and time reversal $\mathcal{T}$ and acts by swapping the
$A$ and $B$ sublattices, taking a complex conjugate, and mapping
$\mathbf{r}\mapsto-\mathbf{r}$. The composite symmetry $\mathcal{Q}$ is often
referred to as the $C_{2z}\mathcal{T}$ symmetry. It is easily checked that
$\\{\mathcal{Z},\mathcal{Q}\\}=0$, $[\mathcal{Z},\mathcal{L}]=0$, and
$[\mathcal{Q},\mathcal{L}]=0$.

We will use these symmetries to fix a specific choice of Bloch eigenbasis for
the flat bands; a procedure we refer to as “gauge fixing” following the
physics terminology.

######  Assumption 3 (Symmetry Assumptions).

We assume that the single particle Hamiltonian $H$ commutes with $\mathcal{Q}$
and anticommutes with $\mathcal{Z}$ and $\mathcal{L}$.

######  Remark 2.

Since we study the nullspace of $H$, we can consider both symmetries that
commute or anticommute with $H$ as both types of symmetries fix the nullspace.

If $\mathcal{Z}$, $\mathcal{Q}$, and $\mathcal{L}$ are symmetries of a
periodic Hamiltonian, it is easily checked that $\mathcal{Z}$ and
$\mathcal{Q}$ map $\mathbf{k}$ to $\mathbf{k}$ and $\mathcal{L}$ maps
$\mathbf{k}$ to $-\mathbf{k}$.

Since the zero energy eigenstates are degenerate within the flat band, once we
fix a basis of orthogonal flat-band eigenfunctions we may perform any $U(2M)$
transformation to this basis and get an alternative choice of eigenbasis which
spans the same space. Any such transformation takes the form of the mapping

$\psi_{n\mathbf{k}}(\mathbf{r})\mapsto\sum_{m\in\mathcal{N}}\psi_{m\mathbf{k}}(\mathbf{r})[U(\mathbf{k})]_{mn}$

where $U(\mathbf{k})$ is a unitary matrix. Under this basis transformation,
the creation and annihilation operators likewise transform as

$\hat{f}^{\dagger}_{n\mathbf{k}}\mapsto\sum_{m\in\mathcal{N}}\hat{f}^{\dagger}_{m\mathbf{k}}[U(\mathbf{k})]_{mn}\qquad\qquad\hat{f}_{n\mathbf{k}}\mapsto\sum_{m\in\mathcal{N}}\hat{f}_{m\mathbf{k}}[U(\mathbf{k})^{*}]_{mn}$

and similarly the form factor transforms as

$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})\mapsto
U(\mathbf{k})^{\dagger}\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})U(\mathbf{k}+\mathbf{q}).$

By fixing a specific choice of eigenbasis, which we refer to as “fixing a
gauge” following physics terminology, we can force the form factor to satisfy
certain properties which will be a key part of our characterization of the
Hartree-Fock ground state. Note that the FBI Hamiltonian only depends on the
form factor $\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ which in turn only
depends on the periodic Bloch functions,
$u_{n\mathbf{k}}(\mathbf{r};\sigma,j)$. As such we will define our gauge
fixing in terms of the periodic Bloch functions,
$u_{n\mathbf{k}}(\mathbf{r};\sigma,j)$, instead of the Bloch functions.

######  Remark 3.

When the Hamiltonian $H$ exhibits band crossings at zero energy, the set
$\\{u_{n\mathbf{k}}(\mathbf{r}):n\in\mathcal{N},\mathbf{k}\in\Omega^{*}\\}$ is
still closed under the action of symmetry operations for
$\mathbf{k}\in\mathcal{K}_{\rm crossing}$ since the basis at the crossing
points is defined through continuity.

Since the zero energy eigenstates $\psi_{n\mathbf{k}}(\mathbf{r})$ are indexed
by $n\in\mathcal{N}$ and $\mathcal{Z}$ is a symmetry of the Hamiltonian, we
may change the gauge and relabel the eigenfunctions so that:

$\begin{split}\mathcal{Z}u_{n\mathbf{k}}(\mathbf{r};\sigma,j)&=(+1)u_{n\mathbf{k}}(\mathbf{r};\sigma,j)\qquad
n>0,\\\
\mathcal{Z}u_{n\mathbf{k}}(\mathbf{r};\sigma,j)&=(-1)u_{n\mathbf{k}}(\mathbf{r};\sigma,j)\qquad
n<0.\end{split}$ (4.4)

In particular, this implies that $u_{n\mathbf{k}}(\mathbf{r})$ is completely
supported on the $A$ sublattice if $n>0$ and completely supported on the $B$
sublattice if $n<0$. We can now properly define the uniformly half-filled,
translation invariant, ferromagnetic Slater determinant states:

######  Definition 4.1 (Ferromagnetic Slater Determinants).

Suppose that the Bloch eigenbasis has been chosen so that it satisfies Eq.
4.4. We define two uniformly half-filled, translation invariant ferromagnetic
Slater determinants, or ferromagnetic Slater determinants for short, to be
many-body states of the form:

$\ket{\Psi_{+}}=\prod_{\mathbf{k}\in\mathcal{K}}\prod_{n>0,n\in\mathcal{N}}\hat{f}_{n\mathbf{k}}^{\dagger}\ket{\mathrm{vac}}\qquad\ket{\Psi_{-}}=\prod_{\mathbf{k}\in\mathcal{K}}\prod_{n<0,n\in\mathcal{N}}\hat{f}_{n\mathbf{k}}^{\dagger}\ket{\mathrm{vac}}.$

That is, $\ket{\Psi_{\pm}}$ fully fills one of the two eigenspaces of
$\mathcal{Z}$.

Since $\mathcal{Q}$ is a symmetry of the Hamiltonian and relates the two
different sublattices, after fixing the Bloch eigenfunctions to satisfy Eq.
4.4, we may fix the span $\\{u_{n\mathbf{k}}(\mathbf{r}):n>0\\}$, and
determine $n<0$ by the sublattice symmetry $\mathcal{Q}$. That is, for all
$n>0$ we define

$u_{(-n)\mathbf{k}}(\mathbf{r};\sigma,j):=\mathcal{Q}u_{n\mathbf{k}}(\mathbf{r};\sigma,j)=\overline{u_{n\mathbf{k}}(-\mathbf{r};-\sigma,j)}.$
(4.5)

Additionally since $\mathcal{L}$ is also a symmetry of the Hamiltonian and
relates $\mathbf{k}$ and $-\mathbf{k}$, after fixing the periodic Bloch
functions to satisfy Eqs. 4.4 and 4.5 we may fix the space
$\\{u_{n\mathbf{k}}(\mathbf{r}):k_{1}\geq 0\\}$ and determine $k_{1}<0$ by the
layer symmetry $\mathcal{L}$. That is,

$u_{n(-k_{1},k_{2})}(\mathbf{r};\sigma,j):=\mathcal{L}u_{n(k_{1},-k_{2})}(\mathbf{r};\sigma,j)=(-1)^{j}u_{n(k_{1},-k_{2})}(-\mathbf{r};\sigma,N-j).$
(4.6)

######  Lemma 4.2.

Suppose that the gauge has been chosen so that periodic Bloch functions
$\\{u_{n\mathbf{k}}(\mathbf{r}):n\in\mathcal{N},\mathbf{k}\in\Omega^{*}\\}$
satisfy Eqs. 4.4, 4.5 and 4.6. Then the form factor
$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ can be written as

$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})=\begin{bmatrix}A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})&\\\
&\overline{A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})}\end{bmatrix}$ (4.7)

where

$[A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})]_{mn}=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}}\sum_{\sigma,j}\overline{\hat{u}_{m,\mathbf{k}}(\mathbf{G};\sigma,j)}\hat{u}_{n(\mathbf{k}+\mathbf{q})}(\mathbf{G}+\mathbf{G}^{\prime};\sigma,j)\qquad
m,n>0.$ (4.8)

Let $\mathcal{K}$ be as in Eq. 2.12, then
$A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ additionally satisfies the following
_sum rule_ :

$\sum_{\mathbf{k}\in\mathcal{K}}\operatorname{Im}{\operatorname{tr}{(A_{\mathbf{k}}(\mathbf{G}))}}=0,\quad\mathbf{G}\in\Gamma^{*}.$
(4.9)

###### Proof.

Recall the definition of the form factor

$[\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})]_{mn}:=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}\in\Gamma^{*}}\sum_{\sigma,j}\overline{\hat{u}_{m\mathbf{k}}(\mathbf{G}^{\prime};\sigma,j)}\hat{u}_{n(\mathbf{k}+\mathbf{q}+\mathbf{G})}(\mathbf{G}^{\prime};\sigma,j).$
(4.10)

Since the definition of the form factor involves a sum over the sublattice
$\sigma$ and $u_{n\mathbf{k}}(\mathbf{r})$ have disjoint sublattice support
for $n>0$ and $n<0$ we immediately see that
$[\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})]_{mn}=0$ if $mn<0$ and hence
$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ can be written as a block
diagonal matrix.

Due to $\mathcal{Q}$ symmetry (Eq. 4.5) we have that

$\begin{split}\hat{u}_{(-n)\mathbf{k}}(\mathbf{G};\sigma,j)&=\int_{\Omega}e^{-i\mathbf{G}\cdot\mathbf{r}}u_{(-n)\mathbf{k}}(\mathbf{r};\sigma,j)\,\mathrm{d}\mathbf{r}\\\
&=\int_{\Omega}e^{-i\mathbf{G}\cdot\mathbf{r}}\overline{u_{n\mathbf{k}}(-\mathbf{r};\sigma,j)}\,\mathrm{d}\mathbf{r}\\\
&=\overline{\int_{\Omega}e^{-i\mathbf{G}\cdot\mathbf{r}}u_{n\mathbf{k}}(\mathbf{r};\sigma,j)\,\mathrm{d}\mathbf{r}}\\\
&=\overline{\hat{u}_{n\mathbf{k}}(\mathbf{G};\sigma,j)}.\end{split}$ (4.11)

From this relation we easily see that

$[\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})]_{(-m)(-n)}=[\overline{\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})}]_{mn}$
(4.12)

which together with the block diagonal structure implies Eq. 4.7. Now we turn
to prove the sum rule Eq. 4.9.

Due to the layer symmetry (Eq. 4.6) we have that

$\begin{split}\hat{u}_{n(-\mathbf{k})}(\mathbf{G};\sigma,j)&=\int_{\Omega}e^{-i\mathbf{G}\cdot\mathbf{r}}(-1)^{j}u_{n\mathbf{k}}(-\mathbf{r};-\sigma,N-j)\,\mathrm{d}\mathbf{r}\\\
&=(-1)^{j}\hat{u}_{n\mathbf{k}}(-\mathbf{G};-\sigma,N-j).\end{split}$ (4.13)

Recall that $\Gamma^{*}$ is closed under the map
$\mathbf{G}\mapsto-\mathbf{G}$ and therefore

$\begin{split}\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})&=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}\in\Gamma^{*}}\sum_{\sigma,j}\overline{\hat{u}_{m\mathbf{k}}(\mathbf{G}^{\prime};\sigma,j)}\hat{u}_{n(\mathbf{k}+\mathbf{q}+\mathbf{G})}(\mathbf{G}^{\prime};\sigma,j)\\\
&=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}\in\Gamma^{*}}\sum_{\sigma,j}\overline{\hat{u}_{m(-\mathbf{k})}(-\mathbf{G}^{\prime};-\sigma,N-j)}\hat{u}_{n(-\mathbf{k}-\mathbf{q}-\mathbf{G})}(-\mathbf{G}^{\prime};-\sigma,N-j)\\\
&=\Lambda_{-\mathbf{k}}(-\mathbf{q}-\mathbf{G})\end{split}$ (4.14)

which for $\mathbf{q}=\bm{0}$ reduces to
$\Lambda_{\mathbf{k}}(\mathbf{G})=\Lambda_{-\mathbf{k}}(-\mathbf{G})=\Lambda_{-\mathbf{k}}(\mathbf{G})^{\dagger}$
where the last equality is due to Eq. 2.18. Note that this immediately implies
that $\Lambda_{\bm{0}}(\mathbf{G})$ is Hermitian.

We recall that elements of $\mathcal{K}$ can be written as
$\mathbf{k}=\frac{i}{n_{k_{x}}}\mathbf{g}_{1}+\frac{j}{n_{k_{y}}}\mathbf{g}_{2}$
for some $i\in[n_{k_{x}}]$ and $j\in[n_{k_{y}}]$. For any $\mathbf{k}$ of this
form, we can find a $\mathbf{k}^{\prime}\in\mathcal{K}$ and
$\mathbf{G}^{\prime}\in\Gamma^{*}$ so that
$-\mathbf{k}=\mathbf{k}^{\prime}+\mathbf{G}^{\prime}$. Therefore, we may
partition the momentum grid $\mathcal{K}$ into three disjoint sets
$\\{\bm{0}\\}$, $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$. The point $\bm{0}$ is
the unique point so that $\mathbf{k}=-\mathbf{k}$ and $\mathcal{K}_{1}$ and
$\mathcal{K}_{2}$ are defined so that for each $\mathbf{k}\in\mathcal{K}_{1}$
there exists a $\mathbf{k}^{\prime}\in\mathcal{K}_{2}$ so that
$\mathbf{k}-\mathbf{k}^{\prime}\in\Gamma^{*}$. Since
$A_{\mathbf{k}}(\mathbf{G})=A_{\mathbf{k}+\mathbf{G}^{\prime}}(\mathbf{G})$
for all $\mathbf{G}^{\prime}\in\Gamma^{*}$ (Eq. 2.19), we have

$\begin{split}\sum_{\mathbf{k}\in\mathcal{K}}A_{\mathbf{k}}(\mathbf{G})&=A_{\bm{0}}(\mathbf{G})+\sum_{\mathbf{k}\in\mathcal{K}_{1}}A_{\mathbf{k}}(\mathbf{G})+\sum_{\mathbf{k}\in\mathcal{K}_{2}}A_{\mathbf{k}}(\mathbf{G})\\\
&=A_{\bm{0}}(\mathbf{G})+\sum_{\mathbf{k}\in\mathcal{K}_{1}}A_{\mathbf{k}}(\mathbf{G})+\sum_{\mathbf{k}\in\mathcal{K}_{1}}A_{-\mathbf{k}}(\mathbf{G})\\\
&=A_{\bm{0}}(\mathbf{G})+\sum_{\mathbf{k}\in\mathcal{K}_{1}}(A_{\mathbf{k}}(\mathbf{G})+A_{\mathbf{k}}(\mathbf{G}))^{\dagger})\\\
\end{split}$ (4.15)

and hence
$\sum_{\mathbf{k}\in\mathcal{K}}\operatorname{Im}{\operatorname{tr}{(A_{\mathbf{k}}(\mathbf{G}))}}=0$.
Here we have used $A_{\bm{0}}(\mathbf{G})$ is Hermitian. ∎

##  5\. Properties of the Flat-Band Interacting Hamiltonian

In this section, we start by reviewing the basics on Hartree-Fock theory for
FBI Hamiltonians (Section 5.1) and then prove some of their important many-
body properties. In particular, we prove that the FBI Hamiltonian (Eq. 2.22)
is positive semidefinite (Section 5.2) and therefore any state which has zero
energy must be a ground state. Using this fact, along with the sum rule (Eq.
4.9), we prove that the ferromagnetic Slater determiant state are ground
states (Section 5.3). Finally, we show that the ferromagnetic Slater
determinants are insulating in the sense of a charge gap, i.e., both adding
and removing an electron costs a finite amount of energy (Section 5.4).

###  5.1. A Review of Hartree-Fock Theory

Recall that Slater determinants defined by the set $\mathcal{S}$ in Eq. 1.5
take the form:

$\ket{\Psi_{S}}=\prod_{i=1}^{M}\prod_{\mathbf{k}\in\mathcal{K}}\hat{b}_{i\mathbf{k}}^{{\dagger}}\ket{\mathrm{vac}},$
(5.1)

where $\ket{\mathrm{vac}}$ is the vacuum state, $M$ is the half the number of
bands, and

$\hat{b}_{i\mathbf{k}}^{{\dagger}}=\sum_{n\in\mathcal{N}}\hat{f}_{n\mathbf{k}}^{{\dagger}}\Xi_{ni}(\mathbf{k})\qquad\sum_{n\in\mathcal{N}}\Xi_{ni}^{*}(\mathbf{k})\Xi_{nj}(\mathbf{k})=\delta_{ij}$
(5.2)

defines the creation operator for the Hartree-Fock orbitals for each
$\mathbf{k}\in\mathcal{K}$.

The Hartree-Fock equations can be expressed in terms of the one-body reduced
density matrix (1-RDM). The 1-RDM associated with a given Slater determinant
$\ket{\Psi_{S}}$ can be written as

$[P(\mathbf{k})]_{nm}=\braket{\Psi_{S}}{\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n\mathbf{k}}}{\Psi_{S}}=\sum_{i=1}^{N_{\rm
occ}}\Xi_{ni}(\mathbf{k})\Xi_{mi}^{*}(\mathbf{k}).$ (5.3)

Due to the orthogonality relation on $\Xi(\mathbf{k})$ given in Eq. 5.2, we
may verify that for each $\mathbf{k}$ the 1-RDM $P(\mathbf{k})$ is an
orthogonal projection onto an $M$-dimensional vector space.

######  Remark 4.

The specific matrix representation of a 1-RDM depends on the choice of basis
used to define the creation and annihilation operators
$\hat{f}_{m\mathbf{k}}^{\dagger}$ and $\hat{f}_{m\mathbf{k}}$. Using the gauge
fixing scheme from Section 4, the two ferromagnetic Slater determinant states
have 1-RDMs as follows (where $\eta\in\\{\pm 1\\}$):

$\braket{\Psi_{\eta}}{\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n\mathbf{k}}}{\Psi_{\eta}}=\begin{cases}\delta_{mn}&\eta
m,\eta n>0\\\ 0&\eta m<0\text{ or }\eta n<0\end{cases}$ (5.4)

That is, the 1-RDM for $\ket{\Psi_{+}}$ is identity for all $\mathbf{k}$ and
all for $m,n>0$ and zero otherwise (and similarly for $\ket{\Psi_{-}}$).

Let us define the manifold of all admissible uniformly-filled 1-RDMs:

$\mathcal{M}:=\\{P\in{\mathbb{C}}^{(2M)\times(2M)}:P^{2}=P,\leavevmode\nobreak\
P^{\dagger}=P,\leavevmode\nobreak\ \operatorname{tr}{(P)}=M\\}.$

Following the standard derivation of Hartree-Fock theory (see e.g., [33, 23]),
the Hartree-Fock energy is given by finding the Slater determinant
$\ket{\Psi_{S}}$ which minimizes the energy of the interacting Hamiltonian:

$\mathcal{E}^{\rm(HF)}=\min_{\ket{\Psi_{S}}}\braket{\Psi_{S}}{\hat{H}_{FBI}}{\Psi_{S}}.$

Since $\hat{H}_{FBI}$ only involves the electron-electron interaction term,
the Hartree-Fock energy can be written as the sum of two functionals acting on
the 1-RDM $P\in\mathcal{M}^{\mathcal{K}}$. In particular,

$\mathcal{E}^{\rm(HF)}=\min_{P\in\mathcal{M}^{\mathcal{K}}}\Big{(}J[P]+K[P]\Big{)}$
(5.5)

where $J[\cdot]$ and $K[\cdot]$ are non-linear functionals in $P$, referred to
as the Hartree and Fock energy functionals, respectively.

Due to the specific form of $\hat{H}_{FBI}$, using Wick’s theorem, the Hartree
and Fock energies (up to a physically irrelevant constant) can be concisely
written in terms the matrix $Q(\mathbf{k}):=2P(\mathbf{k})-I$ as follows:

$\displaystyle J[P]$
$\displaystyle=\frac{1}{|\Omega|N_{\mathbf{k}}}\sum_{\mathbf{k},\mathbf{q}\in\mathcal{K}}\sum_{\mathbf{G}\in\Gamma^{*}}V(\mathbf{G})\operatorname{tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{G})Q(\mathbf{k})\Big{)}}\operatorname{tr}{\Big{(}\Lambda_{\mathbf{k}+\mathbf{q}}(\mathbf{G})^{\dagger}Q(\mathbf{k}+\mathbf{q})\Big{)}},$
(5.6) $\displaystyle K[P]$
$\displaystyle=-\frac{1}{|\Omega|N_{\mathbf{k}}}\sum_{\mathbf{k},\mathbf{q}\in\mathcal{K}}\sum_{\mathbf{G}\in\Gamma^{*}}V(\mathbf{q}+\mathbf{G})\operatorname{tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})Q(\mathbf{k}+\mathbf{q})\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})^{\dagger}Q(\mathbf{k})\Big{)}}.$
(5.7)

###  5.2. Positive Semidefiniteness

In this section, we prove that the many-body Hamiltonian $\hat{H}_{FBI}$ is
positive semidefinite. We start by proving the following

######  Lemma 5.1.

For all $\mathbf{q}^{\prime}\in\mathbb{R}^{2}$, the operator
$\hat{\rho}(\mathbf{q}^{\prime})$ satisfies
$\widehat{\rho}(\mathbf{q}^{\prime})^{\dagger}=\widehat{\rho}(-\mathbf{q}^{\prime})$.

###### Proof.

For any fixed $\mathbf{q}\in\mathcal{K}$ and $\mathbf{G}\in\Gamma^{*}$, we
calculate

$\begin{split}\widehat{\rho}(\mathbf{q}+\mathbf{G})^{\dagger}&=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})^{*}]_{mn}\left(\hat{f}_{n(\mathbf{k}+\mathbf{q})}^{\dagger}\hat{f}_{m\mathbf{k}}-\frac{1}{2}\delta_{\mathbf{q},0}\delta_{mn}\right)\\\
&=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})^{*}]_{nm}\left(\hat{f}_{m(\mathbf{k}+\mathbf{q})}^{\dagger}\hat{f}_{n\mathbf{k}}-\frac{1}{2}\delta_{\mathbf{q},0}\delta_{mn}\right).\end{split}$
(5.8)

Now we would like to make the change of variables
$\mathbf{k}\mapsto\mathbf{k}-\mathbf{q}$. Unfortunately, in general
$\mathbf{k},\mathbf{q}\in\mathcal{K}$ does not imply that
$\mathbf{k}-\mathbf{q}\in\mathcal{K}$. However, due to the definition of
$\mathcal{K}$, we may always find a
$\mathbf{G}_{\mathbf{k},\mathbf{q}}\in\Gamma^{*}$ so that

$\mathbf{k}-\mathbf{q}=\widetilde{\mathbf{k}-\mathbf{q}}+\mathbf{G}_{\mathbf{k},\mathbf{q}}\quad\text{where}\quad\widetilde{\mathbf{k}-\mathbf{q}}\in\mathcal{K}.$
(5.9)

Under the change of variables
$\mathbf{k}\mapsto\widetilde{\mathbf{k}-\mathbf{q}}$ we have that

$\begin{split}\widehat{\rho}(\mathbf{q}+\mathbf{G})^{\dagger}&=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}-\mathbf{q}+\mathbf{G}_{\mathbf{k},\mathbf{q}}}(\mathbf{q}+\mathbf{G})^{*}]_{nm}\left(\hat{f}_{m(\mathbf{k}+\mathbf{G}_{\mathbf{k},\mathbf{q}})}^{\dagger}\hat{f}_{n(\mathbf{k}-\mathbf{q}+\mathbf{G}_{\mathbf{k},\mathbf{q}})}-\frac{1}{2}\delta_{\mathbf{q},0}\delta_{mn}\right)\\\
&=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}-\mathbf{q}}(\mathbf{q}+\mathbf{G})^{*}]_{nm}\left(\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n(\mathbf{k}-\mathbf{q})}-\frac{1}{2}\delta_{\mathbf{q},0}\delta_{mn}\right)\\\
&=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}}(-\mathbf{q}-\mathbf{G})]_{mn}\left(\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n(\mathbf{k}-\mathbf{q})}-\frac{1}{2}\delta_{\mathbf{q},0}\delta_{mn}\right)=\widehat{\rho}(-\mathbf{q}-\mathbf{G})\end{split}$
(5.10)

where in the second line we have used Eq. 2.14 and Eq. 2.19. ∎

Since $\hat{H}_{FBI}$ takes the form

$\hat{H}_{FBI}=\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\widehat{\rho}(\mathbf{q}^{\prime})\widehat{\rho}(-\mathbf{q}^{\prime})=\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\widehat{\rho}^{{\dagger}}(\mathbf{q}^{\prime})\widehat{\rho}(\mathbf{q}^{\prime})$
(5.11)

and $\hat{V}(\mathbf{q}^{\prime})=\hat{V}(-\mathbf{q}^{\prime})>0$, and
immediate corollary is that any many-body state $\ket{\Psi}$ so that
$\hat{H}_{FBI}\ket{\Psi}=0$ must be a ground state.

###  5.3. Ferromagnetic Slater Determinants are Ground States

We can now check that the ferromagnetic Slater determinant states are ground
states by verifying they are zero energy eigenstates.

######  Proposition 5.2 (Proof of Result 1).

Suppose that the single particle Hamiltonian $H$ satisfies the symmetry
assumption (Assumption 3). Then the two ferromagnetic Slater determinant
states (4.1) are exact many-body ground states of the interacting model Eq.
2.22.

###### Proof.

We show that for all $\mathbf{q}^{\prime}=\mathbf{q}+\mathbf{G}$,
$\widehat{\rho}(\mathbf{q}^{\prime})\ket{\Psi_{\pm}}=0$ where
$\ket{\Psi_{\pm}}$ is a ferromagnetic Slater determinant. While we only
consider $\ket{\Psi_{+}}$, the calculation for $\ket{\Psi_{-}}$ follows
similar steps.

$\begin{split}\widehat{\rho}(\mathbf{q}+\mathbf{G})\ket{\Psi_{+}}&=\sum_{\mathbf{k}\in\mathcal{K}}\left[\sum_{m,n}[\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})]_{mn}\left(\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n(\mathbf{k}+\mathbf{q})}-\frac{1}{2}\delta_{mn}\delta_{\mathbf{q},\bm{0}}\right)\right]\ket{\Psi_{+}}\\\
&=\delta_{\mathbf{q},\bm{0}}\sum_{\mathbf{k}\in\mathcal{K}}\left[\sum_{m,n}[\Lambda_{\mathbf{k}}(\mathbf{G})]_{mn}\left(\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n\mathbf{k}}-\frac{1}{2}\delta_{mn}\right)\right]\ket{\Psi_{+}}\\\
&=\delta_{\mathbf{q},\bm{0}}\sum_{\mathbf{k}\in\mathcal{K}}\left[(\sum_{m>0}[\Lambda_{\mathbf{k}}(\mathbf{G})]_{mm})-\frac{1}{2}(\sum_{m}[\Lambda_{\mathbf{k}}(\mathbf{G})]_{mm})\right]\ket{\Psi_{+}}\\\
&=\frac{1}{2}\delta_{\mathbf{q},\bm{0}}\left[\sum_{\mathbf{k}\in\mathcal{K}}(\sum_{m>0}[\Lambda_{\mathbf{k}}(\mathbf{G})]_{mm}-[\Lambda_{\mathbf{k}}(\mathbf{G})^{*}]_{mm})\right]\ket{\Psi_{+}}\\\
&=\frac{1}{2}\delta_{\mathbf{q},\bm{0}}\left[\sum_{\mathbf{k}\in\mathcal{K}}\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}}(\mathbf{G}))\right]\ket{\Psi_{+}}=0\end{split}$
(5.12)

where in the second line, we have use the fact that
$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ is block diagonal and the
ferromagnetic Slater determinant state fully fills the $n>0$ states. The last
expression vanishes due to the sum rule in Eq. 4.9. ∎

###  5.4. Charge Gap at the Thermodynamic Limit

As we will see, adding or removing an electron to a ferromagnetic Slater
determinant always costs a non-zero amount of energy.

######  Proposition 5.3 (Charge Gap of Ferromagnetic Slater Determinants).

Fix $\eta\in\\{\pm 1\\}$ and consider an arbitrary single electron excitation
of the ferromagnetic Slater determinant $\ket{\Psi_{\eta}}$. All such
excitations can be written in terms of creation and annihilation operators of
the following form:

$\hat{c}^{\dagger}=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{\eta\ell<0}\hat{f}^{\dagger}_{\ell\mathbf{k}}c_{\ell\mathbf{k}}\qquad\hat{c}=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{\eta\ell>0}\hat{f}_{\ell\mathbf{k}}\overline{c_{\ell\mathbf{k}}}\qquad\sum_{\mathbf{k}\in\mathcal{K}}\sum_{\eta\ell>0}|c_{\ell,\mathbf{k}}|^{2}=1.$
(5.13)

where $\eta\ell<0$ denotes summation over the set $\\{\ell:\eta\ell<0\\}$ and
similarly for $\eta\ell>0$.

For such excitations, we have

$\begin{split}\braket{\Psi_{\eta}}{\hat{c}^{\dagger}\hat{H}_{FBI}\hat{c}}{\Psi_{\eta}}&=\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k}\in\mathcal{K}}\sum_{\eta\ell>0}\sum_{\eta\ell^{\prime}>0}c_{\ell^{\prime}\mathbf{k}}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})^{\dagger}]_{\ell^{\prime},\ell}\overline{c_{\ell\mathbf{k}}},\\\
\braket{\Psi_{\eta}}{\hat{c}\hat{H}_{FBI}\hat{c}^{\dagger}}{\Psi_{\eta}}&=\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k}\in\mathcal{K}}\sum_{\eta\ell<0}\sum_{\eta\ell^{\prime}<0}\overline{c_{\ell^{\prime}\mathbf{k}}}[\Lambda_{\mathbf{k}}(-\mathbf{q}^{\prime})^{\dagger}\Lambda_{\mathbf{k}}(-\mathbf{q}^{\prime})]_{\ell^{\prime},\ell}c_{\ell\mathbf{k}}.\end{split}$
(5.14)

Hence, since for all $\mathbf{k}$,
$[\Lambda_{\mathbf{k}}(\bm{0})]=\delta_{mn}$, the energy of
$\hat{c}^{\dagger}\ket{\Psi_{\eta}}$ and $\hat{c}\ket{\Psi_{\eta}}$ is
positive for all choices of $\hat{c}^{\dagger}$, $\hat{c}$.

Due to the constraint
$\sum_{\mathbf{k}}\sum_{\ell}|c_{\ell,\mathbf{k}}|^{2}=1$, in the
thermodynamic limit ($N_{\mathbf{k}}\rightarrow\infty$) the energy of the
single electron excitations converges to

$\begin{split}\bra{\Psi_{\eta}}&\hat{c}^{{\dagger}}\hat{H}_{FBI}\hat{c}\ket{\Psi_{\eta}}\\\
&=\int_{\mathbb{R}^{2}}\int_{\Omega}\hat{V}(\mathbf{q}^{\prime})\sum_{\eta\ell>0}\sum_{\eta\ell^{\prime}>0}c_{\ell^{\prime}}(\mathbf{k})[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})^{\dagger}]_{\ell^{\prime},\ell}\overline{c_{\ell}(\mathbf{k})}\,\mathrm{d}{\mathbf{k}}\,\mathrm{d}{\mathbf{q}^{\prime}},\\\
\bra{\Psi_{\eta}}&\hat{c}\hat{H}_{FBI}\hat{c}^{\dagger}\ket{\Psi_{\eta}}\\\
&=\int_{\mathbb{R}^{2}}\int_{\Omega}\hat{V}(\mathbf{q}^{\prime})\sum_{\eta\ell<0}\sum_{\eta\ell^{\prime}<0}\overline{c_{\ell^{\prime}}(\mathbf{k})}[\Lambda_{\mathbf{k}}(-\mathbf{q}^{\prime})^{\dagger}\Lambda_{\mathbf{k}}(-\mathbf{q}^{\prime})]_{\ell^{\prime},\ell}c_{\ell}(\mathbf{k})\,\mathrm{d}{\mathbf{k}}\,\mathrm{d}{\mathbf{q}^{\prime}}.\end{split}$
(5.15)

which is again strictly positive and so the charge gap does not vanish in the
thermodynamic limit. We now proceed to the proof.

###### Proof.

From the canonical anticommutation relations, we have that for all
$\mathbf{q}^{\prime}=\mathbf{q}+\mathbf{G}$

$\begin{split}[\hat{f}_{n\mathbf{k}}^{\dagger}\hat{f}_{m(\mathbf{k}+\mathbf{q}^{\prime})},\hat{f}_{\ell,\mathbf{k}^{\prime}}^{\dagger}]&=\delta_{m\ell}\delta_{\mathbf{k}^{\prime},(\mathbf{k}+\mathbf{q}^{\prime})}\hat{f}_{n\mathbf{k}}^{\dagger},\\\
[\hat{f}_{n\mathbf{k}}^{\dagger}\hat{f}_{m(\mathbf{k}+\mathbf{q}^{\prime})},\hat{f}_{\ell,\mathbf{k}^{\prime}}]&=-\delta_{n\ell}\delta_{\mathbf{k}^{\prime},\mathbf{k}}\hat{f}_{m(\mathbf{k}+\mathbf{q}^{\prime})}.\end{split}$
(5.16)

Therefore,

$\begin{split}[\widehat{\rho}(\mathbf{q}^{\prime}),\hat{f}_{\ell,\mathbf{k}^{\prime}}^{\dagger}]&=\sum_{\mathbf{k}}\sum_{mn}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}\delta_{n\ell}\delta_{\mathbf{k}^{\prime},(\mathbf{k}+\mathbf{q}^{\prime})}\hat{f}_{n\mathbf{k}}^{\dagger},=\sum_{n}[\Lambda_{\mathbf{k}^{\prime}-\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})]_{\ell
n}\hat{f}_{n(\mathbf{k}^{\prime}-\mathbf{q}^{\prime})}^{\dagger},\\\
[\widehat{\rho}(\mathbf{q}^{\prime}),\hat{f}_{\ell,\mathbf{k}^{\prime}}]&=-\sum_{\mathbf{k}}\sum_{mn}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}\delta_{n\ell}\delta_{\mathbf{k}^{\prime},\mathbf{k}}\hat{f}_{m(\mathbf{k}+\mathbf{q}^{\prime})}=-\sum_{m}[\Lambda_{\mathbf{k}^{\prime}}(\mathbf{q}^{\prime})]_{m\ell}\hat{f}_{m(\mathbf{k}^{\prime}+\mathbf{q}^{\prime})}.\\\
\end{split}$ (5.17)

Since the ferromagnetic Slater determinant state $\ket{\Psi_{\eta}}$ fully
fills either the positive or negative bands,
$\hat{f}_{\ell^{\prime}\mathbf{k}^{\prime\prime}}\ket{\Psi_{\eta}}\neq 0$ if
and only if $\eta\ell^{\prime}>0$. Therefore, for a fixed $\eta$, we only
consider $\ell^{\prime}$ so that $\eta\ell^{\prime}>0$. For the ferromagnetic
Slater determinant state $\ket{\Psi_{\eta}}$ for all $\mathbf{q}^{\prime}$ we
have

$\begin{split}\bra{\Psi_{\eta}}&\hat{f}_{\ell\mathbf{k}^{\prime}}^{\dagger}\widehat{\rho}(\mathbf{q}^{\prime})\widehat{\rho}(-\mathbf{q}^{\prime})\hat{f}_{\ell^{\prime}\mathbf{k}^{\prime\prime}}\ket{\Psi_{\eta}}=-\braket{\Psi_{\eta}}{[\widehat{\rho}(\mathbf{q}^{\prime}),\hat{f}_{\ell\mathbf{k}^{\prime}}^{\dagger}][\widehat{\rho}(-\mathbf{q}^{\prime}),\hat{f}_{\ell^{\prime}\mathbf{k}^{\prime\prime}}]}{\Psi_{\eta}}\\\\[4.30554pt]
&=\sum_{mn}[\Lambda_{\mathbf{k}^{\prime}-\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})]_{m\ell}[\Lambda_{\mathbf{k}^{\prime\prime}}(-\mathbf{q}^{\prime})]_{\ell^{\prime}n}\braket{\Psi_{\eta}}{\hat{f}_{m(\mathbf{k}^{\prime}-\mathbf{q}^{\prime})}^{\dagger}\hat{f}_{n(\mathbf{k}^{\prime\prime}-\mathbf{q}^{\prime})}}{\Psi_{\eta}}\\\
&=\sum_{\eta m>0}\sum_{\eta
n>0}[\Lambda_{\mathbf{k}^{\prime}-\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})]_{m\ell}[\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})]_{\ell^{\prime}n}\delta_{mn}\delta_{\mathbf{k}^{\prime}\mathbf{k}^{\prime\prime}}\\\
&=[\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})^{\dagger}]_{\ell^{\prime},\ell}\delta_{\mathbf{k}^{\prime}\mathbf{k}^{\prime\prime}}\end{split}$
(5.18)

where in the last line we have used the identity
$\Lambda_{\mathbf{k}-\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})=\Lambda_{\mathbf{k}}(-\mathbf{q}^{\prime})^{\dagger}$
Eq. 2.18.

Similar calculations show for $\eta\ell^{\prime}<0$

$\bra{\Psi_{\eta}}\hat{f}_{\ell\mathbf{k}^{\prime}}\widehat{\rho}(\mathbf{q}^{\prime})\widehat{\rho}(\mathbf{q}^{\prime})\hat{f}_{\ell^{\prime}\mathbf{k}^{\prime\prime}}^{\dagger}\ket{\Psi_{\eta}}=[\Lambda_{\mathbf{k}^{\prime}}(\mathbf{q}^{\prime})^{\dagger}\Lambda_{\mathbf{k}^{\prime}}(\mathbf{q}^{\prime})]_{\ell^{\prime},\ell}\delta_{\mathbf{k}^{\prime}\mathbf{k}^{\prime\prime}}.$
(5.19)

Therefore, by linearity

$\begin{split}\braket{\Psi_{\eta}}{\hat{c}^{\dagger}\widehat{\rho}(\mathbf{q}^{\prime})\widehat{\rho}(-\mathbf{q}^{\prime})\hat{c}}{\Psi_{\eta}}&=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{\eta\ell>0}\sum_{\eta\ell^{\prime}>0}c_{\ell^{\prime}\mathbf{k}}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})^{\dagger}]_{\ell^{\prime},\ell}\overline{c_{\ell\mathbf{k}}}\\\
\braket{\Psi_{\eta}}{\hat{c}\widehat{\rho}(\mathbf{q}^{\prime})\widehat{\rho}(-\mathbf{q}^{\prime})\hat{c}^{\dagger}}{\Psi_{\eta}}&=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{\eta\ell<0}\sum_{\eta\ell^{\prime}<0}\overline{c_{\ell^{\prime}\mathbf{k}}}[\Lambda_{\mathbf{k}}(-\mathbf{q}^{\prime})^{\dagger}\Lambda_{\mathbf{k}}(-\mathbf{q}^{\prime})]_{\ell^{\prime},\ell}c_{\ell\mathbf{k}}\end{split}$
(5.20)

which implies the proposition. ∎

##  6\. The Hartree-Fock Ground States of the Flat-Band Interacting
Hamiltonian

###  6.1. Rigorous Statement of Result 2

We can now state our main result rigorously:

######  Theorem 4 (Main theorem).

Suppose that the single particle Hamiltonian $H$ satisfies Assumptions 1 and
3. Suppose further that the Monkhorst-Pack grid $\mathcal{K}$ has been chosen
to satisfy Assumption 2. If there exists a $\mathbf{k}\in\mathcal{K}$ so that

  1. (1)

For some $\mathbf{G}$,
$\operatorname{Im}\operatorname{tr}{(A_{\mathbf{k}}(\mathbf{G}))}\neq 0$

  2. (2)

For all non-trivial orthogonal projectors $\Pi$ (i.e. $\Pi$ is not zero or
identity), there exists a $\mathbf{G}^{\prime}$ so that

$\|(I-\Pi)A_{\mathbf{k}}(\mathbf{G}^{\prime})\Pi\|>0$ (6.1)

then the two ferromagnetic Slater determinants are the unique translation-
invariant Hartree-Fock ground states of $\hat{H}_{FBI}$ in Eq. 2.22.

The conditions in Theorem 4 involves checking all orthogonal projectors and
are difficult to verify computationally. For the special case of $2$ and $4$
flat bands, these conditions can be significantly simplified.

######  Corollary 5 (Two Band Case).

Suppose that the single particle Hamiltonian $H$ satisfies Assumptions 1 and
3, and has two flat bands. Suppose further that the Monkhorst-Pack grid
$\mathcal{K}$ has been chosen to satisfy Assumption 2. If there exists a
$\mathbf{k}\in\mathcal{K}$ and $\mathbf{G}$ so that
$\operatorname{Im}(A_{\mathbf{k}}(\mathbf{G}))\neq 0$ then the two
ferromagnetic Slater determinants are the unique translation-invariant
Hartree-Fock ground states of Eq. 2.22.

###### Proof.

This follows immediately from Theorem 4 since in this case
$A_{\mathbf{k}}(\mathbf{G})$ is a scalar ($1\times 1$ matrix) and there are no
nontrivial orthogonal projectors for scalars. ∎

######  Corollary 6 (Four Band Case).

Suppose that the single particle Hamiltonian $H$ satisfies Assumptions 1 and
3, and has four flat bands. Suppose further that the Monkhorst-Pack grid
$\mathcal{K}$ has been chosen to satisfy Assumption 2. If there exists a
$\mathbf{k}\in\mathcal{K}$ so that

  1. (1)

For some $\mathbf{G}$,
$\operatorname{Im}\operatorname{tr}{(A_{\mathbf{k}}(\mathbf{G}))}\neq 0$,

  2. (2)

For some $\mathbf{G}^{\prime},\mathbf{G}^{\prime\prime}$,
$[A_{\mathbf{k}}(\mathbf{G}^{\prime}),A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})]\neq
0,$

then the two ferromagnetic Slater determinants are the unique translation-
invariant Hartree-Fock ground states of Eq. 2.22.

###### Proof.

The first condition is the same as Theorem 4 so we only need to show the
commutator condition implies the condition on projectors.

Since for four bands, $A_{\mathbf{k}}(\mathbf{G})$ is a $2\times 2$ matrix,
the only non-trivial projectors are rank one projectors. Therefore, we only
need to show that for all $\ket{v}\in{\mathbb{C}}^{2}$ with $\|v\|=1$, there
exists a $\mathbf{G}$ so that
$\braket{v^{\perp}}{A_{\mathbf{k}}(\mathbf{G})}{v}\neq 0$ where
$\ket{v^{\perp}}$ is a unit vector orthogonal to $\ket{v}$.

Suppose that $\mathbf{G}^{\prime},\mathbf{G}^{\prime\prime}$ are so that
$[A_{\mathbf{k}}(\mathbf{G}^{\prime}),A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})]\neq
0$ and pick some $\ket{v}\in{\mathbb{C}}^{2}$. We have three cases:

Case 1 $\braket{v^{\perp}}{A_{\mathbf{k}}(\mathbf{G}^{\prime})}{v}\neq 0$:

Take $\mathbf{G}=\mathbf{G}^{\prime}$

Case 2 $\braket{v^{\perp}}{A_{\mathbf{k}}(\mathbf{G}^{\prime})}{v}=0$ but
$\braket{v}{A_{\mathbf{k}}(\mathbf{G}^{\prime})}{v^{\perp}}\neq 0$:

Observe that

$\overline{\braket{v}{A_{\mathbf{k}}(\mathbf{G}^{\prime})}{v^{\perp}}}=\braket{v^{\perp}}{A_{\mathbf{k}}(\mathbf{G}^{\prime})^{\dagger}}{v}\neq
0$ (6.2)

but by Eq. 2.18
$A_{\mathbf{k}}(\mathbf{G}^{\prime})^{\dagger}=A_{\mathbf{k}}(-\mathbf{G}^{\prime})$
so we take $\mathbf{G}=-\mathbf{G}^{\prime}.$

Case 3 $\braket{v^{\perp}}{A_{\mathbf{k}}(\mathbf{G}^{\prime})}{v}=0$ and
$\braket{v}{A_{\mathbf{k}}(\mathbf{G}^{\prime})}{v^{\perp}}=0$:

Since $\ket{v}$ and $\ket{v^{\perp}}$ are an orthogonal basis for
${\mathbb{C}}^{2}$, the assumptions in this case implies that $\ket{v}$ and
$\ket{v^{\perp}}$ are eigenvectors of $A_{\mathbf{k}}(\mathbf{G}^{\prime})$.
Therefore, we can orthogonally diagonalize
$A_{\mathbf{k}}(\mathbf{G}^{\prime})$ as
$A_{\mathbf{k}}(\mathbf{G}^{\prime})=V\Sigma V^{\dagger}$ where $V$ is an
orthogonal matrix whose columns are $\ket{v}$, $\ket{v^{\perp}}$ and $\Sigma$
is a diagonal matrix of eigenvalues.

Since
$[A_{\mathbf{k}}(\mathbf{G}^{\prime}),A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})]\neq
0$, using the eigendecomposition of $A_{\mathbf{k}}(\mathbf{G}^{\prime})$ we
have that $V\Sigma
V^{\dagger}A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})-A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})V\Sigma
V^{\dagger}\neq 0$ which implies

$\Sigma\Big{(}V^{\dagger}A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})V\big{)}-\Big{(}V^{\dagger}A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})V\Big{)}\Sigma\neq
0$ (6.3)

If $V^{\dagger}A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})V$ were a diagonal
matrix then Eq. 6.3 would be zero so it must be that either
$\braket{v^{\perp}}{A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})}{v}\neq 0$ or
$\braket{v}{A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})}{v^{\perp}}\neq 0$. If
$\braket{v^{\perp}}{A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})}{v}\neq 0$ we
can take $\mathbf{G}=\mathbf{G}^{\prime\prime}$. If
$\braket{v}{A_{\mathbf{k}}(\mathbf{G}^{\prime\prime})}{v^{\perp}}\neq 0$,
appealing to Eq. 2.18, we can take $\mathbf{G}=-\mathbf{G}^{\prime\prime}$.

∎

##  7\. Proof of Theorem 4

As we saw in Section 4, with a proper choice of gauge, the sublattice symmetry
$\mathcal{Z}$ implies we can partition the set of flat bands $\mathcal{N}$
into two sets $n>0$ and $n<0$ whose basis functions are supported on the $A$
or $B$ sublattices respectively. The composite symmetry $\mathcal{Q}$ further
implies these two sets are of equal size and are related by an antiunitary
transformation.

To make use of this observation, suppose that we have a 1-RDM $P(\mathbf{k})$
for a model with $2M$ flat bands. If this state is uniformly half-filled, then
for each $\mathbf{k}\in\mathcal{K}$, $P(\mathbf{k})$ can be expressed as a
$(2M)\times(2M)$ projection matrix with rank $M$. Hence, we can write
$P(\mathbf{k})=\Phi(\mathbf{k})\Phi(\mathbf{k})^{\dagger}$ where
$\Phi(\mathbf{k})$ is a $(2M)\times M$ matrix with orthogonal columns.

Since $\Phi(\mathbf{k})$ has orthogonal columns, we may apply the cosine-sine
(CS) decomposition [35] to decompose $\Phi(\mathbf{k})$ so that it respects
the decomposition into $n>0$ and $n<0$:

$\Phi(\mathbf{k})=\begin{bmatrix}U_{1}(\mathbf{k})&\\\
&U_{2}(\mathbf{k})\end{bmatrix}\begin{bmatrix}\tilde{c}(\mathbf{k})&-\tilde{s}(\mathbf{k})\\\
\tilde{s}(\mathbf{k})&\tilde{c}(\mathbf{k})\end{bmatrix}\begin{bmatrix}V(\mathbf{k})^{\dagger}\\\
0\end{bmatrix}$ (7.1)

where

$\begin{split}&\tilde{c}(\mathbf{k})=\operatorname{diag}(\cos{(\theta_{1}(\mathbf{k})/2)},\cos{(\theta_{2}(\mathbf{k})/2)},\cdots,\cos{(\theta_{M}(\mathbf{k})/2)})\\\
&\tilde{s}(\mathbf{k})=\operatorname{diag}(\sin{(\theta_{1}(\mathbf{k})/2)},\sin{(\theta_{2}(\mathbf{k})/2)},\cdots,\sin{(\theta_{M}(\mathbf{k})/2)})\end{split}$
(7.2)

and $U_{1}(\mathbf{k}),U_{2}(\mathbf{k}),V(\mathbf{k})$ are $M\times M$
unitary matrices.

Using this decomposition for $\Phi(\mathbf{k})$, we can express the Hartree-
Fock energy (Eqs. 5.6 and 5.7) in terms of the quantities
$\theta_{i}(\mathbf{k}),U_{1}(\mathbf{k}),U_{2}(\mathbf{k})$ and the blocks of
the form factor $A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ (see Lemma 4.2 for the
definition of $A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$).

Since the Hartree-Fock energy is written in terms of the matrix
$Q(\mathbf{k})$, we begin writing $Q(\mathbf{k})$ in terms of the CS
decomposition Eq. 7.1. By definition we have

$\begin{split}P(\mathbf{k})&=\Phi(\mathbf{k})\Phi(\mathbf{k})^{\dagger}\\\
&=\begin{bmatrix}U_{1}(\mathbf{k})&\\\
&U_{2}(\mathbf{k})\end{bmatrix}\begin{bmatrix}\tilde{c}(\mathbf{k})^{2}&\tilde{c}(\mathbf{k})\tilde{s}(\mathbf{k})\\\
\tilde{c}(\mathbf{k})\tilde{s}(\mathbf{k})&\tilde{s}(\mathbf{k})^{2}\end{bmatrix}\begin{bmatrix}U_{1}(\mathbf{k})^{\dagger}&\\\
&U_{2}(\mathbf{k})^{\dagger}\end{bmatrix}.\end{split}$ (7.3)

Due to our decomposition of $\Phi(\mathbf{k})$, $P(\mathbf{k})$ is a
ferromagnetic Slater determinant if and only if one of the following holds:

  * •

For all $i\in\\{1,\cdots,M\\}$ and all $\mathbf{k}$,
$\theta_{i}(\mathbf{k})=0$, or

  * •

For all $i\in\\{1,\cdots,M\\}$ and all $\mathbf{k}$,
$\theta_{i}(\mathbf{k})=\pi$

These two conditions imply that for the ferromagnetic Slater determinant
states either $\tilde{c}(\mathbf{k})^{2}=I$ or $\tilde{s}(\mathbf{k})^{2}=I$.

Since $Q(\mathbf{k})=2P(\mathbf{k})-I$, we have

$Q(\mathbf{k})=\begin{bmatrix}U_{1}(\mathbf{k})&\\\
&U_{2}(\mathbf{k})\end{bmatrix}\begin{bmatrix}c(\mathbf{k})&s(\mathbf{k})\\\
s(\mathbf{k})&-c(\mathbf{k})\end{bmatrix}\begin{bmatrix}U_{1}(\mathbf{k})^{\dagger}&\\\
&U_{2}(\mathbf{k})^{\dagger}\end{bmatrix}$ (7.4)

where

$\begin{split}&c(\mathbf{k})=\operatorname{diag}(\cos{(\theta_{1}(\mathbf{k}))},\cos{(\theta_{2}(\mathbf{k}))},\cdots,\cos{(\theta_{M}(\mathbf{k}))})\\\
&s(\mathbf{k})=\operatorname{diag}(\sin{(\theta_{1}(\mathbf{k}))},\sin{(\theta_{2}(\mathbf{k}))},\cdots,\sin{(\theta_{M}(\mathbf{k}))}.\end{split}$
(7.5)

Using trigonometric identities, we see that the ferromagnetic Slater states
correspond to having $c(\mathbf{k})=\pm I$.

To prove uniqueness of the ferromagnetic Slater determinant states, we will
show that they are the unique states in $\mathcal{S}$ which achieve the
minimum value for both the Hartree and the Fock energies _simultaneously_. We
will first show that ferromagnetic Slater determinants are minimizers of the
Hartree energy in Section 7.1. Then in Section 7.2, we will show that the
assumptions of Theorem 4 imply that the ferromagnetic Slater determinant is
the unique minimizer of the Fock energy.

###  7.1. Minimizing the Hartree energy

Since $Q(\mathbf{k})=Q(\mathbf{k})^{\dagger}$,

$\overline{\sum_{\mathbf{k}\in\mathcal{K}}\operatorname{tr}{(\Lambda_{\mathbf{k}}(\mathbf{G})Q(\mathbf{k}))}}=\sum_{\mathbf{k}\in\mathcal{K}}\operatorname{tr}{(\Lambda_{\mathbf{k}^{\prime}}(\mathbf{G})^{\dagger}Q(\mathbf{k}^{\prime}))}.$
(7.6)

If $\mathbf{k},\mathbf{q}\in\mathcal{K}$, then
$\mathbf{q}-\mathbf{k}\in\mathcal{K}+\Gamma^{*}$. Since
$\Lambda_{\mathbf{k}+\mathbf{G}^{\prime}}(\mathbf{G})=\Lambda_{\mathbf{k}}(\mathbf{G})$
(Eq. 2.19) we can perform the change of variables
$(\mathbf{k},\mathbf{q})\rightarrow(\mathbf{k},\mathbf{q}-\mathbf{k})$ and
write the Hartree energy as

$J(P)=\frac{1}{|\Omega|N_{\mathbf{k}}}\sum_{\mathbf{G}}V(\mathbf{G})\left|\sum_{\mathbf{k}\in\mathcal{K}}\operatorname{tr}{(\Lambda_{\mathbf{k}}(\mathbf{G})Q(\mathbf{k}))}\right|^{2}.$
(7.7)

Note that since $V(\mathbf{G})>0$, necessarily $J[P]\geq 0$.

Now recall that, due to the symmetry assumptions on the form factor (Lemma
4.2) we can write

$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})=\begin{bmatrix}A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})&\\\
&\overline{A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})}\end{bmatrix}.$ (7.8)

Since we will multiply the form factor by $Q(\mathbf{k})$, it will be
convenient to define

$\begin{split}B^{(1)}_{\mathbf{k}}(\mathbf{G})&=U_{1}(\mathbf{k})A_{\mathbf{k}}(\mathbf{G})U_{1}(\mathbf{k})^{\dagger}\\\
B^{(2)}_{\mathbf{k}}(\mathbf{G})&=U_{2}(\mathbf{k})\overline{A_{\mathbf{k}}(\mathbf{G})}U_{2}(\mathbf{k})^{\dagger}.\end{split}$
(7.9)

The trace in the Hartree energy can then be written in terms of
$B^{(1)}_{\mathbf{k}}(\mathbf{G})$ and $B^{(2)}_{\mathbf{k}}(\mathbf{G})$ as
follows:

$\begin{split}\operatorname{tr}{(\Lambda_{\mathbf{k}}(\mathbf{G})Q(\mathbf{k}))}&=\operatorname{tr}{\Big{(}\begin{bmatrix}B^{(1)}_{\mathbf{k}}(\mathbf{G})&\\\
&B^{(2)}_{\mathbf{k}}(\mathbf{G})\end{bmatrix}\begin{bmatrix}c(\mathbf{k})&s(\mathbf{k})\\\
s(\mathbf{k})&-c(\mathbf{k})\end{bmatrix}\Big{)}}\\\
&=\operatorname{tr}{\Big{(}(B^{(1)}(\mathbf{G})-B^{(2)}(\mathbf{G}))c(\mathbf{k})\Big{)}}.\end{split}$
(7.10)

For the ferromagnetic Slater determinant state, $c(\mathbf{k})=\pm I$ and
hence

$\operatorname{tr}{(\Lambda_{\mathbf{k}}(\mathbf{G})Q(\mathbf{k}))}=\pm\operatorname{tr}{\Big{(}(B^{(1)}_{\mathbf{k}}(\mathbf{G})-B^{(2)}_{\mathbf{k}}(\mathbf{G}))\Big{)}}=\pm
2i\operatorname{Im}{\operatorname{tr}{(A_{\mathbf{k}}(\mathbf{G}))}}$ (7.11)
$J(P)=\frac{2}{|\Omega|N_{\mathbf{k}}}\sum_{\mathbf{G}}V(\mathbf{G})\left|\sum_{\mathbf{k}}\operatorname{Im}{\operatorname{tr}{(A_{\mathbf{k}}(\mathbf{G}))}}\right|^{2}=0$
(7.12)

where the last equality is due to the sum rule Eq. 4.9. Therefore, the
ferromagnetic Slater determinants minimize the Hartree energy.

###  7.2. Minimizing the Fock energy

For these calculations, we will adopt the shorthand
$\mathbf{k}^{\prime}:=\mathbf{k}+\mathbf{q}$ and
$\mathbf{q}^{\prime}=\mathbf{q}+\mathbf{G}$ and generalize the definitions of
$B^{(1)}_{\mathbf{k}}(\mathbf{G})$ and $B^{(2)}_{\mathbf{k}}(\mathbf{G})$ from
the previous section:

$\begin{split}B^{(1)}_{\mathbf{k}}(\mathbf{q}+\mathbf{G})&=B^{(1)}_{\mathbf{k}}(\mathbf{q}^{\prime})=U_{1}(\mathbf{k})A_{\mathbf{k}}(\mathbf{q}^{\prime})U_{1}(\mathbf{k}^{\prime})^{\dagger}\\\
B^{(2)}_{\mathbf{k}}(\mathbf{q}+\mathbf{G})&=B^{(2)}_{\mathbf{k}}(\mathbf{q}^{\prime})=U_{2}(\mathbf{k})\overline{A_{\mathbf{k}}(\mathbf{q}^{\prime})}U_{2}(\mathbf{k}^{\prime})^{\dagger}.\end{split}$
(7.13)

After some lengthy computations (Appendix A) it can be shown that

$\begin{split}K[P]=-\frac{1}{8|\Omega|N_{\mathbf{k}}}\sum_{\mathbf{k}}V(\mathbf{q}^{\prime})\sum_{ij}\bigg{\\{}&|[B^{(1)}_{\mathbf{k}}(\mathbf{q}^{\prime})+B^{(2)}_{\mathbf{k}}(\mathbf{q}^{\prime})]_{ij}|^{2}\cos{(\theta_{i}(\mathbf{k})-\theta_{j}(\mathbf{k}^{\prime}))}\\\
&+|[B^{(1)}_{\mathbf{k}}(\mathbf{q}^{\prime})-B^{(2)}_{\mathbf{k}}(\mathbf{q}^{\prime})]_{ij}|^{2}\cos{(\theta_{i}(\mathbf{k})+\theta_{j}(\mathbf{k}^{\prime}))}\bigg{\\}}.\end{split}$
(7.14)

Let’s take a closer look at each of terms in Eq. 7.14

$\begin{split}-|&[B^{(1)}_{\mathbf{k}}(\mathbf{q}^{\prime})+B^{(2)}_{\mathbf{k}}(\mathbf{q}^{\prime})]_{ij}|^{2}\cos{(\theta_{i}(\mathbf{k})-\theta_{j}(\mathbf{k}^{\prime}))}\\\
&-|[B^{(1)}_{\mathbf{k}}(\mathbf{q}^{\prime})-B^{(2)}_{\mathbf{k}}(\mathbf{q}^{\prime})]_{ij}|^{2}\cos{(\theta_{i}(\mathbf{k})+\theta_{j}(\mathbf{k}^{\prime}))}\end{split}$
(7.15)

From Eq. 7.15, we can understand the fundamental mechanism which forces the
ground state to be a ferromagnetic Slater determinant. Since
$\hat{V}(\mathbf{q}^{\prime})>0$, if we assume that

$|[B^{(1)}_{\mathbf{k}}(\mathbf{q}^{\prime})+B^{(2)}_{\mathbf{k}}(\mathbf{q}^{\prime})]_{ij}|>0\qquad|[B^{(1)}_{\mathbf{k}}(\mathbf{q}^{\prime})-B^{(2)}_{\mathbf{k}}(\mathbf{q}^{\prime})]_{ij}|>0$
(7.16)

then for a state to minimize the Fock energy it must be that

$\begin{split}\cos{(\theta_{i}(\mathbf{k})-\theta_{j}(\mathbf{k}^{\prime}))}=1&\Rightarrow\theta_{i}(\mathbf{k})=\theta_{j}(\mathbf{k}^{\prime})\pmod{2\pi}\\\
\cos{(\theta_{i}(\mathbf{k})+\theta_{j}(\mathbf{k}^{\prime}))}=1&\Rightarrow\theta_{i}(\mathbf{k})=-\theta_{j}(\mathbf{k}^{\prime})\pmod{2\pi}.\end{split}$
(7.17)

as if this weren’t the case then we could decrease the energy further.

The first constraint forces $\theta_{i}(\mathbf{k})$ to be constant
(independent of $i$ and $\mathbf{k}$) and the second constraint forces
$\theta_{i}(\mathbf{k})\in\\{0,\pi\\}$. These two facts combined show that the
ferromagnetic Slater determinants minimize the Fock energy and suggest a
strategy for proving these states are the unique Hartree-Fock minimizers.

By definition

$[B^{(1)}_{\mathbf{k}}(\mathbf{q}^{\prime})\pm
B^{(2)}_{\mathbf{k}}(\mathbf{q}^{\prime})]_{ij}=[U_{1}(\mathbf{k})A_{\mathbf{k}}(\mathbf{q}^{\prime})U_{1}(\mathbf{k}^{\prime})^{\dagger}\pm
U_{2}(\mathbf{k})\overline{A_{\mathbf{k}}(\mathbf{q}^{\prime})}U_{2}(\mathbf{k}^{\prime})^{\dagger}]_{ij}$
(7.18)

While generically, it may be true that that the above quantity does not
vanish, since $U_{1}(\mathbf{k})$ and $U_{2}(\mathbf{k})$ are arbitrary
unitaries, for any $i,j$ we can always find specific choices of
$U_{1}(\mathbf{k})$, $U_{2}(\mathbf{k})$ so that the above vanishes. The
assumptions of Theorem 4 guarantee that enough of these terms do not vanish
for every choice of $U_{1}(\mathbf{k})$, $U_{2}(\mathbf{k})$ to force the
ferromagnetic Slater determinant to be the unique minimizer of the Fock
energy.

To prove the assumptions of Theorem 4 are sufficient, we proceed in two steps.
First we show that for one special $\mathbf{k}$-point, $\mathbf{k}_{*}$,
enough of the entries of $B_{\pm,\mathbf{k}_{*}}(\mathbf{G})$ do not vanish to
force $\theta_{i}(\mathbf{k}_{*})=\theta_{j}(\mathbf{k}_{*})=0$ or
$\theta_{i}(\mathbf{k}_{*})=\theta_{j}(\mathbf{k}_{*})=\pi$ for all
$i,j\in\\{1,\cdots,M\\}$. This implies that the Fock energy maximizing 1-RDM
at $\mathbf{k}_{*}$ is a ferromagnetic Slater determinant. After showing this,
we use the fact that the grid $\mathcal{K}$ has been chosen sufficiently
finely so that Assumption 2 holds. Once this is the case, for any
$\mathbf{k}\in\mathcal{K}$, we can find a path connecting $\mathbf{k}_{*}$ and
$\mathbf{k}$ and we will show that for all momenta along this path, the
minimizing 1-RDM must be the same ferromagnetic Slater determinant as at
$\mathbf{k}_{*}$.

####  7.2.1. Local Uniqueness of Ground State

We focus on the point $\mathbf{k}_{*}$, take $\mathbf{q}=\bm{0}$ and $i=j$. In
this case, Eq. 7.15 reduces to

$\frac{1}{2}|[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G})+B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G})]_{ii}|^{2}+\frac{1}{2}|[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G})-B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G})]_{ii}|^{2}\cos{(2\theta_{i}(\mathbf{k}_{*}))}.$
(7.19)

Now notice that

$\begin{split}\operatorname{Im}{\operatorname{tr}{(A_{\mathbf{k}_{*}}(\mathbf{G}))}}&=\frac{1}{2i}\operatorname{tr}{(A_{\mathbf{k}_{*}}(\mathbf{G})-\overline{A_{\mathbf{k}_{*}}(\mathbf{G})})}\\\
&=\frac{1}{2i}\operatorname{tr}{\Big{(}U_{1}(\mathbf{k}_{*})A_{\mathbf{k}_{*}}(\mathbf{G})U_{1}(\mathbf{k}_{*})^{\dagger}-U_{2}(\mathbf{k}_{*})\overline{A_{\mathbf{k}_{*}}(\mathbf{G})}U_{2}(\mathbf{k}_{*})^{\dagger}\Big{)}}\\\
&=\frac{1}{2i}\operatorname{tr}{(B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G})-B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G}))}.\end{split}$
(7.20)

Since by assumption
$\operatorname{Im}{\operatorname{tr}{(A_{\mathbf{k}_{*}}(\mathbf{G}))}}\neq
0$, it must be there exists an $m$ so that
$[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G})-B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G})]_{mm}\neq
0$. Therefore, to be an optimizer $\cos{(2\theta_{m}(\mathbf{k}_{*}))}=1$
which implies $\theta_{m}(\mathbf{k}_{*})\in\\{0,\pi\\}$. Now we show that
$\theta_{j}(\mathbf{k}_{*})=\theta_{m}(\mathbf{k}_{*})$ for all $j$.

For this part of the proof, we fix the unitary $U_{1}(\mathbf{k}_{*})$ and
show that for this choice of $U_{1}(\mathbf{k}_{*})$ all of the
$\theta_{j}(\mathbf{k}_{*})$ must agree. For this let
$\\{\ket{i}:i\in\\{1,\cdots,M\\}\\}$ denote the standard basis for
${\mathbb{C}}^{M}$ so that for any matrix $A_{ij}=\braket{i}{A}{j}$.

We will prove this result by induction. Let $m_{1}=m$ and suppose that we have
already shown that the angles
$\\{\theta_{m_{i}}(\mathbf{k}_{*}):i\in\\{1,\cdots,n\\}\\}$ are all equal to
$\theta_{m_{1}}(\mathbf{k}_{*})\in\\{0,\pi\\}$. For fixed
$U_{1}(\mathbf{k}_{*})$, consider the orthogonal projector

$\Pi=U_{1}(\mathbf{k}_{*})^{\dagger}\left(\sum_{i=1}^{n}\ket{m_{i}}\bra{m_{i}}\right)U_{1}(\mathbf{k}_{*}).$
(7.21)

For any $\mathbf{G}$ we have

$\begin{split}(I&-\Pi)A_{\mathbf{k}_{*}}(\mathbf{G})\Pi\\\
&=\left(I-U_{1}(\mathbf{k}_{*})^{\dagger}\left(\sum_{i=1}^{n}\ket{m_{i}}\bra{m_{i}}\right)U_{1}(\mathbf{k}_{*})\right)A_{\mathbf{k}_{*}}(\mathbf{G})\left(U_{1}(\mathbf{k}_{*})^{\dagger}\left(\sum_{i=1}^{n}\ket{m_{i}}\bra{m_{i}}\right)U_{1}(\mathbf{k}_{*})\right)\\\\[4.30554pt]
&=U_{1}(\mathbf{k}_{*})^{\dagger}\left(I-\sum_{i=1}^{n}\ket{m_{i}}\bra{m_{i}}\right)U_{1}(\mathbf{k}_{*})A_{\mathbf{k}_{*}}(\mathbf{G})U_{1}(\mathbf{k}_{*})^{\dagger}\left(\sum_{i=1}^{n}\ket{m_{i}}\bra{m_{i}}\right)U_{1}(\mathbf{k}_{*}).\end{split}$
(7.22)

But since the spectral norm is unitarily invariant, the second assumption of
Theorem 4 implies that there exists a $\mathbf{G}^{\prime}$ so that

$\left\|\left(I-\sum_{i=1}^{n}\ket{m_{i}}\bra{m_{i}}\right)U_{1}(\mathbf{k}_{*})A_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})U_{1}(\mathbf{k}_{*})^{\dagger}\left(\sum_{i=1}^{n}\ket{m_{i}}\bra{m_{i}}\right)\right\|>0.$
(7.23)

Let $v,w$ be the top right/left singular vectors of the above operator. Since
$\sum_{i=1}^{n}\ket{m_{i}}\bra{m_{i}}$ is an orthogonal projection and
$\\{\ket{i}:i\in\\{1,\cdots,b\\}\\}$ forms a complete basis we can write

$v=\sum_{i=1}^{n}\alpha_{i}\ket{m_{i}},\qquad
w=\sum_{m\not\in\\{m_{i}:i\in\\{1,\cdots,n\\}\\}}\beta_{m}\ket{m}$ (7.24)

for some constants $\alpha_{i},\beta_{i}\in{\mathbb{C}}$. By definition of
$v,w$ we know that

$\begin{split}0&<\braket{w}{U_{1}(\mathbf{k}_{*})A_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})U_{1}(\mathbf{k}_{*})^{\dagger}}{v}\\\
&=\sum_{i=1}^{n}\sum_{m\not\in\\{m_{i}:i\in\\{1,\cdots,n\\}\\}}\alpha_{i}\overline{\beta_{m}}\braket{m}{U_{1}(\mathbf{k}_{*})A_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})U_{1}(\mathbf{k}_{*})^{\dagger}}{m_{i}}.\end{split}$
(7.25)

Hence, there must exist an $m^{\prime}\not\in\\{m_{i}:i\in\\{1,\cdots,n\\}\\}$
and an $m_{i}$ so that

$[U_{1}(\mathbf{k}_{*})A_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})U_{1}(\mathbf{k}_{*})^{\dagger}]_{m^{\prime},m_{i}}\neq
0.$ (7.26)

Define $m_{n+1}:=m^{\prime}$, since

$\Big{(}B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})+B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})\Big{)}+\Big{(}B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})-B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})\Big{)}=2U_{1}(\mathbf{k}_{*})A_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})U_{1}(\mathbf{k}_{*})^{\dagger}.$
(7.27)

Equation 7.26 implies that either
$[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})+B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})]_{m_{n+1},m_{i}}\neq
0$ or
$[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})-B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})]_{m_{n+1},m_{i}}\neq
0$. For simplicity of discussion suppose
$[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})+B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})]_{m_{n+1},m_{i}}\neq
0$, the other case follows similarly. Now recall the terms in the Fock energy:

$\begin{split}|[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})&+B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})]_{m_{n+1},m_{i}}|^{2}\cos{(\theta_{m_{n+1}}(\mathbf{k}_{*})-\theta_{m_{i}}(\mathbf{k}_{*}))}\\\
&+|[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})-B^{(2)}_{\mathbf{k}_{*}}(\mathbf{G}^{\prime})]_{m_{n+1},m_{i}}|^{2}\cos{(\theta_{m_{n+1}}(\mathbf{k}_{*})+\theta_{m_{i}}(\mathbf{k}_{*}))}.\end{split}$
(7.28)

Since $\theta_{m_{i}}(\mathbf{k}_{*})=\theta_{m_{1}}(\mathbf{k}_{*})$,
$\theta_{m_{1}}\in\\{0,\pi\\}$, and
$[B_{-,\mathbf{k}_{*}}(\mathbf{G}^{\prime})]_{m_{n+1},m_{i}}\neq 0$, to be a
minimizer it must be that
$\theta_{m_{n+1}}(\mathbf{k}_{*})=\theta_{m_{1}}(\mathbf{k}_{*})$. Hence by
induction, $\theta_{i}(\mathbf{k}_{*})\in\\{0,\pi\\}$ and
$\theta_{i}(\mathbf{k}_{*})=\theta_{j}(\mathbf{k}_{*})$ for all $i,j$.

####  7.2.2. Global Uniqueness of Ground State

For simplicity, let us assume that $\theta_{j}(\mathbf{k}_{*})=0$ for all $j$.
The case $\theta_{j}(\mathbf{k}_{*})=\pi$ follows similarly. In this case, for
any $\mathbf{q}$ we have that

$\begin{split}\frac{1}{2}&|[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})+B^{(2)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})]_{ij}|^{2}\cos{(\theta_{j}(\mathbf{k}_{*}+\mathbf{q}))}\\\
&+\frac{1}{2}|[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})-B^{(2)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})]_{ij}|^{2}\cos{(\theta_{j}(\mathbf{k}_{*}+\mathbf{q}))}.\end{split}$
(7.29)

Therefore, to conclude $\theta_{j}(\mathbf{k}_{*}+\mathbf{q})=0$, it is enough
to show that for each $r\in\\{1,\cdots,M\\}$

$\text{one
of}\quad\begin{cases}\max_{m}|[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})+B^{(2)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})]_{rm}|>0&\\\
\max_{m}|[B^{(1)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})-B^{(2)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})]_{rm}|>0&\\\
\end{cases}\quad\text{holds}.$ (7.30)

That is, for each row, we can find a non-zero entry in one of
$B^{(1)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})\pm
B^{(2)}_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})$

Now observe that

$B_{+,\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})+B_{-,\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})=U_{1}(\mathbf{k}_{*})A_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})U_{1}(\mathbf{k}_{*}+\mathbf{q})^{\dagger}$
(7.31)

and recall that a full rank matrix must have a non-zero entry in every row.
Since unitary transformations cannot change the rank of a matrix, if we can
show that $A_{\mathbf{k}_{*}}(\mathbf{q}+\mathbf{G})$ is full rank for some
$\mathbf{G}$, then Eq. 7.30 must be true and hence the fact that the
minimizing 1-RDM at $\mathbf{k}_{*}$ is an FSD state implies that the
minimizing 1-RDM at $\mathbf{k}_{*}+\mathbf{q}$ is the same FSD state.

The following key fact, along with Assumption 2 and the above calculation,
allows us to propagate the ferromagnetic Slater determinant from
$\mathbf{k}_{*}$ to the entire Brillouin zone:

######  Lemma 7.1.

Let $\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}$ and let $\Pi(\mathbf{k})$
denote the flat-band projection at $\mathbf{k}$. If
$\|\Pi(\mathbf{k})-\Pi(\mathbf{k}^{\prime})\|<1$ then there exists a
$\mathbf{G}$ so that
$A_{\mathbf{k}}((\mathbf{k}^{\prime}-\mathbf{k})+\mathbf{G})$ is full rank.

###### Proof.

Proven in Appendix B. ∎

Assuming Lemma 7.1 is true, the argument is as follows. By Assumption 2, for
any $\mathbf{k}\in\Omega^{*}$ we may construct a path
$\\{\mathbf{k}_{i}\\}_{i=1}^{L}$ connecting $\mathbf{k}_{*}$ and $\mathbf{k}$
so that $\|\Pi(\mathbf{k}_{i})-\Pi(\mathbf{k}_{i+1})\|<1$. By the above
calculation, since $\mathbf{k}_{1}=\mathbf{k}_{*}$ and the projector condition
holds between $\mathbf{k}_{1}$ and $\mathbf{k}_{2}$, it must be that the
minimizing 1-RDM at $\mathbf{k}_{2}$ is the same FSD state as at
$\mathbf{k}_{*}$. Since the minimizing 1-RDM at $\mathbf{k}_{2}$ is a FSD
state, we may repeat the same argument to conclude the 1-RDM at
$\mathbf{k}_{3}$ also the same FSD state. Continuing down the path, we
conclude that the minimizing 1-RDM at $\mathbf{k}_{i}$ must be the same FSD
state as at $\mathbf{k}_{*}$ for all $i\in\\{1,\cdots,L\\}$. Since the choice
of $\mathbf{k}$ was arbitrary, we conclude that minimizing 1-RDM must agree
with $\mathbf{k}_{*}$ throughout the whole Brillouin zone completing the
proof.

##  8\. Applications of Theorem 4 to TBG and eTTG

We begin by translating the conditions of Theorem 4 from momentum space to
real space in Section 8.1. We then prove Result 3 by verifying these real
space conditions hold for TBG-2 (Section 8.2), TBG-4 (Section 8.3), and eTTG-4
(Section 8.4).

###  8.1. Conditions of Theorem 4 in Real Space

While the proof of Theorem 4 is stated in momentum space, the conditions also
have natural analogs in real space. We begin by recalling the definition of
$A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$

$[A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})]_{mn}=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}}\sum_{\sigma,j}\overline{\hat{u}_{m\mathbf{k}}(\mathbf{G}^{\prime};\sigma,j)}\hat{u}_{n(\mathbf{k}+\mathbf{q})}(\mathbf{G}+\mathbf{G}^{\prime};\sigma,j).$
(8.1)

We can substitute in the definition of the Fourier transform to conclude

$\begin{split}[A_{\mathbf{k}}&(\mathbf{q}+\mathbf{G})]_{mn}\\\
&=\frac{1}{|\Omega|}\int_{\Omega}\int_{\Omega}\sum_{\mathbf{G}^{\prime}}\sum_{\sigma,j}e^{i\mathbf{G}^{\prime}\cdot\mathbf{r}}e^{-i(\mathbf{G}+\mathbf{G}^{\prime})\cdot\mathbf{r}^{\prime}}\overline{u_{m\mathbf{k}}(\mathbf{r};\sigma,j)}u_{n(\mathbf{k}+\mathbf{q})}(\mathbf{r}^{\prime};\sigma,j)\,\mathrm{d}\mathbf{r}\,\mathrm{d}\mathbf{r}^{\prime}\\\
&=\frac{1}{|\Omega|}\int_{\Omega}\int_{\Omega}\sum_{\mathbf{G}^{\prime}}\sum_{\sigma,j}e^{i\mathbf{G}^{\prime}\cdot(\mathbf{r}-\mathbf{r}^{\prime})}e^{-i\mathbf{G}\cdot\mathbf{r}^{\prime}}\overline{u_{m\mathbf{k}}(\mathbf{r};\sigma,j)}u_{n(\mathbf{k}+\mathbf{q})}(\mathbf{r}^{\prime};\sigma,j)\,\mathrm{d}\mathbf{r}\,\mathrm{d}\mathbf{r}^{\prime}\\\
&=\int_{\Omega}e^{-i\mathbf{G}\cdot\mathbf{r}}\sum_{\sigma,j}\overline{u_{m\mathbf{k}}(\mathbf{r};\sigma,j)}u_{n(\mathbf{k}+\mathbf{q})}(\mathbf{r};\sigma,j)\,\mathrm{d}\mathbf{r}.\end{split}$
(8.2)

We now define the pair product
$\rho_{\mathbf{k},\mathbf{k}+\mathbf{q}}(\mathbf{r})$

$[\rho_{\mathbf{k},\mathbf{k}+\mathbf{q}}(\mathbf{r})]_{mn}:=\sum_{\sigma,j}\overline{u_{m\mathbf{k}}(\mathbf{r};\sigma,j)}u_{n(\mathbf{k}+\mathbf{q})}(\mathbf{r};\sigma,j).$
(8.3)

Since $\rho_{\mathbf{k},\mathbf{k}+\mathbf{q}}$ is periodic with respect to
$\Gamma$ (since $u_{n\mathbf{k}}$ are periodic) we can view the form factor
$A_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ as the Fourier series coefficients of
the pair product. Since the Fourier transform is isometric up to scaling, we
can equivalently state conditions for Theorem 4 in real space.

######  Proposition 8.1 (Theorem 4 in Real Space).

Suppose that the single particle Hamiltonian $H$ satisfies Assumptions 1 and
3. Suppose further that the Monkhorst-Pack grid $\mathcal{K}$ has been chosen
to satisfy Assumption 2. If there exists a $\mathbf{k}\in\mathcal{K}$ so that

  1. (1)

For some $\mathbf{r}$,
$\operatorname{tr}{(\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r}))}\neq\overline{\operatorname{tr}{(\rho_{\mathbf{k},\mathbf{k}}(-\mathbf{r}))}}$

  2. (2)

For all non-trivial projections $\Pi$, there exists a $\mathbf{r}^{\prime}$ so
that

$\|(I-\Pi)\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r}^{\prime})\Pi\|>0$ (8.4)

then the two ferromagnetic Slater determinants are the unique Hartree-Fock
ground states of Eq. 2.22.

We also have the following simple corollaries:

######  Corollary 7 (Two Bands Case in Real Space).

Suppose that the single particle Hamiltonian $H$ satisfies Assumptions 1 and 3
and has two flat bands. Suppose further that the Monkhorst-Pack grid
$\mathcal{K}$ has been chosen to satisfy Assumption 2. If there exists a
$\mathbf{k}\in\mathcal{K}$ and $\mathbf{r}\in\Omega$ so that
$\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r})\neq\overline{\rho_{\mathbf{k},\mathbf{k}}(-\mathbf{r})}$
then the two ferromagnetic Slater determinants are the unique Hartree-Fock
ground states of Eq. 2.22.

######  Corollary 8 (Four Band Case in Real Space).

Suppose that the single particle Hamiltonian $H$ satisfies Assumptions 1 and 3
and has four flat bands. Suppose further that the Monkhorst-Pack grid
$\mathcal{K}$ has been chosen to satisfy Assumption 2. If there exists a
$\mathbf{k}\in\mathcal{K}$ so that

  1. (1)

For some $\mathbf{r}$,
$\operatorname{tr}{(\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r}))}\neq\overline{\operatorname{tr}{(\rho_{\mathbf{k},\mathbf{k}}(-\mathbf{r}))}}$,

  2. (2)

For some $\mathbf{r}^{\prime},\mathbf{r}^{\prime\prime}$

$\det{\begin{bmatrix}\|u_{1\mathbf{k}}(\mathbf{r})\|^{2}-\|u_{2\mathbf{k}}(\mathbf{r})\|^{2}&\braket{u_{1\mathbf{k}}(\mathbf{r}),u_{2\mathbf{k}}(\mathbf{r})}\\\
\|u_{1\mathbf{k}}(\mathbf{r}^{\prime})\|^{2}-\|u_{2\mathbf{k}}(\mathbf{r}^{\prime})\|^{2}&\braket{u_{1\mathbf{k}}(\mathbf{r}^{\prime}),u_{2\mathbf{k}}(\mathbf{r}^{\prime})}\end{bmatrix}}$
(8.5)

then the two ferromagnetic Slater determinants are the unique Hartree-Fock
ground states of Eq. 2.22.

For the proofs of Propositions 8.1, 7 and 8 we refer the reader to Appendix C.

###  8.2. Application of Main Theorem to TBG-2

For our next two propositions, we will use the Jacobi $\theta$ function

$\begin{gathered}\theta_{1}(\zeta|\omega):=-\sum_{n\in\mathbb{Z}}\exp(\pi
i(n+\tfrac{1}{2})^{2}\omega+2\pi
i(n+\tfrac{1}{2})(\zeta+\tfrac{1}{2})),\end{gathered}$ (8.6)

which satisfies

$\theta_{1}(\zeta+m|\omega)=(-1)^{m}\theta_{1}(\zeta|\omega),\ \
\theta_{1}(\zeta+n\omega|\omega)=(-1)^{n}e^{-\pi in^{2}\omega-2\pi i\zeta
n}\theta_{1}(\zeta|\omega).$

Furthermore, $\theta_{1}(\bullet|\omega)$ has simple zeros at $\Lambda$ (and
no other zeros) – see [27]. In the sequel we shall just write
$\theta(\zeta):=\theta_{1}(\zeta|\omega)$ to simplify the notation.

We also use the $\wp$ function $\wp(z):=\wp(z;\tfrac{4}{3}\pi
i\omega,\tfrac{4}{3}\pi i\omega^{2})$. Here $\wp(z;\omega_{1},\omega_{2})$ is
the Weierstrass $\wp$-function – see [27, §I.6]. It is periodic with respect
to ${\mathbb{Z}}\omega_{1}+{\mathbb{Z}}\omega_{2}$ and its derivative has a
pole of order $3$ at $z=0$. It has the property that

$\wp(\omega z)=\omega\wp(z).$ (8.7)

We now turn to verifying the conditions of having a unique many-body ground
state.

######  Proposition 8.2 (Simple magic angle).

Let $\alpha\in{\mathbb{C}}$ be a simple magic angle of TBG-2, then for all
$\mathbf{k}\notin\Gamma^{*}$ there is $\mathbf{G}\in\Gamma^{*}$ such that

$\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}}(\mathbf{G}))\neq 0$

while for $\mathbf{k}\in\Gamma^{*}$ one has

$\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}}(\mathbf{G}))=0\text{ for
all }\mathbf{G}\in\Gamma^{*}.$ (8.8)

###### Proof.

From Corollary 7 we deduce that the existence of some
$\mathbf{G}\in\Gamma^{*}$ such that

$\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}}(\mathbf{G}))\neq 0$

is equivalent to the existence of some $\mathbf{r}$ such that
$\|\mathbf{u}_{\mathbf{k}}(\mathbf{r})\|\neq\|\mathbf{u}_{\mathbf{k}}(-\mathbf{r})\|$,
where we use here the Euclidean norm of the $2$-vector
$\mathbf{u}_{\mathbf{k}}$ that is given by

$\mathbf{u}_{\mathbf{k}}\in\ker_{L^{2}_{2}}(D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{\mathbb{C}^{2}})\setminus\\{0\\}.$

At simple magic angles there is a unique solution $\mathbf{u}_{\mathbf{0}}$
which obey reflection symmetry, see [5, Theo.1] (also Figure 4 for an
illustration).

Figure 4. Log of modulus of flat-band wavefunction with $U_{0}$ as in (3.14)
and $\mathbf{k}=0$, $\alpha\approx 0.58656$ exhibiting reflection symmetry
with zero in the center (left) and $\mathbf{k}=(2,0.7)$ with zero away from
center (right) showing
$\|u_{\mathbf{k}}(\mathbf{r})\|\neq\|u_{\mathbf{k}}(-\mathbf{r})\|.$

Indeed, using symmetry (3.11), we find

$\|\mathbf{u}_{\mathbf{0}}(\mathbf{r})\|=\|\mathcal{L}\mathbf{u}_{\mathbf{0}}(\mathbf{r})\|=\|\mathbf{u}_{\mathbf{0}}(-\mathbf{r})\|$

showing (8.8). By [4, Theo. $3$] the function $\mathbf{u}_{\mathbf{0}}$ has a
simple zero at $\mathbf{r}=\mathbf{0}$ and no other zeros in its fundamental
domain.

From [4, Lemma $3.1$],

$\mathbf{u}_{\mathbf{k}}(\mathbf{r})=F_{\mathbf{k}}(\mathbf{r})\mathbf{u}_{\mathbf{0}}(\mathbf{r})\in\ker(D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{{\mathbb{C}}^{2}}),$
(8.9)

where

$F_{\mathbf{k}}(\mathbf{r}):=e^{\frac{(k_{1}+ik_{2})}{2}(-i(1+\omega)x+(\omega-1)y)}\frac{\theta(\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\frac{k_{1}+ik_{2}}{\sqrt{3}\omega})}{\theta\Big{(}\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}\Big{)}}.$

It follows that $\mathbf{u}_{\mathbf{k}}(\mathbf{r})$ has a unique zero at
$\mathbf{r}=\frac{4\pi}{3\sqrt{3}}(k_{2},-k_{1})^{\top}$ per unit cell (note
that both $\mathbf{u}_{\mathbf{0}}$ and $\theta$ have a simple zero at
$\mathbf{r}=\mathbf{0}$). Thus, $\|\mathbf{u}_{\mathbf{k}}\|$ cannot be an
even function for $\mathbf{k}\notin\Gamma^{*}.$ ∎

Therefore, by Corollary 7, we have the following result

######  Theorem 9.

The ferromagnetic Slater determinants states are the unique ground states of
the corresponding flat-band interacting model of TBG-2.

###  8.3. Application of Main Theorem to TBG-4

A similar result also holds for two-fold degenerate magic angles in chiral
limit TBG.

######  Proposition 8.3 (Two-fold degenerate magic angle).

Let $\alpha\in{\mathbb{C}}$ be a two-fold degenerate magic angle of TBG-4,
then for $\mathbf{k}=\pm\mathbf{q}_{1}$, are invariant points under the
rotation and translation symmetry, there are $\mathbf{G}\in\Gamma^{*}$ such
that

$\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}}(\mathbf{G}))\neq 0,$ (8.10)

while for $\mathbf{k}\in\Gamma^{*}$ one has

$\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}}(\mathbf{G}))=0\text{ for
all }\mathbf{G}\in\Gamma^{*}.$ (8.11)

###### Proof.

For a two-fold degenerate magic angle there is by [5, Theo.$1$], a unique
element $\mathbf{w}_{\mathbf{0}}\in\ker_{L^{2}_{2,1}}(D(\alpha))$ and
$\mathbf{v}_{\mathbf{0}}(\mathbf{r})=\wp(x_{1}+ix_{2})\mathbf{w}_{\mathbf{0}}(\mathbf{r})\in\ker_{L^{2}_{2,0}}(D(\alpha)).$

Since they belong to $L^{2}$-orthogonal subspaces, due to different rotational
symmetry using (8.7), we have
$\operatorname{tr}(A_{\mathbf{0}}(\mathbf{r}))=\|\mathbf{v}_{\mathbf{0}}(\mathbf{r})\|^{2}+\|\mathbf{w}_{\mathbf{0}}(\mathbf{r})\|^{2}.$
Thus, by the symmetry (3.11), we find

$\begin{split}\operatorname{tr}(A_{\mathbf{0}}(\mathbf{r}))&=\|\mathbf{w}_{\mathbf{0}}(\mathbf{r})\|^{2}+\|\mathbf{v}_{\mathbf{0}}(\mathbf{r})\|^{2}\\\
&=\|\mathcal{L}\mathbf{w}_{\mathbf{0}}(\mathbf{r})\|^{2}+\|\mathcal{L}\mathbf{v}_{\mathbf{0}}(\mathbf{r})\|^{2}\\\
&=\|\mathbf{w}_{\mathbf{0}}(-\mathbf{r})\|^{2}+\|\mathbf{v}_{\mathbf{0}}(-\mathbf{r})\|^{2}=\operatorname{tr}(A_{\mathbf{0}}(-\mathbf{r})).\end{split}$
(8.12)

We recall from Corollary 8 the equivalence

$\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}}(\mathbf{G}))\neq 0\text{
for some }\mathbf{G}\text{ if and only if
}\operatorname{tr}(\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r}))\neq\operatorname{tr}(\rho_{\mathbf{k},\mathbf{k}}(-\mathbf{r}))\text{
for some }\mathbf{r}$

showing (8.11), as all $\mathbf{k}\in\Gamma^{*}$ are equivalent to
$\mathbf{k}=\bm{0}.$

We recall from [5, Theo. $3$] and [5, Lemm. $7.1$] that
$\mathbf{w}_{\mathbf{0}}$ has a zero of order $2$ at $\mathbf{r}=\mathbf{0}$
and $\mathbf{v}_{\mathbf{0}}$ has simple zeros at $\pm\mathbf{r}_{S}$ with
$\mathbf{r}_{S}=(\frac{4\pi}{3\sqrt{3}},0)^{\top},$ see Figure 5.

Figure 5. Log of modulus of flat-band wavefunctions two flat bands for
$U_{7/8}$ as in (3.14) with $\mathbf{k}=0$, $\alpha\approx 0.853799$ with
$\mathbf{v}_{\mathbf{0}}$ (left) (with simple zeros at $\pm\mathbf{r}_{S}$ and
$\mathbf{w}_{\mathbf{0}}$(right) with double zero at $\mathbf{r}=0$.

We then define functions
$\mathbf{v}_{\mathbf{k}},\mathbf{w}_{\mathbf{k}}\in\ker(D(\alpha)+(k_{1}+ik_{2})\operatorname{id}_{{\mathbb{C}}^{2}})$
by

$\mathbf{v}_{\mathbf{k}}(\mathbf{r})=\alpha_{\mathbf{k}}F_{\mathbf{k}}(\mathbf{r}-\mathbf{r}_{S})\mathbf{v}_{\mathbf{0}}(\mathbf{r})\text{
and
}\mathbf{w}_{\mathbf{k}}(\mathbf{r})=\beta_{\mathbf{k}}F_{\mathbf{k}}(\mathbf{r})\mathbf{w}_{\mathbf{0}}(\mathbf{r})$

with normalizing constants $\alpha_{\mathbf{k}},\beta_{\mathbf{k}}>0$ that we
allow to change in this proof to simplify the notation and

$F_{\mathbf{k}}(\mathbf{r}):=e^{\frac{(k_{1}+ik_{2})}{2}(-i(1+\omega)x+(\omega-1)y)}\frac{\theta(\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\frac{k_{1}+ik_{2}}{\sqrt{3}\omega})}{\theta\Big{(}\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}\Big{)}}.$ (8.13)

From (8.9), we see that
$\mathbf{v}_{\mathbf{k}_{*}},\mathbf{w}_{\mathbf{k}_{*}}\in\ker(D(\alpha)+i\operatorname{id}_{{\mathbb{C}}^{2}}).$
We observe directly from the function $F_{\mathbf{k}}$ that
$\mathbf{w}_{\mathbf{k}_{*}}$ vanishes to first order at $\mathbf{r}=0$ and
$\mathbf{r}=\mathbf{r}_{S}$, but not at any other point. The function
$\mathbf{v}_{\mathbf{k}_{*}}$ vanishes to second order at
$\mathbf{r}=-\mathbf{r}_{S}$, i.e.
$\mathbf{v}_{\mathbf{k}_{*}}(-\mathbf{r}_{S})=0$ and
$\nabla\mathbf{v}_{\mathbf{k}_{*}}(-\mathbf{r}_{S})=0$.

We conclude that they are $L^{2}$-orthogonal using (8.7), since we can write
for $\gamma\neq 0$

$\begin{split}\mathbf{w}_{\mathbf{k}_{*}}(\mathbf{r}^{\prime})^{=}\gamma\wp(x_{1}^{\prime}+|r_{S}|+ix_{2}^{\prime})\mathbf{v}_{\mathbf{k}_{*}}(\mathbf{r}^{\prime}),\text{
where }\\\
\gamma^{:}=\frac{\mathbf{w}_{\mathbf{k}_{*}}(-\mathbf{r}_{S})}{\|\wp(\bullet+|\mathbf{r}_{S}|)\mathbf{v}_{\mathbf{k}_{*}}(\bullet)\|\lim_{x\to-\mathbf{r}_{S}}\wp(x+|\mathbf{r}_{S}|)\mathbf{v}_{\mathbf{k}_{*}}(x)}\end{split}$

and thus $\mathbf{w}_{x}$ and $\mathbf{v}_{\mathbf{k}_{*}}$ have different
rotational symmetries. We thus have

$\begin{split}\operatorname{tr}(\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}(\mathbf{r}))&=\|\mathbf{w}_{\mathbf{k}_{*}}(\mathbf{r})\|^{2}+\|\mathbf{v}_{\mathbf{k}_{*}}(\mathbf{r})\|^{2}\\\
&=(\alpha_{\mathbf{k}_{*}}|F_{\mathbf{k}_{*}}(\mathbf{r}-\mathbf{r}_{S})\wp(x_{1}+ix_{2})|^{2}+\beta_{\mathbf{k}_{*}}|F_{\mathbf{k}_{*}}(\mathbf{r})|^{2})\|\mathbf{w}_{\mathbf{0}}(\mathbf{r})\|^{2}.\end{split}$

In addition, we record that

$F_{\mathbf{k}_{*}}(\mathbf{r}-\mathbf{r}_{S})=e^{-\frac{2\pi}{3\sqrt{3}}(1+\omega)}e^{\frac{i}{2}(-i(1+\omega)x_{1}+(\omega-1)x_{2})}\frac{\theta(\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\frac{2i}{\sqrt{3}\omega})}{\theta\Big{(}\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\frac{i}{\sqrt{3}\omega}\Big{)}}$

and

$F_{\mathbf{k}_{*}}(\mathbf{r})=e^{\frac{i}{2}(-i(1+\omega)x_{1}+(\omega-1)x_{2})}\frac{\theta(\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\frac{i}{\sqrt{3}\omega})}{\theta\Big{(}\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}\Big{)}}.$

In particular, we have for the modulus of the prefactor of (8.13) as above

$|e^{\frac{i}{2}(-i(1+\omega)x_{1}+(\omega-1)x_{2})}|^{2}=e^{\frac{x_{1}-\sqrt{3}x_{2}}{2}}.$

Combining the above properties, we find

$\begin{split}\operatorname{tr}(\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}(\mathbf{r}))&=e^{\frac{x_{1}-\sqrt{3}x_{2}}{2}}\Bigg{(}\alpha_{\mathbf{k}_{*}}\frac{|\wp(x_{1}+ix_{2})\theta(\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\frac{2i}{\sqrt{3}\omega})|^{2}}{|\theta\Big{(}\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\frac{i}{\sqrt{3}\omega}\Big{)}|^{2}}+\beta_{\mathbf{k}_{*}}\frac{|\theta(\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\frac{i}{\sqrt{3}\omega})|^{2}}{|\theta\Big{(}\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}\Big{)}|^{2}}\Bigg{)}\|\mathbf{w}_{\mathbf{0}}(\mathbf{r})\|^{2}.\end{split}$

We also notice that $2i/(\sqrt{3}\omega)=-i/(\sqrt{3}\omega)+(1-\omega)$ which
allows us together with the identity

$\wp(x_{1}+ix_{2})=-e^{-\frac{2\pi}{\sqrt{3}\omega}}\Bigg{(}\frac{\theta^{\prime}(0)}{\theta(\frac{i}{\sqrt{3}\omega})}\Bigg{)}^{2}\frac{\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}-\tfrac{i}{\sqrt{3}\omega})\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\tfrac{i}{\sqrt{3}\omega})}{\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega})^{2}}$

to reduce the above expression to

$\begin{split}\operatorname{tr}(\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}(\mathbf{r}))&=e^{-(x_{1}-\sqrt{3}x_{2})}\Bigg{(}\alpha_{\mathbf{k}_{*}}\frac{|\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}-\tfrac{i}{\sqrt{3}\omega})|^{4}}{|\theta\Big{(}\frac{3(x_{1}+ix_{2})}{4\pi
i\omega}\Big{)}|^{2}}+\beta_{\mathbf{k}_{*}}|\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\tfrac{i}{\sqrt{3}\omega})|^{2}\Bigg{)}\frac{\|\mathbf{w}_{\mathbf{0}}(\mathbf{r})\|^{2}}{|\theta\Big{(}\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}\Big{)}|^{2}}.\end{split}$

Recalling that $\mathbf{w}_{0}$ and $|\theta|$ are even, we have to analyze
whether for $\gamma:=\beta_{\mathbf{k}_{*}}/\alpha_{\mathbf{k}_{*}}>0$

$\Phi(\mathbf{r})=e^{-(x_{1}-\sqrt{3}x_{2})}\Bigg{(}|\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}-\tfrac{i}{\sqrt{3}\omega})|^{4}+\gamma|\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\tfrac{i}{\sqrt{3}\omega})\theta\Big{(}\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}\Big{)}|^{2}\Bigg{)}$

is even. Evaluating $\Phi$ on the lattice $\Gamma$ we have by [27, Lemma
$4.1$]

$\mathbf{r}=(x_{1},x_{2})^{\top}\in\Gamma\Leftrightarrow\theta\Big{(}\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}\Big{)}=0,$

we find assuming $\Phi$ to be even that

$e^{-(x_{1}-\sqrt{3}x_{2})}|\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}-\tfrac{i}{\sqrt{3}\omega})|^{4}=\Phi(\mathbf{r})=\Phi(-\mathbf{r})=e^{(x_{1}-\sqrt{3}x_{2})}|\theta(\tfrac{3(x_{1}+ix_{2})}{4\pi
i\omega}+\tfrac{i}{\sqrt{3}\omega})|^{4}.$

Choosing, explicitly lattice points $x_{1}=2\pi/\sqrt{3}$ and $x_{2}=2\pi/3$,
which is just $-\mathbf{v}_{1}$, see (2.1), we find

$|\theta(-1-\tfrac{i}{\sqrt{3}\omega})|=|\theta(-1+\tfrac{i}{\sqrt{3}\omega})|$

which is false.

We have thus shown that

$\operatorname{tr}(\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}(\mathbf{r}))\neq\operatorname{tr}(\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}(-\mathbf{r})).$

Either repeating the previous construction for $-\mathbf{k}_{*}$ or observing
that using (3.11)

$\operatorname{tr}(\rho_{-\mathbf{k}_{*},-\mathbf{k}_{*}}(\mathbf{r}))=\|\mathcal{L}\mathbf{w}_{\mathbf{k}_{*}}(r)\|^{2}+\|\mathcal{L}\mathbf{v}_{\mathbf{k}_{*}}(\mathbf{r})\|^{2}=\|\mathbf{w}_{\mathbf{k}_{*}}(r)\|^{2}+\|\mathbf{v}_{\mathbf{k}_{*}}(\mathbf{r})\|^{2}=\operatorname{tr}(\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}(\mathbf{r})),$

we conclude that (8.10) holds. ∎

Using the notation of the previous proof, we also have

######  Proposition 8.4.

Let $\alpha$ be a two-fold degenerate magic angle of TBG,
$\mathbf{k}_{*}:=(0,1)^{\top}$, let
$\mathbf{r}\in\\{\mathbf{0},\pm\mathbf{r}_{S}\\}$, and
$\mathbf{r}^{\prime}\notin\\{\mathbf{0},\pm\mathbf{r}_{S}\\}+\Gamma$. Then the
expression

$\begin{split}D(\mathbf{r},\mathbf{r}^{\prime})&:=\det{\begin{bmatrix}\|\mathbf{w}_{\mathbf{k}}(\mathbf{r})\|^{2}-\|\mathbf{v}_{\mathbf{k}}(\mathbf{r})\|^{2}&\braket{\mathbf{w}_{\mathbf{k}}(\mathbf{r}),\mathbf{v}_{\mathbf{k}}(\mathbf{r})}\\\
\|\mathbf{w}_{\mathbf{k}}(\mathbf{r}^{\prime})\|^{2}-\|\mathbf{v}_{\mathbf{k}}(\mathbf{r}^{\prime})\|^{2}&\braket{\mathbf{w}_{\mathbf{k}}(\mathbf{r}^{\prime}),\mathbf{v}_{\mathbf{k}}(\mathbf{r}^{\prime})}\end{bmatrix}}\\\
&=(\|\mathbf{w}_{\mathbf{k}}(\mathbf{r})\|^{2}-\|\mathbf{v}_{\mathbf{k}}(\mathbf{r})\|^{2})\braket{\mathbf{w}_{\mathbf{k}}(\mathbf{r}^{\prime}),\mathbf{v}_{\mathbf{k}}(\mathbf{r}^{\prime})}\end{split}$
(8.14)

is non-zero.

###### Proof.

From the proof of Proposition 8.3, we know that $\mathbf{v}_{\mathbf{k}}$
vanishes to second order precisely at $-\mathbf{r}_{S}$ and
$\mathbf{w}_{\mathbf{k}}$ vanishes precisely to first order at $0$ and
$\mathbf{r}_{S}.$ This already implies the second line in (8.14) and that
$(\|\mathbf{w}_{\mathbf{k}}(\mathbf{r})\|^{2}-\|\mathbf{v}_{\mathbf{k}}(\mathbf{r})\|^{2})$
is non-zero. Thus, it remains to show that
$\braket{\mathbf{w}_{\mathbf{k}}(\mathbf{r}^{\prime}),\mathbf{v}_{\mathbf{k}}(\mathbf{r}^{\prime})}$
is non-zero. We can write up to a negligible phase
$\mathbf{w}_{\mathbf{k}}(\mathbf{r}^{\prime})=\wp(x_{1}^{\prime}+|r_{S}|+ix_{2}^{\prime})\mathbf{v}_{\mathbf{k}}(\mathbf{r}^{\prime})$
and thus

$\braket{\mathbf{w}_{\mathbf{k}}(\mathbf{r}^{\prime}),\mathbf{v}_{\mathbf{k}}(\mathbf{r}^{\prime})}=\wp(x_{1}^{\prime}+ix_{2}^{\prime}+|r_{S}|)\|\mathbf{v}_{\mathbf{k}}(\mathbf{r}^{\prime})\|^{2}$

and since $\wp(x_{1}^{\prime}+ix_{2}^{\prime}+|r_{S}|)$ has the same zeros as
$\mathbf{w}_{\mathbf{k}}$ and
$\mathbf{v}_{\mathbf{k}}(\mathbf{r}^{\prime})\neq 0$, the claim follows. ∎

Therefore, combining Propositions 8.3 and 8.4 with Corollary 8, we have the
following result

######  Theorem 10.

The ferromagnetic Slater determinants states are the unique ground states of
the corresponding flat-band interacting model for TBG-4.

###  8.4. Application of Main Theorem to eTTG-4

Recall the definition of the $\times$ operation in Eq. 3.16. If we denote by
$\mathbf{u}_{\mathbf{0}}\in L^{2}_{2,2}$ the flat band eigenfunction
associated in the nullspace of $D_{\text{TBG}}(\alpha_{0})$, then

$\mathbf{w}_{\mathbf{0}}:=\mathbf{u}_{\mathbf{0}}\times\mathbf{u}_{\mathbf{0}}\in\ker_{L^{2}_{1,1}}(D_{\text{eTTG}}(\sqrt{2}\alpha_{0})).$

We then have that

$\|\mathbf{w}_{\mathbf{0}}(\mathbf{r})\|=\|\mathbf{w}_{\mathbf{0}}(-\mathbf{r})\|.$

In particular, since $\mathbf{u}_{\mathbf{0}}$ has a simple zero at
$\mathbf{r}=0$, $\mathbf{w}_{\mathbf{0}}$ has a double zero at $\mathbf{r}=0$.

The two flat bands of eTTG-4 are given by the almost everywhere linearly
independent expressions

$w_{\mathbf{k}}(\mathbf{r})=F_{\mathbf{k}}(\mathbf{r})w_{0}(\mathbf{r}),\quad
w_{\mathbf{k}}(\mathbf{r})=F_{\mathbf{k}/2}(\mathbf{r})F_{\mathbf{k}/2}(\mathbf{r})w_{0}(\mathbf{r}).$

Since $\mathbf{w}_{\mathbf{0}}$ has therefore the same properties as
$\mathbf{w}_{\mathbf{0}}$ in the proof of Prop. 8.3. One can then replicate
the proofs of Prop. 8.3 and Prop. 8.4 to find that

######  Proposition 8.5.

Let $\alpha\in{\mathbb{C}}$ be a two-fold degenerate magic angle of eTTG-4,
with $\mathbf{k}_{*}$ as in Prop. 8.4, there are $\mathbf{G}\in\Gamma^{*}$
such that

$\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}_{*}}(\mathbf{G}))\neq 0$
(8.15)

while for $\mathbf{k}\in\Gamma^{*}$ one has

$\operatorname{Im}\operatorname{tr}(A_{\mathbf{k}}(\mathbf{G}))=0\text{ for
all }\mathbf{G}\in\Gamma^{*}.$ (8.16)

######  Proposition 8.6.

Let $\alpha$ be a two-fold degenerate magic angle of eTTG-4 and
$\mathbf{k}=\pm e_{2}$. Moreover, let
$\mathbf{r}\in\\{\mathbf{0},\pm\mathbf{r}_{S}\\}$ and
$\mathbf{r}^{\prime}\notin\\{\mathbf{0},\pm\mathbf{r}_{S}\\}+\Gamma$, then the
expression

$\begin{split}D(\mathbf{r},\mathbf{r}^{\prime})&:=\det{\begin{bmatrix}\|\mathbf{w}_{\mathbf{k}}(\mathbf{r})\|^{2}-\|\mathbf{v}_{\mathbf{k}}(\mathbf{r})\|^{2}&\braket{\mathbf{w}_{\mathbf{k}}(\mathbf{r}),\mathbf{v}_{\mathbf{k}}(\mathbf{r})}\\\