Prominent Cooper Pairing Away From the Fermi Level and its Spectroscopic Signature in Twisted Bilayer Graphene

We investigate phonon-mediated Cooper pairing in flat electronic band systems by solving the full-bandwidth multiband Eliashberg equations for superconductivity in magic angle twisted bilayer graphene using a realistic tight-binding model. We find that Cooper pairing away from the Fermi level contributes decisively to superconductivity by enhancing the critical temperature and ensures a robust finite superfluid density. We show that this pairing yields particle-hole asymmetric superconducting domes in the temperature–gating phase diagram and gives rise to distinct spectroscopic signatures in the superconducting state. We predict several such features in tunneling and angle resolved photoemission spectra for future experiments.

In the conventional theory of phonon-mediated superconductivity, as formulated by Eliashberg [1], Cooper pairing is expressed by the energy dependent function, ∆(ω), which is proportional to the pair's condensation energy [2]. When the electron energy scale dominates over the phonon one, pairing takes place only within a narrow energy window, ω c , around the Fermi surface where ∆(ω) > 0. If in addition the electron-phonon interaction (EPI) is very weak, one recovers the Bardeen-Cooper-Schrieffer (BCS) approximation, ∆(ω) ∆ × θ(ω − |ω c |) [3]. Flat band systems [4] however present a fascinating counterexample to this picture. These materials exhibit extremely narrow electronic bands near the Fermi level that host Van Hove singularities in extended regions of reciprocal space. The flatness of these bands enables a reversed situation where the phonon energy scale dominates over the electron one [5]. It is as yet unclear how the BCS picture applies to flat band superconductors. For example, a BCS approximation implies the peculiar absence of a robust Meissner effect unless extra geometrical terms contribute to the superfluid density [6].
In this work, we unravel the specifics of Cooper pairing in flat band superconductors by providing a full solution to the appropriate Eliashberg equations [7,8]. We focus on magic angle twisted bilayer graphene (TBG), a prototypical flat band superconductor, that has recently attracted tremendous research interest due to the rich variety of physical phenomena that result from a vastly changed electronic structure depending on the twist angle [9][10][11][12][13][14][15]. Most importantly for our present work, two nearly flat bands develop in close vicinity of the Fermi level [9,10,16]. When gated, TBG becomes superconducting at the magic twist angle ∼ 1.1 • with maximum critical temperature T c ranging between 0.5 K and 1.7 K [10,11]. There is a competition with external gating between the superconducting and an insulating state [10,11]. The latter is associated with enhanced correlations that manifest as extended features seen by tunneling experiments in the non-superconducting state of TBG [17][18][19]. Several works have recently emphasized the relevance of the EPI for superconductivity in TBG [20][21][22][23][24]. However, to our knowledge, all previous theories focused on effective descriptions near the Fermi level.
Here we present multiband, full-bandwidth Eliashberg calculations of superconductivity in TBG having as input realistic electron dispersions and EPI. We find that, in contrast to the BCS picture, ∆(ω) > 0 throughout the full bandwidth and show that Cooper pairing stems primarily from electrons away from the Fermi level. By calculating the full temperature-doping superconducting phase diagram, we show that this type of pairing not only contributes to T c , but also does this in a particlehole asymmetric way that agrees with experimental observations. Furthermore, our calculations demonstrate that in flat band superconductors, Cooper pairing away from the Fermi level ensures a finite superfluid density. Finally, we calculate angle resolved photoemission spectroscopy (ARPES) and scanning tunneling spectroscopy (STS) spectra and predict distinct low energy features below T c as signatures of such pairing.
To accurately describe the electronic properties of magic angle TBG, we adopt a ten band tight-binding model [14] that provides a faithful band structure, ξ l (k), that is consistent with findings of the continuum model [25][26][27][28] and ab initio calculations [15,29,30]. Here l is the band index and k the momentum in the mini Brillouin zone (BZ). A perturbation term in the Hamiltonian breaks the particle-hole symmetry and leads to a pair of nearly flat bands very close to the Fermi level, see Fig. 1(a). These two bands are energetically well separated from the rest by over 20 meV and have a very narrow bandwidth, W ≈ 7 meV, so that they harbor flat regions with extremely high DOS near the M-K and (-M)-(-K) symmetry lines [ Fig.1(b)-(c)]. As a result, the rest of the bands are marginally relevant and we hence confine ourselves here to these two flat bands in all following calculations. The inclusion of the complete ten bands has been carefully checked to not alter our results. As shown in Fig. 1(b), the two flat bands form Dirac points at the ±K points of the BZ and exhibit Dirac-like cones at their band extrema, located at Γ. Whereas focus has been put on the role of EPI involving intralayer modes [20][21][22][23], significant interlayer coupling has also been calculated [22]. Driven by recent experiments where superconductivity in TBG has been found to be tunable with the distance between the two graphene sheets [11], here we assume as the mediator of superconductivity the interlayer breathing mode with the characteristic phonon frequency of Ω = 11 meV [22,31]. This mode can provide an attractive interaction between electrons in an energy window that is ∼ 1.6×W. Given this fact and the similarity between the two dispersions, we expect that electrons from these bands couple to the mode similarly, regardless of their band index. Thus, we assume a global electron-phonon coupling strength, g 0 , which we also take as isotropic. By computing T c and requesting it to match to experimental values [10,11] our Eliashberg calculations yield g 0 = 1.5 meV, which is kept fixed henceforth.
To study our electron-phonon coupled system we employ full-bandwidth, multiband Eliashberg theory [7,32]. This theory takes explicitly into account scattering processes involving electrons with energies and momenta that are not restricted to the vicinity of the Fermi surface, therefore it goes beyond Migdal's theorem. As such, it also allows the possibility of Cooper pairing away from the Fermi level to be taken fully into account [7,8]. The matrix self-energy using Pauli matricesρ i , iŝ with Matsubara frequencies ω m =πT (2m + 1), temperature T and matrix Green's function, The chemical potential µ simulates gating effects. From the Dyson equation, we obtain a system of three coupled selfconsistent Eliashberg equations for the mass renormalization Z(iω m ), pairing function φ(iω m ) and chemical potential renormalization χ(iω m ) that are complemented by an equation for the particle filling, n (see Supplemental Material (SM)). The superconducting gap function has s-wave, spin singlet symmetry and is given by ∆(iω m ) = φ(iω m )/Z(iω m ). We first solve selfconsistently the multiband full-bandwidth Eliashberg equations for different values of (T, µ) and subsequently map the results to the corresponding sets of (T, n) values. The obtained Matsubara space results can be used to calculate any thermodynamic property that stems from the system's free energy. Here, this is done for the superfluid weight as will be discussed below. All quantities are then analytically continued from Matsubara to real frequencies by a selfconsistent procedure that is formally exact [33]. This gives us access to the precise real frequency dependent retarded Green's function at all temperatures and therefore to spectroscopic quantities of interest like the ones measured by STS and ARPES (see SM and Refs. [7,32]). All calculations are performed with the Uppsala Superconductivity (UppSC) code [34]. Fig. 2(a) shows our calculated temperature-filling phase diagram with respect to ∆(0). Here, n (0) is the electron filling of the system without gating/doping (µ = 0) and n (h) (n (e) ) are the reference points where the lower (upper) flat band is half-filled and the corresponding DOS peaks (see Fig.1(b)-(c)), which happens for hole or electron doping, i.e. µ h < 0 or µ e > 0. For hole doping (n < 0) we resolve a superconducting dome in qualitative agreement with experiment [9,11]. The counterpart for electron doping (n > 0) exhibits a lower T c as has also been measured [11]. Near half-filling (n (0) ) where insulating states are observed [9,11,35,36], we find T c = 0. This is not surprising since particle-hole pairing is not included in our model. Here we focus on the superconducting domes that have been observed experimentally and are well described by our theory.
With our chosen parameter set we find a maximum gap of ∆ ∼ 163 µeV and T c ∼ 1.05 K for for n (h) , whereas for n (e) , ∆ ∼ 154 µeV and T c ∼ 1 K. Clearly, T c is maximized when the doping level is such that the flat portion of the bands lies on the Fermi level, i.e. when the Fermi level DOS is maximized. This may be expected given that the coupling strength is proportional to the DOS [20,21]. A similar argument could also be conjectured to explain the electron-hole asymmetry in the shape of the two domes of Fig. 2(a), since the particle-hole asymmetry of the TBG bandstructure leads to DOS peaks that differ significantly in height, see Fig.1(c). Notably, such electron-hole asymmetry was not found in previous BCS calculations of the TBG phase diagram, despite the fact that the input bandstructure was particle-hole asymmetric [20,21].
To shed more light on this apparent discrepancy, we repeated our Eliashberg calculations for the same set of parameters but now taking as input only one of the flat bands at a time. In these single-band calculations we first keep only the upper band (red line in Fig.1(b)) and place the Fermi level at µ e and subsequently take only the lower band (blue line in Fig.1(b)) and place the Fermi level at µ h . The obtained results are shown in Fig. 2(b) with dashed (solid) lines for one (two) band calculations. Remarkably, the one band results are identical to each other despite the particle-hole asymmetry between the bands and the pronounced difference in Fermi level DOS. The superconducting gap and the critical temperature are reduced by approximately 40% as compared to the full calculation. It is only when we consider both bands together that T c is enhanced in a particle-hole asymmetric way. These results attest that at a given doping, the band not crossing the Fermi level is always involved in superconductivity, which points to the relevance of Cooper pairing away from the Fermi level in TBG. Such type of pairing has been recently predicted in FeSe/SrTiO 3 [7] where the phonon energy scale is comparable to the electron one.
Having calculated the precise real frequency dependence of all the quantities involved in our theory, we proceed by projecting the superconducting gap function on the underlying renormalized band structure,ξ l (k) [7]. The latter is a solution toξ = [ξ l (k) − µ + χ(ξ)]/Z(ξ) for each k and l. We then map ∆(ω) to ∆(ξ l (k)). Fig. 2(c) shows this projection for n = n (e) along high-symmetry lines. Remarkably, the maximum ∆ value is located at the largest energy away from the Fermi level. This behavior is markedly different from the conventional picture where ∆ < 0 above a given energy, due to the EPI becoming anti-pairing for electrons at such high energies [2]. In our case, ∆ stays positive (i.e. the EPI is pairing) throughout the full bandwidth. Thus, our results show that the Cooper pairing in TBG stems prevalently from electrons away from the Fermi level. Both bands, irrespective of which one crosses the Fermi level, contribute significantly to the Cooper pair formation and hence to T c .
Given the above results, it is now instructive to calculate the so-called London kernel, that describes the local, static current response to an applied transverse vector potential, J α = −Q αβ A β , with α, β = x, y. This is proportional to the superfluid density, n S (or superfluid weight), which in turn is related to the penetration depth via n S ∝ λ −2 . The first (second) summand in Eq. (3) is the diamagnetic (paramagnetic) part. Taking the usual approximations, Eq. (3) reduces to the standard energy-integrated isotropic Eliashberg form [37,38]. However, here we solve Eq. (3) as is. We also note that no extra geometrical (topological) term [6] is included here. The diamagnetic term is usually considered vanishingly small, so that one often focuses on the paramagnetic part which is proportional to ∆ 2 . For a flat band system, with an electronic dispersion of the general form ξ k ∝ ±(k/k f b ) N , where N 2 and k f b controls the extent of the flat band region [4], Eq. (3) would yield Q ≈ 0 if one naively assumes that ∆ → 0 for |ξ k | > ξ k f b . This is because the electron velocity is vanishingly small for |ξ k | < ξ k f b where ∆ would be expected to be large. Such a situation would appear problematic since then n S → 0 and λ → ∞ and it has been proposed that topological terms beyond Eq. (3) would become essential to produce n s = 0 [6]. On the other hand, for momenta away from the flat band region and therefore away from the Fermi level, the rapid increase in velocities would lead to a significant contribution in Eq. (3) if ∆ remains large for |ξ k | > ξ k f b . This is exactly what we find here, as shown in Fig.2(b),(c). To quantify this observation, we introduce the normalized superfluid weight D(ω c ) = |ξ k |≤ωc k in Fig.2(d), for small-ω c , where processes only at the Fermi surface enter into Eq. (3), D(ω c ) ≈ 0. However, for ω c ≥ max|ξ l (k)|, we retrieve the full calculation of Eq. (3) so that D(ω c ) = 1. Thus, there is no obstruction for having a finite superfluid density or a well defined Meissner effect in TBG or any similar flat band system. We now turn to the experimental signatures of our predicted Cooper pairing in TBG by using our obtained real frequency Green function to calculate STS and ARPES spectra (see SM). So far, no ARPES measurements exist for TBG and the existing STS studies have mainly focused on the non-superconducting spectral properties [13,[17][18][19][39][40][41]. Here we present results for n = n (e) . Results for n = n (h) are similar and included in the SM. Fig.3(a) shows the differential conductance, dI/dV ∝ DOS, for several temperatures below and above T c . At T > T c , the shape of the DOS is modified compared to Fig.1(c) due to the included EPI, but the relative height of the peaks remains similar. Below T c , the sharp peak at the Fermi level gives way to two coherence peaks separated by a gap, characteristic of s-wave superconductivity. The peak at negative bias that corresponds to the flat band region of the band below the Fermi level remains sharp. However, closer inspection reveals that it slightly moves to the left as temperature decreases. When we now consider the quasiparticle energy spectrum in the superconducting state, it becomes evident that this peak shift is a manifestation of having a non-decreasing ∆ away from the Fermi level. Fig.3(b) shows the renormalized tunneling spectrum that is the ratio of the spectra below, [dI/dV ] S , and above T c , [dI/dV ] N . Such experimentally accessible spectra are ideal for providing finer details of the modified tunneling due to superconductivity and therefore for identifying the signatures of our predicted Cooper pairing away from the Fermi level. Indeed, Fig.3(b) shows that apart from the coherence peaks around zero bias, there appears a hump-dip and a dip-hump for ω 0. Clearly, these structures occur as a result of the relative displacement between [dI/dV ] S and [dI/dV ] N due to Cooper pairing, see e.g. Eq. (4). Similar arguments apply when we look at frequency regions corresponding to the upper and lower edges of the electronic dispersion, shown in Fig. 3(c)-(d) . At the band edge, the DOS is minimal whereas ∆ acquires its maximum value. Therefore, the renormalized spectra exhibit sharp peaks there that should be experimentally discernible. We also find a less pronounced dip-hump near +0.5 meV [see Fig.3(b)] which could be experimentally detectable via a second-derivative analysis. This feature appears as the particle-hole symmetric replica of the structure near ω = −0.5 meV and is due to the twofold Bogoliubov quasiparticle energy dispersion in the superconducting state [cf. Eq. (4)].
The quasiparticle spectrum can be observed with ARPES that probes directly the momentum and energy resolved spectral function, A(k, ω). Our calculated A(k, ω) is shown in Fig.3(e) for n = n (e) at T T c . In Fig. 3(e), the Bogoliubov bands are clearly resolved to exhibit the characteristic "back-bending" near the Fermi level. Signatures of ∆ remaining large away from the Fermi level are not easy to discern here. However, we observe that the Dirac crossing at K (-K) has a Bogoliubov replica. In the case of hole doped TBG, this replicated Dirac crossing lies below the Fermi level (see SM) and, with sufficient resolution, can in principle be measured.
In conclusion, our full-bandwidth multiband Eliashberg calculations for phonon-mediated superconductivity in magic angle twisted bilayer graphene reveal that the Cooper pairing develops primarily away from the Fermi level and contributes significantly to the T c [7]. The existence of such Cooper pairs provides an alternative solution to the vanishing superfluid density paradox in flat band superconductors [6]. In addition, they give rise to nontrivial spectroscopic features at low energy and temperature that should be distinct from the ones already observed in the non-superconducting state [17][18][19]. Our predictions for the superconducting spectra provide a route for the unambiguous identification of Cooper pairing away from the Fermi level. Our findings should be qualitatively applicable for all phonon-mediated flat band superconductors [4,42,43].
We acknowledge support from the Swedish Research Council (VR), the Röntgen-Ångström Cluster, the K. and A. Wallenberg Foundation (Grant No. 2015.0060) and the Swedish National Infrastructure for Computing (SNIC).