Many-Body Superconductivity in Topological Flat Bands

Jonah Herzog-Arbeitman^1
Aaron Chew^1
Kukka-Emilia Huhtinen^2
Päivi Törmä^2
B. Andrei Bernevig^1,3,4

$^1$Department of Physics, Princeton University, Princeton, NJ 08544
$^2$Department of Applied Physics, Aalto University School of Science, FI-00076 Aalto, Finland
$^3$Donostia International Physics Center, P. Manuel de Lardizabal 4, 20018 Donostia-San Sebastian, Spain
$^4$IKERBASQUE, Basque Foundation for Science, Bilbao, Spain

In a flat band superconductor, bosonic excitations can disperse while unpaired electrons are immobile. To study this strongly interacting system, we construct a family of multi-band Hubbard models with exact eta-pairing ground states in all space groups. We analytically compute their many-body excitations and find that the Cooper pair bound states and density excitations obey an effective single-particle Hamiltonian written in terms of the non-interacting wavefunctions. These bound states possess a unique zero-energy excitation whose quadratic dispersion is determined by the minimal quantum metric. The rest of the bound state spectrum is classified by topological quantum chemistry, which we use to identify Cooper pairs with Weyl nodes, higher angular momentum pairing, and  fragile topology. We also add electron kinetic energy as a perturbation to show that the strongest pairing occurs at half filling and not at the highest density of states. This is similar in spirit to the superconductivity observed in twisted bilayer graphene.

Attractive interactions among electrons can lead to a pairing instability at the Fermi surface and create a superconducting ground state at low temperatures. Both the weakly and strongly interacting limits of these systems can be described by attractive Hubbard models. Such models, however, are rarely exactly solvable. Recent work[Törmä et al., 2022, Kopnin et al., 2011, Julku et al., 2016, Törmä et al., 2018, Tovmasyan et al., 2018, Julku et al., 2020, Tian et al., 2021, Peri et al., 2021, Hu et al., 2019, Huhtinen and Törmä, 2021, Wang et al., 2021, Iskin, 2021, Žitko and Pavešić, 2022, Kitamura et al., 2022] has generated great interest in flat band superconductivity where electrons are completely immobile without interactions. A mean-field calculation in the multiband Hubbard model shows [Peotta and Törmä, 2015] that the superfluid weight can be non-zero due to quantum geometry[Resta, 2011], indicating mobile Cooper pairs despite the infinite effective mass of the unpaired electrons. Lower bounds on the superfluid weight have been obtained from topology [Peotta and Törmä, 2015, Xie et al., 2020, Tovmasyan et al., 2016] and symmetry[Herzog-Arbeitman et al., 2022]. While numerical beyond-mean-field work exists [Julku et al., 2016, Peri et al., 2021, Herzog-Arbeitman et al., 2022, Chan et al., 2022, Chan et al., 2022, Hofmann et al., 2022, Hofmann et al., 2020, Orso and Singh, 2022] and confirms the mean-field picture, the few exact approaches on the attractive Hubbard model governing flat band superconductivity have been limited to ground state properties [Tovmasyan et al., 2016] or one dimension [Tovmasyan et al., 2018]. In this Letter, we show that the many-body problem is analytically tractable in any dimension subject to a simple requirement called the uniform pairing condition [Peotta and Törmä, 2015, Tovmasyan et al., 2016]. Our exact results characterize not only the superconducting ground state, but also its excitations which control the low-temperature behavior.

We consider tight-binding models with electron operators $c^\dag_{\mbf{k},\al,\sigma}$ with spin $\sigma= \u,\d$, orbital $= 1,…, N_orb$, and momentum $k$ in the Brillouin zone (BZ). Throughout, we assume $S_z$ spin conservation, e.g. no spin-orbit coupling. We consider kinetic terms with flat bands at zero energy labeled by $m=1,…,N_f$ separated by a large gap to the other bands (see Fig.~\ref{figure1}b for example). We assume the gap is much larger than the interaction strength $|U|$. Under this approximation, the kinetic term is completely characterized by the single-particle eigenvectors of the flat bands which form the $N_orb ×N_f$ matrix $[U_σ(k)]_m$. By orthonormality, $U^†_σ(k) U_σ(k) = 1_N_f$, and $P^σ_(k) = [U_σ(k) U^†_σ(k)]_$ is a rank $N_f$ projector, or $k$-space correlation function [Bernevig and Hughes, 2013]. The gap to the other bands can be sent adiabatically to infinity so all operators can be projected into the flat bands spanned by
\bea
\label{eq:anticomP}
\bar{c}^\dag_{\mbf{k},\al,\sigma} &= \sum_\be c^\dag_{\mbf{k},\be,\sigma} P^\sigma_{\be \al}(\mbf{k}), \  \{\bar{c}_{\mbf{k}',\al,\sigma}, \bar{c}^\dag_{\mbf{k},\be,\sigma}\} = \delta_{\mbf{k}\mbf{k}'} P^\sigma_{\al \be}(\mbf{k})
\eea
which form an overcomplete basis due to the projection (\App{app:singlepart}). The canonical flat band operators $γ^†_k,m,σ = ∑_ c^†_k,,σ U_,m^σ(k)$ alternatively form a complete basis of the projected Hilbert space. For brevity, we define $P(k) = P^(̆k)$. We assume  $S_z$ conservation and spinful time-reversal symmetry $𝒯$ which yields $P^(̣k) = P(-k)^*$.

\textit{Hamiltonian.} The $N_f$ flat bands have zero kinetic energy and the remaining bands are projected out, so the kinetic term is zero and only interactions remain. We study an attractive Hubbard model proposed in \Ref{PhysRevB.94.245149},
\bea
\label{eq:Ham}
H &= \frac{|U|}{2} \sum_{\mbf{R}\al} (\bar{n}_{\mbf{R},\al,\u}-\bar{n}_{\mbf{R},\al,\d})^2, \ \ \bar{n}_{\mbf{R},\al,\sigma} = \bar{c}^\dag_{\mbf{R},\al,\sigma}\bar{c}_{\mbf{R},\al,\sigma} \\
\eea
where $c̅^†_R,,σ$ is the Fourier transform of $c̅^†_k,,σ$. If we impose the uniform pairing condition (UPC)~[Peotta} and {T{\"o}rm{\"a}}, 2015, Tovmasyan et~al., 2016]
\bea
\label{eq:UPC}
\frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \al}(\mbf{k}) = \frac{N_f}{N_{orb}} \equiv \eps \leq 1\\
\eea
where the $$-independent constant $$ is fixed by $P(k) = N_f$, and $𝒩$ is the number of unit cells, then we find
\bea
H = \frac{\eps |U|}{2} \bar{N} - \sum_{\mbf{R}\al} \bar{n}_{\mbf{R},\al,\u}\bar{n}_{\mbf{R},\al,\d}
\eea
which is the usual Hubbard model up to a chemical potential (note that $N̅ = ∑_Rσ n̅_R,,σ$ is the projected number operator). The UPC thus allows to relate the flat band Hubbard Hamiltonian to the positive-definite Hamiltonian \Eq{eq:Ham}, which is
crucial to analytically solving the ground states and excitations.

Intuitively, if the orbitals are related by symmetry, the UPC should hold [Tovmasyan et~al., 2016]. Indeed, we prove that \Eq{eq:UPC} is enforced by symmetry. Our proof (\App{app:UPC}) uses Schur's lemma and some basic facts from topological quantum chemistry [Bradlyn} et~al., 2017, Cano} and {Bradlyn}, 2020, Bacry et~al., 1988, Aroyo et~al., 2006, Aroyo et~al., 2006]. We show that if the orbitals of the model can all be related by a symmetry operator (technically, if they form an irrep of a Wyckoff position), then the UPC holds. Thus the class of flat bands with the UPC is very large, and encompasses existing examples [Peri} et~al., 2021, Herzog-Arbeitman} et~al., 2022, Bergman et~al., 2008, Rhim and Yang, 2019, Song} et~al., 2020, Huhtinen} and {T{\"o}rm{\"a}}, 2021, T{\"o}rm{\"a}} et~al., 2018, Tovmasyan} et~al., 2018].

There is a generalization[Tovmasyan et~al., 2016] of \Eq{eq:UPC}: if the flat band wavefunction $[U_σ(k)]_m$ vanishes on some orbitals $$, then $n̅_R,α̃,σ= 0$ and the UPC is only required on $≠$. Calugaru et al.[C{\v{a}}lug{\v{a}}ru} et~al., 2021] construct models with $N_f = N_L-N_L̃$ perfectly flat bands on bipartite lattices consisting of two sublattices $L, L̃$ with $[U_σ(k)]_ m = 0$ where $N_L   (N_L̃)$ is the number of orbitals in $L   (L̃)$ per unit cell. Choosing the $L$ lattice to have the UPC, their method constructs exactly flat bands with the UPC in all space groups. In this case, \Eq{eq:UPC} holds with $= N_f/N_L$. In this manner, we construct a 10 band model (see Fig.~\ref{figure1}) with the momentum space irreps of twisted bilayer graphene (TBG)[Bistritzer} and {MacDonald}, 2011, Song} et~al., 2019, Po} et~al., 2018] and the UPC (see Fig.~\ref{figure1}), but with Wilson loop winding [Ahn} et~al., 2019, Bouhon} et~al., 2018, Yu} et~al., 2022] 2 (instead of 1 as in TBG).

\begin{figure*}
\includegraphics[width=\textwidth]{Figure_TBGwilson.pdf}
\caption{Ten band model inspired by TBG. ($a$) A bipartite lattice of $s$ ($p$) orbitals at the 6j (2b) Wyckoff positions reproduces the momentum space irreps [Song} et~al., 2019, Ahn} et~al., 2019, Po} et~al., 2018] of the Bistritzer-MacDonald (BM) model[Bistritzer} and {MacDonald}, 2011] with uniform pairing. ($b$) Band structure of the 10-band model (red) and BM model (black) with Wilson loop inset. ($c$) Cooper pair spectrum with the minimal quantum metric approximation in blue. $p$-wave pairing Leggett modes appear above the zero-energy $s$-wave mode (\App{app:Aaron}). ($d$) Superfluid weight $D_s$ from mean-field theory in units where $e,\hbar,k_B = 1$. The BKT temperature $T_c = \frac{\pi}{8} D_s(T_c)$ is estimated at $T_c \approx 0.06 |U|$.}
\label{figure1}
\end{figure*}

\textit{Superconducting Ground States.} Originally, $H$ in \Eq{eq:Ham} was proposed to exhibit an exact BCS-type groundstate despite the ill-defined Fermi surface of the flat band[Tovmasyan et~al., 2016]. In this work, we instead focus on the more physical ground states at \emph{fixed} particle number. We define the local spin operator $S̅^z_R, = n̅_R,,-n̅_R,,$ appearing in \Eq{eq:Ham}. Because of the projection, $c̅^†_R,,σ$ is extended on the scale of the correlation length, which cannot be adiabatically sent to zero (the atomic limit) if the bands are topological or obstructed [Herzog-Arbeitman} et~al., 2022]. This is crucial for superconductivity in transport since trivial flat bands in the atomic limit would make $H = |U|/2 ∑_R, (S̅^z_R,)^2$ a sum of commuting local terms with trivial dynamics and macroscopic degeneracies. Away from the trivial limit, $H$ retains a single nontrivial symmetry operator $η$
\bea
\label{eq:etaoperator}
\null [\eta^\dag, \bar{S}^z_{\mbf{R},\al}] = 0, \ \eta^\dag=  \sum_{\mbf{k}\al} \bar{c}^\dag_{\mbf{k},\al, \u} \bar{c}^\dag_{-\mbf{k},\al, \d} = \sum_{\mbf{R}\al} \bar{c}^\dag_{\mbf{R},\al, \u} \bar{c}^\dag_{\mbf{R},\al, \d} \ .
\eea
This result relies on the UPC, global $S_z$ conservation, and $𝒯$ (\App{app:Hubermodel}). Here $η^†$ is  simply the creation operator of a zero-momentum $s$-wave Cooper pair, and it forms an $su(2)$ algebra, $[N̅, η^†] = 2 η^†$, which extends $U(1)$ charge conservation [Tovmasyan et~al., 2016].  We define the states $|n⟩ ∝η^†n |0⟩$ which have fixed particle number $N̅ |n⟩ = 2 n |n⟩$ for $n = 0, …, 𝒩 N_f$. Yang first introduced similar eta-pairing states[Yang, 1989] on the square lattice, which is also bipartite. However, his operators carried momentum $π$ and were not ground states. In contrast, $|n⟩$ must be ground states because $H$ is positive semi-definite, and $H |n⟩ ∝η^†n H |0⟩ = 0$.

Because $|n⟩$ differs from the BCS-BEC[Chen} et~al., 2005, Randeria} and {Taylor}, 2013] states $e^z η^†|0⟩$ (which are also exact ground states[Tovmasyan et~al., 2016]), it may not be apparent that they describe a bosonic condensate. To establish this, we show that $|n⟩$ has off-diagonal long-range order in the two-particle density matrix[Yang, 1962]; this is intuitive since the $η$ pairs are superpositions of pairs on \textit{all} sites. In the states $|n⟩$, a generalized Wick's theorem [Peotta}, 2022, Molinari}, 2017, Fetter} and {Walecka}, 1971] holds for the correlators of the Wannier state operators $w^†_Rmσ = 1/√(𝒩) ∑_k e^i R·k γ^†_k,m,σ$ obeying ${w_R,m,σ,w^†_R',m',σ' } = δ_R,R' δ_m,m' δ_σ,σ'$. We find (\App{app:generatingfunctions})
\bea
\label{eq:ODLRO}
\lim_{|\mbf{R}-\mbf{R}'| \to \infty} \braket{n| w^\dag_{\mbf{R}' m' \d}w^\dag_{\mbf{R}' m' \u} w_{\mbf{R}m \u} w_{\mbf{R} m \d} |n} = \nu(1-\nu)
\eea
in the thermodynamic limit where $ν= n/N_f𝒩$ is the filling. Clearly the strongest condensate appears at half-filling of the flat bands. Similar dependence on the filling appears in mean-field results on the multi-band Hubbard model [Peotta} and {T{\"o}rm{\"a}}, 2015] which also indicate half-filling is favorable for the superconductor. We remark that the Wannier operators have a gauge freedom corresponding to the redefinition of eigenvectors $U_σ(k) →U_σ(k) 𝒲_σ(k)$ where $𝒲_σ(k) ∈U(N_f)$ [Marzari} et~al., 2012, Brouder et~al., 2007, Bernevig} et~al., 2006, Monaco} et~al., 2018, Monaco} et~al., 2016, Winkler} et~al., 2016, Soluyanov and Vanderbilt, 2011, Song} et~al., 2020, Song} et~al., 2020]. Nevertheless, \Eq{eq:ODLRO} shows that the resulting correlator is gauge-invariant. This condensate ground state is a universal result in all UPC lattices.

\textit{Electron Excitations.} The presence of an enlarged symmetry algebra and exact ground states suggests that more features of $H$ may be accessible. We now show that the electron (and hole) excitations above the $|n⟩$ ground states are exactly solvable. Our procedure[Bernevig} et~al., 2021, Kang} et~al., 2021, Bernevig} et~al., 2021, Herzog-Arbeitman et~al., 2022, Kang and Vafek, 2019] is to add a single electron $γ^†_k,m,σ$ to the ground state and calculate
\bea
\label{eq:Rham}
\null H \gamma^\dag_{\mbf{k},m,\sigma}\ket{n} = [H, \gamma^\dag_{\mbf{k},m,\sigma}]\ket{n} &= \sum_{m'} \gamma^\dag_{\mbf{k},m',\sigma} \ket{n} [R^\sigma(\mbf{k})]_{m'm}
\eea
since $H|n⟩=0$. Diagonalizing $R^σ(k)$, one obtains the charge $+1$ (electron) excitation spectrum and the eigenstates. Using $H = |U|/2 ∑_q S̅^z_-q, S̅^z_q,$ where $S̅^z_q,$ is the Fourier transform of $S̅^z_R,$ and $S̅^z_q,|n⟩=0$, we find
\begin{align}
\label{eq:Rham2}
\null [H, \gamma^\dag_{\mbf{k},m, \sigma}] \ket{n} &= \frac{|U|}{2} \sum_{\mbf{q}\al}   [\bar{S}^z_{-\mbf{q},\al} , [ \bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, m, \sigma}]] \ket{n}  \\
&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!= \frac{|U|}{2\mathcal{N}}\sum_{\mbf{q} \al,m'l} M^{m'l}_{\sigma,\al}(\mbf{k}+\mbf{q},-\mbf{q}) M^{lm}_{\sigma,\al}(\mbf{k},\mbf{q})  \gamma^\dag_{\mbf{k},m',\sigma} \ket{n} \nonumber
\end{align}
defining $M^mn_,σ(k,q) = U^σ*_m(k+q) U^σ_n(k)$ as the orbital-resolved form factors (\App{app:charge1}). Matching \Eq{eq:Rham2} to \Eq{eq:Rham} and simplifying the form factors yields
\begin{align}
\label{eq:chargegap}
\null [R^\sigma(\mbf{k})]_{m'm} &= \frac{|U|}{2\mathcal{N}}\sum_{\mbf{q} \al} U^{\sigma *}_{\al m'}(\mbf{k}) P^\sigma_{\al\al}(\mbf{k}+\mbf{q}) U^\sigma_{\al m}(\mbf{k}) \\
&= \frac{\eps |U|}{2} \sum_\al U^{\sigma *}_{\al m'}(\mbf{k})U^\sigma_{\al m}(\mbf{k}) = \frac{\eps |U|}{2} \delta_{m'm} \nonumber
\end{align}
where the UPC, \Eq{eq:UPC}, was crucial to compute the $q$ sum. Hence, the electron excitations are independent of $ν$, $N_f$-degenerate per spin, and exactly flat meaning that the single-electron excitations are immobile. Indeed, a mean-field treatment assuming a $k$-independent pairing gap reproduces the $|U|/2$ $s$-wave gap [Huhtinen} et~al., 2022]. The temperature scale $|U|$ describes the onset of pairing, but is expected to be somewhat above the superconducting critical temperature with a pseudogap intervening[Wang} et~al., 2020, Zhang} et~al., 2021, Hashimoto} et~al., 2014, Tovmasyan} et~al., 2018].

The electronic excitations characterize the stability of the superconducting phase, as in BCS theory. When dispersion is added to the flat bands, as in realistic systems, we determine the filling where superconductivity is strongest through the Richardson criterion[Richardson} and {Sherman}, 1964, Richardson}, 1977, Bernevig} et~al., 2021] $E_Δ(N+2) = E(N+2)-2E(N+1) +E(N)$ where $E(N)$ is the ground state energy at particle number $N$. Here $E_Δ$ is an estimate of the Cooper pair binding energy: $E_Δ< 0$ indicates the favorability of pair formation (at $T=0$). We can compute $E_Δ(N)$ perturbatively when dispersion $Ẽ_n(k)$ is added to the flat bands. For convenience, we assume $∑_k,n Ẽ_n(k) = 0$ in which case (\App{app:Richcrit})
\bea
\label{eq:Edelta}
- E_\Delta(\nu) = \eps |U| - |1-2\nu| \tilde{E}_\nu
\eea
where $Ẽ_ν> 0$ is $max_k,n Ẽ_n(k)$ ($|min_k,n Ẽ_n(k)|)$ for $ν≥1/2$ ($ν< 1/2$). Hence the strongest Cooper pair binding occurs at $ν= 1/2$, where the condensate is also strongest (see \Eq{eq:ODLRO}) in agreement with multi-band mean-field results [Peotta} and {T{\"o}rm{\"a}}, 2015], but in contrast to conventional BCS theory where pairing is maximized at the highest density of states. This echoes the phenomenology of twisted bilayer graphene where critical temperature is highest \emph{not} at the van Hove singularities of the flat bands [Cao} et~al., 2018, Oh} et~al., 2021, Park} et~al., 2020, Stepanov} et~al., 2020, Xie} et~al., 2019]. \Eq{eq:Edelta} is not valid when spin-orbit coupling is added, but we provide a general formula in \App{app:Richcrit}.

\textit{Cooper Pair Excitations and Bound States.} We now show that the two-particle excitations, containing the Cooper pair and density excitation spectrum, are analytically solvable. We will also find a set of excitations (Leggett modes [Leggett}, 1966]) in other pairing channels. We find that these modes can carry higher angular momentum, e.g. $d$-wave symmetry, and fragile topology. Our calculation shows in microscopic detail the formation of bosons from the underlying fermion bands.

The two-body excitations are characterized by their global quantum numbers, $S_z$ spin and charge. We start in the charge $+2$ sector and compute the scattering matrix
\begin{align}
\label{eq:Rcharge2}
\null &[H, \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}\gamma^\dag_{-\mbf{k},n,\sigma'}] \ket{n} \\
&=  \sum_{\mbf{k}'m'n'} \gamma^\dag_{\mbf{p}+\mbf{k}',m',\sigma}\gamma^\dag_{-\mbf{k}',n',\sigma'}  \ket{n} [R^{\sigma \sigma'}(\mbf{p})]_{\mbf{k}'m'n',\mbf{k}mn} \nonumber
\end{align}
which is a function of the total momentum $p$. The commutator can be computed since $S̅^z_R,|n⟩=0$. The full scattering matrix can be written (\App{app:charge2})
\bea
R^{\sigma \sigma'}(\mbf{p}) &=  \eps |U|  \mathbb{1} + s^z_\sigma  s^z_{\sigma'}\frac{ |U| }{\mathcal{N}} \mathcal{U}_{\sigma \sigma'}(\mbf{p})\mathcal{U}_{\sigma \sigma'}^\dag(\mbf{p}) \\
\eea
where $s^z_/̆ = ±1$ and $[𝒰_σσ'(p)]_kmn, = U^σ*_m (p+k) U^σ' *_n (-k)$ is a $𝒩N_f^2 ×N_L$ matrix ($𝒩$ is the number of $k$ points). If $σ= σ'$, then $𝒰_σσ(p)$ is symmetric under $p+k,m ↔-k,n$ whereas the fermion operators in \Eq{eq:Rcharge2} are anti-symmetric. Hence the anti-symmetrization of the fermions annihilates the $𝒰_σσ(p)𝒰^†_σσ(p)$ term, so all spin-aligned states have energy $|U|$ for all $p$ and are completely unpaired.

We now consider $σ= -σ'$. Denote $R^̣̆(p) = |U| 1 - 𝒰(p)𝒰^†(p)/𝒩$, dropping spin labels for brevity. First, note that $rank  𝒰 𝒰^†≤rank 𝒰 ≤N_L$ and hence $𝒰𝒰^†$ has a very large kernel of dimension at least $𝒩 N_f^2 - N_L$, corresponding to the particle-particle continuum. Its nonzero eigenvalues are identical to the eigenvalues of the $N_L ×N_L$ matrix $h(p) = 𝒰^†(p)𝒰(p)/𝒩$. This pairing Hamiltonian takes the universal form (\App{app:charge2})
\bea
\label{eq:hpairing}
h_{\al\be}(\mbf{p}) &= \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{k}+\mbf{p})P_{\be \al}(\mbf{k}) \\
\eea
which can be interpreted as the effective single-particle Hamiltonian describing the hopping of a tightly bound Cooper pair (\App{app:spectrumbounds}). However, we emphasize that the \Eq{eq:hpairing} derives from an exact many-body calculation, and its eigenvalues are precisely the binding energies of the various pairing channels.

We denote the eigen-decompositon of \Eq{eq:hpairing} as $h(p) u_μ(p) = _μ(p) u_μ(p)$ for $μ= 0,…,N_L-1$. The many-body energies of $R^̣̆(p)$ are $|U|$ for the $𝒩N_f^2-N_L$ unpaired states and $E_μ(p) = |U|(- _μ(p))$ for the $N_L$ bound states ( \Fig{figure1}c). Returning to \Eq{eq:Rcharge2}, the pair excitations can be written as $η^†_p,μ |n⟩$ where (\App{app:charge2})
\bea
\label{eq:etaeigenstates}
\eta^\dag_{\mbf{p},\mu} &= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k}\al} \frac{u^\al_\mu(\mbf{p})}{\sqrt{\eps_\mu(\mbf{p})}} \bar{c}^\dag_{\mbf{p}+\mbf{k},\al,\u} \bar{c}^\dag_{-\mbf{k},\al,\d}
\eea
which is the many-body analogue of the BCS pairing, with $μ$ enumerating the pairing channels. We will show that $η^†_p,μ$ can carry nontrivial topology protected by $G$.

We remark that the charge-0 density excitations are also exactly solvable. We determine (\App{app:goldstones}) that the spin $±1$ density excitations are gapped with energies above the $|U|$ pairing threshold. The two spin-0 density excitations, of the form $γ^†_p+k,σγ_k,σ$, in fact have a spectrum \emph{identical} to the spin-0 Cooper pairs. We will show momentarily that there is guaranteed to be a gapless Goldstone mode in this spectrum, similar to the Anderson-Bogoliubov phonon [Pitaevskii et~al., 2003, Feynman, 1954, Zwerger, 2011].
However our calculation shows this Goldstone mode is quadratic (\App{app:quantummetric}), whereas
the Anderson-Bogoliubov phonon typically has linear dispersion. The quadratic dispersion can be attributed to the breaking of the $η, η^†$ generators by fixing the particle number[Watanabe}, 2019, Zhang} et~al., 2021].

\textit{Boson Bands and Topology.} We have reduced the two-body scattering matrix to an effective band theory problem given by the Hamiltonian in \Eq{eq:hpairing}. Notably, this result is generic: it is valid for all lattices with the UPC. We now prove a number of universal features of $h(p)$.

First, $_μ(p) ≥0 $. This can be proved from the Schur product theorem[Prasolov, 1994], which states that the Hadamard product $[A ∘B]_ = A_ B_$ of positive semi-definite matrices is positive semi-definite. Therefore $h(p)$ is the sum of positive semi-definite terms, and is also positive semi-definite. Secondly, $_μ(p) ≤$ follows since \Eq{eq:Ham} is the sum of squares and hence $E_μ(p) = |U|(- _μ(p))>0$. Note that the largest eigenvalue of $h(p)$, i.e. the strongest pairing, corresponds to the lowest energy many-body state. We number the bands such that $_0(p) ≥_μ(p)$.

Because $η^†$ is a zero-energy, $s$-wave, two-body excitation with $E_0(0) = 0$, it follows that $_0(0) = $ (\App{app:spectrumbounds}). To study the small $p$ behavior around $_0(0) = $, we apply perturbation theory. The result is $E_0(p) = |U|/N_L g_ij p_i p_j + O(p^4)$ where $g_ij$ is the integrated minimal quantum metric [Huhtinen} et~al., 2022] (\App{app:quantummetric}). This amounts to a many-body proof of the result first reported in mean-field[Peotta} and {T{\"o}rm{\"a}}, 2015] and by a two-body calculation[T{\"o}rm{\"a}} et~al., 2018]. Notably, the bipartite construction[C{\v{a}}lug{\v{a}}ru} et~al., 2021] produces many examples of fragile or obstructed single-particle bands, for which $g_ij$ is lower bounded[Herzog-Arbeitman} et~al., 2022, Yu} et~al., 2022]. \textit{Hence, the low-energy spectrum is universally governed by the quantum geometry of the flat bands.} The presence of tightly bound pairs with quadratic dispersion suggests these models are in the Bose-Einstein condensate regime of superconductivity [Randeria} and {Taylor}, 2013, Chen} et~al., 2005, Julku et~al., 2021, Zwerger, 2011]. Transport in the superconductor can be characterized by the $T=0$ mean-field superfluid weight $D_s$, shown in Fig.~\ref{figure1}(d). For degenerate flat
bands, a mean-field calculation shows (\App{app:meanfield})
\begin{equation}
[D_s]_{ij} = 8 \frac{\eps |U|}{V_c}\nu(1-\nu)g_{ij}, \quad g_{ij} = \frac{1}{\mathcal{N}} \sum_\mbf{k} g_{ij}(\mbf{k})
\label{eq.sfw}
\end{equation}
where $V_c$ is the volume of a unit cell and $g_ij(k)$ is the minimal quantum metric of the degenerate bands [Xie} et~al., 2020, Huhtinen} et~al., 2022, Mera} and {Mitscherling}, 2022, Herzog-Arbeitman} et~al., 2022]. Eq.~\eqref{eq.sfw} predicts $D_s/|U|≈0.277$ for our model, in good agreement with \Fig{figure1}d.

The eigenvalues of $h_(p)$ correspond to the pairing channels, with $s$-wave pairing at low energies, but higher energy modes with nontrivial pairing symmetries appear (see \Fig{figure1}c) and contribute to the finite temperature behavior. Although the zero-energy $s$-wave mode is universal across all space groups $G$, the higher energy spectrum varies. Nevertheless, we can use topological quantum chemistry to determine all possible Cooper pair band connectivities[Bradlyn et~al., 2018, Christensen} et~al., 2021, Elcoro} et~al., 2020, Watanabe} et~al., 2018, Kruthoff} et~al., 2017] and topology for each $G$ (\App{app:Cooperpairspectrum}).

The pairing Hamiltonian $h(p)$ in an effective tight-binding model which inherits \emph{bosonic} single-particle symmetries from the underlying electron space group $G$. \Eq{eq:hpairing} shows $h_(p) = h^*_(-p)$ so $h(p)$ has a spinless time-reversal symmetry. We can also show (\App{eq:CPsymmetry}) that $h(p)$ inherits the space group $G$ but in the real representation $|D_[g]|, g∈G$, where $D_[g]$ is the representation of the spin-$$̆ electrons [Alexandradinata et al., 2020]. The possible band structures with these orbitals are highly constrained by topological quantum chemistry since, by enforcing the UPC within our construction, the orbitals of $h(\mbf{p})$ form a single Wyckoff position[Cano and Bradlyn, 2020] and are symmetry-related. If the orbitals are at low-symmetry (non-maximal Wyckoff) positions, then the bands of $h(\mbf{p})$ are composite and can be deformed into disconnected atomic representations[Cano et al., 2018]. If the orbitals are located at high-symmetry (maximal Wyckoff) positions, the bands are either all connected through high-symmetry points of the BZ or decomposable into topological bands[Bradlyn et al., 2018, Cano et al., 2018]. However, because $h(\mbf{p})$ has a spinless time-reversal symmetry, no stable topology is possible [Brouder et al., 2007, Fang et al., 2012]. Instead, a non-zero symmetry indicator[Bradlyn et al., 2017, Po, 2020, Po et al., 2017, Kruthoff et al., 2017] identifies the presence of topologically protected Weyl nodes [Hughes et al., 2011, Elcoro et al., 2020]. In 2D, we check exhaustively that the only groups with fragile topological Cooper pairs are $G =p6, p6mm$. In $p6mm$, placing $s$ or $p$ orbitals at the kagome (3c) position induces a unique decomposition: the low-energy branch is a topologically trivial $s$-wave Cooper pair, and there is a gap to two connected bands with fragile topology where the Cooper pair wavefunction interpolates between $d$-wave and $p$-wave character across the BZ (app:examples). This decomposition is also possible in $p6$, but there is also a second decomposition with a nonzero symmetry indicator $\theta_2 \in\mathds{Z}_2$. In this case, the Cooper pair bands are connected by Weyl nodes. An exhaustive search of the 3D Wyckoff positions with the UPC yields 22 space groups with decomposable bands (app:bosontopology). We find instances where the low-energy branch is fragile topological. All 22 are kagome/pyrochlore groups where many promising flat band candidates [Regnault et al., 2022] have been identified. Our classification identifies flat band superconductors which host topological Cooper pairs and higher angular momentum spin-singlet pairing due to the multi-orbital lattice.

Conclusions. Although the BCS paradigm is successful in describing weakly-coupled superconductors, the numerous unresolved questions in high-$T_c$ materials[Joynt and Taillefer, 2002, Sigrist and Ueda, 1991, Lee et al., 2006] and moiré heterostructures[Balents et al., 2020, Törmä et al., 2022, Zhang et al., 2021, Park et al., 2021, Lake et al., 2022, Wang et al., 2021, Lian et al., 2019, Zhang et al., 2021]  demonstrate a pressing need for a strong-coupling theory of superconductivity[Zhang et al., 2021, Lapa and Levin, 2020, Lapa, 2020, Richardson, 1963, Ibañez et al., 2009, Claeys et al., 2016]. To this end, we introduced a broad family of flat band Hubbard models constructed within the framework of topological quantum chemistry. Their groundstates are bosonic condensates featuring maximal off-diagonal long-range order at half-filling and a linear-in-$|U|$ pairing gap. Computing the mass of the condensing Cooper pair, we established that symmetry-protected wavefunction topology coheres the superconductor with a nonzero superfluid weight. We also found exact Cooper pair bound states with $p,d,\dots$-wave symmetries below the pairing gap and completely classified Cooper pair bands with fragile topology. Finally, we solved the Goldstone and gapped spin-wave excitations which hint at BEC-like behavior. Further study of these exactly solvable models may provide a new route to high-temperature superconductivity.

B.A.B. and A.C. were supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 101020833),
the ONR Grant No. N00014-20-1-2303,
the Schmidt Fund for Innovative Research,
Simons Investigator Grant No. 404513,
the Packard Foundation,
the Gordon and Betty Moore Foundation through the EPiQS Initiative, Grant GBMF11070 and Grant No. GBMF8685 towards the Princeton theory program.
Further support was provided by the NSF-MRSEC Grant No. DMR-2011750,
BSF Israel US foundation Grant No. 2018226,
and the Princeton Global Network Funds.
JHA is supported by a Hertz Fellowship.
P.T. and K-E.H. acknowledge support by the Academy of Finland under project numbers 303351 and 327293. K-E.H. acknowledges support from the Magnus Ehrnrooth Foundation.

[Törmä et al., 2022]
Päivi Törmä, Sebastiano Peotta, and Bogdan A Bernevig.
Superconductivity, superfluidity and quantum geometry in twisted
multilayer systems.
Nature Reviews Physics, pages 1–15, 2022.

[Kopnin et al., 2011]
N. B. Kopnin, T. T. Heikkilä, and G. E. Volovik.
High-temperature surface superconductivity in topological flat-band
Phys. Rev. B, 83:0 220503, Jun 2011.
URL <https://link.aps.org/doi/10.1103/PhysRevB.83.220503>.

[Julku et al., 2016]
Aleksi Julku, Sebastiano Peotta, Tuomas I. Vanhala, Dong-Hee Kim, and
Päivi Törmä.
Geometric Origin of Superfluidity in the Lieb-Lattice Flat Band.
, 1170 (4):0 045303, July 2016.

[Törmä et al., 2018]
P. Törmä, L. Liang, and S. Peotta.
Quantum metric and effective mass of a two-body bound state in a
flat band.
, 980 (22):0 220511, December 2018.

[Tovmasyan et al., 2018]
Murad Tovmasyan, Sebastiano Peotta, Long Liang, Päivi
Törmä, and Sebastian D. Huber.
Preformed pairs in flat Bloch bands.
, 980 (13):0 134513, October 2018.

[Julku et al., 2020]
A. Julku, T. J. Peltonen, L. Liang, T. T. Heikkilä, and
P. Törmä.
Superfluid weight and Berezinskii-Kosterlitz-Thouless transition
temperature of twisted bilayer graphene.
, 1010 (6):0 060505, February 2020.

[Tian et al., 2021]
Haidong Tian, Shi Che, Tianyi Xu, Patrick Cheung, Kenji Watanabe,
Takashi Taniguchi, Mohit Randeria, Fan Zhang, Chun Ning Lau, and
Marc W. Bockrath.
Evidence for Flat Band Dirac Superconductor Originating from Quantum
arXiv e-prints, art. arXiv:2112.13401, December 2021.

[Peri et al., 2021]
Valerio Peri, Zhi-Da Song, B. Andrei Bernevig, and Sebastian D. Huber.
Fragile Topology and Flat-Band Superconductivity in the
Strong-Coupling Regime.
, 1260 (2):0 027002, January 2021.

[Hu et al., 2019]
Xiang Hu, Timo Hyart, Dmitry I. Pikulin, and Enrico Rossi.
Geometric and conventional contribution to the superfluid weight in
twisted bilayer graphene.
Phys. Rev. Lett., 123:0 237002, Dec 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.123.237002>.

[Huhtinen and Törmä, 2021]
Kukka-Emilia Huhtinen and Päivi Törmä.
Possible insulator-pseudogap crossover in the attractive Hubbard
model on the Lieb lattice.
, 1030 (22):0 L220502, June 2021.

[Wang et al., 2021]
Jie Wang, Jennifer Cano, Andrew J. Millis, Zhao Liu, and Bo Yang.
Exact Landau Level Description of Geometry and Interaction in a
, 1270 (24):0 246403, December

[Iskin, 2021]
M. Iskin.
Two-body problem in a multiband lattice and the role of quantum
Phys. Rev. A, 103:0 053311, May 2021.
URL <https://link.aps.org/doi/10.1103/PhysRevA.103.053311>.

[Žitko and Pavešić, 2022]
Rok Žitko and Luka Pavešić.
Yu-shiba-rusinov states, the bcs-bec crossover, and the exact
solution in the flat-band limit, 2022.
URL <https://arxiv.org/abs/2204.12610>.

[Kitamura et al., 2022]
Taisei Kitamura, Tatsuya Yamashita, Jun Ishizuka, Akito Daido, and
Youichi Yanase.
Superconductivity in monolayer FeSe enhanced by quantum geometry.
Physical Review Research, 40 (2):0 023232,
June 2022.

[Peotta and Törmä, 2015]
Sebastiano Peotta and Päivi Törmä.
Superfluidity in topologically nontrivial flat bands.
Nature Communications, 6:0 8944, November 2015.

[Resta, 2011]
R. Resta.
The insulating state of matter: a geometrical theory.
European Physical Journal B, 790 (2):0
121–137, January 2011.

[Xie et al., 2020]
Fang Xie, Zhida Song, Biao Lian, and B. Andrei Bernevig.
Topology-Bounded Superfluid Weight in Twisted Bilayer Graphene.
, 1240 (16):0 167002, April 2020.

[Tovmasyan et al., 2016]
Murad Tovmasyan, Sebastiano Peotta, Päivi Törmä, and Sebastian D. Huber.
Effective theory and emergent $\text{SU}(2)$ symmetry in the flat
bands of attractive hubbard models.
Phys. Rev. B, 94:0 245149, Dec 2016.
URL <https://link.aps.org/doi/10.1103/PhysRevB.94.245149>.

[Herzog-Arbeitman et al., 2022]
Jonah Herzog-Arbeitman, Valerio Peri, Frank Schindler, Sebastian D.
Huber, and B. Andrei Bernevig.
Superfluid Weight Bounds from Symmetry and Quantum Geometry in Flat
, 1280 (8):0 087002, February 2022.

[Chan et al., 2022]
Si Min Chan, B. Grémaud, and G. G. Batrouni.
Pairing and superconductivity in quasi-one-dimensional flat-band
systems: Creutz and sawtooth lattices.
Phys. Rev. B, 105:0 024502, Jan 2022.
URL <https://link.aps.org/doi/10.1103/PhysRevB.105.024502>.

[Chan et al., 2022]
Si Min Chan, B. Grémaud, and G. G. Batrouni.
Designer Flat Bands: Topology and Enhancement of Superconductivity.
arXiv e-prints, art. arXiv:2206.10651, June 2022.

[Hofmann et al., 2022]
Johannes S. Hofmann, Erez Berg, and Debanjan Chowdhury.
Superconductivity, charge density wave, and supersolidity in flat
bands with tunable quantum metric.
arXiv e-prints, art. arXiv:2204.02994, April 2022.

[Hofmann et al., 2020]
Johannes S. Hofmann, Erez Berg, and Debanjan Chowdhury.
Superconductivity, pseudogap, and phase separation in topological
flat bands.
, 1020 (20):0 201112, November 2020.

[Orso and Singh, 2022]
Giuliano Orso and Manpreet Singh.
Pairs, trimers, and bcs-bec crossover near a flat band: Sawtooth
Phys. Rev. B, 106:0 014504, Jul 2022.
URL <https://link.aps.org/doi/10.1103/PhysRevB.106.014504>.

[Bernevig and Hughes, 2013]
B. Andrei Bernevig and Taylor L. Hughes.
Topological Insulators and Topological Superconductors.
Princeton University Press, student edition edition, 2013.
ISBN 9780691151755.
URL <http://www.jstor.org/stable/j.ctt19cc2gc>.

[Bradlyn et al., 2017]
Barry Bradlyn, L. Elcoro, Jennifer Cano, M. G. Vergniory, Zhijun
Wang, C. Felser, M. I. Aroyo, and B. Andrei Bernevig.
Topological quantum chemistry.
, 5470 (7663):0 298–305, Jul 2017.

[Cano and Bradlyn, 2020]
Jennifer Cano and Barry Bradlyn.
Band Representations and Topological Quantum Chemistry.
arXiv e-prints, art. arXiv:2006.04890, June 2020.

[Bacry et al., 1988]
H Bacry, L Michel, and J Zak.
Symmetry and classification of energy bands in crystals.
In Group Theoretical Methods in Physics, pages 289–308.
Springer, 1988.

[Aroyo et al., 2006]
MI Aroyo, JM Perez-Mato, Cesar Capillas, Eli Kroumova, Svetoslav Ivantchev,
Gotzon Madariaga, Asen Kirov, and Hans Wondratschek.
Bilbao crystallographic server: I. databases and crystallographic
computing programs.
ZEITSCHRIFT FUR KRISTALLOGRAPHIE, 221:0 15–27, 01

[Aroyo et al., 2006]
Mois I. Aroyo, Asen Kirov, Cesar Capillas, J. M. Perez-Mato, and Hans
Bilbao Crystallographic Server. II. Representations of
crystallographic point groups and space groups.
Acta Crystallographica Section A, 620 (2):0
115–128, Mar 2006b.
URL <https://doi.org/10.1107/S0108767305040286>.

[Bergman et al., 2008]
Doron L. Bergman, Congjun Wu, and Leon Balents.
Band touching from real-space topology in frustrated hopping models.
Phys. Rev. B, 78:0 125104, Sep 2008.
URL <https://link.aps.org/doi/10.1103/PhysRevB.78.125104>.

[Rhim and Yang, 2019]
Jun-Won Rhim and Bohm-Jung Yang.
Classification of flat bands according to the band-crossing
singularity of bloch wave functions.
Phys. Rev. B, 99:0 045107, Jan 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevB.99.045107>.

[Song et al., 2020]
Zhi-Da Song, Luis Elcoro, and B. Andrei Bernevig.
Twisted bulk-boundary correspondence of fragile topology.
Science, 3670 (6479):0 794–797, February

[Cǎlugǎru et al., 2021]
Dumitru Cǎlugǎru, Aaron Chew, Luis Elcoro, Yuanfeng Xu,
Nicolas Regnault, Zhi-Da Song, and B. Andrei Bernevig.
General construction and topological classification of crystalline
flat bands.
Nature Physics, 180 (2):0 185–189, December

[Bistritzer and MacDonald, 2011]
Rafi Bistritzer and Allan H. MacDonald.
Moiré bands in twisted double-layer graphene.
Proceedings of the National Academy of Science, 1080
(30):0 12233–12237, Jul 2011.

[Song et al., 2019]
Zhi-Da Song, Zhijun Wang, Wujun Shi, Gang Li, Chen Fang, and
B. Andrei Bernevig.
All Magic Angles in Twisted Bilayer Graphene are Topological.
, 1230 (3):0 036401, Jul 2019.

[Po et al., 2018]
Hoi Chun Po, Liujun Zou, T. Senthil, and Ashvin Vishwanath.
Faithful Tight-binding Models and Fragile Topology of Magic-angle
Bilayer Graphene.
arXiv e-prints, art. arXiv:1808.02482, Aug 2018.

[Ahn et al., 2019]
J. Ahn, S. Park, and B.-J. Yang.
Failure of Nielsen-Ninomiya Theorem and Fragile Topology in
Two-Dimensional Systems with Space-Time Inversion Symmetry: Application to
Twisted Bilayer Graphene at Magic Angle.
Physical Review X, 90 (2):0 021013, April

[Bouhon et al., 2018]
Adrien Bouhon, Annica M. Black-Schaffer, and Robert-Jan Slager.
Wilson loop approach to metastable topology of split elementary band
representations and topological crystalline insulators with time reversal
arXiv e-prints, art. arXiv:1804.09719, Apr 2018.

[Yu et al., 2022]
Jiabin Yu, Ming Xie, Fengcheng Wu, and Sankar Das Sarma.
Euler Obstructed Cooper Pairing in Twisted Bilayer Graphene: Nematic
Nodal Superconductivity and Bounded Superfluid Weight.
arXiv e-prints, art. arXiv:2202.02353, February 2022.

[Yang, 1989]
Chen Ning Yang.
pairing and off-diagonal long-range order in a
hubbard model.
Phys. Rev. Lett., 63:0 2144–2147, Nov 1989.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.63.2144>.

[Chen et al., 2005]
Qijin Chen, Jelena Stajic, Shina Tan, and K. Levin.
BCS BEC crossover: From high temperature superconductors to
ultracold superfluids.
Physics Reports, 4120 (1):0 1–88, June 2005.

[Randeria and Taylor, 2013]
Mohit Randeria and Edward Taylor.
BCS-BEC Crossover and the Unitary Fermi Gas.
arXiv e-prints, art. arXiv:1306.5785, June 2013.

[Yang, 1962]
C. N. Yang.
Concept of off-diagonal long-range order and the quantum phases of
liquid he and of superconductors.
Rev. Mod. Phys., 34:0 694–704, Oct 1962.
URL <https://link.aps.org/doi/10.1103/RevModPhys.34.694>.

[Peotta, 2022]
Sebastiano Peotta.
Superconductivity, generalized random phase approximation and linear
scaling methods.
arXiv e-prints, art. arXiv:2204.12163, April 2022.

[Molinari, 2017]
Luca Guido Molinari.
Notes on Wick's theorem in many-body theory.
arXiv e-prints, art. arXiv:1710.09248, October 2017.

[Fetter and Walecka, 1971]
Alexander L. Fetter and John Dirk Walecka.
Quantum theory of many-particle systems.
Courier Corporation, 2003, 1971.

[Marzari et al., 2012]
Nicola Marzari, Arash A. Mostofi, Jonathan R. Yates, Ivo Souza, and
David Vanderbilt.
Maximally localized Wannier functions: Theory and applications.
Reviews of Modern Physics, 840 (4):0
1419–1475, October 2012.

[Brouder et al., 2007]
Christian Brouder, Gianluca Panati, Matteo Calandra, Christophe Mourougane, and
Nicola Marzari.
Exponential localization of wannier functions in insulators.
Phys. Rev. Lett., 98:0 046402, Jan 2007.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.98.046402>.

[Bernevig et al., 2006]
B. Andrei Bernevig, Taylor L. Hughes, and Shou-Cheng Zhang.
Quantum Spin Hall Effect and Topological Phase Transition in HgTe
Quantum Wells.
Science, 3140 (5806):0 1757, Dec 2006.

[Monaco et al., 2018]
Domenico Monaco, Gianluca Panati, Adriano Pisante, and Stefan Teufel.
Optimal Decay of Wannier functions in Chern and Quantum Hall
Communications in Mathematical Physics, January 2018.

[Monaco et al., 2016]
D. Monaco, G. Panati, A. Pisante, and S. Teufel.
The Localization Dichotomy for gapped periodic quantum systems.
arXiv e-prints, art. arXiv:1612.09557, December 2016.

[Winkler et al., 2016]
Georg W. Winkler, Alexey A. Soluyanov, and Matthias Troyer.
Smooth gauge and Wannier functions for topological band structures
in arbitrary dimensions.
, 930 (3):0 035453, January 2016.

[Soluyanov and Vanderbilt, 2011]
Alexey A. Soluyanov and David Vanderbilt.
Wannier representation of ${\mathbb{z}}_{2}$ topological insulators.
Phys. Rev. B, 83:0 035108, Jan 2011.
URL <https://link.aps.org/doi/10.1103/PhysRevB.83.035108>.

[Song et al., 2020]
Zhi-Da Song, Luis Elcoro, Yuan-Feng Xu, Nicolas Regnault, and B. Andrei
Fragile Phases as Affine Monoids: Classification and Material
Physical Review X, 100 (3):0 031001, July

[Song et al., 2020]
Zhi-Da Song, Luis Elcoro, and B. Andrei Bernevig.
Twisted bulk-boundary correspondence of fragile topology.
Science, 3670 (6479):0 794–797, February

[Bernevig et al., 2021]
B. Andrei Bernevig, Biao Lian, Aditya Cowsik, Fang Xie, Nicolas
Regnault, and Zhi-Da Song.
Twisted bilayer graphene. V. Exact analytic many-body excitations in
Coulomb Hamiltonians: Charge gap, Goldstone modes, and absence of Cooper
, 1030 (20):0 205415, May 2021.

[Kang et al., 2021]
Jian Kang, B. Andrei Bernevig, and Oskar Vafek.
Cascades between Light and Heavy Fermions in the Normal State of
Magic-Angle Twisted Bilayer Graphene.
, 1270 (26):0 266402, December 2021.

[Herzog-Arbeitman et al., 2022]
Jonah Herzog-Arbeitman, Aaron Chew, Dmitri K Efetov, and B Andrei Bernevig.
Reentrant correlated insulators in twisted bilayer graphene at 25 t
(2 $\pi$ flux).
Physical Review Letters, 1290 (7):0 076401,

[Kang and Vafek, 2019]
Jian Kang and Oskar Vafek.
Strong coupling phases of partially filled twisted bilayer graphene
narrow bands.
Phys. Rev. Lett., 122:0 246401, Jun 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.122.246401>.

[Huhtinen et al., 2022]
Kukka-Emilia Huhtinen, Jonah Herzog-Arbeitman, Aaron Chew, Bogdan A.
Bernevig, and Päivi Törmä.
Revisiting flat band superconductivity: Dependence on minimal
quantum metric and band touchings.
, 1060 (1):0 014518, July 2022.

[Wang et al., 2020]
Zhiqiang Wang, Gaurav Chaudhary, Qijin Chen, and K. Levin.
Quantum geometric contributions to the BKT transition: Beyond mean
field theory.
, 1020 (18):0 184504, November 2020.

[Zhang et al., 2021]
Xu Zhang, Kai Sun, Heqiu Li, Gaopei Pan, and Zi Yang Meng.
Superconductivity and bosonic fluid emerging from Moiré flat
arXiv e-prints, art. arXiv:2111.10018, November

[Hashimoto et al., 2014]
Makoto Hashimoto, Inna M. Vishik, Rui-Hua He, Thomas P. Devereaux, and
Zhi-Xun Shen.
Energy gaps in high-transition-temperature cuprate superconductors.
Nature Physics, 100 (7):0 483–495, July 2014.

[Richardson and Sherman, 1964]
R. W. Richardson and N. Sherman.
Exact eigenstates of the pairing-force Hamiltonian.
Nuclear Physics, 52:0 221–238, April 1964.

[Richardson, 1977]
R. W. Richardson.
Pairing in the limit of a large number of particles.
Journal of Mathematical Physics, 180 (9):0
1802–1811, September 1977.

[Cao et al., 2018]
Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi
Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero.
Unconventional superconductivity in magic-angle graphene
, 5560 (7699):0 43–50, April 2018.

[Oh et al., 2021]
Myungchul Oh, Kevin P. Nuckolls, Dillon Wong, Ryan L. Lee, Xiaomeng
Liu, Kenji Watanabe, Takashi Taniguchi, and Ali Yazdani.
Evidence for unconventional superconductivity in twisted bilayer
, 6000 (7888):0 240–245, October 2021.

[Park et al., 2020]
Jeong Min Park, Yuan Cao, Kenji Watanabe, Takashi Taniguchi, and Pablo
Flavour Hund's Coupling, Correlated Chern Gaps, and Diffusivity in
Moiré Flat Bands.
arXiv e-prints, art. arXiv:2008.12296, August 2020.

[Stepanov et al., 2020]
Petr Stepanov, Ipsita Das, Xiaobo Lu, Ali Fahimniya, Kenji Watanabe,
Takashi Taniguchi, Frank H. L. Koppens, Johannes Lischner, Leonid
Levitov, and Dmitri K. Efetov.
Untying the insulating and superconducting orders in magic-angle
, 5830 (7816):0 375–378, July 2020.

[Xie et al., 2019]
Yonglong Xie, Biao Lian, Berthold Jäck, Xiaomeng Liu, Cheng-Li
Chiu, Kenji Watanabe, Takashi Taniguchi, B. Andrei Bernevig, and Ali
Spectroscopic signatures of many-body correlations in magic-angle
twisted bilayer graphene.
, 5720 (7767):0 101–105, Jul 2019.

[Leggett, 1966]
A. J. Leggett.
Number-Phase Fluctuations in Two-Band Superconductors.
Progress of Theoretical Physics, 360 (5):0
901–930, November 1966.

[Pitaevskii et al., 2003]
Lev P Pitaevskii, Sandro Stringari, and Sandro Stringari.
Bose-einstein condensation.
Number 116. Oxford University Press, 2003.

[Feynman, 1954]
R. P. Feynman.
Atomic theory of the two-fluid model of liquid helium.
Phys. Rev., 94:0 262–277, Apr 1954.
URL <https://link.aps.org/doi/10.1103/PhysRev.94.262>.

[Zwerger, 2011]
Wilhelm Zwerger.
The BCS-BEC crossover and the unitary Fermi gas, volume 836.
Springer Science & Business Media, 2011.

[Watanabe, 2019]
Haruki Watanabe.
Counting Rules of Nambu-Goldstone Modes.
arXiv e-prints, art. arXiv:1904.00569, April 2019.

[Prasolov, 1994]
Viktor Vasil_evich Prasolov.
Problems and theorems in linear algebra, volume 134.
American Mathematical Soc., 1994.

[Julku et al., 2021]
Aleksi Julku, Georg M. Bruun, and Päivi Törmä.
Quantum geometry and flat band bose-einstein condensation.
Phys. Rev. Lett., 127:0 170404, Oct 2021.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.127.170404>.

[Mera and Mitscherling, 2022]
Bruno Mera and Johannes Mitscherling.
Nontrivial quantum geometry of degenerate flat bands.
arXiv e-prints, art. arXiv:2205.07900, May 2022.

[Bradlyn et al., 2018]
Barry Bradlyn, L. Elcoro, M. G. Vergniory, Jennifer Cano, Zhijun Wang,
C. Felser, M. I. Aroyo, and B. Andrei Bernevig.
Band connectivity for topological quantum chemistry: Band structures
as a graph theory problem.
Phys. Rev. B, 97:0 035138, Jan 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevB.97.035138>.

[Christensen et al., 2021]
Thomas Christensen, Hoi Chun Po, John D. Joannopoulos, and Marin
Location and topology of the fundamental gap in photonic crystals.
arXiv e-prints, art. arXiv:2106.10267, June 2021.

[Elcoro et al., 2020]
Luis Elcoro, Zhida Song, and B. Andrei Bernevig.
Application of induction procedure and Smith decomposition in
calculation and topological classification of electronic band structures in
the 230 space groups.
, 1020 (3):0 035110, July

[Watanabe et al., 2018]
Haruki Watanabe, Hoi Chun Po, and Ashvin Vishwanath.
Structure and topology of band structures in the 1651 magnetic space
Science Advances, 40 (8):0 eaat8685, August

[Kruthoff et al., 2017]
Jorrit Kruthoff, Jan de Boer, Jasper van Wezel, Charles L. Kane, and
Robert-Jan Slager.
Topological Classification of Crystalline Insulators through Band
Structure Combinatorics.
Physical Review X, 70 (4):0 041069, October

[Alexandradinata et al., 2020]
A. Alexandradinata, J. Höller, Chong Wang, Hengbin Cheng, and Ling Lu.
Crystallographic splitting theorem for band representations and
fragile topological photonic crystals.
Phys. Rev. B, 102:0 115117, Sep 2020.
URL <https://link.aps.org/doi/10.1103/PhysRevB.102.115117>.

[Cano et al., 2018]
Jennifer Cano, Barry Bradlyn, Zhijun Wang, L. Elcoro, M. G.
Vergniory, C. Felser, M. I. Aroyo, and B. Andrei Bernevig.
Building blocks of topological quantum chemistry: Elementary band
, 970 (3):0 035139, Jan

[Bradlyn et al., 2018]
Barry Bradlyn, Zhijun Wang, Jennifer Cano, and B. Andrei Bernevig.
Disconnected Elementary Band Representations, Fragile Topology, and
Wilson Loops as Topological Indices: An Example on the Triangular Lattice.
arXiv e-prints, art. arXiv:1807.09729, July 2018.

[Cano et al., 2018]
Jennifer Cano, Barry Bradlyn, Zhijun Wang, L. Elcoro, M. G.
Vergniory, C. Felser, M. I. Aroyo, and B. Andrei Bernevig.
Topology of Disconnected Elementary Band Representations.
, 1200 (26):0 266401, June

[Fang et al., 2012]
Chen Fang, Matthew J. Gilbert, and B. Andrei Bernevig.
Bulk topological invariants in noninteracting point group symmetric
, 860 (11):0 115112, September 2012.

[Po, 2020]
Hoi Chun Po.
Symmetry indicators of band topology.
Journal of Physics Condensed Matter, 320
(26):0 263001, June 2020.

[Po et al., 2017]
Hoi Chun Po, Ashvin Vishwanath, and Haruki Watanabe.
Complete theory of symmetry-based indicators of band topology.
Nature Communications, 8:0 50, June 2017.

[Hughes et al., 2011]
Taylor L. Hughes, Emil Prodan, and B. Andrei Bernevig.
Inversion-symmetric topological insulators.
, 830 (24):0 245132, Jun 2011.

[Elcoro et al., 2020]
Luis Elcoro, Benjamin J. Wieder, Zhida Song, Yuanfeng Xu, Barry
Bradlyn, and B. Andrei Bernevig.
Magnetic Topological Quantum Chemistry.
arXiv e-prints, art. arXiv:2010.00598, October

[Regnault et al., 2022]
Nicolas Regnault, Yuanfeng Xu, Ming-Rui Li, Da-Shuai Ma, Milena
Jovanovic, Ali Yazdani, Stuart S. P. Parkin, Claudia Felser,
Leslie M. Schoop, N. Phuan Ong, Robert J. Cava, Luis Elcoro, Zhi-Da
Song, and B. Andrei Bernevig.
Catalogue of flat-band stoichiometric materials.
, 6030 (7903):0 824–828, March 2022.

[Joynt and Taillefer, 2002]
Robert Joynt and Louis Taillefer.
The superconducting phases of ${\mathrm{upt}}_{3}$.
Rev. Mod. Phys., 74:0 235–294, Mar 2002.
URL <https://link.aps.org/doi/10.1103/RevModPhys.74.235>.

[Sigrist and Ueda, 1991]
Manfred Sigrist and Kazuo Ueda.
Phenomenological theory of unconventional superconductivity.
Rev. Mod. Phys., 63:0 239–311, Apr 1991.
URL <https://link.aps.org/doi/10.1103/RevModPhys.63.239>.

[Lee et al., 2006]
Patrick A. Lee, Naoto Nagaosa, and Xiao-Gang Wen.
Doping a mott insulator: Physics of high-temperature
Rev. Mod. Phys., 78:0 17–85, Jan 2006.
URL <https://link.aps.org/doi/10.1103/RevModPhys.78.17>.

[Balents et al., 2020]
Leon Balents, Cory R. Dean, Dmitri K. Efetov, and Andrea F. Young.
Superconductivity and strong correlations in moiré flat bands.
Nature Physics, 160 (7):0 725–733, May 2020.

[Zhang et al., 2021]
Yiran Zhang, Robert Polski, Cyprian Lewandowski, Alex Thomson, Yang
Peng, Youngjoon Choi, Hyunjin Kim, Kenji Watanabe, Takashi
Taniguchi, Jason Alicea, Felix von Oppen, Gil Refael, and Stevan
Ascendance of Superconductivity in Magic-Angle Graphene
arXiv e-prints, art. arXiv:2112.09270, December

[Park et al., 2021]
Jeong Min Park, Yuan Cao, Liqiao Xia, Shuwen Sun, Kenji Watanabe,
Takashi Taniguchi, and Pablo Jarillo-Herrero.
Magic-Angle Multilayer Graphene: A Robust Family of Moiré
arXiv e-prints, art. arXiv:2112.10760, December 2021.

[Lake et al., 2022]
Ethan Lake, Adarsh S. Patri, and T. Senthil.
Pairing symmetry of twisted bilayer graphene: a phenomenological
synthesis, 2022.
URL <https://arxiv.org/abs/2204.12579>.

[Wang et al., 2021]
Yuxuan Wang, Jian Kang, and Rafael M. Fernandes.
Topological and nematic superconductivity mediated by ferro-SU(4)
fluctuations in twisted bilayer graphene.
, 1030 (2):0 024506, January

[Lian et al., 2019]
Biao Lian, Zhijun Wang, and B. Andrei Bernevig.
Twisted Bilayer Graphene: A Phonon-Driven Superconductor.
, 1220 (25):0 257002, June 2019.

[Lapa and Levin, 2020]
Matthew F. Lapa and Michael Levin.
Rigorous Results on Topological Superconductivity with Particle
Number Conservation.
, 1240 (25):0 257002, June 2020.

[Lapa, 2020]
Matthew F. Lapa.
Topology of superconductors beyond mean-field theory.
Physical Review Research, 20 (3):0 033309,
August 2020.

[Richardson, 1963]
RW Richardson.
A restricted class of exact eigenstates of the pairing-force

[Ibañez et al., 2009]
Miguel Ibañez, Jon Links, Germán Sierra, and Shao-You Zhao.
Exactly solvable pairing model for superconductors with
p$_{x}$+ip$_{y}$ -wave symmetry.
, 790 (18):0 180501, May 2009.

[Claeys et al., 2016]
Pieter W. Claeys, Stijn De Baerdemacker, and Dimitri Van Neck.
Read-green resonances in a topological superconductor coupled to a
Phys. Rev. B, 93:0 220503, Jun 2016.
URL <https://link.aps.org/doi/10.1103/PhysRevB.93.220503>.

[Wieder et al., 2018]
Benjamin J. Wieder, Barry Bradlyn, Zhijun Wang, Jennifer Cano, Youngkuk
Kim, Hyeong-Seok D. Kim, Andrew M. Rappe, C. L. Kane, and B. Andrei
Wallpaper fermions and the nonsymmorphic Dirac insulator.
Science, 3610 (6399):0 246–251, July 2018.

[Lieb, 1989]
Elliott H. Lieb.
Two theorems on the hubbard model.
Phys. Rev. Lett., 62:0 1201–1204, Mar 1989.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.62.1201>.

[Benalcazar et al., 2017]
Wladimir A. Benalcazar, B. Andrei Bernevig, and Taylor L. Hughes.
Electric multipole moments, topological multipole moment pumping,
and chiral hinge states in crystalline insulators.
, 960 (24):0 245115, Dec 2017.

[Benalcazar et al., 2017]
Wladimir A. Benalcazar, B. Andrei Bernevig, and Taylor L. Hughes.
Quantized electric multipole insulators.
Science, 3570 (6346):0 61?66, Jul 2017.
ISSN 1095-9203.
URL <http://dx.doi.org/10.1126/science.aah6442>.

[Vafek et al., 2017]
Oskar Vafek, Nicolas Regnault, and B. Andrei Bernevig.
Entanglement of exact excited eigenstates of the Hubbard model in
arbitrary dimension.
SciPost Physics, 30 (6):0 043, December 2017.

[Gao et al., 2022]
Jiacheng Gao, Zhaopeng Guo, Hongming Weng, and Zhijun Wang.
Magnetic band representations, Fu-Kane-like symmetry indicators and
magnetic topological materials.
arXiv e-prints, art. arXiv:2204.10556, April 2022.

[Alexandradinata et al., 2014]
A. Alexandradinata, Xi Dai, and B. Andrei Bernevig.
Wilson-Loop Characterization of Inversion-Symmetric Topological
Phys. Rev., B890 (15):0 155114, 2014.

[Song et al., 2021]
Zhi-Da Song, Biao Lian, Nicolas Regnault, and B. Andrei Bernevig.
Twisted bilayer graphene. II. Stable symmetry anomaly.
, 1030 (20):0 205412, May 2021.

[Liang et al., 2017]
Long Liang, Tuomas I. Vanhala, Sebastiano Peotta, Topi Siro, Ari Harju, and
Päivi Törmä.
Band geometry, berry curvature, and superfluid weight.
Phys. Rev. B, 95:0 024515, Jan 2017.
URL <https://link.aps.org/doi/10.1103/PhysRevB.95.024515>.

§ SINGLE-PARTICLE SYMMETRY AND TOPOLOGY

In this section, we define the notation used in the rest of the work for the single-particle projectors, projected many-body operators, and space group symmetries (app:projbands). We prove
that the uniform pairing condition can be symmetry enforced by choosing orbitals to transform as irreps of a single Wyckoff positions (app:UPC), and extend this result to bipartite lattices where perfectly flat bands can be systematically constructed.

§.§ Projected Bands

Our problem projects the main electronic degrees of freedom into the quasi-flat bands of the non-interacting Hamiltonian. In a strong coupling approximation, the quasi-flat bands are set to zero energy and the remaining bands are sent to infinite energy and projected out. As such, the physics depends entirely on the flat band wavefunctions.

We now introduce our notation which differs slightly from PhysRevB.94.245149 where only nondegenerate flat bands (doubly degenerate when including spin) were considered. Let $\mbf{R}$ index the unit cells and $\al = 1, \dots, N_{orb}$ index the orbitals at positions $\mbf{r}_\al$. For now we neglect spin, but it will be introduced later with $\sigma = \u,\d$. The electron operators (obeying canonical anti-commutation relations) are
\bea
\label{eq:cRck}
c^\dag_{\mbf{k},\al} = \frac{1}{\sqrt{\mathcal{N}}} \sum_\mbf{R} e^{- i \mbf{k} \cdot (\mbf{R} + \mbf{r}_\al)} c^\dag_{\mbf{R},\al}, \qquad c^\dag_{\mbf{R},\al} = \frac{1}{\sqrt{\mathcal{N}}}\sum_{\mbf{k}} e^{i \mbf{k} \cdot (\mbf{R} + \mbf{r}_\al)} c^\dag_{\mbf{k},\al} \\
\eea
where $𝒩$ is the total number of unit cells, which is the number of terms in the $R,k$ sums. The general single-particle Hamiltonian reads
\bea
\label{eq:gammasing}
\tilde{H} &=  \sum_{\mbf{k},\al \be}c^\dag_{\mbf{k},\al} \tilde{h}_{\al \be}(\mbf{k}) c_{\mbf{k},\be} = \sum_{\mbf{k},n} \tilde{E}_n(\mbf{k}) \gamma^\dag_{\mbf{k},n}  \gamma_{\mbf{k},n}, \qquad \gamma_{\mbf{k},n} &= \sum_\al [\tilde{U}^\dag(\mbf{k})]_{n,\al} c_{\mbf{k},\al}
\eea
and $h̃(k) = Ũ Ẽ Ũ^†$ where $Ẽ$ is the diagonal matrix of energies, $Ũ$ is the $N_orb ×N_orb$ eigenvector matrix, and the diagonalized electron operators obey
${γ_k,m,γ^†_k',n } = δ_k,k' [Ũ^†Ũ]_mn =  δ_k,k' δ_mn$. For brevity, we often suppress the $k$ dependence.

We assume there is a gap around the $N_f$ quasi-flat bands, which we will refer to as the flat bands; in the strong coupling approximation, we do not require them to be exactly flat. We project all operators into the flat bands by setting $γ_k,n = 0$ if $n$ is not one of the flat bands. Now denote $U(k)$ as the $N_orb ×N_f$ occupied eigenstate matrix which obeys $U^†(k) U(k) = 1_N_f$ and $U(k)U^†(k) = P(k)$ where $P$ is the $N_orb ×N_orb$ projector matrix obeying $P^2 = P, P^†= P, P = N_f$. We recall that in the exactly flat case, the eigenstates are degenerate and hence are only defined up to the gauge freedom $U(k) →U(k) 𝒲(k)$ where $𝒲(k) ∈U(N_f)$ is a unitary matrix with arbitrary $k$ dependence [Herzog-Arbeitman} et~al., 2022]. Both $U^†(k) U(k)$ and $U(k)U^†(k) $ are invariant under this gauge freedom. The projector $P(k)$ encodes the topology and quantum geometry of the bulk bands. It will play an essential role in the many-body problem, highlighting the interplay of topology and interactions.

We now define the projected operators. Because we have set all $γ_k,n$ to zero other than the flat band operators, an (over)complete basis of the projected states is given by
\bea
\label{eq:bargamma}
\bar{c}_{\mbf{k}, \al} \equiv  \sum_{n=1}^{N_{orb}} [\tilde{U}(\mbf{k})]_{\al,n} \gamma_{\mbf{k},n} &= \sum_{n=1}^{N_f} [U(\mbf{k})]_{\al,n} \gamma_{\mbf{k},n} = \sum_{\be} [U(\mbf{k}) U^\dag(\mbf{k})]_{\al \be} c_{\mbf{k},\be} = \sum_{\be} [P(\mbf{k})]_{\al \be} c_{\mbf{k},\be}
\eea
using \Eq{eq:gammasing}. Because the projector $P(k)$ is Hermitian, we have
\bea
\label{eq:cdatPc}
\bar{c}^\dag_{\mbf{k},\al}  &= \sum_\be c^\dag_{\mbf{k},\be}  P^*_{\al \be}(\mbf{k}) = \sum_\be c^\dag_{\mbf{k},\be}  P_{\be \al}(\mbf{k}) \ .
\eea
The $c̅^†_k,$ basis is invariant under the eigenvector gauge freedom (unlike the $γ^†_k,n$ operators). However, the barred operators do not obey canonical anti-commutation relations because they are an overcomplete basis, whereas the $γ^†_k,n$ operators form a complete basis. We compute
\bea
\label{eq:anticom}
\{\bar{c}_{\mbf{k},\al}, \bar{c}^\dag_{\mbf{k}',\be} \} = \sum_{\al'\be'} P_{\al \al'}(\mbf{k}) \{c_{\mbf{k},\al'},c^\dag_{\mbf{k}',\be'} \} P_{\be' \be}(\mbf{k}') = \delta_{\mbf{k}\mbf{k}'} [P(\mbf{k})^2]_{\al \be} = \delta_{\mbf{k}\mbf{k}'} P_{\al \be}(\mbf{k}) \ . \\
\eea
We now need the projected position operators in real space. They are defined in analogy to \Eq{eq:cRck} as
\bea
\label{eq:barcexpansoin}
\bar{c}_{\mbf{R},\al}  &=  \frac{1}{\sqrt{\mathcal{N}}}\sum_{\mbf{k}} e^{-i \mbf{k} \cdot (\mbf{R} + \mbf{r}_\al)} \bar{c}_{\mbf{k},\al} \\
&= \frac{1}{\sqrt{\mathcal{N}}}\sum_{\mbf{k}\be} e^{-i \mbf{k} \cdot (\mbf{R} + \mbf{r}_\al)} P_{\al \be}(\mbf{k}) c_{\mbf{k},\be} \\
&= \frac{1}{\sqrt{\mathcal{N}}}\sum_{\mbf{R}'} \lp \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k}\be} e^{-i \mbf{k} \cdot (\mbf{R} + \mbf{r}_\al - \mbf{R}'-\mbf{r}_\be)} P_{\al \be}(\mbf{k}) \rp c_{\mbf{R}',\be} \\
&= \sum_{\mbf{R}'\be}p_{\al \be}(\mbf{R} - \mbf{R}') c_{\mbf{R}',\be} \\
\eea
where we defined the Fourier transform of the projector, i.e. the real-space correlation function of the non-interacting flat bands, as
\bea
\label{eq:pRandembeddingV}
p_{\al \be}(\mbf{R} -\mbf{R}') = \frac{1}{\mathcal{N}} \sum_{\mbf{k}} e^{-i \mbf{k} \cdot (\mbf{R} + \mbf{r}_\al - \mbf{R}'-\mbf{r}_\be)} P_{\al \be}(\mbf{k})= \frac{1}{\mathcal{N}} \sum_{\mbf{k}} e^{- i \mbf{k} \cdot (\mbf{R} - \mbf{R}')} [V(\mbf{k}) P(\mbf{k}) V^\dag(\mbf{k})]_{\al \be}, \quad  [V(\mbf{k})]_{\al \be} = e^{-i \mbf{k} \cdot \mbf{r}_\al} \delta_{\al \be} \ . \\
\eea
The embedding matrix $V(k)$ defined in \Eq{eq:pRandembeddingV} is a unitary diagonal matrix. Although it is sometimes useful[Herzog-Arbeitman} et~al., 2022] to take the thermodynamic limit $𝒩→∞$, we will keep $𝒩$ finite throughout this work.

We prove some properties about $p(R)$. It is obvious that $p(R) = p^†(-R)$ by Hermiticity of $P(k)$ and $p(R) = N_f δ_R,0$ follows from $P(k) = N_f$. We also have the Fourier inversion
\bea
\label{eq:projdef}
P(\mbf{k}) = \sum_{\mbf{R}} e^{i \mbf{k} \cdot \mbf{R}} V^\dag(\mbf{k})p(\mbf{R})V(\mbf{k}) \ .
\eea
From the projector property $P^2 = P$, we find the convolution identity[Herzog-Arbeitman} et~al., 2022]
\bea
\label{eq:convolution}
p(\mbf{R}) &= \sum_{\mbf{R}'}  p(\mbf{R}-\mbf{R}')p(\mbf{R}') % = \sum_{\mbf{R}'}  p(\mbf{R}+\mbf{R}')p^\dag(\mbf{R}') \ . \\
\eea
which we use to derive the real space anti-commutation relations (and can also be recovered from Eq.~\ref{eq:anticom}). We compute
\bea
\label{eq:realspaceanticom}
\{\bar{c}_{\mbf{R},\al}, \bar{c}^\dag_{\mbf{R}',\be} \} &= \sum_{\mbf{S}\mbf{S}',\al'\be'} p_{\al \al'}(\mbf{R}-\mbf{S})  [p^\dag(\mbf{R}'-\mbf{S}')]_{\be' \be} \{c_{\mbf{S},\al'}, c^\dag_{\mbf{S}',\be'} \} \\
&= \sum_\mbf{S} [p(\mbf{R}-\mbf{S})p^\dag(\mbf{R}'-\mbf{S}) ]_{\al \be} \\
&= \sum_\mbf{S} [p(\mbf{R}-\mbf{R}'+\mbf{S})p^\dag(\mbf{S}) ]_{\al \be} \\
&=p_{\al \be}(\mbf{R}-\mbf{R}') \ . \\
\eea
In the thermodynamic limit where $𝒩 →∞$, the flat bands are gapped so $P(k)$ is continuous on the BZ and the real space correlation function $p(R)$ decays at least exponentially fast.

We now look at the projected density operators $n̅_k, = c̅^†_k, c̅_k,$. We use the anti-commutation relations \Eq{eq:anticom}, finding
\bea
\bar{n}_{\mbf{k},\al}^2 &= \bar{c}^\dag_{\mbf{k},\al} \bar{c}_{\mbf{k},\al} \bar{c}^\dag_{\mbf{k},\al} \bar{c}_{\mbf{k},\al} = \bar{c}^\dag_{\mbf{k},\al} P_{\al \al}(\mbf{k}) \bar{c}_{\mbf{k},\al} - \bar{c}^\dag_{\mbf{k},\al} \bar{c}^\dag_{\mbf{k},\al} \bar{c}_{\mbf{k},\al}  \bar{c}_{\mbf{k},\al} =  P_{\al \al}(\mbf{k})  \bar{n}_{\mbf{k},\al} \\
\eea
where the second term vanishes because $c̅^†_k,$ anti-commutes with itself. In real space, the density $n̅_R, = c̅^†_R, c̅_R,$ obeys
\bea
\label{eq:nR2}
\bar{n}_{\mbf{R},\al}^2 &= \bar{c}^\dag_{\mbf{R},\al} \bar{c}_{\mbf{R},\al} \bar{c}^\dag_{\mbf{R},\al} \bar{c}_{\mbf{R},\al} = \bar{c}^\dag_{\mbf{R},\al}p_{\al \al}(\mbf{R}-\mbf{R}) \bar{c}_{\mbf{R},\al} -\bar{c}^\dag_{\mbf{R},\al} \bar{c}^\dag_{\mbf{R},\al}  \bar{c}_{\mbf{R},\al} \bar{c}_{\mbf{R},\al} = p_{\al \al}(0) \bar{n}_{\mbf{R},\al}  \\
\eea
using \Eq{eq:realspaceanticom}. The correlation function elements $p_(0)$ are onsite and measure local density of each orbital. They will play an important role in the many-body problem.

One may also derive these expressions from the electron operators in the kinetic energy band basis, $γ_,m$.  In this basis, projection simply involves restricting the bands to $m ∈M$, where $M$ labels the low-energy flat bands.  The projected density reads
\begin{align}
{\bar n}_{\kk, \alpha} &= \sum_{m,n} U(\kk)^*_{\alpha, m}
U (\kk)_{\alpha, n} \gamma^\dagger_{\kk,m} \gamma_{\kk,n} \\
{\bar n}_{\kk, \alpha}^2 &= \sum_{m,n} U(\kk)^*_{\alpha, m}
U (\kk)_{\alpha, n} \gamma^\dagger_{\kk,m} \gamma_{\kk,n} \sum_{m',n'} U(\kk)^*_{\alpha, m'}
U (\kk)_{\alpha, n'} \gamma^\dagger_{\kk,m'} \gamma_{\kk,n'} \nonumber \\
&= \sum_{m,n} U(\kk)^*_{\alpha, m}
U (\kk)_{\alpha, n} \sum_{m',n'} U(\kk)^*_{\alpha, m'}
U (\kk)_{\alpha, n'}  \gamma^\dagger_{\kk,m} [\{\gamma_{\kk,n}, \gamma^\dagger_{\kk,m'}\} - \gamma^\dagger_{\kk,m'} \gamma_{\kk,n}] \gamma_{\kk,n'} \nonumber \\
&= \sum_{m,n} U(\kk)^*_{\alpha, m}
U (\kk)_{\alpha, n} \sum_{m',n'} U(\kk)^*_{\alpha, m'}
U (\kk)_{\alpha, n'}  \gamma^\dagger_{\kk,m} [\delta_{n,m'} - 0] \gamma_{\kk,n'} \nonumber \\
&= P_{\alpha\alpha}(\kk) {\bar n}_{\kk,\alpha}.
\label{}
\end{align}
The zero in the second to last line arises from the antisymmetry upon interchange of the $m,m'$ indices.

\subsection{Uniform Pairing Condition from the Space group}
\label{app:UPC}

As Tovmasyan, Peotta, T\"orm\"a, and Huber identified in \Ref{PhysRevB.94.245149}, the solvability of the Hubbard model groundstates and excitations hinges on enforcing a simple criterion called the \textit{uniform pairing condition} (UPC):  $∑_R n̅^2_R, ∝∑_R n̅_R, + N̅$ where $N̅$ is the total (projected) number operator. The UPC was originally proposed to ensure equal mean-field pairing strength on all the orbital sites, but the UPC is simply a property of the single-particle wavefunctions as we will show. Using \Eq{eq:nR2}, we expand $n̅^2_R,$ to find
\bea
\label{eq:UPCMB}
\sum_{\mbf{R}\al} \bar{n}^2_{\mbf{R},\al} = \sum_{\mbf{R}\al} p_{\al\al}(0) \bar{n}_{\mbf{R},\al}
\eea
To enforce \Eq{eq:UPCMB}, we must have
\bea
\label{eq:UPCapps}
p_{\al \al}(0) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al\al}(\mbf{k})  = \eps \text{ if } \bar{n}_{\mbf{R},\al} \neq 0,~~\forall \alpha
\eea
where the constant $$ is independent of $$ and is fixed by $p(0) = ∑_p_(0) = N_f$. We will derive the UPC (\Eq{eq:UPCMB}) from the appearance of the eta-pairing symmetry in the following section \App{app:Hubermodel}.

The case of $n̅_R, = 0$ can occur for a given $$ at all $R$ if the wavefunction of the flat bands is exactly zero on the orbital $$ (see \Eq{eq:cdatPc}). We will see that this is the case in the bipartite crystalline lattices we consider shortly[C{\v{a}}lug{\v{a}}ru} et~al., 2021].

We now show that space group symmetries can enforce the uniform pairing condition \emph{generically}: arbitrary hoppings (preserving $G$) can be chosen as long as the underlying orbitals obey simple symmetry conditions which we make precise shortly. We consider unitary space group symmetries $g ∈G,  g = T_δ g_s$ consisting of a symmorphic part $g_s$ (a rotation or reflection) and non-symmorphic part $T_δ$ (a translation by a fraction of a lattice vector). The vector representation of $g$ on real-space vectors $r$ is defined to be $gr = g_sr + δ$ [Wieder} et~al., 2018]. The representations in momentum space take the form $D[g] = e^i g_s k ·δ D[g_s]$ where $D[g_s]$ is the $N_orb ×N_orb$ representation on the orbitals. On the single-particle Hamiltonian, $g$ enforces
\bea
\label{eq:simhtilde}
e^{i g_s \mbf{k} \cdot \pmb{\delta}} D[g_s] \tilde{h}(\mbf{k}) \lp e^{i g_s \mbf{k} \cdot \pmb{\delta}} D[g_s] \rp^\dag = D[g_s] \tilde{h}(\mbf{k}) D^\dag[g_s] = \tilde{h}(g_s \mbf{k}) , \\
\eea
canceling the non-symmorphic phases. From \Eq{eq:simhtilde}, $D[g_s] P(k) D^†[g_s] =  P(g_sk)$[Herzog-Arbeitman} et~al., 2022, Wieder} et~al., 2018]. It is understood that the representation matrix $D[g]$ only depends on the symmorphic part of $g$, and henceforth we drop the $s$ subscript.

We also need some simple notions from topological quantum chemistry. In a lattice with space group $G$, a Wyckoff position [Cano} and {Bradlyn}, 2020, Bacry et~al., 1988] $x$ of multiplicity $n$ consists of the points $x_1, …, x_n$ in the unit cell which form an orbit of $G$: there exists $g∈G$ such that $x_i = g x_j$ for each $x_i,x_j ∈x$. The site symmetry group of a point $x$ is defined by $G_x = {g∈G| gx = x} ∈G$, and all $G_x$ for $x∈x$ are isomorphic. Hence in any lattice with $G$, we can organize the orbitals into irreps of $G_x$ at some Wyckoff position $x$. All such irreps are tabulated on the \href{https://www.cryst.ehu.es/cgi-bin/cryst/programs/bandrep.pl}{Bilbao Crystallographic Server} [Aroyo et~al., 2006, Aroyo et~al., 2006]. For instance in $G = p2$ generated by $C_2$ rotation and translations, the $1a$ position (of multiplicity 1) is the origin of the unit cell and $G_1a = {1, C_2}$ with two 1D irreps (even and odd under $C_2$). The $2e$ position (of multiplicity 2) is the generic nonmaximal position consisting of a pair of $C_2$-related points, so $G_2e$ is the trivial group ${1}$ (see \Fig{fig:wyckoff}).

\begin{figure*}
\includegraphics[height=.3\textwidth]{1a2c}
\caption{Unit cell of the wallpaper group $p2$ generated by $C_2$ and translations (we take $\mbf{a}_1 = (1,0),\mbf{a}_2 = (0,1)$ for ease). There are four maximal Wyckoff positons 1a = $(0,0)$,1b= $(0,1/2)$,1c= $(0,1/2)$,1d= $(1/2,1/2)$ which all have multiplicity 1 and site symmetry group $G_x = 2$. At 1a for instance, the site symmetry group is $\{1, C_2\}$. The non-maximal Wyckoff position 2e$=\{(x,y),(1-x,1-y)\}$ has multiplicity 2 and consists of two general inversion-related points in the unit cell. Since $C_2$ exchanges the two points in 2e, it is not part of the site symmetry group and thus $G_{2e} = 1$ containing only the identity.}
\label{fig:wyckoff}
\end{figure*}

We can now state our claim.

\textbf{Lemma}. Let $x$ be a Wyckoff position (not necessarily maximal) with site-symmetry group $G_x$ and irreps $χ$. The uniform pairing condition holds between all orbitals forming an irrep $χ$ induced from the Wyckoff position $x$.

\textbf{Proof.}  We define $G̃$ as the group of symmetries of the unit cell, i.e. the space group mod lattice translations. We define the site-symmetry group $H = G_x_1 ≅G_x$ where for concreteness $G_x_1$ is the point group of the site $x_1 ∈x$. For a Wyckoff position of multiplicity $n$ and an irrep $χ$ of $G_x$ with dimension $|χ|$ and (irreducible) representations $D^χ[h]$ for $h∈H$, the Wyckoff position contains $n |χ|$ orbitals per unit cell. We index the orbitals in a tensor product basis $= i,a$ where $i = 1, …, n$ and $a =1,…, |χ|$. It is useful to note that in this basis, the embedding matrices are block diagonal, $V_ia,jb[G] = e^- i r_i ·G δ_ij δ_ab$. It is natural to employ a coset decomposition separating the onsite symmetries of $H$ and the symmetries that connect different points in the Wyckoff position. The formal expression is $G̃ = g̃_1 H + …+ g̃_n H$,  $n = |G̃|/|H|$ such that irreps of $H$ form an induced representation in $G̃$ [Song} et~al., 2020]. Hence the representation matrices of $G$ can be chosen as
\bea
\label{eq:symmetryforms}
D[h] &= \mathbb{1}_n \otimes D^\chi[h], \quad h \in H \\
D[\tilde{g}_i] &= P_{\tilde{g}_i} \otimes \mathbb{1}_{|\chi|}, \qquad \tilde{g}_i \in \{\tilde{g}_1, \dots, \tilde{g}_n\}
\eea
where $[P_g̃]_ij = δ_g̃x_i,x_j$ is a permutation matrix[Song} et~al., 2020] on the $n$ sites of the Wyckoff position, and the Kronecker delta is understood mod lattice vectors. First we prove uniform pairing holds within the $|χ|$ orbitals at a given $x_i$.
Choose an arbitrary diagonal block of $P(k)$ denoted $[P_i(k)]_ab = P_ia,ib(k)$. We have denoted $α= ia$ and $β= ib$.  Note that $P_i(k) = P_i(k+G)$ because the embedding matrix acts as $V_ia,ib[G] = e^-i x_i ·G δ_ab$, proportional to the identity. Then because the BZ is symmetric under $H$, we have
\bea
\frac{1}{\mathcal{N}} \sum_\mbf{k} P_i(\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_i(h\mbf{k}) =  D^\chi[h] \lp \frac{1}{\mathcal{N}} \sum_\mbf{k} P_i(\mbf{k}) \rp D^\chi[h]^\dag
\eea
but because $D^χ[h]$ is irreducible, by Schur's lemma we find $1/𝒩 ∑_k P_i(k) = _i 1_|χ|$ where $_i$ is an arbitrary constant.  That is, the only matrix invariant under an irreducible representation is proportional to the identity.  If there were another invariant matrix, call it $M$, that matrix $M$ would commute with all symmetries and thus will block-diagonalize the representation matrices, implying the representation was reducible.

To prove that $_i = _j$, relating different blocks, we use the fact that there exists $g̃ ∉H$ such that $g̃ x_i = x_j$. Then because $D[g̃] P(k) D^†[g̃] = P(g̃k)$, we find that
\bea
\frac{1}{\mathcal{N}} \sum_\mbf{k} P_i(\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_i(g\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_j(\mbf{k})
\eea
such that $_i = _j$. Thus $1/𝒩 ∑_k P_(k) = $ for all $$ in a the irrep $χ$ of the Wyckoff position $x$. Because $P(k) = N_f$, the number of projected bands, we have
\bea
\label{eq:UPCint}
\frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al\al}(\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} \frac{1}{n |\chi|}\Tr P(\mbf{k}) = \frac{N_f}{n |\chi|} \equiv \eps \leq 1  .
\eea
The value of $$ has physical meaning. From \Eq{eq:projdef}, we see that $= p_(0)$, which is the diagonal of the onsite correlation function. Although we have discussed only unitary symmetry groups in this proof, anti-unitary spatial symmetries, which do not act on spin, also enforce uniform pairing. This is easily seen because $P(k)$ is Hermitian and its diagonal elements must be real, so complex conjugation acts trivially. In \App{app:Aaron}, we use the spinless symmetry $C_6z T$ to enforce uniform pairing between 6 orbitals. \\

%Using \Eq{eq:projdef}, we can rewrite this equation in real space as
%\sum_{\mbf{R}} e^{i \mbf{k} \cdot \mbf{R}} D[g] V^\dag(\mbf{k})p(\mbf{R}) V(\mbf{k}) D^\dag[g] &= \sum_{\mbf{R}} e^{i g\mbf{k} \cdot \mbf{R}} V^\dag(g\mbf{k})p(\mbf{R})V(g\mbf{k}) \\
%Because $V(\mbf{k})$ is a diagonal matrix of phases, $[V^\dag(g\mbf{k})p(\mbf{R})V(g\mbf{k})]_{\al \al} = p_{\al\al}(\mbf{R})$. To treat the lefthand side of \Eq{eq:PRg}, we show that $V(\mbf{k}) D^\dag[g]  = D^\dag[g] W_g(\mbf{k})$ where $W_g(\mbf{k})$ is also diagonal. The proof is direct:
%\null [V(\mbf{k}) D^\dag[g]]_{\al \be} &= e^{- i \mbf{k} \cdot \mbf{r}_\al} D_{\al \be}[g^{-1}] = D_{\al \be}[g^{-1}] e^{- i \mbf{k} \cdot (g^{-1}\mbf{r}_\be - \mbf{R}_\be[g])} = [D[g] W_g(\mbf{k})]_{\al \be}
%where we used the fact that $D_{\al \be}[g^{-1}]$ nonzero only if $g^{-1} \mbf{r}_\be - \mbf{r}_\al$ is a lattice vector. In \Eq{eq:Wprooflaticevector}, we have denoted this lattice vector by $\mbf{R}_\be[g]$. It is possible to omit the $\al$ dependence because, given $\mbf{r}_\be$ and $D[g]$, there is a unique orbital site $\mbf{r}_\al = g^{-1} \mbf{r}_\be $ mod lattice vectors within the unit cell. Now note that $W_g(\mbf{k})$ is also a diagonal matrix of phases, so
%\sum_{\mbf{R}} e^{i \mbf{k} \cdot \mbf{R}} [D[g] V^\dag(\mbf{k})p(\mbf{R}) V(\mbf{k}) D^\dag[g] ]_{\al \al} &= \sum_{\mbf{R}} e^{i g\mbf{k} \cdot \mbf{R}} [ V^\dag(g\mbf{k})p(\mbf{R})V(g\mbf{k}) ]_{\al \al} \\
%\sum_{\mbf{R}} e^{i \mbf{k} \cdot \mbf{R}} [D[g] p(\mbf{R}) D^\dag[g] ]_{\al \al} &= \sum_{\mbf{R}} e^{i g\mbf{k} \cdot \mbf{R}} p_{\al \al}(\mbf{R}) \\
%leading to $[D[g] p(\mbf{R}) D^\dag[g] ]_{\al \al} = p_{\al \al}(g\mbf{R})$ by matching terms. Setting $\mbf{R}=0$, we obtain
%\null [D[g] p(0) D^\dag[g] ]_{\al \al} = p_{\al \al}(0) \ .
%We now use \Eq{eq:P0symmetry} to show that the uniform pairing condition can always be enforced by requiring the orbitals to form the 1D irrep of a single Wyckoff position (not necessarily maximal). To prove this as follows.
%%To do so
%%We also have the identity
%%\null [D[g] V(\mbf{k}) D^\dag[g] ]_{\al \be} &= \sum_{\al'} D_{\al \al'}[g] e^{- i \mbf{k} \cdot \mbf{r}_{\al'}} D_{\be \al'}[g]^*
%%\sum_{\mbf{R}} e^{i \mbf{k} \cdot \mbf{R}} V(g\mbf{k}) D[g] V^\dag(\mbf{k})p(\mbf{R})V(\mbf{k}) D^\dag[g] V^\dag(g\mbf{k}) &= \sum_{\mbf{R}} e^{i \mbf{k} \cdot \mbf{R}} p(g\mbf{R}) \ . \\
%%Let us examine the unitary $V(g\mbf{k}) D[g] V^\dag(\mbf{k})$ in detail. We have
%%\null [V(g\mbf{k}) D[g] V^\dag(\mbf{k})]_{\al \be} &= \sum_{\al'\be'} e^{- i g\mbf{k} \cdot \mbf{r}_\al} \delta_{\al \al'} D_{\al' \be'}[g] e^{i \mbf{k} \cdot \mbf{r}_\be} \delta_{\be' \be} \\
%%&=  D_{\al \be}[g] e^{- i g\mbf{k} \cdot \mbf{r}_\al+ i \mbf{k} \cdot \mbf{r}_\be} \\
%%&=  D_{\al \be}[g] e^{- i \mbf{k} \cdot (g^{-1} \mbf{r}_\al -\mbf{r}_\be)} \ . \\
%%Note that whenever $D_{\al \be}[g] \neq 0$, $g^{-1} \mbf{r}_\al -\mbf{r}_\be = -g^{-1} (g \mbf{r}_\be - \mbf{r}_\al)$ must be a lattice vector. In general this lattice vector will be nonzero and $g$ will relate different entries of $p(\mbf{R})$ to different lattice sites, schematically $p(\mbf{R}'_{\al \be})$. However, we will see that the uniform pairing condition is enforced in the simplest case, when $g^{-1} \mbf{r}_\al -\mbf{r}_\be = 0$ for $D_{\al \be}[g]  \neq 0$.
%%As such we have
%%\null \big[ V(g\mbf{k}) D[g] V^\dag(\mbf{k})p(\mbf{R})V(\mbf{k}) D^\dag[g] V^\dag(g\mbf{k})  \big]_{\al \be} = D_{\al \al'}[g] e^{- i \mbf{k} \cdot (g^{-1} \mbf{r}_\al -\mbf{r}_{\al'})} p_{\al' \be'}(\mbf{R}) D_{\be' \be}[g^\dag] e^{i \mbf{k} \cdot (g^{-1} \mbf{r}_\be -\mbf{r}_\be')}
%%so by matching terms we derive $V(g\mbf{k}) D[g] V^\dag(\mbf{k})p(\mbf{R})V(\mbf{k}) D^\dag[g] V^\dag(g\mbf{k}) =  p(g\mbf{R})$ and thus
%%\null [V(g\mbf{k}) D[g] V^\dag(\mbf{k}),p(0)] = 0 \ . \\
%%We first consider the case where we enforce $p_{\al\al}(0) = \nu$ for all $\al$.
%Consider a Wyckoff position $w$, not necessarily maximal, with multiplicity $m$ and orbitals at $\mbf{r}_{1}, \cdots, \mbf{r}_{m}$ forming the orbit of the point $\mbf{r}_{1}$ within the unit cell. The site-symmetry group $G_w = \{g\in G| g\mbf{r}_{\al} = \mbf{r}_{\al}\}$ consists of symmetries which leave the electrons on the Wyckoff position invariant. By definition there exists $h \notin G_w$ such that $\mbf{r}_\al = h \mbf{r}_1$ for each $\al$. Onsite, the electrons transform in representations of $G_w$, which define the atomic limits in $G$. They are tabulated on the Bilbao crystallographic server. We now require that the orbitals form a 1D irrep of a single Wyckoff position so that there is a single electron at each $\mbf{r}_\al$. In this case, we define the notation $g \be$ to refer to the unique index $\al$ such that $g \mbf{r}_\be = \mbf{r}_\al$ mod lattice vectors. If $g \in G_w$, then $g \al = \al$. In general, we have $D_{\al \be}[g] = \la_{\al}\delta_{\al,g\be}$ where $|\la_\al| = 1$ by unitarity of $D[g]$.
%We now compute
%\null [D[g] p(0) D^\dag[g]]_{\al \al} &= \sum_{\al' \be'} D_{\al \al'}[g] p_{\al' \be'}(0) D^*_{\al \be'}[g] \\
%&= \sum_{\al' \be'}  \la_{\al}\delta_{\al,g\al'}p_{\al' \be'}(0)  \la^*_{\al}\delta_{\al,g\be'} \\
%&= p_{g^{-1} \al g^{-1} \al}(0) \\
%and thus by \Eq{eq:P0symmetry}, $p_{\al \al}(0) = p_{g^{-1} \al g^{-1} \al}(0)$. Because there exists $g \in G$ such that $g^{-1} \al = \be$ for all $\be$, we see that $p_{\al \al}(0)$ is independent of $\al$. \jha{add figures with examples in $p4$, the $2c$ position and the $4d$ position.}

We now extend this result to the case of a bipartite crystalline lattice consisting of the sublattices $L$ and $L̃$. The integers $N_L$ and $N_L̃$ denote the number of orbitals on each sublattice, so $N_orb =  N_L + N_L̃$. Without loss of generality we take $N_L > N_L̃$. We consider single-particle Hamiltonians of the form
\bea
\label{eq:starterham}
\tilde{h}(\mbf{k}) &= \bpm 0 & S_\mbf{k} \\ S_\mbf{k}^\dag & B_\mbf{k}  \epm
\eea
where $S_k$ is a $N_L ×N_L̃$ rectangular matrix and $B_k$ is an $N_L̃ ×N_L̃$ Hermitian matrix. The number of flat bands is $N_f = N_L - N_L̃ > 0$ and their wavefunctions obey [C{\v{a}}lug{\v{a}}ru} et~al., 2021]
\bea
\bpm 0 & S_\mbf{k} \\ S_\mbf{k}^\dag & B_\mbf{k}  \epm \bpm \psi_\mbf{k} \\ 0\epm &= \bpm 0 \\ S^\dag_\mbf{k} \psi_\mbf{k} \epm = 0
\eea
where the existence of $N_f = N_L - N_L̃$ wavefunctions $ψ_k,n, n = 1,…, N_f$ in the nullspace of $S^†_k$ is guaranteed by rank deficiency. Thus the flat band projector
\bea
P(\mbf{k}) &= \sum_{n=1}^{N_f} \bpm \psi_{\mbf{k},n} \\ 0\epm \bpm \psi^\dag_{\mbf{k},n} & 0\epm = \sum_{n=1}^{N_f} \bpm \psi_{\mbf{k},n} \psi^{\dag}_{\mbf{k},n} & 0 \\ 0 & 0 \epm
\eea
has nonzero entries only in the $L$ sublattice. As such, the projected operators $c̅_k, = P_(k)c_k,$ are only nonzero when $$ is an orbital in the $L$ sublattice, and so $n̅_R,∈L̃ = 0$. Hence for \Eq{eq:UPCMB} to hold, we only require $p_(0)$ to be independent of $$ within the $L$ sublattice. As proven in the Lemma above, this is accomplished by requiring the $L$ sublattice to consist of orbitals within a single irrep. More than one irrep would no longer guarantee the uniform pairing condition.  Again because $p(0) = N_f = N_L - N_L̃$, we have
\bea
p_{\al \al}(0) = \frac{N_f}{N_{L}} \equiv \eps < 1 \ .
\eea
For simplicity of notation, we will not distinguish between sums over $∈L ⊕L̃$ and $∈L$. This is because the flat band wavefunctions have no weight in the $L̃$ sublattice, and so the sums are interchangable (any additional terms belonging to $L̃$ are $0$).

The notable feature of our construction is the total freedom to choose the hopping elements once the space group and orbitals are fixed. Using the bipartite crystalline lattice construction, we obtain perfectly flat bands obeying the uniform pairing condition in all space groups. The topology of the flat bands can be determined directly as the formal subtraction of the band representations of the two sublattices, and in many cases is fragile topological or obstructed atomic[C{\v{a}}lug{\v{a}}ru} et~al., 2021]. If the formal subtraction results in a ``negative" irrep at high symmetry momentum $K$, then the flat bands are gapless with the dispersive bands at $K$. In this case, projection into the flat bands is not justified, but this case can be analyzed very generally in mean field theory [Huhtinen} et~al., 2022].

\Ref{2021NatPh..18..185C} also constructs auxiliary models from the bipartite Hamiltonian in \Eq{eq:starterham}. By integrating out the $L̃$ sublattice, one obtains the positive semi-definite $N_L×N_L$ Hamiltonian $h̃_L(k) = S_kS^†_k$ which retains the $N_L - N_L̃$ zero-energy flat bands with wavefunctions $ψ_k,n$, as well as $N_L̃$ dispersive bands. These models $h̃_L(k)$ also possess the uniform pairing condition.

\subsubsection{Spin}
\label{app:spin}

To add spin, we impose $S_z$ conservation and spinful time-reversal symmetry $𝒯$. Together, $S_z$ and $𝒯$ ensure the existence of an eta-pairing operator as we show in \App{app:Hubermodel}. Enforcing these symmetries, the single-particle Hamiltonian takes the form
\bea
\tilde{h}_s(\mbf{k}) &= \bpm \tilde{h}_\u(\mbf{k}) & \\ &  \tilde{h}_\d(\mbf{k}) \\ \epm = \bpm \tilde{h}(\mbf{k}) & \\ &  \tilde{h}^*(-\mbf{k}) \\ \epm
\eea
where we used $D[𝒯] = i σ_y K$, so in particular $𝒯$ acts diagonally on the orbital index $,$.  Here $σ_y$ acts on the spin indices and $K$ is complex conjugation. Because of the global symmetry $S_z$, we can fix the gauge of $𝒯$ on the eigenstates so that $U^(̣k) = U^(̆-k)^*$  and $P^(̣k) = P^(̆-k)^*$. As such we have
\bea
\label{eq:TRSpR}
p_{\al \be}^\d(\mbf{R} -\mbf{R}')
%&= \frac{1}{V} \sum_\mbf{k} e^{-i \mbf{k} \cdot (\mbf{R}+\mbf{r}_\al-\mbf{R}'-\mbf{r}_\be)} P_{\al \be}^\d(\mbf{k})
= \frac{1}{\mathcal{N}} \sum_\mbf{k} e^{-i \mbf{k} \cdot (\mbf{R}+\mbf{r}_\al-\mbf{R}'-\mbf{r}_\be)} P_{\al \be}^\u(-\mbf{k})^* = \lp \frac{1}{\mathcal{N}} \sum_\mbf{k} e^{-i \mbf{k} \cdot (\mbf{R}+\mbf{r}_\al-\mbf{R}'-\mbf{r}_\be)} P_{\al \be}^\u(\mbf{k}) \rp ^*= p^\u_{\al \be}(\mbf{R} -\mbf{R}')^*
\eea
which is just the statement that $𝒯$ acts locally in real space. If the uniform pairing condition is obeyed by $p^(̣0)$, it is also by $p^(̆0)$ using $𝒯$ and $S_z$.

%If instead we impose $SU(2)$ spin symmetry, we have $\tilde{h}_\u(\mbf{k}) = \tilde{h}_\d(\mbf{k}) \equiv \tilde{h}(\mbf{k})$. In this case, $P^\u(\mbf{k}) = P^\d(\mbf{k}) = P(\mbf{k})$ and the uniform pairing condition is obviously obeyed by both spins. As we will show in the next section, $SU(2)$ symmetry alone is not sufficient to protect an eta-pairing operator. To do so, we must additionally require particle-hole symmetry $\mathcal{P}$. On the electron operators, $\mathcal{P}$ is unitary and obeys $\mathcal{P} c_{\mbf{R},\al} \mathcal{P}^\dag = c^\dag_{\mbf{R},\al}$. Using $[\tilde{H}, \mathcal{P}] = 0$, we derive $D[\mathcal{P}] \tilde{h}(\mbf{k}) D[\mathcal{P}]^{-1} = - \tilde{h}(-\mbf{k})$ where $D[\mathcal{P}] = K$. Thus $\mathcal{P}$ imposes a relation on the occupied eigenvector matrix. Note that for a general Hamiltonian with particle-hole symmetry, we can define $U(\mbf{k})$ to be the $N_{orb} \times N_{occ}$ matrix of a particle-hole symmetric set of occupied eigenvectors, i.e. the bands within the energy window $\pm \Delta$ assuming there is a gap in the spectrum at $E= \Delta$. Then $U(\mbf{k})$ obeys
%\tilde{h}(\mbf{k}) U(\mbf{k}) = U(\mbf{k}) \tilde{E}(\mbf{k})
%where $\tilde{E}(\mbf{k})$ is a diagonal matrix of energies sorted in decreasing order. Note that particle-hole symmetry in the spectrum enforces $-\tilde{E}(\mbf{k}) = J \tilde{E}(-\mbf{k})J$ where $J^2=1$ is the exchange matrix with 1 on the anti-diagonal and zero otherwise. Because $\tilde{h}(\mbf{k})^* = - \tilde{h}(-\mbf{k})$, we derive
%U^*(\mbf{k}) \tilde{E}(\mbf{k}) = \tilde{h}(\mbf{k})^* U^*(\mbf{k}) &= - \tilde{h}(-\mbf{k}) U^*(\mbf{k})
%so $\tilde{h}(\mbf{k}) U^*(-\mbf{k})= U^*(-\mbf{k}) (-\tilde{E}(\mbf{k})) = U^*(-\mbf{k}) J \tilde{E}(-\mbf{k})J$. Comparing with \Eq{eq:PHhuUE}, we see that we can pick a particle-hole symmetry eigenvector gauge where
%U^*(-\mbf{k}) J = U(\mbf{k})  \\
%which immediately implies $P(\mbf{k}) = P^*(-\mbf{k})$. Following the same steps as \Eq{eq:TRSpR}, we see that this implies $p_{\al \be}(\mbf{R})^* = p_{\al \be}(\mbf{R})$.

\section{Hubbard Model and Eta-pairing Symmetry}
\label{app:Hubermodel}

In this section, we construct a class of positive semi-definite attractive Hubbard models using the uniform pairing condition generalizing the model of \Ref{PhysRevB.94.245149} (\App{eq:hubbardmodels}). The existence of an eta-pairing symmetry in the models we consider is protected by $S_z$ and $𝒯$. We derive the form of the eta-pairing symmetry (\App{app:etaops}) and show that it expands the $U(1) ×U(1)$ charge-spin symmetry group to a non-abelian $SU(2)$ symmetry group (\App{app:symalg}).

\subsection{Attractive Hubbard Models}
\label{eq:hubbardmodels}

We now add a Hubbard term $H$ to the non-interacting Hamiltonian $H̃$ in \Eq{eq:gammasing}. The form of the Hubbard term was discovered in \Ref{PhysRevB.94.245149} and yields an eta-pairing symmetry, $[H, η] = 0$. The first step is to define a projected Hermitian spin operator
\bea
\bar{S}^z_{\mbf{R},\al} &= \bar{c}^\dag_{\mbf{R},\al,\u}\bar{c}_{\mbf{R},\al,\u}-\bar{c}^\dag_{\mbf{R},\al,\d}\bar{c}_{\mbf{R},\al,\d} = \bar{n}_{\mbf{R},\al,\u}-  \bar{n}_{\mbf{R},\al,\d} \ .\\
\eea
With the uniform pairing condition \Eq{eq:UPCapps}, the form of this operator naturally gives an attractive Hubbard model. We define
\bea
\label{eq:HuberN}
H = \frac{1}{2}|U| \sum_{\mbf{R},\al} (\bar{S}^z_{\mbf{R},\al})^2  \ .
\eea
The crucial feature of $H$ is that it is the sum of squares and thus is positive semi-definite. Hence a state $|ψ⟩$ obeying $H|ψ⟩ = 0$ must be a ground state. Now we show that $H$ is an attractive Hubbard model. Using \Eq{eq:nR2}, we compute
\bea
\label{eq:Hhubbardform}
H &= \frac{1}{2} |U| \sum_{\mbf{R},\al} (\bar{n}_{\mbf{R},\al,\u}-  \bar{n}_{\mbf{R},\al,\d} )^2  \\
&= \frac{1}{2}|U| \sum_{\mbf{R},\al} -2 \bar{n}_{\mbf{R},\al,\u} \bar{n}_{\mbf{R},\al,\d} + p^\u_{\al \al}(0) \bar{n}_{\mbf{R},\al,\u}+ p^\d_{\al \al}(0) \bar{n}_{\mbf{R},\al,\d}  \\
&= -|U| \sum_{\mbf{R},\al} \bar{n}_{\mbf{R},\al,\u} \bar{n}_{\mbf{R},\al,\d} + \frac{\eps |U|}{2} \sum_{\mbf{R} ,\al} \bar{n}_{\mbf{R},\al,\u}+ \bar{n}_{\mbf{R},\al,\d}  \\
&= -|U| \sum_{\mbf{R},\al} \bar{n}_{\mbf{R},\al,\u} \bar{n}_{\mbf{R},\al,\d} + \frac{\eps |U|}{2} \bar{N}\\
\eea
which is a Hubbard interaction with strength $|U|$ and a chemical potential term $|U|/2 N̅$. To get from line 2 to line 3, we have employed the uniform pairing condition $p^_̆(0) = p^_̣(0) = |U|/2$.  We have also used the fact that
\bea
\sum_{\mbf{R},\al,\sigma} \bar{n}_{\mbf{R},\al,\sigma} &= \frac{1}{\mathcal{N}} \sum_{\mbf{k},\mbf{k}'} \sum_{\mbf{R},\al,\sigma}  e^{- i (\mbf{k}-\mbf{k}') \cdot (\mbf{R}+ \mbf{r}_\al)} \bar{c}^\dag_{\mbf{k}' \al,\sigma} \bar{c}_{\mbf{k} \al,\sigma}  = \sum_{\mbf{k}\al,\sigma}  \bar{c}^\dag_{\mbf{k} \al,\sigma} \bar{c}_{\mbf{k} \al,\sigma}  = \sum_{\mbf{k},n,\sigma} \gamma^\dag_{\mbf{k},n,\sigma} \gamma_{\mbf{k},n,\sigma}  = \bar{N} \\
\eea
where the final sum is over the $N_f$ flat bands. Let us make a brief comment about the form of the attractive Hubbard term. The form of $∑_ n̅_R,, n̅_R,,$ is somewhat special because each orbital $$ only couples to itself (but opposite in spin). If there is only one orbital per Wyckoff position (i.e. a 1D irrep of $G_x$), this is the same type of term that one would get from a purely local Coulomb interaction. However if the irrep $χ$ is multi-dimensional, like a pair of $p_x, p_y$ orbitals per site, then a Coulomb interaction coupling the \emph{total} densities would be $∑_i (∑_a n̅_R,ia,)(∑_a n̅_R,ia,)$ in the tensor product notation of \App{app:UPC} where $a$ indexes the orbitals at a given point $x_i ∈x$.  Given that the mechanism of generating attractive Hubbard interactions is typically phonons with a complicated coupling structure, we view \Eq{eq:Hhubbardform} as a simple model capturing the essential physics. Although we focus here on an attractive interaction, a particle-hole transformation[Lieb, 1989] on the $$̣ spins maps all our results to the repulsive case with ferromagnetism.

§.§ Eta Operators

Generically, the $\bar{S}^z_{\mbf{R},\al}$ operators do not commute among themselves because of the projection (see eq:realspaceanticom). However, in the trivial atomic limit case where the orbitals are totally decoupled, $\bar{S}^z_{\mbf{R},\al} \to S_{\mbf{R},\al}$ is commuting for all $\mbf{R},\al$, yielding an infinite number of conserved charges. This is the classical limit where the model contains only strictly local density operators. This limit does not have superconductivity because there is no coherence – every site is decoupled. This is an indication that some obstruction to atomic localization, like topology or symmetry protected obstructed states, is necessary for the projected strong-coupling model to have a superconducting phase.

There is still a nontrivial symmetry operator when the projected bands are not atomic, the eta-pairing symmetry [Tovmasyan et al., 2016]. We use the most general two-body, charge $+2$, spin 0, momentum 0 operator for our ansatz of the $\eta^\dag$ operator. It takes the form

η^†   = ∑_q γ^†_q,m, Ω_mn(q) γ^†_-q,n,

involving the yet undetermined matrix $\Omega_{mn}(\mbf{q})$. We now must discuss the invariance of $\eta^\dag$ under the eigenvector gauge freedom $U^\sigma(\mbf{k}) \to U^\sigma(\mbf{k}) \mathcal{W}^\sigma(\mbf{k})$. Clearly for $\eta^\dag$ to be a physical symmetry, it must be invariant under the eigenvector gauge freedom which severely restricts the form of $\eta$ to a simple ansatz. We will find that with $S_z$ and $\mathcal{T}$, this ansatz is successful. Since $S_z$ and $\mathcal{T}$ are the symmetries relevant for real materials, we only consider this case explicitly in our work. However, we mention that the combination of particle-hole symmetry $\mathcal{P}$ and $SU(2)$ spin symmetry also protects the existence of an $\eta$ operator.

The operators $\gamma^\dag_{\mbf{q},n,\sigma}$ are not gauge invariant. Under the eigenvector gauge freedom, we have $\gamma^\dag_{\mbf{q},n,\sigma} \to \sum_{n'}\gamma^\dag_{\mbf{q},n',\sigma} \mathcal{W}^\sigma_{n'n}(\mbf{q})$ where $\mathcal{W}^\sigma(\mbf{q})$ is a unitary matrix. Hence
η^†   →∑_q,mnm'n' γ^†_q,m', 𝒲_m' m^(̆q) Ω_mn(q)  γ^†_-q,n', 𝒲_n' n^(̣-q)

= ∑_q,m'n' γ^†_q,m', [𝒲^(̆q) Ω(q) 𝒲^†̣(-q)^*]_m' n' γ^†_-q,n',  .

In order for $\eta^\dag$ to be invariant, we must have

𝒲^(̆q) Ω(q) 𝒲^†̣(-q)^*    =   Ω(q)  .
The natural object $\Omega$ transforming in this way is a sewing matrix[Benalcazar et al., 2017, Benalcazar et al., 2017] relating states at $\mbf{q}$ and $-\mbf{q}$. Because of the conjugation in eq:etainvariance, the sewing matrix must be that of an anti-unitary symmetry, like time-reversal.

Fixing the gauge of the $S_z$ and $\mathcal{T}$ sewing matrices so that $U^\u(-\mbf{q})^* = U^\d(\mbf{q})$, we can choose eigenstates such that $\mathcal{W}^\u(\mbf{q}) = \mathcal{W}^\d(-\mbf{q})^*$ and thus $\mathcal{W}^\u(\mbf{q}) \Omega(\mbf{q}) \mathcal{W}^{\d \dag}(-\mbf{q})^* = \mathcal{W}^\u(\mbf{q}) \Omega(\mbf{q}) \mathcal{W}^{\u}(\mbf{q})^\dag$. Thus $\eta^\dag$ is invariant under the gauge freedom if $\mathcal{W}^\u(\mbf{q}) \Omega(\mbf{q}) \mathcal{W}^{\u}(\mbf{q})^\dag = \Omega(\mbf{q}) $ which is satisfied for all arbitrary $\mathcal{W}^\u(\mbf{q})$ iff $\Omega_{mn}(\mbf{q}) \propto \delta_{mn}$.  Thus we arrive at the ansatz

η^†   = ∑_q,m γ^†_q,m, γ^†_-q,m, in the gauge-fixed wavefunctions where $U^\u(\mbf{k}) = U^\d(-\mbf{k})^*$. The $\eta^\dag$ operator takes a simpler form in terms of the projected operators (which are gauge invariant). Recalling $\gamma^\dag_{\mbf{q},m,\sigma} = \sum_\al \bar{c}^\dag_{\mbf{q},\al,\sigma} U^\sigma_{\al, m}(\mbf{q})$, we have

η^†   = ∑_q,m c̅^†_q,, U^_̆m(q) U^_̣m(-q) c̅^†_-q,, = ∑_q, c̅^†_q,, P^_̆(q) c̅^†_-q,, =  ∑_q c̅^†_q,, c̅^†_-q,, =  ∑_R c̅^†_R,, c̅^†_R,, .
In the last equality of eq:etadagniceform, we Fourier transformed to real space. Physically, $\eta^\dag$ can be understood as creating a Cooper pair at zero momentum. Note that the wavefunction of $\bar{c}^\dag_{\mbf{R},\al,\u} \bar{c}^\dag_{\mbf{R},\al,\d}$ decays exponentially outside the correlation length of $p(\mbf{R})$, which describes the “size" of the Cooper pair bound state.

We now proceed to the commutator of $\eta^\dag$ and $\bar{S}_{\mbf{R},\al}$:
[S̅^z_R,,η^†]    = ∑_R' [c̅^†_R,,c̅_R,,-c̅^†_R,,c̅_R,,,c̅^†_R',, c̅^†_R',,]

=   ∑_R' c̅^†_R,, [c̅_R,,,c̅^†_R',, c̅^†_R',,] -c̅^†_R,,[c̅_R,,,c̅^†_R',, c̅^†_R',,]

= ∑_R' c̅^†_R,, p_^(̆R-R') c̅^†_R',, + c̅^†_R,, c̅^†_R',, p_^(̣R-R')

where we used the anti-commutator $\{\bar{c}_{\mbf{k},\al,\sigma}, \bar{c}^\dag_{\mbf{k},\be,\sigma'} \} = \delta_{\mbf{k}\mbf{k}'} P_{\al \be}(\mbf{k}) \delta_{\sigma \sigma'}$ (see eq:realspaceanticom). We expand the $\bar{c}^\dag_{\mbf{R}',\be,\sigma}$ operators using eq:barcexpansoin to find
[S̅^z_R,,η^†]    =  ∑_R'R”,' c̅^†_R,, c^†_R”,', p_^(̆R-R')   p^*̣_'(R'-R”)+ c̅^†_R,, c^†_R”,', p_'^*̆(R'-R”) p_^(̣R-R')

=  ∑_R'R”,' c̅^†_R,, c^†_R”,', p_^(̆R-R')   p^_'(R'-R”)+ c̅^†_R,, c^†_R”,', p_'^(R'-R”) p_^(̣R-R')

= ∑_R'R”,' c̅^†_R,, c^†_R”,', [p^(̆R-R') p^(R'-R”)]_'+ c̅^†_R,, c^†_R”,', [p^(̣R-R') p^(R'-R”)]_'

where we used $p^\u_{\al\be}(\mbf{R})^* = p^\d_{\al\be}(\mbf{R})$ thanks to $\mathcal{T}$ in app:spin. Finally, the convolution identity of $p(\mbf{R})$ in eq:convolution gives
[S̅^z_R,,η^†]    =  ∑_R”,' c̅^†_R,, c^†_R”,', p_'^(R-R”)+ c̅^†_R,, c^†_R”,', p_'^(R-R”)

=  ∑_R”,' c̅^†_R,, c^†_R”,', p_'^*̣(R-R”)+ c̅^†_R,, c^†_R”,', p_'^*̆(R-R”)

= ( c̅^†_R,, c̅^†_R,, + c̅^†_R,, c̅^†_R,, )

= 0
again using eq:barcexpansoin to recover the projected operators in the penultimate line. Thus we have proven $\eta^\dag$ commutes with $\bar{S}^z_{\mbf{R},\al}$, so it is trivially established that $[H, \eta^\dag] = 0$. It follows that the eta-pairing states are all eigenstates:
H η^†n |0⟩    = η^†n H |0⟩= 0
and in particular must be groundstates because $H$ is positive semi-definite. Note that because $[\eta^\dag, \bar{S}^z_{\mbf{R},\al}] = 0$, one could consider much more general (though not necessarily physical) Hamiltonians with an eta-pairing symmetry than the onsite $ H = \frac{|U|}{2} \sum_{\mbf{R},\al}\bar{S}^z_{\mbf{R},\al}\bar{S}^z_{\mbf{R},\al}$ Hamiltonian studied in this work.

One may also calculate the commutators using the $\gamma$ basis. We define the momentum space spin operator

\begin{align}
\bar{S}^z_{\mbf{q},\al} &= \sum_\RR e^{-i\mbf{q}\cdot\RR} \bar{S}^z_{\RR,\alpha} \\
&= \sum_{\RR, \alpha,\sigma} (-1)^\sigma e^{-i\mbf{q}\cdot\RR} {\bar n}_{\RR,\alpha,\sigma} \\
&= \sum_{\kk, \alpha,\sigma} (-1)^\sigma {\bar c}^\dagger_{\kk-\qq,\alpha,\sigma} {\bar c}_{\kk,\alpha,\sigma} \\
&= \sum_{\kk, \alpha,\sigma} (-1)^\sigma [U^*(\kk-\qq)]_{\alpha,m,\sigma} [U(\kk)]_{\alpha,n,\sigma} \gamma^\dagger_{\kk-\qq,m,\sigma} \gamma_{\kk,n,\sigma}.
\label{}
\end{align}

The commutator becomes

\begin{align}
\null [\bar{S}^z_{\pp,\al},\eta^\dag] &= \sum_{\kk,m,n,\sigma} \sum_{\mbf{q},l} [(-1)^\sigma [U^*(\kk-\pp)]_{\alpha,m,\sigma} [U(\kk)]_{\alpha,n,\sigma} \gamma^\dagger_{\kk-\pp,m,\sigma} \gamma_{\kk,n,\sigma} ,\gamma^\dag_{\mbf{q},l,\u} \gamma^\dag_{-\mbf{q},l,\d} ] \\
&=  \sum_{\kk,m,n,\sigma} \sum_{\mbf{q},l} [U^*(\kk-\pp)]_{\alpha,m,\sigma} [U(\kk)]_{\alpha,n,\sigma} \left( \gamma^\dagger_{\kk-\pp,m,\u} \delta_{\kk,\mbf{q}} \delta_{ln} \delta_{\sigma,\u} \gamma^\dag_{-\mbf{q},l,\d} - \gamma^\dagger_{\kk-\pp,m,\d} \delta_{\kk,-\mbf{q}} \delta_{ln} \delta_{\sigma,\d}  \gamma^\dagger_{\mbf{q},l,\u} \right) \\
&= \sum_{\kk,m,n,\sigma} [U^*(\kk-\pp)]_{\alpha,m,\sigma} [U(\kk)]_{\alpha,n,\sigma} \left(  \delta_{\sigma,\u} \gamma^\dagger_{\kk-\pp,m,\u} \gamma^\dag_{-\kk,n,\d} -  \delta_{\sigma,\d} \gamma^\dagger_{\kk-\pp,m,\d} \gamma^\dagger_{-\kk,n,\u} \right) \\
&= \sum_{\kk,m,n,\sigma} \left( [U^*(\kk-\pp)]_{\u,m,\sigma} [U(\kk)]_{\u,n,\sigma} - [U^*(-\kk)]_{\d,n,\sigma} [U(\pp-\kk)]_{\d,m,\sigma} \right) \gamma^\dagger_{\kk-\pp,m,\u} \gamma^\dag_{-\kk,n,\d} = 0.
\label{}
\end{align}

§.§ Symmetry Algebra

The eta-pairing states $\eta^{\dag n}\ket{0}$ are naturally interpreted as zero-momentum states of charge $+2$ Cooper pairs, with $\eta^\dag$ being their creation operator. Because the states $\eta^{\dag n}\ket{0}$ are zero energy for each $n$, we can obtain eigenstates of $H$ by taking superpositions with different numbers of particles. This is the simply the Cooper pair condensate that appears in the mean-field description. However, we work at fixed particle number in this work.

We now show that the eta-pairing operators form an $su(2)$ charge algebra which extends the usual $u(1)$ charge algebra. We compute

[η, η^†]    = ∑_kq [c̅_-k,, c̅_k,, , c̅^†_q,, c̅^†_-q,,]

= ∑_kq δ_k,q c̅_-k,, P^_̆(k)c̅^†_-q,, - δ_k,q c̅^†_q,,  P^_̣(-k) c̅_k,,

=  ∑_k c̅_k,, P^_̆(-k)c̅^†_k,, - c̅^†_k,,  P^_̣(-k) c̅_k,,

= ∑_k P^_̆(-k) P^_̣(k) -c̅^†_k,, P^_̣(k) c̅_k,, - c̅^†_k,,  P^_̆(k) c̅_k,,

= ∑_k P^_̆(-k) P^_̆(-k) -c̅^†_k,, P^_̣(k) c̅_k,, - c̅^†_k,,  P^_̆(k) c̅_k,,

= ∑_k P^(̆-k)  - N̅

= 𝒩N_f - N̅
where we used that $P^\d_{\be \al}(\mbf{k})  = P^\u_{\al \be}(-\mbf{k})$ and $\sum_\be P_{\al \be}(\mbf{k}) \bar{c}_{\mbf{k},\be} = \bar{c}_{\mbf{k},\al}$. $\bar{N}$ is the total number operator so

[N̅, η^†] = 2 η^†,     [N̅, η] = -2 η .
Alternatively, in the $\gamma$ band basis,

\begin{align}
[\eta, \eta^\dag] &= \sum_{\mbf{q},m} \sum_{\mbf{q}',m'} [\gamma_{-\mbf{q}',m',\d} \gamma_{\mbf{q}',m',\u}, \gamma^\dag_{\mbf{q},m,\u} \gamma^\dag_{-\mbf{q},m,\d}] \\
&= \sum_{\mbf{q},m} \sum_{\mbf{q}',m'} \gamma_{-\mbf{q}',m',\d} \{\gamma_{\mbf{q}',m',\u}, \gamma^\dag_{\mbf{q},m,\u}\} \gamma^\dag_{-\mbf{q},m,\d} - \gamma^\dag_{\mbf{q},m,\u} \gamma_{\mbf{q}',m',\u} \{\gamma^\dag_{-\mbf{q},m,\d},\gamma_{-\mbf{q}',m',\d}\} \\
&= \sum_{\mbf{q},m} \sum_{\mbf{q}',m'} \gamma_{-\mbf{q}',m',\d} \delta_{mm'} \delta_{\mbf{q}\mbf{q}'} \gamma^\dag_{-\mbf{q},m,\d} - \gamma^\dag_{\mbf{q},m,\u} \gamma_{\mbf{q}',m',\u} \delta_{mm'} \delta_{\mbf{q}\mbf{q}'} \\
&= \sum_{\mbf{q},m} \gamma_{-\mbf{q},m,\d} \gamma^\dag_{-\mbf{q},m,\d} - \gamma^\dag_{\mbf{q},m,\u} \gamma_{\mbf{q},m,\u} = {\cal N} N_f - {\bar N}.
\label{}
\end{align}

To make the $su(2)$ Lie algebra explicit, the correct normalizations are
[η^a,η^b] = i _abc η^c,     η^x = (η+ η^†)/2,   η^y = i (η-η^†)/2,   η^z = (N̅ - 𝒩 N_f)/2  .

Thus we have shown that the $\eta$ operators form a representation of $SU(2)$. It is simple to find the dimension of the representation. For a dimension $d$ irrep of $SU(2)$, $\eta^{\dag d} = 0$. Indeed, we know that $\eta^{\dag (\mathcal{N} N_f +1)} = 0$ because the projected Hilbert space only contains $2 \mathcal{N} N_f$ electrons, and $\eta^\dag$ adds two electrons to the state. It is possible to use the explicit representations of $SU(2)$ to obtain normalized states $\ket{n} \propto \eta^{\dag n}\ket{0}$, but we will use a more direct method involving a generating function in the next section.

We now briefly address the spin symmetries. In the case of $S_z$ and $\mathcal{T}$, the many-body $S_z$ operator
S̅_z = ∑_k,n γ^†_k,n,γ_k,n, - γ^†_k,n,γ_k,n,
commutes with $H$. To see this, we use $[\bar{S}_z, \gamma^\dag_{\mbf{k},n,\sigma}] = s^z_\sigma \gamma^\dag_{\mbf{k},n,\sigma}$ (where $s^z_\sigma = \pm 1$ for $\sigma = \u/\d$) to check that  $[S̅_z, c̅^†_R,, σ] = s^z_σc̅^†_R,, σ$. Thus $[S̅_z, S̅_R,^z] = 0$ (so $[S̅_z, H] = 0$) and also $[S̅_z, η^†] = 0$. Thus we have a $U(1) ×SU(2)$ symmetry group. We emphasize that $H$ breaks spin $SU(2)$, having only $U(1)$ from conservation of $S_z$. One sign of this, which we will see later, is that the spin $±1$ excitations are gapped. \\

\section{Off-diagonal Long Range Order in the Ground State and Generating Functions}
\label{app:generatingfunctions}

In this section, we show that the eta-pairing ground states have ODLRO and thus display superconductivity at zero temperature. This calculation is tractable using a generating function of the BCS form (\App{app:odlro}). In an effort to be self-contained, we then show that Wick's theorem holds in the ground state (\App{app:wick}) using a similar generating function with Grassman variables.  Our calculations resemble previous calculations involving eta-pairing, for example \Refs{PhysRevLett.63.2144}{2017ScPP....3...43V}.

\subsection{ODLRO}
\label{app:odlro}

To assess the superconducting properties of the groundstate, we will compute the off-diagonals of the two-body density matrix:
\bea
\label{eq:odlroGS}
\lim_{|\mbf{R}- \mbf{R}'| \to \infty} \braket{GS| w^\dag_{\mbf{R}' m \d}w^\dag_{\mbf{R}' m \u} w_{\mbf{R}n \u} w_{\mbf{R} n \d} |GS}
\eea
in an eta-pairing groundstate $|GS⟩$. Here $w^†_R,n,σ$ is the Wannier function obtained from the canonical electron operators
\bea
\label{eq:Wannierspin}
w^\dag_{\mbf{R},n,\sigma} = \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k}} e^{i \mbf{R} \cdot \mbf{k}} \gamma^\dag_{\mbf{k},n,\sigma} .
\eea
Such operators can be chosen to be exponentially decaying in obstructed atomic limit bands, and also in fragile bands if we allow their symmetries to not be represented locally [Winkler} et~al., 2016]. This is because the only obstruction to exponential localization is Chern number [Brouder et~al., 2007]. In the Wannier basis (choosing a time-reversal symmetric gauge for the $γ$ operators where $U^(̆k)=U^(̣-k)^*$), the $η$ operator reads
\bea
\eta^\dag &= \sum_{\mbf{q}m} \gamma^\dag_{\mbf{q},m,\u} \gamma^\dag_{-\mbf{q},m,\d} = \sum_{\mbf{R},\mbf{R}',m}  w^\dag_{\mbf{R},m,\u} w^\dag_{\mbf{R}',m,\d} \frac{1}{\mathcal{N}}\sum_{\mbf{q}} e^{-i (\mbf{R}-\mbf{R}') \cdot \mbf{q}} = \sum_{\mbf{R}}  w^\dag_{\mbf{R},m,\u} w^\dag_{\mbf{R},m,\d} \ . \\
\eea
To calculate the correlator in \Eq{eq:odlroGS}, we will employ a BCS-type ground state as a generating function:
\bea
\label{eq:wangen}
\ket{\{z\}} &= \exp \lp \sum_{\mbf{R},m} z_{\mbf{R},m} w^\dag_{\mbf{R},m,\u} w^\dag_{\mbf{R},m,\d} \rp \ket{0}  \\
\eea
with free parameters $z_R,m$. First we observe that setting $z_R,m = z$ gives
\bea
\label{eq:genorm}
\braket{\{z\}=z|\{z\}=z} &=  \braket{0|e^{\bar{z} \eta} e^{z \eta^\dag}|0} = \sum_n \frac{|z|^{2n}}{n!^2} \braket{0|\eta^n \eta^{\dag n}|0} \\
\eea
which we will use to determine the norms of the eta-pairing states. Second, we note the ODLRO correlators appear when taking derivatives:
\bea
\label{eq:ODLROcorr}
\left. \frac{\del}{\del \bar{z}_{\mbf{R},n}}\frac{\del}{\del z_{\mbf{R}',m}} \braket{\{z\}|\{z\}} \right|_{z_{\mbf{R},m} = z} &=  \braket{0|e^{\bar{z} \eta} w_{\mbf{R},n,\d}  w_{\mbf{R},n,\u} w^\dag_{\mbf{R}',m,\u} w^\dag_{\mbf{R}',m,\d} e^{z \eta^\dag}|0} \\
&= \sum_n \frac{|z|^{2n}}{n!^2} \braket{0|\eta^n w_{\mbf{R},n,\d}  w_{\mbf{R},n,\u} w^\dag_{\mbf{R}',m,\u} w^\dag_{\mbf{R}',m,\d} \eta^{\dag n}|0} \ . \\
\eea
We now must compute $⟨{z}|{z}|$⟩. We use the orthogonality of the Wannier basis, $\{w_{\mbf{R},m,\sigma},w^\dag_{\mbf{R}',n,\sigma'} \}= \delta_{\mbf{R}\mbf{R}'} \delta_{mn} \delta_{\sigma,\sigma'}$ to factor eq:wangen into commuting terms
⟨{z}|{z}|   ⟩= ⟨0| ∏_Rm e^z̅_R,m w_R, m, w_R,m,  ∏_R'n e^z_R',n w^†_R',n, w^†_R', n, |0|
= ⟨0| ∏_Rm (1+z̅_R,m w_R, m, w_R,m,) ∏_R'n (1+ z_R',n w^†_R',n, w^†_R', n, ) |0|
where we used that fact that $(w^\dag_{\mbf{R}',n,\u} w^\dag_{\mbf{R}', n,\d})^2 = 0$. Commuting through the different factors and recalling $w^\dag_{\mbf{R},m}$ obey canonical commutation relations, we see that

⟨{z}|{z}|   ⟩= ⟨0| ∏_Rm (1+z̅_R,m w_R,m, w_R,m, )  (1+ z_R,m  w^†_R,m, w^†_R,m, ) |0|
= ∏_Rm (1+|z_R,m |^2)  .
Comparing to eq:genorm and using the binomial theorem, using eq:genorm we obtain
⟨{z}|{z}||⟩_{z}=z    = (1+|z|^2)^N_f 𝒩 = ∏_Rm (1+|z|^2) = (1+|z|^2)^N_f𝒩
∑_n |z|^2n/n!^2 ⟨0|η^n η^†n|0|   ⟩= ∑_n N_f 𝒩n |z|^2n
so equating coefficients, the normed eta-pairing states are

|n⟩    = 1/n! N_f𝒩n^-1/2 η^†n|0⟩  .
Then to obtain the ODLRO correlators, we observe that
. /z̅_R,n/z_R',m ⟨{z}|{z}|⟩|_{z} = z    = |z|^2 (1+|z|^2)^N_f𝒩-2 = ∑_n N_f 𝒩-2n |z|^2n+2,     R,n ≠R',m  .

Comparing to eq:ODLROcorr and using, we obtain
N_f𝒩n ⟨0|η^n w_R,n,  w_R,n, w^†_R',m, w^†_R',m, η^†n|0|   ⟩=  N_f 𝒩-2n-1  .
For $\mbf{R},n \neq \mbf{R}',m$, we commute $w_{\mbf{R},n,\d}  w_{\mbf{R},n,\u}$ through $ w^\dag_{\mbf{R}',m,\u} w^\dag_{\mbf{R}',m,\d} $ and obtain

⟨n| w^†_R m w^†_R m  w_R' n  w_R n'  |n|   ⟩=  (N_f 𝒩-2)!/(N_f𝒩-n-1)!(n-1)!  (N_f𝒩-n)!n!/(N_f 𝒩)!

=  n(N_f𝒩-n)/N_f 𝒩(N_f 𝒩-1)

=  ν(1-ν) + O(1/𝒩)

where we defined the filling $\nu = n/N_f\mathcal{N} \in (0,1)$. The other notable feature of this calculation is that for $\mbf{R} \neq \mbf{R}'$, there is no dependence on $|\mbf{R}-\mbf{R}'|$ in the correlator. This is reminiscent of the ODRLO in a $T=0$ Bose-Einstein condensate, which also has no $|\mbf{R}-\mbf{R}'|$ dependence and is equal to the boson density. Indeed, comparing to mean field calculations of the superfluid weight[Huhtinen et al., 2022], $\nu (1-\nu)$ is also obtained as the Cooper pair density.

For exponentially decaying Wannier states, eq:ODLROexp is a good local measure of off-diagonal long range order. The ODLRO has the interpretation of correlators between two particles at position $\RR$ and another two at position $\RR'$:  should the Wannier states fail to be localized eq:odlroGS loses this interpretation.  We thus must address the question: are the $w$ degrees of freedom employed in the ODLRO correlation exponentially localized?  In the uniform pairing, where bands are constructed from the $S$-matrix method, the bands are fragile or obstructed atomic and thus exponentially localized Wannier states can be found [Brouder et al., 2007], but the symmetries cannot be represented locally on them if there is fragile topology. We remark that it is possible to obtain perfectly flat bands with a nonzero Chern number $C_\sigma$ with uniform pairing from a single Wyckoff position sublattice, but this requires exponentially decaying (infinite range) hoppings. The total Chern number $C_\u + C_\d$ will vanish due to spin-ful time reversal symmetry, and thus it is possible to find exponentially localized Wannier states by mixing the bands of different spin. However, the calculation in this section relies on the spin-polarized Wannier states \Eq{eq:Wannierspin}, and thus the optimal decay of the Wannier states, for nonzero $C_σ$, is a power law [Monaco} et~al., 2016]. An adaptation of \Eq{eq:ODLROexp} in terms of exponentially localized states when $C_σ≠0$ is left for future work.

\subsection{Wick's theorem}
\label{app:wick}

We now extend the generating function method of \App{app:odlro} to show that a version of Wick's theorem holds in the eta -pairing states. In this section we work in the momentum basis using the $γ^†_k,m,σ$ operators. Wick's theorem shows that arbitrary correlators can be determined by enumerating the contractions, which allows us to compute expressions for the norms of the excitations in App.~\ref{app:excit} and calculate expectation values.

The generating function we introduce is the norm $N(z,ξ)$ of the following state:
\bea
\ket{\xi,z} = \exp \lp \sum_{\mbf{k}n\sigma} \xi_{\mbf{k},n,\sigma} \gamma^\dag_{\mbf{k}, n,\sigma}  \rp e^{z \eta^\dag} \ket{0}, \qquad N(z,\xi) = \braket{\xi,z|\xi,z} \\
\eea
where $ξ_k,n,σ$ are anti-commuting (fermionic) Grassman variables, and bar denotes the complex conjugate $ξ̅ = ξ^†$. Taking $ξ$ derivatives of $N(ξ,z)$ and then setting $ξ= 0$ pulls down $γ^†$ operators which build the correlation function. We observe that
\bea
\label{eq:Wicks}
\gamma^\dag_{\mbf{k}, n,\sigma} \ket{\xi,z} &= \del_{\xi_{\mbf{k},n,\sigma}} \exp \lp \sum_{\mbf{k}'n'\sigma'} \xi_{\mbf{k}',n',\sigma'} \gamma^\dag_{\mbf{k}', n',\sigma'}  \rp e^{z \eta^\dag} \ket{0} =  \del_{\xi_{\mbf{k},n,\sigma}} \ket{\xi,z} \\
\bra{\xi,z} \gamma_{\mbf{k}, n,\sigma} &= -\del_{\bar{\xi}_{\mbf{k},n,\sigma}} \bra{0} e^{\bar{z} \eta} \exp \lp \sum_{\mbf{k}'n'\sigma'} \gamma_{\mbf{k}', n',\sigma'}  \bar{\xi}_{\mbf{k}',n',\sigma'}   \rp  = -\del_{\bar{\xi}_{\mbf{k},n,\sigma}}  \bra{\xi,z} \ .
\eea
We see for instance that
\bea
\label{eq:deldelxixix}
\bra{\xi,z} \gamma_{\mbf{k}', n',\sigma'} \gamma^\dag_{\mbf{k}, n,\sigma} \ket{\xi,z} &=  -\del_{\bar{\xi}_{\mbf{k}',n',\sigma'}}  \bra{\xi,z}  \, \del_{\xi_{\mbf{k},n,\sigma}} \ket{\xi,z} =  \del_{\xi_{\mbf{k},n,\sigma}} \del_{\bar{\xi}_{\mbf{k}',n',\sigma'}}  N(z,\xi) \\
\bra{z} \gamma_{\mbf{k}', n',\sigma'} \gamma^\dag_{\mbf{k}, n,\sigma} \ket{z} &= \left. \del_{\xi_{\mbf{k},n,\sigma}} \del_{\bar{\xi}_{\mbf{k}',n',\sigma'}}  N(z,\xi) \right|_{\xi = 0}\\
\eea
and thus by matching powers of $|z|^2$ in $_ξ_km σ _ξ̅_km σ N(z, ξ)$ and $ ⟨z|γ_k, m,σ γ^†_k m σ |z⟩$, the correlators $ ⟨n|γ_k, m,σ γ^†_k m σ |n⟩$ can be determined. Generally, any correlator with creation operators to the right and annihilation operators to the left can be reduced by repeated application of \Eq{eq:Wicks} to derivatives of $⟨ξ,z|ξ,z|=⟩ N(ξ,z)$ (and then setting $ξ→0$). This defines a normal ordering.

We can compute $N(z, ξ)$ directly. Because $ξ$ and $γ^†$ are all anti-commuting, we have
\bea
\ket{\xi,z} &= \exp \lp \sum_{\mbf{k}n \sigma} \xi_{\mbf{k},n,\sigma} \gamma^\dag_{\mbf{k},n,\sigma} \rp e^{z \eta^\dag} \ket{0} \\
&= \exp \lp \sum_{\mbf{k}n} \xi_{\mbf{k},n,\u} \gamma^\dag_{\mbf{k}, n,\u} + \xi_{-\mbf{k},n,\d} \gamma^\dag_{-\mbf{k}, n,\d} + z \gamma^\dag_{\mbf{k},n,\u} \gamma^\dag_{-\mbf{k},n,\d} \rp \ket{0} \\
&= \prod_{\mbf{k},n} \exp \lp \xi_{\mbf{k},n,\u} \gamma^\dag_{\mbf{k},n,\u}+\xi_{-\mbf{k},n,\d} \gamma^\dag_{-\mbf{k},n,\d} + z \gamma^\dag_{\mbf{k},n,\u} \gamma^\dag_{-\mbf{k},n,\d} \rp \ket{0} \\
&= \prod_{\mbf{k},n}  \big( 1 + \xi_{\mbf{k},n,\u} \gamma^\dag_{\mbf{k}, n,\u}+ \xi_{-\mbf{k},n,\d} \gamma^\dag_{-\mbf{k}, n,\d} + (z+\xi_{-\mbf{k},n,\d} \xi_{\mbf{k},n,\u}) \gamma^\dag_{\mbf{k},n,\u} \gamma^\dag_{-\mbf{k},n,\d}  \big) \ket{0} \\
\eea
using $ξ^2 = γ^†2 = 0$ and $[ξγ^†, η^†]=0$.  The $ξ$ operators must be Grassmanns, else $[ξ_,n,↑ γ^†_,n,↑, ξ_',n',↑ γ^†_',n',↑] ≠0$.  The only term at second order is the cross term from $1/2!(ξ_-k,n, γ^†_-k, n,+γ^†_k, n, ξ_k,n,)^2$. Hence we find
\bea
\label{eq:Nzxiprod}
N(z, \xi) = \braket{\xi,z|\xi,z} &= \prod_{\mbf{k},n} \lp 1+ \bar{\xi}_{\mbf{k},n,\u}  \xi_{\mbf{k},n,\u}  + \bar{\xi}_{-\mbf{k},n,\d}  \xi_{-\mbf{k},n,\d}  +  (\bar{z}+\bar{\xi}_{\mbf{k},n,\u} \bar{\xi}_{-\mbf{k},n,\d} )  (z+\xi_{-\mbf{k},n,\d} \xi_{\mbf{k},n,\u}) \rp  \\
%&= \prod_{\mbf{k},n} \exp \lp \bar{\xi}_{\mbf{k},n,\u}  \xi_{\mbf{k},n,\u} + \bar{\xi}_{-\mbf{k},n,\d}  \xi_{-\mbf{k},n,\d} + z \bar{\xi}_{\mbf{k},n,\u} \bar{\xi}_{-\mbf{k},n,\d} + \bar{z} \xi_{-\mbf{k},n,\d} \xi_{\mbf{k},n,\u} \rp
\eea
which importantly is decoupled in $k$ and $n$. Let us expand out a single factor of \Eq{eq:Nzxiprod}:
\bea
1+ |z|^2  + \bar{\xi}_{\mbf{k},n,\u}  \xi_{\mbf{k},n,\u}  + \bar{\xi}_{-\mbf{k},n,\d}  \xi_{-\mbf{k},n,\d}  + z \bar{\xi}_{\mbf{k},n,\u} \bar{\xi}_{-\mbf{k},n,\d} +   \bar{z} \xi_{-\mbf{k},n,\d} \xi_{\mbf{k},n,\u} +
\bar{\xi}_{\mbf{k},n,\u} \bar{\xi}_{-\mbf{k},n,\d} \xi_{-\mbf{k},n,\d} \xi_{\mbf{k},n,\u}
\eea
which contains only quadratic and quartic terms in $ξ$. Because we take $ξ$ derivatives and then set $ξ→0$, we see that only two- and four- term derivatives (for each $k,m$) are nonzero. From \Eq{eq:Wicks}, we compute
\bea
\label{eq:cor1}
\braket{z|\gamma_{\mbf{k}'m'\sigma'}\gamma^\dag_{\mbf{k}m\sigma} |z} = \del_{\xi_{\mbf{k},m,\sigma}} \del_{\bar{\xi}_{\mbf{k}',m',\sigma'}}  N  |_{\xi=0}&= \frac{\delta_{\mbf{k}\mbf{k}'}\delta_{mm'} \delta_{\sigma \sigma'}}{1+|z|^2} \braket{z|z} \\
\braket{z|\gamma_{\mbf{k}m\d}\gamma_{\mbf{k}'m'\u} |z} =  \del_{\bar{\xi}_{\mbf{k},m,\d}}  \del_{\bar{\xi}_{\mbf{k}',m',\u} }N  |_{\xi=0}&= \frac{z \delta_{\mbf{k},-\mbf{k}'} \delta_{mm'} }{1+|z|^2} \braket{z|z} \\
\braket{z|\gamma^\dag_{\mbf{k}m\u} \gamma^\dag_{\mbf{k}'m'\d} |z}  =\del_{\xi_{\mbf{k},m,\u} }\del_{\xi_{\mbf{k}',m',\d}}  N  |_{\xi=0}&= \frac{\bar{z} \delta_{\mbf{k},-\mbf{k}'} \delta_{mm'}}{1+|z|^2} \braket{z|z} \\
\eea
with all other two derivative terms being zero, and $⟨z|z|=⟩ (1+|z|^2)^N_f 𝒩$. We observe in the second and third line a nonzero anomalous BCS-type contraction between $γγ$ or $γ^†γ^†$. The last nonzero term comes from the four-body derivative:
\bea
\label{eq:cor2}
\braket{z|\gamma_{-\mbf{k}m\d}\gamma_{\mbf{k}m\u} \gamma^\dag_{\mbf{k}m\u} \gamma^\dag_{-\mbf{k}m\d}  |z} =  \del_{\xi_{\mbf{k},m,\u}} \del_{\xi_{-\mbf{k},m,\d}}  \del_{\bar{\xi}_{-\mbf{k},m,\d}} \del_{\bar{\xi}_{\mbf{k},m,\u}}  N  |_{\xi=0}&= \frac{1}{1+|z|^2} \braket{z|z} \ . \\
\eea

So far we have given a recipe to compute all normal ordered correlators. Using the canonical anti-commutation relations of $γ$, all correlators can be manipulated into normal order and then evaluated. However, it is of practical interest to evaluate correlators directly in any order. We now calculate the one-body contractions, \Eqs{eq:cor1}{eq:cor2}, in which a $γ^†γ$ contraction is introduced to be consistent with the fermion anti-commutation relations. This is a standard result (see \Ref{1971qtmp.book.....F}) but we include a proof here to set the notation and be self-contained.  All many-body contractions will be built from these one-body contractions.

We compute the non-normal ordered correlator
\bea
\label{eq:cor3}
\braket{z|\gamma^\dag_{\mbf{k}m\sigma} \gamma_{\mbf{k}'m'\sigma'} |z} &= -\braket{z|\gamma_{\mbf{k}'m'\sigma'} \gamma^\dag_{\mbf{k}m\sigma}  |z} + \delta_{\mbf{k}\mbf{k}'} \delta_{mm'} \delta_{\sigma \sigma'}\braket{z|z} \\
&= \delta_{\mbf{k}\mbf{k}'} \delta_{mm'} \delta_{\sigma \sigma'} \braket{z|z} \lp - \frac{1}{1+|z|^2} + 1 \rp \\
&=  \delta_{\mbf{k}\mbf{k}'} \delta_{mm'} \delta_{\sigma \sigma'} \frac{|z|^2}{1+|z|^2} \braket{z|z} \ . \\
\eea
This correlator will define a two-particle contraction to be included in Wick's theorem:  correlators are equal to the signed sum of all possible two-particle contractions defined in \Eqs{eq:cor1}{eq:cor3}. Because only two- and four-particle contractions are nonzero and the two-particle contractions satisfy Wick's theorem automatically, we only need to check that \Eq{eq:cor2} is consistent with this Wick's theorem. Consider the only non-zero normal ordered four-particle correlator $⟨z|γ_-kmγ_km γ^†_km γ^†_-km  |z|$⟩. There are three possible contractions from Wick's theorem (one vanishing):
⟨z|γ_-kmγ_km γ^†_km γ^†_-km  |z|⟩⟨z|z|⟩    = ⟨z|γ_-kmγ_km|z|⟩⟨z|z|⟩  ⟨z| γ^†_km γ^†_-km|z|⟩⟨z|z|⟩  - ⟨z|γ_-km γ^†_km|z|⟩⟨z|z|⟩  ⟨z|γ_km γ^†_-km|z|⟩⟨z|z|⟩

+ ⟨z|γ_-km γ^†_-km|z|⟩⟨z|z|⟩  ⟨z|γ_km γ^†_km|z|⟩⟨z|z|⟩

= z/1+|z|^2 z̅/1+|z|^2-  (0) (0) + 1/1+|z|^2 1/1+|z|^2

=  |z|^2+1/(1+|z|^2)^2 =  1/1+|z|^2
matching eq:cor2.

For clarity, we list the nonzero elementary contractions of Wick's theorem:

1/⟨z|z|⟩ ⟨z|γ_k'm'σ'γ^†_kmσ |z|   ⟩= 1/1+|z|^2 δ_kk'δ_mm' δ_σσ'

1/⟨z|z|⟩⟨z|γ^†_kmσ γ_k'm'σ' |z|   ⟩= |z|^2/1+|z|^2 δ_kk'δ_mm' δ_σσ'

1/⟨z|z|⟩ ⟨z|γ_kmσγ_k'm'σ' |z|   ⟩= z/1+|z|^2 δ_k,-k'δ_mm' s_z^σ' δ_σ,-σ'

1/⟨z|z|⟩ ⟨z|γ^†_kmσ γ^†_k'm'σ' |z|   ⟩= z̅/1+|z|^2 δ_k,-k'δ_mm' s_z^σ δ_σ,-σ'

and $\braket{z|z} = (1+|z|^2)^{N_f \mathcal{N}}$, $s_z^{\u/\d} = \pm 1$. We will use these results in app:Richcrit to evaluate Richardson's criterion.

§ CHARGE EXCITATIONS

In this section, we compute the spectrum and eigenstates of excitations above the eta-pairing ground states. We show that the fermionic excitations are gapped and perfectly flat (app:charge1), while the bosonic excitations are gapless and dispersive. We evaluate Richardson's criterion to show that kinetic energy generically favors superconductivity at half filling (app:Richcrit). We then study the bosonic (charge +2 and 0) excitations and find that the spin-0 Cooper pair spectrum splits into a set of low-lying bands described exactly by a single-particle pairing Hamiltonian and a flat continuum at the two-electron gap (app:charge2). The spin $\pm 1$ Cooper pair spectrum is entirely flat: all states are at the two-electron gap and are effectively unpaired (app:charge2). We provide explicit examples of the pairing Hamiltonian (app:examples) for obstructed atomic and fragile bands. Finally, we show that the spin-0 density mode (Goldstone mode) has the same spectrum as the spin-0 Cooper pairs. The spin $\pm1$ density excitations are all energetically above the two-electron gap, but have nontrivial dispersion also described by a single-particle Hamiltonian (app:goldstones).

§.§ Charge +1

Given an eta-pairing pair groundstate $\ket{GS}=\ket{n}$, we can compute the charge $+1$, spin $\pm \frac{1}{2}$ excitations from the effective Hamiltonian $R(\mbf{k})$ defined

[H, γ^†_k,n,σ] |GS⟩    =  ∑_m γ^†_k,m,σ R^σ_mn(k)  |GS⟩  .
By diagonalizing $R^\sigma(\mbf{k})$, we obtain exact eigenstates and their energies. Note that $\sigma$ is not summed over. The excitation matrix is diagonal in spin because of $S_z$ conservation.

It is simplest to perform the excitation calculations in momentum space using the canonical $\gamma^\dag_{\mbf{k},n,\sigma}$ operators, which form a complete basis of the projected Hilbert space. First we define momentum space spin operators $\bar{S}^z_{\mbf{R},\al} = \sum_\sigma s^z_\sigma \bar{n}_{\mbf{R},\al,\sigma}$ (recall $s^z_\sigma = \pm1$ for $\sigma = \u,\d$) according to
\bea
\bar{S}^z_{\mbf{k},\al} &= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{R}} e^{-i \mbf{k} \cdot (\mbf{R}+\mbf{r}_\al)} \bar{S}^z_{\mbf{R},\al}, \qquad \bar{S}^z_{\mbf{k},\al} = (\bar{S}^z_{-\mbf{k},\al})^\dag, \qquad \bar{S}^z_{\mbf{R},\al} = \frac{1}{\sqrt{\mathcal{N}}} \sum_\mbf{k} e^{i \mbf{k} \cdot (\mbf{R}+\mbf{r}_\al)}  \bar{S}^z_{\mbf{k},\al} \\
\eea
such that
\bea
H &= \frac{1}{2}|U| \sum_{\mbf{R},\al} (\bar{S}^z_{\mbf{R},\al})^2 = \frac{1}{2}|U| \sum_\mbf{k, \al} \bar{S}^z_{-\mbf{k},\al} \bar{S}^z_{\mbf{k},\al} \ . \\
\eea
$H$ can be written in terms of the $γ$ modes using
\bea
\bar{n}_{\mbf{R},\al,\sigma} &= \bar{c}^\dag_{\mbf{R},\al,\sigma}\bar{c}_{\mbf{R},\al,\sigma} = \frac{1}{\mathcal{N}} \sum_{\mbf{k},\mbf{k}'} e^{- i (\mbf{k} - \mbf{k}') \cdot (\mbf{R}+\mbf{r}_\al)} \bar{c}^\dag_{\mbf{k}',\al,\sigma}  \bar{c}_{\mbf{k},\al,\sigma} \\
&= \frac{1}{\mathcal{N}} \sum_{\mbf{k},\mbf{k}'mn} e^{- i (\mbf{k} - \mbf{k}') \cdot (\mbf{R}+\mbf{r}_\al)} \gamma^\dag_{\mbf{k}',m,\sigma} [U^\sigma_{\al m}(\mbf{k}')^*  U_{\al n}^\sigma(\mbf{k})]\gamma_{\mbf{k},n,\sigma} \\
&\equiv \frac{1}{\mathcal{N}} \sum_{\mbf{k},\mbf{k}',mn} e^{- i (\mbf{k} - \mbf{k}') \cdot (\mbf{R}+\mbf{r}_\al)} \gamma^\dag_{\mbf{k}',m,\sigma} M^{mn}_{\sigma,\al}(\mbf{k},\mbf{k}'-\mbf{k}) \gamma_{\mbf{k},n,\sigma} \\
\eea
where we defined the sublattice form factor
\bea
M^{mn}_{\sigma,\al}(\mbf{k},\mbf{q}) &= U^\sigma_{\al m}(\mbf{k}+\mbf{q})^*  U_{\al n}^\sigma(\mbf{k}) \ .
\eea
Thus the Fourier transform of the spin operator is
\bea
\bar{S}^z_{\mbf{q},\al} &= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{R}} e^{-i \mbf{q} \cdot (\mbf{R}+\mbf{r}_\al)} \bar{S}^z_{\mbf{R},\al} \\
&= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{R}mn} e^{-i \mbf{q} \cdot (\mbf{R}+\mbf{r}_\al)}  \frac{1}{\mathcal{N}} \sum_{\mbf{k},\mbf{k}',\sigma} e^{- i (\mbf{k} - \mbf{k}') \cdot (\mbf{R}+\mbf{r}_\al)} s^z_\sigma \gamma^\dag_{\mbf{k}',m,\sigma} M^{mn}_{\sigma,\al}(\mbf{k},\mbf{k}'-\mbf{k})\gamma_{\mbf{k},n,\sigma} \\
&= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k},\mbf{k}',mn,\sigma}  \frac{1}{\mathcal{N}} \sum_{\mbf{R}}  e^{- i (\mbf{k} +\mbf{q} - \mbf{k}') \cdot (\mbf{R}+\mbf{r}_\al)} s^z_\sigma \gamma^\dag_{\mbf{k}',m,\sigma} M^{mn}_{\sigma,\al}(\mbf{k},\mbf{k}'-\mbf{k})\gamma_{\mbf{k},n,\sigma} \\
&= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k},mn,\sigma} s^z_\sigma \gamma^\dag_{\mbf{q}+\mbf{k},m,\sigma} M^{mn}_{\sigma,\al}(\mbf{k},\mbf{q})\gamma_{\mbf{k},n,\sigma} \\
\eea
and hence it follows that
\bea
\label{eq:formfactorS}
\null [\bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}] &=  \frac{1}{\sqrt{\mathcal{N}}} \sum_m s^z_\sigma \gamma^\dag_{\mbf{k}+\mbf{q},m,\sigma} M^{mn}_{\sigma,\al}(\mbf{k},\mbf{q}) \ .
\eea
To prepare to compute the excitation spectrum, we observe
\bea
\label{eq:SSgamma}
\sum_{\mbf{q}\al} [\bar{S}^z_{-\mbf{q},\al} \bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}] &= \sum_{\mbf{q}\al}   \bar{S}^z_{-\mbf{q},\al}  [ \bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}] +[\bar{S}^z_{-\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}]  \bar{S}^z_{\mbf{q},\al} \\
&= \sum_{\mbf{q}\al}   \bar{S}^z_{-\mbf{q},\al}  [ \bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}] +[\bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}]  \bar{S}^z_{-\mbf{q},\al} \\
&= \sum_{\mbf{q}\al}   [\bar{S}^z_{-\mbf{q},\al} , [ \bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}]] +2 [\bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}]  \bar{S}^z_{-\mbf{q},\al} \\
\eea
Acting on the groundstate where the second term is zero (recall that $S^z_q,$ annihilates $|n⟩$ because $S^z_q, |0⟩ = 0$ and $[S^z_q,,η^†]=0$) and using \Eq{eq:formfactorS}, we derive
\bea
\sum_{\mbf{q}\al} [\bar{S}^z_{-\mbf{q},\al} \bar{S}^z_{\mbf{q},\al}, \gamma^\dag_{\mbf{k}, n, \sigma}] \ket{GS} &=  \frac{1}{\sqrt{\mathcal{N}}} s^z_\sigma \sum_{\mbf{q}\al m'} [\bar{S}^z_{-\mbf{q},\al} , \gamma^\dag_{\mbf{k}+\mbf{q},m',\sigma}] M^{m'n}_{\sigma,\al}(\mbf{k},\mbf{q}) \ket{GS} \\
&=  \frac{1}{\mathcal{N}}\sum_{\mbf{q}\al m m'}  \gamma^\dag_{\mbf{k},m,\sigma}  M^{mm'}_{\sigma,\al}(\mbf{k}+\mbf{q},-\mbf{q}) M^{m'n}_{\sigma,\al}(\mbf{k},\mbf{q}) \ket{GS} \\
\eea
where we used that $(s^z_σ)^2 = 1$. Thus we obtain
\bea
\label{eq:Rsteps}
R^\sigma_{mn}(\mbf{k}) &=  \frac{|U|}{2} \frac{1}{\mathcal{N}} \sum_{\mbf{q} \al m'} M^{mm'}_{\sigma,\al}(\mbf{k}+\mbf{q},-\mbf{q}) M^{m'n}_{\sigma,\al}(\mbf{k},\mbf{q}) \\
&=  \frac{|U|}{2} \frac{1}{\mathcal{N}} \sum_{\mbf{q} \al m'} U^\sigma_{\al m}(\mbf{k})^*  U_{\al m'}^\sigma(\mbf{k}+\mbf{q})  U^\sigma_{\al m'}(\mbf{k}+\mbf{q})^*  U_{\al n}^\sigma(\mbf{k})  \\
&=  \frac{|U|}{2}  \sum_{\al} U^\sigma_{\al m}(\mbf{k})^* \frac{1}{\mathcal{N}} \sum_\mbf{q} P^\sigma_{\al \al}(\mbf{k}+\mbf{q})  U_{\al n}^\sigma(\mbf{k})  \\
\eea
and using the uniform pairing condition \Eq{eq:UPCint}, $1/𝒩 ∑_q P^σ_(k+q) = ϵ$, we have
\bea
\label{eq:Rmn}
R^\sigma_{mn}(\mbf{k})  &=  \frac{|U|}{2}  \sum_{\al} U^\sigma_{\al m}(\mbf{k})^* \eps U_{\al n}^\sigma(\mbf{k}) \\
&=  \frac{\eps |U|}{2}  [U_\sigma^\dag(\mbf{k}) U_\sigma(\mbf{k})]_{mn} \\
&=  \frac{\eps |U|}{2} \delta_{mn} \\
\eea
which is diagonal for all $k$ and has completely flat bands with an energy gap $|U|/2$ above the ground state. \Eq{eq:charge1excit} simply reads
\bea
\label{eq:estateeq}
\null H \gamma^\dag_{\mbf{k},n,\sigma} \ket{GS} &=  \frac{\eps|U|}{2} \gamma^\dag_{\mbf{k},n,\sigma} \ket{GS}
\eea
making use of $H|GS⟩ = 0$. We refer to $|U|/2$ as the single-electron gap, whose analogue in BCS theory is the Bogoliubov quasi-particle gap.

We can repeat this calculation for the hole excitations, defined by
\bea
\null [H, \gamma_{\mbf{k},n,\sigma}] \ket{GS} &=  \sum_m \tilde{R}^\sigma_{nm}(\mbf{k}) \gamma_{\mbf{k},m,\sigma}   \ket{GS} \ .
\eea
Taking the hermitian conjugate of \Eq{eq:formfactorS}, we obtain the requisite formula
\bea
\null [\bar{S}^z_{\mbf{q},\al}, \gamma_{\mbf{k}, n, \sigma}] &=  -\frac{1}{\sqrt{\mathcal{N}}} \sum_m s^z_\sigma \gamma_{\mbf{k}-\mbf{q},m,\sigma} M^{mn}_{\sigma,\al}(\mbf{k},\mbf{q})^* \\
&=  -\frac{1}{\sqrt{\mathcal{N}}} \sum_m M^{nm}_{\sigma,\al}(\mbf{k}+\mbf{q},-\mbf{q}) s^z_\sigma \gamma_{\mbf{k}-\mbf{q},m,\sigma} \ . \\
\eea
All steps in the calculation of $R^σ_mn(k)$ proceed identically for $R̃^σ_nm(k)$ but with the mapping $M^mn_σ,(k,q) →M^nm_σ,(k+q,-q)$, so from \Eq{eq:Rsteps} we find
\bea
\tilde{R}^\sigma_{nm}(\mbf{k}) &= \frac{|U|}{2} \frac{1}{\mathcal{N}} \sum_{\mbf{q} \al m'} M^{m'm}_{\sigma,\al}(\mbf{k},\mbf{q}) M^{nm'}_{\sigma,\al}(\mbf{k}+\mbf{q},-\mbf{q})  \\
&= \frac{|U|}{2} \frac{1}{\mathcal{N}} \sum_{\mbf{q} \al m'} M^{nm'}_{\sigma,\al}(\mbf{k}+\mbf{q},-\mbf{q}) M^{m'm}_{\sigma,\al}(\mbf{k},\mbf{q}) \\
&= R^\sigma_{nm}(\mbf{k}) \\
\eea
and hence single particle and hole excitations have the same energies.  This is not true for systems that do not satisfy the uniform pairing condition.  For example, the flat band systems in magic angle twisted bilayer graphene have different particle and hole excitations [Bernevig} et~al., 2021].

\subsection{Richardson Criterion}
\label{app:Richcrit}

The Richardson criterion[Richardson} and {Sherman}, 1964, Richardson}, 1977] is an approximation of the optimal binding energy of the Cooper pair, defined by $E_Δ(N) = E(N+2) - E(N)   - 2(E(N+1)-E(N))$  where $E(N)$ is the groundstate energy at particle number $N$.  In the flat band case, $E(2n) = 0   ∀n$ due to the eta-pairing symmetry because $|n⟩$ are zero energy ground states. Thus $E_Δ(N) = - |U| < 0$ demonstrating perfect pairing at all densities. Indeed, we will see in \App{app:Cooperpairspectrum} that $|U|$ is the maximum binding energy computed from the charge $+2$ Cooper pair excitations.

When weak dispersion (kinetic energy) is added to the model, we will determine which density shows the strongest pairing: that is, the most negative binding energy $E_Δ(N)$. To evaluate $E(N),E(N+2)$ at $N=2n$ in the eta-pairing groundstates, we need the expectation value of the kinetic energy in the states $|n⟩$, and for $E(N+1)$ we need to compute the expectation value in the charge +1 states atop the $N = 2n$ ground state, which is a four operator correlator and can be computed with Wick's theorem.

We consider a more general perturbation than discussed in the Main Text. We consider a perturbation to the kinetic energy defined in \Eq{eq:gammasing}, namely
\bea
\tilde{H} + \tilde{H}' &= \sum_{\mbf{k},\al \be,\sigma}c^\dag_{\mbf{k},\al,\sigma} \tilde{h}^\sigma_{\al \be}(\mbf{k}) c_{\mbf{k},\be,\sigma} + \sum_{\mbf{k},\al \be,\sigma\sigma'}c^\dag_{\mbf{k},\al,\sigma} \tilde{h}'_{\al \sigma, \be\sigma'}(\mbf{k}) c_{\mbf{k},\be,\sigma'} \\
\eea
where $h̃_(k)$ is the flat band Hamiltonian and we treat $h̃'_σ,σ'(k)$ as a small perturbation (additional hoppings, spin-orbit couplings to break $𝒯$ and $S_z$, etc). When projected into the flat bands (still requiring that the conduction bands are much higher in energy than the scale of $h̃'$ and $|U|$), we find
\bea
\tilde{H}' &\to \sum_{\mbf{k},\al \be,\sigma\sigma'}\bar{c}^\dag_{\mbf{k},\al,\sigma} \tilde{h}'_{\al \sigma, \be\sigma'}(\mbf{k}) \bar{c}_{\mbf{k},\be,\sigma'}\\
&= \sum_{\mbf{k},\al\al' \be\be',\sigma\sigma'}\bar{c}^\dag_{\mbf{k},\al',\sigma} P^{\sigma}_{\al' \al}(\mbf{k})  \tilde{h}'_{\al \sigma, \be\sigma'}(\mbf{k}) P^{\sigma'}_{\be \be'}(\mbf{k}) \bar{c}_{\mbf{k},\be',\sigma'}\\
&= \sum_{\mbf{k},\al\be,\sigma\sigma',mn}\gamma^\dag_{\mbf{k},m,\sigma} U^{\sigma*}_{\al m}(\mbf{k}) \tilde{h}'_{\al \sigma, \be\sigma'}(\mbf{k}) U^{\sigma'}_{\be n}(\mbf{k}) \gamma_{\mbf{k},n,\sigma'} \\
\eea
using \Eq{eq:bargamma} to write the $c̅$ operators in terms of $γ$ operators. Thus the correction to the flat band kinetic energy Hamiltonian takes the form