# Geometric Superfluid Weight of Composite Bands

Guodong Jiang  Department of Physics, University of Nevada, Reno, Reno NV
89502, USA  Department of Applied Physics, Aalto University School of Science,
FI-00076 Aalto, Finland  Yafis Barlas  Department of Physics, University of
Nevada, Reno, Reno NV 89502, USA

###### Abstract

The superfluid weight of an isolated flat band in multi-orbital
superconductors contains contributions from the band’s quantum metric and a
lattice geometric term that depends on the orbital positions in the lattice.
Since the superfluid weight is a measure of the superconductor’s energy
fluctuation, it is independent of the lattice geometry, leading to the minimal
quantum metric of a band [1]. Here, a perturbation approach is developed to
study the superfluid weight and its lattice geometric dependence for composite
bands. When all orbitals exhibit uniform pairing, the quantum geometric term
contains each band’s contribution and an inter-band contribution between every
pair of bands in the composite. Based on a band representation analysis, they
provide a topological lower bound for the superfluid weight of an isolated
composite of flat bands. Using this perturbation approach, an analytical
expression of the lattice geometric contribution is obtained. It is expressed
in terms of Bloch functions, providing a convenient formula to calculate the
superfluid weight for multi-orbital superconductors.

##  I Introduction

The discovery of superconductivity in twisted two-dimensional (2D) crystals
[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15] has brought focus to the
nature of superconductivity in flat or weakly dispersive bands [16, 17, 18,
19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 1, 32]. While the pairing
mechanism remains unresolved, the superconducting fluctuations in flat bands
are exotic enough to warrant a complete understanding [33, 34, 35, 36, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]. In particular, is the
question of superfluid weight (SW) in flat band superconductors, which serves
as an essential criteria to distinguish superconductors from metals [51, 52].
The SW in conventional superconductors, irrespective of the symmetry of order
parameter, is inversely proportional to the effective mass [53, 54], ruling
out flat-band superconductivity. However, in multi-orbital superconductors,
band geometry rescues superconductivity by providing a finite SW, which is
proportional to the quantum metric of the band [17, 18, 19, 20, 22, 21, 23,
24, 25, 26, 28, 29, 30]. This quantum geometric SW dominates in flat band
superconductors like twisted 2D crystals [23, 24, 55]. Moreover, geometric
effects can result in transitions to exotic superconducting states [56, 57,
58], which necessitates the need for further understanding the role of band
geometry in multi-orbital superconductors.

Up to now, most studies of the geometric SW have been restricted to a single
isolated band [17, 19, 20, 22, 21, 23, 24, 25, 26, 28, 29, 30]. This is valid
when the pairing potential is comparable with the bandwidth of the active band
near the Fermi level but much smaller than the bandgap to other remote bands.
However, many examples exist where multiple bands lie within or close to the
pairing interaction window, violating this single-band projection. In these
cases, the inter-band geometric effects [18, 27, 59] can be significant. In
addition, quantum geometric quantities like the Berry curvature and quantum
metric are usually large near the local k-points where the remote bands are
energetically close to the active band. Therefore, it is essential to study
the geometric SW when a composite of bands is within or near the interaction
range.

In this paper, we employ a band projection formalism to study the SW of
composite bands. This projection to the composite transforms the orbital basis
to band basis but projects to a few bands most relevant to the interaction
energy scale determined by the pairing Hamiltonian. We then develop a
perturbation approach to calculate the SW, which is similar to the
$\textbf{k}\cdot\textbf{p}$ method of calculating the effective mass tensor of
bands in solids. This approach provides an alternative derivation of the
general results in Ref. [17, 20] and identifies each term of the quantum
geometric SW from transitions between Bogoliubov bands. The quantum geometric
SW of composite bands contains each band’s contribution and an inter-band
contribution between every pair of bands in the composite. This separation of
SW emphasizes the role of inter-band quantum metric, which originates from the
derivative of inter-band pairing functions. Although the inter-band geometric
terms are negative semidefinite and reduce the SW, the total SW is always
positive semidefinite, implying the local stability of the zero center-of-mass
momentum (CMM) state in the uniform pairing channel.

Analyzing the reduction of superfluid weight due to the inter-band pairing, we
derive a topological lower bound for the SW of an isolated composite of flat
bands (ICFB), Eq. 19. The necessity of such a lower bound arises from a very
fundamental question, as in realistic materials, the scenario of a single
isolated flat band is rare. Instead, many flat or weakly dispersive bands can
coexist in the interaction window. In the simplest case, one can ask if there
is a topological lower bound for the SW when two flat bands of opposite Chern
numbers $\pm 1$ are both within the pairing interaction window.

An analysis of this case from the perspective of band representation (BR)
turns out to be illuminating. If the two bands are constructed from the same
set of Wannier orbital basis of an elementary band representation (EBR) [60,
61], their topological lower bounds tend to cancel, $|+1+(-1)|=0$. In
contrast, when they belong to different EBRs, their topological lower bounds
will add, $|+1|+|-1|=2$. The latter is due to a suppression of the reduction
effect by inter-band pairings. This hints that a larger SW can be achieved if
a multi-orbital superconductor is doped to a set of bands that contain
multiple incomplete EBRs.

In the second half of this paper, we study the relation between lattice
geometry and SW for multi-orbital superconductors. Recently, it was shown that
an additional contribution encoding the modification of order parameters with
orbital positions in a unit cell must be included in the calculation of SW
[31, 1]. This lattice geometric term, combined with the quantum geometric
term, makes the total SW invariant with the choice of orbital positions [1],
therefore lattice geometry-independent [62]. Using the perturbation method, we
derive an analytical expression of the lattice geometric term in terms of
Bloch functions, Eq. 35. This expression provides a convenient way of
computing the lattice geometric term numerically without solving gap
equations.

The rest of this paper is organized as follows. In Sec. II, we develop the
projection procedure for composite bands and present the quantum geometric SW
calculated from the perturbation method. Sec. III contains an analysis of the
inter-band effect on the geometric SW and establishes a topological lower
bound for an ICFB. Sec. IV studies the derivative of gap equation with respect
to CMM using the perturbation method, thereby obtaining the second derivative
of grand potential with order parameters and analyzing its properties. In Sec.
V, we show the geometry independence of the SW and provide a convenient
formula for computing the lattice geometric term. Finally, we conclude our
findings in Sec. VI.

##  II Superfluid weight of Composite bands

We start from the tight-binding lattice Hamiltonian with intra-orbital
attractive interactions ($U_{\alpha}>0$):

$\displaystyle\hat{\mathcal{H}}=\sum_{ij,\alpha\beta,\sigma}t^{\sigma}_{ij,\alpha\beta}c^{\dagger}_{i\alpha\sigma}c_{j\beta\sigma}-\sum_{i,\alpha}U_{\alpha}\hat{n}_{i\alpha\uparrow}\hat{n}_{i\alpha\downarrow},$
(1)

where $i,j$ label the unit cell, $\alpha,\beta=1,...,s$ label the orbitals in
a unit cell, and $\sigma=\uparrow,\downarrow$ labels the spin. The tight-
binding matrix elements, $t^{\sigma}_{ij,\alpha\beta}$ encode the hopping
pattern and amplitudes from orbital $(j,\beta)$ to $(i,\alpha)$. This hopping
graph determines the band dispersion and topology. The second term is the
onsite density-density interaction, with
$\hat{n}_{i\alpha\sigma}=c^{\dagger}_{i\alpha\sigma}c_{i\alpha\sigma}$ the
number operator. We assume time-reversal symmetry (TRS) with zero spin-orbit
coupling, thereby restricting to singlet pairing.

For attractive pairing interactions, we assume a BCS-type ground state and
perform the usual mean-field decoupling of Eq. 1 to arrive at the mean-field
Hamiltonian $\hat{\mathcal{H}}_{MF}$. Throughout this paper, we focus on an
orbital-independent pairing matrix $\hat{\Delta}=\Delta\hat{\mathcal{I}}_{s}$,
referred to as the uniform pairing condition (UPC) in literature [17, 18, 20].
Here $\hat{\mathcal{I}}_{s}$ is the $s$-dim identity matrix in orbital space.
This choice of pairing channel can always be made for any intra-orbital
interactions obeying TRS by tuning parameters $U_{\alpha}$. It corresponds to
a particular solution channel to the BCS gap equation, which always leads to a
positive semidefinite SW tensor.

The superconducting fluctuations to finite CMM can be described within the
same framework by a reduced mean-field Hamiltonian (see App. A),

$\displaystyle\hat{\mathcal{H}}_{MF}(\textbf{q})=$
$\displaystyle\sum_{\textbf{k}}\textbf{C}_{\textbf{k},\textbf{q}}^{\dagger}H_{BdG,\textbf{k}}(\textbf{q})\textbf{C}_{\textbf{k},\textbf{q}}+N\sum_{\alpha=1}^{s}\frac{|\Delta_{\textbf{q},\alpha}|^{2}}{U_{\alpha}}$
$\displaystyle+\sum_{\textbf{k}}\text{Tr}\\{h_{-\textbf{k}+\textbf{q}}^{\downarrow}-\mu_{\textbf{q}}\\},$
(2)

where
$\textbf{C}_{\textbf{k},\textbf{q}}=(c_{\textbf{k}+\textbf{q},\alpha\uparrow},c_{-\textbf{k}+\textbf{q},\alpha\downarrow}^{\dagger})^{T}$
denotes the Nambu spinor, with q half of the CMM, $N$ denotes the number of
unit cells and

$\displaystyle
H_{BdG,\textbf{k}}(\textbf{q})=\begin{pmatrix}h^{\uparrow}_{\textbf{k}+\textbf{q}}-\mu_{\textbf{q}}&\hat{\Delta}_{\textbf{q}}\\\
\hat{\Delta}^{\dagger}_{\textbf{q}}&-(h^{\downarrow
T}_{-\textbf{k}+\textbf{q}}-\mu_{\textbf{q}})\end{pmatrix}$ (3)

is the Bogoliubov-de Gennes (BdG) matrix. TRS imposes that $h^{\downarrow
T}_{-\textbf{k}}=h^{\uparrow}_{\textbf{k}}$. Chemical potential
$\mu_{\textbf{q}}$ and pairing matrix $\hat{\Delta}_{\textbf{q}}$ are attained
self-consistently. Pairing matrix
$\hat{\Delta}_{\textbf{q}}=\text{diag}\\{\Delta_{\textbf{q},1},...,\Delta_{\textbf{q},s}\\}$
in general has its entries complex and q-dependent.

To compute the SW, we adopt the grand potential method [17]. SW is the second
total derivative of the superconducting free energy $F(\textbf{q})$ with
respect to q,

$\displaystyle
D_{s,\mu\nu}=\frac{1}{A}\frac{d^{2}F}{dq_{\mu}dq_{\nu}}\bigg{|}_{\textbf{q}=0}.$
(4)

For simplicity, we set the area $A=1$. The free energy
$F(\textbf{q})=\Omega(\textbf{q},\mu_{\textbf{q}},\hat{\Delta}_{\textbf{q}},\hat{\Delta}_{\textbf{q}}^{*})+\mu_{\textbf{q}}N_{e}$,
where the grand potential

$\displaystyle\Omega(\textbf{q},\mu_{\textbf{q}},\hat{\Delta}_{\textbf{q}},\hat{\Delta}_{\textbf{q}}^{*})=-\frac{1}{\beta}\ln\text{Tr}\\{e^{-\beta\hat{\mathcal{H}}_{MF}(\textbf{q})}\\}$
(5)

is a function of
$\textbf{q},\mu_{\textbf{q}},\hat{\Delta}_{\textbf{q}},\hat{\Delta}_{\textbf{q}}^{*}$
explicitly. $N_{e}$ is the fixed average number of electrons. In the presence
of TRS, $D_{s,\mu\nu}$ can be related to partial derivatives of the grand
potential [17, 31, 1]:

$\displaystyle D_{s,\mu\nu}=\bigg{[}\frac{\partial^{2}\Omega}{\partial
q_{\mu}\partial
q_{\nu}}-\frac{\partial\Delta_{\textbf{q},\alpha}^{I}}{\partial
q_{\mu}}\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{I}\partial\Delta_{\textbf{q},\beta}^{I}}\frac{\partial\Delta_{\textbf{q},\beta}^{I}}{\partial
q_{\nu}}\bigg{]}\bigg{|}_{\textbf{q}=0},$ (6)

where $\Delta_{\textbf{q},\alpha}^{I}$ is the imaginary part of
$\Delta_{\textbf{q},\alpha}$, and we have employed the Einstein summation rule
for $\alpha,\beta=1,...,s$. $U(1)$ symmetry allows one to fix
$\Delta_{\textbf{q},1}$ to be real and positve for all q, so the sum over
$\alpha,\beta$ above can be reduced as from $2$ to $s$ [1]. However, we assume
an arbitrary global phase to keep the discussion general.

The first term of Eq. 6 ($D_{s}^{(1)}$) contains two contributions to the SW
of multi-orbital superconductors—the conventional contribution, which is
inversely proportional to the effective mass, and the quantum geometric
contribution. This quantum geometric SW is proportional to an integral of the
quantum metric [17], which is lattice geometry-dependent [62]. However, as
expected from general considerations, this geometry-dependence is canceled by
the second term of Eq. 6 (the lattice geometric term, $D_{s}^{(2)}$), which
encodes the order parameter’s variation with orbital positions [31, 1, 32].

We consider the scenario that $n\leqslant s$ bands are energetically within or
close to the interaction window [63] around the Fermi level (the case of $n=2$
is shown in Fig. 1 (a)). We refer to the group of $n$ flat or dispersive
bands, whether connected or disconnected from each other, as a “composite”,
and denote it by $\mathcal{C}$. We perform a projection of the full
Hamiltonian Eq. II to the composite, a generalization of the single-band
projection in an earlier paper [58]. The projected mean-field Hamiltonian at
finite q, involving the kinetic energy of the $n$ bands and pairing
interactions between them, reads

$\displaystyle\hat{\mathcal{H}}_{\mathcal{C}}(\textbf{q})=$
$\displaystyle\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}[\xi^{\uparrow}_{l,\textbf{k}+\textbf{q}}c^{\dagger}_{l,\textbf{k}+\textbf{q}\uparrow}c_{l,\textbf{k}+\textbf{q}\uparrow}-\xi^{\downarrow}_{l,-\textbf{k}+\textbf{q}}c_{l,-\textbf{k}+\textbf{q}\downarrow}c^{\dagger}_{l,-\textbf{k}+\textbf{q}\downarrow}]+\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}\xi^{\downarrow}_{l,-\textbf{k}+\textbf{q}}$
$\displaystyle+\sum_{\textbf{k}}\sum_{l,l^{\prime}\in\mathcal{C}}[\Delta_{ll^{\prime},\textbf{k}}(\textbf{q})c^{\dagger}_{l,\textbf{k}+\textbf{q}\uparrow}c^{\dagger}_{l^{\prime},-\textbf{k}+\textbf{q}\downarrow}+h.c.]+N\sum_{\alpha=1}^{s}\frac{|\Delta_{\textbf{q},\alpha}|^{2}}{U_{\alpha}},$
(7)

where
$\xi^{\uparrow}_{l,\textbf{k}+\textbf{q}}=\varepsilon^{\uparrow}_{l,\textbf{k}+\textbf{q}}-\mu_{\textbf{q}}$,
$\xi^{\downarrow}_{l,-\textbf{k}+\textbf{q}}=\varepsilon^{\downarrow}_{l,-\textbf{k}+\textbf{q}}-\mu_{\textbf{q}}$
are the band energies measured from the fluctuated chemical potential. By TRS
we have defined
$\xi_{l\textbf{k}}\equiv\xi^{\uparrow}_{l\textbf{k}}=\xi^{\downarrow}_{l,-\textbf{k}}$.
Projection to the composite naturally defines the inter-band gap functions

$\displaystyle\Delta_{ll^{\prime},\textbf{k}}(\textbf{q})\equiv\langle
u_{l,\textbf{k}+\textbf{q}}|\hat{\Delta}_{\textbf{q}}|u_{l^{\prime},\textbf{k}-\textbf{q}}\rangle,$
(8)

where $l,l^{\prime}$ are band indices restricted in the composite and
$u_{l\textbf{k}}$ is the periodic part of Bloch function of the $l^{th}$
spin-$\uparrow$ band. For the case $l=l^{\prime}$,
$\Delta_{ll^{\prime},\textbf{k}}(\textbf{q})$ is just the intra-band gap
function $\Delta_{l,\textbf{k}}(\textbf{q})$ [58]. Under the UPC, the inter-
band gap functions ($l\neq l^{\prime}$) vanish at $\textbf{q}=0$ by the
orthogonality of Bloch functions but are nonvanishing at $\textbf{q}\neq 0$.

We developed a stationary state perturbation method to diagonalize the
quadratic Hamiltonian of composite bands, Eq. II. The perturbation method has
been commonly used to simplify Kubo formula results [64] but is mainly
restricted to the single-band case without considering band geometry. Here, it
not only reproduces the first term of Eq. 6, $D_{s}^{(1)}$, which agrees with
the general results in Ref. [17, 20], but gives each contribution the meaning
of transitions between Bogoliubov bands.

Figure 1: Interpretation of the band projection, composite, and geometric
terms of Eq. II from Bogoliubov band transitions (see App. D for details). (a)
Two electronic bands (1,2) energetically close to the interaction scale are
projected out from a four-orbital model, which defines the composite. The blue
rectangle indicates the interaction window of $2\Delta$ around the Fermi
level. (b)-(c), zoom-in picture of the composite. (b) Electronic bands mix
with hole bands to form four Bogoliubov bands ($1,\pm$ and $2,\pm$). At
$\textbf{q}=0$, these bands are particle-hole symmetric (PHS); at
$\textbf{q}\neq 0$, the dispersion is deformed and PHS is broken. SW measures
the second derivative of energy with respect to q. (c) For small q, the
deformation can be treated with perturbation, which enables transitions
between Bogoliubov bands. Here, only the transitions from band $1,-$ are
indicated by arrows. The $O(q^{2})$ transition amplitude to band $1,-$ itself,
combined with the $O(q)$ amplitude to band $1,+$ gives the intra-band
geometric SW $\propto g^{1}_{\mu\nu}$; the $O(q)$ transition amplitudes to
band $2,\pm$ give the inter-band terms $\propto
p^{(\pm)}_{12}g^{12}_{\mu\nu}$, respectively.

The first term of Eq. 6 is (see App. D for details)

$\displaystyle D^{(1)}_{s,\mu\nu}=$
$\displaystyle\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}\frac{\Delta^{2}}{E_{l\textbf{k}}^{3}}\bigg{(}\tanh\frac{\beta
E_{l\textbf{k}}}{2}-\frac{\beta E_{l\textbf{k}}}{2}\text{sech}^{2}\frac{\beta
E_{l\textbf{k}}}{2}\bigg{)}\partial_{\mu}\xi_{l\textbf{k}}\partial_{\nu}\xi_{l\textbf{k}}$
$\displaystyle+\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}\tanh\frac{\beta
E_{l\textbf{k}}}{2}\bigg{\\{}\frac{4\Delta^{2}}{E_{l\textbf{k}}}g^{l}_{\mu\nu}(\textbf{k})+\sum_{l^{\prime}\in\mathcal{C},l^{\prime}\neq
l}8\Delta^{2}\bigg{[}\frac{p^{(+)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}+\frac{p^{(-)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\bigg{]}g^{ll^{\prime}}_{\mu\nu}(\textbf{k})\bigg{\\}},$
(9)

where $\Delta$ is the uniform pairing order parameter at finite temperature
$T=1/\beta$, and $E_{l\textbf{k}}=\sqrt{\xi_{l\textbf{k}}^{2}+\Delta^{2}}$ is
the $\textbf{q}=0$ quasiparticle energy. The first line of Eq. II corresponds
to the conventional contribution to SW, and the two separate pieces in the
second line correspond to the intra- and inter-band geometric contributions.
The different transition processes represented by these geometric
contributions can be interpreted from stationary state perturbations and are
illustrated in Fig. 1 (see App. D for details). We have introduced functions

$\displaystyle
p^{(\pm)}_{ll^{\prime}}(\textbf{k})\equiv\frac{1}{2}\bigg{(}1\pm\frac{\xi_{l\textbf{k}}\xi_{l^{\prime}\textbf{k}}+\Delta^{2}}{E_{l\textbf{k}}E_{l^{\prime}\textbf{k}}}\bigg{)},\,\,\,l\neq
l^{\prime},$ (10)

which can be termed as the $\textbf{q}=0$ inter-band coherence factors [65],
in contrast to the intra-band coherence factors ($l=l^{\prime}$) in
conventional BCS theory [53, 64]. $g^{l}_{\mu\nu}=\text{Re}R_{\mu\nu}^{l}$,
$g^{ll^{\prime}}_{\mu\nu}=\text{Re}R_{\mu\nu}^{ll^{\prime}}$ are the real part
of the intra- and inter-band component of the quantum geometric tensor
$R_{\mu\nu}^{l},R_{\mu\nu}^{ll^{\prime}}$, respectively, with

$\displaystyle
R_{\mu\nu}^{l}(\textbf{k})=\langle\partial_{\mu}u_{l\textbf{k}}|(1-|u_{l\textbf{k}}\rangle\langle
u_{l\textbf{k}}|)|\partial_{\nu}u_{l\textbf{k}}\rangle,$ $\displaystyle
R^{ll^{\prime}}_{\mu\nu}(\textbf{k})=-\langle\partial_{\mu}u_{l\textbf{k}}|u_{l^{\prime}\textbf{k}}\rangle\langle
u_{l^{\prime}\textbf{k}}|\partial_{\nu}u_{l\textbf{k}}\rangle,\,\,\,l\neq
l^{\prime}.$ (11)

They together constitute the quantum geometric tensor of the composite bands,

$\displaystyle
R^{\mathcal{C}}_{\mu\nu}(\textbf{k})=\sum_{l\in\mathcal{C}}R^{l}_{\mu\nu}(\textbf{k})+\sum_{l,l^{\prime}\in\mathcal{C},l\neq
l^{\prime}}R^{ll^{\prime}}_{\mu\nu}(\textbf{k}).$ (12)

where
$R^{\mathcal{C}}_{\mu\nu}(\textbf{k})=\text{Tr}\\{\hat{P}_{\mathcal{C}}(\textbf{k})\partial_{\mu}\hat{P}_{\mathcal{C}}(\textbf{k})\partial_{\nu}\hat{P}_{\mathcal{C}}(\textbf{k})\\}$
and
$\hat{P}_{\mathcal{C}}(\textbf{k})=\sum_{l\in\mathcal{C}}|u_{l\textbf{k}}\rangle\langle
u_{l\textbf{k}}|$ is the projection operator.

Tensor $g^{ll^{\prime}}_{\mu\nu}(\textbf{k})$ measures small distances between
Bloch functions of different bands in the momentum space,

$\displaystyle|\langle
u_{l\textbf{k}}|u_{l^{\prime},\textbf{k}+\textbf{q}}\rangle|^{2}\approx-
q_{\mu}q_{\nu}g^{ll^{\prime}}_{\mu\nu}(\textbf{k}).$ (13)

With our convention Eq. II, $g^{ll^{\prime}}_{\mu\nu}(\textbf{k})$ is negative
semidefinite, indicating that the inter-band geometric contributions in Eq. II
are negative and reduce the total SW. However, in App. E we show that the
total SW, which is the sum of the conventional, intra- and inter-band
geometric terms, is always positive semidefinite. This implies that the
$\textbf{q}=0$ BCS state locally minimizes the free energy at any temperature
whenever the UPC is valid for the composite bands.

In addition to the first term $D^{(1)}_{s}$, the perturbation method also
gives an analytical expression of the Hessian matrix
$\partial^{2}\Omega/\partial\Delta_{\textbf{q},\alpha}^{I}\partial\Delta_{\textbf{q},\beta}^{I}|_{\textbf{q}=0}$,
thereby put the lattice geometric term $D_{s}^{(2)}$ into a computationally
convenient form. We address this in Sec. IV and V.

##  III Topological Lower Bound for an Isolated Composite of Flat Bands

The inter-band geometric terms in Eq. II are negative semidefinite and reduce
the SW. In this section, we use this reduction effect along with the relative
topology of the bands in the composite to establish a topological lower bound
for the SW of composite bands. For a single isolated flat band, the absolute
value of the Chern number provides a weak lower bound for the geometric SW
[17]. These bounds stem from the relation

$\displaystyle\text{Tr}g^{m}({\textbf{k}})\geqslant|B^{m}({\textbf{k}})|,$
(14)

for an isolated band $m$, with $B^{m}$ the Berry curvature of the band, and
the fact that the gauge-invariant part of the Wannier localization functional
(WLF) is related to the trace of quantum metric. Additional lower bounds due
to Wilson loop winding numbers or real space invariants in other topological
classifications have also been established [19, 24, 66]. One might ask, for
composite bands whether there is a lower bound for the SW, and what are the
possible generalizations of Eq. 14?

To answer these questions, we restrict our discussion to the zero temperature
limit in the following. At $T=0$, the problem is simplified and Eq. II reads

$\displaystyle D^{(1)}_{s,\mu\nu}=$
$\displaystyle\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}\bigg{[}\frac{\Delta^{2}}{E_{l\textbf{k}}^{3}}\partial_{\mu}\xi_{l\textbf{k}}\partial_{\nu}\xi_{l\textbf{k}}+\frac{4\Delta^{2}}{E_{l\textbf{k}}}g_{\mu\nu}^{l}(\textbf{k})\bigg{]}$
$\displaystyle+\sum_{\textbf{k}}\sum_{l,l^{\prime}\in\mathcal{C},l^{\prime}>l}\frac{16\Delta^{2}}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}p^{(+)}_{ll^{\prime}}(\textbf{k})g^{ll^{\prime}}_{\mu\nu}(\textbf{k}).$
(15)

Here, coherence factors $p^{(-)}_{ll^{\prime}}$ are canceled because they
describe Bogoliubov transitions between two quasiparticle bands or two
quasihole bands [Fig. 1(c)], which are suppressed at $T=0$, a “Pauli blocking”
effect. Following a similar procedure described in App. E, we break the intra-
band quantum metric $g^{l}$ into inter-band ones, then the geometric part of
Eq. III splits into two terms,

$\displaystyle D_{s,\mu\nu}^{(1)\text{geo}}=$
$\displaystyle-\sum_{\textbf{k}}\sum_{l\in\mathcal{C},l^{\prime}\notin\mathcal{C}}\frac{4\Delta^{2}}{E_{l\textbf{k}}}g^{ll^{\prime}}_{\mu\nu}(\textbf{k})$
(16)
$\displaystyle-\sum_{\textbf{k}}\sum_{l,l^{\prime}\in\mathcal{C},l^{\prime}>l}\frac{4\Delta^{2}(\xi_{l\textbf{k}}-\xi_{l^{\prime}\textbf{k}})^{2}}{E_{l\textbf{k}}E_{l^{\prime}\textbf{k}}(E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}})}g^{ll^{\prime}}_{\mu\nu}(\textbf{k}).$

In this expression, the first term contains the inter-band quantum metric
between bands inside and outside the composite, while the second term contains
these between bands inside the composite only.

To establish a topological lower bound, we consider a special class of
composite bands, where each band in the composite is within the interaction
window, and these bands are energetically close to each other—we call this an
“isolated composite of flat bands”( ICFB). Hereafter, our definition of flat
bands encompasses weakly dispersive bands. This requires that the bandwidth of
each band,
$\delta\xi_{l}=\max\\{|\xi_{l\textbf{k}}-\xi_{l\textbf{k}^{\prime}}|\\}$, and
the band gaps between them,
$\delta\xi_{ll^{\prime}}=\max\\{|\xi_{l\textbf{k}}-\xi_{l^{\prime}\textbf{k}^{\prime}}|\\}$
($l,l^{\prime}\in\mathcal{C}$) are all much smaller than the uniform order
parameter $\Delta$. This notion of ICFB gives the separation in Eq. 16 a well-
defined physical meaning.

In Eq. 16, the first term depends on the overall topology of the composite
since it contains the quantum metric between bands inside and outside the
composite, whereas the second term depends on the relative topology between
bands inside the composite. For an ICFB, as long as the first term is nonzero,
it is of order $\Delta$ and dominates over the second term, which is of order
$\delta\xi^{2}/\Delta$.

Furthermore, we assume that the ICFB contains an incomplete set of $n$
disconnected bands ($n<s$), each with Chern number $C_{l}$ ($1\leqslant
l\leqslant n$). When a few bands are connected, the Chern numbers will be
generalized to the integrals of nonabelian Berry curvature over the Brillouin
zone [67]. The total Chern number of an ICFB, $C=\sum_{l=1}^{n}C_{l}$, is
always an integer.

With these clarifications, it is important to note that the inequality Eq. 14
holds for composite bands also [17, 68, 24, 69, 59], which reads

$\displaystyle\text{Tr}g^{\mathcal{C}}({\textbf{k}})\geqslant|B^{\mathcal{C}}({\textbf{k}})|,$
(17)

where $g^{\mathcal{C}}$ and $B^{\mathcal{C}}$ are the quantum metric and Berry
curvature of the composite. Following the definition Eq. 12, we write the
quantum metric of composite bands as

$\displaystyle
g^{\mathcal{C}}_{\mu\nu}=\sum_{l\in\mathcal{C}}g^{l}_{\mu\nu}+\sum_{l,l^{\prime}\in\mathcal{C},l^{\prime}\neq
l}g^{ll^{\prime}}_{\mu\nu}=-\sum_{l\in\mathcal{C},l^{\prime}\notin\mathcal{C}}g^{ll^{\prime}}_{\mu\nu}.$
(18)

This is proportional to the leading term of Eq. 16, since
$\Delta^{2}/E_{l\textbf{k}}$ is approximately a constant for an ICFB. As a
result, the total SW is lowered bounded by $|C|$ (besides a proportional
constant), with $C$ the total Chern number of the composite. Exceptional cases
may exist, e.g., when the composite is complete ($n=s$), such that $C=0$ and
the first term of Eq. 16 vanishes. In these cases, the SW will rely on the
second term of Eq. 16 and the conventional term, which depends on the
dispersion details; then, a simple topological lower bound cannot be obtained.
We will see examples below when the composite has $C=0$ but still provides a
nonzero topological bound.

The result $D_{s}\geqslant|C|$ can be compared to the single isolated band
case. If we add the lower bound $|C_{l}|$ of each band, it would give
$D_{s}\geqslant\sum_{l=1}^{n}|C_{l}|$ for a composite. Therefore, we deduce
that by the reduction effect of the inter-band geometric SW, the lower bound
shrinks from the sum of absolute values to the absolute value of sum,
$|C|=|\sum_{l=1}^{n}C_{l}|$. This poses a contrast between composites with
$\\{C_{l}\\}$ of uniform signs and different signs, e.g., $+,+,...$ and
$+,-,...$, which can be further explained by analyzing the WLF [68, 17].

The absolute value of Chern number of the composite bands, $|C|$, is a lower
bound for the WLF of Wannier orbitals of the composite, providing a measure of
Wannier obstruction. If the composite is Wannier representable, then the bands
inside and outside the composite can be constructed from disjoint sets of
Wannier orbitals. In the language of “topological quantum chemistry” [60, 70,
61, 71], such a composite is said to form a BR. In this case, the inter-band
terms of Eq. III cancel the intra-band terms altogether, giving a zero lower
bound for the leading term of Eq. 16. On the contrary, in the presence of
Wannier obstruction, the composite has to form a BR with some bands outside
the composite. Then the inter-band terms of Eq. III partially cancel the
intra-band terms, leaving a finite lower bound for the leading term of Eq. 16.

Figure 2: An isolated composite of two weakly dispersive bands near the Fermi
level (only spin-up bands are shown), with the blue rectangle indicating the
interaction window of $2\Delta$. (a) The two bands in the composite have Chern
number $C_{1}=C_{2}=1$. The SW is lower bounded by $|C|=2$. (b) The two bands
have $C_{1}=1$ and $C_{2}=-1$, and the SW is lower bound by $|C|=0$. (c), same
as (b), but the two bands belong to different EBRs (shown in red and green),
then the lower bound can be improved to $|C_{1}|+|C_{2}|=2$.

However, the lower bound by $|C|$ in some cases is too weak, and a stronger
lower bound can be derived. Suppose the $s$ bands formed by the $s$ atomic
orbitals are divided into a few EBRs [60, 61], which are constructed from
disjoint sets of Wannier basis. For each EBR, it might be that only part of
the bands belong to the composite. By making an induction from BR to EBR, one
can immediately deduce that the inter-band geometric terms of Eq. III can only
reduce the SW substantially within an EBR, and not between two EBRs [72]. As a
result, the topological lower bound can be improved, with $|C|$ replaced by

$\displaystyle\bigg{|}\sum_{l_{1}=1}^{n_{1}}C^{(1)}_{l_{1}}\bigg{|}+...+\bigg{|}\sum_{l_{J}=1}^{n_{J}}C^{(J)}_{l_{J}}\bigg{|}=\sum_{j=1}^{J}\bigg{|}\sum_{l_{j}=1}^{n_{j}}C^{(j)}_{l_{j}}\bigg{|},$
(19)

where $J$ is the number of EBRs formed by the $s$ atomic orbitals, $n_{j}$ is
the number of bands in the composite that is from the $j^{th}$ EBR (with the
sum $\sum_{j=1}^{J}n_{j}=n$), and $C^{(j)}_{l_{j}}$ is the Chern number of
$l_{j}\,^{th}$ band of the $j^{th}$ EBR. In the extreme case when each band of
the composite belongs to a different EBR, there is no substantial reduction of
SW by the inter-band terms, so the improved lower bound is
$\sum_{l=1}^{n}|C_{l}|$. This simple idea of reducing topological lower bound
within individual EBRs is illustrated in Fig. 2 using the example of a two-
band composite. Similar arguments can be generalized to systems with nonzero
Chern numbers to attain stronger bounds based on real space invariants [66].

##  IV Superfluid Weight and Lattice Geometry

We now study the relationship between SW and lattice geometry using the
perturbation method. Here “lattice geometry” refers to the position of atomic
orbitals in a unit cell of non-Bravais lattices. It was shown in Ref. [1] that
when the atomic orbitals occupy some optimal positions, the second term of SW
$D_{s}^{(2)}$ vanishes and the geometric part of $D_{s}^{(1)}$ becomes the
minimal quantum metric (MQM). Typically, the optimal positions are high-
symmetry points in the real space lattice.

However, if the symmetry of the tight-binding Hamiltonian is lower than that
of the underlying Bravais lattice, the optimal positions can shift along high-
symmetry lines [1]. Since Bravais lattice symmetries no longer protect them,
the optimal positions can be sensitive to temperature, interaction, and the
tight-binding graph $t_{ij,\alpha\beta}^{\sigma}$. This observation suggests
that instead of seeking the optimal positions of MQM, it would be more
straightforward to work with some general positions and find an explicit
expression of the second term $D_{s}^{(2)}$. Notice in Eq. 6, one has to
either numerically compute the derivative of order parameters
$\partial_{\mu}\Delta^{I}_{\alpha}$ (a short-hand notation of
$\partial_{\mu}\Delta^{I}_{\textbf{q},\alpha}|_{\textbf{q}=0}$) or vary the
first term $D_{s}^{(1)}$ with respect to orbital positions to achieve MQM,
both of which may require a large number of calculations.

In this section, we derive an explicit expression of the Hessian matrix
$\partial^{2}\Omega/\partial\Delta_{\alpha}^{I}\partial\Delta_{\beta}^{I}$
(similarly, stands for
$\partial^{2}\Omega/\partial\Delta_{\textbf{q},\alpha}^{I}\partial\Delta_{\textbf{q},\beta}^{I}|_{\textbf{q}=0}$)
in Eq. 6 for general composite bands. It also contains an analysis of its rank
and semidefiniteness, which allows us to write $D_{s}^{(2)}$ into a convenient
form in Sec. V.

Before we proceed, it is essential to clarify and distinguish between
geometry-independent and dependent quantities [62]. Electronic band energies
and topological properties, such as Chern numbers, are determined by the
hopping integrals $t^{\sigma}_{ij,\alpha\beta}$, therefore, are geometry-
independent. Geometry-dependent quantities, like the quantum metric and Berry
curvature, instead depend on the orbital positions. Following Ref. [62], we
adopt the convention that the atomic orbitals are point-like, therefore the
periodic part of the Bloch function $u_{l\textbf{k}}$ depends on the orbital
positions through a componentwise phase,
$u_{l\textbf{k}\alpha}=e^{-i\textbf{k}\cdot\textbf{x}_{\alpha}}u^{0}_{l\textbf{k}\alpha}$.
Here $u^{0}_{l\textbf{k}\alpha}$ are the Bloch components of gauge choice at
some arbitrary reference positions in the unit cell, and $\textbf{x}_{\alpha}$
are the orbital positions measured from these reference points. This explains
the origin of geometric dependence of $D_{s}^{(1)}$, Eq. II.

According to the spirit of Ref. [1], the MQM can be defined in two ways, which
are precisely given below:
Definition 1: MQM corresponds to the orbital positions that make the second
term of Eq. 6, $D_{s}^{(2)}$ vanish.
Definition 2: MQM corresponds to the orbital positions that extremize (in
fact, minimize) some geometric functional of the quantum metric tensor, which
in general takes the form

$\displaystyle
I=\sum_{\textbf{k}}\big{[}\sum_{l\in\mathcal{C}}f_{l}(\textbf{k})\text{Tr}g^{l}(\textbf{k})+\sum_{l,l^{\prime}\in\mathcal{C},l\neq
l^{\prime}}f_{ll^{\prime}}(\textbf{k})\text{Tr}g^{ll^{\prime}}(\textbf{k})\big{]}.$
(20)

$I$ is a function of orbital positions $\textbf{x}_{1},...,\textbf{x}_{s}$
because the quantum metric
$g^{l}_{\mu\nu}(\textbf{k}),g^{ll^{\prime}}_{\mu\nu}(\textbf{k})$ depend on
them explicitly. $f_{l}(\textbf{k})$ and $f_{ll^{\prime}}(\textbf{k})$ are a
set of intra- and inter-band geometry-independent functions. For finite-
temperature SW calculations, they are

$\displaystyle f_{l}(\textbf{k})=\frac{1}{E_{l\textbf{k}}}\tanh\frac{\beta
E_{l\textbf{k}}}{2},$ (21) $\displaystyle
f_{ll^{\prime}}(\textbf{k})=\bigg{[}\frac{2p^{(+)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}+\frac{2p^{(-)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\bigg{]}\tanh\frac{\beta
E_{l\textbf{k}}}{2},$ (22)

such that the functional $I$ is proportional to the trace of the geometric
part of $D^{(1)}_{s}$. Notice that the quasiparticle energy and coherence
factors depend on the band energy only, therefore are geometry-independent.

The equivalence of the two definitions above relies on an assumption made in
Ref. [1] about the rank and semidefiniteness of matrix
$\partial^{2}\Omega/\partial\Delta_{\alpha}^{I}\partial\Delta_{\beta}^{I}$
(see App. H). Below, we prove this assumption by computing the matrix
explicitly.

It was argued in Ref. [1] that the order parameter minimizes the
superconducting free energy. Therefore, the matrix
$\partial^{2}\Omega/\partial\Delta^{I}_{\alpha}\partial\Delta^{I}_{\beta}$ has
to be positive semidefinite. However, this statement may not be true for a
general pairing matrix $\hat{\Delta}$ as a solution to the gap equation of
multi-orbital superconductors. As a solution channel, it guarantees that the
first derivative of grand potential $\Omega$ with respect to $\hat{\Delta}$ is
zero (see App. B), but in general cannot imply information about the second
derivative. The second derivative instead is related to the “stability” of the
channel.

Figure 3: Schematic plot of grand potential in the space of components of
$\hat{\Delta}$, at $\textbf{q}=0$. Multi-orbital superconductors may have
multiple solution channels, $A,B,C$, etc., to some fixed interaction
parameters. A solution channel $\hat{\Delta}$ satisfies the gap equation
$\partial\Omega/\partial\hat{\Delta}=0$, so extremizes $\Omega$, but the
minimization is not guaranteed. A ground state channel (point $A$) ensures the
global minimum, but there may also exist “unstable” solutions (e.g., point
$B$). This instability is irrelevant to convergence for numerically solving
the gap equation.

A multi-orbital superconductor can have multiple channels that are solutions
to some fixed interaction parameters, an example of which can be seen in the
supplementary material of Ref. [58]. If some channel corresponds to the
superconducting ground state, it should be the global minimum of free energy.
Therefore, its Hessian matrix needs to be positive semidefinite (Fig. 3).
However, for solutions that are not the ground state, the Hessian matrix may
not be positive semidefinite. For instance, in App. J, we show for two-orbital
models when the order parameter matrix takes the form
$\hat{\Delta}=\text{diag}\\{\Delta_{1},\Delta_{2}\\}$ with
$\Delta_{1},\Delta_{2}$ both real,
$\partial^{2}\Omega/\partial\Delta^{I}_{\alpha}\partial\Delta^{I}_{\beta}$ is
positive semidefinite if and only if $\text{Sgn}(\Delta_{1}\Delta_{2})=1$. For
this reason, obtaining an explicit expression of this matrix for the uniform
pairing channel is necessary. To achieve this goal, we study the gap equation
for composite bands at finite q and extract the Hessian matrix from its
derivative with respect to q.

Let’s define the finite-q mean-field order parameter as

$\displaystyle\Delta_{\textbf{q},\alpha}=-U_{\alpha}\langle
c_{i\alpha\downarrow}c_{i\alpha\uparrow}\rangle_{\textbf{q}}e^{-2i\textbf{q}\cdot(\textbf{R}_{i}+\textbf{x}_{\alpha})},$
(23)

where $\langle\rangle_{\textbf{q}}$ is the average over the fluctuated pairing
state $\Psi_{\textbf{q}}$, and the phase factor
$e^{-2i\textbf{q}\cdot(\textbf{R}_{i}+\textbf{x}_{\alpha})}$ gets rid of the
superficial q-dependence from $\Delta_{\textbf{q},\alpha}$. Taking the total
derivative of Eq. 23 with respect to $q_{\mu}$ at $\textbf{q}=0$, it yields a
matrix equation [31, 1]

$\displaystyle M_{\alpha\beta}\partial_{\mu}\Delta_{\beta}=V_{\alpha,\mu},$
(24)

where $\partial_{\mu}\Delta_{\alpha}=i\partial_{\mu}\Delta_{\alpha}^{I}$ by
TRS. After choosing a proportional constant, it is easy to show that (see App.
C)

$\displaystyle
M_{\alpha\beta}=\frac{1}{2}\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{I}\partial\Delta_{\textbf{q},\beta}^{I}}\bigg{|}_{\textbf{q}=0},$
(25) $\displaystyle
V_{\alpha,\mu}=-\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}\partial
q_{\mu}}\bigg{|}_{\textbf{q}=0},$ (26)

so the problem of computing matrix
$\partial^{2}\Omega/\partial\Delta^{I}_{\alpha}\partial\Delta^{I}_{\beta}$ is
reduced to determining the gap equation 23.

For example, the zero-temperature gap equation for a single isolated band $m$
can be expressed as [58]

$\displaystyle\Delta_{\textbf{q},\alpha}=$
$\displaystyle\frac{U_{\alpha}}{N}\sum_{\textbf{k}}u_{m,\textbf{k}-\textbf{q},\alpha}^{*}u_{m,\textbf{k}+\textbf{q},\alpha}$
$\displaystyle\times\frac{\Delta_{m,\textbf{k}}(\textbf{q})}{\sqrt{(\xi_{m,\textbf{k}+\textbf{q}}+\xi_{m,\textbf{k}-\textbf{q}})^{2}+4|\Delta_{m,\textbf{k}}(\textbf{q})|^{2}}},$
(27)

from which the derivative can be easily calculated. Here
$\Delta_{m,\textbf{k}}(\textbf{q})=\langle
u_{m,\textbf{k}+\textbf{q}}|\hat{\Delta}_{\textbf{q}}|u_{m,\textbf{k}-\textbf{q}}\rangle$
is the projected intra-band gap function. Unfortunately, for general composite
bands, a closed form like Eq. IV does not exist since one needs to diagonalize
the $2n$-dim mean-field Hamiltonian Eq. II to get the pairing state
$\Psi_{\textbf{q}}$. However, the perturbation method allows us to expand
$\Psi_{\textbf{q}}$ and the gap equation in powers of q. Thus, Eq. 24 as its
first derivative can be exactly calculated.

To proceed, we transform the gap equation 23 from orbital to band basis and
project to the composite, reading

$\displaystyle\Delta_{\textbf{q},\alpha}=-\frac{U_{\alpha}}{N}\sum_{\textbf{k}}\sum_{l,l^{\prime}\in\mathcal{C}}u_{l^{\prime},\textbf{k}-\textbf{q},\alpha}^{*}u_{l,\textbf{k}+\textbf{q},\alpha}\langle
c_{l^{\prime},-\textbf{k}+\textbf{q}\downarrow}c_{l,\textbf{k}+\textbf{q}\uparrow}\rangle_{\textbf{q}}.$
(28)

Then we use the standard perturbation method to calculate the pairing
amplitude $\langle
c_{l^{\prime},-\textbf{k}+\textbf{q}\downarrow}c_{l,\textbf{k}+\textbf{q}\uparrow}\rangle_{\textbf{q}}$
(see App. F). Taking derivative of Eq. 28 with $q_{\mu}$, we finally obtain
Eq. 24 with

$\displaystyle M_{\alpha\beta}=$
$\displaystyle\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}|u_{l\textbf{k}\alpha}|^{2}\delta_{\alpha\beta}\frac{\tanh(\beta
E_{l\textbf{k}}/2)}{2E_{l\textbf{k}}}-\frac{1}{2}\sum_{\textbf{k}}\sum_{l,l^{\prime}\in\mathcal{C}}u_{l^{\prime}\textbf{k}\alpha}^{*}u_{l\textbf{k}\alpha}u_{l\textbf{k}\beta}^{*}u_{l^{\prime}\textbf{k}\beta}$
$\displaystyle\times\bigg{[}\frac{p^{(+)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}\bigg{(}\tanh\frac{\beta
E_{l\textbf{k}}}{2}+\tanh\frac{\beta
E_{l^{\prime}\textbf{k}}}{2}\bigg{)}+\frac{p^{(-)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\bigg{(}\tanh\frac{\beta
E_{l\textbf{k}}}{2}-\tanh\frac{\beta
E_{l^{\prime}\textbf{k}}}{2}\bigg{)}\bigg{]}$ (29)

and

$\displaystyle V_{\alpha,\mu}=$
$\displaystyle\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}(u_{l\textbf{k}\alpha}^{*}\partial_{\mu}u_{l\textbf{k}\alpha}-\partial_{\mu}u_{l\textbf{k}\alpha}^{*}u_{l\textbf{k}\alpha})\frac{\Delta}{2E_{l\textbf{k}}}\tanh\frac{\beta
E_{l\textbf{k}}}{2}-\sum_{\textbf{k}}\sum_{l,l^{\prime}\in\mathcal{C}}u_{l^{\prime}\textbf{k}\alpha}^{*}u_{l\textbf{k}\alpha}\langle
u_{l\textbf{k}}|\partial_{\mu}u_{l^{\prime}\textbf{k}}\rangle$
$\displaystyle\times\bigg{[}\frac{\Delta
p^{(+)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}\bigg{(}\tanh\frac{\beta
E_{l\textbf{k}}}{2}+\tanh\frac{\beta
E_{l^{\prime}\textbf{k}}}{2}\bigg{)}+\frac{\Delta
p^{(-)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\bigg{(}\tanh\frac{\beta
E_{l\textbf{k}}}{2}-\tanh\frac{\beta
E_{l^{\prime}\textbf{k}}}{2}\bigg{)}\bigg{]}.$ (30)

Although matrix $M_{\alpha\beta}$ and vector $V_{\alpha,\mu}$ take a
complicated form, they have some simple properties we shall discuss now. We
assume that the composite $\mathcal{C}$ contains all the $s$ orbitals of the
lattice model; therefore, $M_{\alpha\beta}$ is an $s$-dim square matrix. This
assumption may not hold for some exceptional cases, e.g., the atomic limit of
the single-band projection of two-orbital models (App. I) and the isolated
flat band limit of Lieb lattice, which will be discussed in Sec. V.

From Eq. IV and IV above, one first notices that $M_{\alpha\beta}$ is a real
symmetric matrix, while $V_{\alpha,\mu}$ is a purely imaginary vector.
Moreover, $M_{\alpha\beta}$ contains factors $|u_{l\textbf{k}\alpha}|^{2}$,
$u_{l^{\prime}\textbf{k}\alpha}^{*}u_{l\textbf{k}\alpha}$ and quasiparticle
energy $E_{l\textbf{k}}$, which are all geometry-independent quantities,
therefore is geometry-independent, whereas $V_{\alpha,\mu}$ is geometry-
dependent. These properties are all governed by Eq. 24.

One can also notice that $M_{\alpha\beta}$ of Eq. IV is not invertible, since
$\sum_{\beta=1}^{s}M_{\alpha\beta}=0$ (similarly,
$\sum_{\alpha=1}^{s}V_{\alpha,\mu}=0$), which implies that it has an
eigenvector $\textbf{v}_{0}=(1,1,...,1)^{T}$ with eigenvalue zero. Physically,
this is associated with the $U(1)$ symmetry of the order parameter. The
existence of such a kernel vector implies that the rank of $M_{\alpha\beta}$
cannot exceed $s-1$. In general, one can prove that under the UPC if the
composite bands contain all the $s$ orbitals, then $\text{Rank}(M)=s-1$. We
prove this for two special cases in App. G.

Nevertheless, one can imagine a process in which some orbitals are removed
from the composite, which may arise from topological phase transitions or
changes in the interaction scale. Whenever an $r^{th}$ orbital is removed
($1\leqslant r\leqslant s$), $\textbf{e}_{r}=(...,0,1,0,...)^{T}$ (with the
$r^{th}$ component nonzero only) becomes an additional kernel vector of
$M_{\alpha\beta}$, lowering $\text{Rank}(M)$ by one. The reason is that the
physics will no longer depend on its position as the $r^{th}$ orbital becomes
irrelevant to the composite.

Besides, one can show that $M_{\alpha\beta}$ is positive semidefinite, which
means it has one zero eigenvalue and $s-1$ positive eigenvalues (see App. G).
This vital property, as a consequence of UPC, reveals that the uniform pairing
channel belongs to these “stable” solutions to the gap equation like point $A$
and $C$ in Fig. 3.

##  V Geometry Independence and the Lattice Geometric Term in Terms of Bloch
Functions

With the expression of matrix $M_{\alpha\beta}$ and vector $V_{\alpha,\mu}$,
it is not difficult to show that each component of the SW tensor is geometry-
independent. For this, we write the second term of Eq. 6 as

$\displaystyle
D^{(2)}_{s,\mu\nu}=2M_{\alpha\beta}\partial_{\mu}\Delta_{\alpha}\partial_{\nu}\Delta_{\beta},$
(31)

where $\partial_{\mu}\Delta_{\alpha}=i\partial_{\mu}\Delta_{\alpha}^{I}$.
Under a geometric transformation which translates orbital positions from
$\textbf{x}^{0}_{\alpha}$ to
$\textbf{x}_{\alpha}=\textbf{x}^{0}_{\alpha}+\delta\textbf{x}_{\alpha}$ [73],
it is changed by

$\displaystyle\delta D^{(2)}_{s,\mu\nu}=$ $\displaystyle
2M_{\alpha\beta}(\delta\partial_{\mu}\Delta_{\alpha}\partial_{\nu}\Delta_{\beta}^{0}+\partial_{\mu}\Delta_{\alpha}^{0}\delta\partial_{\nu}\Delta_{\beta}$
$\displaystyle+\delta\partial_{\mu}\Delta_{\alpha}\delta\partial_{\nu}\Delta_{\beta}),$
(32)

where $\partial_{\mu}\Delta_{\alpha}^{0}$ is $\partial_{\mu}\Delta_{\alpha}$
evaluated at the initial positions and
$\delta\partial_{\mu}\Delta_{\alpha}^{0}$ is the change due to such
transformation. Eq. V can be further simplified using Eq. 101 in App. H.
Similarly, the change of the quantum geometric term of SW, $\delta
D^{(1)}_{s}$, can be evaluated by inserting Eq. H into Eq. II. Then, one can
immediately check the geometry independence of the total SW

$\displaystyle\delta D_{s}=\delta D^{(1)}_{s}+\delta D^{(2)}_{s}=0,$ (33)

therefore, we have finished establishing the relation between SW and lattice
geometry.

As mentioned at the beginning of Sec. IV, a shortcoming of Eq. 6 to compute
the SW is that one needs to either solve the gap equation at a few q points to
get the derivative $\partial_{\mu}\Delta_{\alpha}$, or vary the orbital
positions in the whole geometric space to find the MQM positions. However,
since we have treated the gap equation perturbatively to get Eq. IV and IV, it
means the total SW can be computed using simple integrals of Bloch functions
and their derivatives without following these procedures.

In fact, as pointed out in Ref. [31, 1], once we are convinced that
$\text{Rank}(M)=s-1$, it would be convenient to partially project [74] the
matrix $M_{\alpha\beta}$, vector $V_{\alpha,\mu}$ and
$\partial_{\mu}\Delta_{\alpha}$ into an $(s-1)$-dim subspace, in which the
matrix $M$ becomes invertible, so one can write

$\displaystyle\widetilde{\partial_{\mu}\Delta_{\alpha}}=(\widetilde{M}^{-1})_{\alpha\beta}\widetilde{V}_{\beta,\mu},$
(34)

where $\widetilde{}$ denotes the projected quantities. This finally yields a
simple expression for the lattice geometric term:

$\displaystyle
D^{(2)}_{s,\mu\nu}=2(\widetilde{M}^{-1})_{\alpha\beta}\widetilde{V}_{\alpha,\mu}\widetilde{V}_{\beta,\nu}.$
(35)

In this expression, the projected $\widetilde{M}$ and $\widetilde{V}$ are not
unique but can be taken as eliminating the first row and column from matrix
$M$ and eliminating the first component from vector $V$ in Eq. IV and IV. This
corresponds to restricting to the geometric space spanned by the $s-1$
orbitals.

Figure 4: SW calculation of the Lieb lattice model. (a) The gray square
indicates the unit cell of the Lieb lattice and the three orbitals $A,B,C$.
Top: the tight-binding hoppings can be tuned to $(1\pm\delta)t$ to open a gap
$2\sqrt{2}t\delta$ between the flat and dispersive bands; Mid: lattice
geometry assuming the three orbitals sitting on the same site; Bottom:
assuming the orbitals sitting on the regular sites of Lieb lattice. (b)(c) SW
computed using Eq. II and 35 at $T=0$, $t/\Delta=10$, $\mu=0$. (b) assumes the
geometry of (a) mid while (c) assumes the geometry of (a) bottom.

To show the utility of Eq. 35, we take the Lieb lattice model [18, 1] as an
example. A unit cell of the Lieb lattice is shown at the top of Fig. 4 (a),
where parameters $t$ and $\delta$ completely determine the tight-binding
matrix elements. Remember that the hopping graph does not specify the
positions of orbital $A,B$, and $C$, which leaves the degree of freedom for
writing the Bloch Hamiltonian $h^{\sigma}_{\textbf{k}}$. Since the hopping
integrals are usually defined from a given lattice structure, in this case,
one can still imagine the underlying Bravais lattice as a square lattice with
$C_{4v}$ symmetry.

When $\delta\neq 0$, the hopping graph breaks $C_{4}$ symmetry, making the
tight-binding model acquire a lower symmetry than the imagined square lattice.
Then it was found that the optimal positions of MQM of atom $A$ and $C$ no
longer sit at the midpoint of the side but shift towards $B$ atom [1]. As a
contrast, in our calculations using Eq. II and 35, we can place the orbitals
at any general positions in the unit cell. In Fig. 4 (b), we imagine the three
orbitals sitting on the same site, while Fig. 4 (c) assumes they are sitting
on the regular sites of the Lieb lattice. Since we make no approximations and
take into account all the inter-band and lattice geometric effects, the total
SW data of the two geometries (two solid curves in (b) and (c)) completely
match, showing the geometry independence of the SW. One also notices that in
Fig. 4 (b) the $D_{s}^{(1)}$ curve converges to $D_{s}^{(1)}+D_{s}^{(2)}$ in
the large $\delta$ limit, whereas in (c) they converge in the $\delta=0$
limit. This exhibits how the actual MQM geometry approaches the two geometries
in the two limits.

As a last comment, we address a possible issue of defining the composite when
applying Eq. 35 for numerical calculations. For most parameter ranges of the
Lieb lattice, one can take all three bands as the composite, which contains
three orbitals, then matrix $\widetilde{M}$ and vector $\widetilde{V}$ are
2-dim. However, in the isolated flat band limit, $\mu\sim 0$ and $\Delta\ll
t\delta$, the composite effectively contains just the flat band, which
comprises only two orbitals. Then as $\text{Rank}(M)$ is lowered to 1, we need
to project $M$ to a 1-dim subspace to make $\widetilde{M}$ invertible. This
gives a 1-dim matrix $\widetilde{M}$ and vector $\widetilde{V}$, which means
for Eq. IV and IV we only need to evaluate the integral for orbital $A$ or
$C$.

##  VI Conclusion and Outlook

We developed a perturbation approach to treat the center-of-mass momentum
(phase) fluctuation of multi-orbital superconductors. This approach naturally
separates the superfluid weight of composite bands into intra- and inter-band
contributions from Bogoluibov band transitions. Based on the reduction by the
inter-band geometric terms, a topological lower bound was derived for the
superfluid weight of isolated composite flat bands. Applying the theory of
band representation, we found that the elementary band representations should
be considered the minimum “block” for deriving a stronger lower bound.

In studying the geometry independence of the superfluid weight, we emphasized
the role of the derivative of gap equations in getting the Hessian matrix
$\partial^{2}\Omega/\partial\Delta_{\alpha}^{I}\partial\Delta_{\beta}^{I}$.
The perturbation method facilitated the calculation as we expanded the gap
equation in powers of q. We analyzed the properties of this Hessian matrix in
detail and showed how the positive semidefiniteness and rank information stem
from the uniform pairing condition. We finally rewrote the superfluid weight
formula into a convenient form suitable for numerical calculations of any
general tight-binding model.

Since the superfluid weight is the stiffness of the Meissner effect, it
naturally encodes the information about fluctuation to a finite-momentum
pairing state within a given channel, but little information about competition
between different channels, or between superconducting and other correlated
states [75, 76]. However, as we pointed out, the Hessian matrix of grand
potential with respect to order parameters can be nonpositive definite for
multi-orbital superconductors, and contain information about the stability of
the channel. Whether this allows us to determine the ground state is left to
future studies.

###### Acknowledgements.

The authors would like to acknowledge Jonah Herzog-Arbeitman for discussions
on the EBR and the support of UNR/VPRI startup grant PG19012. G. J. has also
been funded by the Jane and Aatos Erkko Foundation and the Keele Foundation as
part of the SuperC collaboration.

##  Appendix A Interacting Hamiltonian with Center-of-Mass Momentum
Fluctuation

By transformation
$c_{i\alpha\sigma}=(1/\sqrt{N})\sum_{\textbf{k}}e^{i\textbf{k}\cdot(\textbf{R}_{i}+\textbf{x}_{\alpha})}c_{\textbf{k}\alpha\sigma}$
(we chose this gauge for $c_{\textbf{k}\alpha\sigma}$ if not otherwise
specified), the onsite intra-orbital interaction takes the form

$\displaystyle\hat{\mathcal{H}}_{int}$
$\displaystyle=-\sum_{i,\alpha}U_{\alpha}c^{\dagger}_{i\alpha\uparrow}c^{\dagger}_{i\alpha\downarrow}c_{i\alpha\downarrow}c_{i\alpha\uparrow}$
(36)
$\displaystyle=-\sum_{\textbf{k}\textbf{k}^{\prime}\textbf{q},\alpha}\frac{U_{\alpha}}{N}c^{\dagger}_{\textbf{k}+\textbf{q},\alpha\uparrow}c^{\dagger}_{-\textbf{k}+\textbf{q},\alpha\downarrow}c_{-\textbf{k}^{\prime}+\textbf{q},\alpha\downarrow}c_{\textbf{k}^{\prime}+\textbf{q},\alpha\uparrow}.$

The conventional BCS theory makes the ansatz that electron
$\textbf{k}\uparrow$ pairs with $-\textbf{k}\downarrow$, therefore the reduced
interaction $\hat{\mathcal{H}}_{int}^{(\text{red})}$ only keeps the
$\textbf{q}=0$ term from above.

We now consider a pairing state of fluctuated CMM $2\textbf{q}$ (i.e. phase
fluctuation), so the reduced interaction keeps a term of fixed and finite q,
reading

$\displaystyle\hat{\mathcal{H}}_{int}^{(\text{red})}(\textbf{q})=-\sum_{\textbf{k}\textbf{k}^{\prime},\alpha}\frac{U_{\alpha}}{N}c^{\dagger}_{\textbf{k}+\textbf{q},\alpha\uparrow}c^{\dagger}_{-\textbf{k}+\textbf{q},\alpha\downarrow}c_{-\textbf{k}^{\prime}+\textbf{q},\alpha\downarrow}c_{\textbf{k}^{\prime}+\textbf{q},\alpha\uparrow}.$
(37)

This form of interaction implies that in the fluctuated state electron
$\textbf{k}+\textbf{q},\uparrow$ pairs with
$-\textbf{k}+\textbf{q},\downarrow$. With the order parameter
$\Delta_{\textbf{q},\alpha}$ defined by Eq. 23, Fourier transform gives

$\displaystyle\Delta_{\textbf{q},\alpha}=-\frac{U_{\alpha}}{N}\sum_{\textbf{k}}\langle
c_{-\textbf{k}+\textbf{q},\alpha\downarrow}c_{\textbf{k}+\textbf{q},\alpha\uparrow}\rangle_{\textbf{q}}.$
(38)

To get Eq. 38, we have used the fact that in the fluctuated state
$\Psi_{\textbf{q}}$ the electron $\textbf{k}+\textbf{q}\uparrow$ pairs with
$-\textbf{k}+\textbf{q}\downarrow$. Using the transformation between orbital
and band basis

$\displaystyle
c_{\textbf{k}\alpha\uparrow}=\sum_{l}u_{l\textbf{k}\alpha}c_{l\textbf{k}\uparrow},$
$\displaystyle
c_{-\textbf{k},\alpha\downarrow}=\sum_{l}u_{l\textbf{k}\alpha}^{*}c_{l,-\textbf{k}\downarrow},$
(39)

one obtains Eq. 28 (the sum over $l$ has been restricted to the composite due
to projection).

With 38, mean-field decoupling yields

$\displaystyle\hat{\mathcal{H}}_{int}^{(\text{red})}(\textbf{q})\simeq\sum_{\textbf{k},\alpha}(\Delta_{\textbf{q},\alpha}c^{\dagger}_{\textbf{k}+\textbf{q},\alpha\uparrow}c^{\dagger}_{-\textbf{k}+\textbf{q},\alpha\downarrow}+h.c.)+N\sum_{\alpha=1}^{s}\frac{|\Delta_{\textbf{q},\alpha}|^{2}}{U_{\alpha}}.$
(40)

For the kinetic part, we shift the dummy momentum k by $\pm\textbf{q}$ for the
spin-$\uparrow$, $\downarrow$ sector, respectively,

$\displaystyle\hat{\mathcal{H}}_{kin}=$
$\displaystyle\sum_{\textbf{k},\alpha\beta}[c^{\dagger}_{\textbf{k}\alpha\uparrow}(h^{\uparrow}_{\textbf{k}})_{\alpha\beta}c_{\textbf{k}\beta\uparrow}+c^{\dagger}_{-\textbf{k},\alpha\downarrow}(h^{\downarrow}_{-\textbf{k}})_{\alpha\beta}c_{-\textbf{k},\beta\downarrow}]-\mu_{\textbf{q}}\sum_{\textbf{k},\alpha}(c^{\dagger}_{\textbf{k}\alpha\uparrow}c_{\textbf{k}\alpha\uparrow}+c^{\dagger}_{-\textbf{k},\alpha\downarrow}c_{-\textbf{k},\alpha\downarrow})$
(41) $\displaystyle=$
$\displaystyle\sum_{\textbf{k},\alpha\beta}\\{c^{\dagger}_{\textbf{k}+\textbf{q},\alpha\uparrow}[(h^{\uparrow}_{\textbf{k}+\textbf{q}})_{\alpha\beta}-\mu_{\textbf{q}}\delta_{\alpha\beta}]c_{\textbf{k}+\textbf{q},\beta\uparrow}-c_{-\textbf{k}+\textbf{q},\alpha\downarrow}[(h^{\downarrow,T}_{-\textbf{k}+\textbf{q}})_{\alpha\beta}-\mu_{\textbf{q}}\delta_{\alpha\beta}]c^{\dagger}_{-\textbf{k}+\textbf{q},\beta\downarrow}\\}+\sum_{\textbf{k}}\text{Tr}\\{h^{\downarrow}_{-\textbf{k}+\textbf{q}}-\mu_{\textbf{q}}\\}.$

In the last line above,
$\sum_{\textbf{k}}\text{Tr}h^{\downarrow}_{-\textbf{k}+\textbf{q}}=\sum_{l,\textbf{k}}\varepsilon^{\downarrow}_{l,-\textbf{k}+\textbf{q}}$
is q-independent because k is summed over the entire Brillouin zone, whereas
$\sum_{\textbf{k}}\text{Tr}\mu_{\textbf{q}}=\mu_{\textbf{q}}sN$ with $N$ the
number of unit cells. Eq. 40 and 41 together give Eq. II.

##  Appendix B Self-consistency Equations

For completeness, we derive the gap equation and electron number equation and
study their properties under TRS, which is similar to Ref. [17]. With the
grand potential Eq. 5, we first derive the gap equation

$\displaystyle\frac{\partial\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}}=0,\,\,\,\forall\,\textbf{q},\alpha.$
(42)

Using Eq. II,

$\displaystyle\frac{\partial\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}}$
$\displaystyle=\frac{\text{Tr}\\{e^{-\beta\hat{\mathcal{H}}_{MF}(\textbf{q})}\sum_{\textbf{k}}c_{-\textbf{k}+\textbf{q},\alpha\downarrow}c_{\textbf{k}+\textbf{q},\alpha\uparrow}\\}}{\text{Tr}\\{e^{-\beta\hat{\mathcal{H}}_{MF}(\textbf{q})}\\}}+\frac{N}{U_{\alpha}}\Delta_{\textbf{q},\alpha}$
$\displaystyle=\sum_{\textbf{k}}\langle
c_{-\textbf{k}+\textbf{q},\alpha\downarrow}c_{\textbf{k}+\textbf{q},\alpha\uparrow}\rangle_{\textbf{q}}+\frac{N}{U_{\alpha}}\Delta_{\textbf{q},\alpha}.$
(43)

This cancels by the definition of order parameter 38; therefore, 42 is proved.

Next, we look at the electron number equation:

$\displaystyle\frac{\partial\Omega}{\partial\mu_{\textbf{q}}}=-N_{e},\,\,\,\forall\,\textbf{q},$
(44)

where the fixed average number of electrons $N_{e}$ is a constraint for the
system; condition $\forall\,\textbf{q}$ means when we compare the fluctuated
state $\Psi_{\textbf{q}}$ of different q, $N_{e}$ is always fixed at the same
value. Using Eq. II again one obtains

$\displaystyle\frac{\partial\Omega}{\partial\mu_{\textbf{q}}}=-\frac{\text{Tr}\\{e^{-\beta\hat{\mathcal{H}}_{MF}(\textbf{q})}\sum_{\textbf{k}\alpha\sigma}c^{\dagger}_{\textbf{k}\alpha\sigma}c_{\textbf{k}\alpha\sigma}\\}}{\text{Tr}\\{e^{-\beta\hat{\mathcal{H}}_{MF}(\textbf{q})}\\}},$
(45)

therefore Eq. 44 is just

$\displaystyle N_{e}=\sum_{\textbf{k}\alpha\sigma}\langle
c^{\dagger}_{\textbf{k}\alpha\sigma}c_{\textbf{k}\alpha\sigma}\rangle_{\textbf{q}}.$
(46)

The two self-consistency equations 38, 46 determine all the properties of the
fluctuated state $\Psi_{\textbf{q}}$, and one can solve them for
$\hat{\Delta}_{\textbf{q}}$, $\mu_{\textbf{q}}$. In the presence of TRS, one
can show that whenever $\\{\hat{\Delta}_{\textbf{q}},\mu_{\textbf{q}}\\}$ is a
solution to the self-consistency equations of q state, then
$\\{\hat{\Delta}_{\textbf{q}}^{\dagger},\mu_{\textbf{q}}\\}$ is a solution to
the equations for $-\textbf{q}$ state (the single isolated band case was
proved in Ref. [58]). Then for the case of diagonal
$\hat{\Delta}_{\textbf{q}}$ matrix,
$\Delta_{-\textbf{q},\alpha}=\Delta_{\textbf{q},\alpha}^{*}$ and
$\mu_{-\textbf{q}}=\mu_{\textbf{q}}$, so there are the following identities:

$\displaystyle\frac{\partial\mu_{\textbf{q}}}{\partial
q_{\mu}}\bigg{|}_{\textbf{q}=0}=0,$ (47)
$\displaystyle\frac{\partial\Delta_{\textbf{q},\alpha}}{\partial
q_{\mu}}\bigg{|}_{\textbf{q}=0}$
$\displaystyle=-\frac{\partial\Delta_{\textbf{q},\alpha}^{*}}{\partial
q_{\mu}}\bigg{|}_{\textbf{q}=0}=i\frac{\partial\Delta_{\textbf{q},\alpha}^{I}}{\partial
q_{\mu}}\bigg{|}_{\textbf{q}=0},$ (48)

where $\Delta_{\textbf{q},\alpha}^{I}$ is the imaginary part of
$\Delta_{\textbf{q},\alpha}$. The importance of TRS will be discussed further
in App. C.

##  Appendix C Superfluid Weight Formula, Derivative of Gap Equation and Time-
reversal Symmetry

In this section, we derive the SW formula Eq. 6. We treat the free energy
$F(\textbf{q})$ as function of q only, while
$\Omega(\textbf{q},\mu_{\textbf{q}},\Delta_{\textbf{q},\alpha},\Delta_{\textbf{q},\alpha}^{*})$
is a function of
$\textbf{q},\mu_{\textbf{q}},\Delta_{\textbf{q},\alpha},\Delta_{\textbf{q},\alpha}^{*}$
explicitly.

Using
$F(\textbf{q})=\Omega(\textbf{q},\mu_{\textbf{q}},\Delta_{\textbf{q},\alpha},\Delta_{\textbf{q},\alpha}^{*})+\mu_{\textbf{q}}N_{e}$,
we have the first derivative

$\displaystyle\frac{dF}{dq_{\mu}}=$
$\displaystyle\frac{\partial\Omega}{\partial
q_{\mu}}+\bigg{(}\frac{\partial\Omega}{\partial\mu_{\textbf{q}}}+N_{e}\bigg{)}\frac{\partial\mu_{\textbf{q}}}{\partial
q_{\mu}}$
$\displaystyle+\sum_{\alpha}\bigg{(}\frac{\partial\Omega}{\partial\Delta_{\textbf{q},\alpha}}\frac{\partial\Delta_{\textbf{q},\alpha}}{\partial
q_{\mu}}+\frac{\partial\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}}\frac{\partial\Delta_{\textbf{q},\alpha}^{*}}{\partial
q_{\mu}}\bigg{)}.$ (49)

The second term cancels by the electron number equation, and the third term
cancels by the gap equation; therefore

$\displaystyle\frac{dF}{dq_{\mu}}=\frac{\partial\Omega}{\partial
q_{\mu}},\,\,\,\forall\,\textbf{q}.$ (50)

An important consequence of TRS is that $F(\textbf{q})$ is an even function of
q, so

$\displaystyle\frac{dF}{dq_{\mu}}\bigg{|}_{\textbf{q}=0}=\frac{\partial\Omega}{\partial
q_{\mu}}\bigg{|}_{\textbf{q}=0}=0.$ (51)

Caution is needed if one wants to prove this like Eq. B and 45, the derivative
should not act on $c_{\textbf{k}+\textbf{q},\uparrow}$ and
$c_{-\textbf{k}+\textbf{q},\downarrow}$ operators since they are to form the
Nambu basis. At the same time, the grand potential depends only on the
eigenvalues of the BdG Hamiltonian. One can calculate $\partial\Omega/\partial
q_{\mu}$ directly using an explicit form like Eq. D.1, and then show Eq. 51
holds as long as $\hat{\Delta}_{\textbf{q}=0}$ is a Hermitian matrix up to a
$U(1)$ phase, which TRS imposes. The case of a single isolated band has been
proved in Ref. [58]. If TRS is broken, SW being the second derivative of free
energy at $\textbf{q}=0$ will lose its meaning as a stiffness tensor, since
the first derivative of free energy is already nonzero.

In the presence of TRS, the second derivative of free energy is

$\displaystyle\frac{d^{2}F}{dq_{\mu}dq_{\nu}}=\frac{\partial^{2}\Omega}{\partial
q_{\mu}\partial q_{\nu}}+\frac{\partial^{2}\Omega}{\partial
q_{\mu}\partial\mu_{\textbf{q}}}\frac{\partial\mu_{\textbf{q}}}{\partial
q_{\nu}}$
$\displaystyle+\sum_{\alpha}\bigg{(}\frac{\partial^{2}\Omega}{\partial
q_{\mu}\partial\Delta_{\textbf{q},\alpha}}\frac{\partial\Delta_{\textbf{q},\alpha}}{\partial
q_{\nu}}+\frac{\partial^{2}\Omega}{\partial
q_{\mu}\partial\Delta_{\textbf{q},\alpha}^{*}}\frac{\partial\Delta_{\textbf{q},\alpha}^{*}}{\partial
q_{\nu}}\bigg{)}$ (52)

To simplify this, we take the total derivative of gap equation 42 with respect
to $q_{\mu}$, yielding

$\displaystyle\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}\partial
q_{\mu}}+\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}\partial\mu_{\textbf{q}}}\frac{\partial\mu_{\textbf{q}}}{\partial
q_{\mu}}$ (53)
$\displaystyle+\sum_{\beta}\bigg{(}\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}\partial\Delta_{\textbf{q},\beta}}\frac{\partial\Delta_{\textbf{q},\beta}}{\partial
q_{\mu}}+\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}\partial\Delta_{\textbf{q},\beta}^{*}}\frac{\partial\Delta_{\textbf{q},\beta}^{*}}{\partial
q_{\mu}}\bigg{)}=0.$

Changing variables
$\Delta_{\textbf{q},\alpha},\Delta_{\textbf{q},\alpha}^{*}\rightarrow\Delta_{\textbf{q},\alpha}^{R},\Delta_{\textbf{q},\alpha}^{I}$
($\Delta_{\textbf{q},\alpha}\equiv\Delta_{\textbf{q},\alpha}^{R}+i\Delta_{\textbf{q},\alpha}^{I}$),
evaluating at $\textbf{q}=0$ and using TRS properties Eq. 47, 48, one can
check $d_{\mu}d_{\nu}F(\textbf{q})|_{\textbf{q}=0}$ gives the SW formula Eq.
6.

Moreover, Eq. 53 at $\textbf{q}=0$ gives

$\displaystyle\frac{1}{2}\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{I}\partial\Delta_{\textbf{q},\beta}^{I}}\bigg{|}_{\textbf{q}=0}\partial_{\mu}\Delta_{\beta}+\frac{\partial^{2}\Omega}{\partial\Delta_{\textbf{q},\alpha}^{*}\partial
q_{\mu}}\bigg{|}_{\textbf{q}=0}=0,$ (54)

which exactly maps to
$M_{\alpha\beta}\partial_{\mu}\Delta_{\beta}-V_{\alpha,\mu}=0$ and the total
derivative of Eq. B. Therefore, we get Eq. 25, 26 in the main text.

##  Appendix D Perturbation Approach to Calculate the Superfluid Weight

In this appendix, we provide the steps of using the perturbation method to
calculate $D_{s}^{(1)}$, Eq. II.

###  D.1 Nondegenerate Perturbation

Starting from the mean-field Hamiltonian projected to the composite bands, Eq.
II, the first and third terms can be organized into the $2n$-dim BdG form

$\displaystyle
H_{\textbf{k}}(\textbf{q})=\begin{pmatrix}\xi_{1,\textbf{k}+\textbf{q}}&\Delta_{1,\textbf{k}}(\textbf{q})&&\Delta_{12,\textbf{k}}(\textbf{q})&..&&\Delta_{1n,\textbf{k}}(\textbf{q})\\\
\Delta_{1,\textbf{k}}(\textbf{q})^{*}&-\xi_{1,\textbf{k}-\textbf{q}}&\Delta_{21,\textbf{k}}(\textbf{q})^{*}&&..&\Delta_{n1,\textbf{k}}(\textbf{q})^{*}&\\\
&\Delta_{21,\textbf{k}}(\textbf{q})&\xi_{2,\textbf{k}+\textbf{q}}&\Delta_{2,\textbf{k}}(\textbf{q})&&&\\\
\Delta_{12,\textbf{k}}(\textbf{q})^{*}&&\Delta_{2,\textbf{k}}(\textbf{q})^{*}&-\xi_{2,\textbf{k}-\textbf{q}}&&&\\\
..&..&&&..&&\\\
&\Delta_{n1,\textbf{k}}(\textbf{q})&&&&\xi_{n,\textbf{k}+\textbf{q}}&\Delta_{n,\textbf{k}}(\textbf{q})\\\
\Delta_{1n,\textbf{k}}(\textbf{q})^{*}&&&&&\Delta_{n,\textbf{k}}(\textbf{q})^{*}&-\xi_{n,\textbf{k}-\textbf{q}}\end{pmatrix},$
(55)

where the empty entries are zeroes. To aid in bookkeeping, in
$H_{\textbf{k}}(\textbf{q})$ we have assigned the $(2l-1)^{th}$ label to
electrons, and the $2l^{th}$ label to holes of the $l^{th}$ band.

Let $\Lambda_{\textbf{k}}(\textbf{q})$ be some generalized Bogoliubov
transformation that exactly diagonalizes $H_{\textbf{k}}(\textbf{q})$,

$\displaystyle\Lambda_{\textbf{k}}(\textbf{q})^{\dagger}H_{\textbf{k}}(\textbf{q})\Lambda_{\textbf{k}}(\textbf{q})=\hat{E}_{\textbf{k}}(\textbf{q}).$
(56)

Typically $\Lambda_{\textbf{k}}(\textbf{q})$ can only be determined
numerically.

$\displaystyle\hat{E}_{\textbf{k}}(\textbf{q})=\text{diag}\\{E_{1\textbf{k}+}(\textbf{q}),E_{1\textbf{k}-}(\textbf{q}),...\\}$
(57)

is a diagonal matrix consisting eigenvalues of $H_{\textbf{k}}(\textbf{q})$
(hat $\,\,\,\hat{}$ is to distinguish from a number). The $2n$ eigenvalues
correspond to the $n$ fluctuated Bogoliubov quasiparticle bands ($+$) and $n$
quasihole bands ($-$). At $\textbf{q}=0$ the spectrum is particle-hole
symmetric (PHS),
$E_{l\textbf{k}\pm}(0)=\pm\sqrt{\xi_{l\textbf{k}}^{2}+\Delta^{2}}$ (we also
use symbol $E_{l\textbf{k}}$ for $\sqrt{\xi_{l\textbf{k}}^{2}+\Delta^{2}}$).

After the diagonalization, $\hat{\mathcal{H}}_{\mathcal{C}}(\textbf{q})$ of
Eq. II reads

$\displaystyle\hat{\mathcal{H}}_{\mathcal{C}}(\textbf{q})=$
$\displaystyle\sum_{\textbf{k}}\big{[}\widetilde{\bm{\gamma}}_{\textbf{k},\textbf{q}}^{\dagger}\hat{E}_{\textbf{k}}(\textbf{q})\widetilde{\bm{\gamma}}_{\textbf{k},\textbf{q}}+\sum_{l\in\mathcal{C}}\xi^{\downarrow}_{l,-\textbf{k}+\textbf{q}}\big{]}+N\sum_{\alpha=1}^{s}\frac{|\Delta_{\textbf{q},\alpha}|^{2}}{U_{\alpha}}$
(58)

where $\widetilde{\bm{\gamma}}_{\textbf{k},\textbf{q}}$ is the Bogoliubov band
spinor that exactly diagonalizes $H_{\textbf{k}}(\textbf{q})$. The projected
grand potential is

$\displaystyle\Omega_{\mathcal{C}}=$
$\displaystyle-\frac{1}{\beta}\ln\text{Tr}\\{e^{-\beta\hat{\mathcal{H}}_{\mathcal{C}}(\textbf{q})}\\}$
$\displaystyle=$
$\displaystyle-\frac{1}{\beta}\sum_{\textbf{k}}\text{Tr}\ln(1+e^{-\beta\hat{E}_{\textbf{k}}(\textbf{q})})+\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}\xi^{\downarrow}_{l,-\textbf{k}+\textbf{q}}$
$\displaystyle+N\sum_{\alpha=1}^{s}\frac{|\Delta_{\textbf{q},\alpha}|^{2}}{U_{\alpha}}$
(59)

In
$D_{s,\mu\nu}^{(1)}=\partial_{\mu}\partial_{\nu}\Omega_{\mathcal{C}}|_{\textbf{q}=0}$,
the derivative acts on explicit q only and does not act on $\mu_{\textbf{q}}$
or $\hat{\Delta}_{\textbf{q}}$, giving [17]:

$\displaystyle D^{(1)}_{s,\mu\nu}=$
$\displaystyle-\frac{1}{2}\sum_{\textbf{k}}\text{Tr}\bigg{\\{}\frac{\partial^{2}\hat{E}_{\textbf{k}}(\textbf{q})}{\partial
q_{\mu}\partial q_{\nu}}\tanh\frac{\beta\hat{E}_{\textbf{k}}(\textbf{q})}{2}$
$\displaystyle+\frac{\beta}{2}\frac{\partial\hat{E}_{\textbf{k}}(\textbf{q})}{\partial
q_{\mu}}\frac{\partial\hat{E}_{\textbf{k}}(\textbf{q})}{\partial
q_{\nu}}\text{sech}^{2}\frac{\beta\hat{E}_{\textbf{k}}(\textbf{q})}{2}\bigg{\\}}\bigg{|}_{\textbf{q}=0}.$
(60)

Eq. D.1 involves the first and second-derivative of
$E_{l\textbf{k}\pm}(\textbf{q})$. If we treat q as perturbation and find the
eigenvalues of $H_{\textbf{k}}(\textbf{q})$ up to the $q^{2}$ order, then
$D_{s}^{(1)}$ can be determined. To proceed, we approximate
$\Lambda_{\textbf{k}}(\textbf{q})$ as

$\displaystyle\Lambda_{\textbf{k}}(\textbf{q})\simeq
U_{\textbf{k}}=U_{1\textbf{k}}\oplus...\oplus U_{n\textbf{k}},$ (61)

where

$\displaystyle
U_{l\textbf{k}}=\begin{pmatrix}w_{l\textbf{k}}&-v_{l\textbf{k}}\\\
v_{l\textbf{k}}&w_{l\textbf{k}}\end{pmatrix}$ (62)

is the Bogoliubov transformation for the $l^{th}$ band at $\textbf{q}=0$, with

$\displaystyle
w_{l\textbf{k}}=\sqrt{\frac{1}{2}\bigg{(}1+\frac{\xi_{l\textbf{k}}}{E_{l\textbf{k}}}\bigg{)}},\,\,\,v_{l\textbf{k}}=\sqrt{\frac{1}{2}\bigg{(}1-\frac{\xi_{l\textbf{k}}}{E_{l\textbf{k}}}\bigg{)}}.$
(63)

Transformation $U_{\textbf{k}}$ makes the off-diagonal elements of
$H_{\textbf{k}}(\textbf{q})$ small therefore splits into two parts,

$\displaystyle
U_{\textbf{k}}^{\dagger}H_{\textbf{k}}(\textbf{q})U_{\textbf{k}}=\hat{E}_{\textbf{k}}(0)+H_{\textbf{k}}^{(1)}(\textbf{q}),$
(64)

where
$\hat{E}_{\textbf{k}}(0)=\text{diag}\\{E_{1\textbf{k}},-E_{1\textbf{k}},...,E_{n\textbf{k}},-E_{n\textbf{k}}\\}$
is the unperturbed PHS Bogoliubov spectrum, and
$H^{(1)}_{\textbf{k}}(\textbf{q})$ is the perturbing matrix.

This is equivalent to expressing $\hat{\mathcal{H}}_{\mathcal{C}}(\textbf{q})$
in the Bogoliubov band basis defined by
$|l\textbf{k}\pm\rangle\equiv\gamma^{\dagger}_{l\textbf{k}\pm}(\textbf{q})|0\rangle$,
where $\gamma_{l\textbf{k}\pm}$ denote the Bogoliubov operators of the
$l^{th}$ band. The 2$n$-component Bogoliubov band spinor
$\bm{\gamma}_{\textbf{k},\textbf{q}}=(\gamma_{1\textbf{k}+}(\textbf{q}),\,\,\gamma_{1\textbf{k}-}(\textbf{q}),\,\,...)^{T}$
is related to the electronic band operator through

$\displaystyle\bm{\gamma}_{\textbf{k},\textbf{q}}=U_{\textbf{k}}^{\dagger}(c_{1,\textbf{k}+\textbf{q}\uparrow},c_{1,-\textbf{k}+\textbf{q}\downarrow}^{\dagger},...)^{T}.$
(65)

As we see below, the separation Eq. 64 results in an excitation Hamiltonian

$\displaystyle\hat{\mathcal{H}}_{\mathcal{C}}^{(1)}(\textbf{q})=\sum_{\textbf{k}}\bm{\gamma}_{\textbf{k},\textbf{q}}^{\dagger}H_{\textbf{k}}^{(1)}(\textbf{q})\bm{\gamma}_{\textbf{k},\textbf{q}},$
(66)

which separately accounts for the SW contributions from intra- and inter-band
transitions.

We first assume that at any k-point, no two bands of the composite are
degenerate ($\xi_{l\textbf{k}}=\xi_{l^{\prime}\textbf{k}}$), whether
accidentally or enforced by symmetry; additionally, we assume no two bands
have opposite energies ($\xi_{l\textbf{k}}=-\xi_{l^{\prime}\textbf{k}}$).
These restrictions mean the $2n$ Bogoliubov bands at $\textbf{q}=0$ are
disconnected throughout the Brillouin zone, so nondegenerate perturbation can
be applied. We find the energy eigenvalues of $H_{\textbf{k}}(\textbf{q})$ up
to $O(q^{2})$ terms give $n$ quasiparticle bands and $n$ quasihole bands

$\displaystyle
E_{l\textbf{k}+}(\textbf{q})=E_{l\textbf{k}}+H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l-1}+\sum_{l^{\prime}\in\mathcal{C},l^{\prime}\neq
l}\frac{|H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}-1}|^{2}}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}+\sum_{l^{\prime}\in\mathcal{C}}\frac{|H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}}|^{2}}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}+O(q^{3}),$
$\displaystyle
E_{l\textbf{k}-}(\textbf{q})=-E_{l\textbf{k}}+H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l}+\sum_{l^{\prime}\in\mathcal{C}}\frac{|H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}-1}|^{2}}{-E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}+\sum_{l^{\prime}\in\mathcal{C},l^{\prime}\neq
l}\frac{|H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}}|^{2}}{-E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}+O(q^{3}).$
(67)

Here $H_{\textbf{k}}^{(1)}(\textbf{q})_{ij}$ denotes the $(i,j)^{th}$ entry of
matrix $H_{\textbf{k}}^{(1)}(\textbf{q})$. In terms of the Bogoluibov basis,
they are

$\displaystyle
H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}-1}\equiv\langle
l\textbf{k}+|H_{\textbf{k}}^{(1)}(\textbf{q})|l^{\prime}\textbf{k}+\rangle,$
$\displaystyle H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}}\equiv\langle
l\textbf{k}-|H_{\textbf{k}}^{(1)}(\textbf{q})|l^{\prime}\textbf{k}-\rangle,$
$\displaystyle
H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}}\equiv\langle
l\textbf{k}+|H_{\textbf{k}}^{(1)}(\textbf{q})|l^{\prime}\textbf{k}-\rangle.$
(68)

which can be explicitly computed from Eq. 55 and 61.

After some algebra, we find the derivatives of Eq. D.1 with respect to q
(which do not act on $\mu_{\textbf{q}}$ or $\hat{\Delta}_{\textbf{q}}$) give

$\displaystyle\partial_{\mu}E_{l\textbf{k}\pm}(\textbf{q})\big{|}_{\textbf{q}=0}=\partial_{\mu}\xi_{l\textbf{k}},$
(69) $\displaystyle\partial_{\mu}\partial_{\nu}E_{l\textbf{k}\pm}(\textbf{q})$
$\displaystyle\big{|}_{\textbf{q}=0}=\pm\bigg{\\{}\frac{\xi_{l\textbf{k}}}{E_{l\textbf{k}}}\partial_{\mu}\partial_{\nu}\xi_{l\textbf{k}}-\frac{4\Delta^{2}}{E_{l\textbf{k}}}g^{l}_{\mu\nu}(\textbf{k})$
$\displaystyle-\sum_{l^{\prime}\in\mathcal{C},l^{\prime}\neq l}8$
$\displaystyle\Delta^{2}\bigg{[}\frac{p^{(+)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}+\frac{p^{(-)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\bigg{]}g^{ll^{\prime}}_{\mu\nu}(\textbf{k})\bigg{\\}},$
(70)

where the coherence factors

$\displaystyle
p^{(+)}_{ll^{\prime}}(\textbf{k})=\frac{1}{2}\bigg{(}1+\frac{\xi_{l\textbf{k}}\xi_{l^{\prime}\textbf{k}}+\Delta^{2}}{E_{l\textbf{k}}E_{l^{\prime}\textbf{k}}}\bigg{)}=(w_{l\textbf{k}}w_{l^{\prime}\textbf{k}}+v_{l\textbf{k}}v_{l^{\prime}\textbf{k}})^{2},$
$\displaystyle
p^{(-)}_{ll^{\prime}}(\textbf{k})=\frac{1}{2}\bigg{(}1-\frac{\xi_{l\textbf{k}}\xi_{l^{\prime}\textbf{k}}+\Delta^{2}}{E_{l\textbf{k}}E_{l^{\prime}\textbf{k}}}\bigg{)}=(w_{l\textbf{k}}v_{l^{\prime}\textbf{k}}-v_{l\textbf{k}}w_{l^{\prime}\textbf{k}})^{2}$
(71)

come from the second derivative
$\partial_{\mu}\partial_{\nu}|H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}}|^{2}$,
$\partial_{\mu}\partial_{\nu}|H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}-1}|^{2}$
and
$\partial_{\mu}\partial_{\nu}|H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}-1}|^{2}$,
$\partial_{\mu}\partial_{\nu}|H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}}|^{2}$,
respectively. Physically, $p^{(+)}_{ll^{\prime}}(\textbf{k})$ accounts for the
transitions between a quasiparticle and quasihole band, while
$p^{(-)}_{ll^{\prime}}(\textbf{k})$ accounts for transitions between two
quasiparticle bands or between two quasihole bands.

Also, one notices that the intra- and inter-band quantum metric comes from the
expansion of intra- and inter-band gap functions. If the fluctuation of
$\hat{\Delta}_{\textbf{q}}$ is ignored, they are (see definition Eq. 8)

$\displaystyle\Delta_{l,\textbf{k}}(\textbf{q})=$
$\displaystyle\Delta\\{1-2\langle
u_{l\textbf{k}}|\partial_{\mu}u_{l\textbf{k}}\rangle q_{\mu}$
$\displaystyle-[\langle\partial_{\mu}u_{l\textbf{k}}|\partial_{\nu}u_{l\textbf{k}}\rangle+c.c.]q_{\mu}q_{\nu}\\}+O(q^{3}),$
(72) $\displaystyle\Delta_{ll^{\prime},\textbf{k}}(\textbf{q})=-2\Delta\langle
u_{l\textbf{k}}|\partial_{\mu}u_{l^{\prime}\textbf{k}}\rangle
q_{\mu}+O(q^{2}),\,\,\,l\neq l^{\prime}.$ (73)

Since all matrix elements $H_{\textbf{k}}^{(1)}(\textbf{q})_{ij}\sim O(q)$,
the intra-band quantum metric $g^{l}_{\mu\nu}$ in Eq. D.1 comes from the
$O(q^{2})$ term of $H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l-1}$ combined
with the $O(q)$ term of $H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l}$ or the
$O(q^{2})$ term of $H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l}$ combined with
the $O(q)$ term of $H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l-1}$. Whereas the
inter-band quantum metric $g^{ll^{\prime}}_{\mu\nu}$ only comes from the
$O(q)$ term of $H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}-1}$,
$H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}}$,
$H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}-1}$ and
$H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}}$ of $l\neq l^{\prime}$. The
transition processes they represent are illustrated in Fig. 1.

Inserting Eq. 69, D.1 into Eq. D.1 one gets Eq. II in the main text.

###  D.2 Degenerate Perturbation

To get Eq. II, we assumed that there are not any two Bogoliubov bands at
$\textbf{q}=0$ being degenerate at any k point, so nondegenerate perturbation
was used. Otherwise, if a few quasiparticle bands are degenerate at some
isolated k points or completely degenerate throughout the momentum space
(which requires the band energy
$\xi_{l\textbf{k}}=\pm\xi_{l^{\prime}\textbf{k}}$), then degenerate
perturbation must be applied.

However, using standard degenerate perturbation, one can prove that Eq. II
remains the same even in the presence of band degeneracy, which agrees with
the general Kubo formula result. In calculating current responses using the
Kubo formula, whenever two bands become degenerate, the only change is to make
the substitution

$\displaystyle\frac{n(E_{1})-n(E_{2})}{E_{1}-E_{2}}\rightarrow\frac{dn(E)}{dE}\bigg{|}_{E=E_{1}}.$
(74)

Similarly, for Eq. II, when $E_{l\textbf{k}}=E_{l^{\prime}\textbf{k}}$, one
only needs to replace

$\displaystyle\frac{\tanh(\beta E_{l\textbf{k}}/2)-\tanh(\beta
E_{l^{\prime}\textbf{k}}/2)}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\rightarrow\frac{\beta}{2}\text{sech}^{2}\frac{\beta
E_{l\textbf{k}}}{2}.$ (75)

Therefore, we conclude that band degeneracy poses no additional mathematical
difficulty for expressing the SW $D_{s}^{(1)}$ as Eq. II.

##  Appendix E Positive Semidefiniteness of the Superfluid Weight

The purpose of this appendix is to show that the first term $D_{s}^{(1)}$, Eq.
II is positive semidefinite. This proof will be independent of lattice
geometry. Since the total SW $D_{s}=D_{s}^{(1)}+D_{s}^{(2)}$ is geometry-
independent and $D_{s}^{(2)}$ can be made zero by choosing the MQM, it also
shows the total $D_{s}$ is positive semidefinite for the uniform pairing
channel, regardless of temperature and any composite bands.

One first notices that the function

$\displaystyle\tanh\frac{\beta E_{l\textbf{k}}}{2}-\frac{\beta
E_{l\textbf{k}}}{2}\text{sech}^{2}\frac{\beta E_{l\textbf{k}}}{2}>0,$ (76)

so the conventional part of Eq. II is positive semidefinite. To analyze the
geometric term, we break the intra-band quantum metric $g^{l}_{\mu\nu}$ into
inter-band ones $g^{ll^{\prime}}_{\mu\nu}$, and collect all the terms with the
same $g^{ll^{\prime}}_{\mu\nu}$ (including the symmetric term
$l\leftrightarrow l^{\prime}$), to arrive at the following split expression,

$\displaystyle D^{(1)\text{geo}}_{s,\mu\nu}=$
$\displaystyle-\sum_{\textbf{k}}\sum_{l\in\mathcal{C},l^{\prime}\notin\mathcal{C}}\frac{4\Delta^{2}}{E_{l\textbf{k}}}\tanh\frac{\beta
E_{l\textbf{k}}}{2}g^{ll^{\prime}}_{\mu\nu}(\textbf{k})$
$\displaystyle-\sum_{\textbf{k}}\sum_{l,l^{\prime}\in\mathcal{C},l^{\prime}>l}4\Delta^{2}\bigg{\\{}\tanh\frac{\beta
E_{l\textbf{k}}}{2}\bigg{[}\frac{1}{E_{l\textbf{k}}}-\frac{2p_{ll^{\prime}}^{(+)}(\textbf{k})}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}-\frac{2p_{ll^{\prime}}^{(-)}(\textbf{k})}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\bigg{]}+l\leftrightarrow
l^{\prime}\bigg{\\}}g^{ll^{\prime}}_{\mu\nu}(\textbf{k}).$ (77)

Since each $g^{ll^{\prime}}_{\mu\nu}$ ($l\neq l^{\prime}$) is negative
semidefinite, all we need to show is that each quantity

$\displaystyle
f(\xi_{l\textbf{k}},\xi_{l^{\prime}\textbf{k}})=\tanh\frac{\beta
E_{l\textbf{k}}}{2}\bigg{[}\frac{1}{E_{l\textbf{k}}}-\frac{2p_{ll^{\prime}}^{(+)}(\textbf{k})}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}-\frac{2p_{ll^{\prime}}^{(-)}(\textbf{k})}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\bigg{]}+l\leftrightarrow
l^{\prime}$ (78)

is positive. Using definition D.1, we find

$\displaystyle f(\xi_{l\textbf{k}},$
$\displaystyle\xi_{l^{\prime}\textbf{k}})=-\bigg{(}\frac{1}{E_{l\textbf{k}}}\tanh\frac{\beta
E_{l\textbf{k}}}{2}-\frac{1}{E_{l^{\prime}\textbf{k}}}\tanh\frac{\beta
E_{l^{\prime}\textbf{k}}}{2}\bigg{)}\frac{\xi_{l\textbf{k}}-\xi_{l^{\prime}\textbf{k}}}{\xi_{l\textbf{k}}+\xi_{l^{\prime}\textbf{k}}}.$
(79)

Note the factor
$\frac{\xi_{l\textbf{k}}-\xi_{l^{\prime}\textbf{k}}}{\xi_{l\textbf{k}}+\xi_{l^{\prime}\textbf{k}}}=\frac{(\xi_{l\textbf{k}}-\xi_{l^{\prime}\textbf{k}})^{2}}{(E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}})(E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}})}$,
and function $(1/x)\tanh x$ has negative slope at $x>0$ (same as Eq. 76),
therefore $f(\xi_{l\textbf{k}},\xi_{l^{\prime}\textbf{k}})>0$ and the positive
semidefiniteness is proved.

##  Appendix F Calculation of the Pairing Amplitude

In this appendix, we give details of using the perturbation method to
calculate the pairing amplitude $\langle
c_{l^{\prime},-\textbf{k}+\textbf{q}\downarrow}c_{l,\textbf{k}+\textbf{q}\uparrow}\rangle_{\textbf{q}}$.

Recall that under transformation $U_{\textbf{k}}$, the mean-field Hamiltonian
can be expressed as

$\displaystyle\hat{\mathcal{H}}_{\mathcal{C}}(\textbf{q})=\sum_{\textbf{k}}\bm{\gamma}^{\dagger}_{\textbf{k},\textbf{q}}[\hat{E}_{\textbf{k}}(0)+H_{\textbf{k}}^{(1)}(\textbf{q})]\bm{\gamma}_{\textbf{k},\textbf{q}},$
(80)

where the perturbation matrix $H_{\textbf{k}}^{(1)}$ mixes Bogoliubov bands of
different indices. Using the standard perturbation method, the annihilation
operator of the band-resolved quasiparticles that exactly diagonalize
$\hat{\mathcal{H}}_{\mathcal{C}}(\textbf{q})$, $\widetilde{\gamma}$, can be
expressed as a linear combination of $\gamma$ operators:

$\displaystyle\widetilde{\gamma}_{l\textbf{k}+}(\textbf{q})=\gamma_{l\textbf{k}+}(\textbf{q})+\sum_{l^{\prime}\in\mathcal{C},l^{\prime}\neq
l}\frac{H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}-1}}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\gamma_{l^{\prime}\textbf{k}+}(\textbf{q})+\sum_{l^{\prime}\in\mathcal{C}}\frac{H_{\textbf{k}}^{(1)}(\textbf{q})_{2l-1,2l^{\prime}}}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}\gamma_{l^{\prime}\textbf{k}-}(\textbf{q})+O(q^{2}),$
$\displaystyle\widetilde{\gamma}_{l\textbf{k}-}(\textbf{q})=\gamma_{l\textbf{k}-}(\textbf{q})+\sum_{l^{\prime}\in\mathcal{C}}\frac{H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}-1}}{-E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\gamma_{l^{\prime}\textbf{k}+}(\textbf{q})+\sum_{l^{\prime}\in\mathcal{C},l^{\prime}\neq
l}\frac{H_{\textbf{k}}^{(1)}(\textbf{q})_{2l,2l^{\prime}}}{-E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}\gamma_{l^{\prime}\textbf{k}-}(\textbf{q})+O(q^{2}).$
(81)

Here, to expend $H_{\textbf{k}}^{(1)}(\textbf{q})_{ij}$ in powers of q, we
must take into account the q-dependence of $\hat{\Delta}_{\textbf{q}}$,
therefore use

$\displaystyle\Delta_{ll^{\prime},\textbf{k}}(\textbf{q})=$
$\displaystyle\Delta\delta_{ll^{\prime}}+(-2\Delta\langle
u_{l\textbf{k}}|\partial_{\mu}u_{l^{\prime}\textbf{k}}\rangle$
$\displaystyle+\sum_{\alpha=1}^{s}u_{l\textbf{k}\alpha}^{*}u_{l^{\prime}\textbf{k}\alpha}\partial_{\mu}\Delta_{\alpha})q_{\mu}+O(q^{2}),$
(82)

instead of Eq. D.1, 73. The normalization of $\widetilde{\gamma}$ operators is
ignored as it leads to $O(q^{2})$ corrections. Eq. F can be viewed as a linear
transformation $\widetilde{\bm{\gamma}}=A\bm{\gamma}$ between the spinor
$\widetilde{\bm{\gamma}}$ and $\bm{\gamma}$, where $A_{ij}=\delta_{ij}+a_{ij}$
is a matrix close to identity, with $a_{ij}\sim O(q)$. The inverse of $A$ has
an approximate form $(A^{-1})_{ij}=\delta_{ij}-a_{ij}+O(q^{2})$, which enables
us to write operators $\gamma$ in terms of $\widetilde{\gamma}$. Since the
electron operators $c$ are related to $\gamma$ through the transformation
$U_{\textbf{k}}$, we can finally write $c$ in terms of $\widetilde{\gamma}$.

Since the $\widetilde{\gamma}$ operators exactly diagonalize the BdG
Hamiltonian of composite bands, we have

$\displaystyle\langle\widetilde{\gamma}_{l\textbf{k}\pm}(\textbf{q})^{\dagger}\widetilde{\gamma}_{l\textbf{k}\pm}(\textbf{q})\rangle_{\textbf{q}}\approx
n_{F}(\pm E_{l\textbf{k}}),$ (83)

to the lowest order of q, where $n_{F}$ is the Fermi-Dirac function. Finally,
one obtains the pairing amplitude

$\displaystyle\langle
c_{l^{\prime},-\textbf{k}+\textbf{q}\downarrow}c_{l,\textbf{k}+\textbf{q}\uparrow}\rangle_{\textbf{q}}=\begin{cases}\displaystyle-\frac{\tanh(\beta
E_{l\textbf{k}}/2)}{2E_{l\textbf{k}}}\bigg{[}\Delta+\big{(}-2\Delta\langle
u_{l\textbf{k}}|\partial_{\mu}u_{l\textbf{k}}\rangle+\sum_{\beta=1}^{s}|u_{l\textbf{k}\beta}|^{2}\partial_{\mu}\Delta_{\beta}\big{)}q_{\mu}\bigg{]}+O(q^{2}),&l=l^{\prime},\\\
\displaystyle\frac{1}{2}\bigg{[}\frac{p^{(+)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}+E_{l^{\prime}\textbf{k}}}\bigg{(}\tanh\frac{\beta
E_{l\textbf{k}}}{2}+\tanh\frac{\beta
E_{l^{\prime}\textbf{k}}}{2}\bigg{)}+\frac{p^{(-)}_{ll^{\prime}}(\textbf{k})}{E_{l\textbf{k}}-E_{l^{\prime}\textbf{k}}}\bigg{(}\tanh\frac{\beta
E_{l\textbf{k}}}{2}-\tanh\frac{\beta
E_{l^{\prime}\textbf{k}}}{2}\bigg{)}\bigg{]}\\\
\displaystyle\times\big{(}2\Delta\langle
u_{l\textbf{k}}|\partial_{\mu}u_{l^{\prime}\textbf{k}}\rangle-\sum_{\beta=1}^{s}u_{l\textbf{k},\beta}^{*}u_{l^{\prime}\textbf{k},\beta}\partial_{\mu}\Delta_{\beta}\big{)}q_{\mu}+O(q^{2}),&l\neq
l^{\prime}.\end{cases}$ (84)

From this to get Eq. IV and IV, we have also used the $\textbf{q}=0$ gap
equation

$\displaystyle
1=\frac{U_{\alpha}}{N}\sum_{\textbf{k}}\sum_{l\in\mathcal{C}}|u_{l\textbf{k}\alpha}|^{2}\frac{\tanh(\beta
E_{l\textbf{k}}/2)}{2E_{l\textbf{k}}},\,\,\,\forall\,\alpha.$ (85)

##  Appendix G Positive Semidefiniteness and Rank of the Hessian Matrix for
Uniform Pairing

This appendix serves to prove the rank and positive semidefiniteness of the
Hessian matrix $M_{\alpha\beta}$ in Eq. IV.

We consider a general class of $s$-dim real symmetric matrices
$M_{\alpha\beta}$, which have a kernel eigenvector
$\textbf{v}_{0}=(1,1,...,1)^{T}$. It implies
$\sum_{\beta=1}^{s}M_{\alpha\beta}=\sum_{\alpha=1}^{s}M_{\alpha\beta}=0$. With
this property, we write $M_{\alpha\beta}$ as

$\displaystyle
M_{\alpha\beta}=\begin{cases}\sum_{\gamma=1,\gamma\neq\alpha}^{s}C_{\alpha\gamma},&\alpha=\beta\\\
-C_{\alpha\beta},&\alpha\neq\beta\end{cases}$ (86)

where $C_{\alpha\beta}$ ($\alpha\neq\beta$) are just the off-diagonal elements
of $M$ with a minus sign. We give three propositions below.

###### Proposition G.1

If $C_{\alpha\beta}>0$ for all $1\leqslant\alpha,\beta\leqslant
s,\,\alpha\neq\beta$, then $M_{\alpha\beta}$ is positive semidefinite and
$\text{Rank}(M)=s-1$.

We use the criterion that all its leading principal minors are nonnegative to
show the positive semidefiniteness. Let’s denote the leading principal minor
of $M$ of order $k$ ($1\leqslant k\leqslant s$) by $D_{k}(M)$, we find

$\displaystyle D_{1}(M)=M_{11}=\sum_{\alpha\neq 1}^{s}C_{1\alpha}\geqslant 0,$
(87) $\displaystyle D_{2}(M)=\det\begin{pmatrix}M_{11}&M_{12}\\\
M_{21}&M_{22}\end{pmatrix}=\sum_{\alpha_{1}\neq 1,\alpha_{2}\neq
2}^{s}C_{1\alpha_{1}}C_{2\alpha_{2}}-C_{12}C_{21},$ (88)

and

$\displaystyle
D_{k}(M)=\sideset{}{{}^{\prime}}{\sum}_{\alpha_{1},...,\alpha_{k}}^{s}C_{1\alpha_{1}}C_{2\alpha_{2}}...C_{k\alpha_{k}}.$
(89)

Here, $\sum^{\prime}$ means the summation imposed by the single rule that
there is not any subset $\\{i_{1},...,i_{j}\\}\subseteq\\{1,...,k\\}$ such
that $\\{\alpha_{i_{1}},...,\alpha_{i_{j}}\\}$ is its permutation. This single
rule leads to the following sub-rules: (i) no $C_{\alpha\alpha}$ appears in
the sum (e.g., 87), since $\\{\alpha\\}$ is a permutation of $\\{\alpha\\}$;
(ii) no squares of any $C$ appears in the sum, i.e. for $\alpha\neq\beta$,
$C_{\alpha\beta}C_{\beta\alpha}$ is not allowed since $\\{\beta,\alpha\\}$ is
a permutation of $\\{\alpha,\beta\\}$ (e.g., 88); (iii) in every term of
$D_{k}(M)$, there is at least an $\alpha_{j}$ ($1\leqslant j\leqslant k$) that
is not in the set $\\{1,...,k\\}$. Sub-rule (iii) is a simple result of the
general permutation rule.

Using (iii), it is easy to show that $\det M=D_{s}(M)=0$, thus
$\text{Rank}(M)\leqslant s-1$. Moreover, since each leading principal minor
$D_{k}(M)\geqslant 0$, $M_{\alpha\beta}$ is positive semidefinite.

Next, we prove $\text{Rank}(M)=s-1$ by contradiction—if
$\text{Rank}(M)\leqslant s-2$, then any minor of order $s-1$ should vanish.
However,

$\displaystyle
D_{s-1}(M)=\sideset{}{{}^{\prime}}{\sum}_{\alpha_{1},...,\alpha_{s-1}}^{s}C_{1\alpha_{1}}C_{2\alpha_{2}}...C_{s-1,\alpha_{s-1}}>0$
(90)

since it contains at least one term $C_{1s}C_{2s}...C_{s-1,s}$, a
contradiction (QED).

We want to emphasize that Prop. G.1 together with Eq. 86 is a statement about
UPC. It can be equivalently stated as “UPC leads to a positive semidefinite
$M$ and $\text{Rank}(M)=s-1$”. We cannot prove the most general case, but will
prove two special cases below, which still contain enough information for
understanding this general property.

###### Proposition G.2

Under the UPC, if an isolated band $m$ contains $s$ orbitals, then the Hessian
matrix $M$ is positive semidefinite and $\text{Rank}(M)=s-1$.

The condition above can be stated as “for each $1\leqslant\alpha\leqslant s$
there exists a neighborhood of some k-point in Brillouin zone such that
$|u_{m\textbf{k}\alpha}|\neq 0$”. Based on Prop. G.1, we only need to show
$C_{\alpha\beta}>0$. From Eq. IV we find for a single isolated band $m$,

$\displaystyle C_{\alpha\beta}=\sum_{\textbf{k}}\frac{\tanh(\beta
E_{m\textbf{k}}/2)}{2E_{m\textbf{k}}}|u_{m\textbf{k}\alpha}|^{2}|u_{m\textbf{k}\beta}|^{2},\,\,\,\alpha\neq\beta,$
(91)

which is positive. Therefore, the proposition is proved.

By induction, we know that whenever an orbital is removed from the band, one
positive eigenvalue of $M$ will approach 0, lowering the rank by 1. This can
be easily seen from that $|u_{m\textbf{k}r}|\rightarrow 0$ if the $r^{th}$
orbital is removed. Then $M$ will have an additional kernel eigenvector
$\textbf{e}_{r}=(...,0,1,0,...)^{T}$ (with its $r^{th}$ component nonzero).

###### Proposition G.3

Under the UPC, for a composite of $n$ bands formed from $s$ orbitals, if (i)
the composite is incomplete ($n<s$), (ii) the $n$ bands are energetically
close to each other compared to the interaction scale $\Delta$, and (iii) the
composite contains all the $s$ orbitals, then the Hessian matrix $M$ is
positive semidefinite and $\text{Rank}(M)=s-1$.

This proposition is based on three conditions, but the first two can be
relaxed. Similarly, the third condition can be stated as “for each
$1\leqslant\alpha\leqslant s$ there exists some band $l\in\mathcal{C}$ and a
neighborhood of some k-point in Brillouin zone such that
$|u_{l\textbf{k}\alpha}|\neq 0$”. The second condition states that the $n$
bands form an ICFB. With this condition one can omit the band index in
$E_{l\textbf{k}}$, $E_{l\textbf{k}}\rightarrow E_{\textbf{k}}$. From Eq. IV,
$C_{\alpha\beta}$ of the composite bands can be expressed as

$\displaystyle C_{\alpha\beta}$
$\displaystyle\approx\sum_{\textbf{k}}\frac{\tanh(\beta
E_{\textbf{k}}/2)}{2E_{\textbf{k}}}\sum_{l,l^{\prime}\in\mathcal{C}}u_{l^{\prime}\textbf{k}\alpha}^{*}u_{l\textbf{k}\alpha}u_{l\textbf{k}\beta}^{*}u_{l^{\prime}\textbf{k}\beta}$
$\displaystyle=\sum_{\textbf{k}}\frac{\tanh(\beta
E_{\textbf{k}}/2)}{2E_{\textbf{k}}}|b_{\textbf{k},\alpha\beta}|^{2},$ (92)

where we have defined

$\displaystyle
b_{\textbf{k},\alpha\beta}=\sum_{l\in\mathcal{C}}u_{l\textbf{k}\alpha}^{*}u_{l\textbf{k}\beta}.$
(93)

Eq. G is the leading term of $C_{\alpha\beta}$ if it is expanded in powers of
$\xi_{l\textbf{k}}-\xi_{l^{\prime}\textbf{k}}$. Since it is required that the
composite is incomplete ($n<s$), we have $|b_{\textbf{k},\alpha\beta}|$
strictly greater than 0. Then $C_{\alpha\beta}>0$, and the proposition is
proved. If instead $n=s$, then
$b_{\textbf{k},\alpha\beta}=\delta_{\alpha\beta}$, making the leading term of
$C_{\alpha\beta}$ vanish.

For cases other than those stated in Prop. G.3, one should be able to
numerically check that $C_{\alpha\beta}>0$ always holds under UPC. Analogous
to the single isolated band case, whenever an $r^{th}$ orbital is removed from
the composite, $\textbf{e}_{r}$ becomes an additional kernel eigenvector and
$\text{Rank}(M)$ is lowered by 1.

##  Appendix H Equivalence of the Two Definitions of Minimal Quantum Metric

In this appendix, we show how the two definitions of MQM given in Sec. IV are
equivalent to each other based on the rank information of matrix
$\partial^{2}\Omega/\partial\Delta^{I}_{\alpha}\partial\Delta^{I}_{\beta}$.

Consider the geometric transformation that each orbital’s position transforms
as $\textbf{x}_{\alpha}=\textbf{x}^{0}_{\alpha}+\delta\textbf{x}_{\alpha}$.
During this transformation, hopping integrals $t^{\sigma}_{ij,\alpha\beta}$
are fixed, while the Bloch components $u_{l\textbf{k}\alpha}$ undergo the
change

$\displaystyle
u_{l\textbf{k}\alpha}=e^{-i\textbf{k}\cdot\delta\textbf{x}_{\alpha}}u^{0}_{l\textbf{k}\alpha}.$
(94)

The intra- and inter-band quantum metric transform accordingly,
$g^{l}_{\mu\nu}=g^{l,0}_{\mu\nu}+\delta g^{l}_{\mu\nu}$,
$g^{ll^{\prime}}_{\mu\nu}=g^{ll^{\prime},0}_{\mu\nu}+\delta
g^{ll^{\prime}}_{\mu\nu}$, with

$\displaystyle\delta g^{l}_{\mu\nu}(\textbf{k})=\langle\delta x_{\mu}\delta
x_{\nu}\rangle^{0}_{l\textbf{k}}-\langle\delta
x_{\mu}\rangle^{0}_{l\textbf{k}}\langle\delta
x_{\nu}\rangle^{0}_{l\textbf{k}}+\big{[}\big{(}\frac{i}{2}\sum_{\alpha}\delta
x_{\alpha,\mu}u_{l\textbf{k}\alpha}^{0*}\partial_{\nu}u_{l\textbf{k}\alpha}^{0}+c.c.-i\langle\delta
x_{\mu}\rangle^{0}_{l\textbf{k}}\langle
u^{0}_{l\textbf{k}}|\partial_{\nu}u^{0}_{l\textbf{k}}\rangle\big{)}+\mu\leftrightarrow\nu\big{]},$
$\displaystyle\delta
g^{ll^{\prime}}_{\mu\nu}(\textbf{k})=\big{[}-\frac{1}{2}\langle\delta
x_{\mu}\rangle_{ll^{\prime}\textbf{k}}^{0}\langle\delta
x_{\nu}\rangle_{l^{\prime}l\textbf{k}}^{0}-\frac{i}{2}\big{(}\langle\delta
x_{\mu}\rangle_{ll^{\prime}\textbf{k}}^{0}\langle
u_{l^{\prime}\textbf{k}}^{0}|\partial_{\nu}u_{l\textbf{k}}^{0}\rangle+\langle\delta
x_{\nu}\rangle_{l^{\prime}l\textbf{k}}^{0}\langle
u_{l\textbf{k}}^{0}|\partial_{\mu}u_{l^{\prime}\textbf{k}}^{0}\rangle\big{)}\big{]}+\mu\leftrightarrow\nu.$
(95)

Here $\langle\rangle^{0}_{l\textbf{k}}$ represents the average over different
orbitals in the initial state $u^{0}_{l\textbf{k}}$, while
$\langle\rangle^{0}_{ll^{\prime}\textbf{k}}$ represents the overlap between
different bands, i.e.,

$\displaystyle\langle\delta
x_{\mu}\rangle^{0}_{l\textbf{k}}=\sum_{\alpha=1}^{s}\delta
x_{\alpha,\mu}|u^{0}_{l\textbf{k}\alpha}|^{2},$ $\displaystyle\langle\delta
x_{\mu}\delta x_{\nu}\rangle^{0}_{l\textbf{k}}=\sum_{\alpha=1}^{s}\delta
x_{\alpha,\mu}\delta x_{\alpha,\nu}|u^{0}_{l\textbf{k}\alpha}|^{2},$
$\displaystyle\langle\delta
x_{\mu}\rangle_{ll^{\prime}\textbf{k}}^{0}=\sum_{\alpha=1}^{s}\delta
x_{\mu,\alpha}u_{l\textbf{k}\alpha}^{0*}u_{l^{\prime}\textbf{k}\alpha}^{0}.$
(96)

With these established, let’s vary the functional $I$ in Eq. 20 around some
fixed positions $\\{\textbf{x}_{\alpha}^{0}\\}$. One can readily show the
following identity,

$\displaystyle\frac{\partial I}{\partial
x_{\alpha,\mu}}\bigg{|}_{\\{\textbf{x}_{\alpha}^{0}\\}}=\bigg{(}\frac{2i}{\Delta}\bigg{)}V^{0}_{\alpha,\mu},$
(97)

where $V_{\alpha,\mu}^{0}$ is Eq. IV evaluated at
$\\{\textbf{x}_{\alpha}^{0}\\}$. Based on Definition 2, if
$\\{\textbf{x}_{\alpha}^{0}\\}$ extremize $I$, then
$V^{0}_{\alpha,\mu}=M_{\alpha\beta}\partial_{\mu}\Delta_{\beta}^{0}=0$,
implying $D_{s}^{(2)}=0$ in Eq. 6. Therefore we have proved that Definition 2
leads to Definition 1.

To show the other way requires the rank of the Hessian matrix
$\partial^{2}\Omega/\partial_{\mu}\Delta_{\alpha}\partial_{\nu}\Delta_{\beta}$
to be $s-1$. The uniform pairing case has been proved in App. G. Assuming
$\text{Rank}(M)=s-1$, then the vanishing of $D_{s}^{(2)}$ in Eq. 6 has two
solutions: either (i) $\partial_{\mu}\Delta_{\alpha}=0$ or (ii)
$\partial_{\mu}\Delta_{1}=\partial_{\mu}\Delta_{2}=...=\partial_{\mu}\Delta_{s}=i\lambda_{\mu}$,
$\mu=x,y$. These two solutions are identical up to a $U(1)$ gauge
transformation [31, 1]. To see this, we write the second solution as

$\displaystyle\Delta_{\textbf{q},\alpha}=\Delta[1+i\bm{\lambda}\cdot\textbf{q}+S_{\alpha}(q^{2})],$
(98)

where $\bm{\lambda}=(\lambda_{x},\lambda_{y})$ is a real vector and
$S_{\alpha}(q^{2})$ is a complex function of $q^{2}$ order. A new order
parameter that has an overall phase difference can be defined,

$\displaystyle\widetilde{\Delta}_{\textbf{q},\alpha}=e^{-i\bm{\lambda}\cdot\textbf{q}}\Delta_{\textbf{q},\alpha},$
(99)

which satisfies $\partial_{\mu}\widetilde{\Delta}_{\alpha}=0$ and recovers the
first solution. Both solutions above lead to the vanishing of extremization
Eq. 97, since $V_{\alpha,\mu}=M_{\alpha\beta}\partial_{\mu}\Delta_{\alpha}$.
We have shown that Definition 1 implies Definition 2. Therefore, the two
definitions are equivalent.

With the knowledge that $\text{Rank}(M)=s-1$, one can further show the minimal
positions exist and are unique, up to fixing the position of one orbital [1].

Starting from some general positions $\\{\textbf{x}_{\alpha}^{0}\\}$, let’s
consider the transformation
$\textbf{x}_{\alpha}=\textbf{x}^{0}_{\alpha}+\delta\textbf{x}_{\alpha}$ again,
under which
$\partial_{\mu}\Delta_{\alpha}=\partial_{\mu}\Delta_{\alpha}^{0}+\delta\partial_{\mu}\Delta_{\alpha}$
and $V_{\alpha,\mu}=V^{0}_{\alpha,\mu}+\delta V_{\alpha,\mu}$. Using Eq. IV
and 94, the change of $V_{\alpha,\mu}$ is found to be

$\displaystyle\delta V_{\alpha,\mu}=-2i\Delta M_{\alpha\beta}\delta
x_{\beta,\mu}.$ (100)

By the invertibility of $M_{\alpha\beta}$ in the $(s-1)$-dim subspace, we
obtain

$\displaystyle\delta\partial_{\mu}\Delta_{\alpha}=-2i\Delta\delta
x_{\alpha,\mu}.$ (101)

Imposing
$\partial_{\mu}\Delta_{\alpha}^{0}+\delta\partial_{\mu}\Delta_{\alpha}=0$,
thus the minimal positions are located at

$\displaystyle\delta
x_{\alpha,\mu}=-\frac{i}{2\Delta}\partial_{\mu}\Delta_{\alpha}^{0},$ (102)

from the initial positions $\\{\textbf{x}_{\alpha}^{0}\\}$.

As discussed at the end of Sec. IV, when some orbitals are removed from the
composite bands, their positions become irrelevant to the geometric functional
$I$. Thus, the MQM orbital positions are hyperlines or hyperplanes in the
geometric space. We use two-orbital models as an example to illustrate this in
App. I.

##  Appendix I Two-orbital Models

Inserting Eq. H into Eq. II, and inserting Eq. 101 into Eq. V, we find both
$\delta D_{s}^{(1)}$ and $\delta D_{s}^{(2)}$ have a quadratic form of the
orbital coordinate change $\delta x_{\alpha,\mu}$, therefore take a parabolic
shape in the $3(s-1)$-dim geometric space. The eigenvalues of Hessian matrix
$M_{\alpha\beta}$ show the steepness of the parabola, while the eigenvectors
tell the directions where the MQM functional undergoes the steepest change.
Whenever an orbital is removed from the composite, one eigenvalue of
$M_{\alpha\beta}$ reaches zero, and its position becomes irrelevant to the
functional $I$. In this appendix, we illustrate this using two-orbital models.

For two-orbital models, one can define the composite to contain either one
(valence or conduction) band or both bands, depending on the Fermi energy and
interaction scale. Since the position of one orbital can permanently be fixed,
we write $I$ as a function $I(\textbf{r})$, where
$\textbf{r}=(x,y)\equiv\textbf{x}_{2}-\textbf{x}_{1}$ is orbital 2’s position
relative to orbital 1. To calculate $I(\textbf{r})$, one can use the method of
the Bloch function, i.e., using Eq. 94 to get H. Here, instead, we use the
Bloch Hamiltonian method.

We choose the lattice geometry where the two orbitals sit on the same site as
the reference geometry, $\textbf{r}=0$. This defines the Bloch Hamiltonian of
a two-orbital model

$\displaystyle
h(\textbf{k})=d_{0}(\textbf{k})\sigma_{0}+\sum_{j=x,y,z}d_{j}(\textbf{k})\sigma_{j}.$
(103)

(we omitted the spin index, so $h(\textbf{k})$ is for spin-$\uparrow$ only and
the spin-$\downarrow$ sector is its time-reversal counterpart) where
$\sigma_{j}$ are Pauli matrices in the orbital space. Then, the quantum metric
is

$\displaystyle g_{\mu\nu}\equiv
g^{v}_{\mu\nu}=g^{c}_{\mu\nu}=-g^{vc}_{\mu\nu}=\frac{1}{4}\partial_{\mu}\hat{\textbf{d}}\cdot\partial_{\nu}\hat{\textbf{d}},$
(104)

where $\textbf{d}=(d_{x},d_{y},d_{z})$ and
$\hat{\textbf{d}}=\textbf{d}/|\textbf{d}|$ is a unit vector.

Under the geometric transformation, $h(\textbf{k})$ is transformed by the
unitary matrix
$\mathcal{U}(\textbf{k})=\text{diag}\\{1,e^{i\textbf{k}\cdot\textbf{r}}\\}$:

$\displaystyle\widetilde{h}(\textbf{k})=\mathcal{U}(\textbf{k})^{\dagger}h(\textbf{k})\mathcal{U}(\textbf{k}),$
(105)

therefore vector d is rotated in the way

$\displaystyle\widetilde{d}_{x}=d_{x}\cos(\textbf{k}\cdot\textbf{r})+d_{y}\sin(\textbf{k}\cdot\textbf{r}),$
$\displaystyle\widetilde{d}_{y}=-d_{x}\sin(\textbf{k}\cdot\textbf{r})+d_{y}\cos(\textbf{k}\cdot\textbf{r}),$
$\displaystyle\widetilde{d}_{z}=d_{z},$ (106)

which gives the quantum metric change

$\displaystyle\widetilde{g}_{\mu\nu}=$ $\displaystyle
g_{\mu\nu}+\frac{1}{4}[(-\hat{d}_{x}\partial_{\mu}\hat{d}_{y}+\hat{d}_{y}\partial_{\mu}\hat{d}_{x})r_{\nu}$
$\displaystyle+(-\hat{d}_{x}\partial_{\nu}\hat{d}_{y}+\hat{d}_{y}\partial_{\nu}\hat{d}_{x})r_{\mu}+(\hat{d}_{x}^{2}+\hat{d}_{y}^{2})r_{\mu}r_{\nu}]$
(107)

where $\mu$ or $r_{\mu}=x$ or $y$. Then $I(\textbf{r})$ of Eq. 20 is

$\displaystyle I(\textbf{r})=I_{0}+I_{1x}x+I_{1y}y+I_{2}(x^{2}+y^{2}),$ (108)

with

$\displaystyle I_{0}=\sum_{\textbf{k}}f(\textbf{k})\text{Tr}g(\textbf{k}),$
$\displaystyle
I_{1\mu}=\frac{1}{2}\sum_{\textbf{k}}f(\textbf{k})(-\hat{d}_{x}\partial_{\mu}\hat{d}_{y}+\hat{d}_{y}\partial_{\mu}\hat{d}_{x}),$
$\displaystyle
I_{2}=\frac{1}{4}\sum_{\textbf{k}}f(\textbf{k})(\hat{d}_{x}^{2}+\hat{d}_{y}^{2}).$
(109)

Here, the form of $f(\textbf{k})$ depends on whether the composite contains
only one band or both bands. For the former,
$f(\textbf{k})=f_{v(c)}(\textbf{k})=\tanh(\beta
E_{v(c)\textbf{k}}/2)/E_{v(c)\textbf{k}}$; for the latter,

$\displaystyle f(\textbf{k})=$
$\displaystyle\frac{1}{E_{v\textbf{k}}}\tanh\frac{\beta
E_{v\textbf{k}}}{2}+\frac{1}{E_{c\textbf{k}}}\tanh\frac{\beta
E_{c\textbf{k}}}{2}$ $\displaystyle-\tanh\frac{\beta
E_{v\textbf{k}}}{2}\bigg{[}\frac{2p_{vc}^{(+)}(\textbf{k})}{E_{v\textbf{k}}+E_{c\textbf{k}}}+\frac{2p_{vc}^{(-)}(\textbf{k})}{E_{v\textbf{k}}-E_{c\textbf{k}}}\bigg{]}$
$\displaystyle-\tanh\frac{\beta
E_{c\textbf{k}}}{2}\bigg{[}\frac{2p_{vc}^{(+)}(\textbf{k})}{E_{v\textbf{k}}+E_{c\textbf{k}}}+\frac{2p_{vc}^{(-)}(\textbf{k})}{E_{c\textbf{k}}-E_{v\textbf{k}}}\bigg{]}.$
(110)

Notice that $I_{0},I_{1\mu},I_{2}$ are all geometry-independent quantities.
Completing the squares in Eq. 108 one get

$\displaystyle
I(\textbf{r})=I_{0}-\frac{I_{1x}^{2}+I_{1y}^{2}}{4I_{2}}+I_{2}\bigg{[}\bigg{(}x+\frac{I_{1x}}{2I_{2}}\bigg{)}^{2}+\bigg{(}y+\frac{I_{1y}}{2I_{2}}\bigg{)}^{2}\bigg{]},$
(111)

indicating that the extremal position is located at
$\textbf{r}_{min}=-(I_{1x},I_{1y})/(2I_{2})$. The sign of the parabola is
determined by $\text{Sgn}(I_{2})=\text{Sgn}(f(\textbf{k}))$. It turns out that
Eq. I is the same as Eq. 78, which we have proved to be positive. Therefore,
for two-orbital models, no matter whether the composite contains one band or
both the valence and conduction band, $I_{2}>0$, and the function is minimized
at the point $\textbf{r}_{min}$. This is another example showing the Hessian
matrix $M$ for a composite of two bands of a two-orbital model is positive
semidefinite and $\text{Rank}(M)=1$, with the first two conditions stated in
Prop. G.3 relaxed.

Next, we consider two prototype models, the Bravais-lattice BHZ model [77, 78]
and the non-Bravais-lattice Haldane model [79] to show how the parabola
evolves under topological phase transitions.

The BHZ model [77, 78] is defined on a square lattice. Choosing the geometric
gauge where both orbitals sit on the same site, the Hamiltonian has

$\displaystyle\textbf{d}(\textbf{k})=(t\sin k_{x},t\sin k_{y},m+t\cos
k_{x}+t\cos k_{y}).$ (112)

Inserting this into Eq. I, one gets $\textbf{r}_{min}=(0,0)$, regardless of
the composite, temperature and value of $m/t$. This is enforced by the $C_{4}$
symmetry of the model hopping graph.

Similarly, the Bloch Hamiltonian of Haldane model has

$\displaystyle
d_{0}(\textbf{k})=2t_{2}\cos\phi\sum_{i=1}^{3}\cos(\textbf{k}\cdot\textbf{a}_{i}),$
$\displaystyle
d_{x}(\textbf{k})=t_{1}[1+\cos(\textbf{k}\cdot\textbf{a}_{2})+\cos(\textbf{k}\cdot\textbf{a}_{3})],$
$\displaystyle
d_{y}(\textbf{k})=t_{1}[\sin(\textbf{k}\cdot\textbf{a}_{2})-\sin(\textbf{k}\cdot\textbf{a}_{3})],$
$\displaystyle
d_{z}(\textbf{k})=m-2t_{2}\sin\phi\sum_{i=1}^{3}\sin(\textbf{k}\cdot\textbf{a}_{i}).$
(113)

Inserting this into Eq. I we find $\textbf{r}_{min}=(0,\sqrt{3})$, regardless
of the composite, temperature, $t_{1},t_{2}$ and $m$ values. This is enforced
by the $C_{6}$ symmetry of the hopping graph.

As $m/t\rightarrow\infty$, both the BHZ and Haldane model transit to the
topologically trivial phase, in which $\hat{d}_{x}^{2}+\hat{d}_{y}^{2}\propto
t^{2}/m^{2}\rightarrow 0$ (during the transition we assume $t,t_{1},t_{2}$ are
fixed but $m$ is varied). Suppose the composite contains the isolated valence
or conduction band only ($\Delta\sim t\ll m$). In that case, we find
$f(\textbf{k})\sim 1/\Delta$ which has no dependence on $m$, therefore
$I_{2}\propto f(\textbf{k})(\hat{d}_{x}^{2}+\hat{d}_{y}^{2})\rightarrow 0$,
implying that the parabola in Eq. 111 becomes flat. This means the functional
$I(\textbf{r})$ becomes a constant in the geometric space as one orbital is
removed from the band.

Conversely, if the composite contains both the valence and conduction band
($\Delta\sim m\gg t$), it always contains two orbitals, even in the
topologically trivial phase. In the atomic limit, from Eq. I, one can show

$\displaystyle f(\textbf{k})\sim\frac{m^{2}}{\Delta^{3}}.$ (114)

Therefore, $I_{2}\sim t^{2}/\Delta^{3}$ remains finite, and the rank of the
Hessian matrix does not decrease. This means $I(\textbf{r})$ of Eq. 111
remains parabolic in the atomic limit.

##  Appendix J Hessian Matrix of Nonuniform Pairing

This appendix studies the properties of the matrix of second derivative of
grand potential with order parameters beyond the uniform pairing. In
particular, we show that this matrix may not be positive semidefinite.

We restrict the interaction to be onsite intra-orbital density-density
interaction. We also assume a time-reversal invariant channel $\hat{\Delta}$
in orbital space, i.e., under TRS, the matrix $\hat{\Delta}$ is mapped to
itself up to a $U(1)$ phase. Then for our interaction type, the pairing matrix
at $\textbf{q}=0$ takes the diagonal form
$\hat{\Delta}=\text{diag}\\{\Delta_{1},\Delta_{2},...,\Delta_{s}\\}$, whose
entries are taken to be real but may be non-uniform.

For simplicity, we consider the case of a single isolated band at $T=0$, with
the gap equation given by Eq. IV. The $\textbf{q}=0$ gap equation reads

$\displaystyle\Delta_{\alpha}=$
$\displaystyle\frac{U_{\alpha}}{N}\sum_{\textbf{k}}|u_{m\textbf{k}\alpha}|^{2}\frac{\Delta_{m,\textbf{k}}(0)}{2E_{m\textbf{k}}},\,\,\,\,\,\,1\leqslant\alpha\leqslant
s.$ (115)

Here $\Delta_{m,\textbf{k}}(0)=\langle
u_{m\textbf{k}}|\hat{\Delta}|u_{m\textbf{k}}\rangle$ is the band-projected gap
at $\textbf{q}=0$, and
$E_{m\textbf{k}}=\sqrt{\xi_{m\textbf{k}}^{2}+\Delta_{m,\textbf{k}}(0)^{2}}$ is
the quasiparticle energy.

Following the discussions at the end of App. C and computing the derivative of
Eq. IV, we find

$\displaystyle
M_{\alpha\beta}=\frac{\delta_{\alpha\beta}}{\Delta_{\alpha}}\sum_{\textbf{k}}\frac{\Delta_{m,\textbf{k}}(0)}{2E_{m\textbf{k}}}|u_{m\textbf{k}\alpha}|^{2}-\sum_{\textbf{k}}\frac{1}{2E_{m\textbf{k}}}|u_{m\textbf{k}\alpha}|^{2}|u_{m\textbf{k}\beta}|^{2},$
(116)

and

$\displaystyle V_{\alpha,\mu}=\sum_{\textbf{k}}\frac{1}{2E_{m\textbf{k}}}[$
$\displaystyle(u_{m\textbf{k}\alpha}^{*}\partial_{\mu}u_{m\textbf{k}\alpha}-\partial_{\mu}u_{m\textbf{k}\alpha}^{*}u_{m\textbf{k}\alpha})\Delta_{m,\textbf{k}}(0)$
$\displaystyle-2|u_{m\textbf{k}\alpha}|^{2}\text{Im}\langle
u_{m\textbf{k}}|\hat{\Delta}|\partial_{\mu}u_{m\textbf{k}}\rangle].$ (117)

Now $U(1)$ symmetry imposes
$\sum_{\beta=1}^{s}M_{\alpha\beta}\Delta_{\beta}=0$ and
$\sum_{\beta=1}^{s}V_{\alpha,\mu}\Delta_{\alpha}=0$, so the kernel vector is
instead $\textbf{v}_{0}=(\Delta_{1},\Delta_{2},...,\Delta_{s})^{T}$.

For the most general case, $\text{Rank}(M)=s-1$. To study the
semidefiniteness, we take two-orbital models as an example. Note that the off-
diagonal elements of $M_{\alpha\beta}$ are all negative, so we write

$\displaystyle M_{\alpha\beta}=\begin{pmatrix}a&-b\\\ -b&c\end{pmatrix},$
(118)

with $b>0$. Kernel eigenvector $\textbf{v}_{0}$ leads to
$a=(\Delta_{2}/\Delta_{1})b$ and $c=(\Delta_{1}/\Delta_{2})b$. Since $a,c$ are
the diagonal entries of $M$, we conclude that for the two-orbital case, the
Hessian matrix is positive semidefinite if and only if
$\text{Sgn}(\Delta_{1}\Delta_{2})=1$. This result largely but not completely
coincides with the pair density-wave transition discovered in Ref. [58] when
one Hubbard interaction turns repulsive, $U_{1}U_{2}<0$. However, as a
reminder, there it was also found that a solution channel with
$\text{Sgn}(\Delta_{1}\Delta_{2})=-1$ exists when both interactions are
attractive, $U_{1},U_{2}>0$.

For models with more than two orbitals, similar but more complicated
conditions will be required for a nonpositive definite Hessian matrix.

## References

  * Huhtinen  _et al._ [2022] K.-E. Huhtinen, J. Herzog-Arbeitman, A. Chew, B. A. Bernevig, and P. Törmä, Phys. Rev. B 106, 014518 (2022).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018a).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018b).
  * Chen  _et al._ [2019] G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, J. Jung, Z. Shi, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Nature 572, 215 (2019).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature 574, 653 (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Burg  _et al._ [2019] G. W. Burg, J. Zhu, T. Taniguchi, K. Watanabe, A. H. MacDonald, and E. Tutuc, Phys. Rev. Lett. 123, 197702 (2019).
  * Shen  _et al._ [2020] C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, T. Taniguchi, R. Yang, Z. Y. Meng, D. Shi, O. V. Yazyev, and G. Zhang, Nature Physics 16, 520 (2020).
  * Cao  _et al._ [2020] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 583, 215 (2020).
  * Liu  _et al._ [2020] X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. Haei Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Nature 583, 221 (2020).
  * Polshyn  _et al._ [2020] H. Polshyn, J. Zhu, M. A. Kumar, Y. Zhang, F. Yang, C. L. Tschirhart, M. Serlin, K. Watanabe, T. Taniguchi, A. H. MacDonald, and A. F. Young, Nature 588, 66 (2020).
  * Chen  _et al._ [2021] S. Chen, M. He, Y.-H. Zhang, V. Hsieh, Z. Fei, K. Watanabe, T. Taniguchi, D. H. Cobden, X. Xu, C. R. Dean, and M. Yankowitz, Nature Physics 17, 374 (2021).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021).
  * Cao  _et al._ [2021] Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 595, 526 (2021).
  * Uchoa and Barlas [2013] B. Uchoa and Y. Barlas, Phys. Rev. Lett. 111, 046604 (2013).
  * Peotta and Törmä [2015] S. Peotta and P. Törmä, Nature Communications 6, 8944 (2015).
  * Julku  _et al._ [2016] A. Julku, S. Peotta, T. I. Vanhala, D.-H. Kim, and P. Törmä, Phys. Rev. Lett. 117, 045303 (2016).
  * Tovmasyan  _et al._ [2016] M. Tovmasyan, S. Peotta, P. Törmä, and S. D. Huber, Phys. Rev. B 94, 245149 (2016).
  * Liang  _et al._ [2017] L. Liang, T. I. Vanhala, S. Peotta, T. Siro, A. Harju, and P. Törmä, Phys. Rev. B 95, 024515 (2017).
  * Törmä  _et al._ [2018] P. Törmä, L. Liang, and S. Peotta, Phys. Rev. B 98, 220511 (2018).
  * Iskin [2018] M. Iskin, Phys. Rev. A 97, 063625 (2018).
  * Hu  _et al._ [2019] X. Hu, T. Hyart, D. I. Pikulin, and E. Rossi, Physical Review Letters 123, 237002 (2019).
  * Xie  _et al._ [2020] F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Phys. Rev. Lett. 124, 167002 (2020).
  * Wang  _et al._ [2020] Z. Wang, G. Chaudhary, Q. Chen, and K. Levin, Phys. Rev. B 102, 184504 (2020).
  * Hofmann  _et al._ [2020] J. S. Hofmann, E. Berg, and D. Chowdhury, Phys. Rev. B 102, 201112 (2020).
  * Julku  _et al._ [2020] A. Julku, T. J. Peltonen, L. Liang, T. T. Heikkilä, and P. Törmä, Phys. Rev. B 101, 060505 (2020).
  * Julku  _et al._ [2021] A. Julku, G. M. Bruun, and P. Törmä, Phys. Rev. Lett. 127, 170404 (2021).
  * Wang  _et al._ [2021] J. Wang, J. Cano, A. J. Millis, Z. Liu, and B. Yang, Phys. Rev. Lett. 127, 246403 (2021).
  * Törmä  _et al._ [2022] P. Törmä, S. Peotta, and B. A. Bernevig, Nature Reviews Physics 4, 528 (2022).
  * Chan  _et al._ [2022] S. M. Chan, B. Grémaud, and G. G. Batrouni, Phys. Rev. B 105, 024502 (2022).
  * Herzog-Arbeitman  _et al._ [2022a] J. Herzog-Arbeitman, A. Chew, K.-E. Huhtinen, P. Törmä, and B. A. Bernevig, Many-body superconductivity in topological flat bands (2022a), arXiv:2209.00007 [cond-mat.str-el] .
  * Peltonen  _et al._ [2018] T. J. Peltonen, R. Ojajärvi, and T. T. Heikkilä, Phys. Rev. B 98, 220504 (2018).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8, 031089 (2018).
  * Xu and Balents [2018] C. Xu and L. Balents, Phys. Rev. Lett. 121, 087001 (2018).
  * Dodaro  _et al._ [2018] J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q. Sun, and C. Wang, Phys. Rev. B 98, 075154 (2018).
  * Xu  _et al._ [2018] X. Y. Xu, K. T. Law, and P. A. Lee, Phys. Rev. B 98, 121406 (2018).
  * Liu  _et al._ [2018] C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Phys. Rev. Lett. 121, 217001 (2018).
  * Isobe  _et al._ [2018] H. Isobe, N. F. Q. Yuan, and L. Fu, Phys. Rev. X 8, 041041 (2018).
  * González and Stauber [2019] J. González and T. Stauber, Phys. Rev. Lett. 122, 026801 (2019).
  * Cea  _et al._ [2019] T. Cea, N. R. Walet, and F. Guinea, Phys. Rev. B 100, 205113 (2019).
  * Chou  _et al._ [2021] Y.-Z. Chou, F. Wu, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 127, 217001 (2021).
  * Phong  _et al._ [2021] V. o. T. Phong, P. A. Pantaleón, T. Cea, and F. Guinea, Phys. Rev. B 104, L121116 (2021).
  * Wu  _et al._ [2018] F. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. 121, 257001 (2018).
  * Wu  _et al._ [2019] F. Wu, E. Hwang, and S. Das Sarma, Phys. Rev. B 99, 165112 (2019).
  * Lian  _et al._ [2019] B. Lian, Z. Wang, and B. A. Bernevig, Phys. Rev. Lett. 122, 257002 (2019).
  * Choi and Choi [2018] Y. W. Choi and H. J. Choi, Phys. Rev. B 98, 241412 (2018).
  * Ojajärvi  _et al._ [2018] R. Ojajärvi, T. Hyart, M. A. Silaev, and T. T. Heikkilä, Phys. Rev. B 98, 054515 (2018).
  * Cea and Guinea [2021] T. Cea and F. Guinea, Proceedings of the National Academy of Sciences 118, e2107874118 (2021).
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Science Advances 7, eabf5299 (2021).
  * Scalapino  _et al._ [1992] D. J. Scalapino, S. R. White, and S. C. Zhang, Phys. Rev. Lett. 68, 2830 (1992).
  * Scalapino  _et al._ [1993] D. J. Scalapino, S. R. White, and S. Zhang, Phys. Rev. B 47, 7995 (1993).
  * Schrieffer [1999] J. Schrieffer,  _Theory Of Superconductivity_ , Advanced Books Classics (Avalon Publishing, 1999).
  * Leggett [1998] A. J. Leggett, Journal of Statistical Physics 93, 927 (1998).
  * Tian  _et al._ [2023] H. Tian, X. Gao, Y. Zhang, S. Che, T. Xu, P. Cheung, K. Watanabe, T. Taniguchi, M. Randeria, F. Zhang, C. N. Lau, and M. W. Bockrath, Nature 614, 440 (2023).
  * Kitamura  _et al._ [2022] T. Kitamura, A. Daido, and Y. Yanase, Phys. Rev. B 106, 184507 (2022).
  * Chen and Huang [2023] W. Chen and W. Huang, Science China Physics, Mechanics & Astronomy 66, 287212 (2023).
  * Jiang and Barlas [2023] G. Jiang and Y. Barlas, Phys. Rev. Lett. 131, 016002 (2023).
  * Mera and Mitscherling [2022] B. Mera and J. Mitscherling, Phys. Rev. B 106, 165133 (2022).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
  * Zak [1980] J. Zak, Phys. Rev. Lett. 45, 1025 (1980).
  * Simon and Rudner [2020] S. H. Simon and M. S. Rudner, Phys. Rev. B 102, 165148 (2020).
  * [63] The concept of interaction window is well-defined self-consistently. Starting from Hubbard models, the window can be first defined as the range of $2|U|$ around the Fermi level, with $U$ the largest Hubbard parameter. By solving the self-consistency gap equation and restricting to a channel, the order parameter $\Delta$ is obtained, which is typically smaller than $|U|$, so the window can be refined to $2\Delta$. After this step, some bands falling between $\Delta$ and $|U|$ may be neglected. For pairing channels that are not uniform, matrix $\hat{\Delta}$ is not proportional to identity, then the window can be defined as $2\text{max}\\{|\hat{\Delta}_{\alpha\beta}|\\}$.
  * Tinkham [2004] M. Tinkham,  _Introduction to superconductivity_ (Courier Corporation, 2004).
  * Kopnin [2011] N. B. Kopnin, JETP letters 94, 81 (2011).
  * Herzog-Arbeitman  _et al._ [2022b] J. Herzog-Arbeitman, V. Peri, F. Schindler, S. D. Huber, and B. A. Bernevig, Phys. Rev. Lett. 128, 087002 (2022b).
  * Böhm  _et al._ [2003] A. Böhm, A. Mostafazadeh, H. Koizumi, Q. Niu, and J. Zwanziger,  _The Geometric phase in quantum systems: foundations, mathematical concepts, and applications in molecular and condensed matter physics_ (Springer, 2003).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Ozawa and Mera [2021] T. Ozawa and B. Mera, Phys. Rev. B 104, 045103 (2021).
  * Cano  _et al._ [2018] J. Cano, B. Bradlyn, Z. Wang, L. Elcoro, M. G. Vergniory, C. Felser, M. I. Aroyo, and B. A. Bernevig, Phys. Rev. B 97, 035139 (2018).
  * Zak [1982] J. Zak, Phys. Rev. B 26, 3010 (1982).
  * [72] Two EBRs may share some atomic orbital components, in which case, since the inter-band quantum metric is negative semidefinite, it will still reduce the superfluid weight. However, this reduction is not substantial, as it is unable to change the lower bound by an integer.
  * [73] This geometric transformation can be understood as physical, which is to shift the actual position of atomic orbitals but keep the hopping integrals fixed, or as artificial, which is to impose a componentwise phase for the Bloch function.
  * [74] This projection is not orthogonal as the resulting vector is not orthogonal to $\textbf{v}_{0}=(1,1,...,1)^{T}$. An orthogonal projection is not necessary since the purpose is only to get a simple expression of $D_{s}^{(2)}$.
  * Hofmann  _et al._ [2022] J. S. Hofmann, D. Chowdhury, S. A. Kivelson, and E. Berg, npj Quantum Materials 7, 83 (2022).
  * Hofmann  _et al._ [2023] J. S. Hofmann, E. Berg, and D. Chowdhury, Phys. Rev. Lett. 130, 226001 (2023).
  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
  * Qi  _et al._ [2006] X.-L. Qi, Y.-S. Wu, and S.-C. Zhang, Phys. Rev. B 74, 085308 (2006).
  * Haldane [1988] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).