Superfluid weight and Berezinskii-Kosterlitz-Thouless transition temperature of twisted bilayer graphene

TL;DR

This paper investigates the superconducting properties of twisted bilayer graphene, focusing on superfluid weight and BKT transition temperature, highlighting differences between local and non-local interactions and emphasizing the role of multiband processes and quantum geometry.

Contribution

It provides a detailed analysis of superfluid weight and BKT transition temperature in twisted bilayer graphene, incorporating both microscopic and low-energy models, and reveals qualitative differences based on interaction types.

Findings

Superfluid weight and BKT temperature depend on interaction schemes.

Quantum geometry influences superfluid properties in the flat band limit.

Distinct experimental signatures predicted for local vs. non-local interactions.

Abstract

We study superconductivity of twisted bilayer graphene with local and non-local attractive interactions. We obtain the superfluid weight and Berezinskii-Kosterlitz-Thouless (BKT) transition temperature for microscopic tight-binding and low-energy continuum models. We predict qualitative differences between local and non-local interaction schemes which could be distinguished experimentally. In the flat band limit where the pair potential exceeds the band width we show that the superfluid weight and BKT temperature are determined by multiband processes and quantum geometry of the band.

Equations

$$D^{s}_{\mu\nu}=\lim_{\mathbf{q}\rightarrow 0}\big{[}\lim_{\omega\rightarrow 0}K_{\mu\nu}(\mathbf{q},\omega)\big{]}.$$

Full text

Recently, it was found that $D^{s}$ has, in addition to the conventional contribution proportional to group velocity, a geometric contribution arising from multiband processes Peotta and Törmä (2015); Julku et al. (2016); Tovmasyan et al. (2016); Liang et al. (2017); Törmä et al. (2018). In a flat band limit the geometric contribution dominates and is bounded from below by the band Berry curvature Liang et al. (2017) and Chern number Peotta and Törmä (2015). Here we show that the geometric contribution dominates $D^{s}$ and $T_{\rm BKT}$ in the flat band regime of TBG. Importantly, we show that including only the few flat bands is not sufficient but one needs also a number of dispersive bands to correctly predict the geometric contribution. Therefore, approximate models of TBG such as those with only flat bands, as used for deriving upper Hazra et al. (2018) and lower Xie et al. (2019) bounds of the superfluid weight and in many other works Guo et al. (2018); Lian et al. (2018); Chen et al. (2019); Xu and Balents (2018); Liu et al. (2018); Tang et al. (2019); Roy and Juričić (2019); Sherkunov and Betouras (2018); Lin and Nandkishore (2018); Classen et al. (2019); Ray et al. (2019); Dodaro et al. (2018); Kang and Vafek (2019), may not be suited for quantitative predictions of TBG superconductivity. Moreover, we predict that, in the flat-band regime, local ($s$-wave) and non-local interactions yield distinct behavior, namely an anisotropic superfluid weight in the latter case. We propose a four-terminal radio frequency spectroscopy experiment that can detect the possible anisotropy and thus distinguish between the two pairing mechanisms.

An outstanding problem in describing the TBG physics theoretically Liu et al. (2018); Huang et al. (2019); Ray et al. (2019); Kozii et al. (2019); Liu et al. (2018); Guo et al. (2018); Kennes et al. (2018); Choi and Choi (2018); Peltonen et al. (2018); Wu et al. (2018); Lian et al. (2018); Wu (2019); Wu and Sarma (2019); Wu et al. (2019a); Sherkunov and Betouras (2018); Su and Lin (2018); Xu and Balents (2018); Roy and Juričić (2019); Po et al. (2018); Isobe et al. (2018); Laksono et al. (2018); Fidrysiak et al. (2018); Tang et al. (2019); Venderbos and Fernandes (2018); Dodaro et al. (2018); Lin and Nandkishore (2018); González and Stauber (2019); Classen et al. (2019); Chen et al. (2019) is the fact that the unit cell of the moiré superlattice with twist angles close to $\theta^{*}$ contains a large amount of carbon atoms [Fig. 1(a)], and therefore TBG theory should take a stand on how to describe the interlayer couplings within this unit cell. Here we use and compare with each other two of the previously used approximation procedures: (1) the renormalized moiré (RM) approach Gonzalez-Arraga et al. (2017); Su and Lin (2018), where we scale some coupling energies by a suitable scaling factor to find the flat bands at a higher $\theta$, resulting into a smaller unit cell, and (2) the Dirac point approximation (DP) Peltonen et al. (2018); Lopes dos Santos et al. (2007, 2012), where we make a low energy approximation near the graphene Dirac points by linearizing the intralayer Hamiltonians and using a cut-off in the superlattice Fourier space. Both of these approaches go beyond those often used in TBG literature, either based on a single-parameter coupling model, or a vastly reduced four-band model Yuan and Fu (2018); Kang and Vafek (2018); Koshino et al. (2018); Hazra et al. (2018); Xie et al. (2019).

Theoretical models In the renormalized moiré lattice method (RM), we deploy the Fermi-Hubbard Hamiltonian as sup $H=H_{\textrm{kin}}-\mu N+H_{\textrm{int}}$, where $H_{\textrm{kin}}=\sum_{i\alpha j\beta\sigma}t_{i\alpha j\beta}c^{\dagger}_{i\alpha\sigma}c_{j\beta\sigma}$ is the kinetic term, $N$ is the total particle number operator and $H_{\textrm{int}}$ is the effective attractive interaction described below. Here $c_{i\alpha\sigma}$ annihilates a fermion in the $\alpha$th lattice site of the $i$th moiré superlattice unit cell with spin $\sigma\in\{\uparrow,\downarrow\}$, $\mu$ is the chemical potential and the hopping $t_{i\alpha j\beta}$ includes both the intra- and interlayer terms.

Since the type of the interaction is not currently known, we consider two different singlet pairing potentials, namely the local pairing $H_{\textrm{int}}=J\sum_{i\alpha}c^{\dagger}_{i\alpha\uparrow}c^{\dagger}_{i\alpha\downarrow}c_{i\alpha\downarrow}c_{i\alpha\uparrow}\equiv H_{\textrm{loc}}$ and the nearest-neighbour (NN) pairing $H_{\textrm{int}}=\frac{J}{2}\sum_{\langle i\alpha j\beta\rangle}h^{\dagger}_{i\alpha j\beta}h_{i\alpha j\beta}\equiv H_{\textrm{RVB}}$, where $h_{i\alpha j\beta}=(c_{i\alpha\uparrow}c_{j\beta\downarrow}-c_{i\alpha\downarrow}c_{j\beta\uparrow})$ and $J<0$ is the interaction strength. The local interaction has been used to study $s$ wave superconductivity, mediated by electron-phonon interaction, both in graphene Zhao and Paramekanti (2006); Uchoa and Castro Neto (2007); Kopnin and Sonin (2008); Hosseini (2015) and in TBG Peltonen et al. (2018); Wu et al. (2018). The non-local, called resonance valence bond (RVB) interaction Anderson (1987); Lee et al. (2006), has also been used both in case of monolayer graphene Baskaran (2002); Black-Schaffer and Doniach (2007); Pathak et al. (2010); Black-Schaffer (2012) and TBG Su and Lin (2018). We keep only the pairing channels by applying mean-field theory to approximate $H_{\textrm{loc}}\approx\Delta_{i\alpha}c_{i\alpha\uparrow}^{\dagger}c_{i\alpha\downarrow}^{\dagger}+\textrm{H.c.}$ and $H_{\textrm{RVB}}\approx\Delta_{i\alpha j\beta}h_{i\alpha j\beta}^{\dagger}+\textrm{H.c.}$, where $\Delta_{i\alpha}=J\langle c_{i\alpha\downarrow}c_{i\alpha\uparrow}\rangle$ and $\Delta_{i\alpha j\beta}=\frac{J}{2}\langle h_{i\alpha\beta}\rangle$ are the superfluid order parameters, respectively.

To reduce the number of lattice sites $M$ within a moiré unit cell (around 12 000 for twist angle $\theta\sim$1\text{,}\mathrm{\SIUnitSymbolDegree}), we apply a rescaling trick Su and Lin ([2018](#bib.bib54)); Gonzalez-Arraga *et al.* ([2017](#bib.bib61)) under which the Fermi velocity of a monolayer graphene and the moiré periodicity remain invariant but $\theta$ becomes larger and thus reduces $M$. In our computations we use the rescaling such that $M=676$ and the rescaled angle is $\theta^{\prime}=$4.41\text{\,}\mathrm{\SIUnitSymbolDegree} sup which reproduces the four narrow bands of the bandwidth of $10\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ found experimentally with $\theta\sim$1\text{,}\mathrm{\SIUnitSymbolDegree}$$ [see Fig. 1(b)].

In the Dirac point continuum method (DP) we employ the low-energy Lopes dos Santos et al. (2012); Peltonen et al. (2018); sup Dirac point approximation for the two graphene layers as $H_{\text{kin}}^{1}=\sum_{\sigma\rho\mathbf{k}\mathbf{G}}c_{\sigma\rho,1}^{\dagger}(\mathbf{k}+\mathbf{G})\hbar v_{F}\bm{\sigma}^{\rho}\dotproduct(\mathbf{k}+\mathbf{G})c_{\sigma\rho,1}(\mathbf{k}+\mathbf{G})$ and $H_{\text{kin}}^{2}=\sum_{\sigma\rho\mathbf{k}\mathbf{G}}c_{\sigma\rho,2}^{\dagger}(\mathbf{k}+\mathbf{G})\hbar v_{F}\bm{\sigma}_{\theta}^{\rho}\dotproduct(\mathbf{k}+\mathbf{G})c_{\sigma\rho,2}(\mathbf{k}+\mathbf{G})$ and couple the layers by $H_{\text{kin}}^{\perp}=\sum_{\sigma\rho\mathbf{k}\mathbf{G}\mathbf{G}^{\prime}}c_{\sigma\rho,1}^{\dagger}(\mathbf{k}+\mathbf{G}+\rho\frac{\Delta\mathbf{K}}{2})t_{\perp}^{\rho}(\mathbf{G}-\mathbf{G}^{\prime})c_{\sigma\rho,2}(\mathbf{k}+\mathbf{G}^{\prime}-\rho\frac{\Delta\mathbf{K}}{2})+\text{H.c.}$ Here $c_{\sigma\rho,l}(\mathbf{k})=(c_{\sigma\rho,lA}(\mathbf{k}),c_{\sigma\rho,lB}(\mathbf{k}))^{\mathsf{T}}$ in the sublattice space, where $c_{\sigma\rho,ls}(\mathbf{k})$ is the annihilation operator for spin $\sigma$, valley $\rho\in\{+,-\}$, layer $l$, sublattice $s$, and wavevector $\mathbf{k}$, $\bm{\sigma}^{\rho}=(\rho\sigma_{x},\sigma_{y})$ is a vector of Pauli matrices in the sublattice space, $\bm{\sigma}_{\theta}^{\rho}=R(\theta)\bm{\sigma}^{\rho}$ is the $\theta$-rotated version of it, $t_{\perp}^{\rho}(\mathbf{G})$ is the Fourier component not of a Slater-Koster Slater and Koster (1954) parametrized interlayer potential (times an exponential factor), $\Delta\mathbf{K}=R(\theta)\mathbf{K}-\mathbf{K}$ is the difference vector from the graphene $\mathbf{K}$ point to its rotated counterpart, and $v_{F}$ is the graphene Fermi velocity. The $\mathbf{k}$ sum is over the the moiré Brillouin zone and the $\mathbf{G},\mathbf{G}^{\prime}$ sums are over the (truncated) reciprocal superlattice.

We then write the total Hamiltonian as $H=H_{\text{kin}}^{1}+H_{\text{kin}}^{2}+H_{\text{kin}}^{\perp}-\mu N+H_{\text{int}}$, where $N$ is the total particle number operator. To describe the superconducting state with a local pairing interaction $\lambda$ we use $H_{\text{int}}=\lambda\sum_{ls}\int\differential{\mathbf{r}}\psi_{\uparrow\rho,ls}^{\dagger}(\mathbf{r})\psi_{\downarrow\bar{\rho},ls}^{\dagger}(\mathbf{r})\psi_{\downarrow\bar{\rho},ls}(\mathbf{r})\psi_{\uparrow\rho,ls}(\mathbf{r})$, which is treated in the mean-field level sup . Here $\bar{\rho}$ is the opposite valley of $\rho$ and $\psi_{\sigma\rho,ls}(\mathbf{r})$ is the continuum electron field operator.

Order parameters, superfluid weight, and pairing symmetry In experiments Cao et al. (2018a); Yankowitz et al. (2019); Lu et al. (2019) superconducting (SC) and correlated insulating states have been observed with the magic angle twist such that insulating states emerge for the flat band fillings $\nu\in\{0,\pm 1,\pm 2,\pm 3\}$ and SC states surround the insulating states close to $\nu\in\{0,\pm 1,\pm 2\}$ with the SC phase near $\nu=-2$ being observed at temperature as high as $\sim$3\text{,}\mathrm{K}$$ Yankowitz et al. (2019); Lu et al. (2019). Here $\nu$ is the electron density per moiré unit cell so that the charge neutrality point (CNP) corresponds to $\nu=0$ and narrow bands are empty (full) when $\nu=-4$ ($\nu=4$).

To determine the superfluid weight $D^{s}$, we first solve order parameters from the BCS gap equations sup . In Figs. 2(a)–(b) we show the spatial profiles for the local and RVB interactions computed with the RM method at $\nu\approx-2$. Here $J$ is chosen such that the maximum value of the order parameter is $\max|\Delta|\approx$3.4\text{,}\mathrm{m}\mathrm{e}\mathrm{V}$$. From Figs. 2(c)–(d) we see $\max|\Delta|$ depending almost linearly on the interaction constant, which is typical for generic flat band systems Kopnin et al. (2011); Peotta and Törmä (2015); Julku et al. (2016); Ojajärvi et al. (2018); Peltonen et al. (2018). From the obtained order parameter values one can compute $D^{s}$. For easier comparison between the RM and DP models, below we use $\max|\Delta|$ as a “parameter”.

To obtain $D^{s}$ we use linear response theory. In the mean-field level Scalapino et al. (1992, 1993) the zero-frequency, long wavelength limit of the current-current response function $K_{\mu\nu}(\mathbf{q},\omega)$ is $D^{s}$, i.e. ($\mu,\nu\in\{x,y\}$)

[TABLE]

In Ref. Liang et al., 2017 this was computed for a generic multi-orbital lattice geometry with the local interaction. The details on how $D^{s}_{\mu\nu}$ is obtained for our different models are discussed in the Supplementary Material (SM) sup .

In Fig. 2(e) we present $D^{s}$ as a function of $\max|\Delta|$ at $\nu=-2$ for both local (obtained with RM and DP) and RVB interactions (only RM). Figs. 2(a)–(b) and 2(e) depict a striking distinction between the local and RVB pairing schemes related to the pairing symmetry and the resulting form of $D^{s}$. The local pairing, yielding an $s$-wave symmetry, conserves the underlying $C_{3}$-symmetry of the TBG lattice [Fig. 2(a)] and $D^{s}$ is isotropic sup , i.e. $D^{s}_{xx}=D^{s}_{yy}$ and $D^{s}_{xy}=D^{s}_{yx}=0$. By contrast, the RVB pairing with strong enough interaction breaks the $C_{3}$-rotational symmetry and yields nematic pairing pattern in real space [Fig. 2(b)] which leads to an anisotropic response, i.e. $D^{s}_{xx}\neq D^{s}_{yy}$ and $D^{s}_{xy}=D^{s}_{yx}\neq 0$. The $s$-wave is gapped, whereas the nematic phase has nodal points in the moiré Brillouin zone [see also the inset of Fig. 2(e)]. The anisotropic $D^{s}$ results into an anisotropic kinetic inductance of TBG, and it can in principle be accessed via radio frequency impedance spectroscopy Chiodi et al. (2011) in a Hall-like four-probe setup.

As seen from Fig. 2(e), $D^{s}$ for the RVB interaction in the weak-coupling regime is still isotropic. This phase has the mixed $(d+id)+(p+ip)$ symmetry with a full energy gap, whereas the nematic phase of the flat band regime is identified as a mixture of $s$, $p$ and $d$-wave components Su and Lin (2018), with the $d$-wave being the dominant symmetry. Our results for the pairing symmetry are in agreement with Ref. Su and Lin, 2018 and they differ from the topological $d+id$ symmetry predicted in many TBG studies Guo et al. (2018); Xu and Balents (2018); Liu et al. (2018); Fidrysiak et al. (2018); Huang et al. (2019); Angeli et al. (2019); Wu (2019); Laksono et al. (2018); Lin and Nandkishore (2018); Classen et al. (2019) and also from other proposed symmetries which include $s$-wave Peltonen et al. (2018); Wu et al. (2018); Choi and Choi (2018); Po et al. (2018); Wu et al. (2019a); Laksono et al. (2018), extended $s$-wave Ray et al. (2019); Liu et al. (2018); Sherkunov and Betouras (2018); Angeli et al. (2019), $p$-wave Isobe et al. (2018); Wu et al. (2019a), $p+ip$-wave Roy and Juričić (2019), $d$-wave Wu et al. (2018); Isobe et al. (2018); Wu et al. (2019a), and $f$-wave Tang et al. (2019); Wu et al. (2019a); Lin and Nandkishore (2018); Classen et al. (2019). Apart from Ref. Su and Lin, 2018, nematic pairing has been predicted only in a few works Dodaro et al. (2018); Wu (2019); Wu and Sarma (2019); Venderbos and Fernandes (2018); Roy and Juričić (2019). The microscopic RM method allows to find the nematic pairing, unlike four-band models of Refs. Yuan and Fu (2018); Kang and Vafek (2018); Koshino et al. (2018).

BKT-transition temperature By computing $D^{s}$, one can determine $T_{\text{BKT}}$. In Fig. 3(a) we show $T_{\text{BKT}}$ as a function of $\max|\Delta|$. We can distinguish two qualitatively different regimes: in the weak-coupling limit the RVB and local interactions yield similar $T_{\text{BKT}}$ whereas for stronger interactions $T_{\text{BKT}}$ depends on the pairing model. Moreover, around $\max|\Delta|\gtrsim$2\text{,}\mathrm{m}\mathrm{e}\mathrm{V} the behaviour of the $T_{\text{BKT}}$ curves is almost linear, in accordance with previous studies Peotta and Törmä ([2015](#bib.bib24)); Julku *et al.* ([2016](#bib.bib25)) where $D^{s}$ of a flat band with the local interaction was shown to depend linearly on the pairing strength. In our case the narrow bands are not exactly flat but slightly dispersive and thus their flat band characteristics manifest only when the interaction strength is sufficiently large Peltonen *et al.* ([2018](#bib.bib49)). Because of this, we call the regime with $\max|\Delta|\gtrsim$2\text{\,}\mathrm{m}\mathrm{e}\mathrm{V} as the flat band limit. In this regime the DP and RM results are in agreement, whereas for weak interactions the results differ due to different band structures.

The difference of the two interaction schemes is further highlighted in Fig. 3(b) which presents the ratio $k_{B}T_{\text{BKT}}/\max|\Delta(T=0)|$. At the flat band limit this ratio approaches a constant whose value depends on the pairing potential. In experiments one can measure $T_{\text{BKT}}$ and in principle also deduce $\Delta$ (from the local density of states) and thus the ratio of these two quantities can be used to characterize the SC pairing observed in experiments.

In Figs. 3(c)–(d) we present $T_{\text{BKT}}$ as a function of $\nu$. The weak-coupling regime shows a dome-shaped structure of $T_{\text{BKT}}$ which reaches its maxima near the half-fillings of the hole- and electron-doped regimes, similar to experiments Lu et al. (2019). In the RM model the hole-doped region is much stronger due to higher density of states at negative energies [see Fig. 1(b)], while the DP model exhibits approximate electron-hole symmetry. Strong asymmetry of RM model is due to the applied rescaling approximation which amplifies the finite but small asymmetry of the unscaled model sup . In the flat band limit, the shape of the one-particle dispersions are, except for the pronounced particle-hole asymmetry of the RM model, completely dissolved.

Geometric contribution One can decompose $D^{s}$ to conventional, $D^{s}_{\textrm{conv}}$, and geometric, $D^{s}_{\textrm{geom}}$, parts, so that $D^{s}=D^{s}_{\textrm{conv}}+D^{s}_{\textrm{geom}}$ Peotta and Törmä (2015); Liang et al. (2017); sup . The conventional term depends on the inverse of the effective mass of the Bloch bands and is thus a single-band contribution, whereas $D^{s}_{\textrm{geom}}$ is a multiband effect depending on the overlap of the Bloch states and their momentum derivatives of the form $\innerproduct{\partial_{\mathbf{k}}n}{m}$, where $\ket{m}$ are the single-particle states of the $m$th Bloch band and $n\neq m$ Liang et al. (2017), i.e. $D^{s}_{\textrm{geom}}=0$ for a single-band system. For a strictly flat band, $D^{s}_{\textrm{conv}}=0$ so its superconductivity is purely a multiband process characterized by a finite $D^{s}_{\textrm{geom}}$. This arises an intriguing question related to TBG system: how much the interband terms between dispersive and narrow bands affect $D^{s}$ via $D^{s}_{\textrm{geom}}$?

We study this question in Fig. 4 for RVB [Fig. 4(a)] and local [Fig. 4(b)] pairing by presenting the total $D^{s}$ and its components. We further show results obtained by taking into account sup either only the $4$ flat bands or the 8 lowest (4 flat, 4 dispersive) bands labeled as $D^{s}_{4}$ and $D^{s}_{8}$, respectively. In both pairing cases the contribution coming from the 4 flat bands only is relatively small for larger interactions. The contribution of 8 bands is larger due to a larger $D^{s}_{\textrm{geom}}$, which is caused by the interband terms between dispersive and flat bands as the terms between the dispersive bands only are negligible. The slight dispersion of the narrow bands results in a finite $D^{s}_{\textrm{conv}}$. Note that $D^{s}_{\textrm{conv}}\approx D^{s}_{4,\textrm{conv}}$, i.e. $D^{s}-D^{s}_{4,\textrm{conv}}$ gives the total $D^{s}_{\textrm{geom}}$. From Fig. 4 we see that, at $\max|\Delta|\sim$1\dots 2\text{,}\mathrm{m}\mathrm{e}\mathrm{V}$$ i.e. when the system enters the flat band regime, $D^{s}_{\textrm{geom}}$ surpasses $D^{s}_{\textrm{conv}}$ for local pairing and becomes significant in the RVB case. An important implication of Fig. 4 is the importance of the dispersive bands when computing $D^{s}$ and the insufficiency of four band models. Also for a non-interacting system, higher bands have been argued to be necessary, but for different (symmetry) reasons Po et al. (2019).

Discussion Our work shows that TBG is characterized by two distinct superconducting regimes. When $\Delta$ is much smaller than the flat band bandwidth, the superfluid weight $D^{s}$ and the BKT transition temperature $T_{\text{BKT}}$ are well described by conventional theory of superconductivity. On the other hand, in this weak coupling regime the results are somewhat different for the RM and DP models. This is consistent with the low energy dispersion in TBG being very sensitive to the details of the model used Walet and Guinea (2019). In the flat band regime where $\Delta$ is larger than the width of the significant density of states in the flat bands, a major contribution to the superfluid weight $D^{s}$ originates from the geometric properties of the bands. The geometric contribution $D^{s}_{\textrm{geom}}$ is proportional to the quantum metric Peotta and Törmä (2015) whose importance in physics has been recently emerging Bengtsson and Życzkowski (2006); Gu (2010); Neupert et al. (2013); Roy (2014); Dobardžić et al. (2013); Peotta and Törmä (2015); Srivastava and Imamoğlu (2015); Gao et al. (2015); Julku et al. (2016); Tovmasyan et al. (2016); Piéchon et al. (2016); Liang et al. (2017); Törmä et al. (2018); Ozawa (2018); Lapa and Hughes (2019); Yu et al. (2018); Gianfrate et al. (2019); Asteria et al. (2019); Tan et al. (2019). Moreover, in the flat band regime, both $D^{s}$ and $T_{\text{BKT}}$ depend sensitively on the pairing mechanism, but not strongly on the employed microscopic model. In particular, for a non-local RVB interaction $D^{s}$ becomes anisotropic, which could be seen in four-terminal radio frequency spectroscopy experiments to reveal information about the pairing mechanism.

Within both of our models, at $\theta^{*}$ the crossover between the two regimes takes place for $\Delta=$1\dots 2\text{,}\mathrm{m}\mathrm{e}\mathrm{V}, implying $T_{\text{BKT}}\approx$1.5\dots 3\text{\,}\mathrm{K}. This is also the ballpark of the experimentally accessed critical temperatures Yankowitz et al. (2019); Lu et al. (2019). Thus the geometric contribution of the superfluid weight and the dependence on the pairing mechanism should be relevant for current experiments. An interesting future direction of research is to include other interaction channels than pairing and explore the insulating states observed in TBG Cao et al. (2018b); Yankowitz et al. (2019); Lu et al. (2019). Based on our results, one can anticipate that quantum geometry and multiband processes are important in superconductivity and correlated states of other twisted multilayer materials Chen et al. (2019); Shen et al. (2019); Liu et al. (2019a); Cao et al. (2019); Zhu et al. (2018, 2018); Zhang and Senthil (2018); Wu et al. (2019b); Zhang et al. (2019, 2019); Liu et al. (2019b); Lee et al. (2019); Schrade and Fu (2019); Classen et al. (2019); Scheurer et al. (2019).

Note added After submission of our manuscript, a related work Hu et al. (2019) appeared at the arXiv preprint server.