Unconventional Superconductivity and Density Waves in Twisted Bilayer Graphene

We study electronic ordering instabilities of twisted bilayer graphene with $n=2$ electrons per supercell, where correlated insulator state and superconductivity are recently observed. Motivated by the Fermi surface nesting and the proximity to Van Hove singularity, we introduce a hot-spot model to study the effect of various electron interactions systematically. Using renormalization group method, we find $d$/$p$-wave superconductivity and charge/spin density wave emerge as the two types of leading instabilities driven by Coulomb repulsion. The density wave state has a gapped energy spectrum at $n=2$ and yields a single doubly-degenerate pocket upon doping to $n>2$. The intertwinement of density wave and superconductivity and the quasiparticle spectrum in the density wave state are consistent with experimental observations.


I. INTRODUCTION
Recently superconductivity was discovered near a correlated insulator state in bilayer graphene with a small twist angle θ ≈ 1.1 • [1,2], where the moiré pattern creates a superlattice with a periodicity of about 13 nm. A correlated insulating state is found at the filling of n = 2 electrons per supercell (n = 0 is the charge neutrality point). Electron or hole doping away from n = 2 by electrostatic gating leads to a superconducting dome, similar to cuprates. Insulating states are also found in trilayer graphene with moiré superlattice [3]. The nature of superconducting and insulating states in graphene superlattices are now under intensive theoretical study [4][5][6][7][8][9][10][11][12][13].
For θ ≈ 1.1 • , the low-energy mini-band of twisted bilayer graphene has a narrow bandwidth of 10 meV scale [14][15][16][17][18][19][20][21][22][23][24][25][26]. However, this energy scale is still much larger than the energy gaps of the superconducting and insulating states, which are on the order of 1 K. Moreover, resistivity shows metallic behavior above 4 K. This is rather different from the case of a Mott insulator in the strong coupling limit, which would become insulating at much higher temperature. Based on these considerations, in this work, we take a weak coupling approach to study ordered states driven by electron correlation in twisted bilayer graphene.
While details of the band structure remain to be fully sorted out, a number of prominent features of the normal state fermiology are robust and noteworthy. First, at small twist angle, the two valleys of graphene have negligibly small single-particle hybridization and give rise to two separate Fermi surfaces that intersect each other in the mini Brillouin zone [26][27][28]. Second, as the carrier density increases, Fermi pockets associated with a given valley first appear around the Dirac points at K and K in the mini Brillouin zone, then these K and K pockets merge at a saddle point on the Γ-M line to become a single pocket centered at Γ. The saddle point associated with this Lifshitz transition has a Van Hove singularity (VHS) with a logarithmic divergence of the density of states (DOS). Third, realistic band structure calculations [28] show that near the Van Hove energy the Fermi surfaces of different valleys contain nearly parallel segments When the Fermi energy crosses the Van Hove energy, a conversion between electron and hole carriers is expected and indeed observed from the sign change of Hall resistance as a function of doping for θ = 2 • [28] and θ = 1.8 • [27]. Remarkably, this sign change occurs at the filling n = 2, indicating that the Fermi energy is very close to VHS. VHS in twisted graphene layers was also detected from pronounced peaks in DOS in STM measurements [29]. As the twist angle becomes smaller, the energy separation between VHS in conduction and valence bands is found to decrease rapidly from 430 meV at θ = 3.4 • to 82 meV at 1.79 • and 12 meV at 1.16 • . The strong reduction of bandwidth is expected to magnify the effects of electron correlation. This understanding is consistent with the fact that correlated insulator and superconducting states are found for θ = 1.1 • , but not for θ = 1.8 • and 2 • . In the following, we take the Fermi surface with good nesting condition and in proximity to VHS as a starting point and study its instabilities in the presence of electron interactions.
Due to the divergent DOS at VHS, electron interaction predominates in patches of the Brillouin zone around saddle points, or "hot spots". When multiple hot spots are present at a given energy, various scattering processes among them may interfere with each other, leading to intertwined density wave and superconducting instabilities. Such hot-spot models were studied with renormalization group (RG) approach in the context of cuprates [30][31][32][33][34], and recently, by Nandkishore, Levitov, and Chubukov in the context of doped monolayer graphene [35].
In this paper, we study interaction-driven ordering instabilities of twisted bilayer graphene around the filling n = 2 using RG by patching the Brillouin zone where the DOS is considerably larger than other parts. Our RG analysis shows how the electron interaction changes as the energy scale is reduced. Nontrivial RG flows of the coupling constants are found as a consequence of the nesting of Fermi surfaces. Susceptibility calculations reveal the possibility of various superconducting and spin/charge density-wave states at low temperature. When Coulomb repulsion is the dominant interaction, d/p-wave superconductivity and charge/spin density wave at a particular nesting wave vector emerge as two leading instabilities. The density-wave state is found to have a gapped energy spectrum at n = 2 and yields a single doubly-degenerate pocket upon doping to n > 2.

II. MODEL
We set up a model to analyze the electron interaction effect in twisted bilayer graphene around the filling n = 2. Our analysis of the interaction-driven instabilities focuses on hot spots, which dominates in the DOS. The hot spots are patches on the Brillouin zone. The Fermi surface nesting in the hot spots may potentially lead to ordering instabilities. In the following, we introduce a hot-spot model for twisted bilayer graphene and then the notion of Fermi surface nesting in hot spots.

A. Hot-spot model
To consider the electron interaction effect and resultant ordering instabilities near the filling n = 2, we focus on hot spots in the Brillouin zone which possess significantly larger electronic spectral weights compared to the other region. Such hot spots are obtained by patching the Brillouin zone around the saddle points. The patches are labeled by τ = ±1 for Fermi surfaces from the two valleys, σ for spins, and α = 1, . . . , 3 for the patches of a given valley (Fig. 1). We denote the position of a patch center by k ατ . The three inequivalent wave vectors connecting the patch centers are defined by Electron interaction is treated as scattering among the patches. By analogy with the g-ology model in onedimensional physics [36][37][38] and Shankar's RG approach [39], we write down the general interaction Hamiltonian compatible with lattice rotational symmetry Here the patch indices satisfy α 1 = α 3 = α 2 = α 4 (i = 1), , and α 1 = α 2 = α 3 = α 4 (i = 4). The valley indices τ 1 , . . . , τ 4 obey the same rule, associated with j. This rule is diagrammatically shown in Fig. 2(a). The interaction describes sixteen independent scattering processes with coupling constants g ij In the following analysis, we consider the momentum-conserving processes depicted in Fig. 2(b). Scattering processes related to these ones by lattice symmetry are not shown. Note that g i3 , g 12 , g 21 , and g 34 do not generally conserve crystal momentum since the patches are located away from the Brillouin zone boundary (see Fig. 12). Umklapp processes are allowed only when the hot spots are located at special momenta. The analysis for that case includes more or all g ij , and is presented in Appendix F. Among the nine momentum-conserving terms, g 11 , g 14 , g 22 , g 24 , g 44 are associated with forward scattering processes, and g 31 , g 32 , g 41 , g 42 are BCS scattering processes involving two electrons with opposite momenta. The nine momentum-conserving terms g ij in the interaction Hamiltonian Eq. (1) can be divided into three groups with different index j. For g i4 terms, all four operators belong to one valley (τ 1 = τ 2 = τ 3 = τ 4 ), thus describing intravalley interactions. g i2 terms are the product of two spin-and valley-conserving fermion bilinear operators that are associated with two different valleys; we call them intervalley density interactions. g i1 terms are the product of two spin-conserving but valley-flipping fermion bilinear operators; we call them intervalley exchange interactions.
We now discuss the microscopic origin of these interaction terms and how the coupling constants g ij should in principle be determined. First, when projected to the lowest mini-band, the long-range part of Coulomb interaction generates intra-and intervalley density interaction g i2 , g i4 , while the short-range part of Coulomb interaction on graphene's lattice scale generates intervalley exchange interaction g i1 involving large momentum transfer. Since Wannier functions in twisted bilayer graphene extend over tens of nanometers, the long-range part of Coulomb interaction is expected to dominate. Based on this factor alone, one would expect density interactions g i2 , g i4 to be orders-of-magnitude larger than the intervalley exchange interaction g i1 . On the other hand, it should be noted that the long-range Coulomb interaction strength is reduced by screening from excited bands that span a wide range of energies from ∼ 10 meV up to the bandwidth of graphene layers ∼ 10 eV. The process of integrating out these excited bands as well as those states of the lowest band outside the patches will significantly renormalize the coupling constants to be used in our patch theory. Moreover, their values are also affected by electron-phonon coupling. Since typical phonon energy in graphene is much larger than the mini-band width, it is reasonable to integrate out the phonons to obtain phonon-mediated electron-electron attraction [10], which renormalizes the values of coupling constants.
In this work, we treat the bare values of g ij as phenomenological parameters and calculate flows of these coupling constants under RG to the one-loop order. Strictly speaking, such a perturbative RG analysis is only legitimate for weak coupling. However, instabilities toward superconductivity and/or density waves are found within the weak-coupling regime (see below), which jus-tifies the one-loop RG analysis.

B. Susceptibilities and Fermi surface nesting
Fermion loops in the RG calculation are associated with bare susceptibilities in the particle-hole and particleparticle channels: where q ± = k ατ ± k βτ , f ( ) is the Fermi distribution, and τ k is the energy dispersion of valley τ with = 0 on the Fermi surface.
There are in total eight susceptibilities in the particleparticle and particle-hole channels at various wave vectors: where s = +, − correspond to intra-and intervalley components, respectively. Among these susceptibilities, χ 2+ (ω = 0) is the DOS within a patch at Fermi energy ρ 0 . χ 0− is the Q = 0 Cooper-pair susceptibility, which involves patches at opposite momenta, belonging to the different valleys. Regardless of Fermi surface geometry, χ 0− exhibits a logarithmic divergence χ 0− (ω) = ρ0 4 ln Λ ω in the presence of time-reversal symmetry, where Λ is the high-energy cutoff and depends on the patch size.
Besides the BCS channel, susceptibilities in other channels may find divergences when Fermi surfaces are perfectly nested: τ k+kατ = − τ k+k βτ for the particle-hole channels and τ k+kατ = τ −k+k βτ for the particle-particle channels. In general, the interplay between BCS and nesting-related interactions can lead to nontrivial RG flows of coupling constants [39,40].
As shown in Fig. 1, the Fermi surfaces associated with different valleys have nearly parallel segments connected by the following nesting vectors: Q − and Q connect occupied states of one valley and unoccupied states of another, while Q + connects occupied states around K and those around K associated with the same valley; see Fig. 1(d). Such Fermi surface nesting strongly enhances three types of bare susceptibilities: intervalley charge or spin density wave susceptibility at Q − and Q (particlehole channel), and intervalley pair density wave at Q + (particle-particle channel).
In the ideal case where the Fermi surface is perfectly nested (see Appendix A for the detailed discussion), these FIG. 3: One-loop corrections to the coupling constants. The solid lines represents the fermion propagators and the wavy lines correspond to interactions. The leftmost diagram involves a particle-particle loop, and the other three diagrams have particle-particle loops.
susceptibilities are logarithmically divergent like the BCS channel. When the Fermi surface is nearly nested, the logarithmic frequency energy dependence still holds approximately within a range Λ 0 < ω ≤ Λ, but the divergence in the ω → 0 limit is cutoff below a smaller energy scale Λ 0 associated with the deviation from perfect nesting.
Finally, when the VHS lies close to the Fermi surface, scatterings among states near VHS points receive particularly large RG corrections because of the large spectral weight. This justifies our use of patch RG approach. When the VHS lies exactly on the Fermi surface, the DOS at ω = 0 is logarithmically divergent and thus leads to an additional log divergence in susceptibilities, see discussion in Appendix B.

III. RG ANALYSIS
Loop corrections to the coupling constants suffer from divergences due to Fermi surface nesting and divergent DOS at the Van Hove energy, which are to be cured with the RG method. We consider corrections to one-loop order (Fig. 3). In the hot-spot model, each loop corrections is associated with the susceptibilities Eqs. (4)- (7). Since the Cooper-pair susceptibility χ 0− always gives the leading divergence regardless of the Fermi surface geometry, we set the RG scale by y ≡ χ 0− ( ). In the Wilsonian RG procedure, we integrate out high-energy modes and rescale the remaining low-energy modes as we increase the RG scale y, which is qualitatively similar to lowering the temperature T down from Λ. The other susceptibilities are measured with respect to y, parameterized by By definition, d 0− = 1 always holds.

A. RG equations
The RG equations for the coupling constants involve the parameters d as (y) defined in Eq. (8) as a function of the RG scale y. In general, d as (y) depends on y, except when the corresponding susceptibility χ as (y) diverges similarly to the BCS susceptibility. This occurs in the presence of Fermi surface nesting. Then the corresponding density-wave channel has a divergent susceptibility, and hence d as (y) = d as is a constant less than or equal to 1 [33][34][35].
In an ideal case for nesting where Fermi surface comprises of corner-sharing triangles (see Fig. 7), perfect nesting occurs simultaneously in three channels: three susceptibilities of the intervalley type, χ 1− (Q − ), χ 2− (Q ), and χ 3− (Q + ), are all logarithmically divergent similar to χ 0− , so that d 1− , d 2− and d 3− are nonzero constants. In contrast, none of the intravalley susceptibilities χ a+ is divergent. Thus d a+ (y) decay as y −1 (away from the Van Hove energy) or y −1/2 (at the Van Hove energy; see Appendix B), and hence become negligible at large y. We neglect the subleading terms d a+ in the following analysis. (RG equations for a generalized model with d a+ and additional interaction terms are presented in Appendix F.) As shown in Fig. 1, the Fermi surface of twisted bilayer graphene is nearly (but not perfectly) nested. In this case, the intervalley susceptibilities χ a− (a = 1, 2, 3) still have a logarithmic dependence on energy from Λ down to a smaller energy Λ 0 . Equivalently, the parameter d a− (y) is approximately constant within the corresponding range of the RG scale 0 ≤ y < y 0 . In the following, we shall analyze the RG flow within this energy range of interest, where the Fermi surface is regarded as nested.
For d a+ = 0, we obtain the RG equations for the nine momentum-conserving coupling constants as follows: . (15) We use the conventionġ ≡ dg/dy. Equations (9)- (15) show how different coupling constants change under RG. Among them, the intervalley interactions g 22 and g 11 involve two patches not related by time-reversal symmetry, hence receive corrections solely from scattering processes related to intervalley nesting. g 32 , g 42 , g 31 , and g 41 involve two patches related by timereversal symmetry, hence receive corrections from both BCS and nesting-related processes. The intravalley interactions g i4 do not flow because they do not participate in either process.
Details of the RG flow in the nine-dimensional parameter space are complicated and can in general be acquired numerically. (See Appendix D for the simplest case without nesting, where the analytic solution is obtained.) Nonetheless, its general feature can be understood easily: BCS corrections decrease repulsive interactions under RG, while nesting-related corrections tend to increase repulsive interactions in the corresponding channels. This important trend is a useful guideline to understand the behavior of the RG flow, which we present later.

B. Ordering instabilities
To analyze various possible instabilities, we consider susceptibilities associated with s-and d-wave spin-singlet superconductivity (s-SC etc.), p-and f -wave spin-triplet superconductivity, CDW, SDW, and pair density wave (PDW). Both p-and d-wave pairings have two degenerate components: (p x , p y ) and (d xy , d x 2 −y 2 ); see Appendix C for detail. Three different wave vectors, Q + , Q − , and Q associated with density-wave orders are indicated by superscripts +, −, and , e.g., SDW − and CDW .
When only the intervalley Fermi surface nesting is considered, the relevant instabilities are superconductivity, CDW/SDW at wave vectors Q − and Q , and PDW at Q + . An occurrence of an instability is indicated by a divergence of the corresponding susceptibility. The susceptibility are obtained by an RPA-like resummation [40] where η is used to label various ordering susceptibilities. χ 0 η (y) is the bare susceptibility in the absence of interaction and V η (y) is the effective interaction strength associated with the ordering.
By a straight-forward diagrammatic calculation, we find V η for various ordering channels as follows: Detailed description and derivation of interaction strengths and susceptibilities are found in Appendix C. When the parameters d as are constant, to the leading order in y, χ η (y) is written as The susceptibility diverges at 1 + V η (y c )χ η (y c ) = 0, leading to an instability at y c . An instability occurs only if the interaction strength is attractive: V η < 0. In meanfield theory, the interaction strength V η (y) is treated as a constant determined by the bare values of the coupling constants g ij . Then the instability temperature is given by T η = Λ exp[−(c p d as |V η |) −1/p ], for y = c p ln p (Λ/ ) (c p > 0) with p = 1 away from the VHS or p = 2 at the VHS. Our analysis here further takes into account the RG-scale dependence of coupling constants g ij and hence the interaction strength V η (y). We shall see that the running coupling constants leads to results beyond mean-field theory.

C. Intertwined superconductivity and density waves
Since the intervalley exchange interaction is likely smaller than the density-density interaction, it is instructive to first analyze cases only with the density-density interactions. For simplicity, we set the strengths of all density-density interactions equal (g 14 = g 24 = g 44 = g 22 = g 32 = g 42 > 0). In the absence of the exchange interactions, some susceptibilities become degenerate as we see from Eqs. (17)-(23): s-SC and f -SC, d-SC and p-SC, and CDW and SDW at each wave vector. Such degeneracy results from the fact that each valley has its own charge conservation and spin rotation symmetry. Similar degeneracy also occurs in exciton insulators when only long-range Coulomb interaction is considered [41].
With this choice of coupling constants, a mean-field analysis does not find any superconducting instability because the pairing interactions shown in Eqs. (17) and (18) are zero. In the presence of Fermi surface nesting, a density wave instability is found in mean-field analysis, whose wave vector is Q − if d 1− (g 22 + g 32 ) > d 2− g 42 , and Q vice versa.
Our RG analysis including the scale dependence of the coupling constants g ij finds qualitatively different results. Figures 4(a) and (b) are the phase diagrams from the one-loop RG analysis on the (d 1− , d 2− ) plane with d 3− = 0 and the (d 1− , d 3− ) plane with d 2− = 0, respectively. First, in both cases, d-or p-wave superconductivity is found for small d 1− , which is absent in mean-field theory. Second, the Q − density-wave state is far more dominant than the Q density-wave state: it already occurs at very small nesting parameter d 1− , even when the Fermi surface nesting condition is much stronger at wave vector Q .
To understand why the Q − density wave and d/p-wave superconductivity emerge as leading instabilities, we examine the flow of intervalley density interactions g 22 , g 32 , and g 42 . Since g 32 and g 42 are associated with BCS scattering processes, they receive renormalization even without nesting and decrease under RG when their initial values are repulsive, as shown in Figs. 4(c), (d). In contrast, since g 22 is associated with a forward scattering process, it is marginal without nesting. According to the RG equation Eq. (10), Fermi surface nesting in the   particle-hole channel (d 1− > 0) increases g 22 under RG. Therefore, in the presence of Fermi surface nesting, only g 22 grows without suppression from the BCS process and thus strongly enhances the Q − density wave fluctuation, making it dominate over the Q density wave. Although g 32 and g 42 both decrease under RG, the former decreases slower because BCS process and density wave nesting at wave vector Q − tend to renormalize g 32 in the opposite way; see Eq. (11). Therefore, a negative g 42 − g 32 < 0 is generated and its magnitude grows under RG. As shown in Eqs. (17) and (18), this attraction provides pairing interaction for both d-SC and p-SC and thus enhances these superconducting susceptibilities, see Fig. 4(f). The attractive pairing interaction should be stronger than that for the Q − density wave. Fermi surface nesting in the particle-particle channel (d 3− > 0) assists superconductivity in that it suppresses the increase of g 22 . Finite nesting in the particle-particle channel yields nonzero PDW susceptibility, but the interaction is repulsive for the PDW + fluctuation; see Eq. (23).
We conclude that in the presence of repulsive intervalley density interactions, the two leading instabilities are charge/spin density wave at wave vector Q − and p/dwave superconductivity. When the Fermi surface nesting in the particle-hole channel is strong (weak), the density wave state (superconductivity) is favored.

IV. ROLE OF INTERVALLEY EXCHANGE INTERACTION
The degeneracies of d/p-wave superconductivity and of CDW/SDW susceptibilities are lifted when intervalley exchange interactions g 11 , g 31 , g 41 are included. Their bare values depend on microscopic details as we have discussed in Sec. II A. For example, such interactions can arise from intervalley scattering mediated by optical phonons. Since the typical phonon frequency is much larger than the mini-band width, intervalley exchange interactions between low-energy electrons may be even attractive. Figure 5 shows the RG flows including both densitydensity and exchange interactions. With a small d 1− While the bare values of intervalley exchange interactions are hard to obtain accurately, we now discuss another important factor in selecting between d-and pwave SC, and between CDW − and SDW − . Both SDW and p-wave SC (which is spin-triplet) breaks the SU(2) spin rotational symmetry, while the CDW and d-wave SC (which is spin-singlet) do not. With SU(2) symmetry, it is known that thermal fluctuations associated with Goldstone modes prevent any true long-range spin order in two dimensions. This argument suggests that d-SC or CDW − can still be realized at nonzero temperature, even when the leading susceptibility above the ordering temperature is p-SC or SDW − .

V. ELECTRONIC STRUCTURE OF DENSITY-WAVE STATES
We now examine the electronic structure in a CDW − state and show that at filling the n = 2, the CDW − state can be insulating. The same conclusion applies to a collinear SDW − because it can be mapped to the CDW − state by performing a sign change on electrons of one spin polarization in one valley.
In density-wave states, the Brillouin zone in momentum space is reduced because of the enlarged unit cell in real space. The previously distinct momentum eigenstates now can hybridize, resulting in a band structure reconstruction. When the number of electrons in the enlarged unit cell is an even integer, the resulting CDW state can be a band insulator.
The CDW − order can occur at three equivalent wave vectors related by the C 3 rotational symmetry: with n φ = (cos φ, sin φ). Below we shall consider a triple-Q CDW state, where the above three wave vectors form the new reciprocal vectors, and hence define the reduced Brillouin zone. Compared to a single-Q state, the triple-Q state is expected to be energetically favorable as it gaps more parts of the Fermi surface especially around the hot spots with large DOS.
In the CDW − state, intervalley order parameter ψ † k+Qi,τ σ ψ kτ σ (τ = −τ ) becomes nonzero. The meanfield Hamiltonian in the CDW − state thus includes the CDW potential in addition to the original electron dispersion: Since the spin structure is irrelevant, we have dropped the spin index σ. Here we assume that the CDW order parameters at Q 1 , Q 2 , Q 3 are equal, so that the resulting state is invariant under the three-fold rotation. For the original Fermi surface shown in Fig. 1, the CDW − wave vector connecting a pair of hot spots is close to the commensurate vector Q − |ΓM|/2 = G/4, where G is the length of the original reciprocal lattice vectors of twisted bilayer graphene. (The analytic expression of the energy dispersion in the normal state τ k is given in Appendix E.) With this choice of CDW wave vector Q − , the reduced Brillouin zone is 4 × 4 smaller than the original Brillouin zone and can be constructed as shown in Fig. 6(a). Since there are two conduction bands (one per valley) in the original Brillouin zone, there are 32 bands in the reduced Brillouin zone. A complete filling of 16 bands corresponds to the filling of n = 2, where correlated insulating behavior was experimentally observed.
When the CDW − order parameter is small, the Fermi surface at the filling n = 2 is not fully gapped due to imperfect nesting. A full gap appears for ∆ ∆ c . For a realistic Fermi surface with good nesting condition, we find the critical value of the order parameter ∆ c = 0.15D, where D is the bandwidth of the original conduction band. The fact ∆ c D justifies our weak coupling approach.
The gapped energy spectrum in the CDW − state with ∆ = 0.16D is presented in Fig. 6(b). Importantly, we note that the direct gap in the CDW state is located at Γ in the reduced Brillouin zone. A single electron pocket (with two-fold spin degeneracy) is present above the gap, while two nearly degenerate hole pockets are present below the gap. The hole pockets have much heavier mass than the electron. These features are consistent with quantum oscillation measurements at densities slightly away from n = 2, as we shall discuss in the next section.
For the commensurate CDW state with Q − = G/4 considered here, the scattering process labeled by g 43 carries momentum 2Q = 4Q − = G, and thus it is allowed. This process corresponds to the intervalley exchange interaction and it is presumably smaller than intravalley and intervalley density interactions. We confirm that inclusion of small g 43 does not alter the RG flow much, and we obtain qualitatively the same result [42].

VI. DISCUSSIONS
In this section, we compare our results with the experiments on twisted bilayer graphene [1]. We have found from RG analysis the intertwining of unconventional superconductivity and density-wave instabilities. We have obtained from band structure calculations the gapped spectrum of density-wave states at the filling n = 2.
On the experiment side, the resistivity measurement at zero magnetic field near n = 2 observes a metallic behavior at high temperatures, then an upturn of resistivity in an intermediate temperature region, before superconductivity appears at the lowest temperature. Furthermore, the in-plane upper critical field of the superconducting state is found to be comparable to the Pauli limit, indicating spin-singlet pairing. The change from insulating to superconducting behaviors is consistent with the intertwined density wave and SC instabilities, shown by the evolution of susceptibility with decreasing energy scale in Figs. 4(e), (f) and also Figs. 5(e)-(f). Finally, when superconductivity is destroyed by the magnetic field, resistivity becomes insulating at the lowest temperature.
We interpret this T = 0 insulating state as a CDW/SDW state at wave vector Q − . We have analyzed a triple-Q − CDW/collinear SDW phase with 4 × 4 periodicity and have shown that a moderate density-wave order parameter fully gaps the energy spectrum at the filling n = 2, consistent with the insulating behavior of resistivity at low temperature. Importantly, at densities slightly above n = 2 (or doping towards complete filling of mini-bands), a single pocket with two-fold degeneracy is found in quantum oscillation measurements. This is consistent with our finding of a single pocket with spin degeneracy above the gap. On the other hand, at densities slightly below n = 2, quantum oscillations have so far not been observed. This is consistent with the fact that the pockets below the gap in our density-wave state have heavy mass.
Fermi surface nesting provides singularities in susceptibilities. The susceptibilities (Lindhard functions) in the particle-hole and particle-particle channels are given by τ and τ denote Fermi surfaces, which correspond to the valley degrees of freedom for the case of twisted bilayer graphene.
The conditions for perfect nesting is given by When the Fermi surfaces are perfectly nested, i.e., one of the above conditions holds in a certain area of the Brillouin zone, singularities in the susceptibilities are found in the static limit ω → 0. One finds a logarithmic divergence in a susceptibility with perfectly-nested Fermi surfaces in two dimensions, so that the susceptibility has ln(Λ/ω) dependence. In order to observe logarithmic dependence in susceptibilities at ω, the nesting condition should hold at the energy scale determined by ω. In other words, if the conditions are approximately met with an accuracy of around ω, we see ln(Λ/ω) behavior at the energy scale ω. It allows to relax the nesting conditions at ω to be for the particle-hole and particle-particle channels, respectively. When Fermi surfaces are perfectly nested, we have δ ph = 0 or δ pp = 0, and ln(Λ/ω) dependence in the susceptibility holds down to the lowest energies. On the other hand, when Fermi surfaces are nearly nested with δ ph/pp (ω) 1 for Λ 0 < ω ≤ Λ, we see a logarithmic enhancement with in the range Λ 0 < ω ≤ Λ, and it is cut off by the lower bound Λ 0 .
Our RG analysis focus on the energy range Λ 0 < ω ≤ Λ, where Fermi surfaces are regarded as nearly nested, and hence assume the parameters d as are constant within the range. For energy scale below Λ 0 , d as can no longer be regarded as constant, and we need to consider the y dependence of d as . Figure 7 shows simplified Fermi surfaces at the Van Hove energy, where the Fermi surfaces are approximated as corner-sharing triangles, to emphasize Fermi surface nesting. For convenience, we regard those Fermi surfaces as large hole Fermi surfaces encircling the Γ point. At a saddle point located along the Γ-M line, a Fermi surface touches another from an adjacent Brillouin zone. Therefore, there are two Fermi surfaces in one patch around a hot spot when we consider the reduced Brillouin zone. For convenience, we call the Fermi surface from the first Brillouin zone as the "inner" Fermi surface and the other from the second Brillouin zone as the "outer" Fermi surface.

Inner and outer Fermi surfaces
From the figure, we find four distinct Fermi surface nestings. First, the Fermi surfaces of patches 1 and 1 (1 and 1 are from different valleys) are perfectly nested in the particle-particle channel, yielding the BCS instability. It always holds for time-reversal-invariant systems. There are also symmetry-related pairs: 2-2 and 3-3 . Second, the inner Fermi surfaces are nested between 1 and 1 , leading to the logarithmic dependence in χ ph (Q , ω) and constant d 2− . In addition, by looking at both inner and outer Fermi surfaces, we can find that the Fermi surfaces from 1 and 2 are nested both in the particle-hole and particle-particle channels. It gives logarithmic dependences in χ ph (Q − , ω) and χ pp (Q + , ω), which results in constant d 1− and d 3− , respectively.
The discussion can be done in parallel near the Van Hove energy. When the Fermi energy is slightly above the Van Hove energy, we can regard the Fermi surfaces as hole pockets around the Γ point similarly as above. If we consider the filling slightly below the Van Hove energy, we should instead look at electron pockets, which surround the K and K points; cf. Fig. 1(c). In this case, a hot spot involves two Fermi surfaces from the two electron pockets. Still, a hot spot contains two Fermi surfaces, and Fermi surface nesting for both of them should be considered.

Appendix B: Susceptibilities near the Van Hove energy
The DOS logarithmically diverges at a saddle point in two dimensions because of the Van Hove singularity. When we approximate the energy dispersion near a saddle point as (δk) = Aδk 2 x − Bδk 2 y , we obtain the DOS (in the vicinity of the saddle point) as where Λ is the high-energy cutoff, which corresponds to the patch size in the present case. We have defined the eight susceptibilities χ as in the particle-hole and particle-particle channels with four different wave vectors in Eqs. (4)- (7). Among the eight susceptibilities, χ 2+ corresponds to the DOS; χ 2+ (ω) = ρ 0 (ω). The divergence in the limit ω → 0 is relaxed by finite chemical potential or away from the Van Hove energy: .
The susceptibility in the BCS channel χ 0− suffers from the Cooper instability in the presence of time-reversal symmetry. An alternative way to state this is that the Fermi surfaces are perfectly nested in this channel. We note that χ 0− is an intervalley susceptibility since the saddle points at k ατ and −k ατ (= k α,τ =τ ) belong to different valleys. Calculating χ 0− (ω), we obtain For µ = 0 (at the Van Hove energy), χ 0− has a double log singularity. It is the leading divergence in the eight susceptibilities and it sets the RG scale y. It is important to note that the two logarithms have different origins: one logarithmic singularity [ln(Λ/ max{ω, µ})] is from the Van Hove singularity and the other [ln(Λ/ω)] from the Fermi surface nesting. The logarithmic divergence from Fermi surface nesting has already been discussed. The divergence of the DOS is suppressed away from the Van Hove energy. However, if µ Λ is satisfied, there still are large spectral weights in the patches around the hot spots, which justifies the use of patch RG.
When Fermi surfaces are nested for a susceptibility in a certain channel and a wave vector (say, χ as ), a logarithmic dependence from Fermi surface nesting is present. (One can find an explicit calculation for the case of monolayer graphene in e.g. Ref. [43].) Then, χ as has the same singularity as χ 0− and the corresponding parameter d as (y) = dχ as /dy becomes constant as we define the RG scale by y ≡ χ pp (0, ). In contrast, Fermi surface nesting is absent for χ as , the corresponding parameter d as (y) decays as y −1/2 for y 1 [33,35].
Appendix C: Susceptibilities for ordering instabilities

Interaction strengths for instabilities
Now we consider an interaction strength that may trigger an ordering instability. We write two-particle interactions in the particle-particle and particle-hole channels asV One can perform a mean-field analysis by considering a mean field composed of the first or last two fermion operators on the right-hand sides. Those two types of interactions are diagrammatically depicted in Fig. 8(a). In the present model, there are 16 distinct scattering processes, and to the lowest order, interaction strengths for ordering instabilities are written as linear combinations of those coupling constants [ Fig. 8(b)]. Various interaction strengths V η corresponding to ordering η with all 16 coupling constants are given as follows: V CDW + =4(g 12 + g 14 + g 32 + g 34 ) − 2(g 21 + g 24 + g 31 + g 34 ), V CDW − =4(g 11 + g 13 + g 31 + g 33 ) − 2(g 22 + g 23 + g 32 + g 33 ), V CDW =2(2g 41 − g 42 + g 43 ) − 2(g 12 + g 13 ) + 4(g 21 + g 23 ), V SDW + = −2(g 21 + g 24 + g 31 + g 34 ), (C10) V SDW − = −2(g 22 + g 23 + g 32 + g 33 ), (C11) V SDW = −2(g 12 + g 13 + g 42 + g 43 ), (C12) V s = −(g 11 + g 14 + g 41 + g 44 ). (C17) In addition to instabilities for superconductivity and density waves, we also write down the interaction strengths for the charge compressibility V c and the uniform spin susceptibility V s . As for superconductivity, we can consider six different pairing symmetries because of the six patches. The six pairing symmetries are distinct in regard to the superconducting order parameters ∆ at the patches (Fig. 9).
Each p-and d-wave pairing has two different pairings, but those two give the same interaction strength V η . We assume that the system respects the point group D 3 , which has the representations A 1 , A 2 , and E. From this viewpoint, s-and f -wave pairings belong to the onedimensional representations A 1 and A 2 , respectively, and each p-and d-wave pairing belong to the two-dimensional representation E. For the other interaction strengths, we assume the A 1 representation.

Derivation of susceptibilities for ordering instabilities
We calculate the susceptibilities for ordering η χ η by summing up RPA-like diagrams, shown in Fig. 8(c); see also Ref. [40]. It yields the susceptibility χ η in the form of The RG equations (9)-(15) involve many variables and parameters, and a full analysis of all nine running coupling constants is rather complicated. Here we analyze the simplest case where there is no Fermi surface nesting other than the BCS channel (d as = 0, except for d 1− = 1 by definition). Then, only g 31 , g 32 , g 41 , and g 42 receive one-loop corrections since they can be regarded as BCS interactions. The RG equations for the four coupling constants are written as where we define g a± = g a2 ± g a1 (a = 3, 4).
We note that g a+ and g a− characterize spin-singlet and spin-triplet SC, respectively. The RG flow for g 4+ and g 3+ is shown in Fig. 10. The coupled RG equations can be rewritten as and thus are readily solved. When we neglect the exchange interactions, g 4± and g 3± are reduced to be g 42 and g 32 , respectively. If we further assume that those two intervalley density-density interactions have the same amplitude g 42 = g 32 , the equality holds under the RG flow from Eq. (D5). Each coupling constant obeys the RG equationġ = −3g 2 , which is in accordance with Shankar's result for BCS interactions without nesting [39].

Appendix E: Energy dispersion of a normal state
We construct an analytic form of the energy dispersion with two valley degrees of freedom, belonging to the point group D 3 and satisfying a filling condition. For the latter, we require that the energy dispersion have VHS points at the filling n = 2. The symmetry conditions are given as follows: The C 3 rotational symmetry is kept within a valley, and hence the energy dispersion satisfies the relation There also exist in-plane C 2 rotations about the Γ-K lines, where one of them is taken along the k y axis. The in-plane C 2 rotation interchanges the two valley and flips k x , thus requiring τ kx,ky = τ −kx,ky .
A function that obeys the symmetry conditions Eqs. (E1) and (E2) can be represented as a Fourier series because of the underlying lattice periodicity. Up to the third order, the symmetry-allowed terms are given by where we define k i y = C i 3 k y and k i x = C i 3 k x . The parameters t 1 , t 2 , t 2 , and t 3 are real. In addition to the D 3 symmetry, we further impose the filling condition, which is fulfilled by tuning the parameters t 1 , t 2 , t 2 , and t 3 . By choosing e.g. t 1 = 1, t 2 = −1, t 2 = 0.3, and t 3 = 0.6, we find saddle points of the energy dispersion, i.e., VHS points, at the filling n = 2, as shown in Fig. 11. The energy bands extend in the range of 0 ≤ E/t 1 ≤ D = 7.2.
We generalize the model for twisted bilayer graphene, consisting of two valleys with SU(2) spin, to a (N v × N s )band model, consisting of N v valleys with SU(N s ) spins. We assume that the valley degrees of freedom are nondegenerate in the kinetic part, and as a result, there are N v separate bands each with N s -fold degeneracy from spin; i.e., the kinetic component has U(1) × SU(N s ) sym-metry. Each Fermi surface is assumed to possess n p hot spots. Such a model is applicable for example to cuprates, monolayer graphene, and graphene superlattices. Our results for graphene superlattices are obtained by setting N v = N s = 2 and n p = 3. A model for cuprates corresponds to N c = 1, N s = 2, and n p = 2 [30][31][32][33][34], and one for monolayer graphene to N c = 1, N s = 2, and n p = 3 [35]. All those cases are in two dimensions; a generalization to other dimensions is straightforward.
The orbital index takes τ = 1, . . . , N v , the spin index σ = 1, . . . , N s , and the saddle point index α = 1, . . . , n p . Redefining the range of the indices, we use the same interaction as Eq. (1). As long as we have the same set of coupling constants, the structure of the RG equations does not depend either on the position of the saddle points in the Brillouin zone or the dimensionality of the system, which merely changes the parameters d as . Note that the dimensionality constrains geometrically possible N v and n s if we consider the same set of coupling constants. Also, note that umklapp processes in the valley space (j = 3) explicitly violate the valley U(1) symmetry.

RG equations for the coupling constants
For the generalized (N v × N s )-band model where each Fermi surface has n p saddle points, we obtain the RG equations for the sixteen coupling constants g ij to oneloop level (Fig. 3): The RG equations for the case of twisted bilayer graphene are obtained with N v = N s = 2 and n p = 3. By setting N v = 1, N s = 2, we can reproduce the previous results for cuprates with n p = 2 [30][31][32][33][34] and for graphene with n p = 3 [35]. Since there are no valley degrees of freedom in such cases, we need only four coupling constants g i and four parameters d a , dropping the second subscripts j from g ij and s from d as .