Critical magnetic fields and electron-pairing in magic-angle twisted bilayer graphene
CCritical magnetic fields and electron-pairing in magic-angle twisted bilayer graphene
Wei Qin, Bo Zou, and Allan H. MacDonald Department of Physics, The University of Texas at Austin, Austin, Texas 78712,USA ∗ (Dated: February 23, 2021)The velocities of the quasiparticles that form Cooper pairs in a superconductor are revealed bythe upper critical magnetic field. Here we use this property to assess superconductivity in magicangle twisted bilayer graphene (MATBG), which has been observed over a range of moir´e bandfilling, twist angle, and screening environment conditions. We find that for pairing mechanisms thatare unrelated to correlations within the MATBG flat bands, minima in an average Fermi velocity v ∗ F ≡ k B T c ‘ c / ~ , where ‘ c is the magnetic length at the critical perpendicular magnetic field, arealways coincident with transition temperature maxima. Both extrema occur near flat-band vanHove singularities. Since no such association is present in MATBG experimental data, we concludethat electronic correlations that yield a band-filling-dependent pairing glue must play a crucial rolein MATBG superconductivity. Introduction—
The observation of superconductingdomes near correlated insulating states in magic-angletwisted bilayer graphene (MATBG) [1, 2] has stimulatedinterest in achieving a full microscopic understanding ofthis relatively simple electronic system [3–19]. The elec-tronic properties of MATBG devices are extremely sen-sitive to multiple tuning knobs, especially electrostaticdoping and twist angle [1, 2, 9], but also interlayer sep-aration [3, 20], vertical displacement field [3], and three-dimensional screening environment [21–23]. The tunabil-ity of MATBG makes it a particularly appealing experi-mental platform for the exploration of strong-correlationsuperconductivity.As often happens in superconductors, observationsthat clearly distinguish between purely electronic pairingmechanisms, possibly related to the correlated insulatingstates [24–35], and conventional phonon-mediated elec-tron pairing [36–38] are sparse. In this Letter, we demon-strate by explicit calculations for representative flat-bandand optical and acoustic phonon-mediated pairing mod-els that if superconductivity were mediated by a pairingglue derived from non-electronic degrees-of-freedom, themaximum critical temperature would always occur whenthe Fermi level is close to the magic-angle flat-band vanHove singularity (VHS). The band filling factor at whichthe VHS appears in a particular system is extremely sen-sitive to electronic structure model details and also tosample-dependence, particularly in connection with bro-ken flavor symmetries and associated Fermi surface re-constructions [39–42]. It is therefore not easily predictedtheoretically and is in all likelihood device dependent. Aswe show here however, it can be determined experimen-tally by combining measured critical temperatures andcritical perpendicular magnetic fields to determine thetypical velocity of the quasiparticles that form Cooperpairs. Using this strategy we argue that there is no ex-perimental correlation between a high transition temper-ature and van-Hove proximity, which argues against con-ventional phonon-mediated pairing.Table I summarizes experimental data for MATBG
TABLE I. Experimental results for the critical temperature T c , critical perpendicular magnetic field H c , Pauli criticalfield H P , and the extracted average Fermi velocity v ∗ F forMATBG in a variety of experiments with different twist angles θ , band fillings ν , and screening environments. θ ( ◦ ) ν T c (K) H c (mT) H P (T) v ∗ F (10 m/s) Refs1.05 -2.02 1.7 - 3.15 - [1]1.16 -2.15 0.5 125 0.93 0.341.27 -2.33 3.1 210 5.74 1.54 [3]1.27 -2.62 2 72 3.7 1.78-2.31 3.1 180 5.74 1.67 [9]1.1 -1 0.14 100 0.26 0.100.67 0.16 300 0.3 0.071.48 0.65 400 1.2 0.231.15 -1.6 0.92 220 1.7 0.46 [21]1.8 0.42 26 0.78 0.561.04 -2.43 1.3 >
50 2.41 < >
60 1.3 < >
50 7.4 < superconductors. Each observation determines an av-erage Fermi velocity defined by v ∗ F ≡ k B T c / ~ q c , where q c = 1 /‘ c is the maximum pairing momentum [43],2 π‘ c = Φ /H c , and the superconducting flux quan-tum Φ = 2 e/hc . The average Fermi velocity definedabove is motivated by the observation that superconduc-tivity is suppressed at finite pairing momentum becausethe electrons that form Cooper pairs differ in energyby ∼ ( d(cid:15) k /dk ) /‘ ∼ ~ v F /‘ , and that it is lost in mean-field theory when this difference is comparable to k B T c .Microscopically, superconductivity is lost in a magneticfield because Landau level Cooper pair states are formedfrom individual electron states that differ in energy by ∼ ~ v F /‘ c [48].We will show by explicit calculation that this experi-mentally available average Fermi velocity very effectivelyidentifies VHS in two-dimensional band structures. Wenote that the experimental values of v ∗ F in Table I aretypically a hundred or more times smaller than the Dirac a r X i v : . [ c ond - m a t . s up r- c on ] F e b -0.2 -0.1 0 0.1 0.2 q x ( m e V ) = -1 -0.2 -0.1 0 0.1 0.2 q x -0.07-0.06-0.05-0.04-0.03-0.02-0.010 F ( m e V ) = -1.0 = -1.5 = -2.0 = -2.5 = -1.0 = -1.5 = -2.0 = -2.5 q
1, where the MBZ is indicated by the dashed hexagon.These results are obtained at zero temperature with the inter-layer tunneling ratio set to η = t AA /t AB = 0 .
85, and Coulombrepulsion u = 40 meV · nm chosen to yield T c ∼ . velocities of isolated graphene sheets, demonstrating thecrucial role of the dramatically flattened moir´e bands[49]. The main point we wish to make here however isthat experimental critical temperature maxima, whichalways occur in a narrow range of filling factor near ν = − . ν . Below we first describe mi-croscopic calculations of T c and H c based on a phonon-mediated pairing model that demonstrate this point, andthen discuss the implications of this observation. Ourmain results are summarized in Fig. 1 and Fig. 2. Critical magnetic field from finite-momentumpairing—
The conclusions in this paper are basedon calculations of the pairing wavevector dependence ofthe condensation energy for MATBG superconductors. We assume that superconducting condensation energycan be calculated using a mean-field approximation,and that the relevant Cooper pairs involve two electronsfrom opposite valleys. Given these assumptions, the the-oretical system properties depend only on the MATBGband structure model and the interaction Hamiltonian.The former is not accurately known at present, mainlybecause bands are renormalized by interactions [39–41],and because these renormalizations are sensitive tothe three-dimensional screening environment [50, 51].We will account for the possible role of band-structurevariations by varying η , the ratio between the intra-sublattice and intersublattice tunneling terms in theBistritzer-MacDonald MATBG model [49].In the phenomenological Ginzburg Landau theory ofan isotropic superconductor [43], the superconductingcondensation energy at finite pairing-momentum δF ( q ) = δF [1 − ( q /q c ) ] , (1)where δF is the condensation energy at zero pairing-momentum and q c is the critical pairing wavevector towhich superconductivity survives. Our microscopic cal-culations are in close agreement with this expression andcan be fit to determine δF and q c . As explained in theSupplemental Material (SM) [43], the critical perpendic-ular magnetic field is related to q c by H c = Φ s q c / π . Inaddition, the supercurrent density as a function of pair-ing wavevector can also be calculated from δF ( q ) using j = (2 e/ ~ )[ ∂F ( q ) /∂ q ] [43]. Mean-field theory for finite-momentum pairing—
Thefinite pairing-momentum mean-field calculations that wehave performed for MATBG superconductors employmodels in which electrons are paired by in-plane opticalphonons [36, 37, 43] competing with repulsive Coulombinteractions. Because the flat band width is small com-pared to the phonon energy, both interactions are essen-tially instantaneous. Optical phonons are a convenientchoice for a putative pairing mechanism because theirinteractions with graphene π -bands are well understood[36, 52]. The case of acoustic-phonon mediated interac-tions is discussed later.The long-range Coulomb interaction is thought tobe responsible for flavor-symmetry-breaking, also dis-cussed later, which plays an important role in MATBGphysics [39–41]. Here we assume that there is neg-ligible Coulomb scattering between valleys and adjustthe strength of repulsive intra-valley Coulomb scatter-ing, which plays a depairing role, to match experimentalcritical temperature. (This adjustment however sets theCoulomb scale to value that is smaller than the expec-tation from a crude estimation. Coulomb repulsion isassumed to have the same strength for electrons fromthe same valley and opposite valleys.) Specifically, themodel electron-electron interaction Hamiltonian H ee =2 u P ll P τsσσ R d r ρ lτsσ ( r ) ρ l ¯ τ ¯ sσ ( r ) , where the densityoperator ρ lτsσ ( r ) = ψ † lτsσ ( r ) ψ lτsσ ( r ), ψ lτsσ ( r ) is the -2 -1 0 1 2 (meV) H c ( T ) -2 -1 0 1 2 (meV) v F * ( m / s ) DO S ( e V - n m - ) VHS -2 -1 0 1 2 (meV) T c ( K ) (a)(b) (c) FIG. 2. (a)-(c) Zero-pairing-momentum critical temperature T c , perpendicular critical magnetic field H c , and averageFermi velocity v ∗ F as functions of chemical potential µ for sev-eral strengths of Coulomb repulsion u in units of meV · nm .These results are calculated for a MATBG band-structuremodel with θ = 1 . ◦ and η = 0 .
7. The dotted lines indi-cate the values of µ for flat-band VHSs in the DOS, which isplotted as a dashed curve in (c). v ∗ F in (c) has a minimumwhere T c in (a) is maximized. real-space annihilation field operator, and l, τ, s, σ = ± H BdG ( q , k ) = (cid:20) H ( k ) − µ q ∆ q ( k )∆ † q ( k ) −H T ( q − k ) + µ q , (cid:21) , (2)where each block acts on four-component sublatticespinors and the pair potential for pairing wavevector q [∆ q ( k )] αα = X α α X k [ V ( k , k )] αα α α [ F q ( k )] α α . (3)Here V ( k , k ) is the total interaction matrix element in-cluding phonon-mediated attraction and Coulomb repul-sion contributions, F q ( k ) is Gorkov’s anomalous Green’sfunction [53], and the summation of k is over the MBZ.The chemical potential µ q is determined self-consistentlyby particle number conservation. The free energy of thefinite pairing-momentum superconducting state F s ( q ) = C q + An µ q + 12 β Tr X k ln f [ − E q ( k )] , (4)where n is the carrier density measured from chargeneutrality, A is the sample area, f ( (cid:15) ) is the Fermi-Diracdistribution function, E q ( k ) are the eigenvalues of theBdG Hamiltonian, and C q = − Tr F † q V F q is a double-counting correction [43]. The superconducting condensa-tion energy is defined as δF ( q ) = F s ( q ) − F n , where F n is the normal-state free energy calculated by Eq. (4) forvanished pair potential. We use numerical results for the condensation energy as a function of pairing wavevector,band filling ν , and model parameters to connect withexperimental observables. Robust correlations between velocity minima and criti-cal temperature maxima–
Because our model interactionis local the real-space pair potential is conveniently para-materized by performing a Fourier expansion. The coef-ficients of reciprocal lattice vectors in this expansion areplotted as a function of pairing wavevector in the MBZ inFig. 1(b). In the insert of Fig. 1(b), the asymmetric be-havior of ∆ P Q is due to rotational symmetry breakingat finite-momentum pairing. Figures 1(c) and (d) de-pict the pairing wavevector dependence of δF ( q ), whichpossesses a minimum at q = 0, indicating that the zero-momentum pairing state is the ground superconductingstate. Moreover, δF ( q ) is nearly isotropic and vanishesat a critical pairing wavevector q c . The small value of q c compared to the MBZ implies that the superconductingcoherence length is many times of the moir´e period.Figures 2 (a)-(b) show mean-field results for the zero-pairing-momentum critical temperature T c , and the criti-cal perpendicular magnetic field H c obtained from finite-pairing-momentum calculations. As expected, T c de-creases monotonously upon strengthening the repulsiveCoulomb interaction between electrons. Dome-like fea-tures are revealed in T c as a function of µ which arereminiscent with the superconducting domes, but are as-sociated with maxima centered on the flat-band VHSs. H c possesses sharp peaks at VHSs consistent with thephenomenological argument that H c is proportional tothe effective mass of paired electrons [43]. Away from theVHSs, H c decreases quickly and approaches zero when T c becomes zero. Figure 2(c) shows the average Fermivelocity extracted from the mean-field results of T c and H c via v ∗ F = k B T c / ~ q c . The average velocities are ahundred or more times smaller than the Dirac velocityof isolated graphene sheets v D ∼ m/s, in agreementwith experiment. In contrast to experiment, v ∗ F possessesminima at the VHSs and therefore at the T c maxima.In Fig. 2(c), we have seen that the average Fermi veloc-ity is almost perfectly a band structure property that isindependent of the pairing strength. In Fig. 3, we studythe correlation between T c and v ∗ F by varying the band-structure parameters η and θ for fixed pairing strength.Over the illustrated parameter range the flat-band widthis larger at smaller η and larger θ [43]. Accordingly de-creasing η or increasing θ tends to reduce T c as shownin Fig. 3 (a) and (d), where dome-like features in T c are peaked around the shifting VHSs. As illustrated inFig. 3 (b) and (e), sharp peaks of H c are also foundat the VHSs, and the magnitude of H c is dramaticallysuppressed for models with larger flat-band width. Theextracted v ∗ F are plotted in Fig. 3 (c) and (f), and areas expected larger for larger flat-band widths. For ev-ery band-structure model, v ∗ F possesses a V-shaped min-imum near the VHS, which is always coincident with the -5 -4 -3 -2 -1 0 1 2 3 4 5 (meV) H c ( T ) -5 -4 -3 -2 -1 0 1 2 3 4 5 (meV) T c ( K ) -5 -4 -3 -2 -1 0 1 2 3 4 5 (meV) v F * ( m / s ) (d) (e) (f) -2 -1 0 1 2 (meV) H c ( T ) -2 -1 0 1 2 (meV) v F * ( m / s ) -2 -1 0 1 2 (meV) T c ( K ) (a) (b) (c) FIG. 3. (a)-(c) T c , H c , and v ∗ F as functions of µ for severalvalues of η for MATBG with θ = 1 . ◦ and u = 30 meV · nm .(d)-(f) T c , H c , and v ∗ F as functions of µ for several differenttwist angles for MATBG with η = 0 . u = 20 meV · nm .Vertical arrows in (c) and (f) highlight the positions of flat-band VHSs. maximum of T c . Discussion—
So far we have neglected both the possi-bility of flavor symmetry breaking and the role of Zeemancoupling to the electronic spin. Indeed, this assumptioncan be questioned since there is strong experimental ev-idence that the strongest superconducting dome occursin a region of filling factor where only two flavors arepartially occupied and the moir´e flat-bands have conse-quently reconstructed [1–9, 11, 14, 15, 39–42]. If we wereto assume that the superconducting state near ν = − . H P ≈ . T c (T) [1], which is much larger than H c for large T c , indicating the orbital effect dominates pair-breaking. Figure 4(a) compares critical magnetic fieldscalculated by including only the Zeeman effect, only theorbital effect, and both effects. Close to the VHS, H P iscomparable to H c because the orbit effect is suppresseddue to small Fermi velocities. Away from the VHSs, how-ever, H P is much larger than H c , consistent with experi-mental observations. The two peaks in H P near the VHSreflect Zeeman splitting. The critical fields calculated in-cluding both effects show that H c is nearly unchangedby Zeeman coupling except when the Fermi level is very -2 -1 0 1 2 (meV) T c ( K ) (meV) F () (b) -2 -1 0 1 2 (meV) B ( T ) (a) FIG. 4. (a) Critical magnetic field calculated at zero pairingwavevector by including Zeeman coupling ( H P ), and calcu-lated from the critical pairing wavevector H c without (w/o)and with (w/t) Zeeman coupling. In the presence of Zeemancoupling, H c = (Φ s / π ) q x,c q y,c , where q x,c and q y,c are thecritical pairing wavevectors along x and y directions. Theseresults are for θ = 1 . ◦ , η = 0 .
7, and u = 30 meV · nm .(b) Electron-acoustic phonon coupling strength λ and theestimated T c for MATBG with θ = 1 . ◦ and η = 0 . α F ( ω ) at µ = − . close to the VHS.We have so far not explicitly included acoustic-phononmediated interactions, which may compete more suc-cessfully with Coulomb interactions because they areretarded. The longitudinal acoustic phonon veloc-ity in an isolated graphene sheet v ph = 2 × m/s is comparable to the flat-band electron veloci-ties. Figure 4(b) illustrates the band-filling depen-dence of the dimensionless electron-acoustic-phonon cou-pling constant λ = 2 R dωα F ( ω ) /ω and the non-trivial frequency dependence of the Fermi-surface aver-aged coupling strength α F ( ω ) for a representative band-structure model [43]. The critical temperatures shownin Fig. 4(b) are estimated by using the formula T c = ~ ω ln . k B exp [ . λ ) λ − µ ∗ (1+0 . λ ) ], where the averaged phonon fre-quency ω ln = exp [2 R dω ln ( ω ) α F ( ω ) / ( λω )] [54], andthe reduced Coulomb coupling strength µ ∗ = 0 .
2. Incontrast to the optical-phonon-mediated interaction, re-tardation does supply a formal justification for reducedCoulomb coupling [55], but since the phonon and elec-tronic energy scales are similar, it still does not justify thelarge reduction needed to match experimental T c scales.That aside, it is clear in Fig. 4(b) that the association of T c maxima with flat-band VHSs applies equally well toacoustic-phonon mediated superconductivity.For a given pairing glue, weak-coupling BCS theorypredicts that T c is positively correlated with the DOS atFermi level, determined by the Fermi surface size and thetypical quasiparticle velocity. In experiment [1, 2, 19],these quantities are often extracted from the frequencyand temperature dependence of magnetic oscillations,which respectively measure Fermi surface area and cy-clotron effective masses. These measurements are nor-mally carried out in the high magnetic field regime, wheresuperconductivity in MATBG is completely suppressed.In this study, we propose extracting average Fermi veloc-ity instead from a combination of the T c and H c , whichis more directly connected to the properties of quasiparti-cles that participate in pairing. We establish the reliabil-ity of this approach by calculating T c and H c for a widevariety of representative MATBG models. In particular,our theoretical calculations show that the definition of v ∗ F in this way are largely independent of pairing model,that they exhibit a V-shaped cusp at the flat-band VHS,and that that v ∗ F is always negatively correlated with T c .Since experimental values of v ∗ F , summarized in TableI and in the SM [43], do not show any indication of asuch correlation between v ∗ F and T c , we conclude thatthe pairing glue is related to electronic fluctuations thatare sensitive to flat-band filling factor and optimized nearthe peaks of the experimental superconducting domes. Acknowledgment—
The authors acknowledge helpfuldiscussions with Andrea Young and Ming Xie. This workwas supported by DOE grant DE-FG02-02ER45958 andWelch Foundation grant TBF1473 ∗ [email protected][1] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, Unconventional super-conductivity in magic-angle graphene superlattices, Na-ture , 43 (2018).[2] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken,J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T.Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-fillingin magic-angle graphene superlattices, Nature , 80(2018).[3] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watan-abe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean,Tuning superconductivity in twisted bilayer graphene,Science , 1059 (2019).[4] Y. Xie, B. Lian, B. J¨ack, X. Liu, C.-L. Chiu, K. Watan-abe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Spec-troscopic signatures of many-body correlations in magic-angle twisted bilayer graphene, Nature , 101 (2019).[5] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian,M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J.Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Maxi-mized electron interactions at the magic angle in twistedbilayer graphene, Nature , 95 (2019).[6] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora,R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. vonOppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge,Electronic correlations in twisted bilayer graphene nearthe magic angle, Nature Physics , 1174 (2019).[7] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule,J. Mao, and E. Y. Andrei, Charge order and broken rota-tional symmetry in magic-angle twisted bilayer graphene,Nature , 91 (2019). [8] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney,K. Watanabe, T. Taniguchi, M. A. Kastner, andD. Goldhaber-Gordon, Emergent ferromagnetism nearthree-quarters filling in twisted bilayer graphene, Science , 605 (2019).[9] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I.Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang,A. Bachtold, A. H. MacDonald, and D. K. Efetov, Su-perconductors, orbital magnets and correlated states inmagic-angle bilayer graphene, Nature , 653 (2019).[10] X. Hu, T. Hyart, D. I. Pikulin, and E. Rossi, Geometricand conventional contribution to the superfluid weight intwisted bilayer graphene, Phys. Rev. Lett. , 237002(2019).[11] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu,K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young,Intrinsic quantized anomalous Hall effect in a moir´e het-erostructure, Science , 900 (2020).[12] A. Julku, T. J. Peltonen, L. Liang, T. T. Heikkil¨a,and P. T¨orm¨a, Superfluid weight and Berezinskii-Kosterlitz-Thouless transition temperature of twisted bi-layer graphene, Phys. Rev. B , 060505(R) (2020).[13] F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Topology-bounded superfluid weight in twisted bilayer graphene,Phys. Rev. Lett. , 167002 (2020).[14] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S.Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, andA. Yazdani, Cascade of electronic transitions in magic-angle twisted bilayer graphene, Nature , 198 (2020).[15] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R.Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. vonOppen, A. Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani,Cascade of phase transitions and Dirac revivals in magic-angle graphene, Nature , 203 (2020).[16] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young,Superconductivity and strong correlations in moir´e flatbands, Nature Physics , 725 (2020).[17] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel,and A. Vishwanath, Charged skyrmions and topologi-cal origin of superconductivity in magic angle graphene,arXiv:2004.00638.[18] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu,A. H. MacDonald, B. A. Bernevig, and D. K. Efetov,Competing zero-field Chern insulators in superconduct-ing twisted bilayer graphene, arXiv:2012.15126.[19] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, andP. Jarillo-Herrero, Tunable phase boundaries and ultra-strong coupling superconductivity in mirror symmetricmagic-angle trilayer graphene, arXiv:2012.01434.[20] S. Carr, S. Fang, P. Jarillo-Herrero, and E. Kaxiras, Pres-sure dependence of the magic twist angle in graphenesuperlattices, Phys. Rev. B , 085144 (2018).[21] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe,T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov,and D. K. Efetov, Untying the insulating and supercon-ducting orders in magic-angle graphene, Nature , 375(2020).[22] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F.Young, Independent superconductors and correlated in-sulators in twisted bilayer graphene, Nature Physics ,926 (2020).[23] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek,and J. I. A. Li, Tuning electron correlation in magic-angle twisted bilayer graphene using coulomb screening, arXiv:2003.11072.[24] H. Guo, X. Zhu, S. Feng, and R. T. Scalettar, Pairingsymmetry of interacting fermions on a twisted bilayergraphene superlattice, Phys. Rev. B , 235453 (2018).[25] J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q.Sun, and C. Wang, Phases of a phenomenological modelof twisted bilayer graphene, Phys. Rev. B , 075154(2018).[26] C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Chiralspin density wave and d + id superconductivity in themagic-angle-twisted bilayer graphene, Phys. Rev. Lett. , 217001 (2018).[27] H. Isobe, N. F. Q. Yuan, and L. Fu, Unconventionalsuperconductivity and density waves in twisted bilayergraphene, Phys. Rev. X , 041041 (2018).[28] M. Fidrysiak, M. Zegrodnik, and J. Spa lek, Unconven-tional topological superconductivity and phase diagramfor an effective two-orbital model as applied to twistedbilayer graphene, Phys. Rev. B , 085436 (2018).[29] D. M. Kennes, J. Lischner, and C. Karrasch, Strong cor-relations and d + id superconductivity in twisted bilayergraphene, Phys. Rev. B , 241407(R) (2018).[30] Y. Sherkunov and J. J. Betouras, Electronic phases intwisted bilayer graphene at magic angles as a result ofVan Hove singularities and interactions, Phys. Rev. B , 205151 (2018).[31] J. Gonz´alez and T. Stauber, Kohn-Luttinger supercon-ductivity in twisted bilayer graphene, Phys. Rev. Lett. , 026801 (2019).[32] B. Roy and V. Juriˇci´c, Unconventional superconductivityin nearly flat bands in twisted bilayer graphene, Phys.Rev. B , 121407(R) (2019).[33] T. Huang, L. Zhang, and T. Ma, Antiferromagneticallyordered Mott insulator and d + id superconductivity intwisted bilayer graphene: a quantum Monte Carlo study,Science Bulletin , 310 (2019).[34] S. Ray, J. Jung, and T. Das, Wannier pairs in super-conducting twisted bilayer graphene and related systems,Phys. Rev. B , 134515 (2019).[35] Y.-Z. You and A. Vishwanath, Superconductivity fromvalley fluctuations and approximate SO(4) symmetry ina weak coupling theory of twisted bilayer graphene, npjQuantum Materials , 16 (2019).[36] F. Wu, A. H. MacDonald, and I. Martin, Theory ofphonon-mediated superconductivity in twisted bilayergraphene, Phys. Rev. Lett. , 257001 (2018).[37] Y. W. Choi and H. J. Choi, Strong electron-phonon cou-pling, electron-hole asymmetry, and nonadiabaticity inmagic-angle twisted bilayer graphene, Phys. Rev. B ,241412(R) (2018).[38] B. Lian, Z. Wang, and B. A. Bernevig, Twisted bilayergraphene: a phonon-driven superconductor, Phys. Rev.Lett. , 257002 (2019).[39] M. Xie and A. H. MacDonald, Nature of the correlatedinsulator states in twisted bilayer graphene, Phys. Rev.Lett. , 097601 (2020). [40] S. Liu, E. Khalaf, J. Y. Lee, and A. Vishwanath, Ne-matic topological semimetal and insulator in magic anglebilayer graphene at charge neutrality, arXiv:1905.07409.[41] J. Liu and X. Dai, Correlated insulating states and thequantum anomalous hall phenomena at all integer fillingsin twisted bilayer graphene, arXiv:1911.03760.[42] M. Xie and A. H. MacDonald, Weak-field hall resistiv-ity and spin/valley flavor symmetry breaking in MAtBG,arXiv:2010.07928.[43] See Supplemental Material at xxxx for detailed de-scriptions of critical magnetic field calculation, the self-consistent mean-field theory for finite-momentum pairingstate, the electron-acoustic phonon interaction, and addi-tional analysis of experimental data, which includes Refs.[44-47].[44] S. Carr, S. Fang, Z. Zhu, and E. Kaxiras, Exact con-tinuum model for low-energy electronic states of twistedbilayer graphene, Phys. Rev. Research , 013001 (2019).[45] K.-H. Bennemann and J. B. Ketterson, Superconductiv-ity: Volume 1: Conventional and Unconventional Super-conductors Volume 2: Novel Superconductors (Springer,Berlin, Heidelberg, 2008).[46] A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani,D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T.Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, andE. Zeldov, Mapping the twist-angle disorder and Landaulevels in magic-angle graphene, Nature , 47 (2020).[47] M. Koshino and Y.-W. Son, Moir´e phonons in twistedbilayer graphene, Phys. Rev. B , 075416 (2019).[48] A. MacDonald, H. Akera, and M. Norman, Quantum me-chanics and superconductivity in a magnetic field, Aus-tralian Journal of Physics , 333 (1993).[49] R. Bistritzer and A. H. MacDonald, Moir´e bands intwisted double-layer graphene, Proceedings of the Na-tional Academy of Sciences , 12233 (2011).[50] Z. A. H. Goodwin, F. Corsetti, A. A. Mostofi, and J.Lischner, Twist-angle sensitivity of electron correlationsin moir´e graphene bilayers, Phys. Rev. B , 121106(R)(2019).[51] Z. A. H. Goodwin, V. Vitale, F. Corsetti, D. K. Efetov,A. A. Mostofi, and J. Lischner, Critical role of device ge-ometry for the phase diagram of twisted bilayer graphene,Phys. Rev. B , 165110 (2020).[52] D. M. Basko and I. L. Aleiner, Interplay of Coulomb andelectron-phonon interactions in graphene, Phys. Rev. B , 041409(R) (2008).[53] G. D. Mahan, Many-Particle Physics (Kluwer Aca-demic/Plenum Publisher, New York, 2000).[54] P. B. Allen and R. C. Dynes, Transition temperature ofstrong-coupled superconductors reanalyzed, Phys. Rev.B , 905 (1975).[55] P. Morel and P. W. Anderson, Calculation of the su-perconducting state parameters with retarded electron-phonon interaction, Phys. Rev. , 1263(1962). upplemental Material: Critical magnetic fields and electron-pairing in magic-angletwisted bilayer graphene Wei Qin, Bo Zou, and Allan H. MacDonald Department of Physics, The University of Texas at Austin, Austin, Texas 78712,USA ∗ (Dated: February 23, 2021) I. CONTINUUM MODEL OF MATBG
In the present study, the moir´e bands of magic angle twisted bilayer graphene (MATBG) are calculated within acontinuum-model Hamiltonian [S1, S2], where Dirac electrons at the two Brillouin zone corners of each graphene layercoupled by periodic sublattice-dependent local tunneling. The K-valley projected Hamiltonian H ,K = (cid:20) − iν F σ − θ/ · ∇ T ( r ) T † ( r ) − iν F σ θ/ · ∇ (cid:21) , (S1)where σ θ = e iθ/ σ z ( σ x , σ y ) e − iθ/ σ z , σ x,y,z are Pauli matrices acting on sublattice. The periodic interlayer localtunneling T ( r ) = P j =1 T j e i q j · r , where T j = t AA σ + t AB [cos( j π/ σ x + sin( j π/ σ y ] . (S2)where q = k θ (0 , q , = k θ ( ±√ / , − / k θ = (8 π/ a ) sin( θ/
2) with graphene lattice constant a = 2 .
46 ˚A. In the main text, we take t AB = 117 meV and t AA = ηt AB with η < H = X αα X k ∈ MBZ ψ † α ( k ) [ H αα ( k ) − µδ αα ] ψ α ( k ) , (S3)where ψ α ( k ) is an annihilation field operator and α = ( lτ sσn ) is a lumped label with l, τ, s, σ = ± n = ( n , n ) specifying the moir´e reciprocal lattice vector G n = n G + n G with n , being integers. As depicted in Fig. S1, G , = k θ ( √ / , ± /
2) are primitive moir´e reciprocal lattice vectors. Inthe momentum space, the moir´e Hamiltonian can be viewed as Dirac cones arranged on a honeycomb lattice andcoupled via nearest-neighbor tunneling matrices, as schematically depicted in Fig. S1. Therefore, the matrix elementof K-valley projected moir´e Hamiltonian[ H K ( k )] lσn,l σ n = ~ ν F e ilθ/ σ z σ · [ k + l ( G n − q )] e − ilθ/ σ z δ ll δ σσ δ nn + X j =1 T j δ l, − l δ σσ δ q j +2 q , G n + n , (S4)where v F ∼ m/s is the Fermi velocity of Dirac electron in isolated graphene. The two valleys in the continuummodel of MATBG are decoupled, indicating H ( k ) = diag[ H K ( k ) , H − K ( k )], where H − K ( k ) = T H K ( k ) with T denotingthe time-reversal operator [S1].Figure S1(b) and (c) illustrate a typical calculation of the density of states at different chemical potentials or bandfillings upon varying the interlayer tunneling ratio η = t AA /t AB for MATBG with θ = 1 . ◦ . When decrease η ,the flat-band width becomes larger and the band fillings of conduction and valence flat-band van Hove singularitiesapproach ν = ±
2, respectively. ∗ Electronic address: [email protected] a r X i v : . [ c ond - m a t . s up r- c on ] F e b T
The Ginzburg-Landau (GL) theory of superconductivity is based on the expansion of free energy of a system inpowers of superconducting order parameter [S4]. In the presence of magnetic field the free energy F s = F n + Z d r (cid:20) ~ m ∗ | ( ∇ − i e ~ c A ) ψ | + α ( T ) | ψ | + β ( T )2 | ψ | + B π (cid:21) , (S5)here ψ = √ n s e iφ is the complex order parameter, φ is the phase of order parameter, n s , m ∗ , and 2 e are the density,effective mass, and total charge of Cooper pair. On the right hand side of Eq. S5, the first term in the integral isresponsible for the kinetic energy of Cooper pair, the second term is the pairing potential energy, the third termresembles the two-body interaction energy, and the last term is magnetic energy. By varying the GL free energy withrespect to magnetic vector potential A and order parameter ψ ∗ , the supercurrent density j = 2 ~ en s m ∗ ( ∇ φ − e ~ c A ) , (S6)and the equation 12 m ∗ ( − i ~ ∇ − ec A ) ψ + α ( T ) ψ + β ( T ) | ψ | ψ = 0 . (S7)The GL coherence length or magnetic length ‘ c characterizes variation of the superconducting order parameter in thereal space, and defined as ‘ c = p ~ / m ∗ | α | . (S8)The upper critical field H c of type-II superconductor can be estimated from Eq. (S7). When the external appliedmagnetic field is close to H c , the superconducting order parameter ψ becomes small, and Eq. (S7) can be linearizedinto 12 m ∗ ( − i ~ ∇ − ec A ) ψ + αψ = 0 , (S9)which resembles the Schr¨odinger equation for a particle with mass m ∗ and charge 2 e subject to magnetic field B = ∇ × A . The solution of Eq. (S9) gives rise to Landau levels of the corresponding particle | α | = ~ ω c ( n + 12 ) , (S10)where n are positive integers and ω c = 2 eB/m ∗ c is the cyclotron frequency. H c is defined as the maximum magneticfield for the solution, namely, H c = m ∗ c | α | / ~ e = Φ s π‘ c , (S11)where Φ s = π ~ c/e ≈ . × − T · m is the magnetic (superconducting) quantum flux.We show that ‘ c and H c can be obtained from the free energy for finite-momentum pairing superconducting state.For pairing wavevector q , the real-space order parameter can be written as ψ = p n s ( q ) e i q · r . In the absence ofexternal magnetic field, the superconducting condensation energy is reduced to δF ( q ) = F s ( q ) − F n = ~ q m ∗ n s ( q ) + αn s ( q ) + 12 βn s ( q ) , (S12)and the relation of Eq. (S7) ~ q m ∗ + α + βn s ( q ) = 0 . (S13)Therefore, Eq. (S12) and Eq. (S13) result in δF ( q ) = − β ( | α | − ~ q m ∗ ) = δF (0)[1 − ( q /q c ) ] , (S14)where q c = p m ∗ | α | / ~ is the critical pairing wavevector defined by δF ( q c ) = 0. Based on Eqs. (S8) and (S11), wehave q c = 1 /‘ c , and H c = Φ s π‘ c = Φ s q c / π. (S15)The supercurrent density of Eq. (S6) is j = 2 ~ en s ( q ) q /m ∗ = 2 en s ( q ) v , (S16)where v = ~ q /m ∗ is the velocity of Cooper pair. For the present study, we can calculate the free energy as a functionof pairing momentum within the microscopic model for MATBG superconductor. Based on Eq. (S12), ∂F ( q ) ∂ ( ~ q ) = ~ q m ∗ n s ( q ) + (cid:20) ~ q m ∗ + α + βn s ( q ) (cid:21) ∂n s ( q ) ∂ ( ~ q ) = n s ( q ) v , (S17)and we have the supercurrent density j = 2 e ~ (cid:20) ∂F ( q ) ∂ q (cid:21) . (S18)In the above derivations, α and β are assumed to be independent of pairing wavevector q . Such a assumption inMATBG can be justified by Fig. 1(c) in the main text, where our numerical calculations of δF ( q ) can be well describedby Eq. (S14). III. INTERACTION HAMILTONIAN
We focus on electron pairing mainly mediated by in-plane optical phonons in graphene, which are thought tocontribute ∼
70% of the total electron-phonon coupling strength [S5]. Based on C v and time-reversal symmetry, thein-plane optical phonons of isolated graphene can be classified into doubly degenerate E modes at BZ center, and A , B modes at BZ corner [S6]. The optical phonon energies ~ ω E = 196 meV and ~ ω A = ~ ω B = 170 meV are muchhigher than the kinetic energy of flat-band electrons in MATBG, indicating the optical phonon-mediated interactionis essentially instantaneous. By further restricting to electron pairing between opposite valley and spin, the opticalphonon-mediated attraction H ep = − g X lτsσ Z d r ψ † lτsσ ( r ) ψ † l ¯ τ ¯ sσ ( r ) ψ l ¯ τ ¯ s ¯ σ ( r ) ψ lτs ¯ σ ( r ) − g X lτsσσ Z d r ψ † lτsσ ( r ) ψ † l ¯ τ ¯ sσ ( r ) ψ lτ ¯ s ¯ σ ( r ) ψ l ¯ τs ¯ σ ( r ) , (S19)where ¯ τ = − τ , ¯ s = − s , ¯ σ = − σ , g and g are estimated to be 52 and 69 meV · nm , giving the electron-electroninteracting strengths mediated by E and A ( B ) phonons, respectively [S2]. The real-space field operator can beexpanded by plane-wave basis via Fourier transformation ψ lτsσ ( r ) = 1 √ A X k ψ lτsσ ( k ) e i ( k + K lτ ) · r , (S20)where A is the area of the system, K lτ denotes the wavevector of the Dirac point labeled by l and valley τ . In themomentum space H ep = − X kk q X lτsσ h g ψ † lτsσ ( k ) ψ † l ¯ τ ¯ sσ ( q − k ) ψ l ¯ τ ¯ s ¯ σ ( q − k ) ψ lτs ¯ σ ( k )+ g X σ ψ † lτsσ ( k ) ψ † l ¯ τ ¯ sσ ( q − k ) ψ lτ ¯ s ¯ σ ( q − k ) ψ l ¯ τs ¯ σ ( k ) i . (S21)The Hamiltonian in the moir´e momentum representation can be obtained by performing the following substitution X k ψ lτsσ ( k ) → X n X k ∈ MBZ ψ lτsσn ( k ) (S22)where k on the right hand side is measured from the γ point of MBZ as depicted in Fig. S1, and ψ lτsσn ( k ) = ψ lτsσ ( k − lτ q + G n ) . (S23)In the moir´e momentum representation, Eq. (S21) is reduced to H ep = − X lττ sσσ X nn n n X qkk ∈ MBZ ( g δ ττ δ σσ + g δ τ ¯ τ ) δ n + n ,n + n ψ † lτsσn ( k ) ψ † l ¯ τ ¯ sσ n ( q − k ) ψ l ¯ τ ¯ s ¯ σ n ( q − k ) ψ lτ s ¯ σn ( k )= − X αα α α g αα α α X qkk ∈ MBZ ψ † α ( k ) ψ † α ( q − k ) ψ α ( q − k ) ψ α ( k ) , (S24)where the collective label α = ( lτ sσn ), δ n + n ,n + n is imposed by the momentum conservation, and the interactionmatrix element g αα α α = ( g δ ττ δ σσ + g δ τ ¯ τ ) δ ll δ l l δ ll δ τ ¯ τ δ τ ¯ τ δ ss δ s s δ s ¯ s δ σ ¯ σ δ σ ¯ σ δ n + n ,n + n . (S25)The model electron-electron interaction of Eq. (2) in the main text becomes H ee = 2 u X ll X τsσσ Z d r ψ † lτsσ ( r ) ψ † l ¯ τ ¯ sσ ( r ) ψ l ¯ τ ¯ sσ ( r ) ψ lτsσ ( r )= 2 u X ll τsσσ X nn n n X qkk ∈ MBZ ψ † lτsσn ( k ) ψ † l ¯ τ ¯ sσ n ( q − k ) ψ l ¯ τ ¯ sσ n ( q − k ) ψ lτsσn ( k )= 2 X αα α α u αα α α X qkk ∈ MBZ ψ † α ( k ) ψ † α ( q − k ) ψ α ( q − k ) ψ α ( k ) , (S26)where u is Coulomb repulsion strength, and the interacting matrix element u αα α α = ( uδ ττ ) δ ll δ l l δ τ ¯ τ δ τ ¯ τ δ ss δ s s δ s ¯ s δ σσ δ σ σ δ n + n ,n + n . (S27)Based on Eqs. (S24) and (S26), the total interacting Hamiltonian can be organized in a compact form as H in = H ep + H ee = 12 X αα α α X qkk ∈ MBZ [ V ( k , k )] αα α α ψ † α ( k ) ψ † α ( q − k ) ψ α ( q − k ) , ψ α ( k ) (S28)with [ V ( k , k )] αα α α = 4( u αα α α − g αα α α ) independent of momentum k and k . IV. BDG HAMILTONIAN WITH FINITE-MOMENTUM PAIRING
Within the mean-field theory, the total interacting Hamiltonian of Eq. (S28) can be decoupled into finite-momentumpairing channel as H in = 12 X αα X k ∈ MBZ (cid:8) ψ † α ( k ) [∆( q )] αα ψ † α ( q − k ) + H.c. (cid:9) + C q , (S29)where q is the pairing momentum with the pair potential defined as[∆ q ( k )] αα = X α α X k [ V ( k , k )] αα α α h ψ α ( q − k ) ψ α ( k ) i (S30) C q = − X αα α α X kk ∈ MBZ [ V ( k , k )] αα α α h ψ † α ( k ) ψ † α ( q − k ) ih ψ α ( q − k ) ψ α ( k ) i . (S31)Given the the anticommutation relation of fermionic operator, the matrix form of the superconducting gap satisfiesthe relation ∆ q ( k ) = − ∆ T q ( k ). In the Nambu spinor representation, the total Hamiltonian H = 12 X k ∈ MBZ Ψ † q ( k ) H BdG ( q , k )Ψ q ( k ) + 12 X k ∈ MBZ
Tr [ H ( q − k ) − µ ] + C q , (S32)where Ψ q ( k ) = (cid:2) ψ ( k ) , ψ † ( q − k ) (cid:3) T , and the BdG Hamiltonian H BdG ( q , k ) = (cid:20) H ( k ) − µ ∆ q ( k )∆ † q ( k ) −H T ( q − k ) − µ (cid:21) . (S33)The self-consistent gap equation can be formulated by using the Green’s function defined as [S7] G ( τ, q , k ) = −h T τ Ψ q ( τ, k )Ψ † q (0 , k ) i = (cid:20) G ( τ, k ) F q ( τ, k ) F † q ( − τ, k ) −G T ( − τ, q − k ) (cid:21) (S34)where τ denotes imaginary time, T τ is the time-order operator, G ( τ, k ) = −h T τ ψ ( τ, k ) ψ † (0 , k ) i and F q ( τ, k ) = −h T τ ψ ( τ, k ) ψ (0 , q − k ) i are the normal and Gorkov’s anomalous Green’s functions, respectively. By performingFourier transformation on τ , we have G ( τ, q , k ) = k B T X n e − iω n τ [ iω n − H BdG ( q , k )] − = k B T X n e − iω n τ U q ( k ) [ iω n − E q ( k )] − U † q ( k ) , (S35)wehre k B is the Boltzmann constant, T denotes temperature, ω n = (2 n + 1) πk B T with integer n , and U q ( k ) is theunitary matrix that diagonalizes H BdG ( q , k ), giving the quasiparticle energy spectrum E q ( k ). Therefore, we have G (0 , q , k ) = U q ( k ) k B T X n [ iω n − E q ( k )] − U † q ( k ) . (S36)Hereafter and in the main text, we use G ( q , k ) to denote G (0 , q , k ) and have G ( q , k ) = (cid:20) G ( k ) F q ( k ) F † q ( k ) −G T ( q − k ) (cid:21) = U q ( k ) f [ E q ( k )] U † q ( k ) , (S37)where f ( (cid:15) ) is the Fermi-Dirac distribution function. Obviously, the normal Green’s function also depends on thepairing momentum q , and the self-consistent gap equation is[∆ q ( k )] αα = X α α X k ∈ MBZ [ V ( k , k )] αα α α [ F q ( k )] α α , (S38)which is independent of k for the present model study. Hereafter, we simply use ∆ q to denote the pair potential. C q = − X αα α α X kk ∈ MBZ [ V ( k , k )] αα α α [ F † q ( k )] α α [ F q ( k )] α α = −
12 Tr F † q V F q , (S39)where V is interacting tensor matrix. To insist the particle number conservation, the chemical potential for finite-momentum pairing state should be determined self consistently by N = X k ∈ MBZ
Tr[ G ( k )] , (S40)with N the total number of the particles in the system.Although the above lumped label α = ( lτ sσn ) is convenient for formulating, the detailed expression for pairpotential ∆ q helps understanding the corresponding physical means. For example,[∆ q ] lτsσn,l ¯ τ ¯ sσ n = 4 X τ σ σ X n n X k [ uδ ττ δ σσ δ σ σ − ( g δ ττ δ σσ + g δ τ ¯ τ ) δ ll δ σ ¯ σ δ σ ¯ σ ] δ n + n ,n + n × h ψ l ¯ τ ¯ sσ ( q − k + l τ q + G n ) ψ lτ sσ ( k − lτ q + G n ) i . (S41)Due to the momentum conservation, we can use a single reciprocal lattice vector Q = G n + G n to label ∆ q insteadof two of ( G n , G n ), namely,[∆ q , Q ] lτsσ,l ¯ τ ¯ sσ = 4 X τ σ σ X n k [ uδ ττ δ σσ δ σ σ − ( g δ ττ δ σσ + g δ τ ¯ τ ) δ ll δ σ ¯ σ δ σ ¯ σ ] × h ψ l ¯ τ ¯ sσ ( q − k + l τ q + Q − G n ) ψ lτ sσ ( k − lτ q + G n ) i . (S42)Therefore the intralayer pair potential[∆ q , Q ] lτsσ,l ¯ τ ¯ sσ = 4 u X n k h ψ l ¯ τ ¯ sσ ( q − k + lτ q + Q − G n ) ψ lτsσ ( k − lτ q + G n ) i− g X n k h ψ l ¯ τ ¯ s ¯ σ ( q − k + lτ q + Q − G n ) ψ lτs ¯ σ ( k − lτ q + G n ) i δ σσ − g X n k h ψ lτ ¯ s ¯ σ ( q − k − lτ q + Q − G n ) ψ l ¯ τs ¯ σ ( k + lτ q + G n ) i . (S43)The interlayer pair potential[∆ q , Q ] lτsσ, ¯ l ¯ τ ¯ sσ = 4 u X n k h ψ ¯ l ¯ τ ¯ sσ ( q − k − lτ q + Q − G n ) ψ lτsσ ( k − lτ q + G n ) i . (S44)In the sublattice space, the matrix form of pairing potential∆ q , Q = (cid:18) [∆ q , Q ] AA [∆ q , Q ] AB [∆ q , Q ] BA [∆ q , Q ] BB (cid:19) . (S45)In the present study, numerical calculations are carried out by taking a cutoff on the reciprocal lattice vector Q since u (meV nm ) ( m e V ) = + 2 u (meV nm ) ( m e V ) = 0 u (meV nm ) ( m e V ) = - 2 u (meV nm ) -I m ( m e V ) = + 2 u (meV nm ) -I m ( m e V ) = - 2 u (meV nm ) -I m ( m e V ) = 0 (a) (b) (c)(d) (e) (f) FIG. S2: Zero-momentum pair potential vs. the strength of electron-electron repulsion u for MATBG with θ = 1 . ◦ and η = 0 .
85. The upper (lower) row show the intralayer (interlayer) intrasublattice pair potentials for different band fillings. the magnitude of pair potential decreases quickly upon increasing Q . In Fig. S1, we depict representative reciprocallattice vectors Q , , , . All of the reciprocal lattice vectors generated by P Q n =1 , , , are included in our calculations,where P are symmetry operations from point group C . V. PAIR POTENTIAL AND CRITICAL TEMPERATURE
Figure S2 (a)-(c) show that the intralayer pair potential ∆ is suppressed by increasing the electron-electronrepulsion u . As suggested by Eq. (S44), the interlayer pair potential ∆ can be labeled by vectors P ( Q − lτ q ).Figure. S2 (d)-(f) show nonmonotonous behaviors of ∆ upon increasing u , attributed to the combined effectsbetween optical phonon-mediated attraction and electron-electron repulsion. In addition to smaller magnitude, ∆ also exhibits a π/ .The mean-field critical temperature T c is obtained identify the temperature where the pair potential vanishes. InFigure S3 shows a typical results of T c as functions of u and chemical potential µ or band filling ν .As complementary results to Fig. 1(b) in the main text, Fig. S4 plot pair potentiala as functions of paring wavevectorfor different band fillings. VI. FREE ENERGY OF THE SUPERCONDUCTING STATE
The free energy of a superconducting system is defined microscopically as F s ( q , T ) = Ω[∆( q ) , µ q , q , T ] + N µ q , (S46) u (meV nm ) T c ( K ) = + 2 = 0 = - 2 -2 -1 0 1 2 (meV) T c ( K ) +2+1-1 = -2 (a) (b) FIG. S3: Mean-field results of the critical temperature T c for MATBG with θ = 1 . ◦ and η = 0 .
85. (a) T c vs. u for differentband fillings. (b) T c vs. chemical potential µ for u = 40 meV · nm . -0.2 -0.1 0 0.1 0.2 q x ( m e V ) = -1 -0.2 -0.1 0 0.1 0.2 q x ( m e V ) = -1.5 -0.2 -0.1 0 0.1 0.2 q x ( m e V ) = -2.5 -0.2 -0.1 0 0.1 0.2 q x ( m e V ) = -2 (a) (b)(d)(c) FIG. S4: (a)-(d) Intralayer intrasublattice pair potentials vs. pairing wavevector q x for MATBG with θ = 1 . ◦ , η = 0 .
85 andband fillings (a) ν = −
1, (b) ν = − .
5, (c) ν = −
2, (d) ν = − . meV
1. The MBZ is highlighted by dashed hexagon. where the thermodynamic grand potential Ω = − β ln Tr( e − βH ) , (S47)with β = 1 /k B T . For MATBG superconductor F s ( q , T ) = − β ln Tr h e − β P k ∈ MBZ Ψ † q ( k ) H BdG ( q , k )Ψ q ( k )+ P k ∈ MBZ
Tr[ H ( q − k ) − µ ( q )]+ C q i + N µ q = C q + N µ q + 12 X k ∈ MBZ
Tr [ H ( q − k ) − µ q ] − β X k ∈ MBZ Tr n ln h e − βE q ( k ) io = C q + An µ q − β X k ∈ MBZ Tr n ln h e − βE q ( k ) io = C q + An µ q + 12 β Tr X k ∈ MBZ ln f [ − E q ( k )] , (S48)where n is the carrier density measured from the charge neutrality of MATBG, A is area of the system. We haveused the relation Tr[ H ( q − k )] = 0 since H ( k ) consists of massless Dirac Hamiltonian. The free energy density f s ( q , T ) = F s ( q , T ) /A = C q /A + n µ q + 18 π β Tr Z MBZ d k ln f [ − E q ( k )] . (S49)Figure. S5 plot the superconducting condensation energy δF ( q ), defined in the main text, for a special case of vanishedelectron-electron repulsion. VII. CRITICAL SUPERCURRENT
Figure S6(a) shows a typical result of supercurrent density j x for MATBG calculated using Eq. (S18). Here j x is an odd function of pairing wavevector q x due to time-reversal symmetry. As illustrated in Fig. 1(d) of the maintext, δF ( q ) is nearly isotropic, suggesting the critical supercurrent density j c can be defined by the maxima of j x asdepicted in Fig. S6(a). Figure S6(b) plots j c and mean-field critical temperature T c as functions of chemical potential µ . Similar to T c , j c exhibits a dome-like feature.Table S1 summarizes experimental data on critical temperature T c and critical suppercurrent I c . The criticalsuppercurrent density can be extracted by j c = I c / √ A c , where the conducting area A c of MATBG sample is ∼ µ m -0.2 -0.1 0 0.1 0.2 q x -7-6-5-4-3-2-10 F ( - m e V ) -6-4-20246 j x ( A m - ) j c j c -2 -1 0 1 2 (meV) j c ( A m - ) T c ( K ) +2+1-1 = -2 (b)(a) FIG. S6: Supercurrent density for MATBG with θ = 1 . ◦ , η = 0 .
85, and u = 40 meV · nm . (a) Superconducting statecondensation energy (circles) and supercurrent density (squares) vs. pairing wavevector q x for ν = −
1. The horizontal dottedlines highlight the critical supercurrent density j c . (b) j c vs. chemical potential µ .The dashed curve gives the correspondingcritical temperature T c .TABLE S1: Experimental results for the superconducting transition temperature T c , and critical current I c for MATBG withdifferent twist angle θ , band filling ν , and screening environment. θ ( ◦ ) ν T c (K) I c (nA) Refs1.05 -1.82 1.2 50 [S8]1.16 -2.15 0.5 171.27 (1.33 Gpa) -2.33 3.1 120 [S9]1.27 (2.21 Gpa) -2.62 2 200-2.31 3.1 70 [S10]1.1 -1 0.14 50.67 0.16 71.48 0.65 161.15 -1.6 0.92 300 [S11]1.8 0.42 121.04 -2.43 1.3 20 [S12]1.08 -2.46 2.4 801.09 -2.79 2.5 201.18 -2.5 0.7 401.12 -2.47 4 160[S8–S12]. I c changes dramatically for MATBG with different twist angles and band fillings. For a given sample, e.g. θ = 1 . ◦ in Table S1, a higher value of T c is associated with a larger value of I c , being qualitatively consistent withour calculation shown in Fig. S6(b). The quantitative value of j c in our calculations is over an order of magnitudelarger than experimental results. There are two possible reasons that lead to suppressions on the experimentallymeasured j c for MATBG. Firstly, the strongest superconducting dome usually occurs in the band filling regime wheretwo flavors are partially occupied associated with reconstructed flat bands [S8–S10, S13]. Our model calculation doesnot include such electronic correlation-induced effect, resulting in an overestimated carrier density and supercurrentdensity. Secondly, the disorders in experimental system is inevitable. For MATBG superconductor, the presence oflocal charge disorder and twist angle disorder, as observed in a recent experimental work [S14], leads to the coexistenceof superconducting and non superconducting regions, which form Josephson links and limit the zero-magnetic-fieldcritical supercurrent. An experimental evidence for the Josephson links is the Fraunhofer-like quantum interference[S8–S12].1 -2 -1.8 -1.6 -1.4 -1.2 -1 v F * ( m / s ) -2 -1.5 -1 T c ( K ) H c ( T ) FIG. S7: Extracted v ∗ F from experimentally reported data of T c and perpendicular H c shown in the insert, where θ = 1 . ◦ and electronic correlation is weaken by metallic screening [S11]. VIII. ELECTRON-ACOUSTIC PHONON INTERACTION
The Fermi surface averaged electron-phonon spectrum function is defined as α F ( ω ) ≡ N µ X nm, kk | g nm ( k , k ) | δ ( (cid:15) n k − µ ) δ ( (cid:15) m k − µ ) δ ( ω − ω ph ) , (S50)where m and n are flat-band indices. Here N µ is the density of states at the Fermi level µ , g nm ( k , k ) is flat-band projected electron-acoustic phonon interacting matrix, (cid:15) n k denotes electron energy for flat band n and moir´emomentum k , and the phonon energy ω ph = v ph | k − k | with v ph = 2 × m/s the sound velocity for isolatedgraphene sheet. In the present study, g nm ( k , k ) is obtained by folding the in-plane acoustic phonon mode of isolatedgraphene into MBZ, where the interlayer potential is ignored. By explicitly considering the interlayer potential, thetwo in-plane acoustic phonon modes from the two layers of graphene recombine into layer symmetric and asymmetricmodes [S15]. We deliver the detailed derivation and expression of g nm ( k , k ) in a separate work. IX. EXPERIMENTAL DATA
Here we show a concrete example of extracting the average Fermi velocity from a recent experimental work onMATBG superconductor. As depicted in Fig. S7, v ∗ F does not show a negative correlation with T c , which is usuallyexpected for weak coupling superconductor with a given pairing glue. Therefore, we argue that the conventional weakpairing mechanism may be insufficient for explaining the exotic superconductivity observed in MATBG. X. PAIRING CHANNEL
In this section, we discuss the potential pairing channels for phonon-mediated superconductivity in MATBG. Asdiscussed in Sec. III, the in-plane optical phonon-mediated electron-electron attraction H op = − g X lτss σ Z d r ψ † lτsσ ( r ) ψ † l ¯ τs σ ( r ) ψ l ¯ τs ¯ σ ( r ) ψ lτs ¯ σ ( r ) − g X lτss σσ Z d r ψ † lτsσ ( r ) ψ † l ¯ τs σ ( r ) ψ lτs ¯ σ ( r ) ψ l ¯ τs ¯ σ ( r ) , (S51)2 TABLE S2: Pairing strengths for different intrasublattice pairing channels.
Pairing channel Valley triplet & Spin singlet Valley singlet & Spin triplet Equal spin h ψ + σ ↓ ψ − σ ↑ i = −h ψ + σ ↑ ψ − σ ↓ i h ψ + σ ↓ ψ − σ ↑ i = h ψ + σ ↑ ψ − σ ↓ i h ψ + σ ↓ ψ − σ ↓ i Pairing strength g + g + g a − u g − g + g a − u g − g + g a − u For pairing symmetry analysis, we ignore the retarded effect of acoustic phonon-mediated electron-electron attraction,which can then be described by H ap = − g a X lτss σ Z d r ψ † lτsσ ( r ) ψ † l ¯ τs σ ( r ) ψ l ¯ τs σ ( r ) ψ lτsσ ( r ) , (S52)where g a represents the interaction strength. The electron-electron repulsion H ee = 2 u X ll X τsσσ Z d r ψ † lτsσ ( r ) ψ † l ¯ τs σ ( r ) ψ l ¯ τs σ ( r ) ψ lτsσ ( r ) . (S53)In contrast to Sec. III, here we release the restriction of electron pairing with opposite spins while still keep theassumption that electron pairing with opposite valleys. Within mean-field theory, the gap equation for real-spaceintrasublattice pair potential∆ + σs, − σs = g h ψ +¯ σs ψ − ¯ σs i − g h ψ +¯ σs ψ − ¯ σs i + ( g a − u ) h ψ + σs ψ − σs i . (S54)The gap equation for intersublattice pair potential∆ + σs, − ¯ σs = − g h ψ + σs ψ − ¯ σs i + ( g a − u ) h ψ + σs ψ − σs i , (S55)indicating that the BZ center E optical phonon modes do not contribute to intersublattice pairing. Table S2 lists thepairing strengths for different intrasublattice pairing channels. Based on these results, we conclude that, for opticalphonon-mediated superconductivity, the spin-singlet intrasublattice pairing is the dominant electron pairing channel.In the main text, we restrict ourselves to consider this dominant electron pairing channel. [S1] R. Bistritzer and A. H. MacDonald, Moir´e bands in twisted double-layer graphene, Proceedings of the National Academyof Sciences , 12233 (2011).[S2] F. Wu, A. H. MacDonald, and I. Martin, Theory of phonon-mediated superconductivity in twisted bilayer graphene,Phys. Rev. Lett. , 257001 (2018).[S3] S. Carr, S. Fang, Z. Zhu, and E. Kaxiras, Exact continuum model for low-energy electronic states of twisted bilayergraphene, Phys. Rev. Research , 013001 (2019).[S4] K.-H. Bennemann and J. B. Ketterson, Superconductivity: Volume 1: Conventional and Unconventional SuperconductorsVolume 2: Novel Superconductors (Springer, Berlin, Heidelberg, 2008).[S5] Y. W. Choi and H. J. Choi, Strong electron-phonon coupling, electron-hole asymmetry, and nonadiabaticity in magic-angletwisted bilayer graphene, Phys. Rev. B , 241412(R) (2018).[S6] D. M. Basko and I. L. Aleiner, Interplay of Coulomb and electron-phonon interactions in graphene, Phys. Rev. B ,041409(R) (2008).s[S7] G. D. Mahan, Many-Particle Physics (Kluwer Academic/Plenum Publisher, New York, 2000).[S8] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional supercon-ductivity in magic-angle graphene superlattices, Nature , 43 (2018).[S9] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean,Tuning superconductivity in twisted bilayer graphene, Science , 1059 (2019).[S10] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A.Bachtold, A. H. MacDonald, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-anglebilayer graphene, Nature , 653 (2019).[S11] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, andD. K. Efetov, Untying the insulating and superconducting orders in magic-angle graphene, Nature , 375 (2020).[S12] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Independent superconductors and correlated insulators intwisted bilayer graphene, Nature Physics , 926 (2020). [S13] M. Xie and A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Phys. Rev. Lett. ,097601 (2020).[S14] A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani, D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T. Taniguchi, P.Moon, M. Koshino, P. Jarillo-Herrero, and E. Zeldov, Mapping the twist-angle disorder and Landau levels in magic-anglegraphene, Nature , 47 (2020).[S15] M. Koshino and Y.-W. Son, Moir´e phonons in twisted bilayer graphene, Phys. Rev. B100