# On the Hartree-Fock Ground State Manifold in Magic Angle Twisted Graphene
Systems

Kevin D. Stubbs<EMAIL_ADDRESS>Department of Mathematics, University of
California, Berkeley, CA 94720, USA  ,  Simon Becker
<EMAIL_ADDRESS>ETH Zurich, Institute for Mathematical Research,
Rämistrasse 101, 8092 Zurich, Switzerland  and  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

###### Abstract.

Recent experiments have shown that magic angle twisted bilayer graphene
(MATBG) can exhibit correlated insulator behavior at half-filling. Seminal
theoretical results towards understanding this phase in MATBG has shown that
Hartree-Fock ground states (with a positive charge gap) can be exact many-body
ground states of an idealized flat band interacting (FBI) Hamiltonian. We
prove that in the absence of spin and valley degrees of freedom, the only
Hartree-Fock ground states of the FBI Hamiltonian for MATBG are two
ferromagnetic Slater determinants. Incorporating spin and valley degrees of
freedom, we provide a complete characterization of the Hartree-Fock ground
state manifold, which is generated by a ${\rm U}(4)\times{\rm U}(4)$ hidden
symmetry group acting on five elements. We also introduce new tools for ruling
out translation symmetry breaking in the Hartree-Fock ground state manifold,
which may be of independent interest.

##  1\. Introduction

Over the past few years, “magic angle” twisted bilayer graphene (MATBG) has
attracted immense attention in the condensed matter physics community due, in
part, to recent experiments which demonstrate correlated insulating [21, 31]
and superconducting phases [10]. The flat band interacting (FBI) Hamiltonian
at the chiral limit, which is an idealized model of MATBG at the moiré
scale111For MATBG, each moiré unit cell contains around $10^{4}$ carbon atoms.
Throughout the paper, details at the atomic scale are suppressed, and the
focus is solely on the periodic structure at the moiré scale. [23], suggests
that the ground state is both simple (at least in some regimes) and has a rich
structure [8, 11, 22, 30, 28, 12, 5, 16, 20, 14, 19, 17, 13, 18, 25, 32, 29,
24]. On one hand, the Hartree-Fock states, which are the simplest quantum
many-body states, can be _exact ground states_ of the FBI Hamiltonian at half-
filling222At other integer fillings, while Hartree-Fock states can be exact
eigenstates of the FBI Hamiltonian, they may not be ground states.. These
quantum states are ferromagnetic, meaning the density matrix behaves
identically at all $\mathbf{k}$ points. On the other hand, the FBI Hamiltonian
with spin and valley333Monolayer graphene features two nonequivalent Dirac
points within the Brillouin zone. Near each Dirac point there is a small
region, referred to as a valley. In MATBG, the wavefunction can be supported
on two valleys simultaneously, referred to as intervalley coherent (IVC)
states. degrees of freedom exhibits a large ${\rm U}(4)\times{\rm U}(4)$
“hidden symmetry” [8, 15]. This implies that the ground state of the FBI
Hamiltonian can, at most, be determined up to this symmetry. To the best of
our knowledge, whether the ferromagnetic ground states (subject to a ${\rm
U}(4)\times{\rm U}(4)$ symmetry) are the only Hartree-Fock ground states at
half-filling remains an open question.

###  1.1. Main Results

In this paper, we characterize the Hartree-Fock ground states of the flat band
interacting (FBI) Hamiltonian at half-filling, which is derived from the
Bistritzer-MacDonald model [7, 26, 9] at the chiral limit [23, 27, 1, 2]. We
refer the readers to Section 2 for a brief review of the FBI Hamiltonian or
[4, Section 2] for a more complete review. Ref. [4] proves that in a
simplified model of MATBG without spin or valley degrees of freedom, among the
Hartree-Fock at half-filling that are (1) uniformly filled, and (2)
translation invariant, there are only two possible ground states called the
ferromagnetic Slater determinant states.

This paper extends the findings of [4] in two ways. First, we prove that in
the absence of spin or valley degrees of freedom in MATBG, the two
ferromagnetic Slater determinant states are the only Hartree-Fock ground
states at half-filling. In particular, the “uniform filling” and “translation
invariant” restrictions have been lifted. Second, when incorporating spin and
valley degrees of freedom in MATBG, the manifold of Hartree-Fock ground states
at half-filling is characterized by a ${\rm U}(4)\times{\rm U}(4)$ symmetry
applied to five unique ferromagnetic Slater determinant states. This provides
a complete characterization of the Hartree-Fock ground states at half-filling
in MATBG.

Furthermore, the FBI Hamiltonian model can be generalized to other models of
twisted bilayer graphene and twisted multilayer graphene. Specifically, the
MATBG model above refers to the chiral limit twisted bilayer graphene
Hamiltonian with two flat bands at zero energy (TBG-2) at a magic angle. Ref.
[4] also analyzes the Hartree-Fock ground state for a general class of FBI
Hamiltonians under certain symmetry and non-degeneracy conditions. The non-
degeneracy conditions in [4] are generally computationally intractable. For
twisted bilayer graphene with four flat bands at a magic angle (TBG-4) and
equal twist angle trilayer graphene with four flat bands (eTTG-4), the non-
degeneracy conditions can be simplified and explicitly verified. As a
byproduct of this study, the non-degeneracy condition is replaced by a new
condition that is computationally more tractable. We also provide a sufficient
condition to rule out translation symmetry breaking for these general systems.

A key property of the FBI Hamiltonian for magic angle twisted graphene systems
is that they can be _frustration-free_. In fact, the Hamiltonians for these
systems can be expressed as a sum of positive semidefinite terms so that the
ground state energy is bounded below by zero. Moreover, there exist ground
states which minimize the energy of each individual term, also resulting in a
value of zero. As a consequence of the form of the FBI Hamiltonian, any many-
body ground state of a frustration-free FBI Hamiltonian is necessarily half-
filled (Lemma 2.1).

Moving to Hartree-Fock theory, our first result is that the Hartree-Fock
ground states of a frustration-free FBI Hamiltonian are always, in a certain
sense, “ferromagnetic”:

######  Result 1 (Informal Version of Proposition 1).

The Hartree-Fock ground states of a frustration-free FBI Hamiltonian are
always “ferromagnetic”, in the sense that knowing the 1-RDM at a single
momentum determines the value of the 1-RDM at all other momenta.

This property can be used to infer that the Hartree-Fock ground states of FBI
Hamiltonians are always uniformly filled. We next give a condition to rule out
translation breaking. In particular, we show that the translation symmetry
breaking of Hartree-Fock ground states of an FBI Hamiltonian is related to the
solution of a family of coupled Sylvester equations (see Section 5 for more
details). This characterization of the Hartree-Fock ground states implies the
following result:

######  Result 2 (Informal Version of Theorems 1, LABEL:, 2, LABEL: and 3).

Let the frustration-free FBI Hamiltonian for TBG-2 be denoted by $\hat{H}_{\rm
FBI}$.

  * •

If $\hat{H}_{\rm FBI}$ is valleyless and spinless, then the two ferromagnetic
Slater determinant states are the only Hartree-Fock ground states.

  * •

If $\hat{H}_{\rm FBI}$ is valleyful and spinless, then all Hartree-Fock ground
states are generated by the orbit of ${\rm U}(2)\times{\rm U}(2)$ on three
elements.

  * •

If $\hat{H}_{\rm FBI}$ is valleyful and spinful, then all Hartree-Fock ground
states are generated by the orbit ${\rm U}(4)\times{\rm U}(4)$ on five
elements.

Finally, we provide an algorithm in Proposition 3 for determining a spanning
set for all potential Hartree-Fock ground states.

###  1.2. Discussion and open question

In this paper, we find conditions which allow us to characterize all Hartree-
Fock ground states of frustration-free FBI Hamiltonians. We then apply these
conditions to the special case of magic angle twisted bilayer graphene with
two flat bands and extend this to include valley and spin degrees of freedom.
It would be interesting to apply the methodology here to other FBI
Hamiltonians corresponding to twisted multilayer graphene systems.

Considering the simplicity of the FBI model and the low amount of electron
correlation of the many-body ground state observed in numerical studies at
half-filling [22, 13], we propose the following conjecture:

######  Conjecture (Many-Body Ground States of Frustration-Free FBI
Hamiltonians).

Any many-body ground states of a frustration-free FBI Hamiltonian can be
written as a linear combination of its Hartree-Fock ground states. In
particular, for TBG-2 (with spin and valley), any many-body ground state can
be written as

(1) $\ket{\Psi}=\int_{{\rm U}(4)\times{\rm
U}(4)}\sum_{i=1}^{5}\alpha_{\mathfrak{g},i}\mathfrak{g}\ket{\Psi_{i}}\,\mathrm{d}\mathfrak{g},\quad\alpha_{\mathfrak{g},i}\in\mathbb{C}.$

Here $\\{\ket{\Psi_{i}}\\}_{i=1}^{5}$ are the ferromagnetic Slater determinant
states which generate the Hartree-Fock ground state manifold.

###  1.3. 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 discuss the Hartree-Fock theory of the FBI Hamiltonian, study
the Hartree-Fock energy, and find conditions for a state to be a Hartree-Fock
ground state.

  * •

In Section 4, we show that the Hartree-Fock ground states of FBI Hamiltonians
always satisfy a modified version of ferromagnetism as compared to the
ferromagnetic Slater determinant ground states discussed in our previous work.

  * •

In Section 5, we show that translation breaking for a Hartree-Fock ground
state can be characterized in terms of solutions to a family of coupled
Sylvester equations [6, Chapter VII.2]. We then show how this Sylvester
equation allows us to understand the appearance of the ${\rm U}(2)\times{\rm
U}(2)$ freedom in valleyful TBG and the ${\rm U}(4)\times{\rm U}(4)$ freedom
in valleyful and spinful TBG.

  * •

Our article contains two appendices. Appendix A contains the computations to
characterize the Hartree-Fock energy and Appendix B contains the proof of a
technical Lemma.

###  1.4. Notation

In this paper, we follow the same conventions as used in [4]. Operators and
matrices acting on the Fock space are represented using the hat notation, such
as $\hat{f}^{{\dagger}},\hat{f},\hat{H}_{\rm 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, such as operators $H$ and $D$. With some
slight abuse of the hat notation, vectors defined in real space are denoted
without the hat, for example, $f(\mathbf{r})$. Their corresponding Fourier
transforms are indicated with the hat notation, as in $\hat{f}(\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. Additionally, the
Frobenius norm of $A$ is denoted $\|A\|_{F}$. The identity is denoted by $I,$
where we occasionally write $I_{n\times n}$ to indicate the matrix size of the
identity.

### Acknowledgments

This work was supported by the Simons Targeted Grants in Mathematics and
Physical Sciences on Moiré Materials Magic (K.D.S., L.L.) and the SNF Grant
PZ00P2 216019 (S.B.). L.L. is a Simons Investigator in Mathematics. We thank
Dumitru Călugăru, Eslam Khalaf, Patrick Ledwidth, Oskar Vafek and Michael
Zaletel for helpful discussions.

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

The Flat-Band Interacting (FBI) Hamiltonian for twisted bilayer graphene is
based on a periodic single particle Hamiltonian which has exactly flat bands.
We highlight the main features of this model and refer the reader to [4,
Section 2] for more details. Following the notation in [4], we let $\Gamma$
and $\Omega$ denote the moiré lattice and real space unit cell in
$\mathbb{R}^{2}$ respectively. Similar, let $\Gamma^{*}$ and $\Omega^{*}$
denote the moiré reciprocal lattice and Brillouin zone in $\mathbb{R}^{2}$
respectively.

Suppose that we are given a single particle Hamiltonian $H$ with flat-bands
which is periodic with respect to $\Gamma$. Given such a Hamiltonian, let
$\\{\psi_{n\mathbf{k}}:\mathbf{k}\in\Omega^{*},n\in\mathcal{N}\\}$ be the
Bloch eigenfunctions where $\mathcal{N}$ is the set indexing the flat bands.
For a system with $N$-layers, each
$\psi_{n\mathbf{k}}(\mathbf{r}):=[\psi_{n\mathbf{k}}(\mathbf{r};\sigma,j)]\in\mathbb{C}^{2N},$
where $\sigma\in\\{A,B\\}$ denotes the sublattice degree of freedom and
$j\in\\{1,\cdots,N\\}$ denotes the layer degree of freedom. Furthermore, let
$u_{n\mathbf{k}}(\mathbf{r})=e^{-i\mathbf{k}\cdot\mathbf{r}}\psi_{n\mathbf{k}}(\mathbf{r})$
denote the $\Gamma$-periodic Bloch functions normalized so that
$\int_{\Omega}\|u_{n\mathbf{k}}(\mathbf{r})\|^{2}\,\mathrm{d}\mathbf{r}=1$,
where $\|u\|:=\sqrt{\sum_{i=1}^{n}|u_{i}|^{2}}$ is the $2$-norm of a vector
$u\in\mathbb{C}^{n}.$ Similarly, we denote by $\langle
v,u\rangle:=\sum_{i=1}^{n}\overline{v_{i}}u_{i}$ the canonical inner-product.

The moiré Brillouin zone $\Omega$ is discretized using a finite grid denoted
by $\mathcal{K}$

(2)
$\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^{*}.$

Here $\mathbf{g}_{1}$, $\mathbf{g}_{2}$ are a pair of generating vectors for
the moiré reciprocal lattice $\Gamma^{*}$ and $(n_{k_{x}},n_{k_{y}})$ are the
number of points in each of the two lattice directions. We also define
$N_{\mathbf{k}}:=\\#|\mathcal{K}|=n_{k_{x}}n_{k_{y}}$.

While in this work $n_{k_{x}},n_{k_{y}}$ are taken to be finite, we will be
interested in the limit where both $n_{k_{x}},n_{k_{y}}$ can become
arbitrarily large towards the thermodynamic limit. We make the following
technical assumption on the grid $\mathcal{K}$ in relation to the single
particle Hamiltonian:

######  Assumption 1.

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}\subseteq\mathcal{K}$ 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\\}$.

This assumption is satisfied for the chiral Bistritzer-MacDonald model [23, 2]
so long as $n_{k_{x}}$ and $n_{k_{y}}$ are both chosen sufficiently large.
This is because the flat bands are gapped from the remaining bands and the
Hamiltonian depends real-analytically on $\mathbf{k}.$

Having fixed the flat-band eigenfunctions and the grid $\mathcal{K}$, we now
define the FBI Hamiltonian. At each $\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}}$ respectively. These many-
body operators satisfy the canonical anti-commutation relation (CAR) as well
as a periodicity condition:

(3)
$\begin{split}\\{\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\\\
\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\forall\mathbf{G}\in\Gamma^{*}.\end{split}$

We also recall the definition of the number operator
$\hat{N}:=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m\in\mathcal{N}}\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{m\mathbf{k}}$.

Next, we define the Fourier coefficients of the periodic Bloch functions
$u_{n\mathbf{k}}$ introduced at the beginning of this section

(4)
$\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{k}+\mathbf{G})\cdot\mathbf{r}}\psi_{n\mathbf{k}}(\mathbf{r};\sigma,j)\,\mathrm{d}\mathbf{r}$

and define the form factor

(5)
$\begin{split}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}&:=\frac{1}{|\Omega|}\sum_{\mathbf{G}^{\prime}\in\Gamma^{*}}\braket{\hat{u}_{m\mathbf{k}}(\mathbf{G}^{\prime}),\hat{u}_{n(\mathbf{k}+\mathbf{q}^{\prime})}(\mathbf{G}^{\prime})}\\\
&=\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}^{\prime})}(\mathbf{G}^{\prime};\sigma,j)\end{split}$

where $\mathbf{k}\in\mathcal{K}$, and
$\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}$.

Finally, let $V$ be radially symmetric function with Fourier transform
$\hat{V}(\mathbf{q})>0$ for all $\mathbf{q}\in\mathbb{R}^{2}$. The FBI
Hamiltonian for $H$ with Coulomb potential $V$ takes the form

(6) $\begin{split}\hat{H}_{\rm
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}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\right).\end{split}$

The operator $\widehat{\rho}(\mathbf{q}^{\prime})$ satisfies
$\widehat{\rho}(-\mathbf{q}^{\prime})=\widehat{\rho}(\mathbf{q}^{\prime})^{\dagger}$
[4, Lemma 5.1], and so $\hat{H}_{\rm FBI}$ is positive semidefinite. Hence any
state with zero energy must be a ground state.

We note that from the form of the FBI Hamiltonian, we can immediately conclude
that any zero energy ground state must be half-filled.

######  Lemma 2.1.

Any many-body ground state $\ket{\Phi}$ of a frustration-free FBI Hamiltonian
is half-filled. That is, it satisfies

(7) $\hat{N}\ket{\Phi}=\frac{\\#|\mathcal{N}|N_{\mathbf{k}}}{2}\ket{\Phi}$

where $\hat{N}$ is the number operator
$\hat{N}=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m\in\mathcal{N}}\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{m\mathbf{k}}$.

###### Proof.

Let $\hat{H}_{\rm FBI}$ be a frustration-free Hamiltonian, any many-body
ground state $\ket{\Phi}$ must satisfy
$\braket{\Phi}{\widehat{\rho}(\mathbf{q}^{\prime})\widehat{\rho}(-\mathbf{q}^{\prime})}{\Phi}=0$.
From the definition of $\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})$, it can be
verified that for all $\mathbf{k}\in\Omega^{*}$,
$\Lambda_{\mathbf{k}}(\bm{0})=I$. In particular,

(8)
$\widehat{\rho}(\mathbf{q}^{\prime}=\bm{0})=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m\in\mathcal{N}}\left(\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{m\mathbf{k}}-\frac{1}{2}\right)=\hat{N}-\frac{\\#|\mathcal{N}|N_{\mathbf{k}}}{2}.$

To be a ground state, $\ket{\Phi}$ must satisfy

(9)
$\braket{\Phi}{\widehat{\rho}(\bm{0})\widehat{\rho}(\bm{0})}{\Phi}=\left\|\left(\hat{N}-\frac{\\#|\mathcal{N}|N_{\mathbf{k}}}{2}\right)\ket{\Phi}\right\|^{2}=0.$

Hence, $\hat{N}\ket{\Phi}=\frac{\\#|\mathcal{N}|N_{\mathbf{k}}}{2}\ket{\Phi}$
as was claimed. ∎

##  3\. Hartree-Fock Theory for Flat-Band Interacting Hamiltonians

We begin this section by first reviewing the Hartree-Fock theory at half-
filling without translation symmetry in Section 3.1. We then derive a simple
formula for the Hartree-Fock energy of the FBI Hamiltonian in Section 3.2. As
a consequence of this formula, we will obtain simple necessary and sufficient
conditions for a Hartree-Fock state to be a ground state of a frustration-free
FBI Hamiltonian.

###  3.1. A Review of Hartree-Fock Theory without Translation Symmetry

For a system with an even number of bands, recall that a general Slater
determinant at half-filling takes the form

(10)
$\ket{\Psi_{S}}=\prod_{i=1}^{MN_{\mathbf{k}}}\hat{b}_{i}^{{\dagger}}\ket{\mathrm{vac}},$

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

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

defines the creation operator for the Hartree-Fock orbitals.

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

(12)
$[P(\mathbf{k}^{\prime},\mathbf{k})]_{nm}=\braket{\Psi}{\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n\mathbf{k}^{\prime}}}{\Psi}=\sum_{i=1}^{MN_{\mathbf{k}}}\Xi_{ni}(\mathbf{k}^{\prime})\overline{\Xi_{mi}(\mathbf{k})}.$

Note that due to the periodicity of the creation and annihilation operators in
Eq. 3, the 1-RDM satisfies the property
$P(\mathbf{k}+\mathbf{G},\mathbf{k}^{\prime}+\mathbf{G}^{\prime})=P(\mathbf{k},\mathbf{k}^{\prime})$
for all $\mathbf{G},\mathbf{G}^{\prime}\in\Gamma^{*}$.

The 1-RDM $[P(\mathbf{k},\mathbf{k}^{\prime})]_{mn}$ may be represented as a
$(2MN_{\mathbf{k}})\times(2MN_{\mathbf{k}})$ matrix as follows

(13)
$P=\begin{bmatrix}P(\mathbf{k}_{1},\mathbf{k}_{1})&P(\mathbf{k}_{1},\mathbf{k}_{2})&\cdots&P(\mathbf{k}_{1},\mathbf{k}_{N_{\mathbf{k}}})\\\
P(\mathbf{k}_{2},\mathbf{k}_{1})&P(\mathbf{k}_{2},\mathbf{k}_{2})&\cdots&P(\mathbf{k}_{2},\mathbf{k}_{N_{\mathbf{k}}})\\\
\vdots&\vdots&\ddots&\vdots\\\
P(\mathbf{k}_{N_{\mathbf{k}}},\mathbf{k}_{1})&\cdots&\cdots&P(\mathbf{k}_{N_{\mathbf{k}}},\mathbf{k}_{N_{\mathbf{k}}})\\\
\end{bmatrix}.$

Due to the orthogonality relation on $\Xi(\mathbf{k})$ given in Eq. 11, one
may verify that $P$ satisfies $P^{2}=P$, $P^{\dagger}=P$, and
$\mathrm{Tr}{(P)}=MN_{\mathbf{k}}$ (i.e. $P$ is an orthogonal projection onto
a $(MN_{\mathbf{k}})$-dimensional vector space). The 1-RDM is translation-
invariant if it takes the form

(14) $\begin{bmatrix}P(\mathbf{k}_{1},\mathbf{k}_{1})&\bm{0}&\cdots&\bm{0}\\\
\bm{0}&P(\mathbf{k}_{2},\mathbf{k}_{2})&\cdots&\bm{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\bm{0}&\cdots&\cdots&P(\mathbf{k}_{N_{\mathbf{k}}},\mathbf{k}_{N_{\mathbf{k}}})\end{bmatrix}.$

###  3.2. Hartree-Fock Energy of Flat-Band Interacting Hamiltonians

The Hartree-Fock energy for FBI Hamiltonians can also be written in terms of
the shifted 1-RDM, $Q(\mathbf{k},\mathbf{k}^{\prime})$, defined as follows:

(15)
$[Q(\mathbf{k},\mathbf{k}^{\prime})]_{mn}:=[P(\mathbf{k},\mathbf{k}^{\prime})]_{mn}-\frac{1}{2}\delta_{mn}\delta_{\mathbf{k}-\mathbf{k}^{\prime}\in\Gamma^{*}}.$

Since $P^{\dagger}=P$ and $P^{2}=P$, one easily checks that $Q^{\dagger}=Q$
and $Q^{2}=\frac{1}{4}I$. Furthermore, $Q(\mathbf{k},\mathbf{k}^{\prime})$
inherits the property
$Q(\mathbf{k}+\mathbf{G},\mathbf{k}^{\prime}+\mathbf{G}^{\prime})=Q(\mathbf{k},\mathbf{k}^{\prime})$
for all $\mathbf{G},\mathbf{G}^{\prime}\in\Gamma^{*}$ from the 1-RDM $P$.

The Hartree-Fock energy can be written in terms of
$Q(\mathbf{k},\mathbf{k}^{\prime})$ as follows (Appendix A) :

(16) $\begin{split}\langle&\Psi|\hat{H}_{\rm
FBI}|\Psi\rangle=\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\left|\sum_{\mathbf{k}\in\mathcal{K}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})Q(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\right|^{2}\\\
&+\frac{1}{4N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k}\in\mathcal{K}}\|\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})\|_{F}^{2}\\\
&-\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})Q(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k}^{\prime})\Lambda_{\mathbf{k}^{\prime}-\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})^{\dagger}Q(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k})\Big{)}}.\end{split}$

We can further simplify this expression by expanding the families of
$(2M)\times(2M)$ matrices
$\\{\Lambda_{\mathbf{k}^{\prime}}(\mathbf{q}^{\prime}):\mathbf{k}^{\prime}\in\mathcal{K}\\}$
and
$\\{Q(\mathbf{k},\mathbf{k}^{\prime}):\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}\\}$
into larger $(2MN_{\mathbf{k}})\times(2MN_{\mathbf{k}})$ matrices similar to
Eq. 13. We start by defining notation for the standard basis on the space of
1-RDMs

(17) $\\{\ket{n,\mathbf{k}}:n\in\mathcal{N},\mathbf{k}\in\mathcal{K}\\}$

where we take the convention that if
$\mathbf{k}+\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}$ then

(18)
$\ket{n,\mathbf{k}+\mathbf{q}^{\prime}}:=\ket{n,\widetilde{\mathbf{k}+\mathbf{q}^{\prime}}}$

where $\widetilde{\mathbf{k}+\mathbf{q}^{\prime}}\in\mathcal{K}$ is the unique
representative of $\mathbf{k}+\mathbf{q}^{\prime}$ in $\mathcal{K}$ modulo
$\Gamma^{*}$ lattice translations. We additionally define the momentum shift
matrix:

(19)
$\pi_{\mathbf{q}^{\prime}}:=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m\in\mathcal{N}}\ket{m,\mathbf{k}}\bra{m,\mathbf{k}+\mathbf{q}^{\prime}}.$

This shifts momentum $\mathbf{k}$ to momentum $\mathbf{k}+\mathbf{q}^{\prime}$
modulo $\Gamma^{*}$ in the basis $\ket{m,\mathbf{k}}$.

Using this basis, for any $\mathbf{q}^{\prime}$ we define the matrix

(20)
$\Lambda(\mathbf{q}^{\prime}):=\sum_{\mathbf{k}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}\ket{m,\mathbf{k}}\bra{n,\mathbf{k}}$

which can be expressed as a direct sum as follows

(21)
$\Lambda(\mathbf{q}^{\prime})=\bigoplus_{n=1}^{N_{\mathbf{k}}}\Lambda_{\mathbf{k}_{n}}(\mathbf{q}^{\prime})=\begin{bmatrix}\Lambda_{\mathbf{k}_{1}}(\mathbf{q}^{\prime})&&&&\\\
&\Lambda_{\mathbf{k}_{2}}(\mathbf{q}^{\prime})&&&\\\
&&\Lambda_{\mathbf{k}_{3}}(\mathbf{q}^{\prime})&&\\\ &&&\ddots&\\\
&&&&\Lambda_{\mathbf{k}_{N_{\mathbf{k}}}}(\mathbf{q}^{\prime})\par\end{bmatrix}.$

Furthermore, the matrix $Q$ can be written in this basis as

(22)
$Q:=\sum_{\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[Q(\mathbf{k},\mathbf{k}^{\prime})]_{mn}\ket{m,\mathbf{k}}\bra{n,\mathbf{k}^{\prime}}.$

Note that in this basis $Q=P-\frac{1}{2}I$. Now observe that

(23)
$\begin{split}\pi_{\mathbf{q}^{\prime}}Q&=\sum_{\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[Q(\mathbf{k},\mathbf{k}^{\prime})]_{mn}\ket{m,\mathbf{k}-\mathbf{q}^{\prime}}\bra{n,\mathbf{k}^{\prime}}\\\
&=\sum_{\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[Q(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k}^{\prime})]_{mn}\ket{m,\mathbf{k}}\bra{n,\mathbf{k}^{\prime}},\end{split}$

where the shift in $\ket{m,\mathbf{k}-\mathbf{q}^{\prime}}$ in the first line
is due to the multiplication on the right by $\pi_{\mathbf{q}^{\prime}}$.
Performing the change of variables
$\mathbf{k}\mapsto\mathbf{k}+\mathbf{q}^{\prime}$ is justified, since
$\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}$ and
$Q(\mathbf{k},\mathbf{k}^{\prime})$ is invariant under shifts by $\Gamma^{*}$
lattice vectors.

Similarly,

(24)
$\begin{split}\pi_{\mathbf{q}^{\prime}}^{\dagger}\Lambda(\mathbf{q}^{\prime})^{\dagger}&=\sum_{\mathbf{k}^{\prime}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}^{\prime}}(\mathbf{q}^{\prime})^{\dagger}]_{mn}\ket{m,\mathbf{k}^{\prime}+\mathbf{q}^{\prime}}\bra{n,\mathbf{k}^{\prime}}\\\
&=\sum_{\mathbf{k}^{\prime}\in\mathcal{K}}\sum_{m,n\in\mathcal{N}}[\Lambda_{\mathbf{k}^{\prime}-\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})^{\dagger}]_{mn}\ket{m,\mathbf{k}^{\prime}}\bra{n,\mathbf{k}^{\prime}-\mathbf{q}^{\prime}}\\\
\end{split}$

where we have used the fact that
$\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})=\Lambda_{\mathbf{k}+\mathbf{G}}(\mathbf{q}^{\prime})$
for any $\mathbf{G}\in\Gamma^{*}$ (see Lemma A.1) to perform the change of
variables $\mathbf{k}^{\prime}\mapsto\mathbf{k}^{\prime}-\mathbf{q}^{\prime}$.

Therefore, the energy of a state $\ket{\Psi}$ can be written

(25) $\begin{split}\braket{\Psi}{\hat{H}_{\rm
FBI}}{\Psi}=&\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\left|\mathrm{Tr}{\Big{(}\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}Q\Big{)}}\right|^{2}\\\
&+\frac{1}{4N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\|\Lambda(\mathbf{q}^{\prime})\|_{F}^{2}\\\
&-\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\mathrm{Tr}{\bigg{(}\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}Q^{\dagger}\pi_{\mathbf{q}^{\prime}}^{\dagger}\Lambda(\mathbf{q}^{\prime})^{\dagger}Q\bigg{)}}.\end{split}$

We now recall the following elementary identity, whose proof is in Appendix B.

######  Lemma 3.1.

Let $A,B$ be arbitrary matrices then we have the following equality

(26)
$\operatorname{Re}{\left(\mathrm{Tr}{\left(ABA^{\dagger}B^{\dagger}\right)}\right)}=\frac{1}{2}\mathrm{Tr}{\Big{(}AA^{\dagger}B^{\dagger}B+A^{\dagger}ABB^{\dagger}\Big{)}}-\frac{1}{2}\|[A,B]\|_{F}^{2}.$

Since the energy of a state is purely real, we can set
$A=\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}$ and
$B=Q^{{\dagger}}$ and use Lemma 3.1 to rewrite Eq. 25 as

(27) $\begin{split}\langle&\Psi|\hat{H}_{\rm FBI}|\Psi\rangle\\\
&=\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left|\mathrm{Tr}{\Big{(}\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}Q\Big{)}}\right|^{2}\\\
&+\frac{1}{4N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\|\Lambda(\mathbf{q}^{\prime})\|_{F}^{2}\\\
&-\frac{1}{2N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\mathrm{Tr}{\bigg{(}\pi_{\mathbf{q}^{\prime}}^{\dagger}\Lambda(\mathbf{q}^{\prime})^{\dagger}\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}QQ^{\dagger}+\Lambda(\mathbf{q}^{\prime})\Lambda(\mathbf{q}^{\prime})^{\dagger}Q^{\dagger}Q\bigg{)}}\\\
&+\frac{1}{2N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\|[Q,\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}]\|_{F}^{2}.\end{split}$

But for Hartree-Fock states $QQ^{\dagger}=Q^{\dagger}Q=\frac{1}{4}I$, the
middle two terms cancel, leaving us with

(28) $\begin{split}\langle&\Psi|\hat{H}_{\rm FBI}|\Psi\rangle\\\
&=\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\left|\mathrm{Tr}{\Big{(}\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}Q\Big{)}}\right|^{2}+\frac{1}{2}\|[Q,\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}]\|_{F}^{2}\right)\\\
&=\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\left|\mathrm{Tr}{\Big{(}\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}\Big{(}P-\frac{1}{2}I\Big{)}\Big{)}}\right|^{2}+\frac{1}{2}\|[P,\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}]\|_{F}^{2}\right)\end{split}$

where in the last line we have used $Q=P-\frac{1}{2}I$ and that the identity
commutes with all matrices. Hence a Hartree-Fock state $\ket{\Psi}$ is a
ground state of $\hat{H}_{\rm FBI}$ if and only if the following two
conditions hold for all $\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}$

(29)
$\displaystyle\mathrm{Tr}{\Big{(}\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}\Big{(}P-\frac{1}{2}I\Big{)}\Big{)}}$
$\displaystyle=0,$ (30)
$\displaystyle[P,\Lambda(\mathbf{q}^{\prime})\pi_{\mathbf{q}^{\prime}}]$
$\displaystyle=0.$

While the trace condition Eq. 29 can be used to eliminate certain states as
ground states, the commutator condition Eq. 30 is a powerful tool for
characterizing possible Hartree-Fock ground states.

##  4\. Generalized Ferromagnetism in Flat-Band Interacting Hamiltonians

In a ferromagnet there are two ground states with equal energy:
$\ket{\uparrow\uparrow\uparrow\uparrow\uparrow\cdots}$ and
$\ket{\downarrow\downarrow\downarrow\downarrow\downarrow\cdots}$. For such
ground states, by knowing the value of the spin at a single site, we can
determine the value at all other sites. Similar to a ferromagnet, if
$\ket{\Psi}$ is a Hartree-Fock ground state of Eq. 6, then by knowing the
value of the 1-RDM at a single momenta there is an explicit method to compute
1-RDM at other momenta. We term this property “generalized ferromagnetism” and
formally have the following proposition:

######  Proposition 1 (Generalized Ferromagnetism).

Suppose that $P$ is the 1-RDM of a Hartree-Fock ground state of a frustration-
free FBI Hamiltonian Eq. 6 and the grid $\mathcal{K}$ satisfies 1. For any
$\mathbf{k}_{1},\mathbf{k}_{2}\in\mathcal{K}$ there exists a family of
invertible matrices
$\\{B_{\mathbf{k}}(\mathbf{q}^{\prime}):\mathbf{k}\in\mathcal{K},\mathbf{q}^{\prime}\in\mathbb{R}^{2}\\}$
so that

(31)
$P(\mathbf{k}_{1}+\mathbf{q}^{\prime},\mathbf{k}_{2}+\mathbf{q}^{\prime})=B_{\mathbf{k}_{1}}(\mathbf{q}^{\prime})^{-1}P(\mathbf{k}_{1},\mathbf{k}_{2})B_{\mathbf{k}_{2}}(\mathbf{q}^{\prime}).$

Hence, if the value of $P(\mathbf{k}_{1},\mathbf{k}_{2})$ is known, then
$P(\mathbf{k}_{1}+\mathbf{q}^{\prime},\mathbf{k}_{2}+\mathbf{q}^{\prime})$ is
determined for all $\mathbf{q}^{\prime}$.

An immediate implication of Proposition 1 is that the Hartree-Fock ground
states must be uniformly filled.

######  Corollary 1.

If $P$ is as in Proposition 1 then for all
$\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}$,
$\mathrm{Tr}{(P(\mathbf{k},\mathbf{k}))}=\mathrm{Tr}{(P(\mathbf{k}^{\prime},\mathbf{k}^{\prime}))}$.

###### Proof of Corollary 1.

Taking $\mathbf{k}_{1}=\mathbf{k}_{2}$ in the statement of Proposition 1 we
see that

(32)
$\begin{split}\mathrm{Tr}{(P(\mathbf{k}_{1}+\mathbf{q}^{\prime},\mathbf{k}_{1}+\mathbf{q}^{\prime}))}=\mathrm{Tr}{(B_{\mathbf{k}_{1}}(\mathbf{q}^{\prime})^{-1}P(\mathbf{k}_{1},\mathbf{k}_{1})B_{\mathbf{k}_{1}}(\mathbf{q}^{\prime}))}=\mathrm{Tr}{(P(\mathbf{k}_{1},\mathbf{k}_{1}))}.\end{split}$

Since the choice of $\mathbf{q}^{\prime}$ was arbitrary, this proves the
result. ∎

The proof of Proposition 1 relies on the following lemma proven in our
previous work:

######  Lemma 4.1 (Lemma 7.1 [4]).

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}\in\Gamma^{*}$ so that
$\Lambda_{\mathbf{k}}((\mathbf{k}^{\prime}-\mathbf{k})+\mathbf{G})$ has full
rank.

We now turn to prove Proposition 1.

###### Proof of Proposition 1.

Let us first fix some $\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}$ and some
$\mathbf{q}^{\prime}\in\mathcal{K}+\Gamma^{*}$. We can then decompose
$\mathbf{q}^{\prime}=\mathbf{q}+\mathbf{G}$ where $\mathbf{q}\in\mathcal{K}$
and $\mathbf{G}\in\Gamma^{*}$. For this choice of $\mathbf{q}+\mathbf{G}$, Eq.
30 implies that for all $\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}$

(33)
$P(\mathbf{k},\mathbf{k}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(\mathbf{q}+\mathbf{G})-\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})P(\mathbf{k}+\mathbf{q},\mathbf{k}^{\prime}+\mathbf{q})=0.$

If $\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ is invertible then we have
that

(34)
$P(\mathbf{k}+\mathbf{q},\mathbf{k}^{\prime}+\mathbf{q})=\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})^{-1}P(\mathbf{k},\mathbf{k}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(\mathbf{q}+\mathbf{G}).$

While we cannot assume that $\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ is
invertible in general, if we choose $\mathbf{q}$ sufficiently small, then by
Lemma 4.1 there exists a $\mathbf{G}$ so that
$\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$ must be invertible.

By 1, we can construct a sequence $\\{\mathbf{k}_{i}\\}_{i=1}^{L}$ which
connects $\mathbf{k}$ to $\mathbf{k}+\mathbf{q}$ so that
$\|\Pi(\mathbf{k}_{i+1})-\Pi(\mathbf{k}_{i})\|<1$ for $i=1,\cdots,L-1$. Along
this path, there always exists a $\mathbf{G}_{i}$ so that the form factor
matrix
$\Lambda_{\mathbf{k}_{i}}(\mathbf{k}_{i+1}-\mathbf{k}_{i}+\mathbf{G}_{i})$ is
invertible. Therefore, we can prove the proposition by recursively applying
Eq. 34.

To make the induction step more clear, we define the momentum difference
$\Delta:=\mathbf{k}^{\prime}-\mathbf{k}$ and the sequence of differences
$\delta_{i}=\mathbf{k}_{i+1}-\mathbf{k}_{i}$. By Lemma 4.1, there exists a
$\mathbf{G}_{1}\in\Gamma^{*}$ so that
$\Lambda_{\mathbf{k}_{1}}(\delta_{1}+\mathbf{G}_{1})$ is invertible.
Therefore, by Eq. 34

(35)
$\begin{split}P(\mathbf{k}_{2},\mathbf{k}_{2}+\Delta)&=P(\mathbf{k}_{1}+\delta_{1},\mathbf{k}_{1}+\Delta+\delta_{1})\\\
&=\Lambda_{\mathbf{k}_{1}}(\delta_{1}+\mathbf{G}_{1})^{-1}P(\mathbf{k}_{1},\mathbf{k}_{1}+\Delta)\Lambda_{\mathbf{k}_{1}+\Delta}(\delta_{1}+\mathbf{G}_{1}).\end{split}$

We can then apply the same argument to the next momentum in the sequence
$\mathbf{k}_{3}$

(36)
$\begin{split}P(\mathbf{k}_{3},\mathbf{k}_{3}+\Delta)&=P(\mathbf{k}_{2}+\delta_{2},\mathbf{k}_{2}+\Delta+\delta_{2})\\\
&=\Lambda_{\mathbf{k}_{2}}(\delta_{2}+\mathbf{G}_{2})^{-1}P(\mathbf{k}_{2},\mathbf{k}_{2}+\Delta)\Lambda_{\mathbf{k}_{2}+\Delta}(\delta_{2}+\mathbf{G}_{2}).\end{split}$

Thus, by repeatedly applying Eq. 34 along the path
$\\{\mathbf{k}_{i}\\}_{i=1}^{L}$ we conclude

(37)
$P(\mathbf{k}_{L},\mathbf{k}_{L}+\Delta)=\left(\prod_{i=1}^{L-1}\Lambda_{\mathbf{k}_{i}}(\delta_{i}+\mathbf{G}_{i})\right)^{-1}P(\mathbf{k}_{1},\mathbf{k}_{1}+\Delta)\left(\prod_{i=1}^{L-1}\Lambda_{\mathbf{k}_{i}+\Delta}(\delta_{i}+\mathbf{G}_{i})\right).$

Therefore, since $\mathbf{k}_{1}=\mathbf{k}$,
$\Delta=\mathbf{k}^{\prime}-\mathbf{k}$, and
$\mathbf{k}_{L}=\mathbf{k}+\mathbf{q}$ this proves the proposition. ∎

##  5\. The Sylvester Equation and Translation Breaking

While we considered Eq. 33 with $\mathbf{q}\neq\bm{0}$ in Section 4, let us
now consider the case when $\mathbf{q}=\bm{0}$. By choosing
$\mathbf{q}=\bm{0}$, from Eq. 33 we see that to be a ground state
$P(\mathbf{k},\mathbf{k}^{\prime})$ must satisfy

(38)
$P(\mathbf{k},\mathbf{k}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(\mathbf{G})-\Lambda_{\mathbf{k}}(\mathbf{G})P(\mathbf{k},\mathbf{k}^{\prime})=0\quad\forall\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K},\leavevmode\nobreak\
\forall\mathbf{G}\in\Gamma^{*}.$

The equation $XA-BX=C$ for some fixed $A,B,C$ is known as the Sylvester
equation. A simple condition ensures that this equation has a unique solution
(see [6, Theorem VII.2.1] with $C=0$).

######  Lemma 5.1.

Consider the Sylvester-type equation with $A,B\in\mathbb{C}^{N\times N}$,

(39) $XA-BX=0\text{ for }X\in\mathbb{C}^{N\times N}.$

The value $X=0$ is always a solution to this equation. Furthermore, $X=0$ is
the unique solution if the eigenvalues of $A$ and $B$ are disjoint.

Combining Eqs. 38 and 5.1 we have the following condition to rule out
translation symmetry breaking.

######  Proposition 2 (No translation symmetry breaking).

Suppose that $P$ is the 1-RDM of a Hartree-Fock ground state of a frustration-
free FBI Hamiltonian Eq. 6. For any
$\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}$, if there exists a
$\mathbf{G}\in\Gamma^{*}$ so that the eigenvalues of
$\Lambda_{\mathbf{k}}(\mathbf{G})$ and
$\Lambda_{\mathbf{k}^{\prime}}(\mathbf{G})$ are disjoint, then
$P(\mathbf{k},\mathbf{k}^{\prime})=0$.

If the eigenvalues of $A$ and $B$ in Eq. 39 are not disjoint then there can be
many possible choices for $X$ which solve $XA-BX=0$. While characterizing all
such possible solutions analytically is difficult in the general case, for the
special case of TBG-2 there is additional structure which allow us to
precisely determine the structure of the degeneracy in solutions to Eq. 38.

We begin by proving that no translation breaking occurs in spinless,
valleyless TBG-2 in Section 5.1 and then characterize the degeneracy in
valleyful and spinful models in Section 5.2. Finally in Section 5.3, we
outline a numerical algorithm which can be used to determine a spanning set
for the vector space of solutions to Eq. 38 and prove its correctness.

###  5.1. Application to Spinless, Valleyless TBG

We begin by converting Eq. 38 from an equation depending on
$\mathbf{G}\in\Gamma^{*}$ to an equation depending on $\mathbf{r}\in\Omega$ by
taking the Fourier transform.

First, we define
$\rho_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{r})\in\mathbb{C}^{(2M)\times(2M)}$
as follows

(40)
$[\rho_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{r})]_{m,n}:=\braket{u_{m\mathbf{k}}(\mathbf{r}),u_{n\mathbf{k}}(\mathbf{r})}.$

Since the periodic Bloch functions, $u_{n\mathbf{k}}(\mathbf{r})$, are
periodic with respect to $\Gamma$,
$\rho_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{r})$ is also periodic with
respect to $\Gamma$; that is,
$\rho_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{r}+\mathbf{R})=\rho_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{r})$
for all $\mathbf{R}\in\Gamma$. A straightforward calculation shows that that
the function $\mathbf{G}\mapsto\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G})$
coincides with the Fourier coefficients of
$\rho_{\mathbf{k},\mathbf{k}+\mathbf{q}}(\mathbf{r})$ when considered as a
periodic function on the unit cell $\Omega$. More specifically,

(41)
$\rho_{\mathbf{k},\mathbf{k}+\mathbf{q}}(\mathbf{r})=\sum_{\mathbf{G}\in\Gamma^{*}}e^{-i\mathbf{G}\cdot\mathbf{r}}\Lambda_{\mathbf{k}}(\mathbf{q}+\mathbf{G}).$

Therefore, by multiplying by $e^{i\mathbf{G}\cdot\mathbf{r}}$ and summing over
$\mathbf{G}\in\Gamma^{*}$ we can transform Eq. 38 to real space as follows

(42)
$P(\mathbf{k},\mathbf{k}^{\prime})\rho_{\mathbf{k}^{\prime},\mathbf{k}^{\prime}}(\mathbf{r})-\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r})P(\mathbf{k},\mathbf{k}^{\prime})=0,\qquad\forall\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K},\leavevmode\nobreak\
\forall\mathbf{r}\in\Omega.$

The Hamiltonian for TBG-2 in the chiral limit takes the form [23, 2]

(43) $H_{\mathbf{k}}=\begin{bmatrix}0&D_{\mathbf{k}}(\alpha)^{\dagger}\\\
D_{\mathbf{k}}(\alpha)&0\end{bmatrix}.$

Here, $D_{\mathbf{k}}(\alpha)=D(\alpha)+(k_{x}+ik_{y})I_{2\times 2}$ where
$D(\alpha)=\begin{bmatrix}2\bar{D}&\alpha V(\mathbf{r})\\\ \alpha
V(-\mathbf{r})&2\bar{D}\end{bmatrix}$ with $\overline{D}=-i\partial_{\bar{z}}$
and
$V(\mathbf{r})=\sum_{n=0}^{2}\omega^{n}e^{-i\mathbf{q}_{n}\cdot\mathbf{r}}$
with $\omega=e^{2\pi i/3}$.

The Hamiltonian $H_{\mathbf{k}}$ satisfies two symmetries $\mathcal{Q}$ and
$\mathcal{L}$ whose actions are defined as follows

(44) $\mathcal{Q}u(\mathbf{r}):=\begin{bmatrix}&&1&\\\ &&&1\\\ 1&&&\\\
&1&&\end{bmatrix}\overline{u(-\mathbf{r})}\qquad\mathcal{L}u(\mathbf{r}):=\begin{bmatrix}&1&&\\\
-1&&&\\\ &&&1\\\ &&-1&\end{bmatrix}\overline{u(-\mathbf{r})}\\\ $

Note that $\mathcal{Q}$ maps $\mathbf{k}$ to $\mathbf{k}$ and $\mathcal{L}$
maps $\mathbf{k}$ to $-\mathbf{k}$ and $[\mathcal{Q},\mathcal{L}]=0$.

Due to $\mathcal{Q}$ symmetry, the zero eigenvectors of $H_{\mathbf{k}}$ can
be written as $\begin{bmatrix}0&v(\mathbf{r})\end{bmatrix}^{\top}$ and
$\begin{bmatrix}\,\overline{v(-\mathbf{r})}&0\end{bmatrix}^{\top}$. Hence
after a basis rotation, $\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r})$ can be
written as follows

(45)
$\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r})=\begin{bmatrix}\|u_{\mathbf{k}}(\mathbf{r})\|^{2}&\\\
&\|u_{\mathbf{k}}(-\mathbf{r})\|^{2}\end{bmatrix}.$

Due to $\mathcal{L}$ symmetry and the fact that $[\mathcal{Q},\mathcal{L}]=0$,
we may also choose an appropriate basis change so that a state
$u_{\mathbf{0}}$ satisfying $D(\alpha)u_{\mathbf{0}}=0$ has the property

$\left\|u_{\mathbf{0}}(\mathbf{r})\right\|=\left\|\mathcal{L}u_{\mathbf{0}}(\mathbf{r})\right\|=\left\|u_{\mathbf{0}}(-\mathbf{r})\right\|,$

where, in this particular case, $\mathcal{L}$ just denotes the upper $2\times
2$ block of the symmetry defined in (44). With these observations, we can now
characterize all Hartree-Fock ground states of the FBI model for TBG-2.

######  Theorem 1.

For TBG-2 at half-filling without the valley or spin degrees of freedom, the
Hartree-Fock ground states take the form

(46)
$P(\mathbf{k},\mathbf{k}^{\prime})=\begin{cases}P_{0}&\mathbf{k}=\mathbf{k}^{\prime}\\\
0&\mathbf{k}\neq\mathbf{k}^{\prime}\end{cases}$

where $P_{0}$ the same for all $\mathbf{k}\in\mathcal{K}$ and equal to either
$\operatorname{diag}(1,0)$ or $\operatorname{diag}(0,1)$.

######  Remark 1.

The states found in Theorem 1 are the same ferromagnetic Slater determinant
states found in our previous work [4].

###### Proof.

By [3, Lemma $3.1$], we know that at a magic $\alpha$ the one-particle
Hamiltonian $H$ exhibits precisely two flat bands at energy zero

(47)
$\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})I_{2\times
2}),$

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{)}}.$

Since $\theta$ is an odd function, $F_{\mathbf{k}}(\mathbf{r})$ satisfies the
reflection symmetry

$F_{\mathbf{k}}(\mathbf{r})=F_{-\mathbf{k}}(-\mathbf{r}).$

The function $\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}$). This implies that

$\|u_{\mathbf{k}}(\mathbf{r})\|=\|F_{\mathbf{k}}(\mathbf{r})u_{\mathbf{0}}(\mathbf{r})\|=\|F_{-\mathbf{k}}(-\mathbf{r})u_{\mathbf{0}}(-\mathbf{r})\|=\|u_{-\mathbf{k}}(-\mathbf{r})\|.$

To show that $P(\mathbf{k},\mathbf{k}^{\prime})=0$ for a disjoint pair
$\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}$, it suffices to show by (42),
that there are $\mathbf{r},\mathbf{r}^{\prime}\in\Omega$ such that

${\|u_{\mathbf{k}}(\mathbf{r})\|\neq\|u_{\mathbf{k}^{\prime}}(\mathbf{r})\|}\text{
and
}\|u_{\mathbf{k}}(\mathbf{r}^{\prime})\|\neq\|u_{\mathbf{k}^{\prime}}(-\mathbf{r}^{\prime})\|.$

To satisfy the first constraint, we may choose the unique point
$\mathbf{r}=\frac{4\pi}{3\sqrt{3}}(k_{2},-k_{1})^{\top}$ at which
$u_{\mathbf{k}}(\mathbf{r})=0,$ then $u_{\mathbf{k}^{\prime}}(\mathbf{r})\neq
0.$

For this choice of $\mathbf{r},$ we also have

$u_{\mathbf{k}^{\prime}}(-\mathbf{r})\neq 0\text{ as long as
}\mathbf{k}^{\prime}\neq-\mathbf{k}.$

This implies that

$P(\mathbf{k},\mathbf{k}^{\prime})=0\text{ for
}\mathbf{k}^{\prime}\neq-\mathbf{k}.$

While we cannot directly conclude that $P(\mathbf{k},-\mathbf{k})=0$ for
$\mathbf{k}\notin\Gamma^{*}$, we may invoke Eq. 31, and use that we have
previously shown that $P(\mathbf{k},\mathbf{k}^{\prime})=0$ whenever
$\mathbf{k}\notin\\{\mathbf{k}^{\prime},-\mathbf{k}^{\prime}\\}$. In
particular, Eq. 31 implies that for any
$\mathbf{q}^{\prime}\in\mathbb{R}^{2}\setminus\Gamma^{*}$,

$\operatorname{rank}(P(\mathbf{k},-\mathbf{k}))=\operatorname{rank}(P(\mathbf{k}+\mathbf{q}^{\prime},-\mathbf{k}+\mathbf{q}^{\prime}))=0.$

This implies that

$P(\mathbf{k},-\mathbf{k})=0\text{ for all
}\mathbf{k}\in\mathcal{K}\setminus\\{\bm{0}\\}.$

Finally, for $\mathbf{k}=\mathbf{k}^{\prime}\neq\mathbf{0}$ it follows from
specializing in (45) to
$\mathbf{r}=\frac{4\pi}{3\sqrt{3}}(k_{2},-k_{1})^{\top}$, in which case
$\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r})=\operatorname{diag}(0,\|u_{\mathbf{k}}(-\mathbf{r})\|^{2})$
with $u_{\mathbf{k}}(-\mathbf{r})\neq 0.$ This implies that
$P(\mathbf{k},\mathbf{k})$ is diagonal. Since $P^{2}=P$ and
$P(\mathbf{k},\mathbf{k}^{\prime})=0$ whenever
$\mathbf{k}\neq\mathbf{k}^{\prime}$, we conclude that
$P(\mathbf{k},\mathbf{k})^{2}=P(\mathbf{k},\mathbf{k})$ and so its eigenvalues
are $0$ and $1$. At half-filling, by Corollary 1, we know
$\mathrm{Tr}{(P(\mathbf{k},\mathbf{k}))}=1$ and so
$P(\mathbf{k},\mathbf{k})=\operatorname{diag}(1,0)$ or
$P(\mathbf{k},\mathbf{k})=\operatorname{diag}(0,1)$ for all $\mathbf{k}$.

Finally, for $\mathbf{k}^{\prime}=\mathbf{k}=\mathbf{0}$, the condition (45)
does not imply any restrictions on $P(\mathbf{0},\mathbf{0}),$ since
$\rho_{\mathbf{0},\mathbf{0}}$ is proportional to the identity matrix.
However, the proof of Proposition 1 shows that one can choose diagonal
$B_{\mathbf{k}}(\mathbf{q}^{\prime})$. This implies by (45) that
$\mathbf{k}\mapsto P(\mathbf{k},\mathbf{k})$ is either
$\operatorname{diag}(1,0)$ or $\operatorname{diag}(0,1)$ for all $\mathbf{k}$.

To complete the proof, we must verify that Eq. 29 also vanishes. This follows
from the sum rule

(48)
$\sum_{\mathbf{k}}\mathrm{Tr}{\left(\Lambda_{\mathbf{k}}(\mathbf{G})\begin{bmatrix}1&\\\
&-1\end{bmatrix}\right)}=0\quad\forall\mathbf{G}\in\Gamma^{*}$

proven in our previous work [4, Lemma 4.3].

∎

###  5.2. Extension to Valleyful and Spinful Models

In the previous section, we showed that in single valley TBG-2 the only
possible Hartree-Fock ground states are the two ferromagnetic Slater
determinants by using the specific properties of the TBG-2 eigenfunctions.
While the final conclusion depended on the specific form of these
eigenfunctions all of the steps up to and including Eq. 42 were done in a
completely general setting (i.e. without reference to an underlying
Hamiltonian). Therefore, to translate these results to the more realistic
valleyful and spinful settings we simply need to substitute in appropriate
definition for $\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r})$.

Following the analysis in [8], for the chiral model of TBG-2 with valley, the
corresponding $\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r})$ takes the form:

(49) $\rho_{\mathbf{k},\mathbf{k}}^{\rm
V}(\mathbf{r})=\begin{bmatrix}\|u_{\mathbf{k}}(\mathbf{r})\|^{2}&&&\\\
&\|u_{\mathbf{k}}(-\mathbf{r})\|^{2}&&\\\
&&\|u_{\mathbf{k}}(-\mathbf{r})\|^{2}&\\\
&&&\|u_{\mathbf{k}}(\mathbf{r})\|^{2}\end{bmatrix}$

and when including valley and spin it takes the form

(50) $\rho_{\mathbf{k},\mathbf{k}}^{\rm
V+S}(\mathbf{r})=\begin{bmatrix}\|u_{\mathbf{k}}(\mathbf{r})\|^{2}&&&\\\
&\|u_{\mathbf{k}}(-\mathbf{r})\|^{2}&&\\\
&&\|u_{\mathbf{k}}(-\mathbf{r})\|^{2}&\\\
&&&\|u_{\mathbf{k}}(\mathbf{r})\|^{2}\end{bmatrix}\otimes I_{2\times 2}.$

We will now show that when including the valley degree of freedom, the ground
state exhibits a ${\rm U}(2)\times{\rm U}(2)$ degeneracy. A similar argument
shows that in the presence of the valley and spin degrees of freedom, the
ground state is unique up to a ${\rm U}(4)\times{\rm U}(4)$ degree of freedom.

######  Theorem 2.

At TBG-2 at half-filling with the valley degree of freedom but without spin,
the Hartree-Fock ground states take the form

(51)
$P(\mathbf{k},\mathbf{k}^{\prime})=\begin{cases}UP_{0}U^{\dagger}&\mathbf{k}=\mathbf{k}^{\prime}\\\
0&\mathbf{k}\neq\mathbf{k}^{\prime}\end{cases}$

where $U$ and $P_{0}$ are the same for all $\mathbf{k}\in\mathcal{K}$.
Furthermore, $U\in{\mathrm{U}}(2)\times{\mathrm{U}}(2)$ and $P_{0}$ is one of
the following matrices $\operatorname{diag}(1,0,0,1)$,
$\operatorname{diag}(1,1,0,0)$, and $\operatorname{diag}(0,1,1,0)$.

###### Proof.

Following the same reasoning used in the proof of Theorem 1, we know that
$P(\mathbf{k},\mathbf{k}^{\prime})=0$ for all
$\mathbf{k}\neq\mathbf{k}^{\prime}$. By Proposition 1, it suffices to choose a
single $\mathbf{k}\in\mathcal{K}$ to fix the ground state. Pick any
$\mathbf{k}_{*}\in\mathcal{K}$ so that
$\|u_{\mathbf{k}_{*}}(\mathbf{r})\|\neq\|u_{\mathbf{k}_{*}}(-\mathbf{r})\|$
for some $\mathbf{r}\in\Omega$. To make the algebra simpler, we define the
permutation matrix $\Pi$:

(52) $\Pi=\begin{bmatrix}1&&&\\\ &&&1\\\ &&1&\\\ &1&&\end{bmatrix}$

which satisfies $\Pi^{\dagger}=\Pi$ and $\Pi^{2}=I$. Applying $\Pi$ to both
sides of Eq. 45 we find that for $P(\mathbf{k}_{*},\mathbf{k}_{*})$ to be the
1-RDM of a ground state it must be that

(53) $\Big{(}\Pi
P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi\Big{)}\Big{(}\Pi\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}(\mathbf{r})\Pi\Big{)}-\Big{(}\Pi\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}(\mathbf{r})\Pi\Big{)}\Big{(}\Pi
P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi\Big{)}=0.$

We will derive a characterization for all possible solutions $\Pi
P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi$ of Eq. 53 at half-filling. Since $\Pi$ is
a unitary, this also characterizes all possible solutions
$P(\mathbf{k}_{*},\mathbf{k}_{*})$ satisfying Eq. 42 after a change of basis.

Notice that

(54) $\begin{split}\Pi\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}^{\rm
V}(\mathbf{r})\Pi&=\begin{bmatrix}\|u_{\mathbf{k}_{*}}(\mathbf{r})\|^{2}&&&\\\
&\|u_{\mathbf{k}_{*}}(\mathbf{r})\|^{2}&&\\\
&&\|u_{\mathbf{k}_{*}}(-\mathbf{r})\|^{2}&\\\
&&&\|u_{\mathbf{k}_{*}}(-\mathbf{r})\|^{2}\end{bmatrix}\\\\[4.30554pt]
&=\begin{bmatrix}\|u_{\mathbf{k}_{*}}(\mathbf{r})\|^{2}I_{2\times 2}&\\\
&\|u_{\mathbf{k}_{*}}(-\mathbf{r})\|^{2}I_{2\times
2}\end{bmatrix}.\end{split}$

So splitting $\Pi P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi$ into $2\times 2$ blocks

(55) $\Pi
P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi=:\begin{bmatrix}\tilde{P}_{11}&\tilde{P}_{12}\\\
\tilde{P}_{12}^{\dagger}&\tilde{P}_{22}\end{bmatrix},$

we can explicitly compute the commutator in terms of $\tilde{P}_{11}$,
$\tilde{P}_{22}$, $\tilde{P}_{12}$:

(56) $[\Pi
P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi,\Pi\rho_{\mathbf{k}_{*},\mathbf{k}_{*}}^{\rm
V}(\mathbf{r})\Pi]=(\|u_{\mathbf{k}_{*}}(-\mathbf{r})\|^{2}-\|u_{\mathbf{k}_{*}}(\mathbf{r})\|^{2})\begin{bmatrix}&\tilde{P}_{12}\\\
-\tilde{P}_{12}^{\dagger}&\end{bmatrix}.$

Since
$\|u_{\mathbf{k}_{*}}(-\mathbf{r})\|^{2}\neq\|u_{\mathbf{k}_{*}}(\mathbf{r})\|^{2}$
for some $\mathbf{r}$, for Eq. 45 to vanish it must be that
$\tilde{P}_{12}\equiv 0$. Hence any solution to Eq. 53 must be block diagonal.

Since $P(\mathbf{k}_{*},\mathbf{k}^{\prime})=0$ for all
$\mathbf{k}_{*}\neq\mathbf{k}^{\prime}$, $\Pi
P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi$ is an orthogonal projection. The set of
all transformations preserving this block structure while maintaining
orthogonality is given by $U\in{\rm U}(2)\times{\rm U}(2)$ where $U$ takes the
form

(57) $U=\begin{bmatrix}U_{1}&\\\ &U_{2}\end{bmatrix}\quad U_{1},U_{2}\in{\rm
U}(2).$

Due to Corollary 1, at half-filling, $\Pi P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi$
has rank 2 and hence all possible solutions to Eq. 53 must lie in one of three
orbits:

(58) $U\begin{bmatrix}1&&&\\\ &1&&\\\ &&0&\\\
&&&0\end{bmatrix}U^{\dagger}\qquad U\begin{bmatrix}1&&&\\\ &0&&\\\ &&1&\\\
&&&0\end{bmatrix}U^{\dagger}\qquad U\begin{bmatrix}0&&&\\\ &0&&\\\ &&1&\\\
&&&1\end{bmatrix}.U^{\dagger}$

Suppose we choose $\Pi P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi$ to be one of the
options in Eq. 58; we now will argue that $\Pi P(\mathbf{k},\mathbf{k})\Pi=\Pi
P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi$ for all $\mathbf{k}\in\mathcal{K}$. Due
to $\mathcal{Q}$ symmetry, the pair product
$\Pi\rho_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{V}(\mathbf{r})\Pi$ has the same
form as in Eq. 54 for all $\mathbf{k},\mathbf{k}+\mathbf{q}$ [8]. Therefore,
the pair product commutes with $\Pi P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi$ and
so following the proof of Proposition 1, we can conclude for a ground state
$\Pi P(\mathbf{k},\mathbf{k})\Pi=\Pi P(\mathbf{k}_{*},\mathbf{k}_{*})\Pi$ for
all $\mathbf{k}\in\mathcal{K}$.

To complete the proof, we must also verify that that Eq. 29 vanishes. Since
$[U,\Pi\rho_{\mathbf{k},\mathbf{k}}^{V}(\mathbf{r})\Pi]=0$ and
$\Pi\rho_{\mathbf{k},\mathbf{k}}^{V}(\mathbf{r})\Pi=I_{2\times
2}\otimes\rho_{\mathbf{k},\mathbf{k}}(\mathbf{r})$, the sum rule Eq. 48
implies Eq. 29 vanishes. The statement of the theorem follows reordering the
basis by the permutation $\Pi$. ∎

######  Remark 2.

We shall compare our result with statements in the physics literature.

  1. (1)

In the physics literature, the structure Eq. 53 is often expressed as
$[\sigma_{z}\tau_{z},P(\mathbf{k},\mathbf{k})]=0$ where $\sigma_{z}$ denotes
the Pauli Z-matrix acting on the sublattice degree of freedom and $\tau_{z}$
denotes the Pauli Z-matrix acting on the valley degree of freedom. The
operator $\sigma_{z}\tau_{z}$ is called the Chern number operator.

  2. (2)

In the physics literature, the ${\rm U}(2)\times{\rm U}(2)$ (or ${\rm
U}(4)\times{\rm U}(4)$) degree of freedom is often stated as a hidden symmetry
of the FBI Hamiltonian. The two perspectives are equivalent. This is because
if $\ket{\Psi}$ is a ground state of $\hat{H}_{\rm FBI}$ and $\mathcal{G}$ is
a symmetry satisfying $[\mathcal{G},\hat{H}_{FBI}]=0$, then
$\mathcal{G}\ket{\Psi}$ is also a ground state.

  3. (3)

Since $U$ is block diagonal, the orbits of $\operatorname{diag}(1,0,0,1)$ and
$\operatorname{diag}(0,1,1,0)$ only contain one element. These two states are
referred to as the quantum hall (QH) states as they occupy the positive
eigenspace and negative eigenspace of the Chern number operator
$\sigma_{z}\tau_{z}$. In contrast, the third orbit contains a manifold of
states which include states which occupy both valleys; such states are
referred to as the intervalley coherent (IVC) states.

Using the same proof strategy, we can easily extend Theorem 2 to include spin:

######  Theorem 3.

At TBG-2 at half-filling with the valley and spin degrees of freedom, the
Hartree-Fock ground states take the form

(59)
$P(\mathbf{k},\mathbf{k}^{\prime})=\begin{cases}UP_{0}U^{\dagger}&\mathbf{k}=\mathbf{k}^{\prime}\\\
0&\mathbf{k}\neq\mathbf{k}^{\prime}\end{cases}$

where $U$ and $P_{0}$ are the same for all $\mathbf{k}\in\mathcal{K}$.
Furthermore, $U\in{\mathrm{U}}(4)\times{\mathrm{U}}(4)$ and $P_{0}$ is equal
to one of the following matrices $\operatorname{diag}(1,0,0,1,1,0,0,1)$,
$\operatorname{diag}(1,1,0,1,1,0,0,0)$,
$\operatorname{diag}(1,1,1,1,0,0,0,0)$,
$\operatorname{diag}(0,0,1,0,0,1,1,1)$, or
$\operatorname{diag}(0,1,1,0,0,1,1,0)$.

###  5.3. An Algorithm for Determining Translation Symmetry Breaking

In cases when an analytical form of the form factor is not known, Eq. 38 can
be used to design an algorithm for identifying when translation symmetry
breaking can occur. From the vectorization identity, we can convert Eq. 38
into a vector equation as follows

(60) $\left(\Lambda_{\mathbf{k}^{\prime}}(\mathbf{G})^{\top}\otimes
I_{(2M)\times(2M)}-I_{(2M)\times(2M)}\otimes\Lambda_{\mathbf{k}}(\mathbf{G})\right)\operatorname{vec}{(P(\mathbf{k},\mathbf{k}^{\prime}))}=0\quad\forall\mathbf{G}\in\Gamma^{*},$

where $\operatorname{vec}{(\cdot)}$ transforms a $(2M)\times(2M)$ matrix into
a vector with $(2M)^{2}$ entries by stacking the columns. Define the
$(2M)^{2}\times(2M)^{2}$ matrix
$M_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{G})$:

(61)
$M_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{G}):=\Lambda_{\mathbf{k}^{\prime}}(\mathbf{G})^{\top}\otimes
I_{(2M)\times(2M)}-I_{(2M)\times(2M)}\otimes\Lambda_{\mathbf{k}}(\mathbf{G}).$

Then $P(\mathbf{k},\mathbf{k}^{\prime})$ satisfies Eq. 38 if and only if its
vectorization lies in the kernel
$M_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{G})$ for all
$\mathbf{G}\in\Gamma^{*}$. We can further reduce our computational effort by
squaring $M_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{G})$ and summing over
$\mathbf{G}$:

(62)
$M_{\mathbf{k},\mathbf{k}^{\prime}}:=\sum_{\mathbf{G}\in\Gamma^{*}}M_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{G})^{\dagger}M_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{G}).$

The following proposition is now immediate

######  Proposition 3.

The matrix $M_{\mathbf{k},\mathbf{k}^{\prime}}$ is positive semidefinite and
its kernel is exactly the intersection of the kernels of all the
$M_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{G})$.

Since the form factors are the Fourier coefficients of a smooth function (Eq.
41) we have that

(63)
$\sum_{\\{\mathbf{G}\in\Gamma^{*}:|\mathbf{G}|>L\\}}\|\Lambda_{\mathbf{k}^{\prime}}(\mathbf{G})\|_{F}^{2}\rightarrow
0\quad\text{as}\quad L\rightarrow\infty$

Hence, in practice we can truncate the sum over $\mathbf{G}$ in Eq. 62 and by
the Weyl bound we can find a spanning set for all
$P(\mathbf{k},\mathbf{k}^{\prime})$ lying in the kernel of
$M_{\mathbf{k},\mathbf{k}^{\prime}}$ by taking a few of the smallest
eigenmodes of the truncated matrix.

######  Remark 3.

In a discretized setting, Proposition 3 also gives us a computationally
efficient method for verifying the “no common invariant subspace” condition in
[4, Theorem 4].

## References

  * [1] Simon Becker, Mark Embree, Jens Wittsten, and Maciej Zworski.  Spectral characterization of magic angles in twisted bilayer graphene.  Phys. Rev. B, 103(16):165113, 2021.
  * [2] Simon Becker, Mark Embree, Jens Wittsten, and Maciej Zworski.  Mathematics of magic angles in a model of twisted bilayer graphene.  Prob. Math. Phys., 3(1):69–103, 2022.
  * [3] Simon Becker, Tristan Humbert, and Maciej Zworski.  Fine structure of flat bands in a chiral model of magic angles.  arXiv:2208.01628, 2022.
  * [4] Simon Becker, Lin Lin, and Kevin D. Stubbs.  Exact ground state of interacting electrons in magic angle graphene, December 2023.
  * [5] B. Andrei Bernevig, Zhi-Da Song, Nicolas Regnault, and Biao Lian.  Twisted bilayer graphene. III. Interacting Hamiltonian and exact symmetries.  Phys. Rev. B, 103(20):205413, 2021.
  * [6] Rajendra Bhatia.  Matrix Analysis, volume 169 of Graduate Texts in Mathematics.  Springer New York, New York, NY, 1997.
  * [7] R. Bistritzer and A. H. MacDonald.  Moire bands in twisted double-layer graphene.  Proceedings of the National Academy of Sciences, 108(30):12233–12237, 2021.
  * [8] Nick Bultinck, Eslam Khalaf, Shang Liu, Shubhayu Chatterjee, Ashvin Vishwanath, and Michael P. Zaletel.  Ground State and Hidden Symmetry of Magic-Angle Graphene at Even Integer Filling.  Phys. Rev. X, 10(3):031034, 2020.
  * [9] Éric Cancès, Louis Garrigue, and David Gontier.  Simple derivation of moiré-scale continuous models for twisted bilayer graphene.  Phys. Rev. B, 107:155403, Apr 2023.
  * [10] Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero.  Unconventional superconductivity in magic-angle graphene superlattices.  Nature, 556(7699):43–50, 2018.
  * [11] Shubhayu Chatterjee, Nick Bultinck, and Michael P. Zaletel.  Symmetry breaking and skyrmionic transport in twisted bilayer graphene.  Phys. Rev. B, 101(16), apr 2020.
  * [12] Ipsita Das, Xiaobo Lu, Jonah Herzog-Arbeitman, Zhi-Da Song, Kenji Watanabe, Takashi Taniguchi, B. Andrei Bernevig, and Dmitri K. Efetov.  Symmetry-broken chern insulators and rashba-like landau-level crossings in magic-angle bilayer graphene.  Nat. Phys., 17(6):710–714, mar 2021.
  * [13] Fabian M Faulstich, Kevin D Stubbs, Qinyi Zhu, Tomohiro Soejima, Rohit Dilip, Huanchen Zhai, Raehyun Kim, Michael P Zaletel, Garnet Kin-Lic Chan, and Lin Lin.  Interacting models for twisted bilayer graphene: A quantum chemistry approach.  Phys. Rev. B, 107(23):235123, 2023.
  * [14] Yuhang Jiang, Xinyuan Lai, Kenji Watanabe, Takashi Taniguchi, Kristjan Haule, Jinhai Mao, and Eva Y. Andrei.  Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene.  Nature, 573(7772):91–95, jul 2019.
  * [15] Biao Lian, Zhi-Da Song, Nicolas Regnault, Dmitri K. Efetov, Ali Yazdani, and B. Andrei Bernevig.  Twisted bilayer graphene. IV. Exact insulator ground states and phase diagram.  Phys. Rev. B, 103(20):205414, 2021.
  * [16] Shang Liu, Eslam Khalaf, Jong Yeon Lee, and Ashvin Vishwanath.  Nematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutrality.  Phys. Rev. Res., 3(1), jan 2021.
  * [17] Shang Liu, Eslam Khalaf, Jong Yeon Lee, and Ashvin Vishwanath.  Nematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutrality.  Phys. Rev. Res., 3(1):013033, 2021.
  * [18] Kevin P Nuckolls, Ryan L Lee, Myungchul Oh, Dillon Wong, Tomohiro Soejima, Jung Pyo Hong, Dumitru Călugăru, Jonah Herzog-Arbeitman, B Andrei Bernevig, Kenji Watanabe, et al.  Quantum textures of the many-body wavefunctions in magic-angle graphene.  Nature, 620(7974):525–532, 2023.
  * [19] Pawel Potasz, Ming Xie, and Allan H. MacDonald.  Exact Diagonalization for Magic-Angle Twisted Bilayer Graphene.  Phys. Rev. Lett., 127(14):147203, 2021.
  * [20] Yu Saito, Jingyuan Ge, Louk Rademaker, Kenji Watanabe, Takashi Taniguchi, Dmitry A. Abanin, and Andrea F. Young.  Hofstadter subband ferromagnetism and symmetry-broken chern insulators in twisted bilayer graphene.  Nat. Phys., 17(4):478–481, jan 2021.
  * [21] Yu Saito, Jingyuan Ge, Kenji Watanabe, Takashi Taniguchi, and Andrea F. Young.  Independent superconductors and correlated insulators in twisted bilayer graphene.  Nature Physics, 16(9):926–930, June 2020.
  * [22] Tomohiro Soejima, Daniel E. Parker, Nick Bultinck, Johannes Hauschild, and Michael P. Zaletel.  Efficient simulation of moire materials using the density matrix renormalization group.  Phys. Rev. B, 102(20):205111, 2020.
  * [23] Grigory Tarnopolsky, Alex Jura Kruchkov, and Ashvin Vishwanath.  Origin of magic angles in twisted bilayer graphene.  Phys. Rev. Let., 122(10):106405, 2019.
  * [24] Chun-Chih Tseng, Xuetao Ma, Zhaoyu Liu, Kenji Watanabe, Takashi Taniguchi, Jiun-Haw Chu, and Matthew Yankowitz.  Anomalous hall effect at half filling in twisted bilayer graphene.  Nature Physics, 18(9):1038–1042, 2022.
  * [25] Glenn Wagner, Yves H Kwan, Nick Bultinck, Steven H Simon, and SA Parameswaran.  Global phase diagram of the normal state of twisted bilayer graphene.  Physical Review Letters, 128(15):156401, 2022.
  * [26] Alexander B. Watson, Tianyu Kong, Allan H. MacDonald, and Mitchell Luskin.  Bistritzer-MacDonald dynamics in twisted bilayer graphene.  Journal of Mathematical Physics, 64(3):031502, March 2023.
  * [27] Alexander B Watson and Mitchell Luskin.  Existence of the first magic angle for the chiral model of bilayer graphene.  J. Math. Phys., 62(9):091502, 2021.
  * [28] Fengcheng Wu and Sankar Das Sarma.  Collective excitations of quantum anomalous hall ferromagnets in twisted bilayer graphene.  Phys. Rev. Lett., 124(4), jan 2020.
  * [29] Fang Xie, Jian Kang, B Andrei Bernevig, Oskar Vafek, and Nicolas Regnault.  Phase diagram of twisted bilayer graphene at filling factor $\nu$=$\pm$3\.  Physical Review B, 107(7):075156, 2023.
  * [30] Ming Xie and A. H. MacDonald.  Nature of the Correlated Insulator States in Twisted Bilayer Graphene.  Phys. Rev. Lett., 124(9):097601, 2020.
  * [31] Yonglong Xie, Biao Lian, Berthold Jäck, Xiaomeng Liu, Cheng-Li Chiu, Kenji Watanabe, Takashi Taniguchi, B Andrei Bernevig, and Ali Yazdani.  Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene.  Nature, 572(7767):101–105, 2019.
  * [32] Geng-Dong Zhou, Yi-Jie Wang, Ninghua Tong, and Zhi-Da Song.  Kondo phase in twisted bilayer graphene.  Physical Review B, 109(4):045419, 2024.

##  Appendix A Calculation of the Hartree-Fock Energy

To perform this calculation, we will recall a few properties of the form
factor which will play a role in the following calculations (see, e.g. [4,
Lemma 4.2] for proofs):

######  Lemma A.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^{*}$

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

Throughout this calculation, we will suppress the summation indexes for
$\mathbf{q}^{\prime},\mathbf{k},\mathbf{k}^{\prime},m,n,m^{\prime},n^{\prime}$.
Expanding the terms in $\widehat{\rho}(\mathbf{q}^{\prime})$ and performing
normal ordering gives

(66) $\begin{split}N_{\mathbf{k}}|\Omega|\hat{H}_{\rm
FBI}=&\sum_{\mathbf{q}^{\prime}}\sum_{\mathbf{k},\mathbf{k}^{\prime}}\sum_{m,n,m^{\prime},n^{\prime}}\hat{V}(\mathbf{q}^{\prime})[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}[\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})]_{m^{\prime}n^{\prime}}\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{m^{\prime}\mathbf{k}^{\prime}}^{\dagger}\hat{f}_{n^{\prime}(\mathbf{k}^{\prime}-\mathbf{q}^{\prime})}\hat{f}_{n(\mathbf{k}+\mathbf{q}^{\prime})}\\\\[4.30554pt]
&+\sum_{\mathbf{q}^{\prime}}\sum_{\mathbf{k}^{\prime}}\sum_{m,n}\hat{V}(\mathbf{q}^{\prime})[\Lambda_{\mathbf{k}+\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})]_{mn}\hat{f}_{m\mathbf{k}^{\prime}}^{\dagger}\hat{f}_{n\mathbf{k}^{\prime}}\\\
&-\frac{1}{2}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{(\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}}\sum_{m,n}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n\mathbf{k}}\right)\\\
&-\frac{1}{2}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}}\mathrm{Tr}{(\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}^{\prime}}\sum_{m^{\prime},n^{\prime}}[\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})]_{mn}\hat{f}_{m^{\prime}\mathbf{k}^{\prime}}^{\dagger}\hat{f}_{n^{\prime}\mathbf{k}^{\prime}}\right)\\\
&+\frac{1}{4}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}}\mathrm{Tr}{(\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{(\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime}))}\right).\end{split}$

By Wick’s theorem, we can evaluate the energy of any Hartree-Fock state,
$\ket{\Psi}$ in terms of the one-body reduced density matrix
$[P(\mathbf{k}^{\prime},\mathbf{k})]_{nm}:=\braket{\Phi}{\hat{f}_{m\mathbf{k}}^{\dagger}\hat{f}_{n\mathbf{k}^{\prime}}}{\Phi}$.
In particular, we have

(67) $\begin{split}N_{\mathbf{k}}|\Omega|&\langle\Psi|\hat{H}_{\rm
FBI}|\Psi\rangle\\\
=&\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}}\sum_{m,n,m^{\prime},n^{\prime}}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}[\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})]_{m^{\prime}n^{\prime}}[P(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k})]_{nm}[P(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k}^{\prime})]_{n^{\prime}m^{\prime}}\\\
&-\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}}\sum_{m,n,m^{\prime},n^{\prime}}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}[\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})]_{m^{\prime}n^{\prime}}[P(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k})]_{n^{\prime}m}[P(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k}^{\prime})]_{nm^{\prime}}\\\
&-\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k}^{\prime}}\sum_{m,n}[\Lambda_{\mathbf{k}^{\prime}+\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})]_{mn}[P(\mathbf{k}^{\prime},\mathbf{k}^{\prime})]_{nm}\\\
&-\frac{1}{2}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{(\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}}\sum_{m,n}[\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})]_{mn}[P(\mathbf{k},\mathbf{k})]_{nm}\right)\\\
&-\frac{1}{2}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}}\mathrm{Tr}{(\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}^{\prime}}\sum_{m^{\prime},n^{\prime}}[\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})]_{mn}[P(\mathbf{k}^{\prime},\mathbf{k}^{\prime})]_{n^{\prime}m^{\prime}}\right)\\\
&+\frac{1}{4}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}}\mathrm{Tr}{(\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{(\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime}))}\right)\end{split}$

which can cleanly be expressed in terms of traces as

(68) $\begin{split}N_{\mathbf{k}}|\Omega|&\langle\Psi|\hat{H}_{\rm
FBI}|\Psi\rangle\\\
=&\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})P(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})P(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k}^{\prime})\Big{)}}\\\
&-\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})P(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})P(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\\\
&-\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}^{\prime}+\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})P(\mathbf{k}^{\prime},\mathbf{k}^{\prime})\Big{)}}\\\
&-\frac{1}{2}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{(\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}}\mathrm{Tr}{(\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})P(\mathbf{k},\mathbf{k}))}\right)\\\
&-\frac{1}{2}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}}\mathrm{Tr}{(\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{(\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})P(\mathbf{k}^{\prime},\mathbf{k}^{\prime}))}\right)\\\
&+\frac{1}{4}\sum_{\mathbf{q}^{\prime}}\delta_{\mathbf{q}^{\prime}\in\Gamma^{*}}\hat{V}(\mathbf{q}^{\prime})\left(\sum_{\mathbf{k}}\mathrm{Tr}{(\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime}))}\right)\left(\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{(\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime}))}\right).\end{split}$

Recalling the definition of the matrix $Q(\mathbf{k},\mathbf{k}^{\prime})$
(Eq. 15), we can combine the first term and the last three terms to get

(69) $\begin{split}N_{\mathbf{k}}|\Omega|\braket{\Psi}{\hat{H}_{\rm
FBI}}{\Psi}=&\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})Q(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})Q(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k}^{\prime})\Big{)}}\\\
&-\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})P(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})P(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\\\
&-\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k}^{\prime}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}^{\prime}+\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})P(\mathbf{k}^{\prime},\mathbf{k}^{\prime})\Big{)}}.\end{split}$

Using Eq. 64 and performing a change of variables
$\mathbf{k}^{\prime}\mapsto\mathbf{k}^{\prime}-\mathbf{q}^{\prime}$ (this
change of variables is valid due to Eq. 65), we can combine the last two terms
in the above equation to get

(70) $\begin{split}N_{\mathbf{k}}|\Omega|\langle\Psi&|\hat{H}_{\rm
FBI}|\Psi\rangle\\\
=&\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})Q(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})Q(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k}^{\prime})\Big{)}}\\\
&-\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})Q(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k}^{\prime})\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})Q(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\\\
&+\frac{1}{4}\sum_{\mathbf{q}^{\prime}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})^{\dagger}\Big{)}}.\end{split}$

Finally, using the identity
$\Lambda_{\mathbf{k}^{\prime}}(-\mathbf{q}^{\prime})=\Lambda_{\mathbf{k}^{\prime}-\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})^{\dagger}$
and performing the change of variables in the first sum
$\mathbf{k}^{\prime}\mapsto\mathbf{k}^{\prime}+\mathbf{q}^{\prime}$, we can
finally write the energy of $\ket{\Psi}$ as

(71) $\begin{split}\langle\Psi&|\hat{H}_{\rm FBI}|\Psi\rangle\\\
=&\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\left|\sum_{\mathbf{k}\in\mathcal{K}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})Q(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\right|^{2}\\\
&+\frac{1}{4N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k}\in\mathcal{K}}\|\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})\|_{F}^{2}\\\
&-\frac{1}{N_{\mathbf{k}}|\Omega|}\sum_{\mathbf{q}^{\prime}\in\Gamma^{*}+\mathcal{K}}\hat{V}(\mathbf{q}^{\prime})\sum_{\mathbf{k},\mathbf{k}^{\prime}\in\mathcal{K}}\mathrm{Tr}{\Big{(}\Lambda_{\mathbf{k}}(\mathbf{q}^{\prime})Q(\mathbf{k}+\mathbf{q}^{\prime},\mathbf{k}^{\prime})\Lambda_{\mathbf{k}^{\prime}-\mathbf{q}^{\prime}}(\mathbf{q}^{\prime})^{\dagger}Q(\mathbf{k}^{\prime}-\mathbf{q}^{\prime},\mathbf{k})\Big{)}}\end{split}$

which completes the calculation.

##  Appendix B Proof of Lemma 3.1

First, we calculate

(72)
$\|[A,B]\|_{F}^{2}=\mathrm{Tr}{([A,B]([A,B])^{\dagger})}=-\mathrm{Tr}{([A,B]([A^{\dagger},B^{\dagger}]))}.$

Now we expand the product of commutators

(73) $\begin{split}[A,B][A^{\dagger},B^{\dagger}]&=(AB-
BA)(A^{\dagger}B^{\dagger}-B^{\dagger}A^{\dagger})\\\
&=ABA^{\dagger}B^{\dagger}-ABB^{\dagger}A^{\dagger}-BAA^{\dagger}B^{\dagger}+BAB^{\dagger}A^{\dagger}\\\
&=\left(ABA^{\dagger}B^{\dagger}+(ABA^{\dagger}B^{\dagger})^{\dagger}\right)-\left(ABB^{\dagger}A^{\dagger}+BAA^{\dagger}B^{\dagger}\right).\end{split}$

Taking the trace of both sides let’s us conclude that

(74)
$-\|[A,B]\|_{F}^{2}=2\operatorname{Re}{\left(\mathrm{Tr}{(ABA^{\dagger}B^{\dagger})}\right)}-\mathrm{Tr}{\left(A^{\dagger}ABB^{\dagger}+AA^{\dagger}B^{\dagger}B\right)}.$

Rearranging we conclude that

(75)
$\operatorname{Re}{\left(\mathrm{Tr}{(ABA^{\dagger}B^{\dagger})}\right)}=\frac{1}{2}\mathrm{Tr}{\Big{(}AA^{\dagger}B^{\dagger}B+A^{\dagger}ABB^{\dagger}\Big{)}}-\frac{1}{2}\|[A,B]\|_{F}^{2}$

which proves the lemma.