Charge-4e superconductivity from multi-component nematic pairing: Application to twisted bilayer graphene
CCharge- e superconductivity from multi-component nematic pairing:Application to twisted bilayer graphene Rafael M. Fernandes and Liang Fu School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455, USA Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139 USA (Dated: January 21, 2021)We show that unconventional nematic superconductors with multi-component order parameterin lattices with three-fold and six-fold rotational symmetries support a charge-4 e vestigial super-conducting phase above T c . The charge-4 e state, which is a condensate of four-electron boundstates that preserve the rotational symmetry of the lattice, is nearly degenerate with a competingvestigial nematic state, which is non-superconducting and breaks the rotational symmetry. Thisrobust result is the consequence of a hidden discrete symmetry in the Ginzburg-Landau theory,which permutes quantities in the gauge sector and in the crystalline sector of the symmetry group.We argue that random strain generally favors the charge-4 e state over the nematic phase, as it actsas a random-mass to the former but as a random-field to the latter. Thus, we propose that two-dimensional inhomogeneous systems displaying nematic superconductivity, such as twisted bilayergraphene, provide a promising platform to realize the elusive charge-4 e superconducting phase. Introduction.
The collective behavior of interactingelectrons in quantum materials can give rise to a plethoraof exotic phenomena. An interesting example is charge-4 e superconductivity [1–7], an intriguing macroscopicquantum phenomena which was theoretically proposedbut is yet to be observed. In contrast to standard charge-2 e superconductors characterized by Cooper pairing, acharge-4 e superconductor is formed by the condensationof four-electron bound states. While a clear manifes-tation of this phase would be vortices with half quan-tum flux, hc e , many of its basic properties, such aswhether its quasi-particle excitation spectrum is gaplessor gapped, remain under debate [6].An interesting question is which systems are promisingcandidates to realize charge-4 e superconductivity. Onestrategy is to consider systems that display two conden-sates, and search for a stable state where pairs of Cooperpairs are formed even in the absence of phase coherenceamong the Cooper pairs. One widely explored option isthe so-called pair-density wave (PDW) state, in which theCooper pairs have a finite center-of-mass momentum [7].An unidirectional PDW is described by two complex gapfunctions ∆ ± Q that have incommensurate ordering vec-tors ± Q . Charge-4 e superconductivity, described by thecomposite order parameter ∆ Q ∆ − Q , is a secondary orderthat exists inside the PDW state. It has been proposedthat the PDW state can melt in two stages before reach-ing the normal state [1], giving rise to an intermediatestate in which there is no PDW order, (cid:104) ∆ ± Q (cid:105) = 0, butthere is charge-4 e superconducting order, (cid:104) ∆ Q ∆ − Q (cid:105) (cid:54) = 0.Such an intermediate phase is called a vestigial phase[8–10], as it breaks a subset of the symmetries bro-ken in the primary PDW state. The main drawbackof this interesting idea is the fact that the occurrenceof PDW states in actual materials seems to be ratherrare [7]. Even from a purely theoretical standpoint, chal-lenges remain in finding microscopic models that give a (a) (b) FIG. 1. A nematic superconducting state in a lattice withthree-fold or six-fold rotational symmetry (here, a honey-comb lattice is shown) is described by a two-component orderparameter (∆ , ∆ ) = ∆ (cos θ, sin θ ), represented here bybound states of electron pairs (red dots). The ellipses repre-sent, schematically, different orientations θ . Two competingvestigial phases are supported: (a) a Potts-nematic phase and(b) a charge-4 e phase. In (a), the angle θ associated with thenematic director is fixed, breaking the three-fold rotationalsymmetry. In (b), the three-fold rotational symmetry is pre-served and a coherent state of bound states of four electronsemerge. In both (a) and (b), (cid:104) ∆ i (cid:105) = 0, i.e. charge-2 e super-conducting order is absent. PDW ground state rather than a uniform superconduct-ing ground state. For these reasons, it is desirable tosearch for other systems that may host vestigial charge-4 e superconductivity.In this paper, we show that nematic superconductorsin hexagonal and trigonal lattices offer a promising alter-native. A nematic superconductor breaks both the gaugesymmetry associated with the phase of the gap functionand the three-fold/six-fold rotational symmetry of thelattice. Importantly, nematic superconductivity has beenexperimentally observed in doped Bi Se [11, 12] and in a r X i v : . [ c ond - m a t . s up r- c on ] J a n twisted bilayer graphene [13], two systems whose latticeshave three-fold rotational symmetry. Superconductingproperties that do not respect the three-fold lattice sym-metry were also observed in few-layer NbSe , although itis unclear whether this is a consequence of a nematic pair-ing state [14, 15]. Unless finite tuning is invoked [16, 17],nematic superconductivity is realized in systems wherethe order parameter transforms as a multi-dimensionalirreducible representation of the relevant point group G [13, 18–22]. Typical examples are two-dimensional repre-sentations (∆ , ∆ ) where ∆ and ∆ correspond to p x -wave/ p y -wave gaps or d x − y -wave/ d xy -wave gaps. Inter-estingly, it has been shown that a secondary compositeorder parameter Φ = (cid:16) | ∆ | − | ∆ | , − ∆ ∆ ∗ − ∆ ∗ ∆ (cid:17) ,corresponding to Potts-nematic order, can onset evenabove the superconducting transition temperature T c [20, 23].In this paper, we show that the very same mechanismthat favors a vestigial nematic phase also promotes a ves-tigial charge-4 e phase characterized by a non-zero com-posite order parameter ψ = ∆ + ∆ , but (cid:104) ∆ i (cid:105) = 0(see Fig. 1). In particular, we find that the effectiveGinzburg-Landau theory obtained after integrating outthe normal-state superconducting fluctuations has thesame form for both the nematic order parameter Φ andthe charge-4 e order parameter ψ . We show that this isa robust result stemming from the existence of a lineartransformation, called a perfect shuffle permutation , thatrelates Φ and ψ in the four-dimensional space spanned by∆ and ∆ . Such a transformation effectively permutesquantities in the “gauge sector” and in the “crystallinesector” of the group U(1) ⊗ G that defines the symmetryproperties of the system.This result implies that there are actually two com-peting vestigial phases that can onset before long-rangesuperconductivity sets in: nematic order, as studied pre-viously [20, 23], and charge-4 e superconductivity. Whilehigher-order terms in the superconducting free-energygenerally favor the former, we show that the couplingto random strain can fundamentally alter the balancebetween them. This is because random strain acts as arandom-field to Φ , but as a random-mass to ψ . Con-sequently, random strain, intrinsically present in actualmaterials, is expected to suppress Potts-nematic ordermuch more strongly than charge-4 e order. We thus con-clude that the most promising candidates to realize ves-tigial charge-4 e superconductivity are relatively inhomo-geneous nematic superconductors with strong supercon-ducting fluctuations, as expected for instance in quasi-2Dsystems. This analysis thus suggests that twisted bilayergraphene [24–32] offers a potentially viable platform torealize this elusive state of matter. Vestigial nematicity: the standard scenario.
We con-sider a nematic superconductor in a lattice with three-fold or six-fold rotational symmetry, described by a two- component order parameter (∆ , ∆ ). For concreteness,hereafter we will focus on the case where the point groupof the lattice is D , and ∆ ≡ (∆ , ∆ ) † transforms asthe E irreducible representation (irrep), correspondingto (cid:0) d x − y , d xy (cid:1) -wave gaps. Importantly, our results aregeneral and hold as long as ∆ transforms as one of thetwo-dimensional E -like irreps of the corresponding pointgroup D , D , C v , etc. The Ginzburg-Landau supercon-ducting action expanded to fourth order in ∆ is given by[16, 20, 23, 33]: S [ ∆ ] = (cid:90) q ∆ ∗ i,q χ − ij ( q ) ∆ j,q + u (cid:90) r (cid:16) | ∆ | + | ∆ | (cid:17) + γ (cid:90) r | ∆ ∆ ∗ − ∆ ∗ ∆ | (1)Here, χ − ij ( q ) is the inverse superconducting suscep-tibility in Fourier space, whereas u > γ areGinzburg-Landau parameters. Furthermore, q = ( q , ω n )and r = ( r , τ ), where q is the momentum, ω n is thebosonic Matsubara frequency, r is the position, and τ isthe imaginary time. Note that S has an enlarged contin-uous rotational symmetry ∆ ± i ∆ → e ± iθ (∆ ± i ∆ ).This emergent continuous rotational symmetry is re-duced to discrete ones only when higher-order terms areincluded, as we discuss later.The superconducting ground state depends on γ : if γ <
0, the action is minimized by ∆ = ∆ (1 , ± i ) † , corre-sponding to a time-reversal symmetry-breaking (TRSB)superconductor that preserves the six-fold rotationalsymmetry of the lattice. On the other hand, for γ > ∆ = ∆ (cos θ, sin θ ) † , witharbitrary θ . Such a pairing state is called nematic, asit preserves time-reversal symmetry but lowers the six-fold rotational symmetry of the lattice to two-fold. It isconvenient to construct the real-valued composite orderparameters ζ ≡ ∆ † σ y ∆ and Φ ≡ (cid:0) ∆ † σ z ∆ , − ∆ † σ x ∆ (cid:1) ,where σ µ is a Pauli matrix that acts on the two-dimensional subspace spanned by ∆ [10, 20, 23]. While ζ transforms as the A irrep of D , and is thus related toTRSB, Φ transforms as the E irrep, being related to six-fold rotational symmetry breaking. Clearly, if the groundstate is ∆ = ∆ (1 , ± i ) † , ζ (cid:54) = 0 but Φ = 0. On the otherhand, if ∆ = ∆ (cos θ, sin θ ) † , ζ = 0 while Φ (cid:54) = 0. Thesign of γ is ultimately determined by microscopic consid-erations. While weak-coupling calculations tend to favor γ < γ > (cid:104) Φ (cid:105) (cid:54) = 0, but super-conducting order is absent, (cid:104) ∆ (cid:105) = 0 (see Fig. 1(a)). Tosee this, we follow the procedure outlined in Ref. [10] andfirst note that the quartic terms in Eq. (1) can be rewrit-ten in terms of the TRSB bilinear ζ = ∆ † σ y ∆ and thetrivial bilinear λ ≡ ∆ † σ ∆ as S (4) = u (cid:82) r λ + γ (cid:82) r ζ .Here, σ is the identity matrix. Now, the Fierz iden-tity (cid:80) µ σ µij σ µkl = 2 δ il δ jk − σ ij σ kl implies a relationshipbetween the bilinears, ζ = λ − Φ . As a result, thequartic term can be rewritten as S (4) = u (cid:82) r λ − γ (cid:82) r Φ ,where u ≡ u + γ and, as defined above, Φ = (Φ , Φ ) = (cid:0) ∆ † σ z ∆ , − ∆ † σ x ∆ (cid:1) is the nematic bilinear. Since γ > S [ ∆ , λ, Φ ] = (cid:90) r Φ γ − (cid:90) r λ u + (cid:90) q ∆ ∗ i,q (cid:2) χ − ij ( q ) + λσ ij − Φ σ zij + Φ σ xij (cid:3) ∆ j,q (2)Note that Φ and λ have been promoted to indepen-dent auxiliary fields. Because the action is quadratic in∆ i , the superconducting fluctuations can be exactly in-tegrated out in the normal state, yielding an effectiveaction for Φ and λ . Since λ does not break any sym-metries, it is always non-zero and simply renormalizesthe static superconducting susceptibility. On the otherhand, Φ is only non-zero below an onset temperature.A large- N calculation [36], as performed in Ref. [23],indicates that (cid:104) Φ (cid:105) (cid:54) = 0 already above T c , implying thatvestigial nematic order precedes the onset of supercon-ductivity (see also the Supplementary Material, SM). In-terestingly, a vestigial nematic phase has been recentlyobserved in doped Bi Se [37, 38]. Competition between nematicity and charge- e super-conductivity. We now show that there is a hidden symme-try between the two-component real-valued nematic or-der parameter Φ and the complex bilinear ψ ≡ ∆ + ∆ .The latter breaks the U(1) gauge symmetry and is pre-cisely the charge-4 e order parameter (see Fig. 1(b)). Im-portantly, ψ (cid:54) = 0 ( ψ = 0) inside the nematic (TRSB)superconducting state.To see the unexpected connection between these twoorder parameters, we need to consider, besides the realbilinears discussed above, complex bilinears formed outof the primary order parameter ∆ , since the lattertransforms as the irrep Γ = e imθ ⊗ E of the groupU(1) ⊗ D . Writing the order parameter explicitly as afour-dimensional vector η ≡ (∆ (cid:48) , ∆ (cid:48)(cid:48) , ∆ (cid:48) , ∆ (cid:48)(cid:48) ) T , wherethe prime (double prime) denotes the real (imaginary)part, the bilinears are generally given by η T ( σ µ ⊗ σ m ) η .Here, the first Pauli matrix (with Greek superscript) inthe Kronecker product σ µ ⊗ σ m refers to the subspaceassociated with the two-dimensional irreducible represen-tation E of the point group D (dubbed the crystallinesector), whereas the second Pauli matrix (with Latin su-perscript) refers to the subspace associated with the U(1) group (dubbed the gauge sector). In this notation, thecomponents of the nematic bilinear become:Φ = η T (cid:0) σ z ⊗ σ (cid:1) η Φ = − η T (cid:0) σ x ⊗ σ (cid:1) η (3)The other real bilinears are given by ζ = η T ( σ y ⊗ σ y ) η and λ = η T (cid:0) σ ⊗ σ (cid:1) η . The charge-4 e bilinear ψ ≡ ψ (cid:48) + iψ (cid:48)(cid:48) , on the other hand, is: ψ (cid:48) = η T (cid:0) σ ⊗ σ z (cid:1) η ψ (cid:48)(cid:48) = η T (cid:0) σ ⊗ σ x (cid:1) η (4)The key point is that, although the Kronecker prod-uct ( M ⊗ N ) is non-commutative, in the case where M and N are square matrices it satisfies the property( M ⊗ N ) = ˜ P T ( N ⊗ M ) ˜ P , where ˜ P is the so-called per-fect shuffle permutation matrix [39]. Here, due to the mi-nus sign in the second equation of (3), a slightly modified2 × P is needed: P = − − (5)Physically, P permutes quantities from the crystallineand the gauge sectors of the four-dimensional spacespanned by η . Note that P is an orthogonal matrix, P − = P T = P . As a result, upon performing the uni-tary transformation ˜ η = P η , we see that while the bilin-ears ζ and λ remain invariant, (Φ , Φ ) → ( ψ (cid:48) , ψ (cid:48)(cid:48) ), i.e.the nematic bilinear is mapped onto the charge-4 e bilin-ear. Consequently, provided that the susceptibility in thequadratic term of Eq. (1) is invariant under the lineartransformation (S19), the effective action in the normalstate has the same functional form with respect to eitherΦ or | ψ | . This is the case if we consider the standardsusceptibility expression χ − ij ( q ) = (cid:0) r + q (cid:1) δ ij , where r ∝ T − T c, is a tuning parameter and T c, is the baresuperconducting transition temperature (see the SM).This is the main result of our paper: for the Ginzburg-Landau action in Eq. (1), which describes a nematicsuperconducting ground state in a lattice with three-foldor six-fold rotational symmetry, an instability towardsa vestigial nematic state at T nem implies an instabilitytowards a vestigial charge-4 e state at the same temper-ature T e = T nem . This degeneracy between nematicityand charge-4 e superconductivity is rooted on the invari-ance of the action upon a perfect shuffle that permuteselements from the crystalline and the gauge sectors. Selecting nematic or charge- e order. While the com-petition between vestigial charge-4 e and nematic ordersis robust, their degeneracy is lifted by additional terms nematic charge4 e T T δ ε FIG. 2. Schematic phase diagram of the vestigial nematic(transition temperature T nem , green) and vestigial charge-4 e ( T e , red) phases. Here, δ ε represents the strength of straininhomogeneity. Because random strain couples as a random-field to the nematic order parameter but as a random-mass tothe charge-4 e order parameter, the former is expected to besuppressed much more strongly than the latter. In the cleansystem, ∆ T ≡ T nem − T e is positive because of the sixth-orderterm in Eq. (6) that restricts the nematic director to threedirections (3-state Potts nematicity) and lifts the emergentdegeneracy between the two vestigial ordered phases. Notethat, as temperature is lowered, a superconducting transi-tion is expected (not shown here). Whether charge-4 e andnematic orders can coexist in the overlapping region of thephase diagram remains to be studied. in the action not considered in the analysis above. Forinstance, additional symmetry-allowed terms in the sus-ceptibility χ ( q ) can favor either the charge-4 e state, inthe case of a hexagonal lattice, or the nematic state, inthe case of a trigonal lattice (see SM). While here we fo-cus on classical phase transitions, where the dynamics of χ ( q ) is not important, the situation changes in the caseof quantum phase transitions, as the couplings betweenthe bosonic fields Φ and ψ and the electrons are expectedto generate different types of bosonic dynamics.More importantly, because the nematic order param-eter Φ is real and transforms as the E irrep of D ,there is a symmetry-allowed cubic term in the nematicaction proportional to (cid:0) Φ + Φ − (cid:1) , where Φ ± = Φ ± i Φ [18, 20, 23, 40, 41]. This term is related to a particularsixth-order term in the superconducting action (1) [33]: S (6) [ ∆ ] ∝ (cid:90) r (∆ + i ∆ ) (∆ ∗ + i ∆ ∗ ) + h . c . (6)In contrast, because ψ is complex and transforms asthe A irrep of D , such a cubic term is not allowed inthe charge-4 e action. This cubic term not only favors thenematic order over the charge-4 e order, but it also lowersthe symmetry of the nematic order parameter from U(1)to 3-state Potts [20, 23, 40, 42]. At first sight, this seemsto suggest that it would be challenging to find a vestigial charge-4 e instability occurring before the onset of vesti-gial nematic order. While it is possible that charge-4 e order could coexist with nematic order and onset at atemperature between T nem and T c (the renormalized su-perconducting transition temperature), this seems to bea rather contrived scenario. However, there is an impor-tant ingredient missing in the analysis: the coupling tolattice degrees of freedom. This is particularly impor-tant for nematic order, as it is known to trigger latticedistortions [41].We thus introduce the strain tensor ε ij = ( ∂ i u j + ∂ j u i ), with u denoting the lattice displace-ment vector. Decomposing it in the irreps of the D group, there are two relevant modes: the longitudinalmode, which transforms as A , ε A ≡ ε xx + ε yy + ε zz ,and the shear mode, which transforms as E , ε E ≡ ( ε xx − ε yy , − ε xy ). The leading-order couplings to thenematic and charge-4 e orders are given, respectively, bythe linear coupling ε E · Φ and by the quadratic coupling ε A | ψ | . While strain can be externally applied, it is in-trinsically present in materials as random strain causedby defects arising in the crystal growth or device fabri-cation. The key point is that random strain acts as arandom-field to the Potts-nematic order parameter, butas a random-mass (also called random- T c ) to the charge-4 e order parameter.This distinction is very important, as random-field dis-order is known to be much more detrimental to long-orderrange order than random-mass disorder. In the specificcase of the 3-state Potts model, random-field is believedto completely kill the Potts transition in two dimensions,and to suppress it in three dimensions [43–45]. Thus,one generally expects random strain to tilt the balancebetween the competing vestigial charge-4 e and nematicorders in favor of the former. The resulting schematicphase diagram is shown in Fig. 2.The condition T e > T nem is not enough to ensurea vestigial charge-4 e phase, as one needs to show alsothat the renormalized superconducting transition tem-perature T c inside the charge-4 e state is split from T e [10]. A large- N analysis indicates that, for sufficientlyanisotropic quasi-two-dimensional systems, T e and T c are indeed split [23, 36]. In this case, while the transi-tion at T e is XY-like, the transition at T c is Ising-likedue to the coupling ψ ∗ (cid:0) ∆ + ∆ (cid:1) between the charge-4 e and the superconducting order parameters [6]. Conclusions.
In this paper, we showed that a nematicsuperconductor in lattices with three-fold or six-foldrotational symmetry supports competing nematic andcharge-4 e vestigial orders. Such a competition is rootedon a perfect shuffle permutation that transforms one or-der parameter onto the other in the four-dimensionalspace spanned by the multi-component superconductingorder parameter. We showed that random strain providesthe most promising tuning knob to favor charge-4 e super-conductivity over nematic order, due to the fact that itacts as a random-field disorder to the latter, but as arandom-mass disorder to the latter. These results estab-lish a new class of systems – nematic superconductors –in which charge-4 e order may be realized.Nematic superconductivity has been now observed indoped Bi Se and in twisted bilayer graphene [11–13].Based on our results, the most favorable conditions forthe observation of charge-4 e superconductivity are sys-tems where superconducting fluctuations are strong (e.g.quasi-2D superconductors) and where random-strain ispresent (e.g. inhomogeneous superconductors). Twistedbilayer graphene seems to satisfy both conditions, giventhe ubiquitous twist angle inhomogeneity [46–49], and isthus a promising place to look for this elusive state ofmatter. Note that the mechanism proposed here, whichrelies on an exact discrete symmetry in Ginzburg-Landautheory for multi-component superconductors in general,is different from a recent proposal for charge-4 e super-conductivity based on an approximate SU(4) symmetryof twisted bilayer graphene [50].We thank A. Chubukov, P. Orth, J. Schmalian, and J.Venderbos for fruitful discussions. This work was sup-ported by the U. S. Department of Energy, Office ofScience, Basic Energy Sciences, Materials Sciences andEngineering Division, under Award No. DE-SC0020045(R.M.F.) and DE-SC0018945 (L.F.). [1] E. Berg, E. Fradkin, and S. A. Kivelson, Nature Phys. ,830 (2009).[2] L. Radzihovsky and A. Vishwanath, Phys. Rev. Lett. , 010404 (2009).[3] E. V. Herland, E. Babaev, and A. Sudbø, Phys. Rev. B , 134511 (2010).[4] D. F. Agterberg, M. Geracie, and H. Tsunetsugu, Phys.Rev. B , 014513 (2011).[5] E.-G. Moon, Phys. Rev. B , 245123 (2012).[6] Y.-F. Jiang, Z.-X. Li, S. A. Kivelson, and H. Yao, Phys.Rev. B , 241103 (2017).[7] D. F. Agterberg, J. S. Davis, S. D. Edkins, E. Fradkin,D. J. Van Harlingen, S. A. Kivelson, P. A. Lee, L. Radz-ihovsky, J. M. Tranquada, and Y. Wang, Annual Reviewof Condensed Matter Physics , 231 (2020).[8] L. Nie, G. Tarjus, and S. A. Kivelson, Proceedings of theNational Academy of Sciences , 7980 (2014).[9] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Rev.Mod. Phys. , 457 (2015).[10] R. M. Fernandes, P. P. Orth, and J. Schmalian, AnnualReview of Condensed Matter Physics , 133 (2019).[11] K. Matano, M. Kriener, K. Segawa, Y. Ando, and G.-q.Zheng, Nature Phys. , 852 (2016).[12] S. Yonezawa, K. Tajiri, S. Nakata, Y. Nagai, Z. Wang,K. Segawa, Y. Ando, and Y. Maeno, Nature Phys. ,123 (2017).[13] Y. Cao, D. Rodan-Legrain, J. M. Park, F. N. Yuan,K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, andP. Jarillo-Herrero, arXiv:2004.04148 (2020).[14] A. Hamill, B. Heischmidt, E. Sohn, D. Shaffer, K.-T. Tsai, X. Zhang, X. Xi, A. Suslov, H. Berger, L. Forr´o,et al., arXiv:2004.02999 (2020).[15] C.-w. Cho, J. Lyu, T. Han, C. Y. Ng, Y. Gao, G. Li,M. Huang, N. Wang, and R. Lortz, arXiv:2003.12467(2020).[16] D. V. Chichinadze, L. Classen, and A. V. Chubukov,Phys. Rev. B , 224513 (2020).[17] Y. Wang, J. Kang, and R. M. Fernandes,arXiv:2009.01237 (2020).[18] L. Fu, Phys. Rev. B , 100509 (2014).[19] J. W. F. Venderbos, V. Kozii, and L. Fu, Phys. Rev. B , 180504 (2016).[20] J. W. F. Venderbos and R. M. Fernandes, Phys. Rev. B , 245103 (2018).[21] Y. Su and S.-Z. Lin, Phys. Rev. B , 195101 (2018).[22] V. Kozii, H. Isobe, J. W. F. Venderbos, and L. Fu, Phys.Rev. B , 144507 (2019).[23] M. Hecker and J. Schmalian, npj Quantum Materials ,26 (2018).[24] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, Nature , 43(2018).[25] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken,J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe,T. Taniguchi, E. Kaxiras, et al., Nature , 80 (2018).[26] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watan-abe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean,Science , 1059 (2019).[27] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir,I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang,et al., Nature , 653 (2019).[28] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney,K. Watanabe, T. Taniguchi, M. A. Kastner, andD. Goldhaber-Gordon, Science , 605 (2019).[29] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian,M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi,J. Hone, C. Dean, et al., Nature , 95 (2019).[30] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule,J. Mao, and E. Y. Andrei, Nature , 91 (2019).[31] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora,R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, et al.,Nature Physics , 1174 (2019).[32] Y. Xie, B. Lian, B. J¨ack, X. Liu, C.-L. Chiu, K. Watan-abe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Na-ture , 101 (2019).[33] M. Sigrist and K. Ueda, Rev. Mod. Phys. , 239 (1991).[34] R. Nandkishore, L. Levitov, and A. Chubukov, NaturePhys. , 158 (2012).[35] R. M. Fernandes and A. J. Millis, Phys. Rev. Lett. ,127001 (2013).[36] R. M. Fernandes, A. V. Chubukov, J. Knolle, I. Eremin,and J. Schmalian, Phys. Rev. B , 024534 (2012).[37] Y. Sun, S. Kittaka, T. Sakakibara, K. Machida, J. Wang,J. Wen, X. Xing, Z. Shi, and T. Tamegai, Phys. Rev. Lett. , 027002 (2019).[38] C.-w. Cho, J. Shen, J. Lyu, O. Atanov, Q. Chen, S. H.Lee, Y. San Hor, D. J. Gawryluk, E. Pomjakushina,M. Bartkowiak, et al., Nature communications , 1(2020).[39] M. Davio, IEEE Transactions on Computers , 116(1981).[40] Y. Xu, X.-C. Wu, C.-M. Jian, and C. Xu, Phys. Rev. B , 205426 (2020).[41] R. M. Fernandes and J. W. F. Venderbos, Science Ad- vances (2020).[42] S. Jin, W. Zhang, X. Guo, X. Chen, X. Zhou, and X. Li,arXiv:1910.11880 (2019).[43] D. Blankschtein, Y. Shapir, and A. Aharony, Phys. Rev.B , 1263 (1984).[44] K. Eichhorn and K. Binder, Journal of Physics: Con-densed Matter , 5209 (1996).[45] M. Kumar, R. Kumar, M. Weigel, V. Banerjee, W. Janke,and S. Puri, Phys. Rev. E , 053307 (2018).[46] A. Uri, S. Grover, Y. Cao, J. Crosse, K. Bagani,D. Rodan-Legrain, Y. Myasoedov, K. Watanabe,T. Taniguchi, P. Moon, et al., Nature , 47 (2020).[47] J. H. Wilson, Y. Fu, S. Das Sarma, and J. H. Pixley,Phys. Rev. Research , 023325 (2020).[48] B. Padhi, A. Tiwari, T. Neupert, and S. Ryu, Phys. Rev.Research , 033458 (2020).[49] C. Tschirhart, M. Serlin, H. Polshyn, A. Shragai, Z. Xia,J. Zhu, Y. Zhang, K. Watanabe, T. Taniguchi, M. Huber,et al., arXiv:2006.08053 (2020).[50] E. Khalaf, P. Ledwith, and A. Vishwanath,arXiv:2012.05915 (2020).[51] J. W. F. Venderbos, V. Kozii, and L. Fu, Phys. Rev. B , 094522 (2016). Supplementary material for “Charge- e superconductivity from multi-componentnematic pairing: Application to twisted bilayer graphene” DERIVATION OF THE EFFECTIVE ACTIONWITHIN LARGE- N Here we derive the explicit form of the effectivenematic/charge-4 e action by performing a large- N calcu-lation, extending the general procedure outlined in Refs.[10, 23, 36]. For a system with D point group symme-try, the most general form of the quadratic part of thesuperconducting action is given by [33, 51]: S (2) [ ∆ ] = (cid:90) q (cid:16) r + q (cid:107) + υq z (cid:17) (cid:16) | ∆ ,q | + | ∆ ,q | (cid:17) + κ (cid:90) q (cid:104)(cid:0) q x − q y (cid:1) (cid:16) | ∆ ,q | − | ∆ ,q | (cid:17) +2 q x q y (cid:0) ∆ ,q ∆ ∗ ,q + ∆ ∗ ,q ∆ ,q (cid:1)(cid:3) (S1)where q (cid:107) ≡ ( q x , q y ) and υ > | κ | < η ≡ (∆ (cid:48) , ∆ (cid:48)(cid:48) , ∆ (cid:48) , ∆ (cid:48)(cid:48) ) T , it can berewritten as: S (2) [ η ] = (cid:90) q η Tq ˆ χ − ( q ) η − q (S2)with: ˆ χ − ( q ) = (cid:16) r + q (cid:107) + υq z (cid:17) ˆ I + κq (cid:107) (cid:16) cos 2 θ ˆ A − sin 2 θ ˆ A (cid:17) (S3)Here, θ ≡ arctan (cid:16) q y q x (cid:17) , the hat denotes a 4 × I ≡ σ ⊗ σ , ˆ A ≡ σ z ⊗ σ , ˆ A ≡ − σ x ⊗ σ . Explicitly:ˆ A = − − (S4)ˆ A = − − − − (S5)We now move on to the quartic terms. The first one isgiven by: S (4)1 [ ∆ ] = u (cid:90) r (cid:16) | ∆ | + | ∆ | (cid:17) = u (cid:90) r (cid:16) η T ˆ I η (cid:17) (S6)Performing a Hubbard-Stratonovich transformation,we introduce the auxiliary field λ and obtain: S (4)1 [ η , λ ] = − (cid:90) r λ u + (cid:90) r λ (cid:16) η T ˆ I η (cid:17) (S7)As for the second quartic term, S (4)2 [ ∆ ] = − γ (cid:90) r (cid:20)(cid:16) | ∆ | − | ∆ | (cid:17) + (∆ ∆ ∗ + ∆ ∗ ∆ ) (cid:105) (S8)there are two different ways to decouple it in terms ofauxiliary fields. In the first case, we introduce the ne-matic field: S (4)2 [ η , Φ ] = (cid:90) r Φ γ − (cid:90) r (cid:104) Φ (cid:16) η T ˆ A η (cid:17) + Φ (cid:16) η T ˆ A η (cid:17)(cid:105) (S9)An alternative way to decouple it is by using the iden-tity: (cid:20)(cid:16) | ∆ | − | ∆ | (cid:17) + (∆ ∆ ∗ + ∆ ∗ ∆ ) (cid:21) = (cid:2) ∆ + ∆ (cid:3) (cid:104) (∆ ∗ ) + (∆ ∗ ) (cid:105) (S10)We then introduce the charge-4 e auxiliary field ψ andobtain: S (4)2 [ η , ψ ] = (cid:90) r | ψ | γ − (cid:90) r (cid:104) ψ (cid:48) (cid:16) η T ˆ B η (cid:17) + ψ (cid:48)(cid:48) (cid:16) η T ˆ B η (cid:17)(cid:105) (S11)where we defined ˆ B ≡ σ ⊗ σ z and ˆ B ≡ σ ⊗ σ x , i.e.ˆ B = − − (S12)ˆ B = (S13)The action can thus be written as: S nem [ η , λ, Φ ] = (cid:90) q η Tq (cid:34) ˆ χ − ( q ) − (cid:88) i Φ i ˆ A i (cid:35) η − q + (cid:90) r Φ γ − (cid:90) r λ u (S14)or, equivalently, S e [ η , λ, ψ ] = (cid:90) q η Tq (cid:34) ˆ χ − ( q ) − (cid:88) i ψ i ˆ B i (cid:35) η − q + (cid:90) r | ψ | γ − (cid:90) r λ u (S15)Here, to simplify the notation, we introduced ψ = ψ (cid:48) and ψ = ψ (cid:48)(cid:48) . Moreover, we defined:ˆ χ − ( q ) = ˆ χ − ( q ) + λ ˆ I (S16) which corresponds to shifting the superconducting massterm to r = r + λ . Integrating out the superconductingfluctuations: (cid:90) D η exp (cid:26) − (cid:90) q η Tq (cid:104) ˆ χ − ( q ) − ˆ V (cid:105) η − q (cid:27) = N exp (cid:26) − (cid:90) q Tr ln (cid:104) ˆ χ − ( q ) − ˆ V (cid:105)(cid:27) gives the effective actions: S (eff)nem [ λ, Φ ] = − ∞ (cid:88) n =1 n (cid:90) q Tr (cid:104) ˆ χ ( q ) ˆ V nem (cid:105) n + (cid:90) r Φ γ − (cid:90) r λ u (S17) S (eff)4 e [ λ, ψ ] = − ∞ (cid:88) n =1 n (cid:90) q Tr (cid:104) ˆ χ ( q ) ˆ V e (cid:105) n + (cid:90) r | ψ | γ − (cid:90) r λ u (S18)where ˆ V nem = (cid:80) i Φ i ˆ A i and ˆ V e = (cid:80) i ψ i ˆ B i . Now,to relate ˆ A i and ˆ B i , we note the identity ( M ⊗ N ) =˜ P T ( N ⊗ M ) ˜ P for 2 × M, N , where ˜ P is theperfect shuffle permutation matrix [39]:˜ P = which is orthogonal. In our case, due to the extra minussign in Eq. (S5), we need a slightly modified orthogonalmatrix: P = − − (S19)which then gives ˆ B i = ˆ P T ˆ A i ˆ P as defined in Eqs. (S4),(S5), (S12), and (S13). Thus, we obtain:Tr (cid:104) ˆ χ ( q ) ˆ V nem (cid:105) n = Tr (cid:104) ˆ˜ χ ( q ) ˆ V e (cid:105) n (S20)upon exchanging ψ i ←→ Φ i . Here, ˆ˜ χ ≡ ˆ P T ˆ χ ˆ P . Thus, aslong as κ = 0 in Eq. (S3), we have ˆ˜ χ ≡ ˆ χ , implying thatthe two actions – nematic and charge-4 e – are identical.We now proceed to investigate the impact of κ (cid:54) = 0. Forsimplicity, we focus on the two-dimensional case, setting υ = 0 in Eq. (S3). We also consider classical finite-temperature phase transitions. Performing the tracesand integrals in Eqs. (S17)-(S18) and assuming uniformorder parameters, we obtain the free-energy densities (theGinzburg-Landau constants γ and u are rescaled by a fac-tor of temperature): F (eff)nem [ r, Φ ] = Φ γ − (cid:90) ∞ qdqπ (cid:0) q + r (cid:1) (cid:104) κ q − ( q + r ) (cid:105) − Φ (cid:90) ∞ qdqπ (cid:0) q + r (cid:1) (cid:104)(cid:0) q + r (cid:1) + 2 κ q (cid:105)(cid:104) κ q − ( q + r ) (cid:105) − ( r − r ) u (S21)and: F (eff)4 e [ r, ψ ] = | ψ | γ − (cid:90) ∞ qdqπ (cid:104)(cid:0) q + r (cid:1) + κ q (cid:105)(cid:104) κ q − ( q + r ) (cid:105) − | ψ | (cid:90) ∞ qdqπ (cid:104)(cid:0) q + r (cid:1) + 6 κ q (cid:0) q + r (cid:1) + κ q (cid:105)(cid:104) κ q − ( q + r ) (cid:105) − ( r − r ) u (S22)The mean-field transitions take place when thequadratic coefficients vanish. This gives the followingcritical values of r , r nem = γJ nem and r e = γJ e , with: J nem ≡ (cid:90) ∞ pdpπ (cid:0) p + 1 (cid:1) (cid:104) κ p − ( p + 1) (cid:105) (S23) J e ≡ (cid:90) ∞ pdpπ (cid:104)(cid:0) p + 1 (cid:1) + κ p (cid:105)(cid:104) κ p − ( p + 1) (cid:105) (S24)An explicit calculation gives: J e = 12 π (1 − κ ) (S25) J nem = J e − π (cid:20) − κ − κ ln (cid:18) κ − κ (cid:19)(cid:21) (S26)which implies that J e > J nem , i.e. r e > r nem . Now, r is generally a decreasing function of temperature, since r ∝ ξ − → ξ is the superconducting correlationlength). Consequently, the κ term in the action favorsthe charge-4 e instability over the nematic instability.Note that this analysis is valid for the case of triangularor hexagonal lattices. For trigonal lattices, there is anadditional allowed term in the susceptibility (S3) thatdepends on q z [51]. Such a term generates the cubicinvariant Φ + Φ − in the nematic free energy [23], withΦ ± = Φ ± Φ , which favors the nematic instability overthe charge-4 ee