Pairing in graphene-based moiré superlattices

We present a systematic classiﬁcation and analysis of possible pairing instabilities in graphene-based moiré superlattices. Motivated by recent experiments on twisted double bilayer graphene showing signs of triplet superconductivity, we analyze both singlet and triplet pairing separately, and describe how these two channels behave close to the limit where the system is invariant under separate spin rotations in the two valleys, realizing an SU(2) + × SU(2) − symmetry. Further, we discuss the conditions under which singlet and triplet can mix via two nearly degenerate transitions, and how the different pairing states behave when an external magnetic ﬁeld is applied. The consequences of the additional microscopic or emergent approximate symmetries relevant for superconductivity in twisted bilayer graphene and ABC trilayer graphene on hexagonal boron nitride are described in detail. We also analyze which of the pairing states can arise in mean-ﬁeld theory and study the impact of corrections coming from ferromagnetic ﬂuctuations. For instance, we show that, close to the parameters of mean-ﬁeld theory, a nematic mixed singlet-triplet state emerges. Our paper illustrates that graphene superlattices provide a rich platform for exotic superconducting states, and allow for the admixture of singlet and triplet pairing even in the absence of spin-orbit coupling.


I. INTRODUCTION
Experiments on twisted bilayers of graphene have recently revealed interaction-induced insulating phases and superconductivity when the relative angle between the layers is fine tuned to yield almost flat moiré bands, which enhances the impact of electronic correlations [1][2][3][4].Due to the strongcoupling nature of the problem, which is corroborated by tunneling spectroscopy measurements [5][6][7][8][9], the form and mechanism of the insulating and superconducting phases are still under debate, despite considerable theoretical effort .Another graphene-based moiré system that displays both superconducting and correlated insulating behavior is ABC-stacked trilayer graphene on hexagonal boron nitride [47,48].In this case, the moiré pattern results from the difference in lattice constants, and it can be controlled by application of a vertical electric field [49,50].
The most recent member of the family of strongly correlated graphene superlattice systems is twisted double bilayer graphene [51][52][53], where two individually aligned ABstacked graphene bilayers are twisted with respect to one another.As theoretical calculations show [54][55][56][57][58][59][60], flat electronic bands can be realized by tuning the twist angle and a vertical electric field.Similar to the above-mentioned graphene moiré systems, both correlated insulating [51][52][53] and superconducting [51,52] phases are observed in experiment.However, in stark contrast to twisted bilayer and trilayer graphene, the Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
superconducting transition temperature is found to increase linearly with a weak in-plane magnetic field [52], which is a strong indication of triplet pairing [42,57].Furthermore, the gap of the correlated insulating phase is seen to increase with an applied magnetic field, indicating ferromagnetic order [51][52][53].There are also clear experimental indications of ferromagnetism in twisted bilayer [4,61,62] and ABC trilayer graphene [63].
In this paper, we study the possible pairing states in graphene moiré superlattices.Motivated by the recent experimental signatures of triplet pairing, we pay special attention to the triplet channel, and possible mixed singlet and triplet phases.While the weak spin-orbit coupling in graphene seems to disfavor the latter class of phases, projections of the Coulomb interaction on the relevant moiré bands evince that the interaction terms that couple the spin degrees of freedom of the two valleys, v = ±, of the system are much weaker than other interaction terms that do not [13,50,57].Together with the nearly valley-diagonal band structure, this indicates that the system is approximately invariant under independent spin rotations in the two valleys.As has been pointed out before [19,44], the associated SU(2) + × SU(2) − symmetry renders the singlet and triplet pairing channels degenerate.This paper will address the questions of (i) under which conditions can singlet and triplet mix when the SU(2) + × SU(2) − symmetry is only weakly broken and (ii) which triplet state transforms into which singlet upon reversing the sign of the symmetrybreaking interactions.In this way, we map out all possible phase diagrams close to the SU(2) + × SU(2) − invariant limit.
In light of the narrow bandwidth and strong correlations of the graphene-based moiré superlattices we are interested in, our analysis will begin with a comprehensive study of exact constraints resulting from symmetry.The symmetry-based classification will then be supplemented with energetics, by FIG. 1. Geometry (top), lattice symmetries (bottom left), and Brillouin zone with symmetries (bottom right) for (a) twisted double bilayer graphene, (b) twisted bilayer graphene, and (c) ABC-stacked trilayer graphene on hexagonal boron nitride.In all three cases, we assume a commensurate superlattice structure for simplicity.For (a) twisted double bilayer, we only show ABAB stacking, which exhibits a π -rotation symmetry C 2x along the x axis.For ABBA stacking, we have a C 2y symmetry instead.Our symmetry analysis of pairing applies to both stacking orders as the applied electric field breaks both of these in-plane rotation symmetries.In (b) twisted bilayer graphene, superconductivity emerges without any applied electric field and, hence, the C 2y symmetry of the lattice has to be taken into account.In addition, the system is approximately invariant under a C 6 symmetry [15] as indicated in gray in the Brillouin zone in (b).In (c), we show only the top boron nitride layer and one of the graphene layers for image clarity in the main panel; the other two graphene layers are indicated in the close-up view of the lattice in light blue.We assume no additional twist such that the moiré pattern is solely due to the lattice mismatch.This leads to the reflection symmetry σ xz , which for the effective two-dimensional description of the system can be viewed as an in-plane rotation C 2y as is done in the main text.Also, this system is believed to exhibit an approximate C 6 symmetry [50] as indicated in gray in the respective Brillouin zone.
studying which of the pairing states can be realized in the weak-coupling limit and what changes in the presence of additional fluctuation corrections.As it has the smallest set of symmetries, we will begin our classification with twisted double bilayer graphene: while the lattice is invariant under threefold rotation, C 3 , perpendicular to the graphene sheets, and under a twofold in-plane rotation [see Fig. 1(a)], the latter is broken due to the vertical electric field that is applied to tune the band structure and to induce superconductivity.It seems currently unclear whether the superconducting state coexists with the likely ferromagnetic correlated insulator and whether, at least in part of the phase diagram, there is a thermal transition directly from the (paramagnetic) normal metal to superconductivity without any ferromagnetic order.For this reason, we will analyze two scenarios separately: (I) there is no ferromagnetic order around the critical temperature, T c , of superconductivity; and (II) there is ferromagnetic order already at T > T c that coexists microscopically with superconductivity for T < T c (or, at the minimum, the associated ferromagnetic moments couple significantly to the superconducting order parameter).We will begin with the analysis of the superconducting states transforming under the irreducible representations (IRs) of the point group C 3 assuming time-reversal symmetry in the high-temperature phase-this is relevant for case I above.In order to capture scenario II, we will later add the coupling to the time-reversal symmetry-breaking magnetic moments and examine how it affects the superconducting transition.This allows us to determine which of the pairing states are compatible with the linear increase of the critical temperature with small magnetic field, B, and to describe the possible phase diagrams in the temperature-B plane.
We also generalize our discussion to twisted bilayer graphene and ABC trilayer graphene on hexagonal boron nitride.Here, we have to take into account an additional twofold rotation symmetry, C 2 , perpendicular to the plane of the system and an in-plane rotation symmetry C 2y along the y axis; these symmetries are either realized as exact microscopic symmetries of the lattice or as approximate emergent symmetries [15,17,50] of those systems [see Figs.1(b) and 1(c)].

A. Brief summary of the main results
Due to the length of the paper, here, we provide a very concise summary of the key results of this paper for the convenience of the reader.
(1) We analyze the consequences of the enhanced SU(2) + × SU(2) − spin symmetry, taking into account the possibility of several consecutive superconducting transitions with their difference in transition temperatures vanishing in the limit where SU(2) + × SU(2) − becomes exact.The resulting complete sets of possible phases for the relevant symmetry groups C 3 and D 3 (or, equivalently, D 6 , see Sec.VI A) are summarized in Tables I, II, and IV.
(2) As follows from these tables, all point groups and all of their IRs allow for singlet-triplet admixed phases in the absence of any spin-orbit coupling.As opposed to the conventional mechanism of singlet-triplet admixture, which is based on a reduced symmetry [64], here, it results from (the proximity to) an enhanced spin symmetry.
(3) To supplement these purely symmetry-based considerations with energetics, we analyze which of those states can be realized in single-band mean-field theory, i.e., whether there exists a form of the effective electron-electron interaction that can stabilize the superconducting state when treated within the mean-field approximation; the result is indicated in the last column in Tables I, II, and IV.This identifies the most important pairing states from a weak-coupling perspective.The presence of any of the remaining pairing phases-as might eventually be established in future experiments-must result from the strong-coupling and/or interband nature of superconductivity.
(4) We also study corrections to mean-field theory coming from ferromagnetic fluctuations, within a simplified phenomenological approach in Sec.V that is justified microscopically in Appendix B1.We analyze two limits.First, we consider the case of weak fluctuations in order to lift the residual degeneracies within mean-field theory.We find that out of the two possible phase diagrams for the IR E close to mean-field theory, shown in Fig. 5, the one in part (b) [part (a)] is favored for spin (orbital) ferromagnetic fluctuations.This reveals that a nematic mixed singlet-triplet phase is a natural candidate pairing phase in graphene moiré superlattices.Second, we analyze which pairing states are favored in the case where the fluctuation corrections dominate over the mean-field contributions (see last column in the tables mentioned above).
(5) We study the coupling of the superconducting states to the magnetic field, B, and examine which states can give rise to a linear increase of the critical temperature for small B: if SU(2) + × SU(2) − is broken significantly, triplet pairing has to dominate for B = 0 and there are only three possible triplet states as leading instabilities for B = 0.In the case where SU(2) + × SU(2) − is an approximate symmetry, even singlet pairing at B = 0 can yield a linear increase.For instance, the possible phase diagrams in the presence of a magnetic field for pairing in the trivial IR A of C 3 are shown in Fig. 3.

B. Relation to other works
Let us briefly comment on the relation of this paper to other works in the literature.While our classification also contains the pure singlet states, which have been subject to intense scrutiny in twisted bilayer graphene, we are mainly interested in elucidating the consequences of the enhanced SU(2) + × SU(2) − spin symmetry with respect to subsequent transitions and the associated nontrivial interplay of singlet and triplet pairing.
In the context of twisted double bilayer graphene, where, recently, signs of triplet pairing have been discovered, [65] mainly focuses on the correlated insulating phase in this system, whereas [57] also discusses pairing.We extend the work of [57] by allowing for momentum-dependent pairing states, contrasting weakly and significantly broken SU(2) + × SU(2) − symmetry, investigating admixed singlet and triplet phases, analyzing fluctuation corrections to mean-field theory, and mapping out the phase diagram in the presence of a magnetic field.In a follow-up work [66], we will complement the analysis of this paper by a microscopic energetic study specifically for twisted double bilayer graphene.In [66], we discuss which IR is expected to be favored, the form of the associated basis functions, and the impact of disorder on superconductivity.

C. Structure of the paper
This paper is organized as follows.As described above, we start with twisted double bilayer graphene.In Sec.II, we introduce the model and the action of the relevant symmetries.We first discuss pairing in the trivial IR of the point group of the system in Sec.III and then generalize to the complex IR E in Sec.IV.Section V demonstrates how strong fluctuations can yield significant corrections to mean-field theory.We extend our analysis to twisted bilayer graphene and ABC trilayer graphene in Sec.VI, and explore the consequences of the additional microscopic and emergent symmetries relevant to those systems.A discussion of our results can be found in Sec.VII.

II. MODEL AND SYMMETRIES
We first focus on the (nearly flat) conduction band of twisted double bilayer graphene which appears to host the superconducting phase observed experimentally [51,52], and later discuss the modifications for the related moiré systems, bilayer and trilayer graphene.Owing to the presence of a gap to other bands in the relevant parameter regime [56][57][58][59], it is reasonable to describe the superconducting instability in a single-band picture.We stress, however, that many of our conclusions are symmetry based and, thus, also apply when several bands are taken into account.Exceptions are provided by the energetic mean-field and fluctuation considerations, where we will specifically comment on the consequences of interband effects that might be present in these systems [8].
Denoting the corresponding electronic creation and annihilation operators by c kσ v , where k is crystal momentum, σ is spin, and v = ± represents the valleys, the general pairing term can be written as Here and in the following, σ j and τ j are Pauli matrices in spin and valley space, respectively, and the 4 × 4 matrix k is the superconducting order parameter.In Eq. ( 1), we have already made the assumptions that only Cooper pairs with zero net momentum form and that superconductivity preserves translational symmetry.Due to the proximity of superconductivity to ferromagnetic order [51,52], relaxing this assumption could be interesting, but we leave this for future work.Consequently, we need not consider IRs of the full space group but rather can concentrate on the point group G of the system and time reversal .
In this regard, we study two different point groups: an approximate point group, where C 3 is the crystalline point group, SU(2) ± is spin rotation in valley v = ±, and U (1) v corresponds to valley charge conservation.As argued in [57], the intervalley "Hund's" coupling J is much smaller than the intravalley-density interaction V .In combination with the fact that the noninteracting band structure only has very small valley mixing, the system is invariant under Eq.( 2) to a good approximation.In the presence of a finite Hund's coupling, Eq. ( 2) is reduced to where SU(2) s is global spin rotation.To define these symmetries more precisely, we specify their representation on the electronic field operators: with P ± = (τ 0 ± τ z )/2 being the valley projection operators.Furthermore, time reversal is represented by the antiunitary operator with To classify superconductivity, we proceed as usual [67] and express k in Eq. ( 1) in terms of the IRs n (with dimension d n ) of the point group as where χ n μ (k), μ = 1, . . ., d n , are partner functions transforming under the IR n.Within the minimal description of pairing in Eq. (1), which only involves one band per valley, χ n μ (k) ∈ C 4×4 are matrices in spin and valley space.
In our case, the point group has the form and G s 2 = SU(2) s .As a consequence, the IRs of G j have the form n = n C 3 × n v × n s where n C 3 , n v , and n s are IRs of C 3 , U (1) v , and G s j , respectively.We can thus rewrite Eq. ( 6) more explicitly as In order to classify superconducting states, we need to consider the different IRs of C 3 , U (1) v , and G s j .Let us begin our discussion of IRs with U (1) v .While it has, in general, countably infinite IRs (one dimensional and with character e im v ϕ , m v ∈ Z), only three are relevant here as all representations with |m v | > 1 cannot be realized with only two valleys.First, there is the trivial representation, m v = 0, with χ m v =0 = aτ 0 + bτ z with a priori unknown a, b.Recalling the extra factor of τ x in Eq. ( 1), this translates to purely intervalley pairing.Secondly, the pair of complex conjugate representations with m v = ±1 has to be considered.Note that due to time-reversal symmetry the complex representations cannot be discussed separately.Here, the basis functions read as χ m v =±1 = τ x ± iτ y ; as such, this corresponds to purely intravalley pairing.
We thus see that U (1) v prohibits the mixing of interand intravalley pairing.As time reversal (5) interchanges the valleys along with sending k → −k and we assume zeromomentum Cooper pairs, we will restrict our discussion to intervalley pairing, i.e., m v = 0 for the rest of the paper.
As is well known [68], C 3 has the following IRs, both of which are one dimensional: the trivial one, A, and the complex representation E (and its complex conjugate partner).We analyze each of these IRs in Secs.III and IV and, in both cases, discuss the differences between G s 1 and G s 2 ; we will also see how the states "connect" once G s 1 is weakly broken to G s 2 due to a small but finite value of Hund's coupling.

III. TRIVIAL REPRESENTATION OF THE CRYSTALLINE POINT GROUP
For simplicity, we begin with the trivial representation A of C 3 , which is real and one dimensional.In fact, the following discussion will not be modified as long as the IR is real and one dimensional and there is no crystalline symmetry relating the two valleys.Interestingly, the last assumption is violated in twisted bilayer graphene and trilayer graphene on boron nitride; see Sec.VI for a detailed discussion of the associated modifications.
As already mentioned, we consider only intervalley pairing which corresponds to a real and one-dimensional IR as well.This means that the order parameter in Eq. ( 7) has the form where χ A (k, v) is invariant under k → gk for all generators g of the crystalline point group (here, we only have g = C 3 ).

A. Limit of exact SU(2) + × SU(2) − symmetry
To proceed further, we have to inspect the scenarios for both G s 1 and G s 2 .We start with the former, i.e., we assume that Hund's coupling is zero.Inserting Eq. ( 8) in the general pairing Hamiltonian (1), we obtain a pairing term of the form where v = ∓ for v = ±, and M kv as well as v = μ η n s μ χ n s μ (v) are matrices in spin space.Fermi-Dirac statistics implies Rewriting pairing in terms of singlet and triplet as We now study the stable superconducting phases in this channel by writing down the most general Ginzburg-Landau expansion constrained by the symmetries Due to the constraint (10) stemming from Fermi-Dirac statistics, we express the free energy in terms of one valley only (say v = +) as F = F[M k+ = χ A (k, +) + ], and the pairing in the other valley just follows from Eq. (10).The most general free energy to quartic order in + , invariant under Eq.( 11) and v → e iϕ v , reads as , so the latter is not an independent term to consider.It further holds that |tr[σ y + σ y + + ], which allows us to set b 3 = 0 in the following without loss of generality.
Using the singular-value decomposition of + , it is straightforward to find all symmetry-inequivalent minima of Eq. (12).There are two different states depending on the sign of b 2 which we label by A m v =0 ( s ; d ), where s and d refer to the singlet and the triplet vector, respectively, A indicates the trivial IR of C 3 , and M k,± = λ ±k σ 0 with λ C 3 k = λ k ; according to the notation introduced above, this state will be labeled as A m v =0 (1; 0, 0, 0).There are (infinitely) many different equivalent representations of this state since, for instance, the transformations in Eq. (11b) mix the singlet and triplet components-as described by the isomorphism SU(2) + × SU(2) − SO(4).However, for the sake of notational clarity, we will henceforth only show one convenient representative of each state.The A m v =0 (1; 0, 0, 0) state preserves time-reversal symmetry and breaks SU(2) + × SU(2) − down to SU(2) s [rotations of the total spin, i.e., ϕ + = ϕ − in Eq. (11b)].

B. Turning on Hund's coupling
In reality, there is, of course, a finite Hund's coupling that reduces G s 1 = SU(2) + × SU(2) − to only global spin rotations, G s 2 = SU(2) s , already in the high-temperature phase.In [57], Hund's coupling J has been estimated to be about 60 times smaller than the intravalley interaction V .Note, however, J might be enhanced due to loop corrections.For this reason, we first classify the possible instabilities in the absence of an approximate SU(2) + × SU(2) − symmetry and then analyze how the different states "connect" for small values of J and whether admixtures of singlet and triplet are possible.
To introduce our notation, we will begin with the classification for the reduced symmetry group G 2 in Eq. (3); in that case, we have either singlet or triplet pairing.

Singlet
This corresponds to the d n s = 1 one-dimensional IR of G s 2 with χ n s = σ 0 in Eq. ( 8).The pairing Hamiltonian simply has the form with λ s kv = λ s −k v and λ s C 3 kv = λ s kv .All symmetries of the high-temperature phase are preserved.We refer to this state as A 1 s m v =0 with the 1 s referring to spin singlet, m v = 0 to intervalley pairing, and A to the trivial representation of C 3 .

Triplet
This pairing channel is associated with the threedimensional IR of G s 2 .A possible choice of basis functions is χ n s μ (v) = σ μ , μ = 1, 2, 3, in Eq. ( 8).As it is a multidimensional representation, the free energy has to be expanded beyond quadratic order.Writing d = (η n s 1 , η n s 2 , η n s 3 ), we have up to quartic order Observe that |d T d| 2 is not an independent quartic term since 2 .The free energy in Eq. ( 14) has two stable minima.For b t 2 > 0, we have d ∝ (1, 0, 0) T and the corresponding pairing term is with As is easily seen, this term preserves time-reversal symmetry and breaks SU(2) s down to spin rotation along a single axis.This state will be referred to as unitary triplet and denoted by the symbol A 3 s m v =0 (1, 0, 0), where the three components just indicate the direction of the triplet vector.If, instead, b t 2 < 0, we obtain d ∝ (1, i, 0) T , whence (16) with λ t kv as above.This is a nonunitary triplet state.It breaks time-reversal symmetry and will be denoted by A 3 s m v =0 (1, i, 0).One might wonder what kind of interaction or band structure would favor (1, 0, 0) and vice versa.In mean-field theory, as detailed in Appendix A, it is straightforward to show by evaluation of a one-loop diagram that Here, ω n are fermionic Matsubara frequencies and ξ k+ is the electronic band energy in valley v = + of the nearly flat band hosting superconductivity.We observe that b t 2 > 0 holds irrespective of microscopic details and, hence, A 3 s m v =0 (1, 0, 0) is generically favored if we neglect corrections beyond meanfield theory (such as residual interactions or frequency dependence of pairing).Intriguingly, there have been experimental reports [69] of intrinsically nonunitary pairing in LaNiC 2 , i.e., nonunitary triplet pairing born out of a paramagnetic normal state.Thus, there is reason to believe that we cannot generically exclude this state, but we do not expect it to show up in any simple mean-field computation.

How do the states connect in the J = 0 limit?
Next, we establish how the three possible states, A 1 s m v =0 , A 3 s m v =0 (1, 0, 0), and A 3 s m v =0 (1, i, 0), connect to the two derived in the previous subsection with enhanced SU(2) + × SU(2) − symmetry, namely, A m v =0 (1; 0, 0, 0) and A m v =0 (1; 1, 0, 0).To this end, we decompose the Ginzburg-Landau expansion (12) into singlet and triplet by writing , singlet and triplet are degenerate at quadratic order in F as a consequence of the enhanced SU(2) + × SU(2) − symmetry.For nonzero J, this degeneracy is lifted and we have where δa can be made arbitrarily small as J → 0. Neglecting, for now, the "back action" of the superconducting order parameter that condenses first on the second one (as described by higher-order terms in the Ginzburg-Landau expansion), we conclude that there are two superconducting transitions at T ± c,0 = T c,0 ± T c with T c = |δa (T c,0 )|/α, taking a(T ) ∼ α(T − T c,0 ) near T c,0 .The extra index 0 in T ± c,0 highlights the fact that the aforementioned higher-order terms in the Ginzburg-Landau expansion can significantly affect the lower transition temperature, T − c = T − c,0 ; of course, this has no effect on the higher transition temperature, T + c = T + c,0 .Before analyzing these corrections, it is useful to estimate the temperature scale T c .Using the expected result, T ± c exp(−1/[(V ± J )ν]) of mean-field theory (from the linearized gap equations)-where is the cutoff and ν the density of states at the Fermi level-leads to The large density of states, taken together with the estimated value of V -which is larger than even the bandwidth [57] of the flat bands-and the relation J V implies that T c T c,0 [70]; the temperature/energy scale T c is most likely too small to be visible in experiments.While the estimate above is only based on mean-field theory, it indicates at least that it is important to study the behavior of superconductivity in the limit of small T c /T c,0 [and, hence, weakly broken SU(2) + × SU(2) − symmetry], accounting for the possibility of two transitions and mixing of singlet and triplet pairing (despite the absence of spin-orbit coupling).Moreover, we will see that nearly degenerate singlet and triplet pairing also has crucial consequences for the behavior of superconductivity in the presence of a magnetic field.
While we postpone the analysis of magnetic fields to Sec.III C, here, we investigate the possibility of an admixture of singlet and triplet in the presence of time-reversal symmetry (relevant to scenario I defined in the introduction).As anticipated above, this requires also considering the quartic terms of Eq. (12).We find (20) neglecting corrections to the quartic terms coming from finite J. Looking at the first transition with the higher transition temperature, we assess which of the two distinct triplet states, A 3 s m v =0 (1, 0, 0) and A 3 s m v =0 (1, i, 0), and the singlet state can be stabilized by starting from A m v =0 (1; 0, 0, 0) or A m v =0 (1; 1, 0, 0) and turning on a finite Hund's coupling J.For this purpose, we can neglect the coupling terms in the third line of Eq. (20).Clearly, if δa < 0 ("anti-Hund's coupling"), we get a singlet state for both A m v =0 (1; 0, 0, 0) and A m v =0 (1; 1, 0, 0).A straightforward way of establishing which of the triplet states is realized when δa > 0 ("conventional" Hund's coupling) proceeds by evaluating their respective free energy in Eq. (20).One finds that the state (1, i, 0) is favored.This brings us to the conclusion that at the first transition (see the schematic phase diagram in Fig. 2).This result is just a consequence of the fact that the form + ∝ σ 0 for the A m v =0 (1; 0, 0, 0) state we had chosen in the previous section can alternatively be written as + ∝ σ x due to the SU(2) + × SU(2) − symmetry and, thus, explicitly assumes the form of the unitary triplet state.Similarly, + ∝ σ 0 + σ z used above for A m v =0 (1; 1, 0, 0) can also be written as . This is why it transitions into the nonunitary triplet state, upon turning on a nonzero Hund's coupling.
In order to determine whether there is a second transition, we have to include the coupling terms between singlet and triplet in the third line of Eq. (20).To illustrate that these terms can be crucial, we consider the case δa > 0 and b 2 > 0, i.e., the triplet state A 3 s m v =0 (1, 0, 0) condenses first.This leads to the coupling between singlet and triplet 2c| in the free energy, where we have made use of the fact that a relative phase of π/2 between singlet and triplet is energetically most favorable.As a result of |d(T )| 2 = (δa − TABLE I. Summary of the different intervalley pairing states transforming under the trivial representation of the point group C 3 in the absence of a magnetic field.For notational convenience, we neglect the extra label m v = 0 to indicate intervalley pairing.λ k is a real-valued and Brillouin-zone-periodic function that is invariant under C 3 .To lowest order, we can take λ k to be independent of k.We also indicate the minimal number of nodes, which state it transforms to when setting J = 0 ["SO (4) parent"] and reversing the sign of J ("Hund's partner"), and whether the state can be found in a single-band mean-field (MF) computation neglecting residual interactions and/or when the ferromagnetic (FM) fluctuation corrections discussed in Sec.V dominate.In the last line, η describes the temperature-dependent strength of admixing of the unitary triplet state.

Pairing
M k+ Nodes SO(4) parent Hund's partner MF/FM a(T ))/(2c), which is valid as long as there is no additional singlet pairing, the growing triplet component induces the extra term which is always larger than the "bare" quadratic term of singlet pairing [in the first line of Eq. ( 20)].Accordingly, there is no second transition (at least close to T c,0 where our Ginzburg-Landau approach is valid) into a state that has a nonzero singlet component.We also checked that Eq. ( 20) does not allow for a first-order transition.
Similarly, all other cases can be scrutinized and one finds that if triplet dominates there is no second transition.However, if singlet has a larger transition temperature (δa < 0), there is a second transition into a phase with singlet and triplet pairing when b 2 < 0. This transition happens at the temperature The stability of the Ginzburg-Landau expansion only requires c > 0 and c > −b 2 , so both T − c < T − c,0 and T − c > T − c,0 are possible.More importantly, unless |b 2 |/c is fine tuned to be of order δ, generically, T − c → T c,0 as J → 0 and the two transitions, if present, are likely too close to be experimentally discernible.Due to the term 2b 2 Re[( s ) 2 d † d * ] in the free energy, we obtain the unitary triplet vector d = d 0 (1, 0, 0) T with s d * 0 ∈ R (same phase).This is to be expected as A summary of these results is provided by the schematic phase diagrams in Fig. 2. We observe that the proximity to the enlarged symmetry in spin space, SU(2) + × SU(2) − , favors the possibility of having a nonzero triplet component: for b 2 < 0, even a negative Hund's coupling (anti-Hund's) allows for d = 0 and leads to the exotic possibility of significant (d 0 s for T − c − T > T c ) singlet-triplet mixing in spite of the absence of spin-orbit coupling.
It is noteworthy that all the states are fully gapped (more precisely, they have no symmetry-enforced nodes) except for the nonunitary triplet A 3 s m v =0 (1, i, 0), which is gapped for one spin species while the other is completely gapless.The admixture of singlet and unitary triplet has two unequal gaps for the two spin species both of which are finite as long as the magnitudes of singlet and triplet are not fine tuned to be equal.All the states, along with their order parameters and properties, are summarized in Table I.
We finally comment on the nature of the thermal phase transition for the different superconducting states once fluctuations of the order parameter are taken into account.Neglecting stray fields, the transition into the singlet phase A 1 s is expected to be a Berezinskii-Kosterlitz-Thouless (BKT) transition with quasi-long-range order of the complex-valued order parameter s below the transition temperature.For the triplet states, it is important to keep in mind that d cannot even have quasi-long-range order as it transforms as a threecomponent vector under spin rotation.For the unitary triplet state [with order-parameter manifold (S 2 × S 1 )/Z 2 ] a BKT transition of the composite charge-4e order parameter d T d is possible and is associated with the (un)binding of half vortices.This is different for the nonunitary state [with order parameter manifold S 3 /Z 2 SO(3)] where d T d = 0 and no BKT transition into a quasi-long-range-ordered superconductor is expected.For the case of the two consecutive transitions in Fig. 2(b) with δ < 0, we first expect a BKT transition into a singlet phase followed by a crossover at which the triplet vector becomes nonzero.
However, we point out that, even in the simplest case of the singlet A 1 s , there are significant corrections to the BKT transition resulting from stray fields and mirror vortices [71], which make the observation of a pristine BKT transition in a (charged) superconductor difficult.We believe that the current status of experiments does not allow one to exclude pairing phases that will not exhibit quasi-long-range order and a BKT transition in the limit of infinite system size.

Expectations within mean-field theory
Lastly, we evaluate what a naïve mean-field computation is expected to yield.In fact, from Eq. ( 17), we already know that the prefactor of the term |d * × d| 2 in Eq. ( 20) must be positive within mean-field theory and, therefore, it holds that b 2 > 0. For completeness, we mention that in the meanfield approximation b 1 = 0, as shown in Appendix A. Consequently, a single-band mean-field computation will generally favor Fig. 2(a) over Fig. 2(b); in other words, only half of the phases proposed in this section can be found in mean field, which we also indicate in the last column of Table I.
However, there is no fundamental mechanism prohibiting the mixing of singlet and triplet via two transitions (see, e.g., [72]) and there are multiple reasons why we can effectively have b 2 < 0 (and b 1 > 0 to ensure stability): for instance, strong residual interactions and fluctuations have been shown to modify the values of the quartic terms in the free energy significantly [43,73], thereby stabilizing phases that are otherwise not possible in the mean-field approximation.Given the small bandwidth and the underlying strong-coupling features of the problem [5][6][7][8][9], it is plausible that there are sizable corrections to mean-field theory.In addition, we recognize that there are other corrections arising from interband pairing, and that disorder can also dress the Ginzburg-Landau expansion.Moreover, it is unclear whether adding frequency dependence to the gap function could be of relevance.
In Sec.V, we will analyze the impact of ferromagnetic fluctuations, which are expected to be relevant for graphene moiré systems [4,[51][52][53][61][62][63], and find that these generically decrease the value of b 2 ; if sufficiently strong, these fluctuations will favor the phase diagram in Fig. 2(b).

C. In the presence of a magnetic field
We now generalize the Ginzburg-Landau expansion to also include the coupling to a Zeeman field Both of these terms can either be due to an applied external magnetic field or due to the correlated insulating state.This enables us to discuss (i) the behavior of the superconducting critical temperature T + c as a function of an external magnetic field in the absence of any ferromagnetic moments associated with the correlated insulating state [case I defined in the introduction].At the same time, we can study (ii) how the transition temperature and the order parameter of superconductivity is affected by the potentially coexisting ferromagnetic order [case II].

Leading superconducting transition
We first turn our attention to the leading superconducting transition with the highest temperature T + c ; potential subsequent superconducting transitions at lower temperatures are addressed later in Sec.III C 2. For the goal of studying the first transition, we can restrict ourselves to quadratic order in the order parameter.Only keeping terms up to quadratic order in the magnetic field as well, we obtain While the prefactors δa, δc 1 , δc 5 , and δc 6 are necessarily zero in the limit J → 0, where the SU(2) + × SU(2) − symmetry becomes exact, all remaining terms can be nonzero (and different in their values) at J = 0. Notice that the third term has not been considered in [57]; this term arises only when both singlet and triplet are allowed for and leads to the admixture of a unitary triplet state with a singlet superconductor.The vanishing of δa and δc 6 at J = 0 is an obvious consequence of the enhanced SU(2) + × SU(2) − symmetry.To see that δc 1 also has to vanish as J → 0, let us take M Z along the z direction; this breaks SU(2) O(2) − , i.e., the system is only invariant under c kv → e iϕ v σ z c kv .Performing this transformation with ϕ + = 0 and ϕ − = π/2, we get ( s , d z ) → (id z , i s ) and, hence, δc 1 → −δc 1 .With the same argument, it can be proven that δc 5 has to go to zero as J → 0. In Appendix A, we show that δc 1 = 0 in meanfield theory within the single-band description, even when SU(2) + × SU(2) − is broken; this results from an emergent valley-exchange symmetry within the single-band mean-field approximation.
In discussing the highest critical temperature and the corresponding order parameter for M Z , M O = 0, it is instructive to first look at the linear-in-field terms in Eq. (24).We find two different cases.If with M Z > 0, the triplet vector is given by d = (0, 1, sign(c 2 )i) T and the critical temperature is Else, if + δa 2 , one finds an admixture of singlet and triplet with order parameter The transition temperature in this case is We see from Eq. ( 26) that there is an approximately equal mixing of singlet and triplet for |δc 1 |M Z |δa| while in the opposite limit, |δc 1 |M Z |δa|, either singlet or triplet dominates depending on whether δa < 0 or δa > 0. The relative phase of π/2 between singlet and triplet makes the pairing state break time-reversal symmetry as is required in order to couple linearly to magnetic moments.
To understand how the approximate SU(2) + × SU(2) − symmetry can naturally explain the linear-in-magnetic-field behavior, we first consider case II, i.e., there is already microscopically coexisting ferromagnetic order (or there is at least a significant coupling between superconductivity and the ferromagnetic moments) at T + c .Then, M Z and M O should be thought of as the combination of the applied external magnetic field and the ferromagnetic order parameter.In this scenario, it is apt to assume |δa| max(δc 1 M Z , c 2 M Z ) and we generically obtain a linear increase of the critical temperature with magnetic field [see Eqs. ( 25) and ( 27)].If c 2 > δc 1 , we obtain the nonunitary triplet state with d ∝ (1, i, 0) T , which we expect close to the J = 0 line, while δc 1 > c 2 leads to the admixture of singlet and triplet with d ∝ (1, 0, 0) T .As δc 1 vanishes for J = 0 and in single-band mean-field theory (even when J = 0), we expect the former scenario to be more likely, which will favor the nonunitary triplet state as the leading instability.
In the case of scenario I, we should view M Z and M O in Eq. ( 24) as resulting entirely from the Zeeman and orbital coupling of the external magnetic field alone.For large magnetic fields where |δa| max (δc 1 M Z , c 2 M Z ), the same conclusions as above will apply and T c will generically vary linearly with the field.However, for sufficiently small magnetic fields, we have |δa| max(δc 1 M Z , c 2 M Z ).In this limit, only δa > 0 favoring the nonunitary triplet pairing is consistent with the transition temperature changing linearly with magnetic field.Alternatively, the system could ultimately be in a singlet state at M Z = 0 (i.e., δa < 0) but the magnitude of δa is sufficiently small such that the "rounding off" of T + c (M Z ) at low M Z cannot be seen in experiment.

Quartic terms and subleading transitions
Having examined the first superconducting transition that takes place upon cooling the system down starting from the normal state, we now assess whether and what type of subsequent superconducting transitions can occur.In this context, we need to include terms quartic in the superconducting order parameter and extend Eq. ( 24) to where we kept only the terms linear in magnetic field, took M Z along the z axis, and reexpressed the triplet in the ).This parametrization is more convenient in the presence of a magnetic field than that used in Eq. (20).Additionally, we have neglected the impact of the magnetic field on the quartic terms.
Taking δc 1 = 0 (as it has to vanish for J = 0), the different possible phase diagrams are summarized in Fig. 3.The possibility illustrated in part (c) of Fig. 3 corresponds to the picture put forward by Ambegaokar and Mermin [74] for He 3 in the presence of a magnetic field, which might very well also apply to twisted double bilayer graphene [52,57].The difference with [74] is that we do not get a third transition since we work with a one-dimensional IR of the spatial point group.
However, there are three other options, depicted in Figs.3(a), 3(b), and 3(d), that we cannot easily exclude given the experimental data: owing to the strong-coupling properties of the problem at hand, a nonunitary triplet state might be dominant at M Z = 0, as seems to be the case in LaNiC 2 [69] and is favored by our fluctuation approach of Sec.V; under this condition, only one transition is expected even when M Z = 0 [see Fig. 3(d)].It could also be that singlet dominates without a magnetic field instead.We can see in Figs.3(a) and 3(b) that, in these two cases, triplet shows up and T c increases linearly when M Z > 2|δa|/c 2 .The small value of T c /T c,0 estimated in Eq. (19) suggests that resolving this initial region, where T c is constant as a function of M Z , is experimentally challenging.28).Thin (thick) black lines correspond to second (first) order transitions.The phases for M Z = 0 are indicated in red and we recover the four different possible temperature dependences of Fig. 2. Recall from Sec. III B 4 that b 2 > 0 and, thus, (a) or (c) is expected in mean-field theory.However, as we will see in Sec.V, strong ferromagnetic fluctuations will favor b 2 < 0 and, hence, (b) or (d).As symmetry requires δc 1 to be proportional to Hund's coupling J, we have set δc 1 = 0 here.For nonzero δc 1 , the singlet superconducting phases will contain an admixture of unitary triplet as described by Eq. ( 26) and a first-order transition into a singlet state (with unitary triplet admixture) will be possible at lower temperatures and nonzero Zeeman field in part (c).Note that the transition temperature from the normal state into the singlet superconductor is constant as we neglect here the nonlinear coupling to the magnetic field.A discussion of the latter can be found in Sec.III C 3.

Nonlinear couplings in a magnetic field
We finally come back to the quadratic couplings to the magnetic field, associated with the terms with prefactors c 3,4 and δc 5,6 in Eq. (24).We first notice that this will lead to an additional quadratic suppression of the leading transition temperatures in Fig. 3; in particular, the transition temperature into the singlet state in part (a) and (b) will not be field independent any more.More interestingly, the suppression of singlet and triplet is enforced to be nearly identical for small J due to the SU(2) + × SU(2) − symmetry, δc 5,6 c 3,4 .Resultantly, if the effective J relevant for superconductivity is indeed small, the nonlinear terms ∝ M 2  Z , M 2 O are not expected to affect the competition between singlet and triplet significantly and the qualitative form of the phase diagrams in Fig. 3 is not modified.

IV. COMPLEX REPRESENTATION OF C 3
In this section, we extend our previous analysis to the complex IR E of the spatial point group C 3 .Time-reversal symmetry necessitates treating the representation and its complexconjugate partner on an equal footing.Alternatively, one can think of a two-dimensional (2D) (reducible) representation with partner functions transforming as x and y under C 3 .
Akin to our discussion earlier, we first study the case of nonzero Hund's coupling, J = 0, with point group G 2 in Eq. ( 3), which enables us to distinguish between singlet and triplet pairing.After discussing all symmetry-allowed singlet and triplet states separately, we will derive the phase diagrams analogous to Fig. 2: we will examine how these states "connect" when adiabatically changing Hund's coupling from negative to positive values, and whether singlet and triplet can mix when J is small and the SU(2) + × SU(2) − symmetry is only weakly broken.

A. Nonzero Hund's coupling
To proceed with singlet pairing, we parametrize M kv in Eq. ( 9) according to while M k− is determined by the Fermi-Dirac constraint (10); X k and Y k are real-valued functions that are continuous on the Brillouin zone and transform as k x and k y under C 3 .A one-parameter family of possible choices for the lowest-order functions (i.e., with minimal number of sign changes in the Brillouin zone) is given by with arbitrary φ ∈ [0, 2π ), where R φ is a 2 × 2 matrix describing rotations by angle φ, R φ = e iφσ y , and Both X k and Y k have to vanish at , K, and K as these momenta are invariant under C 3 .Further, both X k and Y k must have lines of zeros going through these high-symmetry points; the orientation of these lines is, however, not fixed due to the absence of additional reflection or in-plane rotation symmetries-this is different from the situation for twisted bilayer and trilayer graphene in Sec.VI.For Eq. ( 30), the orientation of these zeros changes with φ.
With the parametrization defined in Eq. ( 29), the relevant symmetries act as follows: It readily follows from Eq. ( 31) that the most general free energy up to quartic order reads as The sign of b s 2 therefore distinguishes between two different singlet phases: if b s 2 > 0, we have (η + , η − ) = (1, 0), which corresponds to Exactly as in Sec.III, we always show only one out of the many symmetry-equivalent representations of the order parameter-instead of using a general parametrization of a phase-to make the notation and the discussion of properties of the superconducting state more easily accessible.The state in Eq. ( 33) breaks time-reversal symmetry but preserves C 3 (and spin-rotation symmetry).We refer to this state as a chiral singlet superconductor and denote it by E 1s (1, i) in the following.It is fully gapped (unless the Fermi surfaces go through the , K, or K point) and has been investigated extensively in the recent literature on pairing in twisted bilayer graphene [12,21,27,30,33,36,37,40,42,44].
Conversely, if b s 2 < 0, we find that |η + | = |η − | at the minimum of Eq. (32).As the relative phase ϕ between η + and η − = η + e iϕ is not fixed by Eq. ( 32), one might naively conclude that higher-order terms have to be considered.In fact, in sixth order, there is indeed the contribution and the relative phase ϕ will depend on c 1 /c 2 .However, upon reinserting η − = η + e iϕ into Eq.( 29), we notice that ϕ = 0 simply corresponds to rotating the basis functions X k and Y k into each other, which does not change their transformation behavior under C 3 [ϕ is directly related to φ in Eq. (30a)].Consequently, we can set ϕ = 0 without loss of generality, which implies This state, which we call E 1s (1, 0), breaks C 3 but preserves time-reversal symmetry; this is the nematic singlet phase.Within a single-band mean-field description (see Appendix A), we find b s 1 = b s 2 /2 > 0. As such, mean-field theory generically favors the chiral singlet superconductor over the nematic state E 1s (1, 0); this has been noted before in the context of twisted bilayer graphene [44] and [43] discusses how strong fluctuations can stabilize the nematic phase.
Turning to triplet pairing, we now modify the parametrization (29) to where X k and Y k are defined exactly as before.For simplicity, we introduce the complex-vector notation, d μ = (η μ,1 , η μ,2 , η μ,3 ) T , μ = ±.The representations of the symmetries now read as with R ∈ SO(3) and ω = e i 2π 3 .The most general free-energy expansion is given by up to quartic order, where b t j ∈ R; the different symmetryallowed phases follow from the stable minima of the free energy.When minimizing Eq. ( 38), we take into account that the relative phase between d + and d − can always be absorbed into a redefinition of the basis functions X k and Y k , as for the singlet above.In total, we find eight distinct triplet states which we label by E 3 s (a) through E 3 s (h).Phase diagrams describing which of these phases is realized for a given configuration of the quartic couplings b t j can be found in Appendix C; here, we list all the phases, describe their properties, and refer to Fig. 4 for an illustration of their respective spectra and densities of states.
(a) This state, labeled as E 3 s (a), can be represented by d + = d − = (1, 0, 0) T with the associated order parameter M k+ = X k σ x .More physically, it corresponds to a nematic unitary triplet phase.It preserves time-reversal symmetry, but breaks both SU(2) s spin-rotation symmetry [down to O(2)] and C 3 rotational symmetry.This state has two symmetryenforced nodal points at each Fermi surface around the K, K , or point.Owing to the lack of any reflection symmetry (see the discussion of D 3 in Sec.VI below), the positions of these nodal points are not pinned to any specific direction.
(b) One representative configuration of this phase is given by d + = (1, −i, 0) T /2 and d − = (1, i, 0) T /2; it can thus be seen as a helical triplet, consisting of two time-reversed copies of states with opposite chirality.The order parameter can be more explicitly written as M k+ = X k σ x + Y k σ y , which can alternatively be thought of as a 2D analog of the Balian-Werthamer state of the B phase of superfluid 3 He [75].This state, denoted by E 3 s (b) in the following, only has point nodes at , K, and K , i.e., it is expected to exhibit a full gap for generic Fermi surfaces not going through these highsymmetry points.It preserves time-reversal symmetry.While this state breaks spin-rotation symmetry as well as C 3 , the product of C 3 and a rotation in spin space along σ z with angle 2π/3 is preserved; this can be viewed as the spontaneous formation of spin-orbit coupling.
(c) Here, we can write d + = d − = (1, i, 0) T /2; hence, M k+ = X k (σ x + iσ y ).This is a nematic nonunitary triplet state which breaks time-reversal symmetry and C 3 .One spin species will be gapless while the other will have nodal lines (i.e., point nodes on the Fermi surface).
(d) The triplet vectors in this phase can be written as The Bogoliubov-de Gennes excitation spectrum is calculated using the band structure of the system predicted by the continuum model [58], assuming a pairing term of the form of Eq. ( 9) with an overall scale of 0 = 4 meV.The nodal points/lines are demarcated in dark blue; note that the states (b), (d), (g), and (h), taking α = π/4, are fully gapped; (a) and (f) have nodal points; and (c) and (e) have nodal lines, as is also visible in g(E ).identical to that of the helical triplet E 3 s (b), which is why we group these two states together in Fig. 4.
(e) For this state, we have It consists of only one of the two time-reversed copies with opposite chirality of the E 3 s (b) state discussed above and, thus, is a chiral nonunitary triplet state.This state can be seen as an analog of the A 1 phase of 3 He [75].It preserves C 3 , but breaks SU(2) s spin-rotation symmetry [down to O(2)] and time reversal.Here, one of the spin components will be gapless while the other is fully gapped (as before, except for the high-symmetry points , K, and K which are generically not on the Fermi surface).Note that although the spectrum of this state is not strictly identical to that of the nematic nonunitary triplet E 3 s s (c) we grouped them together in Fig. 4 as their respective plots are practically indistinguishable; this is related to the fact that, in both cases, the low-energy spectrum is dominated by the Fermi surface of one of the spin species.
(f) In this phase, The state can, thus, be thought of as a superposition of two chiral unitary triplets with orthogonal spin polarizations or, when inserted into Eq.( 9), as Cooper pairs of electrons with spin polarization ↓↓ (↑↑) and orbital basis function and spin-rotation symmetry are all broken.The excitation spectrum is given by 2 , so it is characterized by "two gaps," given by |X k ± Y k |, both of which are forced to vanish at two points for each Fermi surface enclosing K, K , and .While the number of nodes of this state and of E 3 s (a) are the same, the spin degrees of freedom on the Fermi surface have nodes at the same two momenta for E 3 s (a).For E 3 s ( f ), however, the two spin species have nodal points at different momenta.
(g) Denoted by T , where the parameter α varies continuously with b t j in the part of the phase diagram where this state is realized.The corresponding order parameter can be written as and can be viewed as a superposition of a chiral nonunitary triplet state and a unitary state with opposite chirality.This state breaks time-reversal symmetry, spin-rotation invariance, and C 3 but preserves the product of C 3 and spin rotation by angle 2π/3 along σ z .So, similar to the state E 3 s (b) above, this state spontaneously entangles rotations in spin and real space and its spectrum [see Fig. 4(g)] is C 3 invariant.It is fully gapped (again, as long as the Fermi surfaces do not go through , K, and K ), with two different gaps , where g α = cos α 1 + sin 2 α.(h) Finally, for the triplet phase E 3 s (h), one has d + = (cos α, 0, i sin α) T , d − = (0, cos α, −i sin α) T , which yields It can be seen as a superposition of the states E 3 s (a) and E 3 s ( f ) to which it reduces for α = π/2 and 0; it will have two nodal points for α close to these limiting cases, but can be fully gapped for other values of α.For α = π/2, this state breaks time-reversal, C 3 , and spin-rotation symmetry.
In Appendix A, we show that b t 5 > 0 and b t 2 = 0 within a single-band mean-field description.Minimizing Eq. ( 38) yields that the phases E 3 s (b) and E 3 s (d ) have the lowest energy and are exactly degenerate for this configuration of quartic couplings.This degeneracy within mean-field theory, which was noted before in [19], will be lifted by corrections resulting, e.g., from residual interactions.In Sec.V, we will find that E 3 s (b) (E 3 s (d )) is favored in the presence of ferromagnetic spin (orbital) fluctuations.We will also see that significant fluctuations can stabilize phases other than the two, E 3 s (b) and E 3 s (d ), favored in mean-field theory.

B. Approximate SU(2) + × SU(2) −
After having classified singlet and triplet separately, we now focus on small Hund's coupling for which SU(2) + × SU(2) − is an approximate symmetry, and singlet and triplet are nearly degenerate at the quadratic level of the free energy.This requires studying them on an equal footing and generalizing the parametrization in Eqs. ( 29) and (36) to include both singlet and triplet, i.e., extending the summation over ν in Eq. ( 36) to ν = 0, 1, 2, 3.In analogy with Sec.III A, we use 2 × 2 matrices and write It is easy to see that the symmetries act according to where, recall, G s 1 ≡ SU(2) + × SU(2) − .Imposing SU(2) + × SU(2) − as an exact symmetry, the most general free energy up to quartic order reads as At first glance, one might think that there are additional terms with extra factors of σ y , similar to the last term in Eq. ( 12).However, as before, all of them can be related to combinations of the terms already present in Eq. ( 41) as outlined in Appendix C. Following the procedure applied in Sec.III to the onedimensional IR A, we now add a small quadratic term, , where s μ and d μ are the singlet and triplet component of μ in Eq. ( 41), i.e., μ = σ 0 s μ + σ • d μ .This term breaks SU(2) + × SU(2) − and, hence, makes singlet and triplet inequivalent.It allows us to study which of the different singlet and triplet states defined above can mix, and to identify "Hund's partners," i.e., which states transform into each other when changing the sign of Hund's coupling J and, accordingly, of δa.This generalizes the phase diagrams in Fig. 2 and Table I to the complex representation.
We find that, out of the eight different triplet states E 3 s (a) to E 3 s (h), only two-E 3 s (a) and E 3 s (d )-do not allow for a singlet-triplet admixture when reversing the sign of J (or δa) so that singlet has the higher transition temperature.The reason for the absence of an admixture is the same as sketched by way of example in Sec.III B: besides pure singlet and pure triplet terms, the quartic terms in Eq. ( 41) also contain couplings between singlet and triplet, as is readily seen by inserting the parametrization μ = σ 0 (the full expansion can be found in Appendix C).At the first transition, one of either singlet or triplet becomes nonzero and, hence, "renormalizes" the quadratic term of the other channel.In some cases, this renormalization can prohibit the presence of a second transition.In the case of phases E 3 s (a) and E 3 s (d ), we just obtain the pure singlets E 1 s (1, 0) and E 1 s (1, i), respectively, without a second transition.The easiest way to interpret why we do not have an admixture in these cases is to look at the associated SO( 4 For all other triplets, Hund's partner is an admixed phase.Specifically, as regards E 3 s (b) and E 3 s (c), Hund's partner is an admixture of a nematic singlet state and a nematic unitary triplet E 3 s (a), with different relative phases and spatial orientations: for the former, the order parameter can be written as i Y k σ 0 + ηX k σ x , where η describes the temperaturedependent strength of mixing, while it is X k (σ 0 + ησ x ) for the latter.On any Fermi surface around one of the highsymmetry points , K, or K , these two states have zero and two nodal points, respectively.Again, the form of the admixed state can be understood from the representation of the triplet state in terms of (η μ ; d μ ).For instance, we have (η Likewise, Hund's partners of E 3 s (e) and E 3 s ( f ) are admixtures of a chiral singlet and a unitary triplet state with the same and opposite chirality, respectively.The associated order parameters can be written as ( While the first of the two states has two fully established gaps, given by (1 ± η) √ X 2 k + Y 2 k (with ± referring to the spin species), the other has two gaps, |X k | and |Y k |, with distinct momentum dependencies; it, thus, exhibits two point nodes per Fermi surface which occur at different positions for the two spin species, similar to the associated triplet phase E 3 s ( f ).
In general, admixing a singlet component at a second transition to a triplet state is less likely to occur as a singlet state has less options to "adapt" (the order parameter comprises two complex numbers for E ) than a triplet state (for which the order parameter comprises six complex numbers).While this is not possible for the one-dimensional representation A (see Fig. 2), the IR E does allow for this scenario but only for the triplet states E 3 s (g) and E 3 s (h): for small δa < 0, we find a second transition where an additional chiral (nematic) singlet component is admixed to E 3 s (g) (E 3 s (h)).As both pure triplet states can be fully gapped, the same holds for the admixed phases.The admixture of the extra singlet component does not change the symmetries of E 3 s (g) and E 3 s (h) listed in Sec.IV A above.Reversing the sign of δa to small positive values, we obtain the same admixed phase.The only difference is that the first transition is a singlet transition into a chiral (nematic) phase and the secondary triplet E 3 s (g) [E 3 s (h)] becomes nonzero at a lower transition temperature.The key results of this section, the pure triplet/singlet states and the possible admixed phases for small J along with their order parameters and properties, are summarized in Table II.As already discussed above, several states are degenerate within single-band mean-field theory.Depending on the form of the corrections to mean-field theory lifting this degeneracy, there are two possible phase diagrams, shown in Fig. 5. Interestingly, we observe that the chiral singlet, E 1 s (1, i), is not the only possible phase close to mean-field theory for anti-Hund's coupling: as can be seen in Fig. 5(b), a secondary phase transition into the nematic mixed singlettriplet state E 1 s (0, i) + E 3 s (a) is predicted.It is a fully gapped state with an anisotropic gap, √ η 2 X 2 k + Y 2 k , breaking rotational symmetry.Note that this route to a nematic superconducting state, indications of which are provided by recent experiments [76], is distinct from that of other works [43,77].Of course, sufficiently large corrections to mean-field theory can in principle yield any of the phases listed in Table II; we will come back to these corrections in Sec.V below.
Let us finally discuss the impact of fluctuations of the order parameter on the thermal phase transitions.As readily follows from the respective order-parameter manifolds, the singlet phases in Table II exhibit a conventional BKT transition, the triplets (a), (b), (d), (f), (g), and (h) will be charge-4e superconductors where only spin-rotation invariant combinations of the triplet vector assume quasi-long-range order at finite temperature, and the triplets (c) and (e) will only display a crossover.However, as pointed out above, none of these three classes of transitions can currently be excluded based on the experimental data.

C. Behavior in a magnetic field
Finally, we turn our attention to the behavior of the pairing states of the complex representation in the presence of a Zeeman field, M Z , and in-plane orbital coupling M O , along the same lines as Sec.III C. From Eqs. ( 31) and (37), it follows that there are three possible coupling terms linear in the field and quadratic in the superconducting order parameter given TABLE II.Summary of possible pairing states transforming under the complex representation E of C 3 .The labeling of the pairing states and their symmetry properties can be found in the main text.The states are ordered by pure singlet, triplet, and admixtures of singlet and triplet.The latter are only expected generically when the SU(2) − × SU(2) + symmetry is weakly broken.We use X k and Y k to denote real-valued continuous functions on the Brillouin zone that transform as k x and k y under C 3 [see, e.g., Eq. ( 30)].The temperature-dependent coefficient η describes the admixture of a triplet/singlet pairing at a second transition to a purely singlet/triplet one.Furthermore, a, b ∈ R vary continuously with system parameters.The minimal number of nodes on any Fermi surface enclosing the , K, or K point is indicated in the column "Nodes."As before, two states are referred to as Hund's partners if they transform into each other under reversing the sign of Hund's coupling (see, e.g., Fig. 5).As singlet and triplet mix for both δ > 0 and δ < 0, there are no Hund's partners for E 3s (g) and E 3s (h); the corresponding mixed phases, contained in the last two lines of the table, are their own Hund's partners.

Pairing
M k+ Nodes Hund's partner MF/FM by Notice that, exactly as for the IR A, there is no linear coupling to the in-plane orbital field, which is prohibited by timereversal and C 3 rotation symmetry.While the first term in Eq. ( 42) is again forced to vanish for J → 0 [for the same reason as δc 1 in Eq. ( 24)], the second singlet-triplet-mixing coupling, c E 2 , is not constrained to be zero for J = 0.However, the emergent symmetry in the single-band mean-field description of Appendix A leads to c E 2 = 0, so it is natural to expect c E

c E
3 such that the last term in Eq. ( 42) describes the dominant linear coupling to the magnetic field-even when J is small.As expounded in Appendix A, the expression for c E 3 is identical in form to that for c 2 in Eq. ( 24).As such, the linear increase of the (first) superconducting transition temperature with small magnetic fields seen in experiment does not permit one to distinguish between the IRs A and E .
There is one difference between the pairing states of the two IRs worth mentioning here: while the form of the leading triplet vector in a magnetic field is completely fixed to be d ∝ (1, i, 0) T for the one-dimensional IR A, the complex IR allows for either the nematic nonunitary E 3 s (c) or the chiral nonunitary E 3 s (e) pairing for nonzero M Z .Which of the two is realized depends on the value of the quartic terms in Eq. ( 38): if b t 2 + b t 3 > 0, the state E 3 s (e) will be preferred while the opposite sign corresponds to E 3 s (c).Within single-band mean-field theory, we find b t 2 = 0 and b t 3 > 0, which leads to phase E 3 s (e).In the next section, we will see that additional ferromagnetic fluctuations will further enlarge the positive value of b t 2 + b t 3 and, consequently, not affect the mean-field prediction that E 3 s (e) is the leading triplet state with the highest transition temperature in the presence of a magnetic field.

V. FLUCTUATION-INDUCED SUPERCONDUCTIVITY
Among the plethora of possible superconducting phases outlined in this paper, only a few can be realized in singleband mean-field theory (see Tables I, II, and IV).This originates from the fact that, within single-band mean-field theory, the ratio of the quartic terms is fixed and only one state or two degenerate states can occur for each IR.However, the presence of sizable correlations in the nearly flat bands of graphene moiré systems is expected to give rise to significant corrections to mean-field theory.This has recently been demonstrated for the case of charge-density-wave fluctuations in twisted bilayer graphene [43], and in the context of nematic fluctuations in the iron-based superconductors [73].In this section, we study how corrections associated with ferromagnetic fluctuations will split the mean-field degeneracies and, if sufficiently strong, realize phases distinct from mean-field theory.
To this end, we will first focus on spin fluctuations.This is prompted by experiments [51][52][53], which indicate a spinpolarized correlated insulating state in twisted double bilayer graphene, and by the fact that the superconducting phase emerges when doping out of this polarized state.Likewise, we also expect ferromagnetic fluctuations to play an important role in twisted bilayer [4,61,62] and trilayer graphene [63].In particular, in the latter two systems, however, these fluctuations will likely not only be of spin but also of orbital origin.This is why we will also discuss orbital fluctuations.
As it is known to capture the essential physics [43,73], we focus in the main text on a phenomenological Ginzburg-Landau-like approach (that does not explicitly take into account fluctuations with nonzero momentum and frequencies), but provide a systematic microscopic derivation in Appendix B1.Representing the ferromagnetic spin moment in valley v = ± by m v , we parametrize its contribution to the free energy as In this expression, χ plays the role of the spin susceptibility (with |δχ| < χ to ensure stability) and we expect δχ > 0 close to a phase where the spin moments in the two valleys are aligned.The ratio δχ/χ controls how strongly the SU(2) + × SU(2) − symmetry is broken down to SU(2) s .

A. Trivial representation
Focusing first on the one-dimensional IR A of C 3 , the magnetic moments couple to the superconducting order parameter in Sec.III according to where we have retained only the couplings invariant under SU(2) + × SU(2) − and assumed that δχ = 0 in Eq. ( 43) is the main symmetry-breaking perturbation.Upon making the association M Z = v m v , we notice that c 2 is the same prefactor as in Eq. ( 24).In the same vein as [43], we integrate out the massive fluctuations of m v .As a consequence of the coupling (44), this yields corrections to the terms quartic in the superconducting order parameters in Eq. ( 20), which can be conveniently split into two categories.First, there are corrections that preserve the SU(2) + × SU(2) − symmetry; these can be restated as renormalizations of the coefficients b 1 and b 2 in Eq. (20).Corrections of the second type break this symmetry, violating the form of the free-energy expansion (20).More explicitly, the renormalization of the free energy F in Eq. ( 20) due to the presence of ferromagnetic spin fluctuations can be compactly stated as where δ 1 = −δ 2 = 2c 2 2 (χ − δχ ) > 0 and δ 3 = 2c 2 2 δχ.As required by symmetry, the contribution δ 3 of the second category breaking the SU(2) + × SU(2) − symmetry is proportional to δχ.
We start with the limit |δχ| χ , where the structure of Eq. ( 20) is asymptotically preserved and the form of the two possible phase diagrams in Fig. 2 is unchanged.Since δ 2 < 0, strong ferromagnetic fluctuations will change the sign of b 2 from its positive mean-field value to negative and, as opposed to mean-field theory, favor the phase diagram in part (b) of Fig. 2 over part (a).We point out that naively taking Eq. ( 45) alone would render the quartic free-energy expansion unstable for large enough χ .However, denoting the mean-field value of b 2 by b 0 2 , there exists a regime, b 0 2 /2 < c 2 s χ < b 0 2 , for which b 2 < 0 due to fluctuation corrections and the free energy in Eq. ( 20) is stable.For larger values of χ , we can imagine adding the sextic term c(tr[ † + + ]) 3 to the free energy to restore stability.
When δχ is of order χ , the ferromagnetic fluctuations described by Eq. ( 43) induce considerable SU(2) + × SU(2) − symmetry-breaking interactions.The presumed sign δχ > 0 brings about a further enhancement of the term −|d * × d| 2  [as is obvious from Eq. ( 45)], which favors nonunitary triplet pairing relative to the SU(2) + × SU(2) − invariant form of the free energy in Eq. (20).Given that δχ < χ, strong ferromagnetic fluctuations are still expected to change the sign of b 2 relative to mean-field theory.The additional effect of δχ lies in effecting an additional first-order transition to a nonunitary triplet state in a third transition at lower temperatures for anti-Hund's coupling in Fig. 2(b).
We have thus shown that significant ferromagnetic fluctuations can reverse the predictions of mean-field theory, and favor the nonunitary triplet state A 3 s (1, i, 0) and the admixed singlet-triplet phase A 1 s + A 3 s (1, 0, 0) in Table I.

B. Complex representation
The same analysis can be performed for the complex IR E of Sec.IV.In this case, the most general SU(2) + × SU(2) − invariant coupling between the superconducting order parameter and the spin fluctuations allows for two independent coupling constants, c ± ∈ R, and has the form Integrating out m v , we again obtain corrections to the free energy which are quartic in the superconducting order parameter.In the limit of SU(2) + × SU(2) − invariance, δχ = 0, these corrections can be represented by renormalizations of the couplings, b j → b j + δb j , in Eq. ( 41) with To study the ramifications of this result, we first consider the limit of weak fluctuations, for which δb j in Eq. ( 47) are much smaller in magnitude than the mean-field value of b 2 .Albeit small, the corrections δb j are crucial here due to the exact degeneracy of the states E 3 s (b) and E 3 s (d ) in mean-field theory observed earlier.From Eq. ( 41) with the replacement b j → b j + δb j , we find the free-energy difference of these two states to be thereby generically favoring E 3 s (b) along with its Hund's partner E 1 s (0, i) + E 3 s (a), defined in Table II; in other words, the phase diagram in Fig. 5(b) is favored over that in part (a).
In the one-band description of Appendix B2, it always holds that c + = c − , which is, in turn, a consequence of an emergent valley-exchange symmetry.However, multiband effects are expected to be present [8] and to lead to nonzero |c + − c − | |c + |, which is enough to lift the degeneracy according to Eq. ( 48).Next, we turn to the limit of strong ferromagnetic fluctuations, where the mean-field values of b j have to be treated as perturbations to the large δb j in Eq. (47).As χ → ∞, we find that, out of the triplet states in Table II, E 3 s (e) has the lowest energy unless c + = c − or c + = −c − .We know that c + c − and, hence, can safely neglect the latter.For the former option, E 3 s (e) is found to be degenerate with E 3 s (c); however, for large but finite χ , the additional contribution to b j from mean-field theory lifts this degeneracy, always selecting E 3 s (e).Out of the multitude of possible pairing states in Table II, strong ferromagnetic fluctuations thus favor the chiral nonunitary triplet state E 3 s (e) and the mixed singlettriplet phase E 1 s (1, i) + E 3 s (d ).Which of these two states is realized depends on whether singlet or triplet has the higher transition temperature (the sign of δa).

C. Orbital fluctuations
Anticipating its relevance for twisted bilayer and trilayer graphene, here, we extend the previous analysis to the case of orbital ferromagnetic fluctuations.Due to the two-dimensional nature of the system, the in-plane orbital moments, , and the out-of-plane moment M z O behave quite differently.Beginning with the complex representation, we already know from Sec. IV C that there is no linear coupling to M O ; however, the superconductor can couple to M z O as For concreteness, one might think of M z O as valley fluctuations, associated with k c † k τ z c k , but our analysis is more general.Taking an energetic contribution quadratic in M z O similar to Eq. ( 43) and integrating over M z O , we obtain a correction to the free energy that can be conveniently expressed as in Eq. (41).It is easily seen that taking this as a small correction to mean-field theory will now favor the phase diagram in Fig. 5(a) over that in part (b).On the other hand, in the limit of strong orbital fluctuations, the chiral unitary, E 3 s (d ), and the chiral nonunitary triplet, E 3 s (e), (along with their Hund's partners) will be favored.This degeneracy will be lifted by the subleading ferromagnetic spin fluctuations, which favor the latter state, E 3 s (e) (and its Hund's partner), as readily follows from Eq. ( 47).In the trivial representation, orbital fluctuations have no impact on which of the two possible phase diagrams in Fig. 2

D. In a magnetic field
Finally, we come back to the impact of fluctuation corrections on the leading triplet phase in the presence of a magnetic field.As we have seen in Sec.IV C, the superconducting state with the highest transition temperature in the presence of a sufficiently strong magnetic field will be a triplet phase due to the linear coupling in the second line of Eq. (42).At the mean-field level, b t 2 + b t 3 > 0, which prefers E 3 s (e) over E 3 s (c) as the order parameter of this phase.Using the relations in Eq. (C6), it is straightforward to rephrase the fluctuation corrections ( 47) and ( 50) of b j in terms of b t j → b t j + δb t j in Eq. ( 38).This yields δb t 1 + δb t 2 = χ (c + − c − ) 2 > 0 and δb t 1 + δb t 2 = 4δb > 0 for spin and orbital fluctuations, respectively.We conclude that, as expected, ferromagnetic fluctuations do not change the mean-field prediction in this case and E 3 s (e) is the dominant triplet order parameter in the presence of a magnetic field, for both strong and weak ferromagnetic fluctuations, and in their absence.

VI. ADDING FURTHER SYMMETRIES
In this section, we will analyze how the results presented above are modified once the additional symmetries, twofold rotation, C 2 , perpendicular to the plane of the system, and in-plane rotation symmetry, C 2y , are added.As shown in Figs.1(b) and 1(c), these symmetries are relevant as either exact microscopic or approximate emergent symmetries of twisted bilayer graphene and ABC trilayer graphene on hexagonal boron nitride, both of which exhibit superconductivity [2,48].

A. Consequences of a C 2 rotation symmetry
One crucial difference in twisted bilayer compared to twisted double bilayer graphene is that the former has an approximate C 2 symmetry [15] that mixes the two valleys, i.e., the system is (approximately) invariant under To relate to our notation used above, we assume that it is sufficient to focus on a single band for describing superconductivity in twisted bilayer graphene as well.This is quite a natural assumption and, unless stated otherwise, we expect our conclusions to hold when additional bands are taken into consideration.This (approximate) symmetry has attracted a lot of attention in the recent theory literature [13,14,17,22] of the system since it, combined with time reversal and C 3 , leads to a C 6 symmetry, which is responsible for not only the presence of (nearly gapless) Dirac cones at K and K but also the (approximate) vanishing of Berry curvature in twisted bilayer graphene.If the twist axis goes through the center of a hexagon, the system has C 6 rotation even as a microscopic symmetry.We note in passing that the (nearly) flat bands obtained in [57,58] for double bilayer graphene do not feature any Dirac cones but have well-separated conduction and valence bands that are characterized by nonzero Chern numbers (at least in some parameter regime); this strongly indicates that C 2 is not an approximate symmetry in twisted double bilayer graphene since C 2 would enforce zero Berry curvature.
In a similar fashion, [50] has argued that the twofold symmetry ( 51) is also an approximate symmetry for ABC trilayer graphene on hexagonal boron nitride, although it is clearly not a microscopic symmetry of the system, as can be seen in Fig. 1(c).
All things considered, it is currently not known whether an approximate C 2 symmetry is relevant for superconductivity in twisted bilayer and ABC trilayer graphene.Therefore, we will now discuss what changes for the possible superconducting instabilities once we assume that the Hamiltonian is also invariant under the transformation in Eq. ( 51).
The C 2 transformation plays a special role in two dimensions as it is equivalent to k → −k and can, thus, significantly affect superconducting instabilities [78].In graphene moiré superlattices, it also relates the two valleys and "interferes" with the Fermi-Dirac constraint (10): decomposing the pairing into singlet and triplet, Eq. ( 10) implies that λ s kv = λ s −k v and λ t kv = −λ t −k v .Consequently, it holds (as long as the pairing matrix elements between different bands can be neglected) that i.e., all representations even (odd) in C 2 must be pure singlet (triplet) states and vice versa.This has a few implications worth mentioning.First, even if C 2 is not a good symmetry (say, it is significantly broken by interactions), SU(2) s spinrotation invariance requires that the first transition must be into a pure singlet or triplet state and, hence, the pairing must be either even or odd under C 2 .In this sense, we can still distinguish between p-wave and d-wave pairing despite the presence of C 2 symmetry-breaking interactions.We emphasize that mixing will only be possible via multiple superconducting transitions (associated with admixtures of singlet and triplet) or interband pairing.The latter is expected to be quite weak given that the typical splitting between the bands at half filling (at least a few meV [5]) is about or more than an order of magnitude larger than the superconducting critical temperature ( 0.15 meV according to [2]).Secondly, if we do have an enhanced SU(2) + × SU(2) − symmetry (or are close to it), singlet and triplet are (nearly) degenerate.This forces the corresponding IRs of the spatial point group D 6 of the system, which behave identically under the subgroup D 3 but are even and odd under C 2 , to be (nearly) degenerate at the quadratic level of the Ginzburg-Landau expansion.For instance, A 1 and B 1 of D 6 have to be degenerate, as summarized in Table III.Without a Zeeman field, an extra C 2 symmetry with action in Eq. ( 53) also has no consequences for the higher-order terms in the free energy since spin-rotation invariance necessitates that all of these terms are even in the triplet vector.The only difference arises in the presence of a Zeeman field or magnetic fluctuations: with C 2 symmetry, it must hold that δc 1 = 0 in Eq. (24) and δc E 1 = 0 in Eq. ( 42) even when the SU(2) + × SU(2) − symmetry is broken.Furthermore, a C 2 symmetry implies c + = c − in Eq. (46).For this reason, weak ferromagnetic spin fluctuations do not lift the degeneracy of mean-field In summary, when classifying superconducting states in twisted bilayer graphene or ABC trilayer graphene on hexagonal boron nitride in the absence of a Zeeman field, it is unimportant whether an approximate C 2 symmetry is relevant or not: singlet and triplet will always be even and odd under it.We can thus work with D 3 (instead of D 6 ) without loss of generality in the following.The only difference with twisted double bilayer graphene (with finite displacement field) is an extra twofold rotation symmetry, C 2y , along the y axes [see Figs.1(b) and 1(c)].Its action on the electronic operators reads as where C 2y k = (−k x , k y ).The upshot of this additional symmetry for the possible superconducting instabilities is clarified in the next subsection.

B. D 3 versus C 3
Due to the additional C 2y symmetry, D 3 is a non-Abelian group and has three IRs-two one-dimensional and one twodimensional representations (refer to the character table in Table III).It is convenient to begin with the one-dimensional IRs A 1 and A 2 and take J = 0. Since C 2y interchanges the valleys, its action on the intervalley pairing order parameter (52) can be written as So, we see that a singlet (triplet) state transforming under A 1 (A 2 ) has no nodes while a singlet (triplet) in the A 2 (A 1 ) channel has symmetry-imposed nodes on the line k y = 0 and along the directions rotated by ±π/3.This creates six nodal points on any surface enclosing the point.We can also readily understand from Eq. ( 55) how the one-dimensional representations "connect" at the SU(2) + × SU(2) − point: at the high-symmetry point, λ s kv = λ t kv , ergo A 1 s  56); for instance, a possible choice for twisted bilayer graphene with Brillouin zone in Fig. 1(b) is given by (X k , Y k ) T = R (π +ϕ)/2 (X (1)  k , Y (1) k ) T with X (1)  k , Y (1)   k in Eq. (30a).To keep the notation short, each line with reference to ϕ 1 or ϕ 2 corresponds to two distinct states with ϕ 1 = 0, π/3 and ϕ 2 = 0, π/2.The indicated number of nodal points refers to a Fermi surface enclosing the point.

Pairing
M k+ Nodes around Hund's partner MF/FM and A 3 s 2 or A 1 s 2 and A 3 s 1 must meet at the J = 0 line in Fig. 2. We summarize these observations in Table IV.
In addition, for the case of the two-dimensional representation E of D 3 , the C 2y symmetry has nontrivial consequences.Once again, we take J = 0 which permits us to study singlet and triplet independently.As singlet pairing has already been analyzed in detail for twisted bilayer graphene (see, e.g., [30]), we are chiefly concerned with the triplet states here.We parametrize the triplet pairing as in Sec.IV A with the sole distinction being that the basis functions X k and Y k are now constrained by the symmetries of D 3 ; we choose them to obey X −C 2y k = −X k and Y −C 2y k = Y k , while transforming as k x and k y under C 3 .A possible choice is given by Eq. ( 30) with φ = π/2 for the Brillouin zone of twisted bilayer graphene in Fig. 1(b).With these conventions, the triplet vector transforms according to (d + , d − ) → (d − , d + ) under C 2y .This does not further constrain the quartic terms in the free energy (38), wherefore we can use the analysis of Sec.IV A for the point group C 3 , bearing in mind the caveat that the relative phase, ϕ, between d + and d − cannot be absorbed in a redefinition of the basis functions X k and Y k any more due to the extra C 2y symmetry.While ϕ has no consequences for E 3 s (d ) or E 3 s (e) and can be absorbed by performing a spin rotation for the phases E 3 s (b) and E 3 s (g), it describes different phases for all other stable minima of Eq. ( 38), and we have to go to higher order in the free-energy expansion to determine its value.Consider E 3 s (a) for instance.Writing d + = (1, 0, 0) T and d − = e iϕ (1, 0, 0) T , it is easy to verify that the most general, ϕ-dependent sextic term to the free energy must have the form c 1 cos(3ϕ) with c 1 ∈ R.This derives from Eq. (34) where the C 2y symmetry forces c 2 to vanish.We thus find ϕ = 2π n/3, n ∈ Z, for c 1 < 0 and ϕ = π/3 + 2π n/3 when c 1 > 0. These two minima correspond to two different states, which can be compactly represented by defining the "rotated" basis functions with X k and Y k as introduced above.The order parameters are for c 1 < 0 and c 1 > 0, respectively.We denote these two states by E 3 s (a) 0 and E 3 s (a) π 3 , respectively.The first state, E 3 s (a) 0 , preserves C 2y , but breaks C 3 rotation symmetry, and has a nodal line which is, as opposed to the states in Sec.IV A, pinned to k y = 0.The other state, E 3 s (a) π 3 , however, breaks C 2y and the nodal line is not pinned to the k x axis.
The remaining triplet states E 3 s (c), E 3 s ( f ), and E 3 s (h) of Sec.IV A can be analyzed in the same way.In all cases, we find two states corresponding to two different discrete values of the relative phase ϕ between d + and d − : for E 3 s (c) and E 3 s (h), we find ϕ = 0 or π/3 as before, whereas E 3 s ( f ) requires even higher-order terms in the free-energy expansion, yielding ϕ = 0 or π/2.In analogy to E 3 s (a) ϕ , we label the states by E 3 s (c) ϕ , E 3 s ( f ) ϕ , and E 3 s (h) ϕ ; their order parameters are the same as those of the corresponding states in Sec.IV A but with the rotated basis functions in Eq. ( 56) using the respective value of ϕ.Taken together, we obtain 12 triplet states for D 3 , which are summarized in Table IV, instead of only eight for the point group C 3 .
Finally, we can also ask how the different states behave for small J, i.e., whether singlet and triplet can mix and which phases are Hund's partners.Exactly as illustrated above for the pure triplet phases, we have to consider higher-order terms that determine the relative phase between the chiral, μ = +, and antichiral, μ = −, basis functions.As this analysis closely parallels our previous discussions, we just present the result in Table IV.In total, there are ten symmetryinequivalent mixed singlet and triplet phases.Seven of them are only possible if δa < 0 (singlet dominates); the remaining three can be realized for either sign of δa.

VII. DISCUSSION AND CONCLUSION
In this paper, we have presented a systematic classification and analysis of superconducting instabilities in graphene moiré systems.To this end, we have focused on zeromomentum Cooper pairs formed out of electrons in different valleys.Intervalley pairing is expected to be the dominant pairing channel as time reversal relates the two valleys.We have first analyzed singlet and triplet pairing separately since spin-orbit coupling is expected to be very weak in graphene.However, theoretical estimates of the interaction terms of twisted bilayer [13], double bilayer [57], and trilayer [50] graphene indicate that these systems are approximately invariant under independent spin rotations in the two valleys, leading to an (approximate) SU(2) + × SU(2) − symmetry and the (near) degeneracy of singlet and triplet pairing.For this reason, we have also classified the pairing instabilities close to this highsymmetry point, analyzing which triplet state transforms into which singlet phase upon changing the sign of the interactions breaking the SU(2) + × SU(2) − symmetry.We have further derived the conditions under which singlet and triplet can mix despite the absence of spin-orbit coupling.
As it has the fewest symmetries, we first considered twisted double bilayer graphene, for which there are also clear experimental indications of triplet pairing [51,52].Here, a displacement field, which is required to stabilize the superconducting state, reduces the point group to C 3 .The pairing states and their properties associated with the real representation A and the complex representation E of C 3 are summarized in Tables I and II, respectively.
Being one dimensional and real, A only allows for one singlet, a unitary and a nonunitary triplet phase, and one mixed phase.The latter is expected to be relevant only if SU(2) + × SU(2) − is weakly broken and the two consecutive transitions in the schematic phase diagram of Fig. 2(b) are very close.Using the values of the coupling constants in [57], we estimate the splitting to be about two orders of magnitude smaller than the critical temperature and, hence, hard to see experimentally [70].Whether renormalization-group corrections could enhance the impact of these weak symmetry-breaking perturbations at energies of order of the transition temperature is an open question, which we leave for future work.The gap structure of the four phases transforming under A is quite different: while the nonunitary triplet is gapless for one of the spin species, the singlet and unitary triplet have a single, fully established gap, and the mixed phase has two finite but distinct gaps for the two spin species.We have further shown that single-band mean-field theory will generically favor the phase diagram in Fig. 2(a) over Fig. 2(b).However, the small bandwidth and strong-coupling nature inherent in the problem makes the applicability of mean-field theory questionable and can lead to significant corrections which might eventually select other phases.We have illustrated these corrections for ferromagnetic fluctuations, expected to be relevant for twisted double bilayer graphene [51][52][53], twisted bilayer [4,61,62], and ABC trilayer graphene [63].We find that the resulting corrections will, as opposed to mean field, generally favor the phase diagram in Fig. 2(b) over that in Fig. 2(a).
The complex representation allows for many more states: two pure singlets, eight triplets, and, if SU(2) + × SU(2) − is only weakly broken, six distinct mixed phases.As compiled in Table II, all of these three classes of states allow for nodal points and fully gapped phases.However, only the triplets can have nodal lines (residual ungapped Fermi surfaces of one spin species).Only one out of the two different triplet states of the IR A allows for an admixture of singlet and triplet for weak anti-Hund's coupling but, in contrast, six out of the eight triplets transforming under E do so.
Out of the possible pairing states in Table II, single-band mean-field theory favors the two triplet states E 3 s (b) and E 3 s (d ) along with their respective Hund's partners-the nematic mixed phase E 1 s (0, i) + E 3 s (a) and the chiral singlet E 1 s (1, i).We show the associated phase diagrams in the vicinity of mean-field theory in Figs.5(a) and 5(b).We have discussed how additional weak ferromagnetic spin (orbital) fluctuations can lift the exact degeneracy of E 3 s (b) and E 3 s (d ), generically favoring the former (latter) and, hence, the phase diagram in Fig. 5(b) [Fig.5(a)].In the limit of strong ferromagnetic fluctuations, we obtain the chiral nonunitary triplet E 3 s (e) or, for weak anti-Hund's coupling, the mixed singlettriplet state E 1 s (1, i) + E 3 s (d ) as the dominant instability.
Motivated by the experimentally observed [52] linear increase of the transition temperature with an in-plane magnetic field in twisted double bilayer graphene and signs of magnetism in bilayer and trilayer graphene, we have also mapped out the possible phase diagrams in the presence of a magnetic field.As expected, if the SU(2) + × SU(2) − symmetry is significantly broken, the linear increase is only consistent with triplet pairing.For pairing in the A channel, there are two possible phase diagrams, shown in Figs.3(c) and 3(d), depending on which triplet state is realized in the absence of a magnetic field.The magnetic field fully determines the form of the leading triplet state to be A 3 s (1, i, 0) in the A channel.For order parameters transforming under E , there are two possibilities for the leading triplet state, E 3 s (c) or E 3 s (e), in a magnetic field; which of the two is realized depends on the value of the quartic couplings in the free energy.Both mean-field theory and ferromagnetic fluctuations favor the E 3 s (e) state.If, however, SU(2) + × SU(2) − is only very weakly broken, singlet pairing as the dominant instability of the system is also consistent with the linear increase of the critical temperature; the two possible phase diagrams for the case of pairing in the IR A are illustrated in Figs.3(a) and 3(b).
We have also derived (within mean-field theory) the key couplings, c 2 in Eq. ( 24) and c E 3 in Eq. ( 42), between the superconducting order parameter and the magnetic field B that determine the slope of the increase, T c , of the critical temperature with magnetic field.We found that they have the exact same mathematical form; as such, the behavior T c 2μ B B, with Bohr magneton μ B , seen in experiment [52], is equally surprising for both pairing channels and does not favor one channel over the other.In both cases, this might either be accidental or due to quantum critical scaling [57].
We have also studied, in Sec.VI, the changes in the classification when there is an extra in-plane rotation symmetry, C 2y , and a twofold rotation, C 2 , perpendicular to the plane.These two symmetries are relevant (either as exact or emergent symmetries) to twisted bilayer graphene and ABC trilayer graphene.We find that, while the C 2 symmetry has no consequences for the classification, C 2y not only pins the nodes of certain pairing states along high-symmetry lines but also leads to more pairing states as summarized in Table IV.
This paper further illustrates that graphene moiré systems provide a very rich playground for novel strongly correlated superconducting phases.We hope that our systematic analysis of pairing in the absence and presence of magnetic fields will help future theoretical and experimental studies to pinpoint the microscopic form of the superconducting state.
The triplet states E 3 s (b) and E 3 s (d ) will have the lowest energy for this configuration of quartic coefficients as argued in the main text.The degeneracy between these two states is lifted by corrections beyond the mean-field approximation, such as the ferromagnetic fluctuations of Sec.V.In the presence of a magnetic field, Eq. (A8) uniquely determines the chiral nonunitary triplet E 3 s (e) as the leading instability (see Sec. IV C).

Coupling to a magnetic field
In this subsection, we will analyze several important coupling terms between the superconductor and the magnetic field from a weak-coupling perspective.The microscopic form of the coupling to the Zeeman, M Z , and in-plane orbital field, M O , reads as where we have absorbed the g factor of the Zeeman coupling into the definition of M Z .This is not possible for the orbital coupling, as its g factor g v (k) depends significantly on momentum.The form of g v (k) is determined by microscopic details such as the Bloch states.All we need here is that g v (k) = −g v (−k), as follows from time-reversal symmetry (5), and we refer to [57] for a microscopic derivation of its momentum dependence.
Let us first note that even when the actual interacting multiband system is not invariant under C 2 the single-band mean-field Hamiltonian, H MF + H B , is left invariant under the action of C 2 in Eq. (51) if we further set d k → −d k in Eq. (A1) and M O → −M O .This emergent symmetry is a consequence of the special role of C 2 in two dimensions as it acts on k in the same manner as time reversal and, as such, can have crucial consequences for superconducting pairing [78].
In the present case, this symmetry implies that the coupling terms δc 1 in Eq. (24) and δc E 1 and c E 2 in Eq. ( 42) will vanish within single-band mean-field theory as is also readily confirmed by explicit calculation; we emphasize, however, that this is not an exact statement and we have checked that a multiband mean-field description allows for nonzero values.Nonetheless, we view the vanishing of these coupling in the weak-coupling single-band limit as an indication that they are likely small in the system.
Finally, the couplings of the Zeeman term to the triplet vector in Eqs. ( 24) and ( 42) are also not constrained by the emergent C 2 symmetry.We find these to be nonzero and given by respectively.Our main observation here is that the forms of c 2 and c E 3 are identical: the nonuniversal part is a momentum integral which, in both cases, is weighted by a function that is invariant under C 3 and has no symmetry-imposed nodes on the Fermi surface.Accordingly, it is not possible to distinguish between the IRs A and E based on the slope of the increase of T + c in small magnetic fields.

APPENDIX B: FLUCTUATION CORRECTIONS TO MEAN FIELD
In this Appendix, we provide further details on Sec.V.

Microscopic derivation
In this first part, we will derive, from a microscopic description of the system, that the prediction of the phenomenological approach of the main text provides the leading correction to the free energy of the superconductor in the limit where the mass of the fluctuations approaches zero.
To this end, we will use the field-theoretical formalism and describe the system by the action describes pairing, similar to the second line of Eq. (A1), where we omitted the term proportional to the order parameter squared, since it is irrelevant for the free-energy contribution at quartic order in the superconducting state.The ferromagnetic fluctuations in valley v with associated bosonic fields φ qv = (φ x qv , φ y qv , φ z qv ) T are described by the action where q ≡ (q, n ) is the bosonic analog of k, i.e., n = 2π T n are bosonic Matsubara frequencies.Physically, χ (q, i n ) plays the role of the (analytic continuation of the) dynamical spin susceptibility [as compared to the static one in Eq. ( 43) of the main text].We will focus here on the SU(2) + × SU(2) − symmetric limit, where [ χ (q)] vv = δ v,v χ (q), and take the conventional Ornstein-Zernike form with velocity v and mass m (related to the correlation length, ξ , as ξ = v/m).Finally, the spin fluctuations couple to the electrons as described by the last contribution, with a coupling constant g that, of course, can be absorbed into χ 0 (or vice versa), but we will keep g explicit here.
We follow [43] and focus on the leading-order correction of the coupling g (second order, ∝ g 2 ) to the free energy, F[ k , d k ], but emphasize that our expressions will differ from those of [43] as we consider a different type of fluctuations (centered around zero rather than finite momenta).These corrections are derived systematically by first integrating out the fermions, expanding the action to quadratic order in the bosonic fields φ vq , which is sufficient to quadratic order in g, and integrating out the massive bosons.
We are interested in terms quartic in the superconducting order parameter, leaving us with the four distinct types of contributions in leading order in g which are represented diagrammatically in Fig. 6.As indicated, the first diagram, in Fig. 6(a), only involves the zero-momentum and zerofrequency, q = 0, fluctuations and it exactly captures the contributions of the simplified approach discussed in Sec.V of the main text.
The remaining diagrams-the self-energy in Fig. 6(b), the ladder in Fig. 6(c), and the vertex correction in Fig. 6(d) to the mean-field box diagram-are fundamentally different: the loop integrals involve integration over finite frequency and momentum of the bosonic fluctuations.As such, it is intuitively clear that they are less singular in the limit m → 0 than the first diagram in Fig. 6(a), which is proportional to χ (q = 0, i n = 0) = m −2 .In fact, these additional contributions can be shown to diverge with log(m) for small m.To illustrate this, let us consider the diagram in Fig. 6(b), which is proportional to where we introduced the function that only depends on the fermionic momenta and frequencies.
Here, λ j k represent the basis functions of the involved superconducting order parameters.In the following, we will cut off the q integral by /v and expand ξ k+q ∼ ξ k + v k • q, allowing us to write with the dimensionless velocity ratio vk := |v k |/v.Since |ω n | π T and we work at finite temperature (given by the critical temperature of superconductivity), the term in the second line of Eq. (B9) is finite in the limit E → 0. The integral, thus, diverges as log(m) at small E (infrared), as stated above.The other two diagrams in Figs.6(c) and 6(d) can be analyzed in the same way and are also found to be subdominant as m → 0 compared to the one in Fig. 6(a).This justifies the approach of Sec.V microscopically.

Enhanced symmetry in the one-band description
In the last part of this Appendix, we discuss why c + c − in Eq. ( 46) is expected.From the previous subsection of this Appendix, we know that the results of the main text on fluctuation-induced superconductivity are captured by zeromomentum and zero-frequency fluctuations.Writing m v := φ q=0v in Eq. (B6) and generalizing to a momentum-and valley-dependent coupling constant, we here consider where g m v (k) = g m v (−k) as a consequence of time-reversal symmetry and we have, as before, assumed that we can focus on a single isolated electronic band.It is easy to see that H MF + H m , with H MF in Eq. (A1), is again invariant under the C 2 symmetry in Eq. ( 51) if we further replace While Eq. ( 44) is automatically invariant under Eq.(B11), the coupling for the two-dimensional representation in Eq. ( 46) is invariant only if c + = c − .Consequently, multiband effects are required for nonzero c + − c − , wherefore we expect its value to be much smaller than c + + c − , as stated in the main text.We also checked by explicit calculation that c + = c − is possible in a multiband description.

APPENDIX C: DETAILS FOR THE COMPLEX REPRESENTATION
In this Appendix, we present additional details of the different phases transforming under the complex representation E of C 3 .FIG. 7. Phase diagram for the free energy in Eq. (38).The different triplet states labeled (a) to (h) are defined in the main text in Sec.IV A.
As a starting point, it is helpful to chart out a phase diagram describing which of the triplet phases E 3 s (a) to E 3 s (h) is realized as a function of the quartic terms b t 1,2,3,4,5 in Eq. (38).Upon recognizing that b t 1 does not affect the form of the order parameter (but is assumed to be chosen so as to guarantee the stability of the expansion), we can conveniently display the phases as a function of b t j /|b t 2 |, j = 3, 4, 5, discussing the two possible signs of b t 2 separately.Such a phase diagram is drawn in Fig. 7.
As the main text contends, there are no independent terms involving σ y to add to the SU(2) + × SU(2) − invariant form of the free energy in Eq. (41).To see this, we note that it suffices to consider terms involving both + and − since terms with only + (or − ) have already been addressed in Sec.III A. Among the terms that mix + and − , the following are consistent with time-reversal and C 3 symmetry: However, all of these terms can be reformulated as It is not difficult to observe that the five different purely triplet quartic terms, β j=3, 4,5,6,7 , are all independent.Consequently, we can parametrize all 12 β j in terms of the five purely triplet terms and realize all of the triplet states of Sec.IV A.

FIG. 3 .
FIG.3.Phase diagram as a function of temperature T and Zeeman field M Z = M Z e x when (a, b) singlet dominates at low fields and (c, d) triplet dominates, which we determine by minimizing Eq.(28).Thin (thick) black lines correspond to second (first) order transitions.The phases for M Z = 0 are indicated in red and we recover the four different possible temperature dependences of Fig.2.Recall from Sec. III B 4 that b 2 > 0 and, thus, (a) or (c) is expected in mean-field theory.However, as we will see in Sec.V, strong ferromagnetic fluctuations will favor b 2 < 0 and, hence, (b) or (d).As symmetry requires δc 1 to be proportional to Hund's coupling J, we have set δc 1 = 0 here.For nonzero δc 1 , the singlet superconducting phases will contain an admixture of unitary triplet as described by Eq. (26) and a first-order transition into a singlet state (with unitary triplet admixture) will be possible at lower temperatures and nonzero Zeeman field in part (c).Note that the transition temperature from the normal state into the singlet superconductor is constant as we neglect here the nonlinear coupling to the magnetic field.A discussion of the latter can be found in Sec.III C 3.
FIG. 4. (i, ii)The lowest lattice harmonics of the basis functions [Eq.(30)].(a)-(h) The momentum dependence of the gap and the density of states, g(E ), for the pairing phases, E 3s (a) through E 3s (h).The Bogoliubov-de Gennes excitation spectrum is calculated using the band structure of the system predicted by the continuum model[58], assuming a pairing term of the form of Eq. (9) with an overall scale of 0 = 4 meV.The nodal points/lines are demarcated in dark blue; note that the states (b), (d), (g), and (h), taking α = π/4, are fully gapped; (a) and (f) have nodal points; and (c) and (e) have nodal lines, as is also visible in g(E ).

FIG. 5 .
FIG.5.The two possible phase diagrams of the complex representation close to mean-field theory, using the labeling of states defined in the main text and TableII.All transitions are second order, except for the one indicated by the thick line, which is first order.In Sec.V, we show that part (b) [part (a)] is favored when taking into account corrections to mean-field theory coming from ferromagnetic spin [orbital] fluctuations.
is realized.This results from the fact that neither M O (see Sec. III C) nor M z O can couple linearly to the superconducting states and their rotational invariant quadratic forms (M z O ) 2 and M 2 O can only couple to | s | 2 + d † d.Consequently, the energetic correction obtained by integrating out the orbital fluctuations will also only depend via | s | 2 + d † d on the superconducting states and, as such, not affect the value of b 2 in Eq. (20) and Fig. 2.

FIG. 6 .
FIG.6.Diagrammatic representation of the four different types of fluctuation corrections quartic in the superconducting order parameter (schematically represented by dashed lines and , * ) to leading order in g in Eq. (B6).The solid black lines with arrows are the bare electronic Green's functions, associated with S c , and the wavy lines denote the bosonic propagator defined in Eq. (B5).The diagram in (a) reproduces the contribution of the "phenomenological approach" of the main text, whereas the remaining diagrams in (b)-(d) are subleading in the limit m → 0.

TABLE III .
Character table of the point group D 3 together with the corresponding basis functions and IRs of D 6 for singlet/triplet pairing.C 2 should only be considered as an approximate symmetry and terms breaking this symmetry will lead to c + = c − , thus lifting the degeneracy of mean-field theory by, e.g., favoring Fig.5(b) over Fig.5(a).Note that the C 2 symmetry also forces c E O in Eq. (49) to vanish for M z O corresponding to valley fluctuations.As such, the approximate C 2 symmetry does not specify whether spin or valley fluctuations are expected to be the dominant source of lifting the mean-field degeneracy.It only indicates that strong ferromagnetic fluctuations are most likely dominated by spin fluctuations.

TABLE IV .
Summary of the different intervalley pairing states classified by the IRs of the point group D 3 .The notation closely parallels that of Table II.Here, we use λ 1 k and λ 2 k to denote continuous functions on the Brillouin zone that are even and odd under (k x , k y ) → (k x , −k y ), respectively, and are both invariant under C 3 , λ j k = λ j C 3 k .Furthermore, X ϕ k and Y ϕ k are rotated basis functions defined in Eq. (