Upper bounds on the superfluid stiffness and superconducting $T_c$: Applications to twisted-bilayer graphene and ultra-cold Fermi gases

Understanding the material parameters that control the superconducting transition temperature $T_c$ is a problem of fundamental importance. In many novel superconductors, phase fluctuations determine $T_c$, rather than the collapse of the pairing amplitude. We derive rigorous upper bounds on the superfluid phase stiffness for multi-band systems, valid in any dimension. This in turn leads to an upper bound on $T_c$ in two dimensions (2D), which holds irrespective of pairing mechanism, interaction strength, or order-parameter symmetry. Our bound is particularly useful for the strongly correlated regime of low-density and narrow-band systems, where mean field theory fails. For a simple parabolic band in 2D with Fermi energy $E_F$, we find that $k_BT_c \leq E_F/8$, an exact result that has direct implications for the 2D BCS-BEC crossover in ultra-cold Fermi gases. Applying our multi-band bound to magic-angle twisted bilayer graphene (MA-TBG), we find that band structure results constrain the maximum $T_c$ to be close to the experimentally observed value. Finally, we discuss the question of deriving rigorous upper bounds on $T_c$ in 3D.

Our work is motivated by the fundamental question: what limits the superconducting (SC) transition temperature T c ? Within BCS mean-field theory, and its extensions like Eliashberg theory, the amplitude of the SC order parameter is destroyed by the breaking of pairs, and T c scales with the pairing gap ∆. The material parameters that control the mean-field T c are the electronic density of states (DOS) at the chemical potential N (0) and the effective interaction, determined by the spectrum of fluctuations that mediate pairing.
Beginning with the pioneering experiments of Uemura [1] and theoretical ideas of Emery and Kivelson [2] on underdoped cuprates, it became clear that the mean field picture of T c scaling with the pairing gap is simply not valid in many novel superconductors. The loss of SC order is then governed by fluctuations of the phase of the order parameter, rather than the suppression of its amplitude, and T c is related to the superfluid stiffness D s . The material parameters that determine D s are rather different from those that determine the pairing gap ∆.
The question of mean field amplitude collapse versus phase fluctuation dominated SC transition is brought into sharp focus by a variety of recent experiments in narrow band and low density systems. One of the most exciting recent developments is the observation of very narrow bands in magic-angle twisted bilayer graphene (MA-TBG) leading to correlation-induced "Mott" insulating states [3] and superconductivity [4] in their vicinity. Flat bands are also also expected to arise in various topological states of matter; see, e.g., [5][6][7][8]. BCS theorybased intuition suggests that narrow bands have a large DOS N (0) and lead to high temperature superconductivity. Is this true or do phase fluctuations limit the T c ?
The extensive compilation of data in Fig. 6 of ref. [4] suggests that all known superconductors have a T c that scales at most like a constant times the "Fermi energy E F ", though there is considerable leeway in defining E F in strongly correlated and multi-band materials. We also note that ultra-cold Fermi gases in the strongly inter-acting regime of the BCS-BEC crossover [9,10] exhibit experimental values [11] of k B T c /E F larger than those observed in the solid state. All of these observations raise the question of ultimate limits on the T c of a superconductor or paired superfluid.
In this paper, we obtain sharp answers to these questions, especially in 2D. First, we derive an upper bound on the superfluid stiffness D s (T ) ≤ D(T ), where D is proportional to the optical conductivity sum rule. This inequality is valid in all dimensions and for arbitrary interactions. We then use the Berezinskii-Kosterlitz-Thouless (BKT) theory in 2D to obtain k B T c ≤ π D(T c )/2.
While the bound on T c is of completely general validity, it is most useful in the strongly correlated regime of narrow-band and low density systems, precisely where conventional mean-field approaches fails. We show that D is necessarily "small" in such systems, and, in many cases of interest, D is essentially determined by the (noninteracting) band structure.
We give several examples that illustrate the usefulness of our bounds for a variety of systems. For a single parabolic band we show that k B T c ≤ E F /8 in 2D. This exact result poses stringent constraints on the T c of the 2D BCS-BEC crossover in ultra cold atoms. We also describe bounds on T c for the 2D attractive Hubbard model, relevant for current optical lattice experiments [12], that demonstrate the tension between breaking of pairs and phase fluctuations, and highlight the connection with a pairing pseudogap [13,14].
Turning to multi-band systems, we use available band structure results [15][16][17][18][19] for MA-TBG to calculate D and thus constrain its T c without any assumptions about the pairing mechanism or order-parameter symmetry. We obtain a rigorous (but weak) bound of 15 K. Using physically motivated approximations, we estimate a bound on T c as low as 6 K.
Finally, we discuss the question of deriving similar bounds in 3D. We show that the presence of non-universal pre-factors in the relation between T c and D s , as well their scaling behavior near a SC quantum critical point, pose challenges in deriving a rigorous bound in 3D.
Results: We first outline our main results and then give a detailed derivation and specific applications. We consider a Fermi system described by the general Hamiltonian where k is crystal momentum, m is a band label, and σ the spin. H K is the kinetic energy and H int describes interactions (electron-phonon, electron-electron, etc.), including those that give rise to superconductivity. The external vector potential A enters H through a Peierl's substitution in the tight-binding representation of H K , but does not affect H int . For now, we ignore disorder and return to it at the end. The macroscopic superfluid stiffness D s determines the free energy cost of distorting the phase of the SC order parameter |∆|e iθ via the Boltzmann factor exp −D s d d r|∇θ| 2 |/2k B T . It is related to the London penetration depth via 1/λ 2 L = (4µ 0 e 2 / 2 )D s in 3D. Microscopically, D s can be calculated as the static, long wavelength limit of the transverse current response [20,21] to a vector potential. (Our results are equally valid for neutral superfluids with rotation playing the role of the magnetic field.) We obtain a rigorous upper bound valid in any dimension where Ω is the volume of the system and M −1 mm (k) is an inverse mass tensor that depends only on the electronic structure of H K ; see eq. (5) below. The temperature and interactions impact D only through c † kmσ c km σ , where the thermal average is calculated using the full H.
We next use D s to provide an upper bound on the SC transition temperature in 2D. We use the Nelson-Kosterlitz [22] universal relation to obtain For a weak coupling superconductor, T c is well described by mean field theory and our result, though valid as an upper bound, may not be very useful. On the other hand, as we show below, for a strongly interacting system the bound gives insight into both the value of T c and on its dependence on parameters.
Bound on superfluid stiffness: The intuitive idea behind D s ≤ D is as follows.
is the optical conductivity spectral weight integrated over the bands in eq. (1), and 4πe 2 / 2 D s is the coefficient of the δ(ω) piece in Re σ(ω) in the SC state; (note: ∞ 0 dω δ(ω) = 1/2). The inequality (2) says that the weight in the SC delta-function must be less than or equal to the total spectral weight.
To derive (2), we use the Kubo formula for D s as a linear response [20,21] to an external vector potential in an arbitrary direction a where D is the diamagnetic response ∼ δ 2 H/δA a 2 , while χ ⊥ is the transverse current-current correlation function. D is given by eq. (2) with Here α, β label orbitals/sites within a unit cell of a Bravais lattice, t αβ (k) is the Fourier transform of the hopping t αβ (r iα − r jα ), and U α,m (k) is the unitary transformation that diagonalizes t αβ (k) to the band basis m (k)δ m,m . The inverse mass tensor in eq. (5) also depends on the direction a = x, y, . . . through the derivative with respect to k a on the right hand side, however, we do not show this a dependence explicitly to simplify the notation. These results are derived in Appendix A, and the relation to the optical sum rule shown in Appendix B; see also ref. [23]. We next turn to the second term in eq. (4). From its Lehmann representation we see that χ ⊥ (q → 0, ω = 0) ≥ 0 at all temperatures; see Appendix C. We thus obtain For a single band system eqs.
(2) and (5) simplify greatly and we get D = (4Ω) −1 k,σ (∂ 2 (k)/∂k 2 a ) n σ (k), where the momentum distribution n σ (k) = c † kσ c kσ . This allows us to recover well-known special cases. (1) With nearest neighbor (NN) hopping on a square or cubic lattice, ∂ 2 (k)/∂k 2 a ∼ (k), and D is proportional to the kinetic energy. (2) A parabolic dispersion (k) = 2 k 2 /2m leads to the simple result D = 2 n/4m, independent of T and of interactions. Here D s (T ) = 2 n s (T )/4m and our bound simply says that the superfluid density n s (T ) ≤ n the total density.
For materials with non-parabolic dispersion and/or multiple bands, D depends on T and interactions. It is thus illuminating to derive a bound for D which depends only on the density. We describe the single band result here, relegating the multi-band generalization to Appendix D. We write H K = − Rδσ t(δ)c † R+δ,σ c R,σ + h.c. with translationally invariant hopping amplitudes t(δ) that depend only the vector δ connecting lattice sites R and R + δ. We couple the system to a vector potential and compute D, which involves terms like . We note that D ≥ 0, since it is the sum rule for Re σ(ω) ≥ 0. We then use the triangle inequality and Cauchy-Schwarz | c † i c j | ≤ n i n j = n to obtain D s ≤ D ≤ n δ δ 2 a |t(δ)|/2. This shows that for small hopping and/or low density, one necessarily has a small D s . T c bound in 2D: For a BKT transition in 2D, the T c and the stiffness D s are related by the universal ratio [22] , we then immediately obtain eq. (3). In an anisotropic system D depends on a = x, y through the ∂ 2 /∂k 2 a in eq. (5). We can use D = max D x , D y to obtain a bound on T c , however, we argue in Appendix H, We emphasize that eq. (3) with D(T c ) on the RHS is sufficient to derive the rigorous results below. However, to obtain the intuitively more appealing result k B T c ≤ π D(0)/2, we need to assume that D s (T ) is a decreasing . 2D Parabolic Dispersion: Consider a single band with (k) = 2 k 2 /2m with density n, so that the Fermi energy E F = π 2 n/m and arbitrary interactions that lead to pairing and superconductivity. Then M −1 (k) = m −1 and Ω −1 k,σ n σ (k; T ) = n independent of T and interactions, so that D = 2 n/4m. Eq. (3) then leads to the simple result which must be obeyed independent of the strength of attraction or order-parameter symmetry, provided the system exhibits a BKT transition. In a weak-coupling superconductor T c will actually be much smaller than E F /8 but, as we discuss next, the bound can be saturated in systems with strong interactions, such as the 2D BCS-BEC crossover experiments in ultra-cold Fermi gases.
2D BCS-BEC crossover: In ultra-cold Fermi gas experiments the two-body s-wave interaction between atoms is tuned using a Feshbach resonance. This has led to deep insights into the crossover [9, 10] from the weak coupling BCS limit with large Cooper pairs all the way to the BEC of tightly bound diatomic molecules. Asymptotically exact results are available in both the BCS and BEC limits, however, the crossover regime between the two extremes is very strongly interacting, with pair size comparable to the inter-particle spacing, and is much less understood. It is precisely here that our exact upper bound (6) is relevant.
The 2D crossover for s-wave pairing is parameterized by the dimensionless interaction [24] log(E b /E F ), where E b is the binding energy of the two-body bound state in vacuum and E F the Fermi energy. In the weak- , with a pre-factor that has been computed including the Gorkov-Melik-Barkhudarov (GMB) correction [25,26]. Clearly T c is much smaller than our bound.
In the BEC limit (E b E F ) the composite bosons have mass 2m, density n/2, and an inter-boson scattering length [27], which is valid in the regime log log 1. This too is smaller than our bound, though our exact result cautions against a naive extrapolation of the BEC limit result into the strong interaction regime. The results of the 2D Fermi gas experiment of ref. [28] seems to violate eq. (6) in the crossover regime. We note, however, that our bound is obtained for a strictly 2D system in the thermodynamic limit, while the experiment is on a quasi-2D system in a harmonic trap, from which it is difficult to accurately determine the BKT T c . The finite size of the trap raises T c ; even the non-interacting Bose gas in a 2D harmonic trap has a non-zero T c .
Magic angle twisted bilayer graphene: Let us next turn to a multi-band system of great current interest. The existence of very narrow bands in MA-TBG was predicted by continuum electronic structure calculations [15,16] that pointed out the crucial role of where θ is the twist angle between the two layers, w is the interlayer tunneling, v 0 F the bare Fermi velocity, and K the Dirac-node location in monolayer graphene. It was predicted that v F in TBG can be tuned to zero [15], with a bandwidth less than 10 meV by choosing certain magic angles θ, the largest of which ≈ 1.1 • has now been achieved in experiments [3,4]. Recently, pressure-tuning of w has also resulted in very narrow bands [29].
Little is known at this time about the nature of the SC state or the pairing mechanism, though the observed nonlinear I-V characteristics [3,4] are consistent with a BKT transition. Proximity to a "Mott" insulator and narrow bandwidth suggest the importance of electron correlations, while the extreme sensitivity of the dispersion to structure suggests that electron-phonon interactions could also be important. We argue here that simply using the available electronic structure information for MA-TBG, and without any prejudice about the interactions responsible for SC, we can put strong constraints on its superconducting T c .
There are two bands for each of the two valleys, one above and the other below the charge neutrality point (CNP) . Each band has a two-fold spin degeneracy, with bands for one valley related to those of the other by timereversal. We include these eight bands in the mm ,σ in eq. (2), while the k is over the moiré Brillouin zone, a hexagon with side 2K sin(θ/2) Kθ. We use the tight-binding model of ref. [17], a multi-parameter fit to the continuum dispersion [15], to calculate M −1 m,m (k) of eq. (5), which is block-diagonal in the valley index, so that there are no cross-valley terms in eq. (2).
We can obtain a more stringent T c bound if use fur- ther physical inputs. The "Mott" gap in the correlated insulator is experimentally [3,4] known be ≈ 0.3 meV, and we expect a superconducting gap which is at most that value. Thus we may assume that, at half-filling away from CNP on the hole doped side, say, the bands above the CNP are essentially empty and unaffected by pairing. Before proceeding, we derive a general result valid for arbitrary interactions which shows that inter-band terms do not contribute to eq. (2) for completely filled or empty bands. To prove this, we again use the Cauchy-Schwarz inequality | c † kmσ c km σ | 2 ≤ n mσ (k)n m σ (k) = 0 when either band m or m is empty. A similar argument works for the filled case after a particle-hole transformation; see Appendix E. Thus c † kmσ c km σ = 0 for m = m , whenever either of the two bands is completely filled or empty, and only m = m terms survive in eq. (2).
To bound T c for MA-TBG near half-filling on the holedoped side of the CNP, we take n m (k) = 0 for the empty bands above the CNP, as explained above. Keeping only band-diagonal terms and using the triangle inequality we obtain D ≤ ( 2 /4Ω) k,m,σ |M −1 mm (k)|n mσ (k). Using n(k) ≤ 1 for the bands below CNP we obtain the bound T c ≤ 14.4 K near half-filling for hole doping using the tight-binding model of ref. [17]. A similar calculation leads to T c ≤ 15.0 K near half-filling for electron doping; see Appendix F. We note that using M −1 and general constraints on n(k) leads to rigorous results, but weakens the bounds.
Finally, we make a physically motivated estimate of D, which yields an improved, but approximate, result. We use the T = 0 band theory result c † kmσ c km σ = δ m,m Θ (µ − m (k)), with the chemical potential µ determined by the density Ω −1 k,m,σ n mσ (k). This, together with M −1 mm (k) calculated from the tight binding model of ref. [17], leads to the density-dependent estimate of D plotted in Fig. 1. We note that using ∂ 2 /∂k 2 versus ∂ 2 /∂k 2 y to calculate M −1 affects our estimates by less than a percent.
The integrated optical spectral weight, given by 2πe 2 / 2 D, vanishes at the band insulators when all bands are either filled or empty. Clearly our bandstructure based estimate does not know about the "Mott" insulating states at half-filling away from CNP. (π/2) times the D plotted in Fig. 1 is an estimated upper bound on the SC T c . The system is not SC over most of the doping range, but our bound is the maximum attainable T c if the system were to exhibit superconductivity. We find the maximum T c to be about 6 K, while the experimental value is 3 K [29].
We note that the T c bounds are sensitive to the precise electronic structure results we use as input for calculating M −1 . As shown in Appendix F, using the tight binding results of ref. [18] for MA-TBG, leads to a T c estimate about 2.5 times higher than the one presented above, based on the band structure of ref. [17]. We emphasize that these differences arise from the fact that the details of the non-interacting band structure of MA-TBG are not very well established. Irrespective of that, our results suggest that MA-TBG is a strongly correlated SC in a phase fluctuation dominated regime. 2D attractive Hubbard model and optical lattices: We next obtain important insights on the value of T c and its interaction-dependence for the 2D attractive Hubbard model, where we can compare our bound with sign problem free Quantum Monte Carlo (QMC) simulations [30]. This system has also been investigated in recent optical lattice experiments [12].
Consider nearest-neighbor (NN) hopping on a square lattice with H = −t i,j σ c † i,σ c jσ + h.c. − |U | i (n i↑ − 1/2) (n i↓ − 1/2). For n = 1 the system has an s-wave SC ground state, exhibiting a crossover from a weak coupling BCS state (|U |/t 1) to a BEC of hard-core on-site bosons (|U |/t 1). The QMC estimate [30] of T c , obtained from the BKT jump in the D s , is a non-monotonic function of |U |/t at a fixed density n; see Fig. 2. The BCS mean field T MFT c correctly describes the weak coupling T c , (For a more accurate estimate, one should take into account the GMB correction [25] which suppresses the numerical pre-factor, but does not alter the functional form of T MFT c .) For |U |/t > 2, T MFT c is the scale at which pairs dissociate and lies well above T c . In the |U |/t 1 limit we see T c ∼ t 2 /|U |, the effective boson hopping.
Our bound permits us to understand T c (|U |/t) in the intermediate coupling regime where there are no other reliable analytical estimates. To estimate D analytically, we need to make an approximation for n(k). If we choose a step-function (as we did for the MA-TBG) we get T c ≤ 0.3t for n = 0.7, independent of |U |/t.
To obtain a better estimate, we note that, as |U |/t increases, the pair-size shrinks and n(k) broadens. In the extreme |U |/t-limit of on-site bosons, n(k) is flat (kindependent), leading to D → 0, since ∂ 2 /∂k 2 x is a periodic function with zero mean whose k-sum vanishes. To model this broadening of n(k), we use the results of the T = 0 BCS-Leggett crossover theory; see Appendix G. This gives us the (approximate) bound plotted in Fig. 2, which has the correct t 2 /|U | asymptotic behavior at large |U |.
In general, we see that For temperatures between the pairing scale T MFT c and T c at which phase coherence sets in, the "normal state" exhibits a pseudogap due to pre-formed pairs [13,14].
Three dimensional systems: Experiments suggest that there may be an upper bound on T c in 3D systems; see, e.g., Fig. 6 of ref. [4]. We have not succeeded in deriving a rigorous bound on the 3D T c , unlike in 2D. There are two challenges that one faces in trying to derive a bound in 3D, one related to rigorous control on numerical pre-factors and the other to the functional form of the relation between T c and D s . Both are related to the fact that in 3D the superfluid stiffness does not have dimensions of energy, unlike in 2D.
Following Emery and Kivelson (EK) [2], we focus on the 3D phase ordering temperature k B T θ = AD s (0) a, which could provide a bound on T c . Here A is a (dimensionless) constant and a is the length-scale up to which one has to coarse-grain to derive an effective XY model. EK use a 2 = πξ 2 , where ξ is the coherence length, and suggest, based on Monte Carlo results for classical XY models, that A 4.4 gave a reasonable account of experiments on underdoped cuprates and other materials.
However, the coefficient A is non-universal and can vary from one system to another. Consider the 3D problem of the BCS-BEC crossover in ultra-cold Fermi gases [10] with 2 k 2 /2m dispersion and interaction, characterized by the s-wave scattering length a s , tuned using a Feshbach resonance. At unitarity (|a s | = ∞), the experimental k B T c 0.17E F [11], while QMC estimates [31,32] range from k B T c 0.15E F − 0.17E F . QMC shows the expected non-monotonic behavior of k B T c /E F as a function of 1/k F a s , with a maximum k B T c /E F 0.22 at a small positive 1/k F a s . The maximum value of k B T c /E F is larger than the non-interacting BEC result, consistent with the rigorous result [33] that repulsive interactions increase the T c of a dilute Bose gas in 3D.
We choose ξ k −1 F near unitarity [34] and try to use k B T θ = A( 2 n/4m)( √ πξ) as a bound on T c . Consistency with the observed k B T c /E F 0.22 then requires A 7.4, quite different from the 4.4 quoted above. We do not know if there is a definite value of A that would give a "phase-ordering" upper bound on T c in 3D.
The following argument suggests that there may, in fact, be no general bound on T c that is linear in D s (0) in 3D. From a practical point of view, one is interested in learning about the highest T c in a class of materials. But, if a general bound were to exist, it should be equally valid in situations where both T c and D s (0) are driven to zero by tuning a (dimensionless) parameter δ → 0 + toward a quantum critical point (QCP). From the action S = 1 2 D s β 0 dτ d d r|∇θ| 2 + . . . describing the phase fluctuations of the SC order parameter, we get the quantum Josephson scaling relation [35] D s (0) ∼ δ (z+d−2)ν . One also obtains, as usual, T c ∼ δ zν , where z and ν are the dynamical and correlation length exponents in d spatial dimensions. Thus T c ∼ [D s (0)] z/(z+d−2) near the QCP. In 2D, this gives a linear scaling between T c and D s (0). However, in 3D we get T c ∼ D s (0) z/(z+1) which, sufficiently close to the QCP, will necessarily violate an upper bound on T c that is conjectured to scale linearly with D s (0). This is not just an academic issue, as experiments see precisely such a deviation from linear scaling with T c ∼ D s (0), consistent with z = 1, both in highly underdoped [36,37] and in highly overdoped [38,39] cuprates.
Concluding remarks: We have thus far ignored disorder. We note that D s of the pure system is necessarily larger than that in the disordered system. This can be seen by generalizing Leggett's bound [40] on the superfluid density (derived in the context of supersolids) to the case of disordered systems [41]. Thus our upper bounds for translationally invariant systems continue to be valid in the presence of disorder, although they can be improved.
Although we have focused on narrow band and low density systems here, our bounds have also important implications for systems close to insulating states, either correlation-driven or disorder-driven. In either case, if there is a continuous superconductor to insulator transition, the superfluid stiffness will eventually become smaller than the energy gap and control the SC T c .
As a design principle, it is interesting to ask if one can have multi-band systems where a narrow band has a large energy gap and large "mean field" T c interacting with a broad band that makes a large contribution to the superfluid stiffness, thus getting the best of both worlds. where we use the notation R = (r iα + r jβ )/2 and r = r iα −r jβ for simplicity. Since we are eventually interested in the long wavelength limit q → 0, we choose a very slowly varying vector potential and write riα r jβ Within linear response theory we can Taylor expand the exponential retaining terms which are linear (paramagnetic) and quadratic (diamagnetic) in A. We transform to Fourier space using t αβ (k) = r t αβ (r)e −ik·r and c iα = Ω −1/2 k e ik·riα d kα . We can then write the current operator j x = δH K /δA x as the sum of the paramagnetic (P ) and diamagnetic (D) current operators given by where we only show the x-component for simplicity. Note that the paramagnetic current operator, when transformed to the band basis, will in general have interband matrix elements [7,8]. The only property of j P x (q) that we will need to use below, however, is that it is a Hermitian operator; see equation (C1).
The superfluid stiffness D s is defined as the static longwavelength limit of the transverse response of the current density j to a vector potential A with q x = 0, q ⊥ → 0, ω = 0 (A6) and ⊥ represents the orthogonal directions to x. Standard linear response theory leads to the Kubo formula where the first term is the diamagnetic term, which is of central interest in this work, and the second is the transverse paramagnetic current-current correlation function. We will focus on the latter in Appendix C, where we show that χ ⊥ jxjx ≥ 0 at all temperatures. Here we focus on the first term that can be read off from the form of the diamagnetic current operator. We find it convenient to write it in the band basis as with the inverse mass tensor given by The unitary transformation U that transforms from the orbital to the band basis is defined by This allows us to write the final result in the band basis using We note several important points about the inverse mass tensor M −1 mm (k). (i) It depends only on the bare band structure, and is independent of temperature and interactions, (ii) it has both diagonal and off-diagonal terms in the band indices. and (iii) it is not simply related to the curvature of the bands ∂ 2 m (k)/∂k 2 x , in contrast to the single-band case in equation (A12).
The standard reference on the formalism for calculating the superfluid stiffness in lattice systems is Scalapino, White and Zhang (SWZ) [21]. Our normalization conventions differ from them and, more importantly, they focus on the special case of a single band model with nearest-neighbor (NN) hopping on a square (or cubic) lattice. Thus it may be useful for us to provide a "dictionary" relating our results to theirs.
In the single-band case our expression for D reduces to where the momentum distribution This result is valid for arbitrary one-band dispersion. For the special case of nearest-neighbor (NN) hopping on a square (or cubic) lattice, it is easy to see that the right hand side of equation (A12) is proportional to the kinetic energy in the x-direction, −K x in the notation of SWZ. Our result thus reduces to Finally, we note that our superfluid stiffness D s is related to that of SWZ by Appendix B: Relation between D and optical spectral weight To see that D is proportional to the optical sum rule spectral weight, we identify the dynamical conductivity σ(ω) as the current response to an electric field E = Using the Kramers-Krönig relation and Reχ jxjx (ω → ∞) → 0, we obtain the sum rule for the optical conductivity as A similar argument for completely filled bands follows from a particle-hole transformation Thus we conclude that for filled and empty bands, the inter-band terms do not contribute to the sum in equation (E1), even in the presence of arbitrary interactions.
Finally, we note the simple fact that within band theory there are no inter-band contributions to D. In the absence of interactions (denoted by subscript 0) we obtain where f is the Fermi function. Magic angles in twisted bilayer graphene were first predicted by the continuum model [15]. Following up on the experimental discovery of correlation-induced insulators and superconductivity in MA-TBG, there has been considerable progress in understanding its electronic structure [17][18][19]. We first focus on the bounds that we obtain from the tight binding model of Koshino et. al. [17], and then at the end of the Appendix compare these with the results we obtain from the tight binding model of Kang and Vafek [18].
The continuum model dispersion [15] is accurately reproduced by the multi-parameter tight binding fit of Koshino et. al. [17] (see Fig. 3) which takes into account hopping over distances up to 9|L M | where L M is the moire lattice vector. We use the hopping integrals presented in the Supplementary Information file eff hopping ver2.dat of ref. [17] to construct the noninteracting Hamiltonian H K of equation (A2). We then identify the unitary matrix U (k) that diagonalizes t αβ (k) (see equation (A10)) and use it together with t αβ (k) to compute the inverse mass tensor Note that we have made explicit here the direction a = x, y as an additional subscript on M −1 .
The inverse mass tensor, obtained from the band structure information as described above, is used to compute D x and D y and bound T c as described in the paper. The additional input needed to determine D using equation (A8) is c † km c km , and we took two different approaches to compute this.
In the first approach, we looked at SC near half-filling on the hole-doped side of the CNP, and argued that the chemical potential was sufficiently far from the CNP that we can take the band above the CNP to be empty. Then using the result of Appendix E we can ignore all interband terms with m = m . For the occupied band we only used the general constraint that n(k) ≤ 1. Using the triangle inequality, we then obtain where the empty bands above the CNP are excluded from the sum. A similar reasoning also works for SC in the vicinity of half-filling on the electron-doped side of the CNP, where we need to use the fact that the bands below CNP are filled to eliminate inter-band terms following Appendix E. We use a particle hole transformation c mk → h † mk , under which t αβ (k) → −t αβ (k) and thus M −1 → −M −1 . We write D in terms of the hole momentum distribution functions n h m (k) = h † mk h mk to get  4. Comparison of (a) the band structure and (b) the integrated spectral weight D for the models in ref. [17] (in black) and ref. [18] (in red).
We have first used m U β,m (k)U † m,α (k) = δ β,α , which follows from the unitarity of U , and then the fact that ∂ 2 t αα (k)/∂k 2 a is a periodic function with zero mean, whose k vanishes. Using the triangle inequality and the general constraint n h (k) ≤ 1, we obtain an expression for electron doping which is similar to the hole-doped case: where now the filled bands below the CNP are excluded from the sum. These bounds, though rigorous, are weak because they involve |M −1 | and only very general constraints on n(k). The second (approximate) approach was to simply use a T = 0 (non-interacting) band-theory estimate. We thus use equation (E4) to obtain with the chemical potential µ determined by the density. We found that D x and D y calculated from the tight binding model of ref. [17] differ by less than a percent. The resulting density-dependent D is shown in Fig. 1 of the main paper. We note that there are many different tight binding models for describing the narrow bands in MA-TBG and our T c bounds depend on this input. We have focused above on the results based on ref. [17] with an electronic structure that has separate charge conservation at the K and K valleys. A rather different model without valley-charge conservation was derived [18] using only time-reversal and point group symmetry. We compare in Fig. 4(a) the band structures of ref. [17] in black and that of ref. [18] in red. The corresponding integrated spectral weights D are shown in Fig. 4(b) using the same color convention. The maximum T c based on the band structure of ref. [18] is 15 K, which is 2.5 times larger than that estimated from ref. [17]. strong coupling regime, and less useful in the weak coupling regime where T c is, in fact, well described by T MFT c , the pair breaking energy scale.