Topology-Bounded Superfluid Weight In Twisted Bilayer Graphene

TL;DR

This paper demonstrates that the nontrivial topology of flat bands in twisted bilayer graphene enhances the superfluid weight, which is crucial for superconductivity, through a topological lower bound related to the Fubini-Study metric.

Contribution

It derives the superfluid weight for TBLG flat bands as an integral of the Fubini-Study metric and links it to topological invariants, revealing a topological lower bound.

Findings

Superfluid weight is expressed as an integral of the Fubini-Study metric.

Topological lower bound of superfluid weight is proportional to the pairing gap.

Trivial flat bands have zero superfluid weight.

Abstract

While regular flat bands are good for enhancing the density of states and hence the gap, they are detrimental to the superfluid weight. We show that the predicted nontrivial topology of the two lowest flat bands of twisted bilayer graphene plays an important role in the enhancement of the superfluid weight and hence of superconductivity. We derive the superfluid weight (phase stiffness) of the TBLG superconducting flat bands with a uniform pairing, and show that it can be expressed as an integral of the Fubini-Study metric of the flat bands. This mirrors results already obtained for nonzero Chern number bands even though the TBLG flat bands have zero Chern number. We further show the metric integral is lower bounded by the topological $C_{2z}T$ Wilson loop winding number of the TBLG flat bands, which renders the superfluid weight has a topological lower bound proportional to the pairing gap. In contrast, trivial flat bands have a zero superfluid weight. The superfluid weight is crucial in determining the BKT transition temperature of the superconductor. Based on the transition temperature measured in TBLG experiments, we estimate the topological contribution of the superfluid weight in TBLG.

Figures

Equations

$$j_i=-[D_s{]}_{ij}{A}_j,$$

$$H_{\mathrm{BdG}}=$$

$$+\frac{1}{2}\sum _{\mathbf{k}}\mathrm{Tr}(\mathcal{H}(\mathbf{k})-\mu ).$$

$$[D_s{]}_{ij}=\frac{8e^2\Delta }{{\hslash }^2}\sqrt{\nu (1-\nu )}\int \frac{d^{2}k}{(2\pi {)}^2}{g}_{ij}(\mathbf{k}),$$

$$g_{ij}(\mathbf{k})=$$

$$+u^{†}(\mathbf{k}){\partial }_{k_i}u(\mathbf{k})u^{†}(\mathbf{k}){\partial }_{k_j}u(\mathbf{k}),$$

$$\Delta (\mathbf{k})=\left[{\Delta }_1{\tilde{U}}_{\mathbf{k}}{\tilde{U}}_{\mathbf{k}}^{†}+{\Delta }_2\left(\text{mathds}1-{\tilde{U}}_{\mathbf{k}}{\tilde{U}}_{\mathbf{k}}^{†}\right)\right]U_T,$$

$$[D_s{]}_{ij}=\frac{2e^2|{\Delta }_1|}{{\hslash }^2}\left(1+\frac{{\Delta }_2}{{\Delta }_1}\right)\sqrt{\nu (1-\nu )}\int \frac{d^{2}k}{(2\pi {)}^2}{g}_{ij}(\mathbf{k})$$

$$g_{ij}(\mathbf{k})=$$

$$+\left({\tilde{U}}_{\mathbf{k}}^{†}{\partial }_{k_i}{\tilde{U}}_{\mathbf{k}}{\tilde{U}}_{\mathbf{k}}^{†}{\partial }_{k_j}{\tilde{U}}_{\mathbf{k}}\right)],$$

$$e_2=\frac{1}{2\pi }\int d^{2}k{f}_{xy}.$$

$$D_s\approx \frac{4e^2{\Delta }_1}{\pi {\hslash }^2}\sqrt{\nu (1-\nu )}|e_2|.$$

$$H_0=\sum _{\mathbf{k}}{\psi }_{\mathbf{k}\alpha }^{†}{\mathcal{H}}_{\alpha \beta }(\mathbf{k}){\psi }_{\mathbf{k}\beta },$$

$$H_0(\mathbf{q})=\sum _{\mathbf{k}}{\psi }_{\mathbf{k}\alpha }^{†}{\mathcal{H}}_{\alpha \beta }(\mathbf{k}-\mathbf{q}){\psi }_{\mathbf{k}\beta }.$$

$$H_{\mathrm{pair}}=\frac{1}{2}\sum _{\mathbf{k}}\left({\Delta }_{\alpha \beta }(\mathbf{k}){\psi }_{\mathbf{k}\alpha }^{†}{\psi }_{-\mathbf{k}\beta }^{†}+h.c.\right),$$

$$H_{\mathrm{MF}}(\mathbf{q})=$$

$$=$$

$$+\frac{1}{2}\sum _{\mathbf{k}\alpha \beta }{\delta }_{\alpha \beta }\left[{\mathcal{H}}_{\alpha \beta }(\mathbf{k}-\mathbf{q})-\mu {\delta }_{\alpha \beta }\right]+\frac{1}{2}\sum _{\mathbf{k}\alpha \beta }\left({\Delta }_{\alpha \beta }(\mathbf{k}){\psi }_{\mathbf{k}\alpha }^{†}{\psi }_{-\mathbf{k}\beta }^{†}+h.c.\right)$$

$$=$$

$$H_{\mathrm{MF}}(\mathbf{q})=$$

$$=\frac{1}{2}\sum _{\mathbf{k}}{d}_{\mathbf{k}}^{†}{\mathcal{H}}_{\mathbf{k}}(\mathbf{q})d_{\mathbf{k}}+\frac{1}{2}\sum _{\mathbf{k}}\mathrm{Tr}({\epsilon }_{\mathbf{k}-\mathbf{q}}-\mu )$$

$${\mathcal{H}}_{\mathbf{k}}(\mathbf{q})$$

$$d_{\mathbf{k}}$$

$$E_{\mathbf{k}}(\mathbf{q})=W_{\mathbf{k}}^{†}(\mathbf{q}){\mathcal{H}}_{\mathbf{k}}(\mathbf{q})W_{\mathbf{k}}(\mathbf{q}),$$

$$W_{\mathbf{k}}(0)=\left(\begin{array}[]{cc}\mathcal{U}_{\mathbf{k}}&\mathcal{V}_{\mathbf{k}}\\ \mathcal{V}^{*}_{-\mathbf{k}}&\mathcal{U}^{*}_{-\mathbf{k}}\end{array}\right)\,,$$

$$\Omega (\mathbf{q})$$

$$=\frac{1}{4}\sum _{\mathbf{k}}\left(\sum _{E_{\mathbf{k}n}(\mathbf{q})<0}{E}_{\mathbf{k}n}(\mathbf{q})+\sum _{E_{\mathbf{k}n}(\mathbf{q})>0}{E}_{\mathbf{k}n}(\mathbf{q})+\sum _{E_{\mathbf{k}n}(\mathbf{q})<0}{E}_{\mathbf{k}n}(\mathbf{q})-\sum _{E_{\mathbf{k}n}(\mathbf{q})>0}{E}_{\mathbf{k}n}(\mathbf{q})\right)+\frac{1}{2}\sum _{\mathbf{k}}\mathrm{Tr}({\epsilon }_{\mathbf{k}-\mathbf{q}}-\mu )$$

$$=\frac{1}{4}\sum _{\mathbf{k},n}{E}_{\mathbf{k}n}(\mathbf{q})-\frac{1}{4}\sum _{\mathbf{k},n}|E_{\mathbf{k}n}(\mathbf{q})|+\frac{1}{2}\sum _{\mathbf{k}}\mathrm{Tr}({\epsilon }_{\mathbf{k}-\mathbf{q}}-\mu )$$

$$=\frac{1}{4}\sum _{\mathbf{k}}\mathrm{Tr}{\mathcal{H}}_{\mathbf{k}}(\mathbf{q})-\frac{1}{4}\sum _{\mathbf{k},n}|E_{\mathbf{k}n}(\mathbf{q})|+\frac{1}{2}\sum _{\mathbf{k}}\mathrm{Tr}({\epsilon }_{\mathbf{k}-\mathbf{q}}-\mu )$$

$$=\frac{1}{4}\left(\sum _{\mathbf{k}}\mathrm{Tr}({\epsilon }_{\mathbf{k}-\mathbf{q}}-\mu )-\sum _{\mathbf{k}}\mathrm{Tr}({\epsilon }_{-\mathbf{k}-\mathbf{q}}-\mu )\right)-\frac{1}{4}\sum _{\mathbf{k},n}|E_{\mathbf{k}n}(\mathbf{q})|+\frac{1}{2}\sum _{\mathbf{k}}\mathrm{Tr}({\epsilon }_{\mathbf{k}-\mathbf{q}}-\mu )$$

$$=-\frac{1}{4}\sum _{\mathbf{k},n}|E_{\mathbf{k}n}(\mathbf{q})|+\frac{1}{2}\sum _{\mathbf{k}}\mathrm{Tr}({\epsilon }_{\mathbf{k}-\mathbf{q}}-\mu ),$$

$${\mathcal{D}}_{\mathbf{k}}(\mathbf{q})=U^{†}(\mathbf{k}-\mathbf{q})\Delta (\mathbf{k})U^{\ast }(-\mathbf{k}-\mathbf{q}),$$

$$[D_{s}]_{ij}=\frac{1}{V}\frac{\partial^{2}\Omega(\mathbf{q})}{\partial{q_{i}}\partial{q_{j}}}\Big{|}_{\mathbf{q}=0}=-\frac{1}{4V}\sum_{\mathbf{k}n}{\rm sgn}(E_{\mathbf{k}n}(\mathbf{q}))\partial_{q_{i}}\partial_{q_{j}}E_{\mathbf{k}n}(\mathbf{q})]\Big{|}_{\mathbf{q}=0}+\frac{1}{2V}\sum_{\mathbf{k}}{\rm Tr}\,\partial_{k_{i}}\partial_{k_{j}}\varepsilon_{\mathbf{k}}\,.$$

$${\partial }_{q_i}{E}_{\mathbf{k}n}(\mathbf{q})$$

$$=\sum _{\alpha ,\beta }{W}_{\mathbf{k}}^{\ast }(\mathbf{q}{)}_{\alpha ,n}{\partial }_{q_i}{\mathcal{H}}_{\mathbf{k}}(\mathbf{q}{)}_{\alpha ,\beta }{W}_{\mathbf{k}}(\mathbf{q}{)}_{\beta ,n}={\left[W_{\mathbf{k}}^{†}(\mathbf{q}){\partial }_{q_i}{\mathcal{H}}_{\mathbf{k}}(\mathbf{q})W_{\mathbf{k}}(\mathbf{q})\right]}_{n,n},$$

$${\partial }_{q_i}{\partial }_{q_j}{E}_{\mathbf{k}n}(\mathbf{q})=[W_{\mathbf{k}}^{†}(\mathbf{q}){\partial }_{q_i}{\partial }_{q_j}{\mathcal{H}}_{\mathbf{k}}(\mathbf{q})W_{\mathbf{k}}(\mathbf{q}){]}_{n,n}+[{\partial }_{q_i}{W}_{\mathbf{k}}^{†}(\mathbf{q}){\partial }_{q_j}{\mathcal{H}}_{\mathbf{k}}(\mathbf{q})W_{\mathbf{k}}(\mathbf{q}){]}_{n,n}+[W_{\mathbf{k}}^{†}(\mathbf{q}){\partial }_{q_j}{\mathcal{H}}_{\mathbf{k}}(\mathbf{q}){\partial }_{q_i}{W}_{\mathbf{k}}(\mathbf{q}){]}_{n,n}.$$

$$\sum _{\alpha \beta }{W}_{\mathbf{k}}^{\ast }(\mathbf{q}{)}_{\alpha ,n}{\mathcal{H}}_{\mathbf{k}}(\mathbf{q}{)}_{\alpha ,\beta }{W}_{\mathbf{k}}(\mathbf{q}{)}_{\beta ,m}=0.$$

Full text

The recently discovered superconducting phase in twisted bilayer graphene has received extensive attention Suárez Morell et al. (2010); Bistritzer and MacDonald (2011); Cao et al. (2018); Yankowitz et al. (2019); Huang et al. (2018); Yuan and Fu (2018); Po et al. (2018); Xu and Balents (2018); Roy and Juričić (2019); Dodaro et al. (2018); Wu et al. (2018); Isobe et al. (2018); Huang et al. (2018); You and Vishwanath (2018); Wu et al. (2018a); Kang and Vafek (2018); Kennes et al. (2018); Guinea and Walet (2018); Peltonen et al. (2018); Fidrysiak et al. (2018); Zou et al. (2018); González and Stauber (2019); Su and Lin (2018); Guo et al. (2018); Song et al. (2018); Po et al. (2018); Ahn et al. (2018); Hejazi et al. (2019); Laksono et al. (2018); Tarnopolsky et al. (2019); Venderbos and Fernandes (2018); Chen et al. (2018); Stauber et al. (2018); Choi and Choi (2018); Jian and Xu (2018); Kozii et al. (2018); Kang and Vafek (2018); Wu et al. (2018b, c); Liu et al. (2018); Lian et al. (2018a); Polshyn et al. (2019). The topology of the lowest two bands (per spin and valley) of twisted bilayer graphene (TBLG) is currently under debate Zou et al. (2018); Song et al. (2018); Po et al. (2018); Ahn et al. (2018); Po et al. (2018). Although theoretical models suggest a nontrivial topological number of these bands, the experimentally measurable effects through which one could prove or falsify the predicted nontrivial topology are scarce. Currently, one viable experimentally observable effect Lian et al. (2018b), predicts that the single-particle magnetic field spectrum of a topologically nontrivial band can cross the single-particle gap, in stark contrast to conventional knowledge and to the in-field spectrum of trivial bands. We here present another effect of a set of topologically nontrivial bands observable at zero field (of the kind present in TBLG) that appears when these bands become superconducting. We show that the superfluid weight in the superconducting state is the sum of two terms: a conventional term, which vanishes when the bands are perfectly flat, and a topological term, bounded from below by the Wilson loop winding number of the $C_{2z}T$ protected topology in TBLG.

This letter is organized as follows. First, we show that by assuming perfectly flat bands and $s$ wave pairing, the superfluid weight can be written as the integral of Fubini-Study metric over the Brillouin zone (BZ), and show that it is lower-bounded by the Wilson loop winding. Secondly, by applying this result to TBLG, we estimate the topological contribution of superfluid weight and explain the relatively high transition temperature.

The two characterizing features of superconductors are the zero DC resistance and Meissner effect. Both of these properties are captured by the celebrated London equation Tinkham (2004). It tells us that the electric current in a superconductor $\mathbf{j}$ is proportional to the gauge potential $\mathbf{A}$ under Coulomb gauge:

[TABLE]

in which the coefficient $[D_{s}]_{ij}$ is called the superfluid weight and it is a tensor in general. Some spacial symmetry, such as $C_{3z}$, requires it to be isotropic in 2D. In some works, such as Ref. Corson et al. (1999), it is called “phase stiffness” (describing energy susceptibility with respect to phase twists). The London equation has two kinds of consequences. One is the perfect diamagnetism feature of superconductor, and the other one is the frequency dependence of AC conductance. In 2D, the superfluid weight is also related with the transition temperature. The phase coherence will disappear at a temperature $T=T_{c}$ given by $\frac{\hbar^{2}D_{s}(T_{c})}{e^{2}k_{B}T_{c}}=\frac{8}{\pi}$ (known as Berezinskii-Kosterlitz-Thouless (BKT) transition Berezinsky (1971); Kosterlitz and Thouless (1973)), because of the creation of vortex-antivortex pairs. Usually $D_{s}(T)$ decreases with increasing temperature, so the transition temperature $T_{c}$ is always lower than $\frac{\pi\hbar^{2}D_{s}(0)}{8e^{2}k_{B}}$. Thus a small $D_{s}$ at zero temperature leads to a low transition temperature.

Two methods, using a static magnetic field or an alternating electric field, are used to measure the superfluid weight (phase stiffness) directly. In 3D, $D_{s}$ can be obtained by measuring the penetration depth $\lambda_{L}=\sqrt{1/\mu_{0}D_{s}}$ of a static magnetic field into the superconductor. This method, however, fails in 2D, as the penetration depth is no longer a simple function of $D_{s}$. Instead, one can measure the AC conductance under an alternating electric field of frequency $\omega$. By taking time derivative on both sides of Eq. (1), we obtain the frequency dependent conductance, and it is inversely proportional to the frequency as $\sigma(\omega)=\frac{iD_{s}}{\omega}$. Therefore the superfluid weight will be $D_{s}=-i\omega\sigma(\omega)$. In experiments, $\sigma(\omega)$ can be measured by the time-domain transmission spectroscopy without any contact with the sample Grüner and Dahl (1998). For example, at zero temperature, the stiffness temperature of $\rm Bi_{2}Sr_{2}CaCu_{2}O_{8+\delta}$ is measured to be $T_{\theta}=55\rm\,K$, and the corresponding superfluid weight is $D_{s}=e^{2}k_{B}T_{\theta}/\hbar^{2}=1.8\times 10^{9}\,\rm H^{-1}$ Corson et al. (1999). As another example, the superfluid weight of MoGe thin film is measured to be $D_{s}=5\times 10^{8}\,\rm H^{-1}$ Misra et al. (2013).

In Landau-Ginzburg (LG) theory of conventional superconductivity, the superfluid weight is given by $D_{s}\approx e^{2}n_{s}/m^{*}$ in which $m^{*}$ is the band effective mass and $n_{s}$ is the superfluid density Scalapino et al. (1992); Tinkham (2004). At zero temperature, all the electrons have contribution to superconducting transport, which means $n_{s}(T=0)$ is equal to the total electron density, and $n_{s}(T)$ usually decreses with increasing temperature. If the band is exactly flat, the band mass will become infinity, and the LG theory tells us the superfluid weight can be zero even when Cooper pairing happens. We use the Bistritzer-MacDonald model to estimate the bandwidth and the conventional contribution of superfluid weight in TBLG. Around the magic angle, the flat band mostly lie in an energy range $|\varepsilon|<W\approx 0.5\,\rm meV$. Hence the effective mass is approximately $m^{*}\approx\hbar^{2}K_{M}^{2}/2W$, where $K_{M}$ is the distance between $\Gamma$ and $K$ in Moiré BZ. Thus the conventional superfluid weight is $[D_{s}]_{\rm tri}\approx e^{2}n_{s}/m^{*}\approx 2e^{2}WN/\hbar^{2}\Omega_{c}K_{M}^{2}=3\sqrt{3}e^{2}WN/4\pi^{2}\hbar^{2}$, where $\Omega_{c}$ is the area of the Moiré unit cell, and $N$ is the number of electrons per Moiré unit cell. Here we assume that the superfluid density is given by the total electron density, which is the upper limit of $n_{s}$. If we consider the case with filling $\nu=1/4$ or equivalently $N=2$, the value of superfluid weight will be $[D_{s}]_{\rm tri}\approx 5\times 10^{7}\,\rm H^{-1}$, and the corresponding BKT transition temperature will not be higher than $0.6\,\rm K$. However, LG theory is valid only when the band is trivial, as the spreading of its Wannier function has a nonzero lower bound, therefore the estimation based on LG theory in this paragraph is not enough Peotta and Törmä (2015); Tovmasyan et al. (2016); Huber and Altman (2010). As a result, we show that even in the exactly flat band limit, the Cooper pairing may acquire nonlocal phase correlations and thus a nonzero superfluid weight, which - as we show - gives rise to a higher transition temperature.

To obtain the contribution of nontrivial band topology to the superfluid weight, we consider a mean-field Bogoliubov-de-Gennes (BdG) Hamiltonian of a superconductor:

[TABLE]

We use $\Omega_{0}$ to denote the ground state energy of $H_{\rm BdG}$, which is also the free energy at zero temperature. We substitute $\mathbf{k}\rightarrow\mathbf{k}-e\mathbf{A}$ by Peierls substitution Luttinger (1951) when a non-zero uniform gauge potential $\mathbf{A}$ is turned on, and the free energy becomes a function of $\mathbf{A}$. We can then expand $\Omega(\mathbf{A})$ to the second order of $A_{i}$ and obtain $\Omega(\mathbf{A})\approx\Omega_{0}+\frac{1}{2}V[D_{s}]_{ij}A_{i}A_{j}$, where $V$ is the area of the sample, and the second order coefficient $[D_{s}]_{ij}$ is the superfluid weight. The first order derivative $\partial_{\mathbf{A}}\Omega(\mathbf{A})$ gives us the electric current, which agrees with the London equation shown in Eq. (1).

The free energy $\Omega(\mathbf{A})$ and thus the superfluid weight can be derived from the BdG Hamiltonian. The general expression of the superfluid weight is the summation of three terms given by Eqs. (S28), (S29) and (S40) in the supplementary material (SM) Sec. II. The first term in Eq. (S28) corresponds to the Landau-Ginzburg contribution, while Eqs. (S29) and (S40) are additional contributions due to the $\mathbf{k}$ dependence of the flat band Bloch wave functions. When the bands are flat, the conventional contribution vanishes, but the wave function contributions Eqs. (S29) and (S40) can be nonzero, and are related to the band topology as we show below.

Before we start our discussion about TBLG, we briefly review the superconductivity in the spin Chern insulator with exactly flat bands studied in Ref. Peotta and Törmä (2015). In this model, $\mathcal{H}(\mathbf{k})$ has both spinful time reversal symmetry and $s_{z}$ conservation, which allows one to define a spin Chern number $C$. The order parameter $\Delta(\mathbf{k})=is_{y}\Delta$ (in which $s_{y}$ is the $y$ direction spin Pauli matrix) is momentum independent, and one can show that the superfluid weight given by sum of Eqs. (S28), (S29) and (S40) can be reduced to the following integral in the BZ:

[TABLE]

where $\nu$ is the filling ratio of the spinful flat bands, and $g_{ij}(\mathbf{k})$ is the Fubini-Study metric evaluated from the Bloch wave function of the spin $\uparrow$ flat band:

[TABLE]

where $u(\mathbf{k})$ is the Bloch wave function at momentum $\mathbf{k}$ of the spin up flat band. This result is derived in the exact flat band limit, so the contribution from the band dispersion (Eq. (S28)) disappears. Thus we discover that LG theory prediction of superfluid weight is not enough even when the flat band is trivial, because we have models which has trivial bands and $\mathbf{k}$ dependent wave function $u(\mathbf{k})$, such as the “topological” phase of SSH model 111The “topological” phase of SSH model is trivial because it has maximally localized Wannier state, although its Wannier function localization functional is bounded by the winding number from below. A nonzero spin Chern number further gives a nonzero lower bound of the Fubini-Study metric, which will be discussed in the following paragraph.

The Fubini-Study metric defines a distance on the BZ torus: two momentum points are close to each other if their wave functions have a large overlap Cheng (2010). The integral of ${\rm tr}\,g=g_{xx}+g_{yy}$ also corresponds to the gauge invariant part of the “Wannier function localization functional”, or the spread functional, which has been studied in detail in previous research Marzari and Vanderbilt (1997); Marzari et al. (2012). The metric is also related to Berry curvature through the quantum geometric tensor defined by $\mathfrak{G}_{ij}=\partial_{k_{i}}u^{\dagger}(\mathbf{k})(1-u(\mathbf{k})u^{\dagger}(\mathbf{k}))\partial_{k_{j}}u(\mathbf{k})$. The real part of $\mathfrak{G}_{ij}$ is the metric $g_{ij}$ and the imaginary part is the Berry curvature. One of the most important properties of $\mathfrak{G}_{ij}$ is its positive definiteness. It can be shown that for arbitrary complex vectors $\{c_{i}\}$, the inequality $\sum_{ij}c^{\dagger}_{i}\mathfrak{G}_{ij}c_{j}\geq 0$ always holds 222See SM Sec. IV. If we choose $c_{x}=1$ and $c_{y}=i$, we will find ${\rm tr}\,g=g_{xx}+g_{yy}\geq-\mathcal{F}_{xy}$; similarly we choose $c_{x}=1$ and $c_{y}=-i$, and we will obtain ${\rm tr}\,g\geq\mathcal{F}_{xy}$. Therefore we prove that the metric is bounded by the absolute value of curvature ${\rm tr}\,g\geq|\mathcal{F}_{xy}|$. From the expression of $D_{s}$ one can easily notice ${\rm tr}\,D_{s}$ is bounded by the spin Chern number $C$. More details of the quantum metric is discussed in supplementary material Sec. IV. In TBLG, the (spin) Chern number is zero, and the system is multi-band, likely with more complicated pairing symmetry. Hence a new bound/limit (if any exists) for the superfluid weight must be obtained.

We first generalize the result of Ref. Peotta and Törmä (2015) to multi-band systems with a more realistic pairing. The free fermion Hamiltonian $\mathcal{H}(\mathbf{k})$ is assumed to be invariant under spinful time reversal transformation, which is represented by $\mathcal{T}=U_{T}\mathcal{K}$, where $U_{T}$ is a real unitary matrix and $\mathcal{K}$ is complex conjugation operator. We do not (any longer) assume momentum independent pairing. $\mathcal{H}(\mathbf{k})$ is diagonalized by $U(\mathbf{k})$ as $\varepsilon_{\mathbf{k}}=U^{\dagger}(\mathbf{k})\mathcal{H}(\mathbf{k})U(\mathbf{k})$ where $\varepsilon_{\mathbf{k}}$ is a diagonal matrix, and we assume that it has $N_{F}$ flat bands at same energy. We also assume the band gap between these flat bands and any other bands is larger than the bandwidth of flat bands, and the interaction between the electrons. In the following discussion we use a time reversal symmetric pairing between Kramers pairs as follows:

[TABLE]

in which $\Delta_{1,2}\in\mathbb{R}$ because of time reversal symmetry. Here $\tilde{U}_{\mathbf{k}}=(u_{1}(\mathbf{k}),u_{2}(\mathbf{k}),\cdots,u_{N_{F}}(\mathbf{k}))$ is the projection of $U(\mathbf{k})$ into these $N_{F}$ flat bands and $u_{i}(\mathbf{k})$ are the eigenvectors of matrix $\mathcal{H}(\mathbf{k})$. This ansatz implies that $s$ wave pairing happens between Kramers pairs and the pairing strength in the flat bands and in all other bands are given by $\Delta_{1}$ and $\Delta_{2}$, respectively. Therefore we have the following three important assumptions in total: 1) the free fermion Hamiltonian with time reversal symmetry has $N_{F}$ flat bands at the same energy near the Fermi level; 2) there is a large band gap between flat bands and other bands; 3) the pairing order parameter satisfies Eq. (8). Because of this large band gap between flat bands and other bands, we can project the BdG Hamiltonian into the flat bands of the free fermion model and then derive the superfluid weight. The result is given by

[TABLE]

[TABLE]

in which ${\rm Tr}(X)=\sum_{n=1}^{N_{F}}(X_{nn})$ stands for the trace over all the flat band indices. Eq. (10) is the the generalization of Fubini-Study metric in Eq. (Topology-Bounded Superfluid Weight In Twisted Bilayer Graphene) to multi-band systems, which is still positive definite 333See SM Sec. III.. The result in Eq. (6) derived in Ref. Peotta and Törmä (2015) is a special case of our result in Eq. (9) when the time reversal transformation is represented by $U_{T}=is_{y}$, the spin $z$ component is conserved, and $\Delta=\Delta_{1}=\Delta_{2}$.

We also notice that if we assume $\frac{\Delta_{2}}{\Delta_{1}}\leq-1$, the superfluid weight will become zero, or even a negative number. A negative superfluid weight is unphysical, denoting that the BCS wave function of such a pairing is not a stable ground state. This also indicates that if there is no constraint on the order parameter $\Delta(\mathbf{k})$, the superfluid weight will not be bounded. However we expect a weaker pairing strength in the bands which are farther away from the Fermi level, or $|\Delta_{2}|<|\Delta_{1}|$. If the pairing in higher bands are much stronger than pairing in the flat bands - a physically impossible situation -, our projection of the BdG Hamiltonian into the flat bands may also become invalid. Hence we later set $\Delta_{2}=0$ in order to estimate the topological contribution of $D_{s}$.

Now we apply Eq. (9) to TBLG, and show that the fragile topology of TBLG flat bands yields a finite lower bound of the superfluid weight although it has zero spin Chern number. The Bistritzer-MacDonald (BM) model Bistritzer and MacDonald (2011) has $C_{2z}T$, $C_{2x}$ and $C_{3z}$ symmetries, in which $T$ stands for the spinless time reversal transformation. If all the spins and valleys are considered here, we will have well-defined time reversal symmetry, although BM model itself does not. The $C_{2z}T$ symmetry is crucial for the flat bands’ topology Song et al. (2018); Ahn et al. (2018); Zou et al. (2018); Po et al. (2018). Because of this symmetry, the two eigenvalues of the non-Abelian Wilson loop have to be complex conjugation to each other Song et al. (2018); Ahn et al. (2018). A winding number can be defined by Wilson loop eigenvalues. $C_{2z}T$ symmetry gives a constraint not only to the Wilson loop but also to the Berry connection and Berry curvature. It can be shown Song et al. (2018); Ahn et al. (2018) that the non-Abelian Berry connection and curvature of the two flat bands can always be written as $\mathbf{A}(\mathbf{k})=-\mathbf{a}(\mathbf{k})\sigma_{2}$ and $\mathcal{F}_{xy}(\mathbf{k})=-f_{xy}(\mathbf{k})\sigma_{2}$ under a proper local gauge choice on a patch in the Brillouin zone (although a global gauge choice which satisfies this condition might not exist 444See SM Sec. V. The transformation $g_{\mathbf{k}}$ can transform $\mathbf{A}(\mathbf{k})$ into $-\mathbf{a}(\mathbf{k})\sigma_{2}$. However we notice it is not a single value function of $\mathbf{k}$. It changes the sign when it goes around a Dirac point..) In SM Sec. V, we prove that the Wilson loop winding number Song et al. (2018), or the “Euler class” in Ref. Ahn et al. (2018) denoted by $e_{2}$, of the two topological bands, which is an integer (if there are more than two bands, the topological classification will be $\mathbb{Z}_{2}$ instead of $\mathbb{Z}$), is given by the integral of $f_{xy}$ over the whole BZ (with Dirac points removed from the integral area):

[TABLE]

The wave function of the two flat bands in TBLG (per spin per valley - which are good quantum numbers for small twist angles) also can be used to define the positive-definite non-Abelian quantum geometric tensor $\mathfrak{G}_{ij}$ (which is a $2\times 2$ complex matrix) and Fubini-Study metric $g_{ij}=\frac{1}{2}{\rm Tr}\,(\mathfrak{G}_{ij}+\mathfrak{G}_{ij}^{\dagger})$. For arbitrary complex vectors $c_{i}\in\mathbb{C}^{2}$, the inequality $\sum_{ij}c^{\dagger}_{i}\mathfrak{G}_{ij}c_{j}\geq 0$ always holds 555See SM Sec. IV.. By choosing vectors $c_{x}$ and $c_{y}$ properly 666We show how we can choose proper $c_{x}$ and $c_{y}$ in SM Sec. V., we find that the metric is bounded by the “Abelian part” $f_{xy}$ of the non-Abelian Berry curvature $\mathcal{F}_{xy}$: ${\rm tr}\,g\geq 2|f_{xy}|$. The derivation of band topology and metric of the TBLG can be found in SM Sec. V.

In small angle TBLG, all bands are 4-fold degenerate with respect to spin $\uparrow,\downarrow$ and with respect to original graphene valley $K$, $K^{\prime}$. The order parameter ansatz in Eq. (8) corresponds to the pairing between opposite spins and valleys, because time reversal transformation in TBLG will flip both spin $\uparrow,\downarrow$ and valley $K,K^{\prime}$. If we assume no pairing in higher bands ($\Delta_{1}=\Delta,\Delta_{2}=0$), which is physically reasonable, the superfluid weight in the exact flat band limit will be the same as the expression shown in Eq. (6). However, both the pairing strength and metric have different meanings. $g_{ij}(\mathbf{k})$ stands for the Fubini-Study metric derived from the wave functions of two flat bands (of with spin $\uparrow$ and valley $K$, all other degenerate spin and valley bands have identical contribution to the superfluid weight). Also $\Delta_{1}$ is no longer the pairing strength in all the bands but only in these flat bands. Because of the two band winding number protected by $C_{2z}T$ symmetry, we have a new lower bound. The inequality ${\rm tr}\,g>2|f_{xy}|$ naturally leads to the lower bound of the trace of the superfluid weight ${\rm tr}\,D_{s}\geq\frac{8e^{2}\Delta_{1}}{\pi\hbar^{2}}\sqrt{\nu(1-\nu)}$ even though the (spin) Chern number here vanishes. Here we used the fact that the winding number of TBLG flat bands is $e_{2}=1$. As we mentioned earlier, $C_{3z}$ symmetry requires the superfluid weight tensor to be isotropic, therefore $D_{s}=\frac{1}{2}{\rm tr}\,[D_{s}]_{ij}$. Then we can use the following equation to estimate the topological contribution of $D_{s}$:

[TABLE]

From this equation, we find that a nonzero superfluid weight is possible even when the bands are exactly flat, as long as Cooper pairing gap is developed ($\Delta_{1}\neq 0$). The parameter $\Delta_{1}$ can be estimated from the measured $T_{c}$. The filling ratio $\nu$ is determined by the carrier density. Below we use experimental data to estimate the value of superfluid weight at zero temperature.

For flat bands with attractive interaction between Kramers pairs, the mean field BCS wave function is a good approximation to the ground state (SM Sec. VIII). To see this, we can perform a partial particle-hole transformation $c_{\mathbf{k},\uparrow}\rightarrow c_{\mathbf{k}\,\uparrow}\,,c_{\mathbf{k}\downarrow}\rightarrow c^{\dagger}_{-\mathbf{k}\downarrow}$, under which the attractive interaction becomes repulsive, and the BCS wave function becomes a ferromagnetic state. Since repulsive interaction favors ferromagnetic state according to the Hund’s rule, we conclude that BCS wave function describes the ground state well.

As previously mentioned, in 2D superconductors, the transition temperature $T_{c}$ measured in experiment is the Berezinskii-Kosterlitz-Thouless temperature (when phase coherence disappears) instead of the BCS mean field transition temperature $T_{c}^{*}\approx\Delta_{1}/2k_{B}$ (when Cooper pairing vanishes). The BKT transition temperature is given by the universal relation $\frac{\hbar^{2}D_{s}(T_{c})}{e^{2}k_{B}T_{c}}=\frac{8}{\pi}$ Kosterlitz and Thouless (1973), and it is generically lower than the BCS mean field transition temperature $T_{c}^{*}$. To derive the BKT temperature of TBLG in the flat band limit, we can generalize the topological superfluid weight expression in Eq. (12) to finite temperatures, which however has no simple analytical expression. By assuming $\Delta_{1}(T)\approx 2k_{B}T_{c}^{*}(1-T/T_{c}^{*})^{1/2}$ Tinkham (2004), we can numerically calculate the temperature dependence of the superfluid weight. As an example, we have plotted $D_{s}$ in FIG. 1 as a function of $T_{c}/T_{c}^{*}$ for filling ratio $\nu=1/4$ (2 electrons per Moiré unit cell), from which we find $T_{c}/T_{c}^{*}=0.35$. See SM Sec. VII for more detailed calculations.

We now estimate the TBLG topological superfluid weight $D_{s}$ at zero temperature. When the bandwidth is small (zero), one expects the superfluid weight to be dominated by the band topology contribution. The experimental transition temperature Cao et al. (2018); Yankowitz et al. (2019) is $T_{c}=1.5\rm\,K$, thus the order parameter can be estimated to be $\Delta_{1}=2k_{B}T_{c}^{*}\approx 0.74\rm\,meV$. By using $\nu=1/4$, the topological superfluid weight is $[D_{s}]_{\rm top}\approx\frac{4e^{2}\Delta_{1}}{\pi\hbar^{2}}\sqrt{\nu(1-\nu)}\approx 1.5\times 10^{8}\rm\,H^{-1}$. One notices that this is an order of magnitude smaller than the superfluid weight in conventional materials, such as BSCCO and MoGe Emery and Kivelson (1995); Basov and Chubukov (2011). However in both cuprates and MoGe thin films, the typical carrier density is around $n\approx 10^{14}\,\rm cm^{-2}$ Uemura et al. (1989) two orders of magnitude larger than that in TBLG (where $n\approx 10^{12}\,\rm cm^{-2}$), hence such a superfluid weight is already large for TBLG.

Moving away from completely flat bands, we find that the conventional term in the superfluid weight mostly depends on the bandwidth $W$, while the topological term mostly depends on the pairing order parameter $\Delta_{1}$ (or the transition temperature $T_{c}$). In TBLG, the bandwidth $W$ and the transition temperature $T_{c}$ have a similar magnitude, and we expect the topological term to have an important contribution. In the strong pairing limit where $\Delta>W$, the superfluid weight will be underestimated if by only the conventional contribution Stauber et al. (2018).

Recently, superconductivity has been observed in other Moiré systems, including twisted double bilayer graphene, and multilayer graphene/boron nitride heterostructure Zhang et al. (2018); Zhang and Senthil (2018); Chen et al. (2018); Moriyama et al. (2019); Chen et al. (2019); Chittari et al. (2019); Shen et al. (2019); Liu et al. (2019); Endo et al. (2019). In these systems, a displacement electric field is necessary for superconductivity. This yields gapped flat bands with a nonzero valley Chern number Zhang et al. (2018); Chittari et al. (2019). The Chern number can be larger than one, which could lead to a larger topological lower bound of superfluid weight. In twisted multilayer graphene, the bandwidth of flat bands is even smaller than that in TBLG. A higher transition temperature (larger $\Delta$) was also observed Shen et al. (2019); Liu et al. (2019). We expect the topological lower bound plays a significant role in the superfluid weight of these systems. Moreover, multi-layer systems with $C_{2z}T$ symmetry can realize larger Euler $e_{2}$ invariants (higher Wilson loop windings), hence increasing the superfluid weight.

In summary, we proved that fragile topology can yield a nonzero superfluid weight (phase stiffness) of a superconductor even when the bands are exactly flat. In TBLG, by assuming exact flat bands and uniform $s$ wave pairing, the superfluid weight can be written as the integral of Fubini-Study metric, which is bounded by the Wilson loop winding number of the two flat bands of TBLG. Based on the continuum model and the superconducting transition temperature measured in TBLG experiments, we estimated the magnitude of the superfluid weight in TBLG, and show that the topological contribution may play a significant role. The same argument can be applied to other Moiré superlattices with flat bands and non-trivial band topology. This shows that topological bands can have $T_{c}$ much higher than their band widths. We also note that our lower bound is consistent with the upper bound studied in Ref. Hazra et al. (2018). For topological flat bands, the upper bound in Ref. Hazra et al. (2018) is around the energy gap $E_{g}$ between the flat bands and other bands, while our lower bound is proportional to the order parameter $\Delta$, which is derived under the assumption that $|\Delta|\ll E_{g}$.

Acknowledgments. We acknowledge helpful discussion with Nai Phuan Ong, Ali Yazdani and Guorui Chen. ZS and BB are supported by the Department of Energy Grant No. DE-SC0016239, the National Science Foundation EAGER Grant No. DMR 1643312, Simons Investigator Grants No. 404513, ONR No. N00014-14-1-0330, NSF-MRSEC No. DMR-142051, the Packard Foundation, the Schmidt Fund for Innovative Research. BL is supported by Princeton Center for Theoretical Science at Princeton University.

Supplementary Material

The supplementary material sections are organized as follows:

• •

Sec. I: MEAN FIELD THEORY, BOGOLIUBOV QUASIPARTICLES, AND THE ZERO TEMPERATURE GRAND POTENTIAL

• •

Sec. II: SUPERFLUID WEIGHT UNDER MEAN-FIELD APPROXIMATION

• •

Sec. III: SUPERFLUID WEIGHT WITHIN FLAT BANDS AND RELATION TO THE QUANTUM METRIC

• •

Sec. IV: QUANTUM GEOMETRIC TENSOR, FUBINI-STUDY METRIC AND BERRY CURVATURE: LOWER BOUNDS

• •

Sec. V: BOUNDS ON THE QUANTUM METRIC FOR DIFFERENT SYSTEMS: SSH CHAIN, CHERN INSULATOR, AND $C_{2z}T$ FRAGILE TOPOLOGY

• •

Sec. VI: SUPERCONDUCTIVITY IN TBLG

• •

Sec. VII: FINITE TEMPERATURE STUDY AND BEREZINSKII-KOSTERLITZ-THOULESS TRANSITION

• •

Sec. VIII: THE VALIDITY OF BCS WAVE FUNCTION IN THE FLAT BAND LIMIT

In the case without specification, we set $\hbar=e=k_{B}=1$ in supplementary material for convenience.

I Mean field theory, Bogoliubov quasiparticles and the zero temperature free energy

In this section, we review the BdG mean field Hamiltonian of a superconductor in the presence of an external gauge field. Before we consider the interaction between electrons, we write down the free fermion Hamiltonian:

[TABLE]

where $\psi_{\mathbf{k}}$ is a spinor with $N_{B}$ components labeled by $\alpha$. Here $\alpha$ can stand for the local orbitals in tight binding models, or reciprocal lattices in continuum models, such as the $\mathbf{Q}$ lattice in TBLG Bistritzer and MacDonald (2011); Song et al. (2018). The spin index has also been included in $\alpha$. The Bloch Hamiltonian $\mathcal{H}(\mathbf{k})$ can be diagonalized by a unitary matrix $U(\mathbf{k})$ as $\varepsilon_{\mathbf{k}}=U^{\dagger}(\mathbf{k})\mathcal{H}(\mathbf{k})U(\mathbf{k})$. If a time dependent uniform gauge field $\mathbf{q}=e\mathbf{A}$ - which corresponds to a uniform electric field -, is applied, we can obtain the Hamiltonian with Peierls substitution by shifting $\mathbf{k}$ to $\mathbf{k}-\mathbf{q}$:

[TABLE]

The attractive interaction between the electrons leads to pairing terms such as $\psi^{\dagger}\psi^{\dagger}$ appear in the Hamiltonian under mean-field approximation. We consider the case without FFLO order, hence only the pairing between electrons with momentum $\mathbf{k}$ and $-\mathbf{k}$ exists. In general, the pairing term takes the form:

[TABLE]

where $\Delta(\mathbf{k})=-\Delta^{\rm T}(-\mathbf{k})$ because of particle hole redundancy. Since we have used Peierls substitution, the effect of gauge field has been included in kinetic term, and the pairing term will not depend on $\mathbf{q}$. The Hamiltonian can be written as a bilinear form of a Nambu spinor $\Psi_{\mathbf{k}}=\left(\psi_{\mathbf{k}1},\psi_{\mathbf{k}2},\cdots,\psi_{-\mathbf{k}1}^{\dagger},\psi_{-\mathbf{k}2}^{\dagger},\cdots\right)^{\rm T}$:

[TABLE]

If the free fermion Hamiltonian can be diagonalized by matrix $U(\mathbf{k})$ as $\varepsilon_{\mathbf{k}}=U^{\dagger}(\mathbf{k})\mathcal{H}(\mathbf{k})U(\mathbf{k})$ where $\varepsilon_{\mathbf{k}}$ is a diagonal matrix, we will be able to use it to diagonalize the diagonal blocks of Eq. (S6):

[TABLE]

By diagonalizing the Bogoliubov-de-Gennes Hamiltonian $\mathscr{H}_{\mathbf{k}}(\mathbf{q})$, the quasiparticle spectrum of the superconductor can be obtained. Assume the matrix $\mathscr{H}_{\mathbf{k}}(\mathbf{q})$ is diagonalized by unitary matrix $W_{\mathbf{k}}(\mathbf{q})$, we will have the following equation:

[TABLE]

where $E_{\mathbf{k}}(\mathbf{q})$ is a diagonal matrix and its eigenvalues $E_{\mathbf{k}n}(\mathbf{q})$ give us the quasiparticle spectrum. When $\mathbf{q}=0$, the matrix $W(0)$ is the Bogoliubov quasiparticle coefficients:

[TABLE]

and because of the unitarity of matrix $W_{\mathbf{k}}(0)$, the blocks satisfy $\mathcal{U}_{\mathbf{k}}\mathcal{U}^{\dagger}_{\mathbf{k}}+\mathcal{V}_{\mathbf{k}}\mathcal{V}^{\dagger}_{\mathbf{k}}=1,\leavevmode\nobreak\ \mathcal{U}_{\mathbf{k}}\mathcal{V}^{\rm T}_{-\mathbf{k}}+\mathcal{V}_{\mathbf{k}}\mathcal{U}^{\rm T}_{-\mathbf{k}}=0$. These matrices will be helpful for us to simplify the expression of superfluid weight later. At zero temperature, the free energy is given by the summation of all the negative eigenvalues of $\mathscr{H}_{\mathbf{k}}(\mathbf{q})$ at every momentum $\mathbf{k}$ in the BZ. More precisely, the free energy at zero temperature can be expressed as

[TABLE]

where $\rm Tr(\cdots)$ represents the trace over the energy bands. Because of the particle hole redundancy, the eigenvalues of matrix $\mathscr{H}_{\mathbf{k}}(0)$ are related with each other by $E_{\mathbf{k}(N_{B}+n)}(0)=-E_{-\mathbf{k}n}(0)<0$ when $1\leq n\leq N_{B}$ in which $N_{B}$ is the number of free electron bands.

II Superfluid weight under mean-field approximation

For completeness, we reproduce the derivation of the superfluid weight. Our procedure is extremely similar to the derivation in Ref. Peotta and Törmä (2015). For simplicity, the off-diagonal block of $\mathscr{H}_{\mathbf{k}}(\mathbf{q})$ is denoted by $\mathcal{D}_{\mathbf{k}}(\mathbf{q})$:

[TABLE]

where $U(\mathbf{k})$ is the unitary matrix which diagonalize the original noninteracting Hamiltonian. $D_{s}$ is defined to be the second order expansion coefficients of free energy, so it can be expressed as the second order derivative of $\Omega(\mathbf{q})$:

[TABLE]

The first order derivative of the quasiparticle energy $E_{\mathbf{k}n}(\mathbf{q})$ is given by:

[TABLE]

here we use Feynman-Hellman theorem to derive the second equality. Further we take the second order derivative of Eq. (II), and it becomes:

[TABLE]

This expression can be further simplified. Suppose $W_{\mathbf{k}}(\mathbf{q})_{\alpha,n}$ is the $n$-th eigenvector of matrix $\mathscr{H}_{\mathbf{k}}(\mathbf{q})_{\alpha,\beta}$ where $\mathbf{q}$ is a parameter, then by definition we have $\sum_{\beta}\mathscr{H}_{\mathbf{k}}(\mathbf{q})_{\alpha,\beta}W_{\mathbf{k}}(\mathbf{q})_{\beta,n}=E_{\mathbf{k}n}(\mathbf{q})W_{\mathbf{k}}(\mathbf{q})_{\alpha,n}$. Because $\mathscr{H}_{\mathbf{k}}(\mathbf{q})$ is a hermitian matrix, we will have the following equation when $n\neq m$:

[TABLE]

Now we take derivative $\partial_{q_{i}}$ on both sides of this equation, and we obtain (for $n\neq m$)

[TABLE]

The unitarity of $W_{\mathbf{k}}(\mathbf{q})$ is used in the third equality. This result can also be written as

[TABLE]

or equivalently

[TABLE]

This theorem (which is nothing else than the classic expression of the non-abelian Berry phase in terms of a sum over eigenstates) can be used to simplify the second and third terms in Eq. (S23):

[TABLE]

and now Eq. (S23) becomes:

[TABLE]

We have hence removed the derivatives of $W_{\mathbf{k}}(\mathbf{q})$, which in principle are not immediately obtainable. Only the derivatives of $\mathscr{H}_{\mathbf{k}}(\mathbf{q})$ are left in the expressions, and these are immediately obtained. The first order derivative of $\mathscr{H}_{\mathbf{k}}(\mathbf{q})$ is given by:

[TABLE]

Similarly the second order derivative of matrix $\mathscr{H}_{\mathbf{k}}(\mathbf{q})$ will be:

[TABLE]

We substitute Eqs. (S26) and (S27) into Eq. (S25) to get the result of superfluid weight. Now only the Bogoliubov coefficients $W_{\mathbf{k}}(0)$ are necessary. By using Eq. (S26) and the first term of Eq. (S25), we obtain the following two terms of superfluid weight:

[TABLE]

The contribution from $\frac{1}{2}\sum_{\mathbf{k}}\partial_{k_{i}}\partial_{k_{j}}\varepsilon_{\mathbf{k}}$ has been included in $[D_{s}]^{(1)}$. Similarly, by substituting Eq. (S24) into Eq. (S25), we get the third term of the superfluid weight:

[TABLE]

We use Eqs. (S17) and (S26) to express $W^{\dagger}_{\mathbf{k}}\partial_{q_{i}}\mathscr{H}_{\mathbf{k}}(\mathbf{q})W_{\mathbf{k}}|_{\mathbf{q}=0}$ by the following block matrix:

[TABLE]

We use these blocks to express the third term of superfluid weight (Eq. (S30)), and it can be written as

[TABLE]

First we study the contribution from these $P_{i}(\mathbf{k})$ matrices. The first term in Eq. (S37) is given by:

[TABLE]

and we notice it has two dummy indices $n$ and $n^{\prime}$. If we simply switch them ($n\leftrightarrow n^{\prime}$), it will become

[TABLE]

and it differs from the second term in Eq. (S37) only by a minus sign. Therefore the summation of the first and second terms will be zero, and all $P_{i}(\mathbf{k})$ matrices will not appear in the final expression. Similarly, we can prove $R_{i}(\mathbf{k})$ will also disappear. Only matrices $Q_{i}(\mathbf{k})$ will be left. By using the particle hole symmetry $E_{\mathbf{k}n}(0)=-E_{-\mathbf{k}(n+N_{B})}$, $[D_{s}]^{(3)}_{ij}$ can be expressed as:

[TABLE]

By adding Eqs. (S28), (S29) and (S40) together, we finally get the general expression of superfluid weight. This result is accurate under mean-field approximation without FFLO order. In the following section we will see that the superfluid weight can be related to the quantum geometric tensor under the flat band approximation. Our result here is very similar to the result shown in Ref. Peotta and Törmä (2015), but it is more general. We used the Nambu spinor with particle-hole redundancy, and our result can also be applied to Hamiltonians with Rashba coupling.

III Superfluid weight within flat bands and relation to the quantum metric

In this section we study the superfluid weight in the flat band limit. The following discussion is based on these assumptions: 1) the free fermion Hamiltonian $\mathcal{H}(\mathbf{k})$ has $N_{F}$ flat bands, with energy $\varepsilon_{0}$, and it can be diagonalized by unitary matrix $U(\mathbf{k})$; 2) the band gap between flat bands and other bands is larger than the band width of flat bands (when the bands are exactly flat, this must be true) and the pairing order parameter. Because of the large band gap, the derivation of free energy and superfluid weight can be simplified dramatically by projecting into the flat bands.

Now we discuss the form of pairing order parameter. It is reasonable to first assume that $\Delta(\mathbf{k})$ obeys spinful time reversal symmetry $\mathcal{T}$, as assumed in Ref. Lian et al. (2018a). This may seem irrelevant with TBLG since BM model does not have time reversal. However if we take spin $\uparrow,\downarrow$ and valley $K,K^{\prime}$ into consideration, the whole system will have time reversal symmetry. The time reversal transformation can always be represented by $U_{T}\mathcal{K}$ in which $\mathcal{K}$ is the complex conjugation, and $U_{T}$ is an real antisymmetric unitary matrix acting on the orbital and spin indices. The fermion operators transform under $\mathcal{T}$ as $\mathcal{T}\psi_{\mathbf{k}i}\mathcal{T}^{-1}=(U_{T})_{ij}\psi_{-\mathbf{k}j}$. If the BdG Hamiltonian has spinful time reversal symmetry, the pairing order parameter will satisfy $U_{T}\Delta^{\dagger}(\mathbf{k})U_{T}=\Delta(\mathbf{k})$. For simplicity we assume the order parameter has the following form:

[TABLE]

where $\tilde{U}_{\mathbf{k}}=(u_{1}(\mathbf{k}),u_{2}(\mathbf{k}),\cdots,u_{N_{F}}(\mathbf{k}))$ stands for the projection of $U(\mathbf{k})$ into flat bands, thus $\tilde{U}_{\mathbf{k}}\tilde{U}_{\mathbf{k}}^{\dagger}$ is projection operator of flat bands, $\left(\mathds{1}-\tilde{U}_{\mathbf{k}}\tilde{U}_{\mathbf{k}}^{\dagger}\right)$ is projector of all the other bands, and $\tilde{U}^{\dagger}_{\mathbf{k}}\tilde{U}_{\mathbf{k}}=\mathds{1}_{N_{F}\times N_{F}}$. Cooper pairing happens between the Kramers pairs, and the pairing strength in flat bands is $\Delta_{1}$, and in other bands it is $\Delta_{2}$. If the BdG Hamiltonian is assumed to be time reversal invariant, then as we mentioned in last paragraph, equation $U_{T}\Delta^{\dagger}(\mathbf{k})U_{T}=\Delta(\mathbf{k})$ must be satisfied. This equation requires $\Delta_{1,2}$ to be real. To show this, we have the following equation

[TABLE]

thus when $\Delta_{1,2}\in\mathbb{R}$, the order parameter $\Delta(\mathbf{k})$ satisfies time reversal symmetry.

As we mentioned earlier, because the band gap between flat bands and other bands is larger than other energy scales in flat bands, it is not necessary to use the tedious method in Sec. II to obtain the free energy and superfluid weight. Therefore, the derivation becomes much simpler if we project the Hamiltonian into these flat bands. The BdG Hamiltonian in Eq. (S13) includes all the energy bands. If we only keep the columns and rows which corresponds to the flat bands, it will become:

[TABLE]

where $\varepsilon_{0}$ is the flat band energy, and $w(\mathbf{k})=\tilde{U}^{\dagger}(\mathbf{k})U_{T}\tilde{U}^{*}(-\mathbf{k})$ is the time reversal sewing matrix of all the flat bands. In order to evaluate the free energy, we have to know the eigenvalues of $\tilde{\mathscr{H}}_{\mathbf{k}}(\mathbf{q})$. Thus we calculate the square of the projected BdG Hamiltonian:

[TABLE]

Note that it is block diagonal. We use $\lambda_{\mathbf{k}n}(\mathbf{q})$ and $\varphi_{\mathbf{k}n}(\mathbf{q})$ to denote the eigenvalues of matrices $\tilde{\mathcal{D}}_{\mathbf{k}}(\mathbf{q})\tilde{\mathcal{D}}^{\dagger}_{\mathbf{k}}(\mathbf{q})$ and $\tilde{\mathcal{D}}^{\dagger}_{\mathbf{k}}(\mathbf{q})\tilde{\mathcal{D}}_{\mathbf{k}}(\mathbf{q})$, respectively. Obviously, both $\tilde{\mathcal{D}}_{\mathbf{k}}(\mathbf{q})\tilde{\mathcal{D}}^{\dagger}_{\mathbf{k}}(\mathbf{q})$ and $\tilde{\mathcal{D}}^{\dagger}_{\mathbf{k}}(\mathbf{q})\tilde{\mathcal{D}}_{\mathbf{k}}(\mathbf{q})$ are Hermitian and semi-positive definite, we find that $\lambda_{\mathbf{k}n}(\mathbf{q})\geq 0$ and $\varphi_{\mathbf{k}n}(\mathbf{q})\geq 0$. Therefore, the eigenvalues of the projected BdG Hamiltonian satisfy $E^{2}_{\mathbf{k}n}(\mathbf{q})=(\varepsilon_{0}-\mu)^{2}+\lambda_{\mathbf{k}n}(\mathbf{q})$ when $1\leq n\leq N_{F}$ and $E^{2}_{\mathbf{k}n}(\mathbf{q})=(\varepsilon_{0}-\mu)^{2}+\varphi_{\mathbf{k}n}(\mathbf{q})$ when $N_{F}+1\leq n\leq 2N_{F}$. The free energy can be written can be written as the following form:

[TABLE]

Then we take the second order derivative of $\Omega(\mathbf{q})$ to get the superfluid weight:

[TABLE]

When $\mathbf{q}=0$, the matrices $\tilde{\mathcal{D}}_{\mathbf{k}}(0)$, $\tilde{\mathcal{D}}_{\mathbf{k}}(0)\tilde{\mathcal{D}}_{\mathbf{k}}(0)^{\dagger}$ and $\tilde{\mathcal{D}}_{\mathbf{k}}(0)^{\dagger}\tilde{\mathcal{D}}_{\mathbf{k}}(0)$ will become

[TABLE]

[TABLE]

Thus we must have $\lambda_{\mathbf{k}n}(0)=\varphi_{\mathbf{k}n}(0)=\Delta_{1}^{2}$, and the denominator in Eq. (S58) will hence become a constant $E_{B}=\sqrt{(\varepsilon_{0}-\mu)^{2}+\Delta_{1}^{2}}$, and the superfluid weight can be written as

[TABLE]

The first order derivative of matrices $\tilde{\mathcal{D}}_{\mathbf{k}}(\mathbf{q})\tilde{\mathcal{D}}^{\dagger}_{\mathbf{k}}(\mathbf{q})$ and $\tilde{\mathcal{D}}^{\dagger}_{\mathbf{k}}(\mathbf{q})\tilde{\mathcal{D}}_{\mathbf{k}}(\mathbf{q})$ are all equal to zero. In order to show this, we can calculate the first order derivative of $\mathcal{D}_{\mathbf{k}}(\mathbf{q})$ and $\mathcal{D}^{\dagger}_{\mathbf{k}}(\mathbf{q})$:

[TABLE]

The first order derivative of the product $\tilde{\mathcal{D}}_{\mathbf{k}}(\mathbf{q})\tilde{\mathcal{D}}^{\dagger}_{\mathbf{k}}(\mathbf{q})$ and $\tilde{\mathcal{D}}^{\dagger}_{\mathbf{k}}(\mathbf{q})\tilde{\mathcal{D}}_{\mathbf{k}}(\mathbf{q})$ will become

[TABLE]

[TABLE]

Then by using Feynman-Hellman’s theorem, we have $\partial_{q_{i}}\lambda_{\mathbf{k}n}(\mathbf{q})|_{\mathbf{q}=0}=\partial_{q_{i}}\varphi_{\mathbf{k}n}(\mathbf{q})|_{\mathbf{q}=0}=0$. Hence the superfluid weight can be further simplified as shown:

[TABLE]

where the integrand of this result is the Fubini-Study metric mentioned in the main text, and it will be discussed in detail in Sec. IV. Because of the pairing term, the chemical potential is no longer at the flat band energy $\varepsilon_{0}$. In fact, it can be determined by the filling ratio $\nu$. Bogoliubov theory tells us that the filling ratio of momentum $\mathbf{k}$ is given by $|v_{\mathbf{k}}|^{2}$:

[TABLE]

and in the flat band limit, $\varepsilon_{\mathbf{k}}=\varepsilon_{0}$ is momentum independent, which leads to a $\mathbf{k}$ independent $v_{\mathbf{k}}$. Therefore, the filling ratio of the flat bands is $\nu=v^{2}_{\mathbf{k}}$. By using this relation, the superfluid weight at $T=0$ can be written as

[TABLE]

In fact the integrand in Eq. (S60) is the Fubini-Study metric. It will be discussed in the following section.

IV quantum geometric tensor, fubini-study metric and berry curvature: lower bounds

Fubini-Study metric mentioned in last section is closely related with Berry curvature through the “quantum geometric tensor”. Suppose we have $n$ orthonormal vectors $u_{m}(\mathbf{k})$, $m=1\ldots n$, in a $N$ dimensional Hilbert space, where $\mathbf{k}$ is some parameter. The quantum geometric tensor can be defined as:

[TABLE]

in which $m,n$ are energy band indices and $i,j$ are spatial direction indices. For convenience we denote $\tilde{u}_{\mathbf{k}}=(u_{n_{1}}(\mathbf{k}),u_{n_{2}}(\mathbf{k}),\cdots,u_{n_{n}}(\mathbf{k}))$ where $n_{1},n_{2},\cdots,n_{n}$ are the indices of the energy bands we are interested in. By using $\tilde{u}_{\mathbf{k}}$, the quantum geometric tensor can be expressed in a more compact form:

[TABLE]

The Hermitian (real) and anti-Hermitian (imaginary) parts of $\mathfrak{G}_{ij}$ are:

[TABLE]

For convenience we use $\mathfrak{g}_{ij}$ to denote ${\rm Re}\,\mathfrak{G}_{ij}$, and $g_{ij}={\rm Tr}\,\mathfrak{g}_{ij}$ is the “Fubini-Study metric” mentioned earlier. It is easy to notice that $g_{ij}$ is actually the integrand appeared in Eqs. (S60). If the Berry connection is defined by $\mathbf{A}=i\tilde{u}^{\dagger}_{\mathbf{k}}\partial_{\mathbf{k}}\tilde{u}_{\mathbf{k}}$, then the imaginary part of $\mathfrak{G}_{ij}$ is proportional to the Berry curvature: ${\rm Im}\,\mathfrak{G}_{ij}=-\frac{1}{2}\mathcal{F}_{ij}=-\frac{1}{2}(\partial_{k_{i}}A_{j}-\partial_{k_{j}}A_{i}-i[A_{i},A_{j}])$. So in conclusion, the quantum geometric tensor can be written as

[TABLE]

By this quantity, the Fubini-Study metric and Berry curvature are related with each other. An important feature of $\mathfrak{G}_{ij}$ is its positive definiteness. Suppose we have several complex vectors $c_{i}\in\mathbb{C}^{n}$, and we will get the following equation:

[TABLE]

Because $\{u_{m}(\mathbf{k})\}$ are orthonormal vectors, the matrix $\left(\mathds{1}-\tilde{u}_{\mathbf{k}}\tilde{u}^{\dagger}_{\mathbf{k}}\right)$ is a projector, and the eigenvalues can only be [math] or $1$. Therefore a scalar product $\varphi^{\dagger}\left(\mathds{1}-\tilde{u}_{\mathbf{k}}\tilde{u}^{\dagger}_{\mathbf{k}}\right)\varphi$ must be non-negative. If we choose the complex vectors $c_{i}$ properly, Eq. (S66) can be used to prove important inequalities between the Fubini-Study metric and Berry curvature, as we will show in Sec. V.

The geometric meaning of Fubini-Study metric can be understood as distance between quantum states. In fact the Bloch wave functions of $N$ bands $\tilde{u}_{\mathbf{k}}$ define a map from the Brillouin zone torus to $\mathbb{CP}^{N-1}$. If we define the distance between two points $\mathbf{k}$ and $\mathbf{k}+d\mathbf{k}$ as shown:

[TABLE]

then by expanding this equation to the second order we will find $d^{2}(\mathbf{k},\mathbf{k}+d\mathbf{k})=g_{ij}(\mathbf{k}_{1})dk_{i}dk_{j}$. A pedagogical introduction to quantum distance and quantum geometric tensor can be found in Ref. Cheng (2010).

A physical interpretation of the Fubini-Study metric is related with the Wannier function localization which is studied in Ref. Marzari et al. (2012); Marzari and Vanderbilt (1997). If the Bloch wave function of a state in band $n$ and momentum $\mathbf{k}$ is $|\psi_{\mathbf{k}n}$, then the Wannier states can be obtained by the discrete Fourier transformation

[TABLE]

The Wannier function localization functional can be defined as follows:

[TABLE]

where $\hat{\mathbf{r}}$ is the position operator. However this functional is represented in the Wannier basis. We will show how to express it under Bloch basis. It is well known that the Bloch wave function can be expressed as the product of a periodic function $u_{n\mathbf{k}}$ (the Bloch function) and a phase $e^{i\mathbf{k}\cdot\hat{\mathbf{r}}}$:

[TABLE]

The overlap between two Bloch functions with different momentum will be:

[TABLE]

The right hand side of this equation is the Bloch states, and we can transform it into Wannier states as shown:

[TABLE]

Because the Wannier functions have a discrete translation symmetry along Bravias lattice, we obtain

[TABLE]

Eq. (S70) now becomes

[TABLE]

Up till now we were using the periodic Bloch wave functions $|u_{\mathbf{k}n}\rangle$ instead of the eigenvectors $u_{m}(\mathbf{k})$ of matrix $\mathcal{H}(\mathbf{k})$. In fact it can be shown that for both continuum models and tight binding models with exactly localized atomic orbitals, we have $\langle u_{m\mathbf{k}}|u_{n\mathbf{k}+\mathbf{q}}\rangle=u^{\dagger}_{m}(\mathbf{k})u_{n}(\mathbf{k}+\mathbf{q})$.

The Bloch wave function is given by the following equation in both tight binding models and continuum models:

[TABLE]

In continuum models, the states $|\mathbf{k},\alpha\rangle$ are given by the plain waves:

[TABLE]

where $\mathbf{Q}_{\alpha}$ is a reciprocal vector in momentum space. By using the plain wave basis, we have

[TABLE]

In tight binding models, the states $|\mathbf{k},\alpha\rangle$ are given by

[TABLE]

where $w_{\alpha}(\mathbf{r}-\mathbf{R}-\tau_{\alpha})$ is the local wave function in the unit cell $\mathbf{R}$ and orbital $\alpha$. The inner product of periodic Bloch wave function with different momenta will be

[TABLE]

If the wave functions of atomic orbitals are assumed to be exactly localized as $w^{*}_{\alpha}(\mathbf{r}-\mathbf{R}_{1}-\tau_{\alpha})w_{\beta}(\mathbf{r}-\mathbf{R}_{2}-\tau_{\beta})=\delta(\mathbf{r}-\mathbf{R}_{1}-\tau_{\alpha})\delta_{\mathbf{R}_{1},\mathbf{R}_{2}}\delta_{\alpha,\beta}$, we will obtain the following equation

[TABLE]

Now we have proved that for both continuum models and tight binding models, the inner product of periodic Bloch wave functions $|u_{n\mathbf{k}}\rangle$ is equal to the inner product of the eigenvectors $u_{n}(\mathbf{k})$ of $\mathcal{H}(\mathbf{k})$. Therefore, Eq. (S70) can be written as

[TABLE]

Then we take the first and second order derivatives of $\mathbf{q}$ on both sides of this equation, and evaluate the result at $q=0$ to obtain:

[TABLE]

The Fourier transformation of these two equations will be

[TABLE]

where $\Omega_{c}$ is the area of the unit cell. For future convenience we divide the localization functional $F$ into the following two parts:

[TABLE]

and let us denote the first and second terms in Eq. (S82) by $F_{I}$ and $\tilde{F}$, respectively. Also $\tilde{F}$ is always positive because of its definition. By using Eqs. (S80) and (S81), the first part $F_{I}$ can be written as the following form:

[TABLE]

The first term is always a real number because it is equal to $\langle 0n|\hat{\mathbf{r}}^{2}|0n\rangle$. So by taking the complex conjugation we obtain

[TABLE]

Also from the normalization condition we have $u^{\dagger}_{n}(\mathbf{k})u_{n}(\mathbf{k})=1$. The second order derivative tells us $u^{\dagger}_{n}(\mathbf{k})\nabla^{2}_{\mathbf{k}}u_{n}(\mathbf{k})+\nabla^{2}_{\mathbf{k}}u^{\dagger}_{n}(\mathbf{k})u_{n}(\mathbf{k})+2\nabla_{\mathbf{k}}u^{\dagger}_{n}(\mathbf{k})\nabla_{\mathbf{k}}u_{n}(\mathbf{k})=0$. Consequently we will obtain

[TABLE]

We use this expression to replace the first term appeared in $F_{I}$, and the result is

[TABLE]

Notice that the integrand is the trace of the Fubini-Study metric ${\rm tr}\,g$ defined by the eigenvectors $u_{n}(\mathbf{k})$. Since the metric is invariant under gauge transformation, $F_{I}$ is gauge invariant. Hence if the integral of ${\rm tr}\,g$ has a nonzero lower bound, the gauge invariant part of Wannier function localization functional will also be bounded. Because $\tilde{F}$ is always positive by definition, the lower bound of the gauge invariant part $F_{I}$ is also the lower bound of the functional $F$ itself. In the following section we will show some examples with bounded Fubini-Study metric.

V bounds on the quantum metric for different systems: SSH chain, Chern insulator and $C_{2z}T$ fragile topology

V.1 SSH Chain

Several examples of metric lower bound are studied in this section. We start our discussion with SSH model. SSH model is a one dimensional model with chiral symmetry and it can be classified by a winding number. In this paragraph we show that the integral of Fubini-Study metric of SSH chain is bounded by the winding number. If the chiral symmetry is represented by $\sigma_{z}$, then the (flat band) Hamiltonian of SSH model can be written as

[TABLE]

in which $\phi(k)$ is a function of $k$. The topological index of this model is the winding number, defined by $\mathcal{W}=\frac{1}{2\pi}\int dk\,\partial_{k}\phi(k)\in\mathbb{Z}$. By diagonalizing $\mathcal{H}(k)$ and finding the wave function of valence band states, we can find the Fubini-Study metric of SSH model is:

[TABLE]

then it is easy to show that the integral of metric is bounded by the winding number:

[TABLE]

in which $a$ is the lattice constant. Hence the integral of Fubini-Study metric in one dimension with chiral symmetry is bounded by the winding number.

V.2 Chern Insulator

The next example will be a band with nonzero Chern number. The integral of Berry curvature of a Chern band must be $2\pi n,n\in\mathbb{Z}$. We consider the case with just one band with a Chern number $C_{1}$. That means $\mathfrak{G}_{ij}=g_{ij}-\frac{i}{2}\mathcal{F}_{xy}$ is a complex number instead of a matrix and $\mathfrak{g}_{ij}=g_{ij}$. If we choose $c_{x}=1,c_{y}=i$, Eq. (S66) will become:

[TABLE]

Similarly, if we choose $c_{x}=1,c_{y}=-i$, we will get ${\rm tr}\,g\geq\mathcal{F}_{xy}$. In conclusion, the trace of the metric is bounded by the local Berry curvature ${\rm tr}\,g\geq|\mathcal{F}_{xy}|$. Obviously, lower bound for the integral of metric is given by the Chern number:

[TABLE]

V.3 Topology of TBLG flat bands

The flat bands in TBLG are not Chern bands. In fact, their Chern number vanishes. The effective continuum model of TBLG per spin per valley is the Bistritzer-MacDonald model Bistritzer and MacDonald (2011) whose two flat bands around the charge neutral point have a nontrivial topology protected by this $C_{2z}T$ symmetry. The free fermion Hamiltonian of valley $K$ and spin up is given by:

[TABLE]

in which $\mathbf{Q},\mathbf{Q}^{\prime}$ take their value in the hexagonal lattice formed by adding $\mathbf{q}_{1,2,3}$ iteratively. Notice that this Hamiltonian is not periodic in momentum space. In fact it changes as $h(\mathbf{k}+\mathbf{G})=V_{\mathbf{G}}h(\mathbf{k})V^{\dagger}_{\mathbf{G}}$ under momentum translation, in which $(V_{\mathbf{G}})_{\mathbf{Q},\mathbf{Q}^{\prime}}=e^{i\alpha}\delta_{\mathbf{Q}^{\prime}-\mathbf{Q},\mathbf{G}}$. The choice of $V_{\mathbf{G}}$ has a phase degree of freedom, and we choose it to be real ($\alpha=0,\pi$). The reason will be explained later. The $C_{2z}T$ transformation and its representation unitary matrix $D(C_{2z}T)$ are defined as follows:

[TABLE]

and the Hamiltonian satisfies $D(C_{2z}T)h^{*}(\mathbf{k})D^{-1}(C_{2z}T)=h(\mathbf{k})$.

V.3.1 $C_{2z}T$ Sewing matrix and off-diagonal Berry connection

The topology of these bands has been studied in previous research where it has been proved that they have a Wilson loop winding number Song et al. (2018); Ahn et al. (2018). Before we prove that the metric is bounded by that number, we briefly review the topology of these flat bands. In momentum space, $T$ sends $\mathbf{k}$ to $-\mathbf{k}$ and then $C_{2z}$ sends it back to $\mathbf{k}$. Since $T^{2}=1$ in this case, there is no Kramers degeneracy, and we can have the sewing matrix of $C_{2z}T$ symmetry. If $u_{n}(\mathbf{k})$ are the eigenstates of the Hamiltonian $h(\mathbf{k})$, the $C_{2z}T$ sewing matrix $B(\mathbf{k})$ will have the following form:

[TABLE]

where $\tilde{u}_{\mathbf{k}}=(u_{1}(\mathbf{k}),u_{2}(\mathbf{k}))$ are the eigenvectors of $h(\mathbf{k})$ which correspond to the flat bands. $C_{2z}T$ symmetry is antiunitary, hence we need a real embedding matrix $V_{\mathbf{G}}$ to make $B(\mathbf{l})$ periodic:

[TABLE]

Because $(C_{2z}T)^{2}=1$, there is no Kramers theorem, then if the two flat bands are not degenerate (away from the Dirac points), the sewing matrix must be diagonal and can be written as

[TABLE]

Because of the periodicity of $B(\mathbf{k})$, we must have $(\theta_{1\mathbf{k}+\mathbf{G}}-\theta_{1\mathbf{k}})=2\pi n,(\theta_{2\mathbf{k}+\mathbf{G}}-\theta_{2\mathbf{k}})=2\pi m$ in which $n,m\in\mathbb{Z}$. By definition, the non-Abelian Berry connection is given by $\mathbf{A}(\mathbf{k})=i\tilde{u}^{\dagger}_{\mathbf{k}}\partial_{\mathbf{k}}\tilde{u}_{\mathbf{k}}$, where $\tilde{u}_{\mathbf{k}}=(u_{1}(\mathbf{k}),u_{2}(\mathbf{k}))$ are the eigenvectors of $h(\mathbf{k})$ which correspond to the flat bands. Due to the property of the sewing matrix Eq. (S94), when $\mathbf{k}$ is not on any Dirac points, the Berry connection satisfies the following equation Song et al. (2018); Ahn et al. (2018):

[TABLE]

and this equation is a constraint on the matrix elements of the Berry connection. It requires the non-Abelian connection to have the following form:

[TABLE]

in which $\mathbf{a}(\mathbf{k})$ is a real vector. It is important to know whether $\mathbf{a}(\mathbf{k})$ is single valued or not. Because the embedding matrix $V_{\mathbf{G}}$ is momentum independent, the Berry connection itself is periodic. Also the winding numbers $n_{1,2}$ defined by phase factors $\theta_{1,2\mathbf{k}+\mathbf{G}}-\theta_{1,2\mathbf{k}}=2\pi n_{1,2}$ must be integers because of periodicity of $B(\mathbf{k})$. Therefore we have $e^{\frac{i}{2}(\theta_{1\mathbf{k}+\mathbf{G}}-\theta_{2\mathbf{k}+\mathbf{G}})}=\pm e^{\frac{i}{2}(\theta_{1\mathbf{k}}-\theta_{2\mathbf{k}})}$. This means the vector $\mathbf{a}(\mathbf{k})$ can be periodic ($\mathbf{a}(\mathbf{k}+\mathbf{G})=\mathbf{a}(\mathbf{k})$) if $(n_{1}-n_{2})$ is an even integer, or antiperiodic ($\mathbf{a}(\mathbf{k}+\mathbf{G})=-\mathbf{a}(\mathbf{k})$) if $(n_{1}-n_{2})$ is an odd integer. . In fact, by using the property of $C_{2z}T$ symmetry, we can prove that $\mathbf{a}(\mathbf{k})$ must be periodic for topological nontrivial bands.

When $\mathbf{k}$ is not on any Dirac points, we can write down the form of non-Abelian Berry curvature by Eq. (S97):

[TABLE]

We will show how the Wilson loop winding number is related with $f_{xy}$ in the following subsection.

V.3.2 Wilson loop winding

It can be shown that the winding number of the Wilson loop of the two active bands can be expressed as the integral of $f_{xy}$ on the whole Brillouin zone. We can label the momentum as $\mathbf{k}=\frac{k_{1}}{2\pi}\mathbf{b}_{1}+\frac{k_{2}}{2\pi}\mathbf{b}_{2}$ in which $\mathbf{b}_{1,2}$ are the primitive vectors of the Moiré reciprocal lattice. The large Wilson loop across the whole BZ along $\mathbf{b}_{2}$ direction with fixed $k_{1}$ is defined as

[TABLE]

where $\mathcal{P}$ stands for the path ordering and $c$ is the straight line path from $(k_{1},0)$ to $(k_{1},2\pi)$. The complex conjugation of this Wilson loop operator will become:

[TABLE]

Since we have proved that the sewing matrix is periodic, we find that $W(k_{1})$ and $W^{*}(k_{1})$ only differ by a unitary transformation. There will be two kinds of possiblities: 1) the two eigenvalues are complex conjugation to each other; 2) the two eigenvalues are real, and because of unitarity of $W(k_{1})$, the two eigenvalues can only be $\{1,1\}$, $\{1,-1\}$ or $\{-1,-1\}$. However, the second case can only give us a trivial band because the Wilson loop is not winding, and hence we only consider the first case in the following discussion. Because the two eigenvalues are complex conjugation to each other, the determinant of this large Wilson loop must be $1$. Therefore, we have

[TABLE]

When there is no Dirac point along this Wilson loop, from Eq. (S97) we know that ${\rm Tr}\,\mathbf{A}(\mathbf{k})=\frac{1}{2}\partial_{\mathbf{k}}(\theta_{1\mathbf{k}}+\theta_{2\mathbf{k}})$, thus $\oint_{c}d\mathbf{k}\cdot{\rm Tr}\,\mathbf{A}(\mathbf{k})=\pi(n_{1}+n_{2})$, where $n_{1,2}$ are the winding numbers of $\theta_{1,2\mathbf{k}}$ along this large Wilson loop. Because of Eq. (S101), $(n_{1}+n_{2})$ must be an even integer, and thus $(n_{1}-n_{2})$ is also an even integer. As we mentioned in last subsection, the off-diagonal element of Berry connection $i\mathbf{a}(\mathbf{k})e^{\frac{i}{2}(\theta_{1\mathbf{k}}-\theta_{2\mathbf{k}})}$ is periodic. Since $e^{\frac{i}{2}(\theta_{1\mathbf{k}+\mathbf{G}}-\theta_{2\mathbf{k}+\mathbf{G}})}=e^{i\pi(n_{1}-n_{2})}e^{\frac{i}{2}(\theta_{1\mathbf{k}}-\theta_{2\mathbf{k}})}=e^{\frac{i}{2}(\theta_{1\mathbf{k}}-\theta_{2\mathbf{k}})}$, we find that $\mathbf{a}(\mathbf{k}+\mathbf{G})=\mathbf{a}(\mathbf{k})$.

The above discussion is about a large Wilson loop $W(k_{1})$ along BZ without going across any Dirac points. However the topology of bands is usually diagnosed by the winding of Wilson loop over the whole BZ, which means $k_{1}$ varies from [math] to $2\pi$ and it must encounter Dirac points during this process. It can be shown that $W(k_{1})$ is continuous when it goes across a Dirac point. Assume that $W_{D}$ is a non-Abelian Wilson loop surrounding a Dirac cone at position $(k_{01},k_{02})$. If we want to shrink this loop to a point with the Dirac point surrounded by it, then the shrinking point must be that Dirac point. The Wilson loop is made out of the projector into the two bands that have the Dirac point, hence it must be the identity. As shown in FIG. S2 (a), the Wilson loop $W(k_{01}-\delta)$ and $W(k_{01}+\delta)$ can be represented as

[TABLE]

and these Wilson line operators $W_{1}^{\pm}$, $W_{2}^{\pm}$, $W_{Dl}$ and $W_{Dr}$ are also shown in FIG. S2 (a). As we just mentioned, the Wilson loop around the Dirac point is the identity when the radius $D$ of that Wilson loop approaches zero, which means $\mathds{1}=W_{Dl}W_{Dr}$ and hence $W_{Dr}=W_{Dl}^{-1}$. Also because $W_{1,2}^{-}$ can be continuously deformed to $W_{1,2}^{+}$ without crossing with any Dirac point, we will obtain that $W_{1,2}^{-}=W_{1,2}^{+}$ when $\delta\rightarrow 0$. One hence easily finds that $W(k_{01}-\delta)$ and $W(k_{01}+\delta)$ are identical at the limit when both $D$ and $\delta$ go to zero. That is our statement that $W(k_{1})$ is continuous when it goes across a Dirac point.

Now we study how the winding number of Wilson loop is related with the off-diagonal Berry connection/curvature which defines the Euler class. To do this we need to know the eigenvalue of the Wilson loop operator. However, the form of its eigenvalue depends on the winding numbers of $\theta_{1\mathbf{k}}$ and $\theta_{2\mathbf{k}}$, i.e., $n_{1},\leavevmode\nobreak\ n_{2}$, along the path of that Wilson loop. As we proved in previous paragraph, $(n_{1}-n_{2})$ is an even integer and there will be two possibilities: 1) both of them are even integers, or 2) both of them are odd integers. If we calculate the Abelian Wilson loop of one band (either of the bands forming the Dirac node) around a Dirac point, the Berry phase must be $\pm\pi$. This means the loop integral of both $\mathbf{A}_{11}(\mathbf{k})$ or $\mathbf{A}_{22}(\mathbf{k})$ around the Dirac point is either $\pm\pi$. From Eq. (S97) we know $\mathbf{A}_{11}(\mathbf{k})=\frac{1}{2}\partial_{\mathbf{k}}\theta_{1\mathbf{k}}$. Thus the winding number of $\theta_{1\mathbf{k}}$ around the Dirac point will be $\pm 2\pi$. Because there is a nonzero winding of $\theta_{1\mathbf{k}}$ around the Dirac point, there will be a branch cut, and the value of $\theta_{1}$ on two sides of the branch cut differs by $\pm 2\pi$. Now we consider two paths $c_{1}$ and $c_{2}$ shown in FIG. S2 (b) and the dashed line stands for the branch cut. Both $c_{1}$ and $c_{2}$ are a large loop which go across the whole BZ. There is a branch cut between the points $(k_{01}-\delta,2\pi)$ and $(k_{01}+\delta,2\pi)$ , but there is no branch between $(k_{01}-\delta,0)$ and $(k_{01}+\delta,0)$. Hence when $\delta\rightarrow 0$, we have $\theta_{1}(k_{01}-\delta,0)=\theta_{1}(k_{01}+\delta,0)$ and $\theta_{1}(k_{01}-\delta,2\pi)=\theta_{1}(k_{01}+\delta,2\pi)\pm 2\pi$. Therefore, the winding number $n_{1}$ changes by one when the path goes across a Dirac point. Similar argument can also be applied to $\theta_{2}$, because $\mathbf{A}_{22}(\mathbf{k})=\frac{1}{2}\partial_{\mathbf{k}}\theta_{2\mathbf{k}}$. In conclusion, $n_{1}$ and $n_{2}$ changes from even (odd) to odd (even) when the path of the Wilson loop goes across a Dirac point.

The eigenvalue of the Wilson loops $W(k_{1})$ can be calculated as shown. To obtain an expression without $\theta_{1,2\mathbf{k}}$, we do the following transformation $u^{\prime}_{\mathbf{k}}=\tilde{u}_{\mathbf{k}}g_{\mathbf{k}}$, in which the transformation $g_{\mathbf{k}}$ is given by

[TABLE]

Then the Wilson loop $W(k_{1})$ can be expressed by $u^{\prime}_{\mathbf{k}}$ as follows:

[TABLE]

When $\theta_{1\mathbf{k}}$ and $\theta_{2\mathbf{k}}$ wind even times, $g(k_{1},0)=g(k_{1},2\pi)$, the eigenvalues of $W(k_{1})$ is the same as $W^{\prime}(k_{1})$. If the winding number is odd, then $g(k_{1},0)=-g(k_{1},2\pi)$, the eigenvalue of $W(k_{1})$ is the same as $-W^{\prime}(k_{1})$. It is hence important to know the expression of $W^{\prime}(k_{1})$. From the expression of $W^{\prime}(k_{1})$, we find that it can be written as

[TABLE]

where $\mathcal{P}$ stands for path ordering. We use Eqs. (S97) and (S102) to derive $\mathbf{A}^{\prime}$, and the result is $\mathbf{A}^{\prime}(\mathbf{k})=-\mathbf{a}(\mathbf{k})\sigma_{2}$. Because $\sigma_{2}$ always commute with itself, $W^{\prime}(k_{1})$ can be expressed without the use of the path ordering:

[TABLE]

By this equation we can get the eigenvalue of Wilson loop $W(k_{1})$ as shown:

[TABLE]

At first glance the Wilson loop spectrum is not a continuous function of $k_{1}$. However, as we will see in the following paragraph, the winding of $\mathbf{a}(\mathbf{k})$ around the Dirac point must be an odd integer, which means $\int_{0}^{2\pi}dk_{2}\,a_{2}(k_{1},k_{2})$ is not continuous, and as a result we obtain a continuous eigenvalue $\xi(k_{1})$. This statement can be proved as follows. By using the transformation in Eq. (V.3.2), the Wilson loop of the two bands around the Dirac point that they form can be written as

[TABLE]

In this equation the Wilson loop starts at point $\mathbf{k}_{i}$, and the it goes around the Dirac point and then ends at $\mathbf{k}_{f}$. Actually $\mathbf{k}_{i}$ and $\mathbf{k}_{f}$ are at the same position in the momentum space but the matrix $g_{\mathbf{k}}$ are different because of the winding of $\theta_{1,2}$. As we proved in the text, the winding numbers of $\theta_{1\mathbf{k}}$ and $\theta_{2\mathbf{k}}$ around the Dirac point are $\pm 2\pi$, which means $g_{\mathbf{k}_{i}}=-g_{\mathbf{k}_{f}}$. Also the non-Abelian Wilson loop is identity when the radius of the loop goes to zero $W_{D}|_{D\rightarrow 0}=\mathds{1}$ , so the following equation holds:

[TABLE]

where $D$ is the infinitesimal loop around a Dirac point. This means the winding of $\mathbf{a}(\mathbf{k})$ is an odd number times $\pi$:

[TABLE]

Now if a Dirac point is assumed to be located at $(k_{01},k_{02})$, we find that the integral $\int_{0}^{2\pi}dk_{2}\,a_{2}(k_{1},k_{2})$ is not continuous and it satisfies:

[TABLE]

where the paths $c_{1,2,3}^{\pm}$ are defined in FIG. S3 (a). When $\delta\rightarrow 0$, we have $\int_{c_{1,3}^{+}}d\mathbf{k}\cdot\mathbf{a}(\mathbf{k})=\int_{c_{1,3}^{-}}d\mathbf{k}\cdot\mathbf{a}(\mathbf{k})$. Therefore, Eq. (S110) becomes the integral of $\mathbf{a}(\mathbf{k})$ along the path $c_{2}^{+}-c_{2}^{-}$, which is actually a loop around the Dirac point. This can tell us that $\int_{0}^{2\pi}dk_{2}\,a_{2}(k_{1},k_{2})$ jumps by $(2n+1)\pi$ at $k_{10}$ if a Dirac point is located at $(k_{01},k_{02})$. Hence, the eigenvalue of $W(k_{1})$ is continuous, because

[TABLE]

That result is consistent with our argument shown in FIG. S2 (a). Now suppose there is no Dirac point between $k_{1}$ and $k_{1}^{\prime}$, then by using Stokes’ theorem, the increment of $\xi$ is given by the integral of $f_{xy}$

[TABLE]

where $\Omega_{k_{1},k_{1}^{\prime}}$ is a rectangle with points $(k_{1},0),\leavevmode\nobreak\ (k_{1},2\pi),\leavevmode\nobreak\ (k_{1}^{\prime},2\pi)$ and $(k_{1}^{\prime},0)$. If there is a Dirac point between $k_{1}$ and $k_{1}^{\prime}$ at $k_{01}$, then we can apply Stokes’ theorem on both sides of Dirac point as shown in FIG. S3 (b), and we obtain the following equations

[TABLE]

Because of the continuity of $\xi(k_{1})$, we know that $\xi(k_{01}+\delta)=\xi(k_{01}-\delta)$ when $\delta\rightarrow 0$. Then by adding the two equations together, we have

[TABLE]

in which $\Omega^{\prime}_{k_{1},k_{1}^{\prime}}$ stands for the same the rectangular region as $\Omega_{k_{1},k_{1}^{\prime}}$, with Dirac points removed. This is consistent with the definition of $f_{xy}$, because it is well defined only when $\mathbf{k}$ is not on any Dirac points. Now we can write down the Wilson loop winding number in terms of the integral of $f_{xy}$:

[TABLE]

where BZ’ is the whole Brillouin zone with Dirac points removed from the integration area.

V.3.3 Fubini-Study metric of $C_{2z}T$ topology

The $C_{2z}T$ symmetry yields a constraint on the quantum geometric tensor. The quantum geometric tensor satisfies the following equation:

[TABLE]

where $B(\mathbf{k})$ is the sewing matrix of the $C_{2z}T$ symmetry. Similar to the Berry connection case, this equation gives a constraint on the matrix elements of the quantum geometric tensor. It can be shown that $\mathfrak{G}_{ij}$ can only take the following form:

[TABLE]

in which $\gamma_{ij}=\mathfrak{g}^{12}_{ij}e^{\frac{i}{2}(\theta_{2}-\theta_{1})}$ is a real symmetric tensor. By definition the Fubini-Study metric of the two flat bands with spin $\uparrow$ and valley $K$ in TBLG is given by $g_{ij}=\mathfrak{g}^{11}_{ij}+\mathfrak{g}^{22}_{ij}$. Similar to the Chern insulator case, if we choose $c_{x}=(1,ie^{\frac{i}{2}(\theta_{2}-\theta_{1})})^{\rm T}$, $c_{y}=(i,-e^{\frac{i}{2}(\theta_{2}-\theta_{1})})^{\rm T}$, then Eq. (S66) can tell us

[TABLE]

similarly, if we choose $c_{x}=(1,-ie^{\frac{i}{2}(\theta_{2}-\theta_{1})})^{\rm T}$, $c_{y}=(i,e^{\frac{i}{2}(\theta_{2}-\theta_{1})})^{\rm T}$, we will have ${\rm tr}\,g\geq-2f_{xy}$. In conclusion, we have ${\rm tr}\,g\geq 2|f_{xy}|$. As we have shown in the last subsection, the integral of $f_{xy}$ gives us the $C_{2z}T$ winding number $e_{2}$. We hence find that the integral of metric is bounded by the winding number:

[TABLE]

VI Superconductivity in TBLG

We now discuss the superconductivity in TBLG. The TBLG system has two spin and two valley quantum number. Thus there are 8 flat bands in total. We denote the fermion annihilation operator by $\psi_{\eta,s,\mathbf{Q},a,\mathbf{k}}$, where $\eta$ stands for graphene valley, $s$ stands for electron spin, $\mathbf{Q}$ stands for the hexagonal lattice site in momentum space formed by $\mathbf{q}_{1,2,3}$ vectors and $a$ stands for graphene sublattice. The BM continuum model in Eq. (S92) is the Hamiltonian of spin $\uparrow$ and valley $K$. To obtain the Hamiltonian of valley $K^{\prime}$, we can use the time reversal transformation $\mathcal{T}$. The fermion operators $\psi$ transform as shown under $\mathcal{T}$:

[TABLE]

in which $s_{i}$, $\tau_{i}$ and $\sigma_{i}$ are Pauli matrices acting on spin, graphene valley and graphene sublattice indices, respectively. The matrix $V$ acting on $\mathbf{Q}$ index is defined by $V_{\mathbf{Q},\mathbf{Q}^{\prime}}=\delta_{\mathbf{Q},-\mathbf{Q}^{\prime}}$. Therefore the representation matrix of spinful time reversal transformation in TBLG is $U_{T}=is_{y}\tau_{x}V$.

By taking all the spin and valley into consideration, the kinetic energy $\mathcal{H}(\mathbf{k})$ and time reversal transformation representation $U_{T}$ can be written as the following diagonal block matrices:

[TABLE]

in which the entries of this matrix stand for $(K\uparrow),(K^{\prime}\uparrow),(K\downarrow)$ and $(K^{\prime}\downarrow)$. Because there is no spin-orbit-coupling in graphene, we have $h_{K\uparrow}(\mathbf{k})=h_{K\downarrow}(\mathbf{k})=h(\mathbf{k})$, where $h(\mathbf{k})$ is the BM continuum model in Eq. (S92), and $h_{K^{\prime}\uparrow}(\mathbf{k})=h_{K^{\prime}\downarrow}(\mathbf{k})=h^{\prime}(\mathbf{k})$. Since the Hamiltonian $\mathcal{H}(\mathbf{k})$ is time reversal invariant, we have $U^{\dagger}_{T}\mathcal{H}^{*}(-\mathbf{k})U_{T}=\mathcal{H}(\mathbf{k})$. This condition tells us that $h^{\prime}(\mathbf{k})=Vh^{*}(-\mathbf{k})V$. Suppose the eigenvectors of $h(\mathbf{k})$ are denoted by $u_{n}(\mathbf{k})$, the flat bands wave function will be

[TABLE]

where $\tilde{u}_{\mathbf{k}}=(u_{1}(\mathbf{k}),u_{2}(\mathbf{k}))$ and $u_{1,2}(\mathbf{k})$ are the eigenvectors of $h(\mathbf{k})$ which corresponds to the two flat bands.

We assume that the pairing order parameter has the form of Eq. (S41). This pairing order parameter means that the pairing happens between electrons which has opposite spin, valley and momentum, because the time reversal transformation $U_{T}=is_{y}\tau_{x}V$ can flip the spin and valley at the same time. Therefore, by using the result in Sec. III, the superfluid weight of TBLG is given by Eq. (S60). However, since $\tilde{U}_{\mathbf{k}}$ are the eigenvectors which corresponds to the all 8 flat bands with all spins and valleys, we cannot directly use the lower bound of Fubini-Study metric of one spin one valley $g_{ij}(\mathbf{k})$. In fact, integrand in Eq. (S60) is related with $g_{ij}(\mathbf{k})$ as shown:

[TABLE]

However, as we have mentioned, Fubini-Study metric defines a positive distance in the BZ, which means it must be a real tensor, hence $g^{*}_{ij}(-\mathbf{k})=g_{ij}(\mathbf{k})$. Therefore we have

[TABLE]

Then by using Eq. (S119) the lower bound of this integral will be

[TABLE]

That naturally leads to the lower bound of the trace of the superfluid weight: ${\rm tr}\,D_{s}\geq\frac{8(\Delta_{1}+\Delta_{2})}{\pi}\sqrt{\nu(1-\nu)}$. Here we used the fact that the winding number of the flat bands in TBLG is one Song et al. (2018). Due to $C_{3z}$ symmetry, $D_{s}$ must proportional to the identity, therefore we obtain $D_{s}\geq\frac{4(\Delta_{1}+\Delta_{2})}{\pi}\sqrt{\nu(1-\nu)}$.

VII Finite temperature study and Berezinskii-Kosterlitz-Thouless transition

As mentioned in the main text, we will assume $\Delta_{2}=0$. Before we estimate the actual value of topological lower bound, we need to know the gap function $\Delta_{1}$ at zero temperature. As we discussed in the main text, the transition temperature measured in experiment is the Berezinskii-Kosterlitz-Thouless temperature $T_{c}$, and it is different from the mean field temperature $T_{c}^{*}$ where Cooper pairing disappears. The condition of BKT transition is universal and it is given by

[TABLE]

where $D_{s}(T)$ is the superfluid weight at temperature $T$. When this transition happens, the superfluid weight jumps from a finite value to zero, which is totally different from the situation in 3D. Hence the BKT transition temperature is always lower $T_{c}<T_{c}^{*}$. Now we use the mean field Hamiltonian to determine the ratio of BKT temperature and mean field temperature. The free energy at finite temperature is given by the following equation:

[TABLE]

in which $E_{\mathbf{k}n}(\mathbf{q})$ are the positive eigenvalues of $\tilde{\mathscr{H}}_{\mathbf{k}}(\mathbf{q})$. Now we take the derivative of this free energy, and the superfluid weight will be

[TABLE]

In Sec. III we have shown that $\partial_{\mathbf{q}}\lambda_{\mathbf{k}n}(\mathbf{q})|_{\mathbf{q}=0}=\partial_{\mathbf{q}}\varphi_{\mathbf{k}n}(\mathbf{q})|_{\mathbf{q}=0}=0$, so all the first order derivative terms in Eq. (S128) will disappear. Also in flat band limit, $E_{n\mathbf{k}}(0)$ is independent with $\mathbf{k}$ and $n$ under the pairing ansatz Eq. (S41), thus the superfluid weight at finite temperature will become

[TABLE]

and we find this result only differs by a hyperbolic tangent factor from the expression at zero temperature in Sec. III. By using Eq. (S58), we can write down the expression of $D_{s}$ at finite temperature. Since it is isotropic due to $C_{3z}$ symmetry, the superfluid weight can be written as

[TABLE]

here $\Delta_{1}(T)$ is the gap function at temperature $T$, and $g(\mathbf{k})$ is the Fubini-Study metric derived from all the flat bands. For convenience we denote $\sqrt{(\varepsilon_{0}-\mu)^{2}+\Delta_{1}(T)}$ by $E_{B}(T)$. Notice that this expression depends on the chemical potential $\mu$ instead of the filling ratio $\nu$. The particle number is no longer simply given by the Bogoliubov quasiparticle coefficient $|v_{\mathbf{k}}|^{2}$ here because of thermal excitations. Because both the free fermion band and Bogoliubov quasiparticle spectrum $E_{B}(T)$ are momentum independent, the filling ratio of flat bands are given by

[TABLE]

here $u^{2},v^{2}=\frac{1}{2}\left(1\pm\frac{\varepsilon_{0}-\mu}{E_{B}(T)}\right)$ are the coefficients of Bogoliubov transformation, $c,c^{\dagger}$ are the electron operators and $\gamma,\gamma^{\dagger}$ are the Bogoliubov quasiparticle operators. When $T\neq 0$, the expectation value of Bogoliubov quasiparticle number is given by Fermi-Dirac distribution, therefore we have $\langle\gamma^{\dagger}_{\mathbf{k},K\uparrow}\gamma_{\mathbf{k},K\uparrow}\rangle=n_{F}(E_{B}(T))$ and $\langle\gamma_{-\mathbf{k},K^{\prime}\downarrow}\gamma^{\dagger}_{-\mathbf{k},K^{\prime}\downarrow}\rangle=1+n_{F}(E_{B}(T))$. Hence the filling ratio will be

[TABLE]

If the temperature $T$, the order parameter $\Delta_{1}(T)$ and the filling ratio $\nu$ are given, the chemical potential and further the Bogoliubov energy $E_{B}(T)$ can be solved from Eq. (S132). Finally we are able to use Eq. (S130) to obtain the superfluid weight with given temperature $T$ and filling ratio $\nu$.

To determine the relationship between BKT temperature and mean field temperature, we have to know the value of $D_{s}(T)$. Here we use flat bands of TBLG as an example. For simplicity, we use the lower bound of the metric to estimate its integral:$\int\frac{d^{2}k}{(2\pi)^{2}}{\rm tr\,}g\approx\frac{4}{\pi}$. Hence the superfluid weight becomes

[TABLE]

in which the chemical potential is related with filling ratio $\nu$ by Eq. (S132). Since it is an order of magnitude estimation here, we can simply use the BCS mean field result $\Delta_{1}(T)\approx\Delta_{1}(0)\left(1-\frac{2T}{\Delta_{1}(0)}\right)^{1/2}$ to determine the temperature dependent superconducting gap. Then we can numerically evaluate the value of Eq. (S133), and the result is shown in FIG. 1 when $\nu=1/4$. The intersection of the solid line and the dotted line is the KT transition point, which satisfies Eq. (S126). We find that when $\nu=1/4$, the transition temperature is $T_{c}\approx 0.34T_{c}^{*}\approx 0.17\Delta_{1}(0)$. Here we were using the lower bound of superfluid weight, which means the ratio $T_{c}/T_{c}^{*}$ is under estimated.

VIII The validity of BCS wave function in the flat band limit

By using Hund’s rule we argue that the BCS mean field ground state is a good approximation of the real ground state. This argument was also applied in the spin Chern number case Peotta and Törmä (2015) and originally this idea comes from quantum Hall ferromagnet Spielman et al. (2000). The BCS wave function can be written as

[TABLE]

in which $n$ is the flat band index. This wave function is equivalent to the following form:

[TABLE]

If we perform the following particle-hole transformation: $c^{\dagger}_{\downarrow,\mathbf{k}}\rightarrow c_{\downarrow,-\mathbf{k}}$. That transformation changes the attractive interaction into repulsive. Originally, the interaction between electrons with opposite spins is attractive. After the transformation, it becomes an attractive interaction between spin up electrons and spin down holes, which is equivalent to a repulsive interaction between electrons with opposite spin. In order to show this clearly, we use the following toy model Hamiltonian

[TABLE]

where $g(\mathbf{q})$ is the Fourier transformation of an attractive potential, and $v$ stand for valley, and $\bar{v}$ is the opposite valley of $v$. After this particle-hole transformation, the toy model Hamiltonian becomes the following form:

[TABLE]

We notice that the sign of $g(\mathbf{q})$ is flipped, which means a repulsive interaction. However, by switching the creation and annihilation operators, the vacuum state is also changed: $\prod_{\mathbf{k},n}c^{\dagger}_{K^{\prime}\downarrow,-\mathbf{k}n}c^{\dagger}_{K\downarrow,-\mathbf{k}n}|0\rangle\rightarrow|0\rangle$. Therefore, under this particle-hole transformation, the BCS wave function becomes

[TABLE]

If we use $(K\uparrow)$ and $(K^{\prime}\downarrow)$ or $(K\downarrow)$ and $(K^{\prime}\uparrow)$ to form a pseudo spin $\frac{1}{2}$ space, then $|\rm BCS\rangle$ transforms into a ferromagnetic state of the pseudo spin. Since the ferromagnetic state is symmetric in the pseudo spin space, its orbital wave function is totally antisymmetric. Therefore, this ferromagnetic state is energetically favored by repulsive interaction, and we conclude that BCS wave function can be a good approximation to the ground state.