\bea
\tilde{H}' = \sum_{\mbf{k}mn\sigma\sigma'} \tilde{E}'_{m\sigma,n\sigma'}(\mbf{k}) \gamma^\dag_{\mbf{k},m,\sigma}\gamma_{\mbf{k},n,\sigma'}, \qquad  \tilde{E}'_{m\sigma,n\sigma'}(\mbf{k})  = \sum_{\al\be}U^{\sigma*}_{\al m}(\mbf{k}) \tilde{h}'_{\al \sigma, \be\sigma'}(\mbf{k}) U^{\sigma'}_{\be n}(\mbf{k}) \\
\eea
or in matrix form, recalling $U_σ(k)$ is an $N_orb ×N_f$ rectangular matrix,
\bea
\tilde{E}'(\mbf{k}) = \bpm U^\dag_\u(\mbf{k}) & \\
& U^\dag_\d(\mbf{k}) \epm \bpm
\tilde{h}'_{\u\u}(\mbf{k}) & \tilde{h}'_{\u\d}(\mbf{k}) \\
\tilde{h}'_{\d\u}(\mbf{k}) & \tilde{h}'_{\d\d}(\mbf{k}) \\
\epm \bpm U_\u(\mbf{k}) & \\
& U_\d(\mbf{k}) \epm
\eea
in the spin $⊗$ orbital tensor product basis. We develop an expression for Richardson's criterion in the general case, but it may be helpful to point out the simplest case: dispersion $Ẽ'_m,σ(k)$ can be added to the flat bands by taking $ Ẽ'_mσ,nσ'(k) = δ_nm δ_σσ' Ẽ'_m,σ(k)$.

We now need to compute expectation values. The ground states at $N=2n$ and $N+2 = 2(n+1)$ are $|n⟩, |n+1⟩$ at zeroth order in $Ẽ'$. When $H̃'$ is added, we determine the shift in energy with degenerate perturbation theory; by computing the matrix elements of $H̃'$ with respect to the unperturbed states:
\bea
\braket{z|\tilde{H}'|z} &= \sum_{\mbf{k}mn\sigma\sigma'} \tilde{E}'_{m\sigma,n\sigma'}(\mbf{k})  \braket{z|\gamma^\dag_{\mbf{k},m,\sigma}\gamma_{\mbf{k},n,\sigma'} |z} \\
&= \sum_{\mbf{k}mn\sigma\sigma'} \tilde{E}'_{m\sigma,n\sigma'}(\mbf{k}) \delta_{mn} \delta_{\sigma\sigma'} \frac{|z|^2}{1+|z|^2} \braket{z|z} \\
&= |z|^2 (1+|z|^2)^{N_f\mathcal{N}-1} 2N_f \mathcal{N} \tilde{E}', \qquad 2N_f \mathcal{N} \tilde{E}' = \sum_{\mbf{k}m\sigma} \tilde{E}'_{m\sigma,m\sigma}(\mbf{k}) \\
\eea
using \Eq{eq:wickcontractions}. We defined $Ẽ' = 1/2N_f 𝒩 ∑_kmσ Ẽ'_mσ,mσ$ as the average of the entire spectrum, which is finite in the thermodynamic limit. Expanding both sides in powers of $z$, we find
\bea
\label{eq:Richarson1}
\sum_n |z|^{2n} \binom{N_f\mathcal{N}}{n} \braket{n|\tilde{H}'|n} &= \sum_n |z|^{2(n+1)} \binom{N_f\mathcal{N}-1}{n} (2N_f \mathcal{N} \tilde{E}') \\
\binom{N_f\mathcal{N}}{n} \braket{n|\tilde{H}|n} &= \binom{N_f\mathcal{N}-1}{n-1}(2N_f \mathcal{N} \tilde{E}') \\
\braket{n|\tilde{H}|n} &= \frac{n}{N_f\mathcal{N}} (2N_f \mathcal{N} \tilde{E}') = \nu (2N_f \mathcal{N} \tilde{E}')
\eea
recalling $ n/N_f𝒩 = ν$ is the filling. Next we need the expectation value of $H̃'$ in the electron excitations \Eq{eq:estateeq}, but first we must normalize them. We compute
\bea
\label{eq:1pnorm}
\braket{z|\gamma_{\mbf{k},m,\sigma} \gamma^\dag_{\mbf{k},m,\sigma}|z} &= \frac{1}{1+|z|^2} \braket{z|z} \\
\binom{N_f \mathcal{N}}{n} \braket{n|\gamma_{\mbf{k},m,\sigma} \gamma^\dag_{\mbf{k},m,\sigma}|n} &= \binom{N_f \mathcal{N} -1}{n} \\
\braket{n|\gamma_{\mbf{k},m,\sigma} \gamma^\dag_{\mbf{k},m,\sigma}|n} &= 1-\nu \ .\\
\eea
Because all $γ^†_k,m,σ|n⟩$ excitations are degenerate, we compute their effective Hamiltonian at first order in degenerate perturbation theory:
\bea
\braket{z|\gamma_{\mbf{k},m',\sigma'} \tilde{H} \gamma^\dag_{\mbf{k},m,\sigma}|z} &= \sum_{\mbf{k}'ijss'} \tilde{E}'_{is,js'}(\mbf{k}')  \braket{z|\gamma_{\mbf{k},m',\sigma'} \gamma^\dag_{\mbf{k}',i,s}\gamma_{\mbf{k}',j,s'} \gamma^\dag_{\mbf{k},m,\sigma}|z} \ . \\
\eea
There are only three contractions:
\bea
\braket{z|\gamma_{\mbf{k},m',\sigma'} \gamma^\dag_{\mbf{k}',i,s}\gamma_{\mbf{k}',j,s'} \gamma^\dag_{\mbf{k},m,\sigma}|z} &= \frac{1}{(1+|z|^2)^2} \Big( \delta_{\mbf{k}\mbf{k}'} \delta_{im'} \delta_{\sigma' s} \delta_{jm} \delta_{\sigma s'} - |z|^2 \delta_{\mbf{k},-\mbf{k}'} \delta_{m'j} \delta_{im} s_z^{s'} \delta_{\sigma',-s'} s_z^{s} \delta_{\sigma,-s} \\
& \qquad\qquad\qquad + |z|^2 \delta_{mm'}\delta_{\sigma \sigma'} \delta_{ij} \delta_{ss'}
\Big) \braket{z|z} \ .
\eea
Plugging in, we find
\bea
\braket{z|\gamma_{\mbf{k},m',\sigma'} \tilde{H}' \gamma^\dag_{\mbf{k},m,\sigma}|z} &= (1+|z|^2)^{N_f \mathcal{N}-2} \Big( \tilde{E}'_{m'\sigma',m\sigma}(\mbf{k}) -|z|^2 s_z^{\sigma'} s_z^{\sigma} \tilde{E}'_{m(-\sigma),m'(-\sigma')}(-\mbf{k}) +  |z|^2 \delta_{mm'} \delta_{\sigma \sigma'} 2N_f \mathcal{N} \tilde{E}' \ .
\Big)
\eea
Now expanding in powers of $|z|^2$ and including the normalization \Eq{eq:1pnorm}, we find
\bea
\label{eq:degepertumatr}
\binom{N_f\mathcal{N}}{n} \frac{\braket{n|\gamma_{\mbf{k},m',\sigma'} \tilde{H}' \gamma^\dag_{\mbf{k},m,\sigma}|n} }{1-\nu} &= \frac{1}{1-\nu} \lp \binom{N_f\mathcal{N}-2}{n} \tilde{E}'_{m'\sigma',m\sigma}(\mbf{k}) - \binom{N_f\mathcal{N}-2}{n-1} s_z^{\sigma'} s_z^{\sigma} \tilde{E}'_{m(-\sigma),m'(-\sigma')}(-\mbf{k}) \right. \\
&\left.\qquad + \binom{N_f\mathcal{N}-2}{n-1} \delta_{mm'} \delta_{\sigma \sigma'}2N_f \mathcal{N} \tilde{E}' \rp \ . \\
\eea
Simplifying the binomal coefficients gives
\bea
\frac{\braket{n|\gamma_{\mbf{k},m',\sigma'} \tilde{H}' \gamma^\dag_{\mbf{k},m,\sigma}|n} }{1-\nu} &= \frac{1}{1-\nu} \lp \frac{(N_f\mathcal{N}-n-1)(N_f\mathcal{N}-n)}{N_f \mathcal{N}(N_f\mathcal{N}-1)} \tilde{E}'_{m'\sigma',m\sigma}(\mbf{k}) - \frac{n(N_f\mathcal{N}-n)}{N_f \mathcal{N}(N_f\mathcal{N}-1)} s_z^{\sigma'} s_z^{\sigma} \tilde{E}'_{m(-\sigma),m'(-\sigma')}(-\mbf{k}) \right. \\
&\left.\qquad + \frac{n(N_f\mathcal{N}-n)}{N_f \mathcal{N}(N_f\mathcal{N}-1)}  \delta_{mm'} \delta_{\sigma \sigma'}2N_f \mathcal{N} \tilde{E}' \rp \\
&=  \lp \frac{N_f\mathcal{N}-n-1}{N_f\mathcal{N}-1} \tilde{E}'_{m'\sigma',m\sigma}(\mbf{k}) - \frac{n}{N_f\mathcal{N}-1} s_z^{\sigma'} s_z^{\sigma} \tilde{E}'_{m(-\sigma),m'(-\sigma')}(-\mbf{k})  + \frac{n}{N_f\mathcal{N}-1}  \delta_{mm'} \delta_{\sigma \sigma'}2N_f \mathcal{N} \tilde{E}' \rp \\
%&=  (1-\nu) \tilde{E}'_{m'\sigma',m\sigma}(\mbf{k}) - \nu  s_z^{\sigma'} s_z^{\sigma} \tilde{E}'_{m(-\sigma),m'(-\sigma')}(-\mbf{k}) + \nu  \delta_{mm'} \delta_{\sigma \sigma'}N_f V \tilde{E}' \\
\eea
To write the effective Hamiltonian more transparently, we define the matrix notation
\bea
\null [\tilde{E}'(\mbf{k})]_{m'\sigma', m\sigma} &= \tilde{E}'_{m'\sigma',m\sigma}(\mbf{k}) \\
\null [\sigma_y \tilde{E}'^T(-\mbf{k}) \sigma_y]_{m'\sigma', m\sigma} &= s_z^{\sigma'} s_z^{\sigma}\tilde{E}'_{m(-\sigma),m'(-\sigma')}(-\mbf{k}) \\
\eea
where $σ_y$ is a Pauli matrix acting only on the spin indices. Because $H̃$ is Hermitian, $E^T = E^*$. Thus we can write the matrix \Eq{eq:degepertumatr} obtained from degenerate perturbation theory as
\bea
\frac{\braket{n|\gamma_{\mbf{k},m',\sigma'} \tilde{H}' \gamma^\dag_{\mbf{k},m,\sigma}|n} }{1-\nu} &= \left[  \frac{N_f\mathcal{N}-n-1}{N_f\mathcal{N}-1}  \tilde{E}'(\mbf{k}) -  \frac{n}{N_f\mathcal{N}-1}  \sigma_y \tilde{E}'^*(-\mbf{k}) \sigma_y + \frac{n}{N_f\mathcal{N}-1} \mathbb{1} \, 2N_f \mathcal{N} \tilde{E}'\right]_{m'\sigma',m\sigma} \ . \\
\eea
The eigen-energies of this matrix split the electron excitations which were degenerate at energy $|U|/2$, forming a band structure for the charge +1 excitations. Let $_min[M(k)]$ denote the smallest eigenvalue of $M(k)$ for all $k$. Then the lowest energy state of the one-particle excitations is
\bea
\label{eq:Richarson2}
E(N+1) = \la_{min}\left[\frac{\eps |U|}{2} \mathbb{1} + \frac{N_f\mathcal{N}-n-1}{N_f\mathcal{N}-1}  \tilde{E}'(\mbf{k}) -  \frac{n}{N_f\mathcal{N}-1}  \sigma_y \tilde{E}'^*(-\mbf{k}) \sigma_y + \frac{n}{N_f\mathcal{N}-1} \mathbb{1} \, 2N_f \mathcal{N} \tilde{E}'\right]
\eea
corresponding to the lowest energy single-particle excitation in perturbation theory.  With \Eqs{eq:Richarson1}{eq:Richarson2}, the Richardson binding energy can be written as
\bea
\label{eq:Delta}
E_\Delta(N) &= E(N+2) + E(N) - 2 E(N+1) \\
&= \frac{n +(n+1)}{N_f\mathcal{N}} (2N_f \mathcal{N} \tilde{E}') \mathbb{1}  - 2 \la_{min} \left[ \frac{\eps|U|}{2}\mathbb{1} + \frac{N_f\mathcal{N}-n-1}{N_f\mathcal{N}-1}  \tilde{E}'(\mbf{k}) -  \frac{n}{N_f\mathcal{N}-1}  \sigma_y \tilde{E}'^*(-\mbf{k}) \sigma_y + \frac{n}{N_f\mathcal{N}-1} \mathbb{1}2 N_f \mathcal{N} \tilde{E}'  \right] \\
&= \lp \frac{2n+1}{N_f\mathcal{N}} -\frac{2n}{N_f\mathcal{N}-1} \rp (2N_f \mathcal{N} \tilde{E}') - \eps|U| -2 \la_{min} \left[ \frac{N_f\mathcal{N}-n-1}{N_f\mathcal{N}-1}  \tilde{E}'(\mbf{k}) -  \frac{n}{N_f\mathcal{N}-1}  \sigma_y \tilde{E}'^*(-\mbf{k}) \sigma_y \right] \\
&= - \eps|U| + 2(1-2\nu) \tilde{E}'  -2 \la_{min} \left[ (1-\nu)  \tilde{E}'(\mbf{k}) -  \nu \sigma_y \tilde{E}'^*(-\mbf{k}) \sigma_y \right] + O(1/\mathcal{N})\\
\eea
where in the last line we took the thermodynamic limit with $ν= n/N_f𝒩$. In going from the second to third lines we have pulled out terms proportional to the identity out of $λ_min$: $λ_min[aI + B] = a + λ_min[B]$. To determine at what filling $ν$ the superconductor is strongest (i.e. has the largest gap), we need to determine when $-E_Δ(N)$ is largest. In general, this would require diagonalizing $(1-ν)  Ẽ'(k) - νσ_y Ẽ'^*(-k) σ_y$ as a function of $ν$. This is easily done for a specific model since $Ẽ'(k)$ is an effective single-particle Hamiltonian. Generically, \Eq{eq:Delta} can reveal material-specific effects on superconductivity due spin-orbit coupling, applied magnetic field, and other small perturbations.

However, we are able to make a general statement when time-reversal symmetry is preserved by $H̃'$, ensuring that $σ_y Ẽ'^*(-k) σ_y = Ẽ'(k)$. We first take $Ẽ'=0$ for convenience; we will restore nonzero $Ẽ'$ later. Then noting that $_min[Ẽ'(k)] ≤0, _max[Ẽ'(k)] ≥0$, we find
\bea
E_\Delta(N) &= - \eps|U| -2 \la_{min} \left[ (1-2\nu)  \tilde{E}'(\mbf{k}) \right]  \\
&= - \eps|U| -2 \begin{cases}
(1-2\nu) \la_{min} [\tilde{E}'(\mbf{k})], & \nu < 1/2\\
(1-2\nu)  \la_{max} [\tilde{E}'(\mbf{k})], & \nu > 1/2\\
\end{cases} \\
&= - \eps|U| + 2|1-2\nu| \begin{cases}
|\la_{min} [\tilde{E}'(\mbf{k})] |, & \nu < 1/2\\
\la_{max} [\tilde{E}'(\mbf{k})], & \nu > 1/2\\
\end{cases} \\
\eea
so because the second term is non-negative with a minimum at $ν= 1/2$, we see that $E_Δ(N)$ is minimized (the binding energy is largest and negative) at $ν= 1/2$. This holds for arbitrary time-reversal symmetric $Ẽ'(k)$ and is independent of the position of the peak in the density of states, which is very different from standard BCS theory.

Restoring nonzero $E'$ gives
\bea
E_\Delta(N) &= - \eps|U| +2(1-2\nu) {\tilde E}' -2 \begin{cases}
(1-2\nu) \la_{min} [\tilde{E}'(\mbf{k})], & \nu < 1/2\\
(1-2\nu)  \la_{max} [\tilde{E}'(\mbf{k})], & \nu > 1/2\\
\end{cases} \\
&= - \eps|U| -2 \begin{cases}
(1-2\nu) (\la_{min} [\tilde{E}'(\mbf{k}) -{\tilde E}' I]), & \nu < 1/2\\
(1-2\nu)  (\la_{max} [\tilde{E}'(\mbf{k}) -{\tilde E}' I]), & \nu > 1/2\\
\end{cases}.
\eea  Because we have subtracted the average $Ẽ'$ from the matrix $Ẽ'()$, $_min [Ẽ'(k) -Ẽ' I] ≤0$ and $_max [Ẽ'(k) -Ẽ' I] ≥0$.  The energy $E_Δ(N)$ is still minimized at $ν= 1/2$.

\subsection{Charge $+2$ Excitations}
\label{app:charge2}

We now study the charge $+2$ excitations. To do so, we will compute the scattering matrix $R^σσ'_k'm'n',kmn(p)$ on the complete basis of charge $+2$ excitations at many-body momentum $p$ defined by
\bea
\label{eq:Hcharge2}
\null [H, \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}\gamma^\dag_{-\mbf{k},n,\sigma'}] \ket{GS} &=  \sum_{\mbf{k}'} \gamma^\dag_{\mbf{p}+\mbf{k}',m',\sigma}\gamma^\dag_{-\mbf{k}',n',\sigma'}  \ket{GS} [R^{\sigma \sigma'}(\mbf{p})]_{\mbf{k}'m'n',\mbf{k}mn} \ .
\eea
There are three flavors of Cooper pair. The $σ= σ' = ±1/2$ cases lead to spin $±1$ Cooper pairs, while the $σ= - σ'$ case leads to a spin 0 Cooper pair. We will find that only the spin 0 Cooper pair has low-energy eigenstates below the two-particle continuum. For instance, we know that there must be a zero-energy excitation $η^†$ at $p=0$ since $η$ is a symmetry.

This calculation is made simpler by using the results from the charge $+1$ electron excitation matrix. Using \Eqs{eq:formfactorS}{eq:SSgamma}, the single electron commutator is
\bea
\null [H, \gamma^\dag_{\mbf{k}, n, \sigma}] &= \frac{\eps |U|}{2} \gamma^\dag_{\mbf{k}, n, \sigma} + \frac{|U|}{2} \sum_{\mbf{q} \al m} 2  \frac{1}{\sqrt{\mathcal{N}}} s^z_\sigma \gamma^\dag_{\mbf{k}+\mbf{q},m,\sigma} M^{mn}_{\sigma,\al}(\mbf{k},\mbf{q}) \bar{S}^z_{-\mbf{q},\al}
\eea
and the second term annihilates the groundstate. Now it is straightforward to compute
\bea
\null [H, \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}\gamma^\dag_{-\mbf{k},n,\sigma'}] \ket{GS} &=  \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}  [H, \gamma^\dag_{-\mbf{k},n,\sigma'}] \ket{GS}  + [H, \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}] \gamma^\dag_{-\mbf{k},n,\sigma'} \ket{GS}  \\
&= \eps |U| \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}  \gamma^\dag_{-\mbf{k},n,\sigma'} \ket{GS}  + |U| \frac{1}{\sqrt{\mathcal{N}}}  \sum_{\mbf{q} \al m'} s^z_\sigma \gamma^\dag_{\mbf{p}+\mbf{q}+\mbf{k},\sigma,m'} M^{m'm}_{\sigma,\al}(\mbf{k}+\mbf{p},\mbf{q}) [\bar{S}^z_{-\mbf{q},\al}, \gamma^\dag_{-\mbf{k},n,\sigma'}] \ket{GS} \\
&= \eps |U| \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}  \gamma^\dag_{-\mbf{k},n,\sigma'} \ket{GS}  \\
&\qquad + |U| \frac{1}{\mathcal{N}}  s^z_\sigma  s^z_{\sigma'} \sum_{\mbf{q} \al m'} \gamma^\dag_{\mbf{p}+\mbf{q}+\mbf{k},\sigma,m'} \gamma^\dag_{-\mbf{k}-\mbf{q},\sigma',n'} M^{m'm}_{\sigma,\al}(\mbf{k}+\mbf{p},\mbf{q})  M^{n'n}_{\sigma',\al}(-\mbf{k},-\mbf{q})  \ket{GS} \\
\eea
where in the second line we introduced the commutator $[S̅^z_-q,, γ^†_-k,n,σ']$ since $S̅^z_-q,$ annihilates the vacuum. Resumming gives
\bea
\null  [H, \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}\gamma^\dag_{-\mbf{k},n,\sigma'}] \ket{GS}   &= \eps |U| \gamma^\dag_{\mbf{p}+\mbf{k},m,\sigma}  \gamma^\dag_{-\mbf{k},n,\sigma'} \ket{GS}  \\
&\qquad + |U| s^z_\sigma  s^z_{\sigma'}  \frac{1}{\mathcal{N}} \sum_{\mbf{k}' \al} \gamma^\dag_{\mbf{p}+\mbf{k}',\sigma,m'} \gamma^\dag_{-\mbf{k}',\sigma',n'}  \ket{GS}  M^{m'm}_{\sigma,\al}(\mbf{k}+\mbf{p},\mbf{k}'-\mbf{k})  M^{n'n}_{\sigma',\al}(-\mbf{k},\mbf{k}-\mbf{k}') \\
\eea
and thus we find
\bea
\label{eq:Rscater}
\null [R^{\sigma \sigma'}(\mbf{p})]_{\mbf{k}'m'n',\mbf{k}mn} &= \eps |U| \delta_{\mbf{k}\mbf{k}'} \delta_{mm'} \delta_{nn'} + s^z_\sigma  s^z_{\sigma'} |U| \frac{1}{\mathcal{N}} \sum_\al M^{m'm}_{\sigma,\al}(\mbf{k}+\mbf{p},\mbf{k}'-\mbf{k})  M^{n'n}_{\sigma',\al}(-\mbf{k},\mbf{k}-\mbf{k}') \ .
\eea
The second term describes the nontrivial scattering processes. We observe that
\bea
\label{eq:nontrivtermMM}
\sum_\al M^{m'm}_{\sigma,\al}(\mbf{k}+\mbf{p},\mbf{k}'-\mbf{k})  M^{n'n}_{\sigma',\al}(-\mbf{k},\mbf{k}-\mbf{k}') &= \sum_\al U^{\sigma*}_{m' \al}(\mbf{p}+\mbf{k}') U^{\sigma}_{m \al}(\mbf{p}+\mbf{k}) \, U^{\sigma'*}_{n' \al}(-\mbf{k}') U^{\sigma'}_{n \al}(-\mbf{k}) \\
&= \sum_\al \mathcal{U}^{\sigma \sigma'}_{\mbf{k}'m'n',\al}(\mbf{p}) \mathcal{U}^{\sigma \sigma' *}_{\mbf{k}mn,\al}(\mbf{p}), \quad \mathcal{U}^{\sigma \sigma'}_{\mbf{k}mn,\al}(\mbf{p}) = U^{\sigma *}_{m \al}(\mbf{p}+\mbf{k}) U^{\sigma' *}_{n \al}(-\mbf{k})  \ .
\eea
Let us consider the case where $σ= σ'$. Note that if $σ= σ'$, $𝒰^σσ'_k'm'n',(p)$ is symmetric under $k' →-(p+k'), m' ↔n'$. But the electron operators $γ^†_p+k',m',σγ^†_-k',n',σ$ which contract against $𝒰^σσ'_k'm'n',(p)$ in \Eq{eq:Hcharge2} are anti-symmetric under $k' →-(p+k'), m' ↔n'$. Thus many-body states which are \emph{not} in the kernel of \Eq{eq:nontrivtermMM} all vanish by fermion anti-symmetry when $σ= σ'$. Hence all the spin $±1$ Cooper pairs all have energy $|U|$, which is twice the single-electron gap. This indicates they are effectively unpaired. Going forward, we take $σ= ,̆ σ' = $̣.

We can think of $\mathcal{U} \equiv \mathcal{U}^{\u \d}$ as a large rectangular matrix of dimension $\mathcal{N} N_f^2 \times N_{orb}$. (In the bipartite case, the matrix is $\mathcal{N} N_f^2 \times N_{L}$ because $U^\sigma_{m\al}(\mbf{k}) = 0$ for $\al \notin L$.) In matrix notation, we have

R^̣̆(p)    = |U| 1 - |U| 1/𝒩 𝒰(p) 𝒰^†(p)  .

Note that $ \mathcal{U}(\mbf{p}) \mathcal{U}^\dag(\mbf{p})$ is positive semi-definite, so the spin 0 Cooper pair spectrum is below $\eps |U|$, indicating attractive pairing. Because $\mathcal{U}$ is rectangular, it has low rank and $\frac{1}{\mathcal{N}} \mathcal{U}\mathcal{U}^\dag$ has many zero eigenvalues. Its nonzero eigenvalues are the same as the eigenvalues of $h(\mbf{p}) = \frac{1}{\mathcal{N}}\mathcal{U}^\dag(\mbf{p}) \mathcal{U}(\mbf{p})$ which is a matrix of dimension $N_{orb} \times N_{orb}$. To see this, define the eigenvector $u$ with energy $\eps(\mbf{p})$ by $h(\mbf{p}) u = \eps(\mbf{p}) u$. Then $\mathcal{U} u$ is an eigenvector of the scattering matrix: $\frac{1}{\mathcal{N}} \mathcal{U}\mathcal{U}^\dag (\mathcal{U} u) =  \mathcal{U} h(\mbf{p}) u = \eps(\mbf{p}) \mathcal{U} u$. This gives $\dim h(\mbf{p}) = N_{orb}$ nontrivial eigenvalues, which are nonnegative because $ \mathcal{U}\mathcal{U}^\dag $ is positive-semi-definite. Note that the rank of $\mathcal{U}\mathcal{U}^\dag$ is at most $N_{orb}$, so all other eigenvalues are zero. Hence we are led to study the spectrum of

h_(p)    = 1/𝒩 [𝒰^†(p) 𝒰(p)]_

= 1/𝒩 ∑_kmn U^_m (p+k) U^_n (-k)  U^*̆_m (p+k) U^*̣_n (-k)

= 1/𝒩 ∑_k P^_(p+k) P^_(-k)

=1/𝒩 ∑_k P_(p+k) P_(k)

where we used $P^\d(\mbf{k}) = P^\u(-\mbf{k})^*$, and set $P^\u(\mbf{k}) = P(\mbf{k})$.

A number of extremely useful statements will be proven about this effective single-particle Hamiltonian in the next section. For now, we simply remark that the eigenvectors $u^\be_\mu(\mbf{p})$ and eigenvalues $\eps_\mu(\mbf{p})$ are defined as
∑_h_(p) u^_μ(p) = _μ(p) u^_μ(p)  .

We normalize the eigenvectors such that $\sum_\be u^\be_\mu(\mbf{p}) u^{\be *}_\nu(\mbf{p}) = \delta_{\mu \nu}$ where $\mu,\nu \in 0,\dots, N_{orb}-1$ label the boson bands. It is convenient to decreasingly order the eigenvalues of $h(\mbf{p})$ such that $\mu = 0$ corresponds to the largest eigenvalue. The full eigenstates of the scattering matrix eq:Rscater not in the kernel of $\mathcal{U}\mathcal{U}^\dag$ are $\sum_\al \mathcal{U}_{\mbf{k}mn,\al}(\mbf{p}) u^\al_\mu(\mbf{p})$ with eigenvalue
R^̣̆(p) 𝒰 u_μ= |U| (- _μ(p) )  𝒰 u_μso the many-body energies are $|U| (\eps - \eps_\mu(\mbf{p}) )$; $\mu=0$ corresponds to the lowest energy band of many-body states.

The eigenvectors $\mathcal{U} u_\mu$ have norm squared $u^\dag_\mu(\mbf{p}) \mathcal{U}^\dag(\mbf{p}) \mathcal{U}(\mbf{p}) u_\mu(\mbf{p}) = \mathcal{N} u^\dag_\mu(\mbf{p}) h(\mbf{p}) u_\mu(\mbf{p}) = \mathcal{N} \eps_\mu(\mbf{p})$. Thus we define the many-body operators

η^†_p,μ    = 1/√(𝒩 _μ(p)) ∑_kmn γ^†_p+k,m,γ^†_-k,n, 𝒰_kmn,(p) u^_μ(p)

= 1/√(𝒩 _μ(p)) ∑_kmn γ^†_p+k,m,γ^†_-k,n, U^*̆_m (p+k) U^*̣_n (-k)  u^_μ(p)

= 1/√(𝒩) ∑_k u^_μ(p)/√(_μ(p)) c̅^†_p+k,,c̅^†_-k,,

= 1/√(𝒩) ∑_R u^_μ(p)/√(_μ(p))  e^-i p ·(R+r_) c̅^†_R,,  c̅^†_R,,

which obey
[H, η^†_p,μ] |GS⟩    = |U| (- _μ(p) ) η^†_p,μ |GS⟩  .

Comparing with eq:etadef, we see that the Cooper pair operators only show onsite pairing in $\mbf{R},\al$ in the $\bar c$ degrees of freedom, indicating closely bound states. We will show that, up to a normalization, $\eta^\dag_{\mbf{p}=0,\mu=0} \propto \eta^\dag$ and that $\eps_{\mu=0}(\mbf{p}=0) = \eps$. Thus the Cooper pair operators $\eta^\dag_{\mbf{p},\mu}$ are deformations of the commuting $\eta^\dag$ symmetry with the same $\bar{c}^\dag_{\mbf{R},\al,\u}  \bar{c}^\dag_{\mbf{R},\al,\d}$ orbital structure. They are gapless at $\mbf{p}=0$ as a consequence of Goldstone's theorem. The spectrum of the $\eta^\dag_{\mbf{p},\mu}$ operators contains the low lying Cooper pair (which condenses) as well as higher energy Cooper pair excitations (also known as Leggett modes).

We calculate the normalization of the operators in eq:etap.  First, building the excitations above the vacuum gives norm
⟨0|η_q,ν  η^†_p,μ |0|   ⟩= 1/𝒩 ∑_kk', u^_μ(p) u^*_ν(q)/√(_μ(p) _ν(q)) ⟨0|c̅_-k', c̅_q+k',, c̅^†_p+k,,c̅^†_-k,, |0|
= 1/𝒩 ∑_kk', u^_μ(p) u^*_ν(q)/√(_μ(p) _ν(q)) δ_kk' δ_pq P_^(̆p+k) P_^(̣-k)

= δ_pq/√(_μ(p) _ν(q))  ∑_  u^*_ν(q)   1/𝒩 ∑_kP_^(̆p+k) P_^(̣-k) u^_μ(p)

= δ_pq/√(_μ(p) _ν(q))  ∑_  u^*_ν(q) h_(p) u^_μ(p)

= δ_pq/√(_μ(p) _ν(q)) _μ(p) δ_μν

= δ_pq δ_μν  .

To compute the normalization of the states $\eta^\dag_{\mbf{p},\mu}\ket{n}$, we use Wick's theorem. First we need the correlator

⟨z|γ_-k',n', γ_q+k',m', γ^†_p+k,m,γ^†_-k,n, |z⟩    = 1/(1+|z|^2)^2 |z|^2 δ_q,0 δ_n'm' δ_p,0 δ_nm - 0 + δ_kk' δ_nn' δ_p,qδ_mm' ⟨z|z|
=  δ_p,q/(1+|z|^2)^2 |z|^2 δ_p,0 δ_n'm'  δ_nm  + δ_kk' δ_nn'δ_mm' ⟨z|z| ⟩.

We will also need the identity (see eq:nontrivtermMM)

1/𝒩∑_km 𝒰^σ,-σ_kmm,(0)    =1/𝒩 ∑_km U^σ*_m (k) U^-σ*_m (-k) = 1/𝒩 ∑_km U^σ*_m (k) U^σ_m (k) = 1/𝒩∑_k P^σ_(k) = using the uniform pairing condition. Thus we obtain

⟨z|η_q,ν η^†_p,μ|z⟩    = δ_p,q (1+|z|^2)^N_f 𝒩-2 /𝒩 √(_μ(p)_ν(q)) ∑_kk',  ∑_mn,m'n' 𝒰^*_k'm'n',(q) u^*_ν(q) |z|^2 δ_n'm' δ_p,0 δ_nm + δ_kk' δ_nn' δ_mm' 𝒰_kmn,(p) u^_μ(p)

= δ_p,q (1+|z|^2)^N_f 𝒩-2 /𝒩 √(_μ(p)_ν(p)) ∑_ u^*_ν(p) |z|^2 δ_p,0 ∑_k'm'𝒰^*_k'm'm',(0) ∑_km𝒰_kmm,(0)  + [𝒰^†(p)𝒰(p)]_ u^_μ(p)

= δ_p,q (1+|z|^2)^N_f 𝒩-2 /√(_μ(p)_ν(p)) ∑_ u^*_ν(p) |z|^2 δ_p,0 𝒩 ^2 +h_ (p)u^_μ(p)

= δ_p,q (1+|z|^2)^N_f 𝒩-2 /√(_μ(p)_ν(p)) |z|^2 δ_p,0 𝒩 ^2  ∑_ u^*_ν(0) u^_μ(0) +δ_μν _μ(p) .

In app:Cooperpairspectrum, we will prove that $\eps_0(\mbf{p}=0) = \eps$ and $ \sum_{\al} u^{\al}_\mu(0) = \delta_{\mu,0} \sqrt{N_{L}}$ where $N_L$ is the number of orbitals where $P_{\al \be}(\mbf{k})$ is nonzero. With this simple result, we obtain
⟨z|η_q,ν η^†_p,μ|z⟩    = δ_p,q (1+|z|^2)^N_f 𝒩-2 /√(_μ(p)_ν(p)) |z|^2 δ_p,0 𝒩 N_L^2 δ_μ,ν δ_μ,0  +δ_μν _μ(p)

= δ_p,q δ_μ,ν (1+|z|^2)^N_f 𝒩-2 |z|^2 δ_p,0 𝒩 N_Lδ_μ,0 +1

= δ_p,q δ_μ,ν (1+|z|^2)^N_f 𝒩-2 |z|^2 δ_p,0 δ_μ,0 N_f 𝒩  +1

making use of $\eps = N_f/N_{L}$. We see that there is an anomalous term when $\mbf{p}=0,\mu=0$. This is because $\eta^\dag_{\mbf{p},0} \propto \eta^\dag$, which produces a zero-energy eta-pairing state. These states are the ground states we have already computed the normalization for in eq:normedetastates. Expanding in powers of $z$ and picking $\mbf{p}\neq 0, \mu \neq 0$, we obtain
∑_n |z|^2n N_f𝒩n ⟨n|η_q,ν η^†_p,μ|n⟩    = ∑_n |z|^2n δ_p,q δ_μ,ν  N_f𝒩-2n ,     p≠0, μ≠0

⟨n|η_q,ν η^†_p,μ|n⟩    = (1-ν)^2 δ_p,q δ_μ,ν

and in the last line we took the thermodynamic limit $\mathcal{N}\to \infty$.

§.§ Examples of the Cooper Pair Spectrum

To illustrate the Cooper pair spectrum, we consider three examples: an inversion-protected SSH chain in 1D, an obstructed atomic limit in 2D, and an $S$-matrix construction of fragile bands in 2D. We use the same notation $\tilde{h}(\mbf{k})$ for the single-particle electron Hamiltonian, $P(\mbf{k})$ for the flat band projector, and $h(\mbf{p})$ for the bosonic pairing Hamiltonian for all the models used in this section. We recall that in general, the projection into the flat bands is valid when $W \ll |U| \ll t$ where $W$ is the bandwidth of the flat band ($W=0$ if the band is perfectly flat), $|U|$ is the interaction strength, and $t$ is the single-particle gap above the flat band.

§.§.§ SSH Chain

The 1D Hamiltonian we consider was studied in 2020Sci...367..794S as a minimal model of a flat band obstructed atomic insulator. We place $s$ and $p$ orbitals at the 1a position (the center of the unit cell) in the 1D space group $pm$. The electron Hamiltonian is
h̃(k)    = sink   σ_2 + cosk   σ_3
with exactly flat bands at energies $\pm 1$. The eigenstate of the lower band is
U(k) = 1/2 1 - e^-ik
1+ e^- i k

and one can calculate the Berry phase $U^\dag i \del_k U = \frac{1}{2}$, indicating an obstruction to the trivial atomic phase where the Wannier functions are onsite. Indeed, $h(k)$ obeys inversion symmetry with $D[\mathcal{I}] h(k)D^\dag[\mathcal{I}] = h(-k)$ where $D[\mathcal{I}] = \sigma_3$, and real space invariants and/or symmetry eigenvalues diagnose $h(k)$ as being in an obstructed atomic limit [Song et al., 2020].

Projecting into the single flat band, we find that the uniform pairing condition is obeyed:

∫dk/2π P_(k) = N_f/N_orb = = 1/2 ,     P(k) = U(k)U^†(k) = (

sin^2k/2     i/2 sink

-i/2 sink     cos^2 k/2

)  .
Note that the uniform pairing condition is not enforced by inversion because the $s$ and $p$ orbitals are different irreps of the 1a Wyckoff position and hence are not interrelated by crystalline symmetry. In other words, adding $\mathcal{I}$-preserving hoppings could destroy the uniform pairing condition.  However, there may be non-crystallographic symmetries that enforce UPC.  In this model, two anti-commuting symmetries: ${\cal C} = \sigma_1$ and ${\cal Z} = \kk \rightarrow \kk + \pi$, combine into a commuting unitary symmetry ${\cal CZ}$ that enforces UPC:

\begin{align}
{\cal CZ} P(\kk) {\cal CZ} = \sigma_1 \begin{bmatrix}
\cos ^2\frac{k+\pi}{2} & -\frac{i}{2} \sin (k +\pi) \\
\frac{i}{2} \sin (k +\pi) & \sin ^2 \frac{k+\pi}{2} \\
\end{bmatrix} \sigma_1 = P(\kk) \\
P_{\alpha\alpha}(\kk) = \sum_\beta [\sigma_{1}]_{\alpha\beta} P_{\beta\beta}(\kk+\pi).
\label{}
\end{align}

The form of the interaction is still density-density, with the spin up $s$-orbital density multiplying the spin down $s$-orbital density and the same for the $p$-orbitals:
H_int    = -|U| ∑_R,= s,p n̅_R,, n̅_R,, .

Using the formula for the projector in eq:SSHproj, we find the pairing Hamiltonian
h_(p) = ∫dk/2π P_(k+p) P_(k)  = 1/8 (

2 + cosp     cosp

cosp     2 + cosp

which has eigenvalues $\eps_\mu(p) = \frac{1}{4}(1+\cos p), \frac{1}{4}$. The many-body energies are $E_\mu(p)/|U| = \eps - \eps_\mu(p) = \frac{1}{4}(1- \cos p), \frac{1}{4}$. The unpaired continuum of states appears at $E/|U| = \eps = \frac{1}{2}$. The spectrum is depicted in fig:SSHOAL(a).

As a final comment, we note that in this symmetry group $pm$ , one may imagine deforming the $s$ and $p$ orbitals at 1a into symmetric and antisymmetric combinations of orbitals at non-maximal Wyckoff positions of $s$ orbitals at 2b. This deformation involves a unitary transformation that combines $s$ and $p$ orbitals at position ${\bf R}$ into symmetric and antisymmetric combinations.  While the orbitals now obey the uniform pairing condition, the density-density interactions are not invariant under this transformation.  However, the following interaction
H    = -|U|2 ∑_R,σ= ,̆ [n̅_R,s,σ + n̅_R,p,σ]^2
is invariant, though we are currently unable to treat this interaction with the eta-pairing operator.

§.§.§ $p3$ Compact Obstructed atomic Insulator

A second compact OAI model with the uniform pairing condition and flat bands was first presented in 2022PhRvL.128h7002H in the wallpaper group $p3$. The Hamiltonian is built from $s$ and $p_x, p_y$ orbitals on the 1a position (the unit cell center) which form the $C_3$ irreps $A, {}^1E, {}^2E$. As constructed in 2022PhRvL.128h7002H, the single valence band of the Hamiltonian forms a compact obstructed atomic limit, much like the obstructed phase of the SSH chain in app:SSHchain, and is indicated by real space invariants. Thus $\eps = \frac{N_f}{N_{orb}} = \frac{1}{3}$. The valence band eigenstate is
U(k) = 1/3 1
1 +  1/3 1
e^- 2πi/3
e^ 2πi/3 e^i (k_1 + k_2) +  1/3 1
e^-2πi/3 e^i k_2  .
We can calculate the pairing Hamiltonian analytically from $P(\mbf{k}) = U(\mbf{k}) U^\dag(\mbf{k})$:
h_(p)    = ∫dk_1 dk_2/(2π)^2 P_(p+k) P_(k) = 1/9 δ_ + 2/81 (cos(p_1 + p_2) + cosp_1 + cosp_2)

h(p)    = 2/81 (cos(p_1 + p_2) + cosp_1 + cosp_2) 1     1     1
1     1     1
1     1     1
+ 1/91

which has eigenvalues $\eps_\mu(\mbf{p}) = \frac{1}{27}(3 + 2 \cos (p_1+p_2) + 2 \cos p_1 + 2 \cos p_2), \frac{1}{9},  \frac{1}{9}$. The many-body spectrum $E_\mu(\mbf{p})/|U| = \frac{1}{3} - \eps_\mu(\mbf{p})$ is shown in fig:SSHOAL(b) along with the un-paired continuum at $E/|U| = \eps = \frac{1}{3}$.

The many-body Cooper pair/density excitation spectrum $E_\mu(\mbf{p})= |U|(\eps - \eps_\mu(\mbf{p}))$ is shown for $(a)$ the 1D SSH-like model in app:SSHchain and $(b)$ the 2D compact obstructed atomic limit in app:p3OAL. The blue band is the quadratic gapless mode guaranteed by the $su(2)$ eta-pairing symmetry. The red band is perfectly flat for both models and is singly degenerate at $\frac{E}{|U|} = \frac{1}{4}$ in $(a)$ and doubly degenerate at $\frac{E}{|U|} = \frac{2}{9}$ in $(b)$. The dashed black line shows the macroscopically degenerate continuum of unpaired states at $E = \eps |U|$. The flatness of the Leggett modes (in red) is due to the simple, low-harmonic structure of the eigenstates. See the Cooper pair spectrum in the Main Text for nontrivially dispersing Leggett modes.

§.§.§ $p6mm$ Fragile $S$-matrix Bands

Finally, we construct a perfectly flat fragile band in symmetry group $p6mm$ using the $S$-matrix method introduced in 2021NatPh..18..185C. To this end, we choose a 1D irrep (only one to enforce UPC) and place there an $s$-orbital, as we show in App. <ref> the Cooper pair transforms like an electronic $s$-orbital under the crystalline symmetries.  As the set of bands induced from an atomic position is trivial, to get a fragile set of bands one requires the induced band representation decomposes into two sets, one of them fragile.  As shown in tab:2Dcooperpairs, the Cooper pair band structure we desire can be induced by placing any 1D irrep (for instance an $s$ orbital) at the 3c (kagome) positions. The uniform pairing condition is guaranteed because the orbitals form an irrep of a single Wyckoff position. To obtain flat bands, we use the $S$-matrix construction define the two sublattices as $L = A_1^{3c}$ with $N_L = 3$ orbitals per unit cell and  $\tilde{L} = A_1^{1a}$ with $N_{\tilde{L}} = 1$ orbital per unit cell. The $A_1$ irrep denotes the $s$ orbital representation at the 1a and 3c positions (see the Bilbao Crystallographic Server <https://www.cryst.ehu.es/cgi-bin/cryst/programs/mbandrep.pl>). The momentum space irreps of the two sublattices are
ℬ_L    = A_1^3c ↑6mm = Γ_1 ⊕Γ_5 + K_1 ⊕K_3 + M_1 ⊕M_3 ⊕M_4

ℬ_L̃    = A_1^1a ↑6mm = Γ_1 + K_1 + M_1  .

In the bipartite $S$-matrix construction, hoppings are only allowed between the $L$ and $\tilde{L}$ sublattices resulting in $N_f = N_L - N_{\tilde{L}} = 3- 1  = 2$ flat bands. The momentum space irreps of the flat bands are[Cǎlugǎru et al., 2021]
ℬ_L ⊖ℬ_L̃ = Γ_5 + K_3  + M_3 ⊕M_4
which is fragile topological [Song et al., 2020]. To be concrete, we choose hoppings such that the single-particle electron Hamiltonian is

h̃(k) = 0     S^†(k)
S(k)    0_3×3 ,     S(k) = t cosk ·a_1/2
cosk ·(a_1+a_2)/2
cosk ·a_2/2 = t cosk_1/2

where $\mbf{a}_1 = (0,1), \mbf{a}_2 = C_3 \mbf{a}_1$ are the lattice vectors and $\mbf{a}_1/2, (\mbf{a}_1+\mbf{a}_2)/2, \mbf{a}_2/2$ are the sites of the 2c position. The projector into the flat bands of $\tilde{h}(\mbf{k})$ is given by

P(k)    =  0     0     0     0
0     2 + cosk_2 + cos(k_1 +k_2)     - 2 cosk_1/2 cosk_1+k_2/2     - cosk_1/2 cosk_2/2
0     - 2 cosk_1/2 cosk_1+k_2/2     2 + cosk_1 + cosk_2     - 2 cosk_2/2 cosk_1+k_2/2
0     - 2 cosk_1/2 cosk_2/2     - 2 cosk_2/2 cosk_1+k_2/2     2 + cosk_1 + cosk_1 + k_2
/ (3 + cosk_1 + cos(k_1 + k_2) + cosk_2 )
which is identically zero on the $1a$ position and obeys
∫dk_1 dk_2/(2π)^2 P_(k) = 0

which shows that the uniform pairing condition is satisfied on the 3c position, as guaranteed by app:UPC, with $\eps = \frac{N_f}{N_L} = \frac{2}{3}$.

We are unable to compute $h_{\al \be}(\mbf{p})$ analytically from the projectors in eq:projfrag, but it is very simple to compute the spectrum numerically.  The lattice model and Cooper pair band structure are depicted in Fig. <ref>.  The irreps of these bands are precisely those irreps induced by $s$ orbitals at $3c$:

ℬ_Cooper    = A_1^3c ↑6mm = Γ_1 ⊕Γ_5 + K_1 ⊕K_3 + M_1 ⊕M_3 ⊕M_4.
The Cooper pair band structure is decomposable, and as shown in Fig. <ref>(b), the blue band and red bands carry irreps
ℬ_1    = Γ_1 + K_1 + M_1

ℬ_2    = Γ_5 + K_3 + M_3 ⊕M_4  .

The first band carries the $\Gamma_1$ irrep that is the zero-energy Cooper pair, and its global topology is trivial because $\Gamma_1+K_1+M_1$ is the band representation of an $s$ orbital at the 1a position. The second set of bands is fragile topological as predicted from topological quantum chemistry, and cannot induced from local orbitals.

$(a)$ Bipartite hoppings $t$ of the single-particle electron Hamiltonian eq:tildeh6mm resulting in two-fold degenerate fragile flat bands. The three $s$ orbitals at the 3c position are blue, and the one $s$ orbital at the 1a position is marked in red. $(b)$ The many-body Cooper pair/Goldstone mode spectrum $E_\mu(\mbf{p})= |U|(\eps - \eps_\mu(\mbf{p}))$ is shown. The blue band is the quadratic gapless mode guaranteed by the $su(2)$ eta-paiting symmetry. The two red bands are fragile topological (see Eq. <ref>) and are gapped from the low-energy quadratic band. The the dashed black line shows the macroscopically degenerate continuum of unpaired states at $E = \eps |U|$.

§.§ Density Excitations

We now study the (spin) density modes which are in the form

[H, γ^†_p+k,m,σγ_k,n,σ'] |GS⟩    =  ∑_k' γ^†_p+k',m',σγ_k',n',σ'  |GS⟩ [ℛ^σσ'(p)]_k'm'n',kmn  .
The eigenstates obtained from eq:densitymodes are charge-neutral, so they contribute to the partition function at fixed particle-number.

To compute $\mathcal{R}^{\sigma \sigma'}(\mbf{p})$, we need the identities from app:charge1:
[S̅^z_q,,γ_k,n,σ]    = -1/√(𝒩) s^z_σ∑_m M^nm_σ,(k-q,q)  γ_k -q,m,σ

[H, γ_k, n, σ]    = |U|/2 γ_k, n, σ - |U|/2 ∑_q m 2  1/√(𝒩) s^z_σM^nm_σ,(k-q,q)  γ_k-q,m,σ S̅^z_-q,
and as usual, the second term in $[H, \gamma_{\mbf{k}, n, \sigma}]$ annihilates the groundstate. The scattering matrix $\mathcal{R}$ is obtained from a direct calculation:
[H, γ^†_p+k,m,σγ_k,n,σ'] |GS⟩    =  γ^†_p+k,m,σ  [H, γ_k,n,σ'] |GS⟩  + [H, γ^†_p+k,m,σ] γ_k,n,σ' |GS⟩

= |U| γ^†_p+k,m,σ  γ_k,n,σ' |GS⟩  + |U| 1/√(𝒩)  ∑_q m' s^z_σγ^†_p+q+k,m',σ M^m'm_σ,(k+p,q) [S̅^z_-q,, γ_k,n,σ'] |GS⟩

= |U| γ^†_p+k,m,σ  γ_k,n,σ' |GS⟩

- |U| 1/𝒩  s^z_σs^z_σ' ∑_q m' n' γ^†_p+q+k,σ,m' γ_k+q,σ',n' M^m'm_σ,(k+p,q)  M^nn'_σ',(k+q,-q)  |GS⟩

from which we identify the scattering matrix
[ℛ^σσ'(p)]_k'm'n',kmn    = |U| δ_kk' δ_mm' δ_nn' - s^z_σs^z_σ' |U| 1/𝒩 ∑_M^m'm_σ,(k+p,k'-k)  M^nn'_σ',(k',k-k')

= |U| δ_kk' δ_mm' δ_nn' - s^z_σs^z_σ' |U| 1/𝒩 ∑_U^σ*_m' (p+k') U^σ_m (p+k) U^σ' *_n (k) U^σ' _n' (k')

= |U| δ_kk' δ_mm' δ_nn' - s^z_σs^z_σ' |U| 1/𝒩 ∑_U^σ*_m' (p+k') U^σ' _n' (k')  U^σ_m (p+k) U^σ' *_n (k)  .

Let us recall the definition from eq:nontrivtermMM,
𝒰^σσ'_kmn,(p)    = U^σ*_m (p+k) U^σ' *_n (-k)  .
In the Cooper pair case (see discussion after eq:nontrivtermMM), the anti-commutation of the fermions, schematically $\gamma^\dag_{i\u} \gamma^\dag_{j\u} + \gamma^\dag_{j\u} \gamma^\dag_{i\u} =0$, prevented nontrivial excitations of nonzero spin. In the density excitation case, $\gamma^\dag_{i\u} \gamma_{i \d} + \gamma^\dag_{j\u} \gamma_{i \u} \neq 0$, so nontrivial spin $\pm1$ excitations are allowed. Enumerating the possible choices of $\sigma, \sigma'$ in eq:densitymodes gives the effective Hamiltonians
γ^†_p+k,m,γ_k,n,:     ℛ^̆̆(p)    = |U| 1 - |U| 1/𝒩 𝒰_̣̆(p) 𝒰_̣̆^†(p)

γ^†_p+k,m,γ_k,n,:     ℛ^̣̣(p)    = |U| 1 - |U| 1/𝒩 𝒰_̣̆(p) 𝒰_̣̆^†(p)

γ^†_p+k,m,γ_k,n,:     ℛ^̣̆(p)    = |U| 1 + |U| 1/𝒩 𝒰_̆̆(p) 𝒰_̆̆^†(p)

γ^†_p+k,m,γ_k,n,:     ℛ^̣̆(p)    = |U| 1 + |U| 1/𝒩 𝒰_̣̣(p) 𝒰_̣̣^†(p)  .

Note that the first two lines are spin 0 operators and the $\mathcal{U}\mathcal{U}^\dag$ term comes with an overall minus, while the last two lines are spin $\pm 1$ and the $\mathcal{U}\mathcal{U}^\dag$ term comes with an overall plus. Because $\mathcal{U}\mathcal{U}^\dag$ is positive semi-definite, we learn that the spin $0$ excitations are always at or below the two-electron gap $\eps |U|$, whereas the spin $\pm 1$ excitations are at or above the $\eps |U|$. This leads to a simple phenomenological prediction: the spin susceptibility is exponentially suppressed at temperatures below the spin gap $\eps|U|$.

By the same argument in eq:lowrank, the nonzero eigenvalues of $\mathcal{U}\mathcal{U}^\dag$ are determined by the effective single-particle Hamiltonians

ℛ^̆̆:   h_̆̆(p)    = h(p)

ℛ^̣̣:   h_̣̣(p)    = 1/𝒩 ∑_k P^_(p+k) P^_(-k) =  1/𝒩 ∑_k P^_(p+k) P^_(-k)

= h(p)

ℛ^̣̆:   h_̣̆(p)    = 1/𝒩 ∑_kmn U^_m (p+k) U^_n (-k)  U^*̆_m (p+k) U^*̆_n (-k)  =1/𝒩 ∑_k P^_(p+k)  P^_(k)

=  1/𝒩 ∑_k P_(p+k) P_(-k)^*

ℛ^̣̆:   h_̣̆(p)    = 1/𝒩 ∑_k P^_(p+k)  P^_(-k) = h_̣̆(-p)^*

so we see that the two branches of spin $0$ Goldstone modes have the same spectrum as the spin 0 Cooper pair, $h(\mbf{p})$ in eq:CPham. This is evidence of an enlarged symmetry group, which we will explain in forthcoming work. In eq:harrows, the Hamiltonians $h^{\u\d}(\mbf{p})$ of the spin $\pm1$ density excitations are related by $\mathcal{T}$ which flips the spins. If the single-particle bands have $SU(2)$ spin symmetry, then $P(\mbf{k}) = P(-\mbf{k})^*$ and $h_{\u \d}(\mbf{p}) = h(\mbf{p})$. This symmetry will manifest itself in the excitation spectrum, where there will be particle-hole symmetric excitations about the flat band continuum at $E = \eps |U|$.  Otherwise, the spectrum of $h_{\u \d}(\mbf{p})$ is generically different from that of $h(\mbf{p})$. Because our interest in this work is the low-energy excitations, we postpone a detailed study of $h_{\u \d}(\mbf{p})$ to future work.

We now write down the spin-0 density operators as in eq:etap:
ρ_p,μ,σ    = 1/√(𝒩 _μ(p)) ∑_kmn γ^†_p+k,m,σγ_k,n,σ 𝒰^σ,-σ_kmn,(p) u^_μ(p)

= 1/√(𝒩 _μ(p)) ∑_kmn γ^†_p+k,m,σγ_k,n,σ U^σ*_m (p+k) U^-σ, *_n (-k)  u^_μ(p)

= 1/√(𝒩) ∑_k u^_μ(p)/√(_μ(p)) c̅^†_p+k,,σc̅_k,,σ

= 1/√(𝒩) ∑_R u^_μ(p)/√(_μ(p))  e^-i p ·(R+r_) n̅_R,,σ  .

To compute the normalized states, we use Wick's theorem. To start, we require the correlator
⟨z| γ^†_k',n',σ γ_q+k',m',σ γ^†_p+k,m,σγ_k,n,σ|z|   ⟩= δ_pq/(1+|z|^2)^2 |z|^4 δ_p,0 δ_n'm' δ_nm - 0 + |z|^2 δ_k,k' δ_nn' δ_mm' ⟨z|z| ⟩.

The presence of the $\delta_{\mbf{p},0}$ term is similar to eq:CPnormcorr. Following the same steps as in eq:etapnorm and using eq:calUupc, we find

⟨z|ρ^†_q,ν ρ_p,μ|z⟩    = δ_p,q (1+|z|^2)^N_f 𝒩-2 /𝒩 √(_μ(p)_ν(q)) ∑_kk',  𝒰^σ,-σ_k'm'n',(q)^* u^*_ν(q) |z|^4 δ_n'm' δ_p,0 δ_nm + |z|^2 δ_kk' δ_nn' δ_mm' 𝒰^σ,-σ_kmn,(p) u^_μ(p)

= δ_p,q (1+|z|^2)^N_f 𝒩-2 /√(_μ(p)_ν(p)) ∑_ u^*_ν(p) |z|^4 δ_p,0 𝒩 ^2 +|z|^2 h_ (p)u^_μ(p)

=  δ_p,q δ_μ,ν (1+|z|^2)^N_f 𝒩-2 |z|^4 δ_p,0 δ_μ,0 N_f 𝒩  + |z|^2  .

The case of $\mbf{p}=0,\mu=0$ again corresponds to the symmetry $\bar{N}_\u$ and $N̅_$̣ operators at zero energy. Otherwise, we compute
∑_n |z|^2n N_f𝒩n ⟨n|ρ^†_q,ν ρ_p,μ|n⟩    = ∑_n |z|^2n+2 δ_p,q δ_μ,ν  N_f𝒩-2n ,     p≠0, μ≠0

⟨n| ρ^†_q,ν ρ_p,μ |n⟩    = ν(1-ν) δ_p,q δ_μ,ν

which fixes the normalization of the states.

§ SPECTRUM AND TOPOLOGY OF THE COOPER PAIRS

In this section, we study the effective single-particle boson Hamiltonian $h(\mbf{p})$ (eq:hpairing of the Main Text) in detail, first showing clearly the invariance of all many-body observables under the momentum space gauge freedom corresponding to the orbital positions (eq:conventons). We then prove a number of generic features of the spectrum of $h(\mbf{p})$ (app:spectrumbounds), demonstrating a unique low-energy branch corresponding to the $s$-wave Cooper pairing channel. Its effective mass is determined by the minimal quantum metric (app:quantummetric). We then prove the symmetry properties of $h(\mbf{p})$ and determine the many-body symmetry transformations of the Cooper pair (eq:CPsymmetry). Finally, we use topological quantum chemistry to enumerate the possible band connectivities and discover fragile topological Cooper pair bands (app:bosontopology).

§.§ Momentum Space Conventions

Before discussing the spectrum of $h(\mbf{p})$, we discuss the dependence of $h(\mbf{p})$ on the choice of momentum space operator convention in eq:cRck. We start by recalling the single-particle Hamiltonian
H̃    = ∑_RR',,σ c^†_R,,σ t^σ_(R - R') c_R',,σ = ∑_k,σ c^†_k,,σ h̃^σ_(k) c_k,,σ

where the momentum operators are $c^\dag_{\mbf{k},\al,\sigma} = \frac{1}{\sqrt{\mathcal{N}}} \sum_\mbf{R} e^{- i \mbf{k} \cdot (\mbf{R} + \mbf{r}_\al)} c^\dag_{\mbf{R},\al,\sigma}$. The choice of the $\mbf{r}_\al$ term in the momentum space operators is arbitrary but convenient for representations of the space group symmetries [Wieder et al., 2018, Herzog-Arbeitman et al., 2022]. We now show that the spectrum of $h(\mbf{p})$ and the many-body operators $\eta^\dag_{\mbf{p},\mu}$ are all invariant under $\mbf{r}_\al \to \mbf{r}_\al + \mbf{x}_\al$. However, under $\mbf{r}_\al \to \mbf{r}_\al + \mbf{x}_\al$, we see nontrivial transformations
c^†_k,,σ    →e^- i k ·x_ c^†_k,,σ

h̃^σ_(k)     →e^i k ·x_ h̃^σ_(k) e^-i k ·x_  .
If we define the local embedding matrix $[V_\mbf{x}(\mbf{k})]_{\al \be} = e^{-i \mbf{k} \cdot \mbf{x}_\al} \delta_{\al \be}$, then the transformation can be written $\tilde{h}(\mbf{k}) \to V_\mbf{x}^\dag(\mbf{k}) \tilde{h}(\mbf{k}) V_\mbf{x}(\mbf{k})$. Thus the electron eigenvectors $U_\sigma(\mbf{k})$ obey $U_\sigma(\mbf{k}) \to V_\mbf{x}^\dag(\mbf{k}) U_\sigma(\mbf{k})$ and the projectors transform like $\tilde{h}$: $P_\sigma(\mbf{k}) \to V^\dag_\mbf{x}(\mbf{k}) P_\sigma(\mbf{k})V_\mbf{x}(\mbf{k})$. As such, we have
P^σ_(k) →e^i k ·(x_-x_) P^σ_(k)  .

We now consider the transformation of the Cooper pair Hamiltonian:
h_(p)    →1/𝒩 ∑_k e^i (p+k) ·(x_-x_) e^i k ·(x_-x_) P_(p+k) P_(k)

= e^i p ·(x_-x_) 1/𝒩 ∑_k P_(p+k) P_(k)

= [V^†_x(p) h(p) V_x(p)]_  .
Hence the Cooper pair eigenvectors transform just like the electron eigenvectors: $u_\mu(\mbf{p}) \to V^\dag_\mbf{x}(\mbf{p}) u_\mu(\mbf{p})$. We see that the Cooper pair operators are invariant under the embedding:
η^†_p,μ    = 1/√(𝒩) ∑_k u^_μ(p)/√(_μ(p)) c̅^†_p+k,,c̅^†_-k,,

→1/√(𝒩) ∑_k u^_μ(p) e^i p·x_/√(_μ(p)) e^-i (k+p)·x_ c̅^†_p+k,, e^i k·x_ c̅^†_-k,,

= η^†_p,μ

It is also obvious that the spectrum of the Cooper pairs, $\eps_\mu(\mbf{p})$, is invariant because the $V_\mbf{x}(\mbf{k})$ conjugate the Cooper pair Hamiltonian. Hence, as expected, all of the many-body observables are invariant under the arbitrary choice of $\mbf{x}_\al$.

§.§ Cooper Pair Hamiltonian and Spectrum

The fact that the Cooper pair and density excitations are determined by an effective single-particle Hamiltonian $h(\mbf{p})$ is very evocative. It shows that the tightly-bound Cooper pairs behave much like a free boson whose hoppings arise from the quantum geometry, encoded in $P(\mbf{k})$, of the flat bands. We will study the effective single-particle boson Hamiltonian and prove a number of properties about its spectrum and topology.

Let us restate the effective Hamiltonian derived in eq:CPham,
h_(p) = 1/𝒩 ∑_k P_(p+k) P_(k)
where $P(\mbf{k})= U(\mbf{k})U^\dag(\mbf{k})$ is the Hermitian projector into the spin-$\u$ electronic flat bands. The $,$ orbital indices range from $1,…, N_orb$. In the bipartite crystalline lattice construction, $N_orb →N_L$ is understood to be the number of orbitals in the $L$ sublattice. ($P_(k) = 0$ for $$ or $∉L$.)

First we show $h(p)$ is Hermitian:
\bea
h^*_{\be \al}(\mbf{p}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P^*_{\be \al}(\mbf{p}+\mbf{k}) P^*_{\al \be}(\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{p}+\mbf{k}) P_{\be \al}(\mbf{k}) = h_{\al \be}(\mbf{p}) \ .
\eea
The eigenvalues of $h(p)$ are denoted by $_μ(p)$ where $μ= 0, …, N_orb - 1$ will denote the boson band index. We will now show that the spectrum of $h(p)$ is bounded by $0 ≤_μ(p) ≤$. Importantly, the many-body energies (see \Eq{eq:lowrank}) of the boson excitations are given by $|U| (- _μ(p) )$ so the bounds on the eigenvalues of $h(p)$ show that the many-body spectrum of the boson bands is between $0$ and $|U|$. Note that the largest eigenvalue of $h(p)$ corresponds to the lowest energy in the many-body spectrum. Thus we can think of $|U|_μ(p)$ as the binding energy gained by the pairing relative to the unpaired continuum of states at energy $|U|$.

To show $h(p)$ is positive semi-definite, we observe that for an arbitrary vector $v$,
\bea
\sum_{\al \be} v^*_\al h_{\al \be}(\mbf{p}) v_\be &= \frac{1}{\mathcal{N}} \sum_\mbf{k} \sum_{mn,\al\be} v^*_\al U_{\al m}(\mbf{p}+\mbf{k}) U^*_{\be m}(\mbf{p}+\mbf{k})  U_{\be n}(\mbf{k}) U^*_{\al n}(\mbf{k})v_\be \\
&= \frac{1}{\mathcal{N}} \sum_\mbf{k} \sum_{mn}  [\sum_\al U^*_{\al m}(\mbf{p}+\mbf{k})  v_\al U_{\al n}(\mbf{k})]^* \sum_\be U^*_{\be m}(\mbf{p}+\mbf{k}) v_\be U_{\be n}(\mbf{k})  \\
&= \frac{1}{\mathcal{N}} \sum_\mbf{k} \sum_{mn} V^*_{mn}  V_{mn}  \\
&= \frac{1}{\mathcal{N}} \sum_\mbf{k} ||V||_F^2 \geq 0 \ .
\eea
where we defined the $N_f ×N_f$ matrix $V_mn = ∑_U^*_m(p+k) v_U_n(k)$ and used the Frobenius norm $||V||_F^2 = ∑_mn |V_mn|^2$. This establishes $_μ(p) ≥0$. Alternatively, we observe that $P_(p+k) P_(k) = [P(p+k) ∘P^T(k)]_$ is the Hadamard product of two positive semi-definite matrices $P(p+k)$ and $P^T(k)$. By the Schur product theorem, $P_(p+k) P_(k)$ is positive semi-definite, and so $h_(p)$ is the sum of positive semi-definite matrices with positive coefficients, and hence is positive-semi definite itself.

To prove $_μ(p) ≤$, we show that $1 - h(p)$ is positive semi-definite. Using the uniform pairing condition, note that
\bea
\label{eq:zerobound}
\null [\eps \mathbb{1} - h(\mbf{p})]_{\al \be} &= \eps \delta_{\al \be} \,- h_{\al \be}(\mbf{p})  = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \al}(\mbf{k}) \delta_{\be\al} \ - h_{\al \be}(\mbf{p}) \\
&= \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{p}+\mbf{k}) \delta_{\be\al} \ - h_{\al \be}(\mbf{p}) \\
&= \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{p}+\mbf{k}) (\delta_{\be\al} - P_{\be \al}(\mbf{k})) \\
&= \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{p}+\mbf{k}) Q_{\be \al}(\mbf{k}) \\
\eea
where we defined $Q(k) = 1 - P(k)$. Note that $Q(k)$ is also a positive semi-definite matrix (it is a projector), and so $1- h(p)$ is positive semi-definite by the Schur product theorem.

We have shown that $_μ(p) ≤$. We now prove that this bound is tight by exhibiting the corresponding eigenvector at $p=0$. Claim: there is an eigenvector $u_0^(p=0) = 1/√(N_orb)$ with eigenvalue $$. This is simply the normalized uniform eigenvector corresponding to the eta-pairing symmetry $η^†$. The claim is proven using the uniform pairing condition:
\bea
\label{eq:u0mathcaln}
\null [h(0) u_0]_\al &= \frac{1}{\mathcal{N}} \sum_\mbf{k} \sum_\be P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) u_0^\be(0) \\
&= \frac{1}{\sqrt{N_{orb}}} \frac{1}{\mathcal{N}} \sum_\mbf{k} \sum_\be P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) \\
&= \frac{1}{\sqrt{N_{orb}}} \frac{1}{\mathcal{N}} \sum_\mbf{k} [P(\mbf{k})^2]_{\al \al} \\
&= \frac{1}{\sqrt{N_{orb}}} \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \al}(\mbf{k}) \\
&= \eps u^\al_0 \ . \\
\eea
The existence of this state is tied to the $η^†$ symmetry. The $p=0$ Cooper pair excitation spectrum includes the $η^†$ operator in its Hilbert space, and  $η^†$ commutes with the Hamiltonian and therefore has zero energy. Because the many-body energy is $|U|(- _μ(p))$, there must be a band with $_μ(0) = $, and as such this excitation is gapless. We refer to this as the $μ=0$ band. Indeed, the excitation operator defined in \Eq{eq:etap} yields
\bea
\eta^\dag_{\mbf{p}=0,\mu=0} &= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k}\al} \frac{1/\sqrt{N_{orb}}}{\sqrt{\eps}} \bar{c}^\dag_{\mbf{k},\al,\u}\bar{c}^\dag_{-\mbf{k},\al,\d} = \frac{1}{\sqrt{N_f\mathcal{N}}} \sum_{\mbf{k}\al} \bar{c}^\dag_{\mbf{k},\al,\u}\bar{c}^\dag_{-\mbf{k},\al,\d} = \frac{1}{\sqrt{N_f \mathcal{N}}} \eta^\dag
\eea
using $= N_f/N_orb$. (In the bipartite case, $N_orb →N_L$.)

We now show that the gapless Cooper pair band is unique, i.e. that there are no other zero-energy excitations for all $p$. To do so, we first Fourier transform into real space using \Eq{eq:projdef}, which reads
\bea
P_{\al \be}(\mbf{k}) = \sum_{\mbf{R}} e^{i \mbf{k} \cdot (\mbf{R}+\mbf{r}_\al-\mbf{r}_\be)} p_{\al \be}(\mbf{R}) \ .
\eea
A direct calculation gives
\bea
\label{eq:hprealspace}
h_{\al \be}(\mbf{p}) &= \sum_{\mbf{R}\mbf{R}'} \frac{1}{\mathcal{N}} \sum_\mbf{k} e^{i (\mbf{p}+\mbf{k}) \cdot (\mbf{R}+\mbf{r}_\al-\mbf{r}_\be) + i \mbf{k} \cdot (\mbf{R}'+\mbf{r}_\be-\mbf{r}_\al)} p_{\al \be}(\mbf{R}) p_{\be \al}(\mbf{R}') \\
&= \sum_{\mbf{R}} e^{i \mbf{p} \cdot (\mbf{R}+\mbf{r}_\al-\mbf{r}_\be)} p_{\al \be}(\mbf{R}) p_{\be \al}(-\mbf{R}) \\
&= \sum_{\mbf{R}} e^{i \mbf{p} \cdot (\mbf{R}+\mbf{r}_\al-\mbf{r}_\be)} |p_{\al \be}(\mbf{R})|^2 \\
\eea
using $p^†(R) = p(-R)$. We recognize the last line of \Eq{eq:hprealspace} as an expression for the matrix elements of the real space hopping matrix $t_(R-R') = |p_(R-R')|^2$. This shows that the effective hoppings felt by the bosons are due to the correlation function $p(R)$ of the flat bands. Note that $spec  t_(R-R') = spec |p_(R-R')|^2 = {spec  h(p) ∀p }$. To show that $_0(p=0)$ is the unique largest eigenvalue of the whole band structure of $h(p)$, the crucial feature is that the hopping elements are non-negative: $ |p_(R-R')|^2 ≥0$. Hence we can apply the Perron-Frobenius theorem, which states than a non-negative irreducible matrix has a \emph{unique} largest positive eigenvalue, which we have shown explicitly is $$. Because $|p_(R-R')|^2$ is non-negative, we must also check that $t_(R-R')$ is irreducible to apply the theorem. An irreducible matrix $A$ is one for which $A e ∉e ∀e$ where $e$ is a basis of any coordinate subspace. Physically, irreducibility is equivalent to the lattice being strongly connected, meaning any two orbitals are connected by a sequence of hoppings $|p_(R-R')|^2$. This is \emph{not} the case in a trivial atomic limit where $p_(R-R') ∝δ_RR'$. \Ref{2022PhRvL.128h7002H} proved nonzero lower bounds on $||p(R≠0)||^2$ in terms of real space invariants [Song} et~al., 2020] which are the indices diagnosing obstructed atomic insulators and fragile phases. Excluding certain pathological cases (like a 2D system formed from decoupled 1D chains, or an occupied trivial atomic limit band added with a direct sum into the projector $P(k)$), the results of \Ref{2022PhRvL.128h7002H} show that $t_(R-R')$ is strongly connected if the flat bands are in an obstructed, fragile, or stable phase.

We prove a further result on the spectrum of $h(p)$. So far we have shown that $0 ≤_μ(p) ≤$. We now show that the spectral average of the band structure is $^2$:
\bea
\frac{1}{\mathcal{N}N_{L}} \sum_{\mbf{p}\mu} \eps_\mu(\mbf{p}) = \frac{1}{\mathcal{N}N_{L}} \sum_{\mbf{p}} \Tr h(\mbf{p}) = \frac{1}{\mathcal{N}N_{L}} \sum_{\mbf{p} \al} h_{\al \al}(\mbf{p}) = \frac{1}{\mathcal{N}^2N_{L}} \sum_{\mbf{p}\mbf{k} \al} P_{\al \al}(\mbf{p}+\mbf{k})  P_{\al \al}(\mbf{k}) =  \frac{1}{N_{L}} \sum_{\al}  \eps^2 = \eps^2 \ .  \\
\eea
Recalling the many-body energies $E_μ(p) = |U|(- _μ(p))$, the many-body density of states $g(E)$ obeys
\bea
\label{eq:dosint}
\int_{0}^{\eps|U|} E g(E) dE = \frac{1}{\mathcal{N}} \sum_{\mbf{p}\mu} E_\mu(\mbf{p})  = |U| N_L (\eps - \eps^2) =  |U| N_f (1 - \eps) \ . \\
\eea
The upper bound of integral in \Eq{eq:dosint} is understood to be just below $|U|$, where there is a thermodynamically large degeneracy of flat density excitation bands (analogous to the unpaired continuum in the Cooper pair sector).

Our last result concerns a simple transformation of the model under $P(k) →1-P(k)$, which is equivalent to exchanging the projected and unprojected bands. For simplicity, we assume that all orbitals form an irrep of a single Wyckoff position, e.g. we do not use the bipartite construction in what follows. Note that under $P(k) →1-P(k)$, we have $N_f →N_orb - N_f$ and hence $→1-$. Now we compute
\bea
h_{\al \be}(\mbf{p}) &\to \frac{1}{\mathcal{N}} \sum_\mbf{k} (\delta_{\al \be} - P_{\al \be}(\mbf{k}+\mbf{p}))(\delta_{\al \be} - P_{\be \al}(\mbf{k})) \\
&=\delta_{\al \be} - \delta_{\al \be}  \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al\al}(\mbf{k}+\mbf{p}) - \delta_{\al \be}  \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al\al}(\mbf{k}) + h_{\al \be}(\mbf{p}) \\
&=\delta_{\al \be} - 2 \eps \delta_{\al \be}+ h_{\al \be}(\mbf{p}) \ . \\
\eea
Thus the many-body energies transform as
\bea
E_\mu(\mbf{p}) = |U|(\eps - \eps_\mu(\mbf{p})) \to  |U|\big(1 - \eps - (1 - 2 \eps +\eps_\mu(\mbf{p})) \big) = |U|(\eps - \eps_\mu(\mbf{p})) \\
\eea
and hence are invariant. The energy of the unpaired continuum states does change, however, $|U| →(1-) |U|$. Because we proved above (\Eq{eq:zerobound}) that the bound state energies $E_μ(p) ≤|U|$ are below the unpaired continuum, we see in fact that $E_μ(p) ≤min{|U|, (1-) |U|}$. This follows by taking $P(k) →1-P(k)$ and using $E_μ(p) →E_μ(p)$, $→1-$.

\subsection{Cooper Pair mass and Minimal Quantum Metric}
\label{app:quantummetric}

We proved in \App{app:spectrumbounds} that the many-body spectrum contains a gapless mode $E_μ=0(p=0) = 0$. In this section, we use perturbation theory to compute the small $p$ expansion around this point, yielding the low energy behavior of the bosonic excitations. The result we will prove is $_0(p) = - g_ij/N_L p_i p_j$ where $g_ij$ is the minimal quantum metric (introduced in \Ref{2022arXiv220311133H} and explicitly constructed below) and $N_L$ is the number of orbitals where $P(k)$ is nonzero. We sum over repeated spatial indices $i,j$ in this section. As a result, the many-body energy is
\bea
E_0(\mbf{p}) = |U|(\eps - \eps_0(\mbf{p})) = |U| \frac{g_{ij}}{N_{L}} p_i p_j + \dots
\eea
which shows the low-lying Goldstone and Cooper pair modes are quadratic. This amounts to a many-body calculation of the Cooper pair mass, which has been previously approximated from the two-body problem[T{\"o}rm{\"a}} et~al., 2018] and is directly related to the mean-field superfluid weight [Peotta} and {T{\"o}rm{\"a}}, 2015, T{\"o}rm{\"a} et~al., 2022].

Because we have proven that $_0(p=0)$ is the unique maximal eigenvalue in \App{app:spectrumbounds}, its small $p$ corrections can be obtained using non-degenerate perturbation theory. At zeroth order, $h(0)$ has eigenvectors $u_μ$, with $u_0 = 1/√(N_L)$ corresponding to the maximal eigenvalue $_μ= $. We determine the corrections to $h(0)$ by expanding in $p$:
\bea
h_{\al \be}(\mbf{p}) &= \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) + p_i \frac{1}{\mathcal{N}} \sum_\mbf{k} \del_i P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) + \frac{1}{2} p_i p_j \frac{1}{\mathcal{N}} \sum_\mbf{k} \del_{ij} P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) + \dots  \ . \\
\eea
Denote $h^i_ = _i h_ |_p=0$. By time-reversal $h^*_(p) = h_(-p)$, we expand to first order and find
\bea
(p_i h^i_{\al \be} )^* &= - p_i h^i_{\al \be}, \\
\eea
so $h^i_ = - h^i_$ is anti-symmetric because $(h^i_)^* = h^i_$ using Hermiticity. This can also be seen with an integration by parts. This immediately shows that the first order correction to $_0(p)$ vanishes because
\bea
u^\dag_0 h^i u_0 &= \frac{1}{N_L} \sum_{\al \be} h^i_{\al \be} = 0 \ . \\
\eea
We now consider the second order correction. There are contributions to the second order correction:
\bea
\label{eq:pertubtwoerms}
\frac{1}{2} p_ip_j \del_{ij} \eps_0(\mbf{p}) &= \frac{1}{2} p_ip_j \lp  2 \Re \sum_{\mu \neq 0} \frac{u^\dag_0 h^i u_\mu u_\mu^\dag h^j u_0}{\eps - \eps_\mu(0)} + \frac{1}{\mathcal{N}} \sum_{\mbf{k} \al \be} \frac{1}{N_L}  \del_{ij} P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) \rp \equiv \frac{1}{2} p_i p_j (A_{ij} + B_{ij}) \ .
\eea
The first term $A_ij$ is the usual second order correction from a linear perturbation $p_i h^i$, and the second term $B_ij$ is simply the expectation value of the second order correction in $u_0$. The second term is related to the quantum metric after an integration by parts:
\bea
B_{ij} &= \frac{1}{N_{L}} \frac{1}{\mathcal{N}} \sum_{\mbf{k}\al\be} \del_{ij} P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) =  -\frac{1}{N_{L}} \frac{1}{\mathcal{N}} \sum_\mbf{k}  \Tr \del_i P(\mbf{k}) \del_j P(\mbf{k}) = - \frac{2}{N_{L}} g_{ij} \preceq 0,  \\
g_{ij} &= \frac{1}{\mathcal{N}} \sum_\mbf{k} g_{ij}(\mbf{k}) , \qquad  g_{ij}(\mbf{k}) = \frac{1}{2} \Tr \del_i P(\mbf{k})\del_j P(\mbf{k})
\eea
using the convention $𝒢_ij(k) = P(k)_i P(k) _j P(k) = g_ij(k) - i/2 f_ij(k)$. Note that \Ref{2022PhRvL.128h7002H} proves $P(k) {_i P(k) , _j P(k) } = _i P(k)_j P(k)$. Next, we observe that $A_ij ≽0$ because $-_μ(0) ≥0$. Thus $A$ and and $B$ have competing contributions to $_ij _0(p)$.

Because the spectrum of $h(p)$ is invariant under invariant under $r_→r_+ x_$,  $1/2 p_ip_j _ij _0(p)$ is invariant. However, $A_ij$ and $B_ij$ are not individually invariant. We now show that there is a unique choice of $x_$, up to a constant, where $A_ij = 0$ and the quantum metric alone determines $_ij _0(p)$. Because $A_ij + B_ij$ is independent of $x_$, taking $A_ij →0$ corresponds to the maximal $B_ij$ and hence the minimal $g_ij$ for all $x_$.

First we show that there exists $x_$ where $A_ij = 0$, and then we show such a choice respects the symmetries of the system. Note that we keep $r_$, the physical location of the Wyckoff position, constant.

Returning to \Eq{eq:pertubtwoerms}, we recall that $- _μ> 0$ and $u^†_0 p_i h^i u_μu_μ^†p_j h^j u_0$ is nonnegative. Thus $A_ij = 0$ iff $u_μ^†h^j u_0 = 0$ for all $μ≠0$. Because the $u_μ$ eigenvectors form a complete and orthogonal basis, $u_μ^†h^j u_0 = 0$ iff $ h^j u_0 ∝u_0$. Thus $u_0$ must be an eigenvector of $h^i$. For a given choice of orbital locations, this is not guaranteed. However, $h^i$ is not invariant under shifts of orbital location $r_→r_+ x_$. Using $P_(k) →e^i k ·(x_-x_) P_(k)$, the transformation is
\bea
h^i_{\al \be} =  \frac{1}{\mathcal{N}} \sum_\mbf{k} \del_i P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) &\to \frac{1}{\mathcal{N}} \sum_\mbf{k} \del_i P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) + i (x^i_\al - x^i_\be)  \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{k}) P_{\be \al}(\mbf{k}) \\
&= h^i_{\al\be} + i  (x^i_\al - x^i_\be) h_{\al \be}(0)
\eea
and $x^i_$ denotes the $i$th spatial component of $x_$. Note that $h_(0)$ is Hermitian and real, so $ i (x^i_- x^i_) h_(0)$ is Hermitian and anti-symmetric, as is $h^i_$. If $x_= x$ is a constant shift, then $h^i_ →h^i_$.

Our task is to look for $x_$ such that $u_0$ is an eigenvector of $h^i$ (which guarantees that $A_ij=0$). The eigenvalue equation is
\bea
\label{eq:eigevalcond}
\sum_\be ( h^i_{\al\be} + i  (x^i_\al - x^i_\be) h_{\al \be}(0) ) u^\be_0 = \la u^\al_0 \\
\eea
where $$ is the eigenvalue. First we show that if \Eq{eq:eigevalcond} condition holds, then $= 0$. Contracting both sides with $u^_0 = 1/√(N_L)$ gives
\bea
\la &= \frac{1}{N_{L}} \sum_{\al \be} h^i_{\al\be} + i (x^i_\al - x^i_\be) h_{\al \be}(0) = 0 \\
\eea
because of anti-symmetry in $$. Then \Eq{eq:eigevalcond} reads
\bea
0 &= \sum_\be h^j_{\al\be} + i  \sum_\be (x^j_\al - x^j_\be) h_{\al \be}(0) \\
0 &= \sum_\be h^j_{\al\be} + i x^j_\al \sum_\be h_{\al \be}(0) - i  \sum_\be  h_{\al \be}(0)  x^j_\be \\
\sum_\be  h_{\al \be}(0)  x^j_\be &= -i \sum_\be h^j_{\al\be} + \eps x^j_\al  \\
\null [ (h(0) -\eps) x^j]_\al &= -i \sum_\be h^j_{\al\be}
\eea
where we used the uniform pairing condition $∑_h_(0) = $. We now define the vector $-i ∑_h^j_  = h^j_$ which is real and obeys $∑_h^j_= 0$, so it is orthogonal to the uniform eigenvector $u_0$. Because $h(0)$ has a unique eigenvector $u_0$ with eigenvalue $$ (see \Eq{eq:u0mathcaln}), $(h(0) -)$ is invertible up to the $u_0$ eigenspace. This means $x^j_$ is uniquely determined up to constant shifts because $u_0^∝1$. Thus we obtain $x^j_$ up to constant shifts by applying the pseudo-inverse[Prasolov, 1994] $(h(0) -)^+$:
\bea
\label{eq:xalphafixed}
x^j_\al &= [-i (h(0) -\eps)^+ h^j u_0\sqrt{N_L}]_\al, \qquad [h^j u_0\sqrt{N_L}]_\al = \sum_\be h^j_{\al \be} \\
\mbf{x}_\al &= [i (\eps-h(0))^+ (\pmb{\nabla}h(\mbf{p})|_{\mbf{p}=0}) u_0\sqrt{N_L}]_\al  \ . \\
\eea
The pseudo-inverse of a Hermitian matrix is very simple. For a hermitian matrix $h$ with eigenvalues $_μ$ and eigen-decomposition $h = ∑_μ_μp_μ$, the pseudo-inverse is $h^+ = ∑_μ^+_μp_μ$ where $p_μ$ is the projector onto the $_μ$ eigenspace, $^+ = 1/$ if $≠0$, and $^+ = 0$ otherwise. Hence $hh^+ = h^+h$ acts as the identity on the nonzero eigenspaces and 0 otherwise.

Having proven in \Eq{eq:xalphafixed} that orbitals shifts $x_$ always exist for which $A_ij=0$, we now show that $x_$ preserve the space group symmetries $g ∈G$. In particular, this will show that orbitals at high-symmetry Wyckoff positions, whose location is fixed by symmetry, are necessarily the choices where $A_ij = 0$. In this case, the minimal quantum metric is the quantum metric evaluated in the $r_$ convention \Eq{eq:cRck}.

To prove this result, we need the action of the space group symmetries on $h(p)$. We prove in \App{eq:CPsymmetry} that
\bea
\Gamma[g] h(\mbf{p}) \Gamma^\dag[g] = h(g^{-1} \mbf{p}), \qquad g \in G
\eea
where the representation $Γ_[g]=|D_[g]|$ is a permutation matrix. It follows that
\bea
\label{eq:ghp}
\Gamma[g] \pmb{\nabla}  h(\mbf{p}) \Gamma^\dag[g] = \pmb{\nabla}  h(g^{-1} \mbf{p}) = g^{-1} \pmb{\nabla}h|_{g^{-1} \mbf{p}} \ .
\eea
Evaluating \Eq{eq:ghp} at $p = 0$ gives $Γ[g] ∇  h|_p=0Γ^†[g]  = g^-1 ∇h|_p=0$. Then
\bea
\label{eq:gxgG}
g^{-1} \mbf{x}_\al &= [i (\eps-h(0))^+ (g^{-1}\pmb{\nabla}h(\mbf{p})|_{\mbf{p}=0}) u_0]_\al \\
&= [i (\eps-h(0))^+ \Gamma[g] (\pmb{\nabla}h(\mbf{p})|_{\mbf{p}=0}) \Gamma^\dag[g] u_0]_\al \\
&= \sum_\be \Gamma_{\al \be}[g] [i (\eps-h(0))^+  (\pmb{\nabla}h(\mbf{p})|_{\mbf{p}=0}) [g] u_0]_\be \\
&= \sum_\be \Gamma_{\al \be}[g] \mbf{x}_\be \\
\eea
where we used $Γ[g] h(0) Γ^†[g] = h(0)$ and $Γ[g] u_0=u_0$ because $Γ[g]$ is a permutation matrix and $u_0$ is the uniform eigenvector. Thus if two orbitals are related by a symmetry $g$, their $x_$ that guarantee the minimal metric are also related by $g$. With \Eq{eq:gxgG}, we can prove that $x_=0$ if the site symmetry group $G_x$ of the Wyckoff position $x$ contains a symmetry $h$ such that $h v ≠v$ for all $v≠0$ and $Γ[h] = 1$. As an example, $h$ can be any nontrivial rotation at a maximal Wyckoff position. In this case
\bea
h \mbf{x}_\al =  \sum_\be \Gamma_{\al \be}[h] \mbf{x}_\be = \mbf{x}_\al,
\eea
but by assumption, the only solution to $h x_= x_$ is $x_= 0$. We give two examples in space group $p4mm$ (see \Fig{fig:wyckoff4}). (1) Consider an $s$ (or $p$) orbital at the 2c $= {1/2 x̂,1/2 ŷ}$ position. The site symmetry group is $G_2c = 2mm$ consisting of $C_2$ reflection and mirrors. The $s$ orbitals transform in a 1D irrep, so $Γ[C_2] = 1$. Thus $-x_= C_2 x_ = x_$, so $x_= 0$. (2) Consider $p_x,p_y$ orbitals at the 1b $=1/2 x̂+1/2 ŷ$ position. The representation of $C_4$ is $D[C_4] = i σ_2$ (the $p_x,p_y$ orbitals are rotated into each other), so $Γ[C_4]= σ_1$. Then $Γ[C_2]= 1$ and $-x_= C_2 x_= x_$ as before. We remark that if the orbitals are not at a maximal Wyckoff position (such as the 6j position in \App{app:Aaron}), the symmetries allow $x_$ to be nonzero.

\begin{figure*}
\includegraphics[height=.3\textwidth]{p4wyckoff.pdf}
\caption{Unit cell of the wallpaper group $p4mm$ generated by $C_4$ and translations (we take $\mbf{a}_1 = (1,0),\mbf{a}_2 = (0,1)$ for ease). There are three maximal Wyckoff positons 1a = $(0,0)$,1b= $(1/2,1/2)$,2c= $\{(0,1/2),(1/2,0)\}$. The site-symmetry groups are $G_{1a} = 4mm, G_{1b} = 4mm, G_{2c} = 2mm$. In all cases, the rotations of the site-symmetry groups are enough to pin $\mbf{x}_\al = 0$. }
\label{fig:wyckoff4}
\end{figure*}

%\subsection{Features of the higher-energy Cooper Pairs}

%We have shown that the lowest energy Cooper pair (and Goldstone) excitation is at $E=0$, is of $s$-wave character, and has a quadratic dispersion given by the minimal quantum metric. We now prove some similar features of the higher-energy Cooper pairs, which have other pairing symmetries ($p$-wave, $d$-wave, etc) as we will see in \Eq{eq:CPsymmetry}. We focus on features near the $\Gamma$ point.
%Let us consider the boson Hamiltonian at the Gamma point:
%h_{\al \be}(\mbf{p}=0) &= \frac{1}{V} \sum_\mbf{k} P_{\al \be}(\mbf{k})P_{\be \al}(\mbf{k}) = \frac{1}{V} \sum_\mbf{k} |P_{\al \be}(\mbf{k})|^2 > 0
%where the strict inequality occurs because $\sum_\mbf{k} |P_{\al \be}(\mbf{k})|^2 = 0$ iff $P_{\al \be}(\mbf{k}) = 0$ for all $\mbf{k}$, but this impossible if the flat bands are not in a trivial atomic limit. (On a bipartite lattice, $\al,\be \in L$ is assumed.) Hence $h(0)$ has strictly positive entries and the Perron-Frobenius theorem requires there to be a unique largest eigenvalue $\eps$. We denote $0 \leq \eps_1 < \eps$ by the second largest eigenvalue. The gap between $\eps_1$ and $\eps$ is characterized by
%\kappa = \frac{\eps_1}{\eps} \leq \frac{1-m/M}{1+m/M}
%where $M = \max_{\al \be} h_{\al \be}(0), m = \min_{\al \be} h_{\al \be}(0)$. This fundamental inequality on the Perron-Frobenius spectral gap is due to Hopf \cite{hopf}. We now provide a bound on $\kappa$ to study the gap (at zero momentum) to the second Cooper pair branch.
%First we give an upper bound on $M$. Let $|| \cdot ||$ denote the operator norm, which is sub-multiplicative. Then for any normalized eigenvectors
%|| u^\dag h(0) v || \leq ||u^\dag|| \, ||h(0)||  \, ||v|| = \eps
%where we used that $||h(0)|| = \eps$ is the largest eigenvalue of a Hermitian matrix. Choosing $u,v$ to be elementary vectors with zeros everywhere except $\al$, $\be$ respectively, we find $h_{\al \be}(0) \leq \eps$. Thus $M \leq \eps$.
%Finding a lower bound on $m$ is more subtle. This is because in the trivial atomic limit where $P(\mbf{k})$ is diagonal and independent of $\mbf{k}$, clearly $m=0$.
%To find a lower bound on $m$, we use Cauchy-Schwartz to show
%\frac{1}{V} \sum_\mbf{k} |P_{\al \be}(\mbf{k})|^2 \geq \lp \frac{1}{V} \sum_\mbf{k} |P_{\al \be}(\mbf{k})| \rp^2 \geq \left| \frac{1}{V} \sum_\mbf{k} P_{\al \be}(\mbf{k}) \right|^2

\subsection{Space Group Symmetry and Representations of the Boson Bands}
\label{eq:CPsymmetry}

Because $h(p)$ is a single-particle Hamiltonian, its bands can display nontrivial topology. In fact, we will show that topological quantum chemistry can be applied to $h(p)$ to prove a number of nontrivial statements about the possible band decompositions and symmetry representations. To bring this formalism to bear, we must first work out the symmetries of $h(p)$. We denote the $N_L×N_L$ representations of the symmetries of $h_(p)$ by $Γ[g]$.

First we show that $h(p)$ has a bosonic time-reversal symmetry with representation $Γ[𝒯] = K$. We check
\bea
\label{eq:trsboson}
h^*_{\al \be}(\mbf{p}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P^*_{\al \be}(\mbf{p}+\mbf{k}) P^*_{\be \al}(\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{k})  P_{\be \al}(\mbf{p}+\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{k}-\mbf{p})  P_{\be \al}(\mbf{k})  = h_{\al \be}(-\mbf{p}) \ .
\eea
We also need to check the embedding matrices. Using $P(k+G) = V[G] P(k) V^†[G]$ where $V_[G] = δ_ e^ - i r_·G$ is the embedding matrix and $G$ is a reciprocal lattice vector (see the Appendix of \Ref{2022PhRvL.128h7002H}), we compute
\bea
\label{eq:orbitallocation}
h_{\al \be}(\mbf{p}+\mbf{G}) &= \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(\mbf{p}+\mbf{G}+\mbf{k}) P_{\be \al}(\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} e^{-i \mbf{G} \cdot \mbf{r}_\al} P_{\al \be}(\mbf{p}+\mbf{k}) e^{i \mbf{G} \cdot \mbf{r}_\be}P_{\be \al}(\mbf{k}) = [V(\mbf{G}) h(\mbf{p}) V^\dag(\mbf{G})]_{\al \be} \\
\eea
which is the same as the transformation of the electron Hamiltonian. The embedding matrix contains the location $r_$ of the orbital basis in the lattice. Since the embedding matrix $V[G]$ of $h(p)$ is the same as the initial electron Hamiltonian (see \Eq{eq:pRandembeddingV}), we find that the orbitals of the effective boson tight-binding model are located at the same position as the electron orbitals. The representations of the orbitals, however, may be different.

To study the spatial symmetries, we recall that $D[g] P(k) D^†[g] = P(gk)$ (see \Eq{eq:simhtilde}). In indices, we have
\bea
\label{eq:Hreps1}
h_{\al \be}(g\mbf{p}) &= \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(g\mbf{p}+\mbf{k}) P_{\be \al}(\mbf{k}) = \frac{1}{\mathcal{N}} \sum_\mbf{k} P_{\al \be}(g\mbf{p}+g\mbf{k}) P_{\be \al}(g\mbf{k}) \\
&= \frac{1}{\mathcal{N}} \sum_{\mbf{k}, \kappa\kappa'\la\la'} D_{\al \kappa}[g] P_{\kappa \kappa'}(\mbf{p}+\mbf{k})  D_{\be \kappa'}^*[g] D_{\be \la}[g] P_{\la \la'}(\mbf{k})  D_{\al \la'}^*[g] \\
&= \frac{1}{\mathcal{N}} \sum_{\mbf{k}, \kappa\kappa'\la\la'} D_{\al \kappa}[g] D_{\al \la'}^*[g] \ P_{\kappa \kappa'}(\mbf{p}+\mbf{k}) P_{\la \la'}(\mbf{k}) \  D_{\be \kappa'}^*[g] D_{\be \la}[g] . \\
\eea
We now need the following fact which we prove momentarily: in every point group, there exists a choice of basis where
\bea
\label{eq:DDdD}
D_{\al \kappa}[g] D_{\al \la'}^*[g] = \delta_{\la \kappa}|D_{\al \la}[g]|^2, \qquad g \in G \ .
\eea
To understand \Eq{eq:DDdD}, we see that each column of $D[g]$ can only have a single non-zero element. Because $D[g]$ is unitary, this element must be a complex phase. Thus \Eq{eq:DDdD} is equivalent to the statement that $D[g]$ is a complex permutation matrix. This is trivially satisfied if there is a single orbital per unit cell. More generally, a group with this property, that there exists a basis where all irreps are complex permutation matrices, is called a monomial group. It was recently proven that all spinless point groups in 2D and 3D are monomial [Alexandradinata et~al., 2020]. Because $D[g]$ are spinless representations of the point group (or little group) of the $Γ$ point $p=0$, \Eq{eq:DDdD} holds. Conveniently, all the representation matrices given on the Bilbao crystallographic server \url{https://www.cryst.ehu.es/cgi-bin/cryst/programs/representations_point.pl?tipogrupo=spg} are in the complex permutation matrix form so there is only one nonzero entry per row. However, recall that our convention in \App{app:spin} is to take $D[𝒯]$ to act as the identity on the orbital index. It is now convenient to define
\bea
\label{eq:defgamma}
D_{\al \kappa}[g] D_{\al \la'}^*[g] = \delta_{\la \kappa}|D_{\al \la}[g]|^2 = \delta_{\la \kappa}|D_{\al \la}[g]| \equiv \delta_{\la \kappa} \Gamma_{\al \la}[g], \qquad \Gamma_{\al \be}[g] = |D_{\al \be}[g]|
\eea
where $Γ[g]$, which are real permutation matrices, are the representations of $g∈G$ on $h(p)$. This follows from \Eqs{eq:Hreps1}{eq:defgamma}:
\bea
\label{eq:gammag}
h_{\al \be}(g\mbf{p}) &= \frac{1}{\mathcal{N}} \sum_{\mbf{k}, \la\la'} D_{\al \la}[g]D_{\al \la}^*[g]  P_{\la \la'}(\mbf{p}+\mbf{k}) P_{\la' \la}(\mbf{k}) D_{\be \la'}^*[g] D_{\be \la'}[g] \\
&= \frac{1}{\mathcal{N}} \sum_{\mbf{k}, \la\la'} \Gamma_{\al \la}[g]  P_{\la \la'}(\mbf{p}+\mbf{k}) P_{\la' \la}(\mbf{k}) \Gamma_{\be \la'}[g]\\
&= [\, \Gamma[g] h(\mbf{p}) \Gamma^\dag[g] \,]_{\al \be} .
\eea
In this way, we have shown that $h(p)$ possess all the space group symmetries of the electron Hamiltonian, but in a ``bosonized" form. Because $Γ[g]$ are all real permutation matrices, we can think of the boson orbitals as those obtained from replacing all electron orbitals with $s$ orbitals. We emphasize that this result holds on the UPC lattices we have constructed, where the orbitals form an irrep of a single Wyckoff position (\App{app:UPC}).

Because $Γ[g]$ is real, it is obvious that $Γ[𝒯]$ commutes with all spatial operators, so the topology of the Cooper pair bound states is classified by time-reversal symmetric space groups with $Γ[𝒯]^2 = +1$. There is an additional constraint on the spectrum: we proved that $u_0^(p=0) = 1/√(N_L)$. But $Γ[g]$ is a permutation matrix, so $Γ[g] u_0(p=0)  = u_0(p=0)$ and hence the irrep of the $_0(p)$ band at $p=0$ is always the trivial irrep. We now show that in the higher Cooper pair bands and at nonzero momenta, the Cooper pair wavefunction can (and will) exhibit other pairing symmetries.

Using the representation matrices $Γ[g]$, we can compute the action of $g$ on the many-body Cooper pair creation operators. From \Eq{eq:etap}, we first derive
\bea
g \eta^\dag_{\mbf{p},\mu} g^\dag &= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k}\al\al'\al''} \frac{u^\al_{\mu}(\mbf{p})}{\sqrt{\eps_\mu(\mbf{p})}} \bar{c}^\dag_{g\mbf{p}+g\mbf{k},\al',\u} D^\u_{\al' \al}[g] \bar{c}^\dag_{-g\mbf{k},\al'',\d}D^\d_{\al'' \al}[g] \\
&= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k}\al\al'\al''} \frac{1}{\sqrt{\eps_\mu(\mbf{p})}} \bar{c}^\dag_{g\mbf{p}+\mbf{k},\al',\u} D_{\al' \al}[g]  D^*_{\al'' \al}[g] u^\al_{\mu}(\mbf{p}) \bar{c}^\dag_{-\mbf{k},\al'',\d} \\
&= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k}\al\al'\al''} \frac{1}{\sqrt{\eps_\mu(g\mbf{p})}} \bar{c}^\dag_{g\mbf{p}+\mbf{k},\al',\u} \delta_{\al' \al''} \Gamma_{\al' \al}[g] u^\al_{\mu}(\mbf{p}) \bar{c}^\dag_{-\mbf{k},\al'',\d} \\
\eea
where we used that $D[g] = D^[̆g] = D^[̣g]^*$ by $𝒯$ and $S_z$. Now we use the space group symmetries of $h(p)$ such that $_μ(gp) = _μ(p)$ and $∑_Γ_' [g] u^_μ(p) = ∑_μ' u^'_μ'(gp) B_μ' μ^g(p)$ where $[B^g(p)]_μ'μ = u^†_μ'(gp) Γ[g] u_μ(p)$ is the sewing matrix of the bands [Bernevig and Hughes, 2013]. If $_μ(p) ≠_μ'(p)$, then $B^g_μ' μ(p) = 0$ and thus
\bea
\label{eq:etagtrans}
g \eta^\dag_{\mbf{p},\mu} g^\dag &= \frac{1}{\sqrt{\mathcal{N}}} \sum_{\mbf{k}\al \mu'} \frac{u^{\al}_{\mu'}(g\mbf{p})}{\sqrt{\eps_\mu(g\mbf{p})}} B^g_{\mu' \mu }(\mbf{p}) \bar{c}^\dag_{g\mbf{p}+\mbf{k},\al,\u}  \bar{c}^\dag_{-\mbf{k},\al,\d} = \sum_{\mu'} \eta^\dag_{g\mbf{p},\mu'} B^g_{\mu' \mu}(\mbf{p}) \ . \\
\eea
At the high symmetry points, $B_g(p)$ is the representation matrix of the bands of $h(p)$, and hence is determined by the irreps of the Cooper pair bands. Hence, the symmetry and topology of the single-particle boson Hamiltonian $h(p)$ is detectable through $B_g(p)$ in the many-body states. We emphasize that although the bound states Cooper pairs are always spin-singlets, they can have nontrivial symmetry representations ($p$-wave, $d$-wave, etc) due to their multi-orbital nature. This is perhaps unfamiliar because \emph{with only spin indices}, a Cooper pair wavefunction which is odd under $p →- p$ must be a spin triplet because the spin indices are anti-symmetrized, for example. In contrast, we have shown here that the low-energy Cooper pair bound states, which are spin-singlets, can have be odd under $p →- p$ since their orbital indices are anti-symmetrized even though their spin indices are symmetric.

\subsection{Topology of the Boson Bands}
\label{app:bosontopology}

Let us review the constraints placed on the bands of $h(p)$. First, we chose the electron orbitals to be irreps of a single Wyckoff position in order to guarantee the UPC (\App{app:UPC}), and we showed that the orbitals of the boson Hamiltonian $h(p)$ are at the same Wyckoff position determined by $V[G]$ in \Eq{eq:orbitallocation}. (Note that in the bipartite crystalline lattice construction, the Wyckoff position of the boson Hamiltonian is that of the (larger) $L$ sublattice.) However, the electron orbitals transform in the representation $D[g]$ whereas the boson orbitals transform in the representation $Γ[g]$ (\Eq{eq:gammag}). Since $Γ_[g] = |D_[g]|$, the orbitals of $h(p)$ are effectively $s$ orbitals. $h(p)$ also possesses spinless time-reversal (see \Eq{eq:trsboson}). Lastly, we proved that the lowest energy state at the gamma point $p = 0$ is always the trivial ($s$) irrep.

We now use Topological Quantum Chemistry [Bradlyn} et~al., 2017] to determine the momentum space irreps which are induced from the boson orbitals satisfying the constraints in the previous paragraph. Because the symmetry of the Cooper pair wavefunction at a given $p$ is determined by the induced momentum space irreps, we can completely classify the possible pairing symmetries and Cooper pair band connectivities in every space group.

First we make the momentum space irreps explicit. Consider a high symmetry momentum $p$ satisfying $g p = p + G$ where $G$ is a reciprocal lattice vector. Then  $h(p)$ commutes with $V^†[G] Γ[g]$ and the Cooper pair wavefunctions at $p$ transform in an irrep of the little group operators $V^†[G] Γ[g]$ (see the Appendix of \Ref{2022PhRvL.128h7002H}). By computing the irreps $V^†[G] Γ[g]$ for a given model, one can determine the possible pairing symmetries. The Bilbao Crystallographic server has tabulated all such irreps for every Wyckoff position and space group. For instance, in the kagome model in \App{eq:kagomecooper} (see \Fig{fig:6mmfrag2}), we computed the irreps at the gamma point to be $Γ_1$ (a $s$-wave 1D irrep) and $Γ_5$ (a 2D $d$-wave irrep, angular momentum $±2$). This is because the little group of the gamma point contains the operator $Γ[C_6]$ which is a $3×3$ permutation matrix with eigenvalues $1, e^±2pi i/3$. From the Bilbao Crystallographic server, we observe that the $K$ point has irreps $K_1$ ($s$-wave) and $K_3$ ($d$-wave), and the $M$ point has irreps $M_1$ ($s$-wave) and $M_3 ⊕M_4$ ($p$-wave). Note that the $s$-wave irrep at $Γ$ is at zero energy, so physically this is the mode that condenses, leading to the eta-pairing groundstates over zero momentum and zero angular momentum. However at higher temperature, the higher energy branch is also populated.

\begin{figure*}
\includegraphics[height=.3\textwidth]{p6mm_hoppings}
\includegraphics[height=.3\textwidth]{6mm_fragile}
\caption{Cooper pair band structure of the Hamiltonian \Eq{eq:tildeh6mm} on the kagome sites. ($(a)$) in \App{eq:kagomecooper}. $(b)$ The band representation computed in \App{eq:kagomecooper} is $\mathcal{B}_1 = \Gamma_1 + K_1 + M_1$ for the low-energy blue band, and $\mathcal{B}_2 = \Gamma_5 + K_3 + M_3\oplus M_4$ for the higher energy pair of red bands.}
\label{fig:6mmfrag2}
\end{figure*}

The global structure of the Cooper pair bands in momentum space is also determined by topological quantum chemistry. In general, a set of bands in momentum space can only be gapped if the irreps obey a set of compatibility relations [Bradlyn et~al., 2018]. To apply these compatibility relations to the Cooper pair bands, we recall that the Cooper pair orbitals are induced from a single Wyckoff position, the Cooper bands form an elementary band representation (EBR). If the Wyckoff position is non-maximal, then in general the Cooper pair bands will be trivially decomposable: they can be separated into gapped atomic bands. If the Wyckoff position is maximal (like the $3c$ kagome position in $G = p6mm$), then there are two possibilities: the EBR could be indecomposable meaning that all Cooper pair bands are connected, or the EBR could be decomposable. If the EBR is decomposable, it can be split into disconnected groups of bands, and necessarily at least one band in this decomposition is topological [Cano} et~al., 2018, Bradlyn} et~al., 2018].
The Bilbao crystallographic server tabulates these properties for all EBRs, and explicitly writes the possible irreps of the split EBR. Finally, there is an additional constraint due to time-reversal symmetry which prevents stable topological bands [Brouder et~al., 2007] in one, two, and three dimensions. Thus if a branch of a split EBR has a nonzero symmetry indicator (which would diagnose stable topology if the band were gapped everywhere), there must be a Weyl node in the spectrum which prevents the bands from being gapped everywhere [Hughes} et~al., 2011, Elcoro} et~al., 2020, Bradlyn} et~al., 2017]. However, the Cooper pair bands can be gapped with nontrivial topology if the split EBR has a decomposition into fragile topological bands.

Using the Bilbao Crystallographic server, we exhaustively tabulate every EBR in the 230 space groups which is decomposable and obeys the constraints in the first paragraph of \App{app:bosontopology} (notably, that there are $s$ orbitals on the Wyckoff position). We leave a complete study of the magnetic space groups to other work.

\begin{figure*}
\includegraphics[height=.3\textwidth]{wyckoffhex.pdf}
\caption{We show the maximal Wyckoff positions in $p6$ and $p6mm$. The 1a position forms a triangular lattice, the 2b position forms a hexagonal lattice, and the 3c position forms a kagome lattice. The site symmetry groups of the positions are $G_{1a} = 6, G_{2b} = 3, G_{3c} = 2$. }
\label{fig:wykhex}
\end{figure*}

To illustrate the method, we first consider the 2D wallpaper groups with spinless time-reversal symmetry. Almost every Wyckoff position is indecomposable, meaning that the bands are all connected through high symmetry points in the BZ. Of the 17 spinless wallpaper groups, only $p61'$ and $p6mm1'$ have decomposable bands induced from the trivial irrep at a maximal Wyckoff position. In both cases, this Wyckoff position is the 3c position (see \Fig{fig:wykhex}) which has site symmetry group $G_3c = 2 = {1,C_2}$. Putting $s$ or $p$ electron orbitals at the 3c position yields the Cooper pair representations $Γ[C_2] = 1$ and $Γ[C_6] = P_3$ where $P_3$ is the permutation matrix exchanging the 3 sites. This is the representation of an $s$-orbital at the 3c position, denoted $A_3c ↑G$.  From the Bilbao crystallographic server we find that in $G=p61'$ (the $1'$ denotes time-reversal), there are two possible band decompositions
\bea
\label{eq:banddecomp}
A_{3c} \uparrow p61' &= \Gamma_1 \oplus \Gamma_3\Gamma_5 + K_1 \oplus K_2K_3 + M_1 \oplus M_2 \oplus M_2 \\
&= \begin{cases}
(\Gamma_1 + K_1 + M_1) +  (\Gamma_3\Gamma_5 + K_2K_3 + 2M_2) & \\
(\Gamma_1 + K_1 + M_2) +  (\Gamma_3\Gamma_5 + K_2K_3 + M_1\oplus M_2) & \\
\end{cases}
\eea
where the parentheses denote which bands are connected. Both decompositions contain a single band with a $Γ_1$ irrep (the trivial) irreps, which is the quadratic band at low energy. The two higher bands are connected and gapped from the low energy band. The momentum space irreps are defined
\bea
\begin{array}{c|ccc}
61' &1& C_6   \\
\hline
\Gamma_1 &1&1 \\
\Gamma_3\Gamma_5 &2 & -1   \\
\end{array}, \qquad
\begin{array}{r|rr}
3 &1& C_3 \\
\hline
K_1 & 1 & 1 \\
K_2 K_3 & 2 & -1 \\
\end{array},\qquad
\begin{array}{c|cc}
m&1& M_x \\
\hline
M_1 &1&1\\
M_2 &1 & -1 \\
\end{array} \ .
\eea
We now consider the topology of the two branches in \Eq{eq:banddecomp}. We observe that the $Γ_1 + K_1 + M_1$ band is $A_1a ↑p61'$ and hence is a trivial atomic band. Its complement $Γ_3Γ_5 + K_2K_3 + 2M_2 = A_3c ↑G ⊖A_1a ↑G$ is fragile topological. This is an instance of the general result that gapping band induced from a maximal Wyckoff position must result in topology [Cano} et~al., 2018, Bradlyn} et~al., 2018]: since the lower band is trivial, the higher two bands must be topological. The other allowed decomposition has a $Γ_1 + K_1 + M_2$ band. Using the Smith normal form [Elcoro} et~al., 2020], we compute the symmetry indicator $þ= m(Γ_2) + m(M_2) 2$, which is the number of odd inversion eigenvalues mod 2. We observe that $þ= 1$ for $Γ_1 + K_1 + M_2$, indicating is must have Weyl nodes connecting it to the other bands.

We systematize this information in the following table. W.P. stands for the Wyckoff position of the orbitals. The $η$ irreps column denotes the orbital of $h(p)$ of the site symmetry group of the Wyckoff position (physically it is always an $s$ orbital, but is denoted with different symbols in different space groups), and $e^-$ irrep denotes the possible electron orbitals that yield the bosonic orbital through $Γ[g] = |D_[g]|$. $#$ numbers the possible decompositions and $þ$ is the symmetry indicator, which results in Weyl nodes due to $Γ[𝒯]$. Our classification enumerates the topology of the Cooper pair bands (the first entry corresponds to the gapless quadratic band). When the band is trivial, its atomic representation is given. Otherwise fragile topology or topologically protected Weyl nodes are indicated.

\begingroup
\label{tab:2Dcooperpairs}
\renewcommand\arraystretch{1.3}
%\begin{longtable}{| C{.04\textwidth} | p{.25\textwidth} | p{.25\textwidth} || p{.23\textwidth} | p{.23\textwidth}|}
\begin{longtable}{|c|c|c|c|c|c|c|c|}
\hline
$G$ & \text{name} & W.P. &\text{$\eta^\dag$ irreps ($e^-$ irreps)} & \# &$\theta$ &\text{Classification}  \\
\hline
16 & $p6$ &  & & &$\th = m(\Gamma_2) + m(M_2) \mod 2 \in \mathds{Z}_2  $ & \\
\hline
& & 3c & $A \  (A,B) $ & 1 &$0,0 $ & $A^{1a} \uparrow G$, fragile \\
& & 3c&  $A \  (A,B) $& 2 &$1,1$ & Weyl, Weyl\\
\cline{3-7}
\hline
17 & $p6mm$ &  & & & none & \\
\hline
& & 3c & $A_1 \  (A_1,A_2,B_1,B_2) $ & 1 & & $A_1^{1a} \uparrow G$, fragile \\
\hline
\end{longtable}
\endgroup
We see that in 2D, the only cases with decomposable bands at maximal Wyckoff positions are in $p6$ and $p6mm$. In all other groups and maximal Wyckoff positions, the Cooper pair bands are connected. In the decomposable cases listed above, only the higher energy bands can be fragile topological.

Performing an exhaustive search 3D, we find 22 space groups containing decomposable bands with a trivial irrep at the $Γ$ point. In the case where the symmetry indicator is zero, we check whether the band is fragile or atomic. They are enumerated in the following table. There are many cases where the low-lying Cooper pair band has topologically protected Weyl nodes. We also highlight $G = Fd3̅c$, space group 228, where the only allowed decomposition is into fragile bands.

%SG 191 $P6/mmm$ with candidate CaNi_5.  (8 others, but lowest band is trivial)
%SG 166 $R\bar{3}m$ with candidate Ca_2NCl.  (8 others, but lowest band is trivial)
%SG 163 $P\bar{3}1c$ with candidate $KAg\[CN\]_2$.
%SG 176 $P6_3/m$ with candidate $RbMo_3S_3$.

%G = 147, 148, 162, 163, 164, 165,166,167,168,175,176,177,183,184,191,192,193,194,202,209,225,228 \\

\begingroup
\renewcommand\arraystretch{1.3}
%\begin{longtable}{| C{.04\textwidth} | p{.25\textwidth} | p{.25\textwidth} || p{.23\textwidth} | p{.23\textwidth}|}
\begin{longtable}{|c|c|c|c|c|c|c|c|}
\hline
$G$ & \text{name} & W.P. &\text{$\eta$ irreps ($e^-$ irreps)} & \# &$\theta$ &\text{Classification}  \\
\hline
147 & $P\bar{3}$ &  & & &$(\mathds{Z}_2, \mathds{Z}_4) $ & \\
\hline
& &  \multirow{4}{*}{3e} &  \multirow{4}{*}{$A_g \  (A_g,A_u) $}  & 1 &$(0,0), (0,0) $ & $A^{1a}_g \uparrow G$, fragile\\
& & &  & 2 &$(0,1), (0,3) $ & Weyl, Weyl\\
& & &  & 3 &$(1,3), (1,1) $ & Weyl, Weyl\\
& & &  & 4 &$(1,0), (1,0) $ & Weyl, Weyl\\
\cline{3-7}
& & \multirow{4}{*}{3f} & \multirow{4}{*}{$A_g \  (A_g,A_u) $}  & 1 &$(0,3), (0,1) $ & Weyl, Weyl\\
& & &  & 2 &$(0,0), (0,0) $ & $A^{1b}_g \uparrow G$, fragile\\
& & &  & 3 &$(1,2), (1,2) $ & Weyl, Weyl\\
& & &  & 4 &$(1,3), (1,1) $ & Weyl, Weyl\\
\hline
148 & $R\bar{3}$ &  & & &$(\mathds{Z}_2, \mathds{Z}_4) $ & \\
\hline
& & \multirow{4}{*}{9d} &  \multirow{4}{*}{$A_g \  (A_g,A_u) $}  & 1 &$(1,3), (1,1) $ & Weyl, Weyl\\
& & &  & 2 &$(0,0), (0,0) $ & $A^{3b}_g \uparrow G$, fragile\\
& & &  & 3 &$(1,2), (1,2) $ & Weyl, Weyl\\
& & &  & 4 &$(0,3), (0,3) $ & Weyl, Weyl\\
\cline{3-7}
& & \multirow{4}{*}{9e} &  \multirow{4}{*}{$A_g \  (A_g,A_u) $} & 1 &$(0,0), (0,0) $ & $A^{3a}_g \uparrow G$, fragile\\
& & &  & 2 &$(1,1), (1,3) $ & Weyl, Weyl\\
& & &  & 3 &$(1,3), (1,1) $ & Weyl, Weyl\\
& & &  & 4 &$(1,0), (1,0) $ & Weyl, Weyl\\
\hline
162 & $P\bar{3}1m$ &  & &  & $\mathds{Z}_2$  & \\
\hline
& & 3f & $A_g \  (A_g,A_u,B_g,B_u) $ & 1 & 0 & $A_{1g}^{1a} \uparrow G$, fragile \\
\cline{3-7}
& & 3g & $A_g \  (A_g,A_u,B_g,B_u) $ & 1 & 0 & $A_{1g}^{1b} \uparrow G$, fragile \\
\hline
163 & $P\bar{3}1c$ &  & &  & $\mathds{Z}_2$   & \\
\hline
& & \multirow{4}{*}{6g} &  \multirow{4}{*}{$A_g \  (A_g,A_u) $} & 1 & $1,1$ & Weyl, Weyl \\
& &  &  & 2 & $0,0$ & fragile, fragile \\
& &  & & 3 &$1,1$ & Weyl, Weyl \\
& &  &  & 4 &$0,0$ & $A_g^{2b} \uparrow G$, fragile \\
\hline
164 & $P\bar{3}m1$ &  & &  & $\mathds{Z}_2$ & \\
\hline
& & 3e & $A_g \  (A_g,A_u, B_g, B_u) $ & 1 & 0,0 & $A_{1g}^{1a} \uparrow G$, fragile \\
\cline{3-7}
& & 3f & $A_g \  (A_g,A_u, B_g, B_u) $ & 1 & 0,0 & $A_{1g}^{1b} \uparrow G$, fragile \\
\hline
165 & $P\bar{3}c1$ &  & &  & $\mathds{Z}_2$ & \\
\hline
& & \multirow{2}{*}{6e} &  \multirow{2}{*}{$A_g \  (A_g,A_u) $}  & 1 & $1,1$ & Weyl, Weyl \\
& &  & & 2 & $0,0$ & $A_{g}^{2b} \uparrow G$, fragile \\
\hline
166 & $R\bar{3}m$ &  & &  & $\mathds{Z}_2$ & \\
\hline
& & 9d & $A_g \  (A_g,A_u,B_g,B_u) $ & 1 & $0,0$ & $A_{1g}^{3b} \uparrow G$, fragile  \\
\cline{3-7}
& & 9e & $A_g \  (A_g,A_u,B_g,B_u) $ & 1 & $0,0$ & $A_{1g}^{3a} \uparrow G$, fragile \\
\hline
167 & $R\bar{3}c$ &  & &  & $\mathds{Z}_2$ & \\
\hline
& &  \multirow{2}{*}{18d} &  \multirow{2}{*}{$A_g \  (A_g,A_u) $}  & 1 & $1,1$ & Weyl, Weyl \\
& &  &  & 2 & $0,0$ & $A_g^{6b} \uparrow G$, fragile\\
\hline
168 & $P6$ &  & &  & $\mathds{Z}_2$ & \\
\hline
& &  \multirow{2}{*}{3c} &  \multirow{2}{*}{$A \  (A,B) $}  & 1 & $0,0$ & $A^{1a}\uparrow G$, fragile \\
& &  & & 2 & $1,1$ & Weyl, Weyl \\
\hline
175 & $P6/m$ &  & &  & $(\mathds{Z}_2,\mathds{Z}_2,\mathds{Z}_2)$ & \\
\hline
& &  \multirow{2}{*}{3f} & \multirow{2}{*}{$A_g \  (A_g,A_u,B_g,B_u)$} & 1 & $(0,0,0),(0,0,0)$ & $A_g^{1a}\uparrow G$, fragile \\
& &  & & 2 & $(0,1,1),(0,1,1)$ & Weyl, Weyl \\
\cline{3-7}
& &  \multirow{2}{*}{3g} & \multirow{2}{*}{$A_g \  (A_g,A_u,B_g,B_u)$}  & 1 & $(0,0,0),(0,0,0)$ & $A_g^{1b} \uparrow G$, fragile \\
& &  &  & 2 & $(0,1,1),(0,1,1)$ & Weyl, Weyl \\
\hline
176 & $P6_3/m$ &  & &  & $\mathds{Z}_2$ & \\
\hline
& &  \multirow{2}{*}{6g} & \multirow{2}{*}{$A_g \  (A_g,A_u) $}   & 1 & $1,1$ & Weyl, Weyl \\
& &  & & 1 & $0,0$ & $A^{2b}_g \uparrow G$, fragile \\
\hline
177 & $P622$ &  & &  & none & \\
\hline
& & 3f & $A_1 \  (A_1,B_1,B_2,B_3)$ & 1 &  & $A_1^{1a}\uparrow G$, fragile \\
\cline{3-7}
& & 3g & $A_1 \  (A_1,B_1,B_2,B_3)$ & 1 &  & $A_1^{1b} \uparrow G$, fragile \\
\hline
183 & $P6mm$ &  & &  & none & \\
\hline
& & 3c & $A_1 \  (A_1,A_2,B_1,B_2) $ & 1 &  & $A_1^{1a}\uparrow G$, fragile \\
\hline
184 & $P6cc$ &  & &  & $\mathds{Z}_2$ & \\
\hline
& &  \multirow{2}{*}{3c} & \multirow{2}{*}{$A \  (A,B) $}  & 1 & 0,0 & $A_1^{2a}\uparrow G$, fragile \\
& &  &  & 2 & 1,1 & Weyl, Weyl \\
\hline
191 & $P6/mmm$ &  & &  & none & \\
\hline
& & 3f & $A_g \  (A_g,A_u,B_{1g},B_{1u},B_{2g},B_{2u}) $ & 1 &  & $A_{1g}^{1a}\uparrow G$, fragile \\
\cline{3-7}
& & 3g & $A_g \ (A_g,A_u,B_{1g},B_{1u},B_{2g},B_{2u}) $ & 1 &  & $A_{1g}^{1b}\uparrow G$, fragile \\
\hline
192 & $P6/mcc$ &  & &  & $\mathds{Z}_2$ & \\
\hline
& & 6f & $A_1 \  (A_1, B_1, B_2, B_3) $ & 1 & 0,0 & $A_{1}^{2a}\uparrow G$, fragile \\
\cline{3-7}
& &  \multirow{2}{*}{6g} & \multirow{2}{*}{$A_g \ (A_g, A_u, B_g, B_u) $}  & 1 &  0,0 & $A_{1}^{2a}\uparrow G$, fragile \\
& &  &  & 2 & $1,1$ & Weyl, Weyl \\
\hline
193 & $P6_3/mcm$ &  & &  & none & \\
\hline
& & 6f & $A_g \  (A_g, A_u, B_g, B_u) $ & 1 &  & $A_{1g}^{2b}\uparrow G$, fragile \\
\hline
194 & $P6_3/mmc$ &  & & & none & \\
\hline
& & 6g & $A_g \  (A_g, A_u, B_g, B_u) $ & 1 & & $A_{1g}^{2a}\uparrow G$, fragile \\
\hline
202 & $Fm\bar{3}$ &  & &  &  none& \\
\hline
& &  \multirow{2}{*}{24d} & \multirow{2}{*}{$A_g \  (A_g, A_u) $ } & 1 & & $A_{g}^{4a}\uparrow G$, fragile \\
& &  &  & 2 & & $A_{g}^{4b}\uparrow G$, fragile \\
\hline
209 & $F432$ &  & &  &  none& \\
\hline
& &  \multirow{2}{*}{24d} & \multirow{2}{*}{$A_1 \  (A_1, B_1) $} & 1 & & $A_{1}^{4a}\uparrow G$, fragile, $T_2^{4b} \uparrow G$ \\
& &  &  & 2 & & $A_{1}^{4b}\uparrow G$, fragile, $T_2^{4a} \uparrow G$  \\
\hline
225 & $Fm\bar{3}m$ &  & &  & none& \\
\hline
& &  \multirow{2}{*}{24d} & \multirow{2}{*}{$A_g \  (A_g,A_u, B_{1g},B_{1u}) $}   & 1 & & $A_{1g}^{4a}\uparrow G$, fragile \\
& &  &  & 2 & & $A_{1g}^{4b}\uparrow G$, fragile \\
\hline
228 & $Fd\bar{3}c$ &  & &  & none& \\
\hline
& & 48d & $A \  (A,B) $ & 1 & & fragile, fragile \\
\hline

\end{longtable}
\endgroup

We use the explicit formulae in \Ref{2022arXiv220410556G} for the symmetry indicators [Po}, 2020, Po} et~al., 2017, Kruthoff} et~al., 2017] in the above groups using the Smith normal form[Elcoro} et~al., 2020]. The results are below. The notation is the same as on the \href{https://www.cryst.ehu.es/cgi-bin/cryst/programs/bandrep.pl}{Bilbao Crystallographic Server}.
\begingroup
\renewcommand\arraystretch{1.3}
%\begin{longtable}{| C{.04\textwidth} | p{.25\textwidth} | p{.25\textwidth} || p{.23\textwidth} | p{.23\textwidth}|}
\begin{longtable}{|c|c|c|}
\hline
$G$ & \text{name} & Symmetry Indicator $\th$  \\
\hline
\multirow{3}{*}{147} & \multirow{3}{*}{$P\bar{3}$}  & $(\mathds{Z}_2, \mathds{Z}_4) $ \\
\cline{3-3}
&& $m(M_1^+)+m(\Gamma_1^+) \mod 2 $ \\
&& $m(A_1^+)-2m(A_2^+A_3^+)-m(L_1^+)+m(M_1^+)-m(\Gamma_1^+)+2m(\Gamma_2^+\Gamma_3^+) \mod 4 $ \\
\hline
\multirow{3}{*}{148} & \multirow{3}{*}{$R\bar{3}$}  & $(\mathds{Z}_2, \mathds{Z}_4) $ \\
\cline{3-3}
&& $m(T_1^+)+m(L_1^+) \mod 2 $ \\
&& $m(F_1^+)-m(L_1^+)+m(T_1^+)+2m(T_2^+T_3^+)-m(\Gamma_1^+) + 2m(\Gamma_2^+ \Gamma_3^+) \mod 4 $ \\
\hline
\multirow{2}{*}{162} & \multirow{2}{*}{$P\bar{3}1m$}  & $\mathds{Z}_2$ \\
\cline{3-3}
&& $ m(A_1^-) + m(A_3^-) + m(\Gamma_1^-) + m(\Gamma_3^-) + m(L_2^+) + m(M_2^+) \mod 2 $ \\
\hline
\multirow{2}{*}{163} & \multirow{2}{*}{$P\bar{3}1c$}  & $\mathds{Z}_2$ \\
\cline{3-3}
&& $m(A_3) +m(\Gamma_1^-) + m(\Gamma_3^+) +m(M_2^+) \mod 2 $ \\
\hline
\multirow{2}{*}{164} & \multirow{2}{*}{$P\bar{3}m1$}  & $\mathds{Z}_2$ \\
\cline{3-3}
&& $m(A_1^-)+m(A_3^-) +m(\Gamma_1^-) + m(\Gamma_3^+) +m(L_2^+)+m(M_2^+) \mod 2 $ \\
\hline
\multirow{2}{*}{165} & \multirow{2}{*}{$P\bar{3}c1$}  & $\mathds{Z}_2$ \\
\cline{3-3}
&& $m(A_3)+m(\Gamma_1^-)+m(\Gamma_3^+) + m(M_2^+) \mod 2 $ \\
\hline
\multirow{2}{*}{166} & \multirow{2}{*}{$R\bar{3}m$}  & $\mathds{Z}_2$ \\
\cline{3-3}
&& $m(F_2^+)+m(L_2^+) +m(T_1^-) + m(T_3^+) +m(\Gamma_1^-)+m(\Gamma_3^+) \mod 2 $ \\
\hline
\multirow{2}{*}{167} & \multirow{2}{*}{$R\bar{3}c$}  & $\mathds{Z}_2$ \\
\cline{3-3}
&& $m(F_2^+)+m(\Gamma_2^+)+m(\Gamma_3^-) \mod 2 $ \\
\hline
\multirow{2}{*}{168} & \multirow{2}{*}{$P6$}  & $\mathds{Z}_2$ \\
\cline{3-3}
&& $m(A_1)+m(M_1) \mod 2 $ \\
\hline
\multirow{4}{*}{175} & \multirow{4}{*}{$P6/m$}  & $(\mathds{Z}_2, \mathds{Z}_2,\mathds{Z}_2) $ \\
\cline{3-3}
&& $m(A_1^+)+m(A_1^-)+m(\Gamma_1^+)+m(M_1^-) \mod 2 $ \\
&& $m(A_1^+)+m(A_1^-)+m(M_1^+)+m(M_1^-) \mod 2 $ \\
&& $m(A_1^+)+m(A_1^-)+m(A_2^+)+m(\Gamma_1^+)+m(L_2^+)+m(M_1^-) \mod 2 $ \\
\hline
\multirow{2}{*}{176} & \multirow{2}{*}{$P6_3/m$}  & $\mathds{Z}_2$ \\
\cline{3-3}
&& $m(M_1^+)+m(\Gamma_1^+) \mod 2 $ \\
\hline
\multirow{1}{*}{177} & \multirow{1}{*}{$P622$}  & none \\
\hline
\multirow{1}{*}{183} & \multirow{1}{*}{$P6mm$}  & none \\
\hline
\multirow{2}{*}{184} & \multirow{2}{*}{$P6cc$}  & $\mathds{Z}_2$  \\
\cline{3-3}
&& $m(A_3A_4)+m(M_3) \mod 2 $ \\
\hline
\multirow{1}{*}{191} & \multirow{1}{*}{$P6/mmm$}  & none \\
\hline
\multirow{2}{*}{192} & \multirow{2}{*}{$P6/mcc$}  & $\mathds{Z}_2$  \\
\cline{3-3}
&& $m(A_5) +m(\Gamma_1^-)+m(\Gamma_3^+)+m(\Gamma_5^-)+m(\Gamma_6^-)+m(M_2^+)+m(M_3^+) \mod 2 $ \\
\hline
\multirow{1}{*}{193} & \multirow{1}{*}{$P6_3/mcm$}  & none \\
\hline
\multirow{1}{*}{194} & \multirow{1}{*}{$P6_3/mmc$}  & none \\
\hline
\multirow{1}{*}{202} & \multirow{1}{*}{$Fm\bar{3}$}  & none \\
\hline
\multirow{1}{*}{209} & \multirow{1}{*}{$F432$}  & none \\
\hline
\multirow{1}{*}{225} & \multirow{1}{*}{$Fm\bar{3}m$}  & none \\
\hline
\multirow{1}{*}{228} & \multirow{1}{*}{$Fd\bar{3}c$}  & none \\
\hline
\end{longtable}
\endgroup

%\subsection{Particle-hole symmetric Kinetic Energy}
%As proved in \App{app:etakinetic}, particle-symmetric kinetic energy does not break the eta-pairing symmetry so we can employ a more sophisticated analysis. Because the $\ket{n}$ eta-pairing states remain exact eigenstates in this case, we can compute the excitations on top of them, now obtaining new terms from the kinetic energy.
%The charge +1 excitations are trivial, we simply find
%R_{mn}(\mbf{p}) = \frac{\eps|U|}{2} \delta_{mn} + \tilde{\eps}_m(\mbf{k}) \delta_{mn}
%where we dropped the spin label because of the $SU(2)$ symmetry required when we have $\mathcal{P}$. We see that the condition for the charge gap to remain open is $\max |\eps_n(\mbf{k})| < \nu |U|/2$. We assume the kinetic energy is always small enough to satisfy this bound. The charge +2 spectrum is more interesting. Adding kinetic energy changes the scattering matrix by adding kinetic energy on the diagonal:
%R^{\u \d}(\mbf{p}) &= \eps |U| \mathbb{1} + \mathcal{K} - |U| \frac{1}{V} \mathcal{U}(\mbf{p}) \mathcal{U}^\dag(\mbf{p}), \qquad \mathcal{K}_{\mbf{k}'m'n',\mbf{k}mn}(\mbf{p}) = \delta_{\mbf{k}\mbf{k}'} \delta_{mm'} \delta_{nn'} (\tilde{\eps}_m(\mbf{p}+\mbf{k})+\tilde{\eps}_n(-\mbf{k}) ) \ . \\
%To analyze this problem, we use the matrix determinant lemma:
%\det (D + U V^\dag) &= \det (1+ V^\dag D^{-1} U)\det D \\
%where $D$ is an invertible diagonal matrix. In our application the diagonal elements of $D$ are $\eps|U| + \tilde{\eps}_m(\mbf{p}+\mbf{k})+\tilde{\eps}_n(-\mbf{k})$. Clearly if the charge $+1$ gap is nonzero, then $D$ has non-zero eigenvalues and is invertible.
%We can now use the matrix determinant lemma to obtain an effective zero-body problem for the Cooper pair bands. The eigenvalues $\la$ of the scattering matrix obey
%0 &= \det (D - \la - \frac{|U|}{V} \mathcal{U}\mathcal{U}^\dag)  \\
%&=\det(D- \la) \, \det (1 - \frac{|U|}{V}\mathcal{U}^\dag (D- \la)^{-1} \mathcal{U}) \ .  \\
%We assume that $\det(D- \la) \neq = 0$ so that the inverse matrix is well-defined. This is equivalent to looking at small $\la$, the low-lying many-body energies where we expect the Cooper pair bands to emerge. Another way of thinking about this is to assume $\tilde{\eps}$ is small, so that $D$ is the matrix of two non-interacting charge +1 excitation energies. Assuming $\det(D- \la) \neq = 0$  is then equivalent to looking for the new nontrivial bound states. In this case, we must look for zeros of the other determinant.
%This is governed by the spectrum of the effective single-particle Hamiltonian
%\null [\mathcal{H}_\la(\mbf{p})]_{\al \be} &= |U| \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{P^m_{\al \be}(\mbf{p}+\mbf{k}) P^{-n}_{\be \al}(\mbf{k})}{\eps |U| + \tilde{\eps}_m(\mbf{p}+\mbf{k})+\tilde{\eps}_n(-\mbf{k}) - \la} \\
%&= |U| \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{P^m_{\al \be}(\mbf{p}+\mbf{k}) P^{n}_{\be \al}(\mbf{k})}{\eps |U| + \tilde{\eps}_m(\mbf{p}+\mbf{k})-\tilde{\eps}_{n}(\mbf{k}) - \la} \\
%where $P^m(\mbf{k})$ is the projector into the $m$th electron band with energy $\tilde{\eps}_m(\mbf{k})$. Numerically at least, we can determine the eigenvalues of $\mathcal{H}_\la(\mbf{p})$ as a function of $\la$; whenever such an eigenvalue equals 1, $\la$ is an eigenvalue of the many-body scattering matrix. As a simple check, take $\tilde{\eps} = 0$ so that the denominator is constant and can be factored out of the integral. Then $\mathcal{H}_\la(\mbf{p}) = (\eps - \la/|U|)^{-1} h(\mbf{p})$ whose eigenvalues are $\eps_\mu(\mbf{p})/(\eps - \la/|U|)$. Thus when the eigenvalues of $\mathcal{H}_\la(\mbf{p})$ are equal to 1, we have $\la = |U|(\eps - \eps_\mu(\mbf{p}))$ which is exactly the many-body energy in the flat band case. In this sense, $\mathcal{H}_\la(\mbf{p})$ is a very natural generalization.
%As a start to understanding this matrix, let's check that the uniform eigenvector $u_0$ is still a good approximation of the Cooper pair at the $\Gamma$ point. We will check later that symmetries should enforce this, although it is possible for other eigenvalues to switch in the spectrum. Then we have
%\mathcal{H}_\la(0) u_0 &= \frac{|U|}{\sqrt{N_{orb}}} \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{[P^m(\mbf{k}) P^n(\mbf{k})]_{\al\al}}{\eps|U| + \tilde{\eps}_m(\mbf{k})-\tilde{\eps}_n(\mbf{k}) - \la} \\
%&= \frac{|U|}{\sqrt{N_{orb}}} \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{\delta_{mn}[P^n(\mbf{k})]_{\al\al}}{\eps|U| + \tilde{\eps}_n(\mbf{k})-\tilde{\eps}_n(\mbf{k}) - \la} \\
%&= \frac{|U|}{\sqrt{N_{orb}}} \frac{1}{V} \sum_\mbf{k} \sum_{n} \frac{[P^n(\mbf{k})]_{\al\al}}{\eps|U| - \la} \\
%&= \frac{|U|}{\sqrt{N_{orb}}} \frac{1}{V} \sum_\mbf{k}\frac{[P(\mbf{k})]_{\al\al}}{\eps|U| - \la} \\
%&= \frac{|U|}{\sqrt{N_{orb}}} \frac{\eps}{\eps|U| - \la} \\
%where we used the orthogonality of the projectors. We see that $u_0$ is an eigenvector with eigenvalue $\frac{\eps|U|}{\eps|U| - \la} = 1$ when $\la = 0$. Thus we confirm that the Cooper pair spectrum is still gapless, as it must be.
%We also see that $\mathcal{H}$ is positive semi-definite as long as $\la < \min \ \eps |U| + \tilde{\eps}$ which is the stability criterion on $\det D$. This is because the denominator is positive definite and the numerator $P^m \circ (P^n)^T$ is positive semi-definite due to Schur's theorem. We now check that there is again a bosonic time reversal
%\mathcal{H}_\la(\mbf{p})^* &= |U| \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{P^m_{\al \be}(\mbf{p}+\mbf{k})^* P^{n}_{\be \al}(\mbf{k})^*}{\eps |U| + \tilde{\eps}_m(\mbf{p}+\mbf{k})-\tilde{\eps}_{n}(\mbf{k}) - \la} \\
%&= |U| \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{P^{-m}_{\al \be}(-\mbf{p}-\mbf{k}) P^{-n}_{\be \al}(-\mbf{k})}{\eps |U| + \tilde{\eps}_m(\mbf{p}+\mbf{k})-\tilde{\eps}_{n}(\mbf{k}) - \la} \\
%&= |U| \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{P^{m}_{\al \be}(-\mbf{p}+\mbf{k}) P^{n}_{\be \al}(\mbf{k})}{\eps |U| + \tilde{\eps}_{-m}(\mbf{p}-\mbf{k})-\tilde{\eps}_{-n}(-\mbf{k}) - \la} \\
%&= |U| \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{P^{m}_{\al \be}(-\mbf{p}+\mbf{k}) P^{n}_{\be \al}(\mbf{k})}{\eps |U| - (\tilde{\eps}_{m}(-\mbf{p}+\mbf{k})-\tilde{\eps}_{n}(\mbf{k})) - \la} \\
%To understand the dispersion corrections to the mass, we expand $\null [\mathcal{H}_\la(\mbf{p})]_{\al \be} $ in $\tilde{\eps}$:
%\null [\mathcal{H}_\la(\mbf{p})]_{\al \be} &= \frac{|U|}{\eps |U| - \la}\lp h_{\al \be}(\mbf{p}) + \frac{1}{V} \sum_\mbf{k} \sum_{mn} \frac{P^m_{\al \be}(\mbf{p}+\mbf{k}) P^{n}_{\be \al}(\mbf{k})}{\eps |U| - \la}(\tilde{\eps}_m(\mbf{p}+\mbf{k})-\tilde{\eps}_{n}(\mbf{k}))  \rp\\
%&= \frac{|U|}{\eps |U| - \la}\lp h_{\al \be}(\mbf{p}) + \frac{1}{V} \sum_\mbf{k} \frac{\tilde{h}_{\al \be}(\mbf{p}+\mbf{k}) P_{\be \al}(\mbf{k}) - P_{\al \be}(\mbf{p}+\mbf{k}) \tilde{h}_{\be \al}(\mbf{k})}{\eps |U| - \la} \rp\\
%where $\tilde{h}(\mbf{k}) = \sum_m \eps_m(\mbf{k}) P_m(\mbf{k})$ is the projected single-particle Hamiltonian. At first order in $\eps$, we can take the expectation value:
%\frac{1}{N_{orb}}\frac{1}{V} \sum_\mbf{k} \frac{\Tr \tilde{h}(\mbf{p}+\mbf{k}) P(\mbf{k}) - P(\mbf{p}+\mbf{k}) \tilde{h}(\mbf{k})}{\eps |U| - \la} &= \frac{1}{N_{orb}}\frac{1}{V} \sum_\mbf{k} \frac{\Tr (\tilde{h}(\mbf{k}+\mbf{p}) - \tilde{h}(\mbf{k}-\mbf{p}) ) P(\mbf{k}) }{\eps |U| - \la}

\section{Ten Band Model for Twisted Bilayer Graphene}
\label{app:Aaron}

In this section, we put forward a variation of the ten band model discussed in Ref.~\onlinecite{2018arXiv180802482P}. Our ten band model captures the flat band irreps of twisted bilayer graphene via an $S$-matrix construction consisting of $6$ orbitals in the $L$ sublattice and $4$ orbitals in the $L̃$ sublattice, yielding two flat bands. However, we compute the Wilson loop[Song} et~al., 2019, Alexandradinata et~al., 2014] of the flat bands in our model and find it has winding number $2$, instead of winding number $1$ as in the Bistritzer-MacDonald Hamiltonian [Song} et~al., 2019, Po} et~al., 2018, Bouhon} et~al., 2018]. In contrast, the model of \Ref{2018arXiv180802482P}, however, in the larger lattice has orbitals at both the $1a$ and $3f$ positions, and thus does not satisfy uniform pairing because the orbitals at the different Wyckoff positions are not related by symmetry.

To remedy this issue we propose a modified $S$-matrix which instead has in the $L$ lattice $s$ orbitals in the 6j position, which are related by $C_6z 𝒯$ and hence have uniform pairing.  The $L̃$ lattice possesses $p_x, p_y$ orbitals at the $2c$ positions (the hexagon).  In SSG $177.151$, the irreps induced by the two sublattices read
\begin{align}
E_{2c} \uparrow G &= 2 \Gamma_3 + (2K_1 \oplus K_2K_3) + (2M_1 \oplus 2M_2) \\
A_{6j} \uparrow G &= (\Gamma_1 \oplus \Gamma_2 \oplus 2\Gamma_3) + (2K_1 \oplus 2K_2K_3) + (3M_1 \oplus 3M_2) \\
{\cal B}_{\text{FB}} &= (\Gamma_1 \oplus \Gamma_2) + K_2K_3 + (M_1 \oplus M_2).
\label{}
\end{align}  The irreps of the flat band correctly match the irreps carried by the flat bands in twisted bilayer graphene, and thus these two bands must be topological.  Note that as this is a lattice model, one cannot correctly implement the particle-hole symmetry that protects the strong topology in TBG [Song} et~al., 2021].  The real-space Hamiltonian involves only NN hopping and NNN hopping between $s$ and $p$ orbitals, and there are only four complex independent parameters:

\begin{align}
H_\text{kin}^{\nu,\sigma} &= \sum_{\RR} t_{NN,x}^{\nu,\sigma} c^\dagger_{\RR,1,\sigma,\nu} d_{\RR,p_x,1,\sigma,\nu} + t_{NN,y}^{\nu,\sigma} c^\dagger_{\RR,1,\sigma,\nu} d_{\RR,p_y,1,\sigma,\nu} \nonumber \\
&+ t_{NNN,x}^{\nu,\sigma} c^\dagger_{\RR,1,\sigma,\nu} d_{\RR-\aaa_2,p_x,2,\sigma,\nu} + t_{NNN,y}^{\nu,\sigma} c^\dagger_{\RR,1,\sigma,\nu} d_{\RR-\aaa_2,p_y,2,\sigma,\nu} + (p6') + (H.c.)
\label{}
\end{align} where $σ$ denotes the spin and $ν$ denotes the valley, giving four identical copies of the kinetic Hamiltonian.  Here the $c$ orbitals belong to the $L$ sublattice, and are labeled $α= 1 …6$, labeling each of the $6j$ positions, while the $d$ orbitals belong to the $L̃$ sublattice and are labeled $β= {p_x, 1}, {p_y, 1}, {p_x, 2}, {p_y, 2}$, labeling the $2c$ positions.  In addition to the spin degree of freedom, the TBG model has a valley degree of freedom, $ν$, from the modes arising from the Dirac cones from the $K$ and $K'$ points in the BZ.  $p6'$ indicates that the appropriate symmetry-related terms are added to yield the full tight-binding Hamiltonians.  We find that the choice of parameters
\begin{align}
t_{NN,x}^{\nu=+1,\sigma=\u} &= 76.50+2.69i, t_{NN,y}^{\nu=+1,\sigma=\u} = -14.64+35.81i, \nonumber \\
t_{NNN,x}^{\nu=+1,\sigma=\u} &= 12.17-7.44i, t_{NNN,y}^{\nu=+1,\sigma=\u} = 5.47-5.45i
\label{}
\end{align} yields an excellent fit of the low-lying bands to the BM model, as depicted in Fig.~\ref{fig:TBGUPC}(b).  The lattice vectors are $_1 = (-√(3)/2,-1/2)a_0, _2 = (√(3)/2,-1/2)a_0$, where $a_0$ is the unit cell length.

\begin{figure*}
\includegraphics[width=.95\textwidth]{Figure_TBGwilson}
\caption{We summarize the 10-band UPC model of TBG. $(a)$ shows the electron orbitals of the model in the unit cell, $(b)$ shows the single-particle bands with the Wilson loop over the flat bands inset, $(c)$ shows the Cooper pair spectrum, and $(d)$ shows the mean-field superfluid weight at finite temperature.}
\label{fig:TBGUPC}
\end{figure*}

The main difference between this ten-band model and previous models arises from the addition of the valley degree of freedom.  The most natural Hubbard-type interaction would be one that includes all spin and valleys summed together in the density:

\begin{align}
H = \dfrac{|U|}{2} \sum_{\RR,\alpha} \left(\sum_{\nu,\sigma} {\bar n}_{\RR,\alpha,\sigma,\nu} \right)^2.
\label{}
\end{align}

However, we find that this interaction is not tractable with our eta-pairing formalism.  Instead, however, the slightly different interaction
\begin{align}
H = \frac{|U|}{2} \sum_{\RR,\alpha} ({\bar S}^z_{\RR, \alpha, \nu\sigma = +1})^2 + ({\bar S}^z_{\RR, \alpha, \nu\sigma = -1})^2,
\label{}
\end{align} where the spin-valley operators $S̅^z_, α, νσ$ are defined as
\begin{align}
{\bar S}^z_{\RR, \alpha, \nu\sigma = +1} &= {\bar n}_{\RR, \alpha, \nu = +1, \sigma = \uparrow} - {\bar n}_{\RR, \alpha, \nu = -1, \sigma = \downarrow} \\
{\bar S}^z_{\RR, \alpha, \nu\sigma = -1} &= {\bar n}_{\RR, \alpha, \nu = +1, \sigma = \downarrow} - {\bar n}_{\RR, \alpha, \nu = -1, \sigma = \uparrow},
\label{eq:TBG_interaction}
\end{align}  can be recast into an eta-pairing problem.  So long as $U(1)$ valley, $U(1)$ spin, and $T$ are preserved, there are two eta-pairing operators that commute with the interaction Hamiltonian
\begin{align}
\eta_{+1}^\dagger = \sum_{\kk,m} \gamma^\dagger_{\kk,m,\nu = +1, \sigma = \uparrow} \gamma^\dagger_{-\kk,m,\nu = -1, \sigma = \downarrow} \\
\eta_{-1}^\dagger = \sum_{\kk,m} \gamma^\dagger_{\kk,m,\nu = +1, \sigma = \downarrow} \gamma^\dagger_{-\kk,m,\nu = -1, \sigma = \uparrow}.
\label{}
\end{align} This spin-valley-polarized model splits into two eta-pairing models, each giving the Cooper pair spectrum shown in Fig.~\ref{figure1}(c).  The two Hamiltonians are labeled by spin-valley, $νσ= +1$ and $νσ= -1$, and are
\begin{align}
H^{\nu\sigma = +1} &= \sum_{\sigma\nu = +1} H^{\nu,\sigma}_\text{kin} + \frac{|U|}{2} \sum_{\RR,\alpha} ({\bar S}^z_{\RR, \alpha, \nu\sigma = +1})^2 \\
H^{\nu\sigma = -1} &= \sum_{\sigma\nu = -1} H^{\nu,\sigma}_\text{kin} + \frac{|U|}{2} \sum_{\RR,\alpha} ({\bar S}^z_{\RR, \alpha, \nu\sigma = -1})^2.
\label{}
\end{align}

We now deduce the irreducible representations of the Cooper pair bands.  Let us focus on the $νσ= +1$ sector.  As argued and rigorously proved in App.~\ref{app:Cooperpairspectrum}, the Cooper pairs themselves transform like $s$ orbitals induced from the larger sublattice $L$, which in this case is the 6j position.  We offer here an intuitive argument for why this must be so.  Consider an \emph{unprojected} Cooper pair operator at site $$ and flavor $α$, in the larger sublattice $L$:
\begin{align}
\eta_{\RR, \alpha}^\dagger = c^\dagger_{\RR,\alpha,\uparrow} c^\dagger_{\RR,\alpha,\downarrow}.
\label{}
\end{align}  Let us assume there is one orbital per position (at least in the larger sublattice $L$), so that the crystalline symmetry $g$ is a permutation matrix up to a phase: the orbital $α$ is located at position $x_α$, and under $g$ this is mapped to orbital $β$ located at position $x_β= g x_α$.  Under a crystalline symmetry $g$, this operator transforms as
\begin{align}
g \eta_{\RR, \alpha}^\dagger g^{-1} &= g c^\dagger_{\RR,\alpha,\uparrow} c^\dagger_{\RR,\alpha,\downarrow} g^{-1} \\
&= \sum_{\beta,\beta'} D[g]_{\alpha\beta} D^*[g]_{\alpha\beta'} c^\dagger_{g\RR,\beta,\uparrow} c^\dagger_{g\RR,\beta',\downarrow} \\
&= \eta_{g\RR, \beta}^\dagger,~~{\bf x}_\beta = g {\bf x}_\alpha.
\label{}
\end{align}  Any phases introduced in the permutation matrix $D$ is canceled by the complex conjugate from the time-reversed $D$ arising from the opposite spin.  Thus, these $η$ operators transform as $s$ orbitals.  When projecting into the flat bands, the orbitals belonging to the smaller sublattice $L̃$ possess no wavefunction weight, and so those bands do not participate in the Cooper pair Hamiltonian, leaving $s$ orbitals in the larger $L$ sublattice.

In our TBG model, these irreps read
\begin{align}
A_{6j} \uparrow G &= (\Gamma_1 \oplus \Gamma_2 \oplus 2\Gamma_3) + (2K_1 \oplus 2K_2K_3) + (3M_1 \oplus 3M_2).
\label{}
\end{align}
The $6j$ position is non-maximal, so the bands form a composite representation which can in principle be separated into gapped bands. The $Γ_1$ irrep is the zero-energy Cooper pair which has $s$-wave symmetry. The $Γ_2$ irrep is odd under $C_2x$ and is of $p_z$ orbital character. The $Γ_3$ irrep is formed from angular momentum $±1$ states under the $C_3z$ operator, and so is of $p_x,p_y$ orbital character (see \Tab{tab:32irreps} below).
\begin{center}
\label{tab:32irreps}
\begin{table}
\begin{tabular}{c | c c c }
& $1$ & $C_{3z}$ & $C_{2x}$ \\
\hline
$\Gamma_1$ & $1$ & $1$ & $1$ \\
$\Gamma_2$ & $1$ & $1$ & $-1$ \\
$\Gamma_3$ & $2$ & $-1$ & $0$ \\
\end{tabular}
\caption{Irreps of the group $32$ which is the unitary subgroup of $6'22'$, the little group of the Gamma point in twisted bilayer graphene.}
\end{table}
\end{center}

\section{Mean-field superfluid weight in degenerate flat bands}
\label{app:meanfield}

In this section, we derive the mean-field superfluid weight of a set of $N_f$ degenerate isolated flat bands with uniform pairing. We take $e,ħ,k_B=1$. In this section only, we use bra-ket notation for the single-particle eigenstates since we are working at the mean-field level. Our calculation is straightforward generalization of the $N_f = 1$ case demonstrated in \Ref{2022arXiv220311133H}. The result agrees with prior work (\Ref{2020PhRvL.124p7002X}) on degenerate flat bands, but correctly accounts for the minimal quantum metric in the self-consistent mean-field solution [Huhtinen} et~al., 2022].

We describe our system with the attractive mean-field approximated Hubbard model
\begin{align}
H &= \sum_{\mbf{k}} \crea{\vec{c}_{\mbf{k}}} H_{\rm BdG}(\mbf{k})
\ani{\vec{c}_{\mbf{k}}} + \sum_{\mbf{k}} {\rm
Tr }\, \tilde{h}^{\d}(\mbf{k}) - N_{orb} \mathcal{N}\mu -
\mathcal{N}\sum_{\alpha}\frac{|\Delta_{\alpha}|^2}{U}, \\
H_{\rm BdG}(\mbf{k}) &= \begin{pmatrix}
\tilde{h}^{\u}(\mbf{k}) - \mu \mathbb{1}_{N_{orb}} & \vec{\Delta}\\
\vec{\Delta}^{\dag} & -\tilde{h}^{\d}(\mbf{k})+\mu \mathbb{1}_{N_{orb}}
\end{pmatrix}, \label{eq.hbdg}
\end{align}
where $μ$ is the chemical potential, $h̃^σ$ is the
single-particle hopping matrix for spin $σ$ and
$c⃗_k =
c_-k,α=1,,…,c_-k,α=N_orb,)$. The
matrix $Δ⃗ = diag(Δ_1,…,Δ_N_orb)$
contains the mean-field order parameters $Δ_α =
Δ_iα =
Uc_Rαc_Rα$, where
$c_Rασ$ destroys a particle with spin $σ$ at
orbital $α$ in the unit cell $R$. We assume the order parameters
are independent of the unit cell. Since the we have shown in \App{app:excit} that the lowest energy gap is $s$-wave, zero momentum, and has the wavefunction $∝∑_R, c̅^†_R,,c̅^†_R,,$, we expect our mean-field ansatz to accurately capture this Cooper pair branch.

The superfluid weight for this Hamiltonian is given
by~[Huhtinen} et~al., 2022] (where $∂_i=∂/∂k_i$)
\begin{align}
[D_s]_{ij} = \frac{1}{\mathcal{N}V_c}
\sum_{\mbf{k},ab}\frac{n_F(E_a)-n_F(E_b)}{E_b-E_a}\big[
&\bra{\psi_a}\partial_{i}\widetilde{H}_{\mbf{k}} \ket{\psi_b}
\bra{\psi_b} \partial_{j}\widetilde{H}_{\mbf{k}} \ket{\psi_a} \nonumber \\
\ket{\psi_b} \bra{\psi_b}
\ket{\psi_a}
\big] - \frac{1}{V_c}C_{ij}, \label{eq.linrespresult}
\end{align}
\bea
\partial_{i}\widetilde{H_{\mbf{k}}} &= \begin{pmatrix}
\frac{\partial \tilde{h}^{\u}(\mbf{k}')}{\partial
k_{i}'}\bigg|_{\mbf{k}'=\mbf{k}} & 0 \\
0 & \frac{\partial \tilde{h}^{\d}(\mbf{k}')^*}{\partial
\end{pmatrix}, \
\delta_{i}\Delta = \begin{pmatrix}
0 & \frac{{\rm d}\vec{\Delta}}{{\rm d}q_{i}}\bigg|_{\mbf{q}=\vec{0}} \\
\frac{{\rm d}\vec{\Delta}^{\dag}}{{\rm d}q_{i}}\bigg|_{\mbf{q}=\vec{0}} & 0
\end{pmatrix}, \
C_{ij} = \frac{1}{U}\sum_{\alpha} \frac{{\rm
d}\Delta_{\alpha}}{{\rm d}q_{i}}\frac{{\rm
d}\Delta_{\alpha}^*}{{\rm d}q_{j}}\bigg|_{\mbf{q}=\vec{0}} +
{\rm H.c.}
\eea
The vector $q$ in the derivatives $dΔ_α/dq_i$ corresponds to the insertion of an electromagnetic field with $A∼q$, which changes the self-consistent order parameters as  $Δ_iα→Δ_iαe^2iq·(R+r_)$. The Fermi distribution is $n_F(E)=1/(1+e^βE)$, where $β=1/T$ is the inverse temperature. The eigenvalues
and eigenvectors of $H_BdG$ are $E_a$ and $|ψ_a⟩$
respectively, and $γ_z = σ_z⊗1_N_orb$,
where $σ_i$ are the Pauli matrices. The coefficient
$[n_F(E_a)-n_F(E_b)]/(E_b-E_a)$ should be understood as $-∂n_F(E)/∂E$ when $E_a=E_b$. The total volume of the system is
$𝒩V_c$ where $V_c$ is the volume of a unit cell.

We will focus on the superfluid weight for a set of
degenerate flat bands $ℬ$ with energy
$ϵ_m$ ($m̅$ denotes a flat band) in a system with time-reversal symmetry
($h̃^(k) = [h̃^(-k)]^*$)
and uniform pairing. Due to time-reversal symmetry, there always
exists a choice of intra-unit-cell orbital positions which guarantees
that the derivatives of the order parameters at $q=0⃗$
vanish~[Huhtinen} et~al., 2022]. With the uniform pairing condition, these positions are given by \Eq{eq:xalphafixed}. Let us assume that this is the
chosen set of positions, so that the superfluid weight of
Eq.~\ref{eq.linrespresult} simplifies to
\begin{equation}
[D_s]_{ij} = \frac{1}{V_c\mathcal{N}}
\sum_{\mbf{k},ab}\frac{n_F(E_a)-n_F(E_b)}{E_b-E_a}\big[
\bra{\psi_a}\partial_{i}\widetilde{H}_{\mbf{k}} \ket{\psi_b}
\bra{\psi_b} \partial_{j}\widetilde{H}_{\mbf{k}} \ket{\psi_a}
\ket{\psi_b} \bra{\psi_b}
\partial_{j}\widetilde{H}_{\mbf{k}}\gamma^z
\ket{\psi_a}
\big], \label{eq.simpl_lin_resp}
\end{equation}
This can be expressed in terms of the dispersion relations and Bloch
functions of the bands by expressing the eigenvectors of $H_BdG$
as $|ψ_a⟩ =
∑_m=1^N_orb(w_+,am|+⟩⊗|m_k⟩ +
where $|m_k⟩$ is the eigenvector of
$h̃^(k)$ with eigenvalue $ϵ_m,k$. We
denote the eigenvectors of $σ_z$ with eigenvalues $±1$ by
$|±⟩$. With these definitions, Eq.~\eqref{eq.simpl_lin_resp} can
be rewritten as
\begin{align}
[D_s]_{ij} &= \frac{1}{V_c\mathcal{N}}\sum_{\mbf{k}}\sum_{mn} C_{pq}^{mn}
[j_{i}(\mbf{k})]_{mn}[j_{j}(\mbf{k})]_{pq}, \\
C_{pq}^{mn} &= 4\sum_{ab} \frac{n_F(E_a)-n_F(E_b)}{E_b-E_a}
[j_{i}]_{mn} &=
\bra{m_{\sigma,\mbf{k}}}\partial_{i}\tilde{h}^{\sigma}(\mbf{k})\ket{n_{\sigma,\mbf{k}}}\nonumber
\\
&= \delta_{mn}\partial_{i}\epsilon_{m,\mbf{k}} +
(\epsilon_{m,\mbf{k}}-\epsilon_{n,\mbf{k}})\brakett{\partial_{i}m_{\mbf{k}}}{n_{\mbf{k}}} . \label{eq:currentmatrix}
\end{align}

In a system with uniform pairing, we can make the ansatz
$Δ = Δ1$,
and $H_BdG =
∑_k∑_m=1^N_orb[(ϵ_m-μ)^2σ^z+Δ^2σ^x]⊗|m⟩⟨m|$, where the momentum dependence of $_m(k)$ and $|m_k⟩$ is suppressed. The
eigenvalues of $H_BdG$ are then $±E_m = ±√((ϵ_m-μ)^2+Δ^2)$, and the corresponding
eigenfunctions are $|ψ_m^+⟩ = (u_m|+⟩+v_m|-⟩)⊗|m⟩$ and $|ψ_m^-⟩ = (-v_m|+⟩ + u_m|-⟩)⊗|m⟩$, where (again suppressing the momentum dependence)
\begin{equation}
\label{eq:vnok}
u_m = \frac{1}{\sqrt{2}}\sqrt{1 + \frac{\epsilon_{m}-\mu}{E_m}}, \:\:
= \frac{1}{\sqrt{2}}\sqrt{1-\frac{\epsilon_{m}-\mu}{E_n}}.
\end{equation}
As Eq.~\ref{eq:currentmatrix} shows, all contributions to the superfluid weight which involve
$[j_i^σ]_mn$ with $m,n∈ℬ$ vanish (because $_m(k)$=$_n(k)$=$_m$ for the degenerate flat bands) vanish. Thus the
superfluid weight reduces to the geometric contribution
\begin{align}
[D_s]_{ij} &=
\frac{2|\Delta|^2}{V_c\mathcal{N}}\sum_{\mbf{k}}\sum_{m\in\mathcal{B}}\sum_{n\notin\mathcal{B}}
\left[ \frac{{\rm tanh}(\beta E_{\overline{m}}/2)}{E_{\overline{m}}}
- \frac{{\rm tanh}(\beta E_n/2)}{E_n}\right]
\frac{\epsilon_n-\epsilon_{\overline{m}}}{\epsilon_n+\epsilon_{\overline{m}}
- 2\mu} \left(
\brakett{\partial_{i}m}{n}\brakett{n}{\partial_{j}m} + {\rm H.c.}
\right) \nonumber \\
&+ \frac{|\Delta|^2}{V_c\mathcal{N}}\sum_{\mbf{k}}\sum_{\substack{m\neq
n\\m,n\notin\mathcal{B}}} \left[ \frac{{\rm tanh}(\beta
E_m/2)}{E_m} - \frac{{\rm tanh}(\beta E_n/2)}{E_n} \right]
\frac{\epsilon_n-\epsilon_m}{\epsilon_n+\epsilon_m-2\mu}\left(
\brakett{\partial_{i}m}{n}\brakett{n}{\partial_{j}m} + {\rm H.c.}
\right).
\end{align}
Let $W_n$ be the band gap between the bands in $ℬ$ and the
$n$th band. We now assume that $W_n≫|ϵ_m-μ|$
and $W_n≫|Δ|$ for all $n∉ℬ$, and the superfluid
weight reduces to~[Liang et~al., 2017]
\begin{equation}
[D_s]_{ij} = \frac{4\Delta^2}{V_c} \frac{{\rm tanh}(\beta
\end{equation}
where $g_ij =
(1/𝒩)∑_k∑_m∈ℬ∑_n∉ℬRe∂_imnn∂_jm$ is the
integral of the quantum metric of the set of bands $ℬ$. Since we assumed
above that the derivatives of the order parameters are zero, this
is in fact the minimal quantum metric~[Huhtinen} et~al., 2022].

At $T=0$, the filling and the pairing gap are given by
\begin{align}
N_f\mathcal{N}\nu &= \sum_{\mbf{k}}v^2{\rm Tr} \, P_{\mbf{k}} \, =
N_f\mathcal{N}v^2 \\
\Delta_{\alpha} = \Delta &=
\frac{|U|}{\mathcal{N}}\sum_{\mbf{k}}u v
[P_{\mbf{k}}]_{\alpha\alpha} = \eps|U|\sqrt{\nu(1-\nu)},
\end{align}
with $=N_f/N_L$, where $N_L$ is the number of orbitals where the flat band
states have a nonzero weight. We have used that $v = v_m̅,
u = u_m̅$ are independent of $k$ and $m̅$ for a flat band (see Eq. \Eq{eq:vnok}) and that
$Tr  P = N_f$. Noticing that
$|Δ|/E_m=2u v=2√(ν(1-ν))$,
we get the zero-temperature superfluid weight
\begin{equation}
[D_s]_{ij} = 8\frac{|U|\eps}{V_{c}}\nu(1-\nu)g_{ij}.
\end{equation}

\end{document}