Nematic superconductivity in twisted bilayer graphene
NNematic superconductivity in twisted bilayer graphene
Dmitry V. Chichinadze, Laura Classen, and Andrey V. Chubukov School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
Twisted bilayer graphene displays insulating and superconducting phases caused by exceptional flatteningof its lowest energy bands. Superconductivity with highest T c appears at hole and electron dopings, near half-filling for valence or conduction bands. In the hole-doped case, the data show that three-fold lattice rotationsymmetry is broken in the superconducting phase, i.e., a superconductor is also a nematic. We perform a com-prehensive analysis of superconductivity in twisted-bilayer graphene within an itinerant approach and present amechanism for nematic superconductivity. We take as an input the fact that at dopings, where superconductivityhas been observed, the Fermi energy lies in the vicinity of twist-induced Van Hove singularities in the densityof states. We argue that the low-energy physics can be properly described by patch models with six Van Hovepoints for electron doping and twelve Van Hove points for hole doping. We obtain pairing interactions for thepatch models in terms of parameters of the microscopic model for the flat bands, which contains both local andtwist-induced non-local interactions. We show that the latter give rise to attraction in di ff erent superconductingchannels. For electron doping, there is just one attractive d − wave channel, and we find chiral d ± id super-conducting order, which breaks time-reversal symmetry, but leaves the lattice rotation symmetry intact. Forhole-doping, we find two attractive channels, g and i -waves, with almost equal coupling constants. We showthat both order parameters are non-zero in the ground state, and explicitly demonstrate that in this co-existencestate the threefold lattice rotation symmetry is broken, i.e., a superconductor is also a nematic. We find twopossible nematic states, one is time-reversal symmetric, the other additionally breaks time-reversal symmetry.Our scenario for nematic superconductivity is based on generic symmetry considerations, and we expect it to beapplicable also to other systems with two (or more) attractive channels with similar couplings. I. INTRODUCTION
Twisted hexagonal heterostructures recently joined the family of condensed matter systems which display superconduc-tivity with yet unsettled pairing mechanism. Superconductivity (SC) has been observed in twisted bilayer graphene [1–4],twisted double-bilayer graphene [5–7], and trilayer graphene on boron nitride [8, 9] under various tuning conditions con-trolled by twist angle, pressure, filling, or external field. This high degree of tunability holds the promise to enlighten thepairing problem from many di ff erent angles and thereby improve our understanding of superconductivity in correlated elec-tron systems. Like in several other correlated systems, SC borders insulated phases, in which fermions either get localizedby Mott physics, or develop a competing order, which gaps excitations near the Fermi surface.The occurrence of SC and insulating phases is ascribed to an exceptional band flattening, which comes along with a verylarge hexagonal moir´e pattern in real space [10, 11]. The small bandwidth of the resulting isolated flat band increases therelative strength of both electron-electron [12–15] and electron-phonon [16, 17] interaction and promotes correlation e ff ects[18]. In twisted bilayer graphene (TBG), insulating and superconducting phases have been induced by adjusting the bandflatness either by changing the twist angle between the two graphene layers to the so-called magic angle [11] close to 1.1 ◦ ,or upon applying pressure in the vicinity of the magic angle [1–4]. The most prominent insulating states occur near halffilling of either conduction or valence bands at the density of two electrons or two holes per moir´e unit cell (near n = ± n = n = − T c ∼ K [4] wasdetected close to half-filling of the valence band ( n = − T c ≤ . K have been foundnear the half-filled conduction band and in between other integer fillings. The suppression of superconductivity by a smallmagnetic field points to spin-singlet pairing [1].The data indicate [19] that TBG with a twist angle near the magic one lies in a regime of moderate coupling, where theratio of Coulomb interaction (reduced by the large spatial scale of the moir´e pattern) and the width of the flat band is of orderone (both are in the range of 10 −
20 meV). Consequently, arguments have been made for both moderate coupling itinerantapproach and strong coupling Mott-type approach. Arguments for the strong coupling approach have been rationalized bythe fact that insulating states have been detected not only near n = ±
2, but also near other integer fillings n = ± ± T c of 3K is muchsmaller than the bandwidth, and that insulating behavior is rather fragile – it disappears already for small T ≥ K and forsmall fields H ≤ T , Ref. [2–4]. Within the itinerant approach, insulating states are viewed as competing states with sometype of order in the particle-hole channel, and superconductivity and competing orders are largely based on the notion thatclose to half-filling [19, 20] and, possibly, other integer fillings, the chemical potential nearly coincides with twist-induced a r X i v : . [ c ond - m a t . s up r- c on ] J un Van Hove (VH) singularities in the density of states, Ref. [21]. This generally amplifies the e ff ect of interactions both inparticle-hole and particle-particle channels [22].In our study, we analyze superconductivity within an itinerant approach. Earlier studies considered both phonon [23–26]and purely electronic pairing mechanisms [27–45]. The works on the electronic mechanism often explored a scenario,where the enhancement of the pairing in the doubly degenerate d − wave channel at densities near the VH singularities leadsto chiral d ± id superconductivity. A similar analysis had been previously performed for a single layer of graphene at highdoping [46, 47].The primary goal of our study is to take a new look at superconductivity in the presence of VH singularities near theFermi level in view of the recent experimental observations that a discrete C rotational symmetry is likely broken in thesuperconducting state of hole-doped TBG [48]. Strong nematic fluctuations were earlier reported in the STM measurementsin the normal state [19, 49]. The authors of [48] argued that for some hole dopings C breaking extends to the normal state,but for other dopings, where superconductivity has been observed, C breaking is only present below T c . The direction ofthe nematic order is di ff erent at dopings where nematicity is only seen in the superconducting state and where it extendsto the normal state. This likely indicates that the nematicity in the normal state is not the source of the nematic order atdopings where it emerges below T c . We take these results as a motivation for our study and analyze a possibility to obtain C breaking only in the superconducting state.The breaking of a lattice rotational symmetry is usually associated with nematic order, and a state in which this symmetryis broken below the superconducting T c is called nematic superconductor. Nematic superconductivity has been earlier dis-cussed for LiFeAs (Ref. [50]) and doped topological insulator Bi Se (Refs. [51–58]). We argue that the scenarios proposedfor these materials do not apply to TBG. The electronic scenario, discussed in earlier works, yields d ± id superconductingorder [32–35, 38, 39, 41, 44, 45]. This order breaks time-reversal symmetry, but preserves C rotational symmetry. Phonon-mediated pairing with s − wave superconductivity does not break C symmetry either. Here, we propose a novel scenario forpairing in TBG, which gives rise to nematic superconductivity. We argue that a combination of the geometry of VH pointsnear n = −
2, which double in number compared to electron doping, and the form of the e ff ective interaction, which wasargued to possess both on-site and nearest-neighbor components [13–15], gives rise to an attraction in two pairing channels.We note that, as soon as the two pairing channels are attractive, our analysis is based on symmetry and independent on thepairing mechanism leading to the attraction. In the classification of the lattice rotation group D , appropriate for TBG, oneof the attractive channels is the doubly degenerate E channel and the other is the single-component A channel. The d -wavestate, discussed in earlier works, belongs to E channel. Based on the number of nodes along the Fermi surface, the SCstates that we found here correspond to doubly degenerate ” g -wave” and ” i − wave”, respectively (each gap component inthe E channel changes sign eight times under a 2 π rotation around the center of the Brillouin zone, while the gap in the A channel changes sign twelve times under a full rotation). The presence of higher lattice harmonics in the superconductinggap structure of TBG was also discussed in [59]. We argue that the coupling constants in the two channels are quite close, soin a wide temperature range below T c the system is in the coexistence state, where both SC orders are present. Taken alone,each of the two states does not break C symmetry: the order parameter in the E channel E ± iE (the g -wave analogueof d ± id ) breaks U (1) phase and Z time-reversal symmetry, and the order parameter in the A channel breaks U (1) phasesymmetry. We show that in the coexistence state, C symmetry is broken, along with the overall U (1) phase symmetry. Thebreaking of C symmetry is due to the presence of special coupling terms in the Landau free energy, which are linear in the A order parameter and cubic in the E order parameter. We found two coexistence states: in one time-reversal symmetry isadditionally broken, in the other it is preserved. The time-reversal-symmetric state develops if the special coupling betweenthe two order parameters is su ffi ciently strong. Otherwise, time-reversal is broken in the coexistence state. We show thecorresponding phase diagrams in Fig. 1 as function of a tuning parameter α T , which determines the non-local interactionstrength and regulates how close the two di ff erent SC states are in energy.It is instructive to compare our findings with the earlier proposals for nematic superconductivity in TBG. For spin-singlet pairing, earlier works [27, 60] focused on the two-component E state, without an admixture of the A state. Inthis situation, nematic superconductivity can develop if the solution for the gap is non-chiral, ( ∆ E , ∆ E ) = ∆ E (cos γ, sin γ ),and the minima of the free energy are at three values of γ . We looked into this possibility, but found that for parametersextracted from the microscopic model [13, 14, 61] that we use, the solution for the E state is the chiral ∆ E ± i ∆ E , whichbreaks time-reversal, but preserves C lattice rotational symmetry. It was suggested [60] that fluctuation corrections canpotentially change the free energy and make the nematic configuration energetically favorable, if the nematic component ofdensity wave fluctuations is large in the normal state. In a similar spirit, it was argued in Ref. [37] that fluctuation-inducednematic superconductivity can develop in the vicinity of a transition into a nematic orbital ferromagnet. We did not analyzefluctuation corrections or preformed nematic phases in our model. Instead, we focus on the superconductivity comingalready from the bare interaction. In Ref. [62] nematic superconductivity in the triplet channel has been explored. Thiswork is likely applicable to twisted double-bilayer graphene, where data suggest spin-polarized pairing [6]. In TBG, which T � T C ×U(1)× Z U(1) U(1)× Z C is broken , T is broken Coexistence phase ( ) ,, + 1 2 C manifold { + 1 2 { Z T � T U(1) U(1)× Z C × U ( ) × Z C × U ( ) × Z C ×U(1) C i s b r ok e n T i s b r ok e n C i s b r ok e n T i s b r ok e n C i s b r ok e n T i s p r e s e r v e d Coexistence phase , , ( ) { C m a n i f o l d FIG. 1. A schematic phase diagram for the case when the pairing channels with E and A symmetry are nearly degenerate (i.e., theattractive interaction in the two channels has almost the same magnitude). We use the strength of the non-local interactions α T as aparameter to distinguish between the cases when the interaction in one of the channels is stronger than in the other. Immediately belowthe onset of the pairing, SC develops in one of the two channels, and the gap symmetry is either g ± ig ( E -channel, green shaded region)or i ( A -channel, blue shaded region), both states are rotationally symmetric. At a lower temperature, a coexistence state develops (redshaded region). Depending on the parameters in the Landau free energy, this state is either type I or type II. State of type I (left panel)evolves continuously between pure A and pure E states. This state breaks C rotational symmetry, i.e., it is a nematic superconductor,and also breaks time-reversal symmetry. For state of type II (right panel), there is an additional intermediate phase inside the coexistenceregion (purple region). In this intermediate phase, C rotational symmetry is broken, but time-reversal symmetry is preserved. In bothpanels, the transition between the E state and the co-existence state is first order (red dashed line), and the transition between the co-existence state and the A state is second order (solid red line). The three states in the ”bubble” in the middle of each panel are symmetrypartners of the C manifold; time-reversal symmetry is realized by exchanging E ↔ E . The directions of blue, green, and brown arrowscorrespond to the phases of A , E , and E order parameters, ∆ i = | ∆ i | e i φ i , counted from the y axis. For definiteness, we set the A orderparameter to be real. we consider, experiments point to spin-singlet pairing [1].The doubling of the number of Van Hove points as function of twist angle or pressure, and the di ff erence in the numberof VH points in valence and conduction bands has been considered for twisted bilayer graphene in Refs. [30, 63] andfor monolayer jacutingaite in Ref. [64]. In particular, Ref. [30] analyzed Kohn-Luttinger superconductivity within themodel with twelve VH points and Hubbard interactions. They found attraction in several channels, with the dominant onebeing spin-triplet and C symmetric. Our analysis di ff ers from Ref. [30] in two aspects. First, we argue that the non-localcomponent of the interaction gives rise to an attraction in more that one channel already at the bare level, and, second, weargue that, to find nematic superconductivity, one needs to move below the highest T c and analyze the coexistence phase.We also argue that the di ff erence in the number of VH points in twisted bilayer graphene (six for electron doping vs twelvefor hole doping) leads to di ff erent SC states, and that nematic superconductivity develops only for hole doping.The structure of our paper is as follows. In Sec. II we introduce the patch models with six and twelve VH points andextract the model parameters from the underlying microscopic tight-binding model with local and non-local interactions. InSec. III we solve the corresponding linearized gap equation and determine the symmetries of its solutions. We use the resultto derive the Landau free energy for spin-singlet superconductivity in TBG in Sec. IV. We analyze the free energy in detailand present possible SC phase diagrams in the end. We conclude in Sec. VI. Several details of the derivations are presentedin the Supplementary material.Before we proceed, we present a brief summary of our results. Summary of the results
We use as an input for our study the e ff ective tight-binding model for the flat bands introduced in Refs. [13, 14, 61, 65].We consider fermionic densities at which VH singularities are located near the chemical potential, and introduce patchmodels for fermions in hot regions near the VH points. We argue that the proper patch model near n ≈ n = − D , appropriate for TBG. This point group has three irreduciblerepresentations: two one-dimensional representations A , A and one two-dimensional representation E [66]. We find thatthe interaction in some pairing channels is attractive, because of the non-local component.For the six-patch model near n = E -channel ( d -wave) is attractive, while the A channel ( s -wave) isrepulsive. The A channel ( f -wave) does not contribute to spin-singlet pairing. This is in accordance with earlier studies ofthe six-patch model for TBG [28] and single-layer graphene [46]. The corresponding SC order parameter can be representedby a vector ∆ E p with number of components equal to the number of patches. Because the E channel is two-dimensional,there are two independent order-parameter vectors ∆ E p and ∆ E p , and the full SC gap is a linear combination ∆ SC p = ∆ E ∆ E p +∆ E ∆ E p . The type of the superconducting order depends on which linear combination minimizes the free energy. It isdetermined by the sign of the coupling term between ∆ E and ∆ E : β | ∆ E + ∆ E | . When β >
0, a chiral ∆ E ± i ∆ E statedevelops, when β < ∆ E , ∆ E ) = ∆ E (cos γ, sin γ ) . For the microscopic model that we use forTBG, we find β >
0, i.e. the SC state is chiral ∆ E ± i ∆ E . This order breaks time-reversal symmetry, but preserves latticerotational symmetry (the gap amplitude is the same at all six VH points).For the twelve patch model near n = −
2, we find that two channels, E and A , are attractive with nearly equal couplingconstants. Analogously to the six-patch case, the corresponding order parameters are twelve-component vectors, which canbe expressed as ∆ SC , E = ∆ E ∆ E p + ∆ E ∆ E p and ∆ SC , A = ∆ A ∆ A p . The minimization of the free energy for the E state takenalone again yields the chiral ∆ E ± i ∆ E state that breaks time-reversal symmetry, but preserves C lattice rotation symmetry.The A state with a single gap amplitude ∆ A preserves both time-reversal and lattice rotation symmetries. However, weshow that new states emerge at low temperatures, when both E and A gaps are non-zero. To study the order parameterin the coexistence state, we derive the Landau functional F [ ∆ E , ∆ E , ∆ A ]. The functional, taken to quartic order in ∆ i ,contains regular mixed terms, quadratic in ∆ E , ∆ E , and in ∆ A : γ ( | ∆ E | + | ∆ E | ) | ∆ A | + γ (cid:16) ( ∆ E + ∆ E ) ¯ ∆ A + c . c (cid:17) , andthe asymmetric term δ ¯ ∆ A (cid:104) ( ∆ E − i ∆ E ) ( ¯ ∆ E − i ¯ ∆ E ) + ( ∆ E + i ∆ E ) ( ¯ ∆ E + i ¯ ∆ E ) (cid:105) + c.c. The coe ffi cients γ , , β , and δ areall expressed via the parameters of the underlying microscopic model. The asymmetric term is special in the sense that itis linear in ∆ A and qubic in ∆ E , , yet it is invariant under all symmetry transformations from the D space group on thehexagonal lattice, as well as under time-reversal and U (1) gauge transformations. We argue that because of this term, theorder parameter in the coexistence state breaks C lattice rotational symmetry.To illustrate the root of the C breaking, we analyze separately the special case when the asymmetric term is absent,and the generic, proper case when it is present. In the special case, we found that there are two coexistence states. Bothare highly degenerate, with order parameter manifold U (1) × U (1) × Z in one phase, and U (1) × U (1) in the other. Thepresence of two U (1)’s implies that there is an additional continuous degeneracy besides the conventional U (1) overallphase degeneracy. The extra Z in one phase is associated with time-reversal. In the other phase, time-reversal operationis a part of the extra U (1) symmetry. In a generic case, when δ is finite, we find that the additional U (1) gets discretized.For small δ , we find that there exists a single coexistence phase with order parameter manifold U (1) × C × Z , where U (1)is phase degeneracy, C is a discrete symmetry with respect to lattice rotations, and Z is associated with time-reversal.The superconducting order breaks all three symmetries, including C symmetry of lattice rotations. This implies that thecoexistence state is a nematic superconductor. For larger δ we find that there appears a region within the coexistence state,where the order parameter manifold is U (1) × C . A SC order in this range is again nematic, but it does not break time-reversal symmetry. For the parameters of the microscopic model of Refs. [13, 14, 61], the value of δ is close to the boundarywhere the state with broken C and unbroken time-reversal symmetry develops. We therefore cannot rigorously argue for oragainst time-reversal breaking in the superconducting state of TBG. Still, we emphasize that for any δ >
0, the SC state inour twelve-patch model near n = − C lattice rotational symmetry, i.e., the SC state is also a nematic state. This isconsistent with the experiments, which near n = − L L a a a b b b A B
FIG. 2. The moir´e honeycomb superlattice and translation vectors relevant for the tight-binding Hamiltonian. Vectors a i correspond tonearest-neighbor hopping, b i correspond to hopping between fifth nearest neighbors (denoted by (cid:104) i j (cid:105) in Eq. (1)), and L i are primitivelattice vectors. Black and gray circles show A and B sublattice sites, respectively. II. EFFECTIVE PATCH MODELS FROM TIGHT-BINDINGA. Fermiology of twisted bilayer graphene
While a brute-force microscopic description of TBG is obstructed by the huge unit cells of the moir´e superlattice, low-energy continuum [11, 20, 67–74] and tight-binding [10, 75–78] models that couple both layers have been very successful inanalyzing the electronic properties of TBG – including the theoretical prediction of flat bands itself [11, 67, 70, 76]. Morerecent works derived e ff ective tight-binding models for the superlattice based on localized Wannier states exclusively for theisolated flat bands [12, 14, 61, 65]. These Wannier states have a three-peak structure centered around sites of the honeycomblattices, which is dual to the triangular moir´e lattice where the local charge density is concentrated. Such a structure givesrise to hopping between further neighbors. In this context, the ability to write down a tight-binding model exclusively forTBG flat bands has been discussed. To do so, one has to overcome Wannier obstructions [79, 80]. The obstruction occursif one implements symmetries at incommensurate twist angles which are not exact, but assumed to emerge. To avoid theobstruction, but still construct a tight-binding model for the flat-bands only, one can either consider commensurate structuresnear the magic angle with well-defined bands [65], or sacrifice one of the approximate symmetries [14]. This is why themodel we use below has a three-fold symmetry instead of a six-fold one.In this work, we employ the model for the dispersion proposed by Yuan and Fu [61] (see also [14]). We start from writingdown the tight-binding Hamiltonian for the moir´e superlattice in real space in terms of Wannier states H T B = H (0) T B + H (1) T B (1) H (0) T B = − (cid:88) i µ c † i c i + (cid:88) (cid:104) i j (cid:105) t (cid:104) c † i c j + h . c . (cid:105) + (cid:88) (cid:104) i j (cid:105) t (cid:104) c † i c j + h . c . (cid:105) , (2) H (1) T B = (cid:88) (cid:104) i j (cid:105) t (cid:20)(cid:16) c † i × c j (cid:17) z + h . c . (cid:21) . (3)Here, the sums go over the sites of the honeycomb lattice, which are centered on the AB or BA regions of the moir´e patternin TBG. The operators c i = (cid:16) c i , x , c i , y (cid:17) T annihilate electrons with p -wave-like orbital index x and y , µ denotes the chemicalpotential, t , t are real hopping amplitudes between nearest- and fifth-nearest-neighbors, and (cid:104) i j (cid:105) denotes fifth-nearestneighbor (see Fig. 2). A fifth-nearest neighbor is equivalent to a second-nearest neighbor within the same sublattice. Forsimplicity, we suppressed a spin index.The Hamiltonian H (1) T B possesses an orbital and spin U (1) × S U (2) symmetry, D space symmetry of the TBG lattice, andis symmetric under time reversal. It yields four spin-degenerate bands with dispersions E ±± = T d ± T sd ± (cid:112) | T sd | , (4) E n e r g y FIG. 3. Twisted bilayer graphene band structure calculated for t = , t = − . , t = .
08. The dashed lines show the positions of VanHove singularities for electron-doping ( µ (cid:39) . µ (cid:39) − . µ are di ff erent in the two cases. where T d = − µ + t (cos b k + cos b k + cos b k ) , (5) T sd = t exp( ik x ) + − i k x √ k y , (6) T sd = t (sin b k + sin b k + sin b k ) . (7)In Fig. 3 we show the calculated band structure for a particular choice of hopping magnitudes. This band structure is ingood agreement with previously published results [11, 61, 65, 71, 74]. In particular, it reproduces the splitting of the bandsalong the Γ M -line, obtained in first-principles calculations [2, 78, 81, 82] and e ff ective low-energy models [11, 71, 74]. Thebands are orbitally-polarized in terms of chiral orbitals c ± = ( c x ± ic y ) / √ µ down from the chargeneutrality point µ =
0, one first reaches the Lifshitz transition at which isolated hole pockets centered on the Γ M line appear(see Fig. 5). At such a transition, there is no VH singularity in the DOS. The reason is that VH points in this case are notsaddle points, but local maxima of the band spectrum. However, decreasing µ further, one reaches the value of µ at whichanother Lifshitz transition occurs. This time, the corresponding VH points are saddle points of the dispersion, and the DOSis logarithmically singular (this is what we earlier called a Van Hove singularity). The dashed lines in Fig. 3 mark thevalues of chemical potential at which saddle-type VH points are located on the Fermi surface. We will focus on these pointsbecause the large DOS increases the tendency towards superconductivity and competing orders. As we will be interestedonly in the states near the saddle-type VH points, we avoid a subtle issue whether in the presence of all symmetries ofTBG, the tight-binding model of Eq. (4), based on localized Wannier states exclusively for the isolated flat bands (Refs.[14, 61, 65]), is adequate everywhere in the Brillouin zone, or if there exist special k − points away from VH regions, whereone needs to invoke other bands to properly describe excitations [79].In Fig. 5 we show the evolution of the Fermi surface for either hole or electron doping (negative or positive µ ). We seedi ff erent behavior in the two cases. Upon electron doping, the system reaches a Lifshitz transition with six VH singularities,located away from the Brillouin zone boundary, but along high symmetry Γ M directions. Upon hole doping, the first Lifshitztransition creates additional pockets along the Γ M lines, but does not give rise to VH singularities. As µ decreases further,the system does undergo a Lifshitz transition accompanied by VH singularities. In this last case, there are twelve VH pointsin the Brillouin zone, each located away from the Brillouin zone boundary and also away from high symmetry directions.We show the Fermi surfaces with six and twelve Van Hove singularities separately in Fig. 4.For other values of hopping integrals we found three other scenarios: (a) Lifshitz transitions with six VH singularitiesfor both electron and hole doping, (b) six VH singularities at a Lifshitz transition for hole doping and twelve for electrondoping, and (ii) Lifshitz transitions with twelve VH singularities for both electron and hole doping. In our analysis belowwe focus on the Fermi-surface geometry in Figs. 3, 4 , and 5, because it appropriately describes the observed electron-holeasymmetry of the superconducting states in TBG [1–4]. We also note in passing that as the twelve VH singularities atthe Lifshitz transition upon hole doping form six sets of pairs with small separation within a pair, there is the intriguingpossibility [63] that for fine-tuned hopping parameters, VH points within each pair merge and create a set of six VHsingularities, each leading to a stronger (power-law) divergence of the DOS [83]. We, however, will not study this special � ��� �� K' K'
FIG. 4. Fermi surfaces for the two bands (shown by di ff erent colors) for µ = .
716 (left) and µ = − .
329 (right). Circles indicatethe positions of Van Hove singularities. For µ (cid:39) .
716 (electron doping) there are six Van Hove singularities, located along Γ M andsymmetry-related directions in the Brillouin zone. For µ (cid:39) − .
329 (hole doping) there are twelve VH singularities; neither is locatedalong a high-symmetry direction. case. � �� K' FIG. 5. Evolution of the Fermi surfaces with doping for the same hopping parameters as in Fig. 3. Top: evolution upon electron dopingfrom charge neutrality (Dirac) point on the left to Van Hove doping on the right. Bottom: same evolution upon hole doping.
B. An e ff ective low-energy patch model Below we analyze superconductivity near Lifshitz transitions accompanied by singularities in the DOS. For this we focuson states near the Van Hove points and introduce e ff ective patch models with either six or twelve patches. We first expandthe energies Eq. (4) around the VH points and approximate the hopping Hamiltonian by H = N p (cid:88) a = (cid:88) τ = ± (cid:88) σ = ↑ , ↓ (cid:15) a ( k ) f † a τσ ( k ) f a τσ ( k ) , (8)where f a τσ ( k ) annihilates an electron with momentum k in the vicinity of patch a in band τ with spin σ . The patch indexruns through a = . . . N p with N p = ff erent patch points (cid:15) a ( k ) = α a k x − β a k y with sgn( α a ) = sgn( β a ) are related by D and time-reversal symmetry,inherited from the microscopic Hamiltonian of Eq. (1). Thereby, half of the patch points belong to one of the two bands FIG. 6. Sketch of the interactions in six-patch and twelve-patch models. Blue and red dots mark patches around Van Hove points. Weshow only the interactions relevant for superconductivity. crossing the Fermi energy at VH doping, while the other half belongs to the other band. The patch points belonging to thesame band are related by a threefold rotation symmetry, while the patches from di ff erent microscopic bands are related byinversion. We show the location of patches in Figs. 4, 5, and 6.We next consider all symmetry-allowed coupling terms between fermions in the patches. In general, there are fourtypes of allowed interactions. These are intra-patch and inter-patch density-density and exchange interactions. Umklappprocesses are not allowed because VH singularities do not appear at momenta connected by a reciprocal lattice vector.A simple bookkeeping analysis shows that there are 6 (18) symmetry-allowed couplings for the six-patch (twelve-patch)model without orbital-mixing terms, and 9 (27), when these terms are included. The orbital mixing terms were found tobe very small numerically in the microscopic model [13, 14], and we do not include them. The most general interactingHamiltonian for the six-patch model without orbital mixing is [28] H Int p = (cid:88) a = (cid:88) τ = ± (cid:20) u f † a τ f a τ f † a τ f a τ + v f † a τ f a τ f † a ¯ τ f a ¯ τ + u f † a τ f a τ f † a + τ f a + τ + v f † a τ f a τ f † a + τ f a + τ + j f † a τ f a + τ f † a + τ f a τ + (cid:18) g f † a τ f a + τ f † a ¯ τ f a + τ + h.c. (cid:19)(cid:21) (9)We introduced ¯ τ = − τ and a labels half of the patches. We omitted spin indices for simplicity – the spin structure of eachterm is (cid:80) σ,σ (cid:48) f † σ f σ f † σ (cid:48) f σ (cid:48) .For the twelve-patch model, the most general interaction Hamiltonian is H int p = (cid:88) a = (cid:88) τ = ± (cid:20) u f † a τ f a τ f † a τ f a τ + v f † a τ f a τ f † a ¯ τ f a ¯ τ + u f † a τ f a τ f † a + τ f a + τ + v f † a τ f a τ f † a + τ f a + τ + u f † a τ f a τ f † a + τ f a + τ + v f † a τ f a τ f † a + τ f a + τ + j f † a τ f a + τ f † a + τ f a τ + g (cid:16) f † a τ f a + τ f † a ¯ τ f a + τ + h.c. (cid:17) + j f † a τ f a + τ f † a + τ f a τ + g (cid:16) f † a τ f a + τ f † a ¯ τ f a + τ + h.c. (cid:17) + u + f † a τ f a τ f † a + ( − a τ f a + ( − a τ + u − f † a τ f a τ f † a − ( − a τ f a − ( − a τ + v + f † a τ f a τ f † a + ( − a ¯ τ f a + ( − a ¯ τ + v − f † a τ f a τ f † a − ( − a ¯ τ f a − ( − a ¯ τ + j + f † a τ f a + ( − a τ f † a + ( − a τ f a τ + g + (cid:16) f † a τ f a + ( − a τ f † a ¯ τ f a + ( − a ¯ τ + h.c. (cid:17) + j − f † a τ f a − ( − a τ f † a − ( − a τ f a τ + g − (cid:16) f † a τ f a − ( − a τ f † a ¯ τ f a − ( − a ¯ τ + h.c. (cid:17)(cid:21) (10)with the patch index being defined modulo 6.Below, we will need the subset of interactions relevant to pairing, which are between fermions with opposite momenta. Inour model these fermions belong to di ff erent bands. The pairing interactions then only involve fermions with patch indices( a , τ ) and ( a , ¯ τ ) (see Fig. 6). This reduces the number of interaction terms relevant for SC to two, v , g , for the six-patchmodel and to five, v , g , g , g + , g − , for the twelve-patch model. We sketch the interactions which will be important for thepairing problem in Fig. 6.To estimate the values of the couplings, we need to compare Eqs. (9) and (10) with the corresponding interaction terms inthe microscopic model. A typical approximation for the four-fermion interaction term for a system with screened Coulomb QQTT TQQT
FIG. 7. Graphic representation of three types of interactions in the Hamiltonian of Eq. (11). QQ (orange) are density-density interactions, T Q and QT (blue) describe processes, in which an electron interacts with a local density while hopping to a neighboring site, and T T (green) describes pair hopping processes. interaction is to keep it local, i.e., approximate the interaction by the on-site Hubbard density-density interaction. The caseof TBG was argued to be di ff erent, because there is substantial overlap between Wannier states localized at neighboring sites[13, 14]. This peculiar property leads to a new form of the interaction Hamiltonian [13, 14], in which local density-densityinteractions and terms describing assisted nearest neighbor hopping are of the same order and have to be considered onequal grounds. We follow Kang and Vafek [13] and write the interaction Hamiltonian in real space as H int = V (cid:88) R (cid:88) o = x , y (cid:88) σ = ↑ , ↓ O o ,σ ( R ) , (11)where O o ,σ ( R ) = Q o ,σ ( R ) + α T T o ,σ ( R ) , (12) Q o ,σ ( R ) = (cid:88) n = (cid:16) c † o σ A ( R + a n ) c o σ A ( R + a n ) + c † o σ B ( R − a n ) c o σ B ( R − a n ) (cid:17) , (13) T o ,σ ( R ) = (cid:88) n = (cid:16) c † o σ B ( R − a n ) c o σ A ( R + a n ) − c † o σ A ( R + a n + ) c o σ B ( R − a n ) + h . c . (cid:17) . (14)Here, R runs over the centers of the honeycomb lattice corresponding to the triangular moir´e pattern, A , B denote thesublattice indexes, a n are three translation vectors of the honeycomb sites (see Fig. 2) and o is the Wannier orbital index,inherited from the original valley degrees of freedom. There are three types of interaction terms: QQ , T Q , QT , and T T terms. The QQ term describes local (Hubbard) density-density interactions within a single honeycomb, the T Q , QT termsdescribe the processes, in which an electron interacts with a local density while hopping to a neighboring site, and the T T term describes the pair-hopping processes, in which electrons interact with each other, while hopping to a neighboring site(see Fig. 7). In momentum space, H int becomes H int = (cid:88) k , q , k (cid:48) , q (cid:48) (cid:88) σ,σ (cid:48) , o , o (cid:48) δ ( k − q + k (cid:48) − q (cid:48) ) (cid:20) (cid:88) v , v (cid:48) F vv (cid:48) v (cid:48) v c † o σ vk c † o (cid:48) σ (cid:48) v (cid:48) k (cid:48) c o (cid:48) σ (cid:48) v (cid:48) q (cid:48) c o σ vq + α T (cid:88) v (cid:44) v (cid:48) (cid:18) F vvvv (cid:48) c † o σ vk c † o (cid:48) σ (cid:48) vk (cid:48) c o (cid:48) σ (cid:48) vq (cid:48) c o σ v (cid:48) q + F vvv (cid:48) v c † o σ vk c † o (cid:48) σ (cid:48) vk (cid:48) c o (cid:48) σ (cid:48) v (cid:48) q (cid:48) c o σ vq + F vv (cid:48) vv c † o σ vk c † o (cid:48) σ (cid:48) v (cid:48) k (cid:48) c o (cid:48) σ (cid:48) vq (cid:48) c o σ vq + F v (cid:48) vvv c † o σ v (cid:48) k c † o (cid:48) σ (cid:48) vk (cid:48) c o (cid:48) σ (cid:48) vq (cid:48) c o σ vq (cid:19) + α T (cid:88) v (cid:44) v (cid:48) (cid:16) F vv (cid:48) vv (cid:48) c † o σ vk c † o (cid:48) σ (cid:48) v (cid:48) k (cid:48) c o (cid:48) σ (cid:48) vq (cid:48) c o σ v (cid:48) q + F vvv (cid:48) v (cid:48) c † o σ vk c † o (cid:48) σ (cid:48) vk (cid:48) c o (cid:48) σ (cid:48) v (cid:48) q (cid:48) c o σ v (cid:48) q (cid:17)(cid:21) , (15)where v , v (cid:48) label the sublattice indexes A and B , F vv (cid:48) v (cid:48)(cid:48) v (cid:48)(cid:48)(cid:48) are coupling functions, which we present in Supplementary ma-terial, and α T measures the strength of non-local interactions (the QQ , T Q , QT , and T T terms are O (1) , O ( α T ) and O ( α T )0 v g i-a ii = ii iiii aa + FIG. 8. Diagrammatic representation of a system of coupled gap equations. Gray triangle is a fully renormalized superconducting vertex,red and blue lines are Green’s functions of electrons from the two bands. Summation over a is implied. terms, respectively). Kang and Vafek estimated α T to be around 1 /
4. We will use α T as a parameter, but keep it close to1 / H int to the band basis and project it onto the patches around VH points. We show the details in theSupplementary material and here present the results. For the interactions relevant to SC, we obtain, in units of V fromEq. (11), v = g = . + . α T (16)for the six-patch model and v = g = . + . α T g = . + . α T g + = . + . α T g − = . + . α T (17)for the twelve-patch model. Observe that the prefactors of the α T terms are large numbers.The interactions in Eqs. (16) and (17) are the bare ones. The true interactions relevant to superconductivity are thee ff ective, fully irreducible ones, which includes all corrections from particle-hole bubbles and also all renormalizations inthe particle-particle channel from fermions with energies above the characteristic scale, which is, roughly, the largest of T c and the Fermi energy at the VH points. The flow of the couplings upon integrating out fermions with higher energies is oftendescribed within renormalization group (RG) approach [84–87]. These renormalizations are particularly relevant when thebare interaction is repulsive in all pairing channels, as renormalizations may overcome the bare repulsion and make theinteraction attractive in one or more pairing channels, below certain energies. Physically, these renormalizations make thee ff ective interaction non-local, and the growing non-local component eventually gives rise to a sign change of the pairinginteraction in certain channels. In our case, the bare interaction is already non-local, and we show in the next section that itis already attractive in one channel for the six-patch model and in two channels for the twelve-patch model. In this situation,the RG-type renormalization of the bare interaction may a ff ect the magnitudes of the attractive couplings, but will unlikelychange qualitatively the results obtained with the bare interactions. We therefore proceed without including the RG flowof the couplings. We emphasize that here we only focus on the pairing channel and do not address the issue of competingorders. To study the interplay between superconductivity and competing orders, RG-type calculations are required.We also note in passing that previous works did apply RG to both six-patch models [28, 33, 34] and a twelve-patchmodel [30] for TBG. However, these works considered the cases when the RG flow of the couplings (or, at least, Kohn-Luttinger renormalizations from particle-hole bubbles) is necessary to overcome a bare repulsion and induce an attractivepairing interaction. III. GAP EQUATION
To study superconductivity, we introduce the gap function as ∆ στσ (cid:48) τ (cid:48) a = (cid:104) f a τσ f a τ (cid:48) σ (cid:48) (cid:105) , where, we remind, a labels thepatches, and τ and σ are band and spin indices. In the absence of orbital mixing, spin-singlet and spin-triplet channelsare degenerate in TBG because direct exchange between patches related by time inversion is absent [28]. A finite orbitalmixing slits spin-singlet and triplet channels. Depending on the sign of the orbital mixing term, either spin-singlet or spin-triplet SC will be favored [28, 32]. Experimentally, superconductivity in TBG is destroyed by small magnetic fields [1],which is consistent with spin-singlet pairing. We therefore will focus on spin-singlet pairing. We assume that the orbitalmixing term is smaller than the other interaction terms. In this situation, the Cooper pairs with zero total momentum arepredominantly made by fermions from di ff erent bands, see Figs. 6 and 8. Accordingly, wee set τ (cid:48) = ¯ τ in ∆ στσ (cid:48) τ (cid:48) a and expressit as ∆ στσ (cid:48) ¯ τ a = ∆ a i σ y ( ∆ a is the same for the two choices of τ ). The matrix gap equation then reduces to a set of three (six)1coupled equations for the six-patch (twelve-patch) model: ∆ p = − Π v g g g v g g g v ∆ p , (18) ∆ p = − Π v g − g g g g + g − v g + g g g g g + v g − g g g g g − v g + g g g g g + v g − g + g g g g − v ∆ p . (19)Here ∆ p = ( ∆ , ∆ , ∆ ) T , ∆ p = ( ∆ , ∆ , ∆ , ∆ , ∆ , ∆ ) T , and Π is the particle-particle polarization bubble (the same forall pairs of fermions). Diagonalizing the matrix gap equation, we obtain eigenvalues and eigenfunctions in di ff erent pair-ing channels. We classify the solutions of the gap equation according to the irreducible representations of the point group D = C × C , whose elements are rotations along the z-axis by ± π/ C ) and twofold rotations along the y- and symmetry-equivalent axes ( C ), Ref. [66]. The D group has two one-dimensional irreducible representations, called A and A ,and one two-dimensional representation, called E (the corresponding eigenvalue is doubly degenerate). Each representationcontains an infinite set of eigenfunctions, some describe spin-singlet and some spin-triplet order. The generic form of eigen-functions in A is cos(6 n θ k ) for spin-singlet pairing ( n = , , .. ) and sin((6 n + θ k ) for spin-triplet pairing with the polarangle θ k = arctan k y / k x counted from the k x axis. For A , the eigenfunctions are cos((6 n + θ k ) for spin-triplet and sin(6 n θ k )for spin-singlet pairing. For E , the eigenfunctions are (cos((6 n + θ k ) , sin((6 n + θ k )) and (cos((6 n + θ k ) , sin((6 n + θ k ))for spin-singlet pairing and (cos((6 n + θ k ) , sin((6 n + θ k )) and (cos((6 n + θ k ) , sin((6 n + θ k )) for spin-triplet pairing.The gap equation decouples between di ff erent representations, but not between di ff erent eigenfunctions within the same rep-resentation. In a generic case, when all Fermi surface points are relevant to pairing, all partial components get coupled in thegap equation. In patch models, the gap equation simplifies because only a limited number of harmonics is distinguishable.For simplicity, we will use the lowest harmonics to describe our solution of the gap equation.In our sign convention, a specific channel becomes attractive when the eigenvalue turns from negative to positive. Weshow below that if we keep only local QQ terms in the interaction (i.e., set α T = α T exceeds somecritical value, specific to a given channel.Solving the gap equation for the six-patch model, we find that only one eigenfunction from A and one from E contributeto spin-singlet pairing. To express the corresponding eigenfunctions, we note that the patches are centered along highsymmetry directions. In this situation, the polar angles of the patch locations θ i , i = , ,
3, are related by θ = θ − π/
3, and θ = θ − π/
3. We can then write the eigenfunctions as ∆ A p = (1 , , ∆ E p = (cos(2 θ ) , cos(2 θ − π/ , cos(2 θ + π/ ∆ E p = (sin(2 θ ) , sin(2 θ − π/ , sin(2 θ + π/ . (20)The eigenfunction in the A representation has the same sign in all patches and is analogous to an s − wave. The eigenfunc-tions in the E representation change sign four times as one makes the full circle along the Fermi surface, and in this respectare analogous to d − wave.The eigenvalues in the A and E channels are λ A p = − Π ( v + g ) (21) λ E p = − Π ( v − g ) . (22)Substituting the values of couplings from (16), we obtain λ A p = − Π V (1 . + . α T ) (23) λ E p = − Π V (0 . − . α T ) . (24)2We plot the eigenvalues as functions of α T in the left panel of Fig. 10. We see that the coupling in the A channel is negative(repulsive) for all values of α T , but the one in the E channel becomes attractive for α T > . θ m + = θ − m π/ θ = π − θ , and θ m + = θ − m π/ ∆ A p = (1 , , , , , ∆ A p = ( − , , − , , − , ∆ E + p = (cos(2 θ ) , cos(2 θ ) , cos(2 θ − π/ , cos(2 θ + π/ , cos(2 θ + π/ , cos(2 θ − π/ ∆ E + p = (sin(2 θ ) , − sin(2 θ ) , sin(2 θ − π/ , − sin(2 θ + π/ , sin(2 θ + π/ , − sin(2 θ − π/ ∆ E − p = (cos(4 θ ) , cos(4 θ ) , cos(4 θ + π/ , cos(4 θ − π/ , cos(4 θ − π/ , cos(4 θ + π/ ∆ E − p = (sin(4 θ ) , − sin(4 θ ) , sin(4 θ + π/ , − sin(4 θ − π/ , sin(4 θ − π/ , − sin(4 θ + π/ A eigenfunction corresponds to s − wave. The A eigenfunction is proportionalto sin(6 θ k ) at the patch points and corresponds to i -wave (12 nodes along the Fermi surface). Eigenfunctions E + and E − from E correspond to d − wave and g − wave, respectively (4 nodes and 8 nodes, see Fig. 9). Because the two states with E symmetry have di ff erent number of nodes, they decouple in the gap equation (this does not hold beyond the patch model).The eigenvalues in the four decoupled channels are λ A / p = − Π (cid:2) v + g ± ( g − + g + + g ) (cid:3) , (26) λ E ± p = − Π (cid:20) v − g ± (cid:113) g − + g + + g − g − g + − g − g − g + g (cid:21) . (27)Substituting the values for the couplings from (17), we obtain λ A / p = − Π V (cid:104) . + . α T ± (0 . + . α T ) (cid:105) , (28) λ E ± p = − Π V (cid:20) . − . α T ± (cid:113) . + . α T + . α T (cid:21) . (29)We plot the eigenvalues as functions of α T in the right panel of Fig. 10. We see that λ A p and λ E + p are repulsive for all valuesof α T , but λ A p and λ E − p become attractive for α T > .
19 and 0 .
21, respectively. Note that the values of α T needed forattraction are smaller than in the six-patch model and also smaller than the estimate for α T ∼ .
23 presented in Ref. [13].Furthermore, over some range of α T ∼ /
4, the couplings λ A p and λ E − p are almost identical, i.e., the critical temperatures T A c and T Ec are approximately the same.We emphasize that this observation represents a qualitative di ff erence to previous studies of superconductivity withinpatch models for Van Hove filling in TBG [30, 33, 88] because in our case no Kohn-Luttinger type corrections (or higherorder corrections associated with spin or charge fluctuations) are needed to induce attractive pairing interactions. As aconsequence, we anticipate a higher critical temperature than typically expected for Kohn-Luttinger type superconductivity.Fluctuation corrections will increase non-local interactions, i.e. shift α T to a larger value and somewhat increase T A c and T Ec from already non-zero values.In the next section we derive the Landau free energy and analyze the SC state below T c . We show that the near-degeneracybetween the eigenfunctions in the A and E − channels leads to a highly non-trivial phase diagram with a large region ofthe coexistence state, where superconductivity breaks not only the U (1) phase symmetry, but also three-fold C latticerotational symmetry, i.e., the SC state is a nematic superconductor. We show that two types of nematic superconductorsemerge, depending on system parameters. One additionally breaks time reversal symmetry, the other preserves it.Before we proceed, we make an adjustment of the E − state in the 12-patch model for further convenience. Namely, weuse the fact that any rotation of the two components of E − is still an eigenfunction and rotate them by 3 π/
4. The newcomponents E − = E and E − = E , are ∆ E p = (cid:0) cos(4 θ + π/ , . . . , cos(4 θ + π/ (cid:1) ∆ E p = (cid:0) sin(4 θ + π/ , . . . , sin(4 θ + π/ (cid:1) (30)3 ����������� �������� � � � � FIG. 9. The eigenfunctions for E and A states in the twelve-patch model. Numbering of Van Hove points is shown on top. Middle:components ∆ E p = cos(4 θ i + π/
4) (red) and ∆ E p = sin(4 θ i + π/
4) (green). Bottom: ∆ A p = sin(6 θ i ) (blue). Circles, squares, anddiamonds denote the values of gap function at di ff erent VH points. Lines are continuous functions of θ , obtained using symmetryreasoning. Viewed as continuous functions, ∆ E p and ∆ E p have eight nodes, and ∆ A p has twelve nodes. With this choice of basis eigenfunctions, ∆ E p and ∆ E p exchange as ∆ E , p → − ∆ E , p under the twofold rotations around k y and symmetry related axes. IV. LANDAU FREE ENERGY AND THE STRUCTURE OF THE SUPERCONDUCTING STATE
We express the gap function in the superconducting state as a linear combination of the eigenfunctions for the attractivepairing components. In the six-patch model we have ∆ SC p = ∆ E ∆ E p + ∆ E ∆ E p , (31)4 (a) (b) (c) FIG. 10. The eigenvalues for the six-patch model λ E and λ A (a) and for the twelve-patch model, λ E + , λ E − , and λ A ((b) and (c)), asfunctions of α T ( λ A is irrelevant and not shown). When an eigenvalue turns positive, the interaction in the corresponding pairing channelbecomes attractive. For the six patch model, λ E > α T > .
98. For the twelve patch model, λ E − , and λ A become positive when α T exceeds certain values. Panel (b) is for the interactions, extracted from the microscopic model (see text). For panel (c), we increasedexchange interactions by a factor of two. Observe that λ E − and λ A are nearly degenerate over some range of α T . where ∆ E and ∆ E are complex numbers. In the twelve-patch model we have ∆ SC p = ∆ A ∆ A p + ∆ E ∆ E p + ∆ E ∆ E p , (32)where ∆ E , ∆ E , and ∆ A are complex numbers.To analyze superconducting ground states, we derive the Landau free energy, F p = F p ( ∆ E , ∆ E ) and F p = F p ( ∆ E , ∆ E , ∆ A ), find their minima and obtain the magnitudes and phases of ∆ E and ∆ E for the six-patch modeland of ∆ A , ∆ E , and ∆ E for the twelve-patch model. The functional form of the Landau free energy for each model isdetermined by D and U (1) symmetries [89], however which superconducting state is realized depends on the parameters ofthe Landau free energy. We obtain these parameters by applying a Hubbard-Stratonovich decomposition to the underlyingfermionic model and integrating out fermions (see Supplementary material for details). A. Six-patch model
The Landau functional for the six-patch model to order ∆ E , has the form F p = α (cid:16) | ∆ E | + | ∆ E | (cid:17) + β (cid:16) | ∆ E | + | ∆ E | (cid:17) + β | ∆ E + ∆ E | . (33)As usual, near a superconducting instability, α ∝ ( T − T c ) and β >
0. The coupling β can be of any sign, as long as β + β >
0. Minimizing with respect to amplitudes and phases of ∆ E and ∆ E we find that for β > F p is minimized by( ∆ E , ∆ E ) = ∆ E e i φ (1 , ± i ) (Refs. [46, 47]). This state breaks U (1) phase symmetry and additionally breaks Z time-reversalsymmetry. For β < F p is minimized by ( ∆ E , ∆ E ) = ∆ E (cos γ, sin γ ), where γ is arbitrary. To fix γ , one needs toinclude terms of sixth order in ∆ E , . The relevant sixth-order term is [27, 33, 46, 57, 58, 89] F (6)6 p = λ (cid:104) ( ∆ E − i ∆ E ) ( ¯ ∆ E − i ¯ ∆ E ) + c.c (cid:105) . (34)For our six-patch model for electron-doped TBG, we derived β from the underlying microscopic model and found β >
0, i.e., the SC state is a nodeless chiral superconductor. Such a state, dubbed d ± id , has been found in several earlierstudies of superconductivity in TBG [32–35, 38, 39, 41, 44, 45]. It breaks time-reversal symmetry, but does not break C lattice rotational symmetry.It was argued that nematic fluctuations [27, 60] in the normal state (more accurately, nematic components of charge orspin density wave fluctuations) do a ff ect β , and if these fluctuations are strong, they can, in principle, reverse the sign of β and convert the SC state in the six-patch model into a nematic SC. We did not analyze the strength of nematic fluctuations inour six-patch model. Instead we show how a nematic SC state can still develop in the twelve-patch model for hole doping,even if β >
0, due to the presence of another superconducting component.5
B. Twelve-patch model
As we demonstrated in Sec. III, there are two attractive pairing channels for the twelve-patch model – one-component A and two-component E channels. Up to fourth order in the gap function, the Landau free energy is F p = α (cid:16) | ∆ E | + | ∆ E | (cid:17) + α | ∆ A | + β (cid:16) | ∆ E | + | ∆ E | (cid:17) + β | ∆ E + ∆ E | + β | ∆ A | + γ (cid:16) | ∆ E | + | ∆ E | (cid:17) | ∆ A | + γ (cid:104)(cid:16) ∆ E + ∆ E (cid:17) ¯ ∆ A + (cid:16) ¯ ∆ E + ¯ ∆ E (cid:17) ∆ A (cid:105) + δ (cid:104)(cid:16) ∆ E | ∆ E | + ¯ ∆ E ∆ E − ∆ E | ∆ E | + ∆ E | ∆ E | + ¯ ∆ E ∆ E − ∆ E | ∆ E | (cid:17) ¯ ∆ A + c.c. (cid:105) . (35)where bar on top of ∆ means complex conjugation, and α ∝ T − T Ec , α ∝ T − T A c change sign at the critical temperaturesfor the pairing in E and A channels. We find that all prefactors for the fourth-order terms – β , β , β , γ , γ and δ , arepositive.Immediately below the largest of T Ec and T A c , the system develops either E or A superconducting order. When T Ec islarger, we have for the order in the E channel F E p = α (cid:16) | ∆ E | + | ∆ E | (cid:17) + β (cid:16) | ∆ E | + | ∆ E | (cid:17) + β | ∆ E + ∆ E | . (36)This F E p has the same form as F p in the six-patch model. Like there, we found β >
0. Then the state immediately below T Ec is a nodeless SC, which breaks time-reversal symmetry, but does not break C lattice rotational symmetry.When T A c is larger, we have for the order in the A channel F A p = α ∆ A + β ∆ A . (37)The A order is odd under C rotations, but it does not break C lattice rotation symmetry.We now consider coexistence states, in which both A and E order parameters are non-zero. We see from Eq. (35) thatthere are two types of terms in F p , which contain products of ∆ A and ∆ E , . The terms with coe ffi cients γ and γ are”conventional” bi-quadratic terms, which in a generic case set relative magnitudes and phases of A and E gap components.However, there is the additional term in Eq. (35) with prefactor δ , which is linear in ∆ A and cubic in ∆ E . Such a termis allowed by all symmetries. Indeed, one can explicitly verify that it is symmetric with respect to an overall U (1) phaserotation and does not change under C and C rotations. For the invariance under C , it is essential that our choice ofeigenfunctions ∆ E and ∆ E transform under C as ∆ E ↔ − ∆ E (see Eq. (30) and discussion after it). The structure of this δ term is similar to that of the sixth-order term in Eq. (34) of the six-patch model. Indeed, the δ term can be re-expressed as − δ ∆ A (cid:32) (cid:104) ( ∆ E − i ∆ E ) ( ¯ ∆ E − i ¯ ∆ E ) + ( ∆ E + i ∆ E ) ( ¯ ∆ E + i ¯ ∆ E ) (cid:105) − i (cid:104) ( ∆ E − i ∆ E ) ( ¯ ∆ E − i ¯ ∆ E ) − ( ∆ E + i ∆ E ) ( ¯ ∆ E + i ¯ ∆ E ) (cid:105) (cid:33) + c.c. (38)We will show that the δ term in the twelve-patch model and the sixth-order term in the six-patch model will play a similarrole regarding the breaking of lattice rotation symmetry. We note in passing that the term cubic in one SC order parameterand linear in the other was recently proposed in Ref. [90] in the context of chiral p - and f -wave pairing states on the squarelattice, with application to Sr RuO .To understand the role played by the δ term, it is instructive to first consider the structure of the coexistence state withoutthis term, and then add it. This is what we do next.
1. The structure of the coexistence state for δ = Without loss of generality, we choose the phase of complex ∆ A to be zero, i.e., set ∆ A to be real. We parametrizecomplex ∆ E and ∆ E as (cid:32) ∆ E ∆ E (cid:33) = ∆ E e i φ + (cid:32) e i φ − cos γ e − i φ − sin γ (cid:33) , (39)6 T � T U(1)×U(1)× Z U(1)× U(1)
Coexistence phase
U(1) U(1)× Z ������ ( ) , , Z =1 2Z =1 2 � � � U(1) �� � �� � �� U(1) � � cos � � � sin � FIG. 11. A schematic phase diagram for the special case δ = α T as tuning parameter. The orange circle marks the point where A and E channels are exactly degenerate. In the greenshaded region the superconducting state is pure E ( E ± iE ), in the blue shaded region the SC is pure A . In the red shaded region, both E = ( E , E ) and A gap components are non-zero. The dashed red line marks a phase transition between two coexistence states withdi ff erent order parameter manifolds. On the left of this line the manifold is U (1) × U (1) × Z , on the right it is U (1) × U (1). In bothcases, the order parameter manifold contains an additional continuous U (1) symmetry. The insets illustrate symmetry operations. Thedirections of blue, green, and brown arrows correspond to the phases of A , E , and E order parameters. where γ ∈ [0 , π/ δ = F δ = p = α ∆ E + α ∆ A + β ∆ E + β ∆ A + γ ∆ E ∆ A + ˜ β (cid:16) cos γ + sin γ cos φ − (cid:17) + γ (cos 2 φ + cos 2 φ − − sin 2 φ + sin 2 φ − cos 2 γ ) , (40)where ˜ β = β ∆ E , ˜ γ = γ ∆ E ∆ A . Minimizing the functional, we find two types of solutions, one for ˜ γ < ˜ β , another for˜ γ > ˜ β . The first solution is realized when the coexistence state emerges out of the E state, the second is when it emergesfrom the A state.For ˜ γ < ˜ β we obtain from minimization cos 2 φ − = − ˜ γ ˜ β cos 2 φ + cos 2 γ = ˜ γ sin 2 φ + ˜ β sin 2 φ − . (41)At ∆ A =
0, ˜ γ =
0, and Eq. (41) yields γ = π/ φ − = ± π/
4, as expected for the pure E state.Substituting cos 2 φ − and cos 2 γ from Eq. (41) into the Landau free energy, we find that F δ = p = α ∆ E + α ∆ A + β ∆ E + β ∆ A + γ ∆ E ∆ A − ˜ γ ˜ β = α ∆ E + α ∆ A + β ∆ E + β ∆ A + γ ∆ E ∆ A − γ ∆ A β , (42)does not depend on φ + . This implies that the order parameter manifold contains, in addition to U (1) total phase symmetry,another, extra U (1), associated with the freedom to rotate the common phase of ∆ E and ∆ E with respect to ∆ A . In addition,Eq. (41) for fixed φ + allows two solutions ( φ − , γ ) and ( − φ − , π/ − γ ). One solution transforms into the other if we interchange ∆ E into ∆ E . The full order parameter manifold is then U (1) × U (1) × Z . One can verify that this Z is associated withtime-reversal symmetry.For ˜ γ > ˜ β , the solution Eq. (41) disappears. The new minima are at φ + = φ − = ± π/ φ + = ± π/ φ − = γ : F δ = p = α ∆ E + α ∆ A + β ∆ E + β ∆ A + γ ∆ E ∆ A + ˜ β − γ = α ∆ E + α ∆ A + β ∆ E + β ∆ A + γ ∆ E ∆ A + β ∆ E − γ ∆ E ∆ A (44)This means that the order parameter manifold again has an additional continuous U (1) symmetry. To obtain the full orderparameter manifold in this case, we note that the four solutions in Eq. (43) can be re-expressed as (cid:32) ∆ E ∆ E (cid:33) = i ∆ E (cid:32) cos γ sin γ (cid:33) , (45)if we allow γ to vary between zero and 2 π . This implies that the order parameter manifold for ˜ γ > ˜ β is U (1) × U (1). Thereis no additional Z , because the phase of ∆ E and ∆ E in (45) is either the same or di ff ers by π , in which case phase reversaldoes not create a distinct SC state. Put di ff erently, the interchange ∆ E ↔ ∆ E can be absorbed into a variation of γ .We show the phase diagram for δ = α T as a tuning parameter.Along the transition line at ˜ γ = ˜ β , one of the E components vanishes, and the order parameter manifold reduces to U (1) × Z .The existence of the continuous U (1) symmetry in the order parameter manifold is highly unusual. In general, one wouldexpect only one U (1) to be present, associated with the symmetry with respect to rotations of the common phase. We willsee below that in the presence of the δ term, the continuous U (1) symmetry is replaced by a discrete C symmetry.
2. The structure of the coexistence state for nonzero δ Next, we consider the full Landau free energy, Eq. (35), with the δ -term. We use the same parametrization as in Eq. (39).The full analysis of Eq. (35) is rather cumbersome, but the outcome can be understood by just expanding near the boundariesof the coexistence phase. Near the left boundary, where ∆ A (cid:28) ∆ E , we have at δ = γ = π/ φ − = ± π/
4. Accordingly,at finite δ , we set γ = π/ + (cid:15) γ and φ − = ± ( π/ + (cid:15) − ), where (cid:15) γ , (cid:15) − ∼ δ ∆ A / ∆ E (cid:28)
1. We will see below that this expansionis valid for δ < β γ . The solutions with opposite sign of φ − transform into each other under ∆ E ↔ ∆ E , i.e., the orderparameter manifold contains Z associated with time reversal, like for δ = (cid:15) γ , (cid:15) − , we obtain to leading order in ∆ A / ∆ E : (cid:15) γ = δ β ∆ A ∆ E sin φ + (cid:15) − = δ β ∆ A ∆ E cos φ + (46)Substituting these expressions back into the Landau free energy we obtain F p = F δ = p − δ β ∆ E ∆ A − δ β γ − δ β cos 3 φ + ∆ E ∆ A + δ β γ − δ β ∆ A . (47)We see that the free energy now depends on φ + via the cos 3 φ + term. For δ < β γ , minimization with respect to φ + yieldsthree solutions φ + = (0 , π/ , − π/ δ term reduces the additional continuous U (1) symmetry to a discrete C symmetry. The system spontaneously chooses one out of three allowed values of φ + , and thereby breaks lattice rotationalsymmetry and becomes a nematic superconductor. Note that one of the states has φ + = γ = π/
4. For thisstate, the magnitudes of ∆ E and ∆ E are equal, only the relative angle 2 φ − varies with ∆ A . However, the two E componentsof the gap are not equal in any given patch, as one gets multiplied by cos(4 θ i + π/ θ i + π/ θ i specify the directions towards VH points. For the other two solutions ( φ + = ± π/ E components of the gap are the same as for the first solution if we rotate θ i by ± π/
3. For δ > β γ , the φ + -dependentterm in the free energy Eq. (47) changes sign. In this case, another solution, with φ − approximately ± π/ ∆ A (cid:28) ∆ E ,becomes energetically favorable.We also note that (i) the prefactor for the term quadratic in ∆ A in Eq. (47) is negative, i.e., for non-zero δ the transitiontemperature into the coexistence state is larger than the original T A c , where α in Eq. (35) changes sign and (ii) the freeenergy (47) has a term proportional to ∆ A . This term renders the transition between the pure E state and the coexistencestate first order.8We consider next the situation near the right boundary of the coexistence phase, where ∆ E (cid:28) ∆ A . Let us assume fordefiniteness that without the δ -term, φ − = π/ φ + = ∆ E = i cos γ, ∆ E − = − i sin γ ), cf. Eq. (43). When δ is non-zero,we expand φ − = π/ + ε − and φ + = + ε + . Minimizing with respect to ε ± , we obtain ε + = δ γ ∆ E ∆ A sin γ − cos γ sin γ cos γε − = δ γ ∆ E ∆ A (cos 3 γ − sin 3 γ ) , (48)and at the minimum F p = α ∆ E + α ∆ A + (cid:32) β + β − δ γ (cid:33) ∆ E + β ∆ A + ( γ − γ ) ∆ E ∆ A + β δ (sin(6 γ ) − ∆ E γ ∆ A . (49)Contrary to the previous case, there is no U (1) breaking term at order O ( δ ). However, such a term appears at order δ withthe structure δ sin 6 γ . Minimizing with respect to γ , we obtain γ = π/ + π n /
3, where n = , . . . , γ within a 2 π interval. One can verify that out of these six solutions, three are time-reversal partners of the other three, i.e., time-reversal symmetry is broken. One can understand this on physical grounds,because once the phase di ff erence 2 φ − between ∆ E and ∆ E becomes di ff erent from π , φ − and − φ − describe non-identicalgap configurations, hence under time-reversal the system transforms into a physically di ff erent state. The remaining threesolutions transform into each other under elements of C , i.e., the order parameter manifold is U (1) × C × Z , the same thatwe obtained near the left boundary of the coexistence phase. Note that in Eq. (49) the correction to α vanishes, and thereis no ∆ E term. As a consequence, the transition from the pure A state into the coexistence state is second order as long as4 (cid:16) β + β − δ γ (cid:17) β > ( γ − γ ) (see Refs. [91–93]).We verified that near the left boundary of the coexistence state, φ − increases with ∆ A (cf. Eq. (46)), and near the rightboundary φ − decreases as ∆ E increases (cf. Eq. (48)), i.e., ∆ E and ∆ E rotate towards each other. This strongly suggeststhat the gap structure in the coexistence state evolves continuously for small, but non-zero δ . We solved numerically for thegap at arbitrary ratio of ∆ A / ∆ E and found that this is indeed the case if δ < β γ . Specifically, for the ”symmetric” statewith φ + = γ = π/
4, we found a continuous change of φ − inside the coexistence phase from φ − ∼ π/ ∆ A (cid:28) ∆ E to φ − (cid:39) π/ ∆ A (cid:29) ∆ E . We show the phase diagram in Fig. 1 along with the structure of the pure and coexistence states.
3. The case of large δ We now show that a new state emerges at δ > β γ , which breaks C symmetry, but preserves time-reversal symmetry.To see this, we look again at the solutions close to the left and right boundaries. We found before that one of the solutionsfrom the C manifold is a symmetric one: φ + = γ = π/
4, i.e., | ∆ E | = | ∆ E | . Let us keep these values of φ + and γ , butnot assume that ∆ A / ∆ E is small and treat φ − as parameter. We will use this as an ansatz for the ground state for larger δ andthen verify that it is a stable minimum.Substituting into Eq. (35), we obtain F p = α ∆ E + α ∆ A + β ∆ E + β ∆ A + γ ∆ E ∆ A + cos(2 φ − ) (cid:16) γ ∆ E ∆ A + β ∆ E + √ δ ∆ E ∆ A cos( φ − ) (cid:17) . (50)One can check that at large enough δ , the free energy has smallest value when φ − = ± π . For this φ − , Eq. (50) reduces to F = α ∆ E + α ∆ A + (4 β + β ) ∆ E + β ∆ A + γ + γ ) ∆ E ∆ A − √ δ ∆ E ∆ A . (51)In such a state the phase of the two E components of the gap is opposite to the phase of the A component, i.e., all three gapcomponents, viewed as vectors, are directed along the same axis. Such a state preserves Z time-reversal symmetry.Rotational C symmetry requires that there must be two other states with the same energy. In total, we find φ + = γ = π/ φ − = ± πφ + = π/ γ = π/ φ − = π/ φ + = π/ γ = π/ φ − = − π/ . (52)9We now analyze where this ”collinear” state is located in the phase diagram. For this we assume that it is present forsome ∆ A and ∆ E and check its stability. For definiteness we choose the ”symmetric” state with φ + = , φ − = π , γ = π/ φ − = π + ε − , φ + = ε + and γ = π/ + ε γ . Substituting this into the free energy, we obtain to secondorder in ε i F p = α ∆ E + α ∆ A + (4 β + β ) ∆ E + β ∆ A + γ + γ ) ∆ E ∆ A − √ δ ∆ E ∆ A + (cid:104) √ δ ∆ E ∆ A − γ ∆ E ∆ A − β ∆ E (cid:105) ε − + √ δ ∆ A ∆ E ε γ + (cid:104) √ δ ∆ E ∆ A − γ ∆ E ∆ A (cid:105) ε + . (53)The stability conditions are then √ δ ∆ E ∆ A − γ ∆ E ∆ A ≥ √ δ ∆ E ∆ A − γ ∆ E ∆ A − β ∆ E ≥ , (54)These conditions set the boundaries of the collinear phase at √ δ + (cid:112) δ − β γ β ∆ A ≥ ∆ E ≥ √ δ − (cid:112) δ − β γ β ∆ A (55)and ∆ E ≥ γ √ δ ∆ A . (56)The first boundary is where fluctuations near φ − = π become unstable, the second is where fluctuations near φ + = γ do not give an additional constraint. Combining these two conditions, we obtain that the phasewith unbroken time-reversal symmetry exists once δ exceeds 2 β γ (we need √ δ + (cid:112) δ − β γ ] / β ≥ γ / √ δ ).It starts as a line in the phase diagram at δ = β γ and expands into the coexistence phase for larger δ . We show the phasediagram at large δ in Fig. 1 along with the states from the C manifold inside the collinear phase, and in Fig. 12 we presentthe plots of the total gap function ∆ SC p = ∆ A ∆ A p + ∆ E ∆ E p + ∆ E ∆ E p for the three regions within the coexistence phase inthe right panel of Fig. 1.The condition δ ≥ β γ coincides with the condition that φ − near the left boundary of the coexistence phase jumps from π/ π/
4. It then further increases with ∆ A and reaches π at the left boundary of the state with unbroken time-reversalsymmetry. The evolution of the gap between the right boundary of the coexistence state and the collinear phase is moreinvolved and we refrain from discussing it in detail. We note in passing that there is a certain analogy between the phasediagram and excitations in our case and for a 2D Heisenberg antiferromagnet in a magnetic field, whose phase diagram alsocontains an intermediate up-up-down phase with collinear ordering of spins in the three sublattices [94]. V. GAP STRUCTURE ALONG THE FULL FERMI SURFACE AND EXPERIMENTAL CONSEQUENCES
The gap structure along the full Fermi surface is shown in the three panels (a)-(c) in Fig. 12 for the three regions of thephase diagram in the right panel of Fig. 1 (the gap structure for the phase diagram in the left panel of Fig. 1 is the same, butwithout middle panel (b)). The gap function in panel (b) is for a SC state which preserves time-reversal symmetry, and hasnodes. The variation of this gap function with the angle θ along the Fermi surface is ∆ ( θ ) = ∆ A (cid:32) sin 6 θ + ∆ E ∆ A √ θ (cid:33) . (57)The number of nodes depends on the ratio ∆ E / ∆ A : for ∆ E / ∆ A (cid:46) / (2 √
2) there are twelve nodes, for ∆ E / ∆ A (cid:38) / (2 √ θ = , π/ , π, π/ θ and sin 4 θ components vanish). The location of the other nodal pointsdepends on the ratio ∆ E / ∆ A .0 (a) (b) (c) FIG. 12. The magnitude of the total gap function ∆ SC p = ∆ A ∆ A p + ∆ E ∆ E p + ∆ E ∆ E p along the Fermi surface, when ∆ A p , ∆ E p , and ∆ E p are viewed as functions of continuous θ rather than of θ i at Van Hove points. The three panels correspond to three coexistence statesin the right panel of Fig. 1. We have chosen the symmetric state with | ∆ E | = | ∆ E | (one of the states in C manifold). Panel (b) is forthe ”collinear” state in the middle of the right panel of Fig. 1, and panels (a) and (c) are for the states to the left and to the right of thecollinear state, respectively. We used ∆ A / ∆ E = . , φ − = .
86 in panel (a), ∆ A / ∆ E = .
71 in panel (b), and ∆ A / ∆ E = . , φ − = .
26 inpanel (c). The gap functions in panels (a) and (c) are complex numbers, and | ∆ SC p | has no nodes. The gap function in the collinear phaseis real and has nodes, because ∆ E p , ∆ E p and ∆ A p have nodes (see Fig. 9). The gap structure in the left panel of Fig. 1 is the same as inpanels (a) and (c). The gap functions in panels (a) and (c) are for the states with broken time-reversal symmetry. These gap functions arenodeless by obvious reasons. They can be parameterized by ∆ ( θ ) = ∆ A (cid:34) sin 6 θ − √ ∆ E ∆ A e i φ + (cid:0) cos φ − (cid:2) cos(4 θ − γ ) + sin(4 θ − γ ) (cid:3) + i sin φ − (cid:2) cos(4 θ + γ ) + sin(4 θ + γ ) (cid:3)(cid:1)(cid:35) , (58)where φ + , φ − and γ evolve as functions of the ratio ∆ E / ∆ A (see Eqs. (46) and (48)). The magnitude of the angle variationof ∆ ( θ ) depends on the ratio ∆ A / ∆ E , and is larger in panel (c) [i.e., for the state at a larger α T in the right panel of Fig. 1].The gap structures, presented in Fig. 12, can be probed experimentally, by, e.g., QPI analysis of STM data, and ARPESexperiments. Therefore, they are testable predictions of our theory. The states with and without nodes can also be distin-guished by other techniques, e.g., by measuring the flux penetration depth. The ratio of ∆ A and ∆ E likely can be variedby, e.g., changing the twist angle or adding uniform strain, which changes the degree of the non-locality of the interactionsand hence a ff ects our parameter α T in Fig. 1. A discrete C symmetry breaking was reported in Ref. [48] and motivatedour study. It can also be detected in STM studies and in angle-resolved photoemission spectroscopy with nanoscale reso-lution [95]. A time-reversal symmetry breaking can be detected via a broad range of probes [96], including measurementsof Kerr rotation [97] and zero-field muon-spin relaxation, which detects weak internal magnetic fields produced by spon-taneous currents, generated around impurities by time-reversal breaking superconducting order [98–101]. Domain walls insuch superconductors have magnetic signatures that could be detected in scanning SQUID and Scanning Hall probe mi-croscope measurements [102]. It was also proposed [103] that a nematic superconductor possesses topological skyrmions(bound states of two spatially separated half-quantum vortices), which can be detected by STM.On a qualitative level, a high T c / T F ratio, observed in magic-angle twisted bilayer graphene [1], is more consistentwith the existence of attractive pairing interactions at the bare level rather than with Kohn-Luttinger scenario, in which theattraction develops at second order in the interaction and is likely much weaker. VI. CONCLUSIONS
In this paper we performed a comprehensive analysis of superconductivity near Van Hove (VH) filling in twisted bilayergraphene (TBG) within an itinerant approach. The key motivation for our study has been the recent experimental finding [19,48] that the superconducting order in hole-doped TBG near n = − C lattice rotational symmetry, i.e., the SC stateis also nematic.We used as an input the e ff ective tight-binding Hamiltonian for the moir´e superlattice, which describes flat bands [14, 61].We argued that there are at least two VH fillings, one for hole doping, the other for electron doping. At VH filling for electrondoping, there are six VH points, located along high symmetry directions in the Brillouin zone, but away from the zoneboundary. At VH filling for hole doping, there are twelve VH points. They are symmetry related, but each is located awayfrom symmetry directions and the zone boundary. We derived e ff ective six-patch and twelve-patch models for fermions near1VH points and projected the interactions into the pairing channel. For the six-patch model, there are two symmetry-allowedpairing interactions in the spin-singlet channel. For the twelve-patch model this number is five. We obtained the values ofthe interactions by matching the patch models with the microscopic model of Kang and Vafek [13], which contains bothlocal (Hubbard) and non-local interactions. The relative strength of the non-local interaction is measured by the parameter α T , which was estimated to be around 0.23 in TBG. We argued that for this α T , the non-local interactions give rise toattraction in certain channels already at the ”bare” level, i.e., without including corrections to the pairing interaction fromthe particle-hole channel. In other words, the reconstruction of the band structure due to the twist and projection onto thenearly flat bands leads to attractive pairing interactions in hole-doped TBG. The attraction exists for both spin-singlet andspin-triplet channels. We concentrate on the first because experiments on TBG point to spin-singlet pairing [1].The symmetry of the superconducting order parameter can be classified based on the irreducible representations of thelattice rotation symmetry group D . They include two one-dimensional representations, A and A , and one two-dimensionalrepresentation, E . Each representation contains an infinite set of di ff erent eigenfunctions, but most become indistinguishablewithin patch models. For the six-patch model, we found that the relevant eigenfunctions are a constant ( s -wave) from A and d − wave-like (cos 2 θ i , sin 2 θ i ) from E , where θ i set the directions towards six VH points. We found that the interactionin the E channel is attractive and gives rise to d ± id SC order. It breaks time-reversal symmetry, but preserves C latticerotational symmetry. This agrees with earlier results for the six-patch model [28, 34] and with earlier studies of single-layergraphene around VH filling [46, 47].For the twelve-patch model, we found four di ff erent pairing channels: one in A , with a constant eigenfunction, one in A , with an eigenfunction changing signs between neighboring patches, and two in E with eigenfunctions (cos 2 θ i , sin 2 θ i )and (cos 4 θ i , sin 4 θ i ). We found that A and the 4 θ E channel are attractive, and that for realistic α T the coupling constantsin the two channels have near-equal magnitudes. We showed that pure A order breaks only U (1) phase symmetry, and pure E order is similar to that in the six-patch model, i.e., it breaks U (1) and Z time-reversal symmetry, but preserves C .Our key result is that in the coexistence state, where both E and A order parameters are non-zero, C symmetry is broken.We argued that this happens due to two reasons: (i) conventional biquadratic couplings between E and A order parametersdo not specify the coexistence state, and the order parameter manifold has an extra U (1) symmetry, in addition to the U (1)total phase symmetry, and (ii) the Landau free energy to quartic order contains a symmetry allowed term, which is linear in ∆ A and qubic in ∆ E , . This term breaks the extra U (1) symmetry down to threefold C . The system spontaneously choosesone of three equivalent states from C manifold and by doing this breaks C . As a result, the coexistence state turns out tobe a nematic superconductor. We found two phases with broken C . In one, time-reversal symmetry is also spontaneouslybroken. In the other, it is preserved.Our results present a scenario for the breaking of threefold lattice rotation symmetry in the superconducting state ofhole-doped TBG near n = −
2, where nematic superconductivity has been observed [19, 48]. We also consider it as ageneric, symmetry-based mechanism how a superconductor can break lattice rotational symmetry. We also emphasize thatalthough in our scenario the nematic long-range order emerges only in the coexistence superconducting phase, nematicorder generally survives in some range outside the coexistence phase, and nematic fluctuations are strong in the wholeregion where the pairing susceptibility is enhanced in both E and A channels. VII. ACKNOWLEDGMENTS
We thank E. Andrei, M. Christensen, R. Fernandes, L. Fu, D. Goldhaber-Gordon, P. Jarillo-Herrero, J. Kang, A. Klein,L. Levitov, M. Navarro Gastiasoro, J. Schmalian, D. Sha ff er, O. Vafek, J. Venderbos, and A. Vishwanath for fruitful discus-sions. The work was supported by U.S. Department of Energy, O ffi ce of Science, Basic Energy Sciences, under Award No.de-sc0014402. A.V.C. is thankful to Aspen Center for Physics (ACP) for hospitality during the completion of this work.ACP is supported by NSF grant PHY-1607611. [1] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature , 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, et al. , Nature , 80 (2018).[3] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science ,1059 (2019), https: // science.sciencemag.org / content / / / et al. , Nature , 653(2019). [5] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature (2020),10.1038 / s41586-020-2260-6.[6] X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, arXiv e-prints , arXiv:1903.08130(2019), arXiv:1903.08130 [cond-mat.mes-hall].[7] C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, T. Taniguchi, R. Yang, Z. Y. Meng, D. Shi,O. V. Yazyev, and G. Zhang, Nature Physics , 520 (2020).[8] G. Chen, L. Jiang, S. Wu, B. Lyu, H. Li, B. L. Chittari, K. Watanabe, T. Taniguchi, Z. Shi, J. Jung, Y. Zhang, and F. Wang, NaturePhysics , 237 (2019).[9] G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, J. Jung, Z. Shi,D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Nature , 215 (2019).[10] E. Su´arez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B , 121407 (2010).[11] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences , 12233 (2011),https: // / content / / / , 031089 (2018).[13] J. Kang and O. Vafek, Phys. Rev. Lett. , 246401 (2019).[14] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X , 031087 (2018).[15] K. Seo, V. N. Kotov, and B. Uchoa, Phys. Rev. Lett. , 246402 (2019).[16] F. Wu, E. Hwang, and S. Das Sarma, Phys. Rev. B , 165112 (2019).[17] F. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. , 257001 (2018).[18] T. M. R. Wolf, J. L. Lado, G. Blatter, and O. Zilberberg, Phys. Rev. Lett. , 096802 (2019).[19] 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, Nature , 95 (2019).[20] Y. Kim, P. Herlinger, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, and J. H. Smet, Nano Letters , 5053 (2016), pMID:27387484, https: // doi.org / / acs.nanolett.6b01906.[21] G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong, and E. Y. Andrei, Nature Physics , 109 (2010).[22] K. L. Hur and T. M. Rice, Annals of Physics , 1452 (2009), july 2009 Special Issue.[23] R. Samajdar and M. S. Scheurer, “Microscopic pairing mechanism, order parameter, and disorder sensitivity in moir superlattices:Applications to twisted double-bilayer graphene,” (2020), arXiv:2001.07716 [cond-mat.supr-con].[24] B. Lian, Z. Wang, and B. A. Bernevig, Phys. Rev. Lett. , 257002 (2019).[25] Y. W. Choi and H. J. Choi, Phys. Rev. B , 241412 (2018).[26] T. J. Peltonen, R. Ojaj¨arvi, and T. T. Heikkil¨a, Phys. Rev. B , 220504 (2018).[27] J. W. F. Venderbos and R. M. Fernandes, Phys. Rev. B , 245103 (2018).[28] H. Isobe, N. F. Q. Yuan, and L. Fu, Phys. Rev. X , 041041 (2018).[29] S. Ray, J. Jung, and T. Das, Phys. Rev. B , 134515 (2019).[30] J. Gonz´alez and T. Stauber, Phys. Rev. Lett. , 026801 (2019).[31] Y. Sherkunov and J. J. Betouras, Phys. Rev. B , 205151 (2018).[32] Y.-Z. You and A. Vishwanath, npj Quantum Materials , 1 (2019).[33] Y.-P. Lin and R. M. Nandkishore, Phys. Rev. B , 214521 (2018).[34] Y.-P. Lin and R. M. Nandkishore, Phys. Rev. B , 085136 (2019).[35] W. Chen, Y. Chu, T. Huang, and T. Ma, Phys. Rev. B , 155413 (2020).[36] X.-C. Wu, K. A. Pawlak, C.-M. Jian, and C. Xu, arXiv e-prints , arXiv:1805.06906 (2018), arXiv:1805.06906 [cond-mat.str-el].[37] J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q. Sun, and C. Wang, Phys. Rev. B , 075154 (2018).[38] C. Xu and L. Balents, Phys. Rev. Lett. , 087001 (2018).[39] M. Fidrysiak, M. Zegrodnik, and J. Spałek, Phys. Rev. B , 085436 (2018).[40] E. Laksono, J. N. Leaw, A. Reaves, M. Singh, X. Wang, S. Adam, and X. Gu, Solid State Communications , 38 (2018).[41] C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Phys. Rev. Lett. , 217001 (2018).[42] Y. Su and S.-Z. Lin, Phys. Rev. B , 195101 (2018).[43] Z. Liu, Y. Li, and Y.-F. Yang, Chinese Physics B , 077103 (2019).[44] D. M. Kennes, J. Lischner, and C. Karrasch, Phys. Rev. B , 241407 (2018).[45] L. Classen, C. Honerkamp, and M. M. Scherer, Phys. Rev. B , 195120 (2019).[46] R. Nandkishore, L. Levitov, and A. Chubukov, Nature Physics , 158 (2012).[47] M. L. Kiesel, C. Platt, W. Hanke, D. A. Abanin, and R. Thomale, Phys. Rev. B , 020507 (2012).[48] P. Jarillo-Herrero, Talk at KITP Rapid Response Workshop .[49] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Nature , 91 (2019).[50] Y. S. Kushnirenko, D. V. Evtushinsky, T. K. Kim, I. V. Morozov, L. Harnagea, S. Wurmehl, S. Aswartham, A. V. Chubukov, andS. V. Borisenko, arXiv e-prints , arXiv:1810.04446 (2018), arXiv:1810.04446 [cond-mat.str-el].[51] Y. S. Hor, A. J. Williams, J. G. Checkelsky, P. Roushan, J. Seo, Q. Xu, H. W. Zandbergen, A. Yazdani, N. P. Ong, and R. J. Cava,Phys. Rev. Lett. , 057001 (2010).[52] L. Fu and E. Berg, Phys. Rev. Lett. , 097001 (2010).[53] L. Fu, Phys. Rev. B , 100509 (2014). [54] K. Matano, M. Kriener, K. Segawa, Y. Ando, and G.-q. Zheng, Nature Physics , 852 (2016).[55] J. W. F. Venderbos, V. Kozii, and L. Fu, Phys. Rev. B , 180504 (2016).[56] S. Yonezawa, K. Tajiri, S. Nakata, Y. Nagai, Z. Wang, K. Segawa, Y. Ando, and Y. Maeno, Nature Physics , 123 (2017).[57] M. Hecker and J. Schmalian, npj Quantum Materials , 26 (2018).[58] C.-w. Cho, J. Shen, J. Lyu, S. H. Lee, Y. San Hor, M. Hecker, J. Schmalian, and R. Lortz, arXiv e-prints , arXiv:1905.01702(2019), arXiv:1905.01702 [cond-mat.supr-con].[59] X. Wu, W. Hanke, M. Fink, M. Klett, and R. Thomale, Phys. Rev. B , 134517 (2020).[60] V. Kozii, H. Isobe, J. W. F. Venderbos, and L. Fu, Phys. Rev. B , 144507 (2019).[61] N. F. Q. Yuan and L. Fu, Phys. Rev. B , 045103 (2018).[62] M. S. Scheurer and R. Samajdar, “Pairing in graphene-based moir superlattices,” (2019), arXiv:1906.03258 [cond-mat.supr-con].[63] N. F. Q. Yuan, H. Isobe, and L. Fu, Nature Communications , 5769 (2019).[64] X. Wu, M. Fink, W. Hanke, R. Thomale, and D. Di Sante, Phys. Rev. B , 041117 (2019).[65] J. Kang and O. Vafek, Phys. Rev. X , 031088 (2018).[66] M. Hamermesh, Group theory and its application to physical problems (Courier Corporation, 2012).[67] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. , 256802 (2007).[68] E. J. Mele, Phys. Rev. B , 161405 (2010).[69] E. J. Mele, Phys. Rev. B , 235439 (2011).[70] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B , 155449 (2012).[71] P. Moon and M. Koshino, Phys. Rev. B , 195458 (2012).[72] P. Moon and M. Koshino, Phys. Rev. B , 205404 (2013).[73] M. Koshino and P. Moon, Journal of the Physical Society of Japan , 121001 (2015), https: // doi.org / / JPSJ.84.121001.[74] N. N. T. Nam and M. Koshino, Phys. Rev. B , 075311 (2017).[75] S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov, Phys. Rev. B , 165105 (2010).[76] G. Trambly de Laissardire, D. Mayou, and L. Magaud, Nano Letters , 804 (2010), pMID: 20121163,https: // doi.org / / nl902948m.[77] J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Phys. Rev. B , 205414 (2014).[78] S. Fang and E. Kaxiras, Phys. Rev. B , 235153 (2016).[79] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Phys. Rev. B , 195455 (2019).[80] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Phys. Rev. B , 085435 (2018).[81] G. Trambly de Laissardi`ere, D. Mayou, and L. Magaud, Phys. Rev. B , 125413 (2012).[82] Y. Cao, J. Y. Luo, V. Fatemi, S. Fang, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero,Phys. Rev. Lett. , 116804 (2016).[83] H. Isobe and L. Fu, Phys. Rev. Research , 033206 (2019).[84] S. Maiti and A. V. Chubukov, AIP Conference Proceedings , 3 (2013), https: // aip.scitation.org / doi / pdf / / , 129 (1994).[86] W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Sch¨onhammer, Rev. Mod. Phys. , 299 (2012).[87] C. Platt, W. Hanke, and R. Thomale, Advances in Physics , 453 (2013), https: // doi.org / / , 125434 (2013).[89] M. Sigrist and K. Ueda, Rev. Mod. Phys. , 239 (1991).[90] B. Zinkl, M. H. Fischer, and M. Sigrist, Phys. Rev. B , 014519 (2019).[91] K.-S. Liu and M. E. Fisher, Journal of Low Temperature Physics , 655 (1973).[92] P. Calabrese, A. Pelissetto, and E. Vicari, Phys. Rev. B , 054505 (2003).[93] A. Eichhorn, D. Mesterh´azy, and M. M. Scherer, Phys. Rev. E , 042141 (2013).[94] A. V. Chubukov and D. I. Golosov, Journal of Physics: Condensed Matter , 69 (1991).[95] M. I. B. Utama, R. J. Koch, K. Lee, N. Leconte, H. Li, S. Zhao, L. Jiang, J. Zhu, K. Watanabe, T. Taniguchi, P. D. Ashby,A. Weber-Bargioni, A. Zettl, C. Jozwiak, J. Jung, E. Rotenberg, A. Bostwick, and F. Wang, “Visualization of the flat electronicband in twisted bilayer graphene near the magic angle twist,” (2019), arXiv:1912.00587 [cond-mat.mes-hall].[96] J. Carlstr¨om, J. Garaud, and E. Babaev, Phys. Rev. B , 134518 (2011).[97] A. Kapitulnik, Physica B: Condensed Matter , 151 (2015), special Issue on Electronic Crystals (ECRYS-2014).[98] Z. L. Mahyari, A. Cannell, C. Gomez, S. Tezok, A. Zelati, E. V. L. de Mello, J.-Q. Yan, D. G. Mandrus, and J. E. Sonier, Phys.Rev. B , 020502 (2014).[99] X. Hu and Z. Wang, Phys. Rev. B , 064516 (2012).[100] J. Garaud, M. Silaev, and E. Babaev, Physica C: Superconductivity and its Applications , 63 (2017), ninth internationalconference on Vortex Matter in nanostructured Superdonductors.[101] S.-Z. Lin, S. Maiti, and A. Chubukov, Phys. Rev. B , 064519 (2016).[102] J. Garaud and E. Babaev, Phys. Rev. Lett. , 017003 (2014).[103] A. A. Zyuzin, J. Garaud, and E. Babaev, Phys. Rev. Lett. , 167001 (2017). Supplemental Material
S1. DIAGONALIZATION OF SINGLE-PARTICLE HAMILTONIAN AND THE INTRODUCTION OF PATCHOPERATORS
We introduce the following lattice vectors (see Fig. 2) a =
12 ( − , √ , a =
12 ( − , − √ , a = (1 , , (S1) b =
32 ( − , √ , b =
32 ( − , − √ , b = (3 , , (S2) L =
12 (3 , √ , L = (0 , √ . (S3)and Fourier-transform the Hamiltonian Eq. (1). As a result, we obtain H T B = (cid:88) k (cid:16) c † xA c † xB c † yA c † yB (cid:17) T d T ∗ sd iT sd T sd T d iT sd − iT sd T d T ∗ sd − iT sd T sd T d c xA c xB c yA c yB , (S4)which we diagonalize to determine the single-particle bands given in Eq. (4) in the main text. Note, that we used spinlessfermions in (S4) reflecting the spin degeneracy of bands in the absence of spin-orbit coupling. To diagonalize Eq. (S4) weuse the transformation U − H T B U = Diag( E −− , E + − , E − + , E ++ ) with U = √| T sd | i (cid:112) T ∗ sd − i (cid:112) T ∗ sd − i (cid:112) T ∗ sd i (cid:112) T ∗ sd − i √ T sd i √ T sd − i √ T sd i √ T sd − (cid:112) T ∗ sd − (cid:112) T ∗ sd (cid:112) T ∗ sd (cid:112) T ∗ sd √ T sd √ T sd √ T sd √ T sd (S5)and the eigenvalues are given in Eq. (4).With the orbtial-to-band transformation U we transform the fermionic spinor c = (cid:16) c † xA c † xB c † yA c † yB (cid:17) written in the Wannierorbital basis to the spinor d = (cid:16) d † d † d † d † (cid:17) written in band basis by a standard linear transformation c = U d . (S6)Note, that spinor c gains momentum dependence from the transformation U . With this now we have an explicit relationbetween the orbital and the band operators, hence we can rewrite the interaction term in the band basis. As U depends onmomentum, this yields the so-called orbital makeup, i.e. additional momentum-dependent factors (which are sometimescalled “coherence factors”) in the interaction from the transformation.In our model we are interested in fermions located at patches near the VH points. So as a next step, we introduce patchband operators, which “live” at the positions of the VH singularities. More specifically, we introduce patch operators f i , f † i with i = ..
12 being the patch index defined as f i = d j ( k i ) , f † i = d † j ( k i ) , (S7)where k i are the positions of VH singularities in the Brillouin zone and j is the band index. Note, that j is di ff erent fordi ff erent patches, since the VH singularities are made of di ff erent bands. The coupling constants can then be obtained byan explicit calculation in which the coupling functions are translated to momentum space and multiplied by four coherencefactors related to fermions on the involved patches. S2. EXPLICIT FORM OF COUPLING FUNCTIONS
In this section we present the explicit form of coupling functions from Eq. (15). In the equations below we use the samevector notation as in Eqs. (S1),(S2),(S3). QQ − term coupling functions, which depend only on transferred momentum, aregiven by F AAAA = (cid:16) + e i L ( k − q ) + e i L ( k − q ) (cid:17) (cid:16) + e i L ( k (cid:48) − q (cid:48) ) + e i L ( k (cid:48) − q (cid:48) ) (cid:17) , F ABBA = (cid:16) + e i L ( k − q ) + e i L ( k − q ) (cid:17) (cid:16) e i a ( k (cid:48) − q (cid:48) ) + e i ( L + a )( k (cid:48) − q (cid:48) ) + e i ( L + a )( k (cid:48) − q (cid:48) ) (cid:17) , F BAAB = (cid:16) e i a ( k − q ) + e i ( L + a )( k − q ) + e i ( L + a )( k − q ) (cid:17) (cid:16) + e i L ( k (cid:48) − q (cid:48) ) + e i L ( k (cid:48) − q (cid:48) ) (cid:17) , F BBBB = (cid:16) e i a ( k − q ) + e i ( L + a )( k − q ) + e i ( L + a )( k − q ) (cid:17) (cid:16) e i a ( k (cid:48) − q (cid:48) ) + e i ( L + a )( k (cid:48) − q (cid:48) ) + e i ( L + a )( k (cid:48) − q (cid:48) ) (cid:17) . (S8)Couplings of T QQT − and T T − terms can also be written in an explicit form. For T QQT we get F AABA = (cid:16) + e i L ( k − q ) + e i L ( k − q ) (cid:17) (cid:16) e − i δ q (cid:48) + e i L k (cid:48) e − i ( L + a ) q (cid:48) + e i L k (cid:48) e − i ( L + a ) q (cid:48) − e i L k (cid:48) e − i a q (cid:48) − e i L k (cid:48) e − i ( L + a ) q (cid:48) − e − i ( L + a ) q (cid:48) (cid:17) , F ABAA = (cid:16) + e i L ( k − q ) + e i L ( k − q ) (cid:17) (cid:16) − e − i a k (cid:48) e i L q (cid:48) − e i ( L + a ) k (cid:48) e − i L q (cid:48) − e − i ( L + a ) k (cid:48) + e i a k (cid:48) + e i ( L + a ) k (cid:48) e − i L q (cid:48) + e i ( L + a ) k (cid:48) e − i L q (cid:48) (cid:17) , F AAAB = F AABA ( k ↔ k (cid:48) , q ↔ q (cid:48) ) , F BAAA = F ABAA ( k ↔ k (cid:48) , q ↔ q (cid:48) ) , F BABB = (cid:16) e i a ( k − q ) + e i ( L + a )( k − q ) + e i ( L + a )( k − q ) (cid:17) (cid:16) e − i a q (cid:48) + e i L k (cid:48) e − i ( L + a ) q (cid:48) + e i L k (cid:48) e − i ( L + a ) q (cid:48) − e i L k (cid:48) e − i a q (cid:48) − e i L k (cid:48) e − i ( L + a ) q (cid:48) − e − i ( L + a ) q (cid:48) (cid:17) , F BBAB = (cid:16) e i a ( k − q ) + e i ( L + a )( k − q ) + e i ( L + a )( k − q ) (cid:17) (cid:16) − e − i a k (cid:48) e i L q (cid:48) − e i ( L + a ) k (cid:48) e − i L q (cid:48) − e − i ( L + a ) k (cid:48) + e i a k (cid:48) + e i ( L + a ) k (cid:48) e − i L q (cid:48) + e i ( L + a ) k (cid:48) e − i L q (cid:48) (cid:17) , F ABBB = F BABB ( k ↔ k (cid:48) , q ↔ q (cid:48) ) , F BBBA = F BBAB ( k ↔ k (cid:48) , q ↔ q (cid:48) ) , (S9)while T T − couplings are F ABAB = (cid:16) e − i a q + e i L k e − i ( L + a ) q + e i L k e − i ( L + a ) q (cid:48) − e i L k e − i a q − e i L k e − i ( L + a ) q − e − i ( L + a ) q (cid:17) × (cid:16) − e − i a k (cid:48) e i L q (cid:48) − e i ( L + a ) k (cid:48) e − i L q (cid:48) − e − i ( L + a ) k (cid:48) + e i a k (cid:48) + e i ( L + a ) k (cid:48) e − i L q (cid:48) + e i ( L + a ) k (cid:48) e − i L q (cid:48) (cid:17) , (S10) F AABB = (cid:16) e − i a q + e i L k e − i ( L + a ) q + e i L k e − i ( L + a ) q (cid:48) − e i L k e − i a q − e i L k e − i ( L + a ) q − e − i ( L + a ) q (cid:17) × (cid:16) e − i a q (cid:48) + e i L k (cid:48) e − i ( L + a ) q (cid:48) + e i L k (cid:48) e − i ( L + a ) q (cid:48) − e i L k (cid:48) e − i a q (cid:48) − e i L k (cid:48) e − i ( L + a ) q (cid:48) − e − i ( L + a ) q (cid:48) (cid:17) , (S11) F BABA = (cid:16) − e − i a k e i L q − e i ( L + a ) k e − i L q − e − i ( L + a ) k + e i a k + e i ( L + a ) k e − i L q + e i ( L + a ) k e − i L q (cid:17) × (cid:16) e − i a q (cid:48) + e i L k (cid:48) e − i ( L + a ) q (cid:48) + e i L k (cid:48) e − i ( L + a ) q (cid:48) − e i L k (cid:48) e − i a q (cid:48) − e i L k (cid:48) e − i ( L + a ) q (cid:48) − e − i ( L + a ) q (cid:48) (cid:17) , (S12) F BBAA = (cid:16) − e − i a k e i L q − e i ( L + a ) k e − i L q − e − i ( L + a ) k + e i a k + e i ( L + a ) k e − i L q + e i ( L + a ) k e − i L q (cid:17) × (cid:16) − e − i a k (cid:48) e i L q (cid:48) − e i ( L + a ) k (cid:48) e − i L q (cid:48) − e − i ( L + a ) k (cid:48) + e i a k (cid:48) + e i ( L + a ) k (cid:48) e − i L q (cid:48) + e i ( L + a ) k (cid:48) e − i L q (cid:48) (cid:17) . (S13) S3. HUBBARD-STRATONOVICH DERIVATION AND THE SYMMETRY PROPERTIES OF THE LANDAU FREEENERGY
In this section we provide an explicit derivation of the Landau free energy using the Hubbard-Stratonovich transformationand discuss symmetry properties of the functional. We first derive the free energy. The derivation for the six-patch modelis analogous to the calculation presented in the Supplementary material for [46]. Here, we present the derivation for thetwelve-patch model. We begin by writing down the Lagrangian L = (cid:88) i f † i ( ∂ τ − (cid:15) ik ) f i − (cid:16) f † f † f † f † f † f † f † f † f † f † f † f † (cid:17) v g − g g g g + g − v g + g g g g g + v g − g g g g g − v g + g g g g g + v g − g + g g g g − v f f f f f f f f f f f f , (S14)where we absorbed the chemical potential µ into the dispersion (cid:15) k , i - is the patch index running from 1 to 6 and theother indices were omitted for shortness. The spin and band structure of each term is f † a τσ f † a ¯ τ ¯ σ f a ¯ τ ¯ σ f a τσ , where the barlabels opposite spin or band, i.e. we consider spin-singlet pairing with zero total momentum. Note, that due to the D symmetry of our low-energy model fermionic dispersions are identical for all patches (upon rotations). We decompose thepairing interaction in Eq. (S14) into its eigenvalues and eigenvectors. A pairing instability can develop when at least oneeigenvalue is positive and the dominant pairing channel is determined by the largest eigenvalue. We find that two eigenvaluescan become positive. They belong to either the A or E representation of D . As A is a one-dimensional irreduciblerepresentation, the corresponding eigenvalue is unique, while the eigenvalue of the two-dimensional E representation is two-fold degenerate. We discard subleading channels in the following, which is justified close to the largest critical temperaturebecause critical temperatures of subleading channels are much smaller. In contrast, the critical temperatures of A and E order are very similar. We can classify the eigenvectors corresponding to the leading pairing instabilities according to thelattice harmonics. Interestingly, the eigenvectors appear in higher “angular” momentum channels in our case following g E ( k x , k y ) = k x − k x k y + k y = cos 4 θ , g E ( k x , k y ) = k x k y ( k x − k y ) = sin 4 θ and g A ( k x , k y ) = k x k y − k x k y + k x k y = sin 6 θ with θ = arctan k y / k x . That is we find the eigenvectors ∆ A p = (cid:0) g A ( k ) , . . . , g A ( k ) (cid:1) ∆ E − p = (cid:0) g E ( k ) , . . . , g E ( k ) (cid:1) ∆ E − p = (cid:0) g E ( k ) , . . . , g E ( k ) (cid:1) (S15)with the six VH points k , . . . , k . In the following, we will use the orthonormal expressions ∆ Γ p → ∆ Γ p / (cid:107) ∆ Γ p (cid:107) . Letus also note that we can choose any linear, orthonormal combination of ∆ E − p and ∆ E − p in the two-dimensional subspacecorresponding to the degenerate eigenvalue, i.e. ∆ E − p → cos( α ) ∆ E − p − sin( α ) ∆ E − p , ∆ E − p → sin( α ) ∆ E − p + cos( α ) ∆ E − p . Thiscorresponds to rotation within the E -subgroup and the resulting eigenvectors with proper normalization can be representedas ∆ A p = N − A (cid:0) sin(6 θ ) , . . . , sin(6 θ ) (cid:1) T ∆ E − p = N − E (cid:0) cos(4 θ + α ) , . . . , cos(4 θ + α ) (cid:1) T ∆ E − p = N − E (cid:0) sin(4 θ + α ) , . . . , sin(4 θ + α ) (cid:1) T , (S16)with N − A , E being the normalization coe ffi cients. Explicitly, for a particular choice of hopping parameters (see caption forFig. 3), by numerical diagonalization of the linearized gap equation we obtain ∆ A p = √ ( − , , − , , − , , ∆ E − p = ( − . , . , . , − . , − . , . ∆ E − p = (0 . , − . , − . , . , − . , . . To an excellent ap-proximation these gap values are fitted by single harmonics given in Eq. S16 with α = π/
4. We use these three vectorsfurther to perform the Hubbard-Stratonovich transformation.In the bulk text we have provided a solution which shows that at a critical value of the parameter α T , all three channelsbecome degenerate. Close to this value, the pairing function can be expressed as a linear combination of the eigenvectors ∆ SC p = ∆ A ∆ A p + ∆ E ∆ E − p + ∆ E ∆ E − p . (S17)We are interested in the Landau action in the vicinity of exactly this point, which can be constructed as an expansion in ∆ A , ∆ E and ∆ E . As we said in the main text, we can anticipate the form of the Landau action from symmetry properties;it has to be invariant under U (1) and C transformations. At the same time, we can derive it explicitly from the fermionicLagrangian, which also gives us the bare values of the couplings in the Landau action.For the following derivation, it is convenient to introduce three 6 × d A , d E , and d E in patch space. Thesematrices are diagonal and given by d A = diag( ∆ A p ) d E = diag( ∆ E − p ) d E = diag( ∆ E − p ) (S18)The order parameters then can be defined via these matrices ∆ A = λ (cid:68) f T d A f (cid:69) , ∆ E = λ (cid:68) f T d E f (cid:69) , ∆ E = λ (cid:68) f T d E f (cid:69) , (S19)where f T = ( f , . . . , f ) and λ is the triply degenerate eigenvalue. Using these bosonic fields we can now apply the Hubbard-Stratonovich transformation to Eq. (S14) and obtain the Lagrangian in terms of bosonic and fermionic fields written in theNambu-Gor’kov notation: L = (cid:16) f † f (cid:17) (cid:32) G − + ∆ A d A + ∆ E d E + ∆ E d E ∆ ∗ A d A + ∆ ∗ E d E + ∆ ∗ E d E G − − (cid:33) (cid:32) ff † (cid:33) + | ∆ A | + | ∆ E | + | ∆ E | λ , (S20)where G ± are particle and hole Green’s functions with G ± = diag( G , ± , . . . , G , ± ). The Green’s functions are given by G − i , ± = i ω ∓ (cid:15) ik , where ω is the Matsubara frequency and i is the patch number index. Now, we can integrate out the fermionsand obtain the Lagrangian written in terms of only bosonic fields L = Tr ln (cid:32) G − + ∆ A d A + ∆ E d E + ∆ E d E ∆ ∗ A d A + ∆ ∗ E d E + ∆ ∗ E d E G − − (cid:33) + | ∆ A | + | ∆ E | + | ∆ E | λ , (S21)where the trace is going over the patch and Nambu-Gor’kov spinor spaces (the overall matrix dimension in (S21) is twelve).The trace over the patch space is an analog of the momentum integration in the continuous notation. We exploit the propertythat the Green’s functions commute with the order parameter matrices and expand Eq. (S21) in small ∆ A , ∆ E , ∆ E up tothe fourth order. As a result, we obtain F = F + F F = A (cid:16) | ∆ E | + | ∆ E | (cid:17) + A | ∆ A | F = K (cid:20) (cid:32) | ∆ E | + | ∆ E | + | ∆ A | (cid:33) + (cid:16) | ∆ E | | ∆ E | + | ∆ A | | ∆ E | + | ∆ A | | ∆ E | (cid:17) + (cid:16) ∆ E ¯ ∆ E + ∆ E ¯ ∆ E + ∆ A ¯ ∆ E + ∆ E ¯ ∆ A + ∆ A ¯ ∆ E + ∆ E ¯ ∆ A (cid:17) + δ (cid:104)(cid:16) − ∆ E | ∆ E | − ¯ ∆ E ∆ E + ∆ E | ∆ E | (cid:17) ¯ ∆ A + c . c (cid:105) + δ (cid:104)(cid:16) ∆ E | ∆ E | + ¯ ∆ E ∆ E − ∆ E | ∆ E | (cid:17) ¯ ∆ A + c . c (cid:105) , (S22)where the upper bar means complex conjugation and K = (cid:82) d k / (2 π ) G i , + G i , − G i , + G i , − is the fermionic box diagram withmomentum integration constrained to the vicinity of the patches. Note that because of rotation and time-reversal symme-try the constant K is the same for all patches. Coe ffi cients in Eq. (S22) correspond to the ones mentioned in the maintext in the following way: β = / , β = / , β = / , γ = / , γ = / . To obtain Eq. (S22) we used the iden-tities Tr( d E ) = Tr( d E ) = / , Tr( d E d E ) = Tr( d E d E d E d E ) = / , Tr( d A d E ) = Tr( d A d E ) = Tr( d A d E d A d E ) = Tr( d A d E d A d E ) = Tr( d A ) = / , and the fact that most traces over the patch space of the cube of one matrix times theother matrix vanishes. However, there are exceptions for special C - and U (1)-symmetric combinations, which lead tothe last two lines in Eq. (S22). We obtain for their prefactor δ = − δ ≈ − .
128 for the choice of the basis for the two-dimensional representation E as in Eq. (S16) with α = π/
4. Note that in contrast to the other coupling constants, δ and δ depend on a basis change of the two-dimensional representation, because the polynomials cubic in ∆ E , ∆ E also transformnon-trivially under a basis rotation. Of course, in total, the action is invariant. Explicitly, if we rotate our basis via (cid:32) ∆ (cid:48) E ∆ (cid:48) E (cid:33) = R α (cid:32) ∆ E ∆ E (cid:33) R α = (cid:32) cos α − sin α sin α cos α (cid:33) (S23)the polynomials p = ∆ E | ∆ E | − ∆ E | ∆ E | − ¯ ∆ E ∆ E , p = − ∆ E | ∆ E | + ∆ E | ∆ E | + ¯ ∆ E ∆ E and δ, δ transform accordingto (cid:16) δ (cid:48) δ (cid:48) (cid:17) = (cid:16) δ δ (cid:17) R α (cid:32) p (cid:48) p (cid:48) (cid:33) = R T α (cid:32) p p (cid:33) (S24)so that their product in F is invariant. This allows us to choose our basis so that δ = − δ2