Superconductivity from valley fluctuations and approximate SO(4) symmetry in a weak coupling theory of twisted bilayer graphene

The recent discovery of the Mott insulating and superconducting phases in twisted bilayer graphene has generated tremendous research interest. Here, we develop a weak coupling approach to the superconductivity in twisted bilayer graphene, starting from the Fermi liquid regime. A key observation is that near half filling, the fermiology consists of well nested Fermi pockets derived from opposite valleys, leading to enhanced valley fluctuation, which in turn can mediate superconductivity. This scenario is studied within the random phase approximation. We find that inter-valley electron pairing with either chiral (d + id mixed with p−ip) or helical form factor is the dominant instability. An approximate SO(4) spin-valley symmetry implies a near degeneracy of spin-singlet and triplet pairing. On increasing interactions, commensurate inter-valley coherence wave (IVCW) order can arise, with simultaneous condensation at the three M points in the Brillouin Zone, and a 2 × 2 pattern in real space. In simple treatments though, this leads to a full gap at fillings ± (1/2 + 1/8), slightly away from half-filling. The selection of spin-singlet or spin triplet orders, both for the IVCW and the superconductor, arise fcase corresponds to the Hundsrom SO(4) symmetry breaking terms. Mott insulators derived from phase fluctuating superconductors are also discussed, which exhibit both symmetry protected and intrinsic topological orders.


INTRODUCTION
There has been considerable interest in studying artificial lattices induced by a long wavelength Moiré potential in graphene and related materials. These experiments have recently gathered momentum with the observation of superconductivity and correlated Mott insulators in bilayer graphene twisted to a particular "magic angle". The Moiré superlattice induced in bilayer graphene twisted by a small angle leads to isolated bands near charge neutrality, whose bandwidth can be tuned by twist angle. [1][2][3][4][5][6][7][8][9][10][11][12] On approaching certain magic angles, the largest being 1.1°, the bandwidth is significantly reduced allowing for correlation physics to take hold. Indeed, recent studies on twisted bilayer graphene (tBLG) near the magic angle have revealed the presence of Mott insulators 13 at fractional filling of the bands, as well as superconductivity 14 in close proximity to some of the Mott insulators. While Mott physics has also been observed in a different Moiré superlattice system, induced by a boron nitride substrate on ABC trilayer graphene, 15 here we will focus on the tBLG system, which has already generated a significant amount of theoretical interest. [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] The band structure of tBLG at small twist angles can be understood from a continuum model 1,4,6 that couples the Dirac points in the individual graphene layers via the interlayer tunneling. Due to the small twist angles involved, there is a separation of scales between the atomic lattice and the Moiré superlattice which implies that commensuration effects can be neglected. 6 The opposite Dirac points in each layer are then essentially decoupled, leading to a valley quantum number n v = ±1 for each electron (n v = +1 for K valley and −1 for K′ valley), which is reversed under time reversal symmetry (as valleys are exchanged). Including both spin and valley degrees of freedom it takes 8 electrons (per Moiré unit cell) to completely fill the nearly flat bands that appear near neutrality. The additional factor of two in the filling appears due to band contacts present at neutrality and protected by symmetry. Charge neutrality then corresponds to four filled and four empty bands, which meet at Dirac cones at the K points of the Moiré Brillouin zone (MBZ). Measuring the electron charge density n from neutrality, the fully filled and fully empty bands occur at ±n s (~2.7 × 10 12 cm −2 for magic angle tBLG). In ref. 13,14 , an insulating state was also observed at f ¼ n=n s ¼ ∓ 1=2, i.e., at half filling both below and above neutrality (hence the term Mott insulator), where there were two (six) electrons per Moiré unit cell. Furthermore, superconductivity was observed around the f = −1/2 Mott insulator, i.e., around the Mott insulator on the hole doped side of neutrality.
Although interactions and the band width are both estimated to be comparable in magic angle tBLG, here we consider approaching the problem from the weak coupling limit, i.e., we imagine moving slightly away from the magic angle, which is motivated by the following considerations. First, although the energy scale of the bandwidth W 2 and interactions U 13 were estimated to be of order 10~20 meV, the Mott gap observed in transport experiments is much smaller~0.4 meV, and could be closed with an in-plane Zeeman field of roughly the same strength. Therefore, the system is not deep in the Mott regime, where the Mott gap would be of the same order as U. Next, doping the Mott insulator towards neutrality results very quickly in a metal with a big Fermi surface, where superconductivity is also observed. This regime could be approached from weak coupling. On the other hand, the side facing the band insulator (i.e., hole doping the f = −1/2 Mott insulator or electron doping the f = +1/ 2 Mott insulator) behaves like a "doped" Mott insulator, with both Hall conductivity and quantum oscillations pointing to a small Fermi surface composed of just the doped carries.
Finally, in both iron-pnictides [34][35][36] and overdoped cuprates, 37-39 weak coupling approaches have been relatively successful at least in predicting the gap symmetry. However, both these calculations relied on band structures with some degree of nesting. Does the fermiology of tBLG support such a nesting driven scenario? Interestingly, on moving slightly away from the magic angle, multiple band structure calculations 6,10,13,40 for small angle tBLG bands reveal a relatively strong nesting feature in the vicinity of half filling, albeit at wave-vectors not simply related to the filling. Such nesting is not expected in a single orbital model on the triangular lattice, but appears here quite generally from having opposite valleys that give rise to a pair of Fermi surfaces related by time reversal symmetry, each of which is constrained by the microscopic symmetries C 3 , M y and C 2 T as defined in ref. 19 . Within a random phase approximation (RPA), we show that nesting-enhanced valley fluctuations give rise to an inter-valley pairing in the "d/p-wave" channel (d-wave and p-wave are generally mixed under C 3 symmetry). An important ingredient is the presence of an approximate SO(4) symmetry. Although four component electrons (spin and valley) might suggest an SU (4) symmetry, this is strongly broken by the valley-dependent band structure. Instead, we obtain separate spin SU(2) rotation symmetries for the two valleys SU(2) K × SU(2) K′~S O(4) with interactions that only depend on the slowly varying part of the electron density. This symmetry ensures a degeneracy of the spin singlet and triplet inter-valley pairing (with valley indices adjusted to ensure the antisymmetry of the pair wave function). Further weak symmetry breaking terms are expected to split this degeneracy, the experimentally reported Pauli limiting behavior 14 suggests a spin singlet superconductor. This would require invoking a weak anti-Hunds coupling, leading to a inter-valley spin-singlet superconductor with chiral d + id and p−ip mixed pairing, while the more conventional Hunds coupling would favor a spin triplet superconductor with chiral d + id/p−ip pairing. 16 Note, in this setting, there is no symmetry distinction between d + id and p−ip pairing. However, depending on their relative strengths, a topological phase transition occurs characterized by different quantized thermal Hall conductivities (chiral central charge c = 4 vs c = −2). At strong coupling, or with explicit rotation symmetry breaking, nematic superconductivity with two or four nodes may also be stabilized.
Our general picture is illustrated the phase diagram in Fig. 1, in the vicinity of f = −1/2, which is obtained based on a mean-field model Eq. (12) to be discussed in details later. Tuning the twist angle θ towards the magic angle θ mag effectively decreases the ratio W/U between the band width W and the interaction U, which pushes the system towards strong coupling. Superconductivity will first emerge in the weak coupling regime. At stronger coupling, a simple nesting based picture predicts a inter-valley coherence wave order, with ordering wave vector at the three M points of the MBZ, although a full gap is opened only at filling f = −(1/2 + 1/8) or at 25% hole doping. A gap at half filling can open if interactions also modify the electronic dispersion, but this is outside the scope of the present treatment. The RPA approach does not apply to the strong coupling regime (as indicated in the phase diagram by the fading-out color), but we will also comment on alternative approaches that can tackle the strong correlation physics. We should also keep in mind that apart from the ratio W/ U, the twist angle θ also influence the band structure especially when θ gets close to the magic angle θ mag . Since the band structure becomes very sensitive to all kinds of perturbations near the magic angle, it is hard to draw universal features right at the magic angle. Thus we will stay a little bit away from the (first) magic angle θ mag by considering 1:2θ mag tθt2θ mag , which can provide us a relatively robust band structure and also place us closer into the weak coupling region in the phase diagram Fig. 1.
Given the volume of recent theoretical output we have to restrict our comments to a few selected references that are closest to this work. ref. 16 starts with an SU(4) Mott insulator, and predicted a topological superconductor on doping the Mott insulator. Our conclusions are similar, although we have an SO(4) (rather than SU(4)) symmetry, and we adopt a weak coupling approach which avoids conflict with localizing electrons in the narrow bands of tBLG. 19 As in reference 19 we favor a spin-singlet inter-valley ordering, albeit at a finite wave vector, and inter-valley fluctuations drive pairing of a spin-singlet superconductor. Finally, adding strong SO(4) symmetry breaking terms to our model reproduces the s-wave pairing in ref. 21 . Although 27,28,33 also predicts topological superconductor from weak coupling/quantum Monte Carlo, their models differ significantly from ours. Our proposed pairing mechanism based on the fluctuation of incipient order is similar to ref. 17 , while we identify the leading incipient order to be the valley fluctuation, which differs from the spin fluctuation in. 17 This paper is organized as follows. We start by proposing an effective model for tBLG, deriving the low-energy band structure in Sec. IIA and formulating the symmetry-allowed generic interaction in Sec. IIB. We then analyze the instabilities in all fermion-bilinear channels within the RPA approach in Sec. IIC and find a leading instability in the inter-valley coherence channel. We study valley fluctuation mediated pairing in Sec. IID and identify the dominant superconducting order parameter. We sketch two descriptions for the insulating phase adjacent to the superconducting phase: a Slater insulator with inter-valley coherence wave order in Sec. IIE and a topologically ordered Mott insulator obtained by projecting out charge fluctuations of the superconductor in Sec. IIF. Finally, we study the SO(4) symmetry breaking effects in Sec. IIG and close with a discussion in Sec. III.

RESULTS
Band structure and fermi surface nesting We first formulate an effective Hamiltonian that describes the electrons in the Moiré band near the Fermi surface. Our starting point is the continuum model of the tBLG proposed in ref. 4,5 which first focuses on the band structure around one valley (say Fig. 1 Schematic phase diagram in the vicinity of f = −1/2, which is obtained by self-consistent mean-field calculation according to Eq. (12) in the low temperature limit (details will be discussed later). TSC topological superconductor, IVCW inter-valley coherence wave. The strong coupling regime (closer to the magic angle) is not captured by this approach. We will mainly focus on the weak coupling regime in this work. The superconductivity is slightly stronger on approaching the van-Hove singularity which is on the electron doped side (neutrality is on the right) the K valley) where c Kl k denotes the K valley electron originated from the layer l (with l = ±1 labeling the top and the bottom layers respectively). h lk = v F (k−K l )⋅σ l captures the Dirac dispersion of the electron near the K l valley, where 0Þ by an angle φ l = lθ/2 determined by the twist angle θ, and accordingly σ l ¼ e Àiφ l σ 3 =2 σe iφ l σ 3 =2 is also rotated from the standard Pauli matrices σ = (σ 1 , σ 2 ) by the same angle. T q a ¼ w 0 þ w 1 ðq a σÞ Á b z þ iw 3 σ 3 describes the interlayer tunneling to the lowest-order of the momentum transfers, as specified by q 1 = K − −K + , q 2 = R 2π/3 q 1 and q 3 = R −2π/3 q 1 in Fig. 2a. In general, T q a depends on three real parameters w 0 , w 1 and w 3 (a typical setting is w 0 ≈ w 1 |q a | ≫ w 3 ). 1,4,6 Such a generic form of T q a can be pinned down by symmetry arguments given in ref. 19 . The Hamiltonian H K′ around the K′ valley can be obtained from H K simply by a time-reversal operation T : c Klk ! Kc K 0 l ;Àk (with K being the complex conjugation operator). Putting together, H 0 = H K + H K′ provides a full description of the low-energy electronic band structure of the tBLG in the continuum limit.
By diagonalizing the Hamiltonian H K (with an appropriate momentum cutoff), we obtain the single-particle band structure as shown in Fig. 3a. The bands are defined in the Moiré Brillouin zone (MBZ), as depicted in Fig. 2a with high symmetry points labeled.
The K + and K − valleys from either layers rest on the Moiré K M and K 0 M point respectively. We focus on the middle band around the charge neutrality, which will become flat as the twist angle θ approaches to the magic angles θ mag . A prominent feature of this band is that its energy contours (Fermi surfaces) around the −1/2 filling typically take triangular shapes around the Γ M point in the MBZ, as shown in Fig. 2a, which was observed in several band theory calculations for small twist angles. 6,10,13,40 The triangular distortion of the Fermi surface is generic on symmetry ground, as it is the lowest order (in terms of angular momentum) distortion that is consistent with all the valley-preserving lattice symmetries C 6 T and M y . 19,41 Indeed it is a rather robust feature for a range of twist angles θgt; rsim1:2θ mag and is also stable against perturbations like lattice relaxation, 42 as long as we are not too close to the magic angle. We assume that such triangular shape Fermi surface is relevant to the low-energy physics of the tBLG near the magic angle at −1/2 filling and base our analysis on this assumption. The key idea is that the almost parallel (well nested) sides of the triangular Fermi surfaces between K and K′ valleys could lead to strong valley fluctuations, which further provides the pairing glue for the superconductivity.
We describe a systematic procedure to extract an low-energy effective band structure from the continuum model described above. Briefly, the end result is a single band model with the dispersion ϵ K;k ¼ k 2 À μ þ αðk 3 x À 3k x k 2 y Þ around the K valley and ϵ K 0 ;k ¼ ϵ K;Àk around the K′ valley. In more detail, we proceed as follows. To model the triangular Fermi surface around the Γ M point, we first derive the effective band theory near Γ M . One systematic and unbiased approach is to first collect the singleparticle wave vectors |mk〉 around Γ M in the middle band (including both its upper and lower branches), and then construct a density matrix ρ / P mk jmkihmkj out of these states (note that | mk〉 are not orthogonal in the orbital space). By diagonalization ρ ¼ P i jψ i ip i hψ i j, we can identify the leading natural orbitals |ψ i 〉 (orbitals with largest weights pi ). The number n of the leading orbitals to be involved in the effective theory can be set by the desired fidelity level. To retain above 95% fidelity, s.t. P n i¼1 p i > 0:95, we typically need to take up to six orbitals (i.e., n = 6). Projecting the continuum model Eq. (1) to the six orbitals leads to the effective Hamiltonian where κ ± k ¼ v 1 ðk x σ 0 ± ik y σ 3 Þ and λ k ¼ ϵ 2 þ v 2 k Á σ are set by four real parameters ϵ 1;2 and v 1,2 . The band structure of the six-orbital model is shown in Fig. 3b. We can see that the features around Γ M is well captured compared to the continuum model in Fig. 3a, but the Dirac dispersions around K M and K 0 M can not be described by the six-orbital model (as expected). The six-orbital model provides a simpler and more flexible description of the near-Γ M band structure compared to the continuum model. 43 Its parameters can be determined by fitting to the first-principle calculations or experimental observations towards a more realistic modeling.
One can further simplify the six-orbital model by integrating out the high-energy electrons in the top and bottom bands, reducing the 6 × 6 Hamiltonian h Kk in Eq. (2) to its first 2 × 2 block: h 0 Kk ¼ ðϵ 1 À bk 2 Þσ 1 þ a Re k 3 þ σ 0 þ O½k 4 , which describes both branches of the middle band, where k ± ≡ k x ± ik y and the coefficients are If we only focus on the lower branch, the effective band theory boils down to a single-orbital model where we have chosen to rescaled the energy such that the single-orbital depends on only one tuning parameter α ¼ a=b ¼ 2ϵ 2 v 2 =ðϵ 2 2 À ϵ 2 1 Þ that characterizes the strength of the triangular Fermi surface anisotropy. The band structure of ϵ k is plotted in Fig.  3c. In Eq. (3), the K′ valley Hamiltonian is also included, which can be inferred from that of the K valley by the time-reversal symmetry T : c Kk ! Kc K 0 ;Àk . The Fermi surfaces in both valleys are drawn in Fig. 2b with μ = 1, α = 1/3 for example. One can see that the model essentially captures the triangular shape of the Fermi surface. There are three nesting vectors between K and K′ pockets, which are set by the chemical potential μ: p ; 0Þ and Note that the electronic spin degrees of freedom can be included in Eq.
In this single-orbital model, the notions of filling fraction and nesting commensurability are lost, but by going back to the original continuum model, we can identify the commensurate wavevector that has a high degree of nesting, which is found to be the M M points, i.e., Q 1 ≃ q 2 −q 1 /2. A commensurate perfect nesting will be achieved at the filling −5/8, which is hole-doped by 25% from the half-filling. We will show later in Sec. IIE that including a commensurate inter valley ordering with a period corresponding to the M M point of the MBZ, we can induce a full gap for relatively small order parameters, and obtain an insulating state when we are at the filling −(1/2 + 1/8) in the microscopic model given by the continuum theory Eq. (1).

Interactions and SO(4) symmetry
We now introduce interactions into the single-orbital model in Eq.
(3). As the electron c = (c K↑ , c K↓ , c K′↑ , c K′↓ ) in the MBZ carries both the spin (σ = ↑, ↓) and the valley (v = K, K′) degrees of freedom, one may expect an emergent U(4) symmetry at low energy that rotates all four components of the electron, as pointed out in ref. 16,19,20,29 . However, the electron kinetic energy (the band structure) strongly breaks this U(4) symmetry. For example, the triangular Fermi surface anisotropy α in the band Hamiltonian Eq.
(3) explicitly breaks the symmetry as the Fermi surface deformations are opposite between the two valleys as shown in Fig. 2. The U(4) symmetry is broken down to U(1) c × U(1) v × SO(4), where U (1) c is the charge U(1) symmetry generated by n c ¼ c y σ 00 c, U(1) v denotes the emergent valley U(1) symmetry generated by n v ¼ c y σ 30 c and SO(4)~SU(2) K × SU(2) K′ stands for the two independent SU(2) spin rotation symmetries in both valleys generated by S v ¼ c y v σc v (for v = K, K′ separately). The original SU(4) generators that are broken by the Fermi surface anisotropy α form a (complex) SO(4) vector, which corresponds to the intervalley coherence (IVC) order I μ ¼ c y K s μ c K 0 as proposed in ref. 19 , where s μ are defined to be (s 0 , s 1 , s 2 , s 3 ) ≡ (σ 0 , −iσ 1 , −iσ 2 , −iσ 3 ) for σ μ being the Pauli matrices. The pairing channels can also be classified by the SO(4) symmetry. There are only two possibilities: ) that transforms as SO(4) (pseudo)scalar. These operators are summarized in Table 1, which exhaust all fermion bilinear operators that can be written down on a local Wannier orbital. (See the Supplementary Material I for a more detailed classification of fermion bilinear operators.) Therefore any U(1) c × U(1) v × SO(4) symmetric local interaction should be mediated by one of these fermion bilinear channels. Further taken into account the time-reversal symmetry T (that interchanges valleys), it turns out that there are only two linearly independent and symmetric interactions (see the Supplementary Material I for the derivation of independent local interactions), which can be written purely in terms of density-density interactions as where n vq P k;σ c y vσkþq c vσk is the density operator of each valley. Since the density-density interaction is generally repulsive, we expect both parameters U 0 and U 1 to be positive (typically U 0 ≈ U 1 > 0). At the special point of U 0 = U 1 = U, the U(4) symmetry is restored for the interaction H int . However, even if H int is tuned to the U(4) symmetric point, when combined with the kinetic energy H 0 in Eq. (3), the symmetry of the full Hamiltonian (4). Later in Sec. IIG, we will further discuss the effect of adding small interaction terms to finally break the emergent SO(4) symmetry down to the microscopic SO(3) spin rotation symmetry.
In summary, by putting together Eqs. (3) and (4), we propose an effective model H = H 0 + H int for the tBLG with Fermi level resting in the lower branch of the nearly-flat band but not too close to the charge neutrality (such that the Fermi surface is still within the control of Γ M point expansion). More specifically, we assume that the Fermi level does not go beyond the van Hove singularity that separates Fermi pockets around the K M points near charge neutrality from those centered around Γ M , see also Fig. 2a. Our remaining goal is to analyze the model within a weak coupling approach.
Random phase approximation We calculate the renormalized interactions within the random phase approximation (RPA) [34][35][36] to analyze the electron instabilities in all six fermion bilinear channels as enumerated in Table 1. We will first restrict our analysis within the s-wave channels for simplicity. For each fermion bilinear operator A q ¼ 1 2 P k χ Àkþq Aχ k generally expressed in the Majorana basis χ k , we evaluate its bare static (zero frequency) susceptibility on the ground state of the single-orbital model H 0 . Then we rewrite the interaction H int ¼ g 0 P q A y q A q þ Á Á Á in the same channel to extract the bare coupling g 0 . The RPA corrected coupling is then given by g RPA (q) = g 0 (1 + g 0 χ 0 (q)) −1 . Admittedly, the RPA approach may not capture the interwind fluctuations in different channels. More systematic and unbiased approaches such as the function renormalization group 44 could be implemented to improve the result in the future.
The largest (in magnitude) value of g RPA (q) is taken and plotted in Fig. 4a as a function of U 0 = U 1 = U for various channels. The most attractive coupling appears in the IVC channel, which is associated with the operator I μ q ¼ Figure 4b shows the bare susceptibility of the IVC fluctuation and Fig. 4c is its RPA corrected coupling, which peaks strongly around three momentums that exactly correspond to the nesting momentums Q 1,2,3 . So Electrical charge is labeled by q c , thus q c = 0 corresponds to charge neutral (particle-hole) operators, while q c = 2 corresponds to Cooper pair (particleparticle) operators. The valley quantum number of the U(1) v symmetry is labeled by q v , hence inter-valley coherence order is obtained on condensing q v = 2 operators. Non-Abelian symmetry representations are labeled by the dimension (with a prime to denote the pseudorepresentation) Y.-Z. You and A. Vishwanath as the bare interaction is strong enough, I μ will condense at these momentums, leading to a finite-momentum IVC order, which we called the inter-valley coherence wave (IVCW). Suppose the nesting vector is pinned by the Moiré pattern to M M . Upon doping, the nesting condition will quickly deteriorate and the IVCW order will cease to develop. Nevertheless the low-energy valley fluctuations can play the role of the pairing glue, mediating an effective pairing interaction between electrons. A hint that can already be observed from Fig. 4 in which the attractive coupling diverges in the I μ channel, while at the same time a repulsive coupling in the s-wave inter-valley pairing Δ μ channel also diverges. This implies that if the pairing form factor is allowed to change sign along the Fermi surface (which goes beyond s-wave), the repulsive coupling in this pairing channel can be effectively converted to an attractive one, leading to a strong pairing instability based on the Kohn-Luttinger mechanism. 45,46 The details will be discussed in the following.

Superconductivity
To pin down the pairing instability mediated by the valley fluctuations, we take the RPA corrected interaction in the IVCW channel I μy I μ and recast it in the inter-valley pairing channel Δ μy Δ μ (restricting to the zero momentum pairing c k c −k ) X q;μ where Kk σ 2 s μ c K 0 Àk is the inter-valley pairing operator, recall that (s 0 , s) = (σ 0 , −iσ). The attractive interaction (g RPA < 0) in the IVCW channel implies the repulsive interaction (−g RPA > 0) between Δ μ k and Δ μ Àkþq . So the pairing can gain energy only if there is a relative sign change between the pairing form factors connected by the nesting momentums Q a (at which the scattering is the strongest), i.e., Δ μ k ¼ ÀΔ μ ÀkþQa , as illustrated in Fig. 5a. By solving the linearized gap equation, X the leading gap function (i.e., the eigen function Δ μ k with the largest eigenvalue λ) is found to be of the form where u μ and v μ are complex vectors, and the form factor w k ¼ w d k 2 þ þ w p k À is a linear combination of the d + id and the p−ip waves with real coefficients w d and w p , as shown in Fig. 5b. The mixing between the d + id and the p−ip pairing is generic, because in the presence of the triangular Fermi surface distortion α, the angular momentum is only mod 3 conserved, meaning that there is no distinction between the d + id and the p − ip wave on symmetry ground. The ratio |w p /w d | carries the dimension of momentum and sets a momentum scale k Q = |w p /w d |, which is expected to be associated with the nesting momentum k Q ≃ |Q a |/ 2. The form factor w k has three zeros (vortices) on the circle of k Q in the momentum space. If the Fermi surface circumvents the zeros from outside (or inside), the pairing will be dominated by d + id (or p − ip) wave.
To determine the coefficients u μ and v μ in Eq. (7), we can write down the Landau-Ginzburg (LG) free energy F within the meanfield theory, 16 (see also the Supplementary Material I,I for the derivation of Landau-Ginzburg free energy and a more detailed analysis of competing orders) As studied in ref. 16 , the free energy admits two types of minimum, which are degenerated in energy, chiral : where ϕ, ϕ 1 , ϕ 2 are arbitrary phases and n μ ; n μ 1 ; n μ 2 are real O(4) vectors with n μ 1 n μ 2 ¼ 0. The chiral solution preferentially choose the form factor of one chirality (either w k or w Ã k ), which corresponds to four copies of the d + id or the p−ip superconductors (or its timereversal partners). The helical solution is a superposition of w k (in one spin sector) and w Ã k (in the other spin sector), which corresponds to two copies of the d ± id or the p ∓ ip superconductors.
In the valley and spin space, Δ μ k transforms as a (complex) SO(4) vector, whose four components corresponds to the spin-singlet pairing Δ 0 k and the spin-triplet pairing In the presence of the emergent SO(4) symmetry, the singlet and triplet pairings are degenerated. This can be considered as an SO(4) generalization of the SO(3) pairing Δ k proposed in ref. 16 , such that the singlet pairing is also included as a possible option in our discussion. However, the SO(4) symmetry is not exact in the tBLG. Any inter-valley spin-spin interaction will break the SO(4) symmetry down to the global (valley-locked) SO(3) spin rotation symmetry, and thus splits the degeneracy between singlet and triplet pairings. If the singlet pairing is favored, then only the chiral gap function is possible, because there is no room for two perpendicular O(4) vectors n μ 1 and n μ 2 to coexist just in the singlet channel. If the triplet pairing is favored, then both the chiral and helical gap functions are allowed. We will discuss the effective of explicit SO(4) symmetry breaking in more details later.  In general, the superconductor will be a topological superconductor (TSC) with fully gapped Fermi surface. [47][48][49] The chiral TSC breaks the time-reversal symmetry and also breaks the U (1) c × U(1) v × SO(4) symmetry to Z F 2 Uð1Þ v SOð3Þ. The topological classification for the chiral TSC is Z. If the d + id (or p − ip) pairing is stronger, the topological index will be ν = 8 (or ν = −4), which admits 8 (or 4) chiral Majorana edge modes. The helical (non-chiral) TSC preserves the (spin-flipping) time-reversal symmetry Z T 2 (under which c Kk ! Kiσ 2 c K 0 ;Àk ; c K 0 k ! Kiσ 2 c K;Àk ) and breaks the U(1) c × U(1) v × SO(4) symmetry to Z F 2 Uð1Þ v SOð2Þ. The SO(2) symmetry may be loosely called a spin U(1) s symmetry since it corresponds to a joint spin rotation for both valleys (in either the same or the opposite manner). In the presence of both U(1) v and U(1) s , the topological classification of the helical TSC is also Z. If the d ± id (or p ∓ ip) pairing is stronger, the topological index will be ν = 4 (or ν = −2), which admits 4 (or 2) helical Majorana edge modes. It is also possible to fine tune the ratio w d / w p to the topological phase transition between the d-wave and pwave TSC, then superconducting gap will close at the nodal points on the Fermi surface resulting in 12 Majorana cones in the bulk.
Finally, we would like to briefly comment on the possibility of the nematic d-wave or p-wave pairing. We could go beyond the mean-field theory by considering more general momentumdependent quartic terms in the LG free energy X k;k 0 If κ kk′ satisfies P k;k 0 κ kk 0 ðw Ã k w k 0 Þ 2 <0, the LG free energy will have only one type of minimum, (see Supplemental Information I,I for justifications of the assumption and the solution) nematic : where ϕ 1 , ϕ 2 are arbitrary phases and n μ is a real O(4) vector. This solution corresponds to the nodal d-wave or p-wave pairing, as Δ μ k $ Reðe iðϕ 1 Àϕ 2 Þ w k Þn μ , which preserves the time-reversal symmetry and breaks the U(1) c × U(1) v × SO(4) symmetry down to Z F 2 Uð1Þ v SOð3Þ. The nodal line lies along the direction set by ϕ 1 −ϕ 2 , which breaks the C 3 rotational symmetry. So the nodal superconductor also has a "nematic" (orientational) order. 50,51 As the Fermi surface is not fully gapped, the nematic superconductor is not topological and has no protected edge mode. Apart from strong coupling, explicit breaking of C 3 rotation symmetry could also favor nematic superconductivity.

Slater insulator and valley order
When the Fermi surface is tune to optimal nesting, the strong nesting instability could lead to the condensation of the IVC order parameter I μ at the nesting momentums, which drives the system into the IVCW phase. In the weak coupling approach, the IVCW and the TSC order compete for the Fermi surface density of state. Here we provide a mean-field theory calculation that captures both competing orders and gives a rough estimate of the overall structure of the phase diagram. We start with the mean-field Hamiltonian H MF that incorporates the order parameters of both the IVCW I 0 Q and the TSC Δ 0 k (which are restricted to the singlet channel without loss of generality given the SO(4) symmetry), where H 0 is taken to be the single-orbital model Eq. (3) and Q is summed over the three nesting vectors Q 1,2,3 . g I = g RPA (Q) and g Δ = avg k,k′∈FS g RPA (k + k′) are the effective couplings in the IVC and the pairing channels respectively. Both of them originate from the RPA corrected coupling g RPA (q) and are expected to scale together with the interaction strength U = U 0 = U 1 . By tracing out the electron, we obtain the free energy F ¼ Àβ À1 ln Tre ÀβHMF for the order parameters I 0 Q and Δ 0 k . (See the Supplementary Material I, I,I for details about the self-consistent mean-field calculation.) We find the free energy saddle point solution in the low temperature limit for different W=U $ g À1 I ; g À1 Δ (where W is the band width) and different chemical potentials μ. This allows us to map out the mean-field phase diagram (in the zero temperature limit) as shown in Fig. 1. As we tune the twist angle towards the magic angle, the band gets flatten and the effective coupling increases. The TSC phase will first appear at low temperature. With stronger coupling, the IVCW phase will emerge around the optimal nesting and gradually expand to a wider range of chemical potential.
As we fix the couplings at g I = 0.8 and g Δ = 0.4 (the energy unit is set by the band dispersion in H 0 ), assume that the optimal nesting is around μ = 1 (such that the nesting momentum is jQj ¼ ffiffiffiffiffi 3μ p % 1:73), and take the anisotropy parameter to be α = 1/3, we can obtain a mean-field phase diagram (for finite temperature) as shown in Fig. 6 (by solving the free energy saddle point equations). The fermilogy at different representative points in the phase diagram are shown in Fig. 6. In the metallic phase, the Fermi surface consists of electron pockets around K and K′ valleys (drawn together). In the TSC phase, the Fermi surface is gapped by the inter-valley pairing with the pairing form factor shown in color (following Fig. 5b). The pairing can be either chiral or helical within the mean-field theory. In the IVCW phase, the K′ pocket (in light red) is shifted away from the K pocket (in light blue) by the three nesting vectors Q 1,2,3 . Deep in the IVCW phase, the Fermi surface can be fully gapped. In between T IVCW and T ins , small (reconstructed) hole or electron pockets remain on the Fermi level. However, using the single-orbital model Eq. (3) as the starting point, we have lost track of the notion of the Moiré Brillouin zone (MBZ) and we can not tell if the nesting vector Q is commensurate with the Moiré lattice or not.
To further investigate the commensurability of the nesting vector and the corresponding filling of IVCW state, we have to fall back on the continuum model Eq. (1) for H 0 , such that the MBZ can be referred. We would like to explore if the IVCW order can fully gap out the Fermi surface and lead to an insulator. We will first focus on the commensurate IVCW order. From the shape of the Fig. 6 Mean-field phase diagram in the vicinity of f = −1/2 and at finite temperature. TSC topological superconductor, IVCW intervalley coherence wave. The TSC appears below T c around the IVCW insulator on both the hole and electron doped sides, with a d + id and p−ip mixed inter-valley pairing. The IVCW order on sets at the temperature T IVCW and becomes strong enough to full gap out the Fermi surface below T ins . On the hole doping side, the metallic IVCW phase has a single hole pocket with twofold spin degeneracy. The transition temperatures T c and T IVCW are correlated since they arise from the same interaction g RPA Fermi surfaces in Fig. 2a, the nesting vectors are most likely to be commensurate if they connect the Γ M point to the M M points in the MBZ. With this, we consider the IVCW order where the valley fluctuations simultaneously develops at the three M M points in the MBZ (corresponding to the nesting vector Q 1 = q 2 −q 1 /2 and its C 3 related partners Q 2,3 ).
The commensurate IVCW order breaks both the U(1) v × SO(4) symmetry and the translation symmetry. It leads to a 2 × 2 modulation on the Moiré lattice as demonstrated in Fig. 7a. As the unit-cell is enlarges to four Moiré sites, the Brillouin zone will be reduced to 1/4 of the MBZ, as illustrated in Fig. 7b. The lower branch of the band (from charge neutrality to the band bottom) will be folded to eight bands in the reduced Brillouin zone (rBZ), which consist of four folded bands for each valley. As we turn on the IVCW order to mix the K and K′ valleys together, a full gap opens between the third and the fourth bands (counting from bottom) as shown in Fig. 7c, d. Counting from the charge neutrality, this corresponds to the filling f = −5/8, but not the filling f = −1/2 as one may expect. In fact, the −1/2 level lies in the continuum above the IVCW gap, as indicated in Fig. 7c, d. At the filling −5/8, the system becomes an IVCW ordered band insulator, which may be called a Slater insulator (to be distinguished from the Mott insulator). There is a simple geometric picture to explain the seemly strange −5/8 filling. In the ideal case, if we consider the K and K′ pockets to be straight triangles connecting the M M points, illustrated as the dashed lines in Fig. 7b, the nesting will be perfect at the desired M M momentum and the corresponding filling is indeed −5/8 by counting the areas of the triangles. Therefore, although the commensurate IVCW order can lead to a fully gapped insulator, but the filling of the insulator has a 1/8 deficit from the −1/2 filling. We also checked that if the ordering momentum is changed to the Γ M or K M point momentum, no gap opening is observed with weak to medium IVCW order. While the −5/8 filling sounds peculiar, we note that in a recent experiment 52 of tBLG, separate quantum oscillations (Landau fans) are observed to emerge from f = −1/2 and f ≈−0.6, the later of which is closer to f = −5/8 = −0.625, although more evidences are still needed to verify or falsify this insulating state as a commensurate IVCW state.
However, if we go beyond the commensurate nesting and relax the nesting vector from the M M momentum, it is possible to obtain an incommensurate IVCW insulator for a range of fillings around −5/8, including the −1/2 filling typically, as long as the nesting condition is good. Another possibility is that the band structure may receive self-energy corrections from the interaction in such a way that the −1/2 filling Fermi surface turns out to admit good commensurate nesting. But in either picture, the −1/2 filling is not special compared to other fillings in terms of forming a Slater insulator, which still does not provide a natural explanation for the specific filling of the Mott insulator. This suggests that the Mott insulator in the tBLG might be a strongly correlated phase beyond the weak coupling picture like Fermi surface nesting. In this case, a strong coupling approach is required to understand the observed Mott insulator at precisely −1/2 filling. Below we discuss a scenario of Mott insulator that naturally arise from quantum disordering the adjacent superconducting phase by double-vortex condensation. [53][54][55][56][57] Mott insulator and topological order One approach towards a strong-coupling Mott state is to start from the adjacent superconducting state and then suppress the U (1) c charge fluctuation by proliferating double vortices of the superconductivity (SC) order parameter (or equivalently 2π fluxes seen by the electron). [53][54][55][56][57] Single vortices of the SC order parameter become anyonic excitations in the resulting Mott state, such that the Mott phase acquires intrinsic topological order. 58,59 In this approach, the nature of the topological order in the Mott phase will be closely related to the nature of the SC order in the adjacent SC phase. Here we assume that the nature of the SC order will remain qualitatively the same as we increase the interaction strength from the weak-coupling to the strongcoupling regime. This assumption is consistent with the past experience of unconventional superconductors including cuprates and iron-pnictides. 60 Assuming this, we can take the SC orders obtained from the weak-coupling approach as input to provide us with more insights about the possible orders in the Mott phase.
On the field theory level, this amounts to first fractionalizing the electron c vσ into a bosonic parton b and a fermionic parton f vσ as c vσ = bf vσ following a slave-boson approach, [61][62][63][64] where v = K, K′ labels the valley and σ = ↑,↓ labels the spin. Both bosonic and fermionic partons couple to the emergent gauge field. We assign the U(1) c symmetry charge to the bosonic parton and the U(1) v × SO(4) symmetry charge to the fermionic parton, in close analogy to the spin-charge separation in cuprates. [65][66][67] The fermionic parton is assumed to be in one of the SC state, such that once the bosonic parton condenses, the electronic SC state will be recovered. As we go from the (electronic) SC phase to the Mott phase, the bosonic parton is expected to acquire a gap across the transition, such that the charge fluctuations will be gapped and the U(1) c symmetry will be restored in the Mott phase. Then the fermionic parton SC state essentially becomes a (generalized version of) quantum spin liquid with intrinsic topological order and symmetry fractionalization 68-72 of valley and spin quantum numbers. Hence such a Mott state may be called a valley-spin liquid (VSL). Different types of SC states correspond to different types of Mott states, as summarized in Table 2. On the other hand, charge doping the VSL states will drive the bosonic parton condensation 〈b〉 ≠ 0, which identifies the fermionic parton f vσ  with the electron c vσ = 〈b〉f vσ , and converts the topological order back to the corresponding SC order. So the correspondence between the SC states and the Mott states in Table 2 are mutually consistent.
The chiral VSL sate can be viewed as the d + id (or p−ip) chiral TSC state of the fermionic parton, which enjoys the SO(8) 1 (or SO (4) −1 ) topological order. 73 They admit Abelian Chern-Simon theory 74-78 descriptions L CS ¼ 1 4π K IJ a I^d a J with the K matrices given by where A v and A s are the background fields that probe the U(1) v × U (1) s symmetry. The Z 2 topological order have four anyon sectors: 1, e, m and ε, where e and m are bosons with mutual-semionic statistics, and they fuse to the fermionic parton ε. For the p-wave helical VSL, e and m must separately carry either the U(1) v or the U (1) s symmetry charge, and ε carries both charges. The mutualsemionic statistics between e and m implies that the p-wave helical VSL state will have a fractionalized valley-spin Hall conductance σ vsH = −1/2. Moreover, because the fermionic spinon ε is a Kramers doublet T 2 ¼ À1 À Á under the timereversal symmetry, it must be the case that one of e or m is a Kramers doublet and the other one is a Kramers singlet T 2 ¼ þ1 À Á , such that the time-reversal anomaly vanishes. 87,88 So the p-wave helical VSL state is a Uð1Þ v Uð1Þ s Z T 2 symmetry enriched topological (SET) state. [89][90][91] For the d-wave helical VSL, m can be charge neutral and Kramers singlet, whereas e and both carry the U(1) v × U(1) s charge and are Kramers doublet. This can be viewed as a trivial Z 2 topological order on top of a U(1) v × U(1) s bosonic symmetry protected topological (BSPT) state. 77,[92][93][94][95][96][97] The Z 2 topological order can be removed by condensing the charge neutral boson m. Then the Mott insulator simply realizes a U(1) v × U(1) s BSPT state with quantized valley-spin Hall conductance σ vsH = 1.
Finally, if we start with the nematic superconductor, the corresponding Mott state will be a gapless Z 2 VSL with nodal fermionic partons and gapped visons. 50 The symmetry of this VSL state is Uð1Þ c Uð1Þ v SOð3Þ Z T 2 . Like the nematic superconductor, the C 3 rotation symmetry is still broken in the VSL state, so there will be a coexisting nematic order in this Mott insulator.
In all cases, the emergent SO(4) symmetry is broken in the Mott phase. But the remaining symmetry is still sufficient to protect a two-fold degeneracy of the electron. For the chiral VSL, the electron transforms (projectively) as spin-1/2 (spinor representation) of the SO(3) symmetry. For the helical VSL, the electron forms Kramers doublet under the time-reversal symmetry. For the nematic VSL, both SO(3) and time-reversal protections are present. The symmetry protected two-fold degeneracy in the valley-spin space is consistent with the experimentally observed Landau fan 14 near the Mott phase with the filling-factor sequence 2; 4; 6; … Consider for example, the spin singlet VSL phase, which is connected to the spin singlet chiral superconductor. Here, spin degeneracy is present, and although valley remains a good quantum number, since the phase itself breaks time reversal symmetry, the degeneracy between opposite valleys is lost. Although it is hard to estimate the strength of this effect, the symmetry dictated degeneracy is just twofold.
Breaking SO(4) symmetry Both the IVCW and the TSC phases break the emergent SO(4) symmetry, as their order parameters I μ and Δ μ are SO(4) vectors. The four (complex) components of the order parameters correspond to the orderings in the spin-singlet and the spintriplet channels, which are degenerated in the presence of the SO (4) symmetry. However, the SO(4) symmetry is never exact in reality. The explicit SO(4) symmetry breaking can split the degeneracy. We will analyze the effects of the SO(4) symmetry breaking in the following.
We first consider the Heisenberg spin-spin interaction between valleys, where S vq ¼ P k c y vkþq σc vk (for v = K, K′) is the spin operator. The J (q) < 0 (or J(q) > 0) case corresponds to the Hunds (or anti-Hunds) interaction. It belongs to the (1,1) representation (the symmetric rank-2 tensor) of the SO(4) ≃ SU(2) K × SU(2) K′ group, which locks the two SU(2) subgroups together and breaks the SO(4) symmetry down to SO(3). The interaction H J admits decompositions in the IVC and the pairing channel as where Δ 0 and I 0 are the spin-singlet orderings (as SO(3) scalar), and Δ and I are the spin-triplet orderings (as SO(3) vector).
Depending on the sign of the inter-valley Heisenberg interaction J (q), the spin-triplet (or spin-singlet) pairing is favored if the interaction is Hunds (or anti-Hunds) like. If we assume an anti-Hunds interaction (i.e., J(q) > 0), then according to Eq. (16), the interactiont will provide attractive interactions for both the IVC and the pairing in the spin-singlet channel. The anti-Hunds interaction could arise from the renormalized Hubbard interaction by integrating out high energy electrons as proposed in ref. 21,98 .
Note that the spin-singlet TSC can only be a chiral TSC as discussed in Sec. IID previously. However, if the inter-valley interaction turns out to be Hunds like, the spin-triplet pairing could also be favored, which admits both chiral and helical TSC. The possibilities are summarized in Table 3. However, if the SO(4) symmetry breaking is implemented in the IVC channel, the result can be very different. Suppose we consider the following enhanced attraction (i.e., g(q) < 0) in the I 0 channel, so as to single out the spin-singlet IVCW order. The same interaction would be translated into the pairing channel as which is also an attractive interaction in the spin-singlet pairing channel Δ 0 (note that g < 0). In contrast to Eq. (5), only the I 0y I 0 interaction is involved in Eq. (17), which completely changes the interaction sign in the singlet pairing channel. Under the RPA correction, g(q) peaks strongly around the nesting momentums q = Q 1,2,3 , thus the attractive interaction between Δ 0 k and Δ 0 ÀkþQa effectively reduces the energy gain in the singlet channel, due to the sign-changing TSC pairing form factor (i.e., Δ 0 k ¼ ÀΔ 0 ÀkþQa ). Therefore a slightly enhanced attractive interaction in the spinsinglet IVCW channel will actually suppresses the spin-singlet TSC pairing and favors the spin-triplet TSC pairing, as summarized in Table 3. The spin-triplet TSC can be either chiral or helical as discussed Sec. IID previously. Although H J in Eq. (16) and H g in Eq. (17) are both SO(4) symmetry breaking terms in the (1,1) representation that favor the singlet IVCW order, yet their effects on splitting the singlet-triplet degeneracy in the TSC channel is completely opposite. This has to do with the fact that under the RPA correction, the interaction H J in the spin channel is not sensitive to the nesting effect, but the interaction H g in the valley channel exhibit a strong nesting effect. This results in very different momentum-dependence of their coupling functions (J(q) or g(q)), which finally divide the fate of the singlet-triplet splitting. The competition between these two symmetry breaking effects demands further analysis by more refined approach such as the functional renormalization group, 99,100 which will be left for future works. 101 Finally, we would like to comment on the connection to ref. 21 , where the valley XY interaction H g in Eq. (17) was considered to be the dominant interaction in the tBLG. In this case, the emergent SO(4) symmetry is strongly broken. The effective attraction in the spin-singlet pairing channel can simply drive the s-wave valleysymmetric spin-singlet pairing, which then leads to a nontopological superconductor as in, 21. Therefore, whether the superconductivity in the tBLG is topological or not could sensitively depend on the form and the strength of the SO(4) symmetry breaking interactions, as summarized in Table 3.

Effect of electric field
Within the framework of the weak coupling theory, we can further consider the effective of a vertical electric field. In the continuum model, turning on the electric field amounts to introducing a potential difference between the layers, As the time-reversal symmetry T remains unbroken under the electric field, the K′ valley Hamiltonian H K0 ¼ T H K T À1 is still related to that of the K valley H K by the time-reversal operation. We will focus on the band structure around the K valley. Figure 8 shows the effect of the electric field on the band structure and the Fermi surfaces for the cases of (a) V = 0.2v F |q a | and (b) V = 0.6v F | q a |. One can see that the Fermi surface is distorted as the electric field shifts the Dirac cones at K M and K 0 M relative to each other in energy (as they originated from the top and the bottom layers respectively).
We can follow the procedure described in Sec. IIA to extract the effect of the vertical electric field in the single-orbital model. However a symmetry analysis already suffices to determine the more relevant deformation of the Fermi surface. Given that the electric field breaks the M y :k + →k − mirror symmetry and preserves the C 3 :k + →e 2πi/3 k + rotational symmetry, new terms can be added to the single-orbital model Eq. (3) as The α′ and α″ terms describes the rotation and deformation of the Fermi surface as shown in Fig. 8 a for weak field. If the electric field is of the same order of the band width, the Fermi surface could be strongly deformed as in Fig. 8b, which goes beyond the perturbative description of Eq. (19). As a consequence of the Fermi surface deformation, the Fermi surface nesting between K and K′ valley will be suppressed by the electric field, therefore both the IVCW and the SC instability should reduce with the electric field. However as the deformation effect α″ enters at a higher order perturbation, one expects that nestingdriven orders remains insensitive to weak electric field, until the interlayer electric potential difference V reaches the order of the band width. Additional effects of interlayer electric field such as modulation of substrate effects due to the vertical displacement of the 2D electron gas can also play a role, and were not included in this analysis.
IVCW inter-valley coherence wave, TSC (inter-valley) topological superconductivity (d + id/p−ip-wave), s-SC (inter-valley) s-wave superconductivity. I 0 and Δ 0 are in the spin-singlet channel, I and Δ are in the spin-triplet channel

DISCUSSION
In summary, we presented a weak coupling analysis of valley fluctuation mediated superconductivity in twisted bilayer graphene. We started with a momentum space formalism of the lowenergy effective Moiré band structure, so as to circumvent the obstruction to constructing valley symmetric Wannier tight binding models. We identified the triangular (three-fold) anisotropy of the Fermi surface is a universal feature of the Moiré band structure around the charge neutrality, as it is the lowest-order distortion that is consistent with all the lattice symmetries. The Fermi surface anisotropy has important implications. The triangular shape of the Fermi surface allows a unique nesting between the parallel triangle sides of opposite valley Fermi pockets. This leads to enhanced valley fluctuations near half-filling, which in turn can provide the pairing glue which is demonstrated using the RPA. By solving the pairing gap equation with the RPA corrected interaction, we obtain the leading pairing instability in the intervalley channel with a d + id and p−ip mixed pairing form factor. The mixing between the d-wave and p-wave pairing is generic, because with triangular anisotropy, and the remaining C 3 symmetry, the angular momentum of the electron is only conserved modulo three, so there is no distinction between d + id and p−ip on symmetry ground. Further taking spin into account, one obtains both spin singlet and triplet chiral superconductors, parameterized by a fourvector n μ , where the μ = 0 component corresponds to the spinsinglet. Additionally, helical pairing orders were also discussed, parameterized by two orthogonal n μ vectors.
We emphasized the approximate SO(4) spin-valley symmetry. The naive SU(4) symmetry of four component electrons in valleyspin space is broken by the Fermi surface distortion which is opposite between the two valleys, leading to U(1) v × SO(4) symmetry. The SO(4) symmetry allows us to discuss the spinsinglet and spin-triplet pairings on equal footing, and n μ transforms as an SO (4) vector. This degeneracy is lifted by SO(4) breaking perturbations and we argued that a Hunds (anti-Hunds) interaction, i.e., an inter-valley ferromagnetic (antiferromagnetic) spin interaction favors spin-triplet (spin-singlet) pairing, which can be probed by studying the response to a Zeeman field. The presence of an approximate SO(4) symmetry could still have observable consequences which would be interesting to explore further. For example, if the SO(4) breaking is not too strong, a Zeeman field would tune a transition between singlet and triplet superconductors at low temperatures. Thus twisted bilayer graphene may provide an opportunity to study different SC phases and the phase transitions between them.
We propose two scenarios for the insulating phases. First, pushing to stronger interactions we see that inter valley coherence order can develop at the nesting vectors. The close commensurate wavevectors are the three M points corresponding to (0, π), (π, 0) and (π, π) at the midpoint of the triangular lattice Brillouin Zone edges. The simultaneous condensation of IVC order at these three wavevectors leads to a Slater insulator, although in our model a full gap obtains slightly below half filling at f ¼ À 1 2 À 1 8 . Future work should establish if a more complete treatment of interactions changes this conclusion. Nevertheless, other aspects of the phenomenology appear promising. For example, on hole doping the IVCW insulator, a single Fermi pocket appears, with two fold degeneracy (see Fig. 6). This agrees with the observed quantum oscillation experiments on the hole doped side of the Mott insulator below neutrality, where a Landau fan degeneracy in multiples of two was observed. As with superconductivity, the SO(4) symmetry implies a degeneracy between spin singlet and spin triplet IVCW orders, the latter being a kind of spin density wave. The same Hunds (anti-Hunds) SO(4) breaking interaction also picks out the spin-triplet (spin-singlet) IVCW order.
In ref. 13 , superconductivity was found to coexist with the insulating phase, i.e., superconducting puddles form even at half filling and establish a phase coherent state at very low temperatures. Assuming the orders are not spatially segregated, this places constrains on the possible pairs of order parameters which are likely not to destroy each other immediately. 102,103 In fact, the spin singlet IVCW order parameter and the spin singlet TSC order parameter anticommute with each other, therefore they are allowed to coexist in general (although they may compete for Fermi surface density of state on the level of energetics). In contrast, a triplet IVCW order will serve as a pair breaking order with respect to the singlet TSC order, such that it will rapidly destroy superconductivity. As with superconductivity, a Zeeman field may stabilize the spin triplet IVCW at reduced temperatures, which would be interesting to explore in future experiments. Similarly, the spin triplet IVCW and superconducting orders are mutually compatible. The common origin of superconductivity and IVCW order implies that their transition temperatures should scale together if interactions are enhanced. Both orders should also be experimentally testable.
Finally, we have considered in detail topologically ordered Mott insulators arising from freezing the charge fluctuations in the candidate superconducting states. We show that condensing double-vortices in the spin-singlet chiral TSC leads to a chiral valley-spin liquid state in the Mott phase, where the time-reversal symmetry is broken spontaneously. Thus the valley degeneracy is lifted in the Mott insulator, consistent with the two-fold Landau level degeneracy in the quantum oscillation experiment. The coexistence of such an insulator with the superconductivity is also natural, as the two phases only differ by chargeon condensation, which can form puddles in the presence of inhomogeneity.
Our study already reveals a plethora of orders and their interrelations on the basis of approximate symmetries as well as quantum interference effects. Undoubtedly, this just the scratches the surface of an even richer set of exciting phenomena made possible in this new experimental platform.
We notice that several related works appear around the same time. ref. 104 focuses on the nesting among hotspots at the van Hove energy where the (CDW′, SDW′) and (singlet SC, triplet SC) orders correspond to the IVCW order I μ and the inter-valley SC order Δ μ in our notation. ref. 105 points out that a fluctuating O(n) vector order (with n > 2) is crucial in explaining the emergence of SC inside the insulating phase.

METHOD
The band structure is calculated by exact diagonalization of the continuum model Eq. (1), where 120 bands are kept in the calculation. Only the middle bands are taken to build the effective model. The symmetry analysis of the interaction is carried out in the Majorana basis. As the electron carries valley and spin degrees of freedom, the Majorana basis is 8 dimensional. We systematically classify all the 28 fermion bilinear operators by enumerating all antisymmetric 8 × 8 matrices. Each representation is combined with its conjugate representation to form symmetric fourfermion interactions. We collect all interaction terms constructed in this way and simplify them to the form of Eq. (4).
The self-consistent mean-field calculation is carried out based on the singleorbital model. Given the reflection and rotation symmetry, we only need to solve the mean-field Hamiltonian in Eq. (12) a triangle region 0; 0 ð ÞÀðL; 0Þ À ðΛ=2; ffiffi ffi 3 p Λ=2Þ with the momentum cut-off Λ = 1.7. The order parameters in other regions in the momentum space are related by symmetry transformations. Within the triangular region, we take a momentum grid of 72 momentum points and solve the mean-field Hamiltonian at each point. The mean-field saddle point in found by iteratively update mean-field expectation values. The parameters are taken to be g I = 0.8, g Δ = 0.56. To increase the stability of the mean-field iteration, we adopt a trick to blur the pairing order parameter Δ 0 k in the momentum space by a Lorentzian kernel with a characteristic radius of 3 units on the momentum grid.

DATA AVAILABILITY
The dataset generated during and/or analyzed during the current study are available in the GitHub repository https://github.com/EverettYou/WeakCouplingTBG.