Superconductivity and Charge Density Wave in Twisted Bilayer Graphene
We study electronic ordering instabilities of twisted bilayer graphene with two 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 intertwined spin-singlet (triplet) superconducting and charge (spin) density-wave orders, consistent with the experimental observation of resistivity changing from metallic to insulating and finally to superconducting behavior upon decreasing temperature.
Recently superconductivity was discovered near a correlated insulator state in bilayer graphene with a small twist angle exp1 (); exp2 (), where the moiré pattern creates a superlattice with a periodicity of about 13 nm. The insulating state is found at the filling of electrons per supercell ( is the charge neutrality point). Electron or hole doping away from by electrostatic gating leads to a superconducting dome, similar to cuprates. Mott insulating states are also found in trilayer graphene with moiré superlattice Wang (). The nature of superconducting and insulating states in graphene superlattices are now under intensive theoretical study Yuan (); Xu (); Volovik (); Po (); Zhang (); Kivelson (); Philips (); Baskaran (); Lee (); Yang ().
For , the low-energy mini-band of twisted bilayer graphene has a narrow bandwidth of 10 meV scale Nam (); TB1 (); TB2 (); TB3 (); TB4 (); TB5 (); DFT1 (); DFT2 (); Nori2 (); Neto (); Neto2 (); MacDonald (); Mele (), which provides the setting for superconducting and insulating states with gaps on the order of 1 K. While details of the band structure remain to be fully sorted out, several essential features of the fermiology appear to be robust and deserve attention. First, at small twist angle, the two valleys of graphene have negligibly small single-particle hybridization and thus give rise to two sets of intersecting Fermi surfaces MacDonald (); PRL2016 (); Kim (). Second, as the carrier density increases, Fermi pockets associated with a given valley first appear around the Dirac points at and in the mini Brillouin zone, then these and pockets merge at a saddle point on the – 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, as shown by realistic band structure calculations Kim (), near the Van Hove energy the Fermi surfaces of different valleys contain nearly parallel segments and hence are nearly nested, see Fig. 1. Such Fermi surface nesting strongly enhances density wave fluctuations.
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 Kim () and PRL2016 (). Remarkably, this sign change occurs at the filling , indicating that the Fermi energy is very close to VHS.
The VHS in twisted graphene layers was also detected from pronounced peaks in DOS in STM measurements Andrei (). 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 to 82 meV at and 12 meV at . The strong reduction of band width 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 , but not for and . We regard the metallic state at and featuring VHS near Fermi energy as an important guide to the normal state of the correlated insulator and superconductor at .
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 Schulz (); Dzyaloshinskii (); Lederer (); Furukawa (); LeHur (), and recently, by Nandkishore, Levitov, and Chubukov in the context of doped monolayer graphene Nandkishore ().
In this paper, we study interaction-driven ordering instabilities of twisted bilayer graphene around the filling 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.
The patches are labeled by for Fermi surfaces from the two valleys, for spins, and for the patches of a given valley (Fig. 1). We denote the position of a patch center by . The three inequivalent wave vectors connecting the patch centers are defined by (intravalley), and (intervalley) (, ), see Fig. 1(e). Electron interaction is treated as scattering among the patches. By analogy with the -ology model in one-dimensional physics 1d_1 (); 1d_2 () and Shankar’s RG approach Shankar (), we write down the general interaction Hamiltonian compatible with lattice rotational symmetry
Here the patch indices satisfy , , , and . The valley indices obey the same rule, associated with . This rule is diagrammatically shown in Fig. 2(a). The interaction describes sixteen independent scattering processes with coupling constants , sketched in Figs. 2(b) and (c). Note that , , , and do not conserve crystal momentum since the patches are located away from the Brillouin zone boundary. We exclude those seven momentum-nonconserving terms from the following analysis. The analysis including all is described in Supplemental Material (SM) SM (). Among the nine momentum-conserving terms, , , , , are forward scattering interactions, while , , , are BCS interactions involving two electrons with opposite momenta.
The bare values of coupling constants are given by the effective electron interaction at the cutoff energy scale of the patch theory, which in principle can be obtained by integrating out high-energy states outside the patches or in excited bands. In the following, we treat the bare values of as input 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 justifies the one-loop RG analysis.
Fermion loops in the RG calculation are associated with bare susceptibilities in the particle-hole and particle-particle channels:
where , is the Fermi distribution, and is the energy dispersion of valley with on the Fermi surface.
There are in total eight susceptibilities in the particle-particle and particle-hole channels at various wave vectors:
where correspond to intra- and intervalley components, respectively.
Among these susceptibilities, is the DOS within a patch at Fermi energy . is the Cooper-pair susceptibility, which involves patches at opposite momenta, belonging to different valleys. Regardless of Fermi surface geometry, exhibits a logarithmic divergence in the presence of time-reversal symmetry, where is the ultraviolet cutoff and depends on the patch size.
Besides the BCS channel, susceptibilities in other channels may find divergences when Fermi surfaces are perfectly nested: for the particle-particle channels and for the particle-particle channels. In general, the interplay between BCS and nesting-related interactions can lead to nontrivial RG flows of coupling constants Shankar (); Chubukov ().
As shown in Fig. 1, the Fermi surfaces associated with different valleys have nearly parallel segments connected by the following nesting vectors: and connect occupied states of one valley and unoccupied states of another, while connects occupied states around and those around associated with the same valley, see Fig. 1(e). Such Fermi surface nesting strongly enhances three types of bare susceptibilities associated with an intervalley charge or spin density wave (particle-hole channel) at and , and intervalley pair density wave (particle-particle channel) at respectively. In the ideal case where the Fermi surface is perfectly nested (see SM SM ()), these 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 , but the divergence in the limit is cutoff below a smaller energy scale 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 is logarithmically divergent and thus leads to an additional log divergence in susceptibilities, see discussion in SM SM ().
Since the Cooper-pair susceptibility always gives the leading divergence regardless of the Fermi surface geometry, we set the RG scale by . In the Wilsonian RG procedure, we integrate out high-energy modes and rescale the remaining low-energy modes as we increase the RG scale , which is qualitatively similar to lowering the temperature down from . The other susceptibilities are measured with respect to , parameterized by
With one-loop calculations, we obtain the RG equations for the momentum-conserving nine coupling constants
We use the convention .
In order to solve the RG equations for the coupling constants, we need to know the parameters defined in Eq. (7) as a function of the RG scale . In general, depends on , except when the corresponding susceptibility diverges logarithmically like the BCS susceptibility. This occurs in the presence of Fermi surface nesting. Then the corresponding density-wave channel has a divergent susceptibility, and hence is a constant less than or equal to Furukawa (); Nandkishore ().
In the ideal case where Fermi surface comprises of corner-sharing triangles (see SM SM ()), perfect nesting occurs simultaneously at three different wave vectors: , and . Therefore, three susceptibilities of intervalley type, , and , are all logarithmically divergent, so that and are nonzero constants. In contrast, the intravalley susceptibilities are not divergent, and hence decay as and become negligible at large . We neglect the subleading terms hereafter.
As shown in Fig. 1, the Fermi surface of twisted bilayer graphene is nearly (but not perfectly) nested. In this case, the intervalley susceptibilities () still have a logarithmic dependence on energy from down to a smaller energy . Equivalently, the parameter is approximately constant within the corresponding range of the RG scale . In the following, we shall analyze the RG flow within this energy range of interest, where the Fermi surface is regarded as nested.
The values of , , depend on the shape of the Fermi surfaces and are generally different due to the different lengths of nested Fermi surface segments at the corresponding wave vectors. For example, for the Fermi surface at an energy slightly above VHS shown in Figs. 1(c) and (d), the strongest nesting is expected to occur in the particle-hole channel at wave vector , making larger than and .
A full analysis of all nine running coupling constants is rather complicated. To gain insight, we find it useful to consider the scenario where certain interactions are dominant and drive strong density-wave fluctuation associated with one particular nesting wave vector. As an example, we study the case of intervalley density wave fluctuations at wave vector , driven by and . Setting the bare values of all other couplings to zero, we obtain two coupled RG equations involving and only:
Figure 3 shows the RG flows for and . For , when is satisfied, flows to large negative values, driving the BCS instability in either spin-singlet or triplet channel. When , the system flows to weak coupling. This phase diagram, shown in Fig. 3(a), can also be obtained from mean-field theory by noting that is the pairing interaction in the spin-singlet (triplet) channel. The absence of Fermi surface nesting rules out density-wave instabilities.
For , the RG flow of , is significantly modified, and the weak-coupling region shrinks [Fig. 3(b)]. This means that attractive pairing interaction can be generated at low energies from repulsive bare interaction, a result that goes beyond mean-field theory. Furthermore, the Fermi surface nesting allows density-wave instabilities in addition to the BCS instability. As a result, superconductivity and density-wave instabilities can be strongly intertwined and evolve nontrivially with decreasing temperature (see below).
To analyze various possible instabilities, we consider susceptibilities associated with - and -wave spin-singlet superconductivity (-SC etc.), - and -wave spin-triplet superconductivity, CDW, SDW, and pair density wave (PDW). Both - and -wave pairings are doubly degenerate SM (). Three different wave vectors, , , and associated with density-wave orders are indicated by superscripts , , and , e.g., SDW and CDW.
The interaction strengths for the instabilities are given by SM ()
where is the susceptibility in the absence of interaction. Here the subscript denotes various superconducting and density-wave channels. When the parameters are constant, to the leading order in , is written as
A susceptibility diverges at , leading to an instability. An instability occurs for . To the leading order in , the instability temperature is obtained as , for (), with away from the VHS or at the VHS.
We analyze the susceptibilities again for the simplest nontrivial case where only two coupling constants and are nonzero and the parameter characterizes the nesting condition at wave vector . In this case, there are four possible instabilities: spin-singlet SC, spin-triplet SC, CDW, and SDW. Figure 4(a) shows instabilities as a function of the nesting parameter and the bare values of , at . When is constant, only the ratio matters, and thus we parametrize the coupling constants by . All four possible instabilities are found in the phase diagram, characterized by a divergent susceptibility. We confirm that divergences of susceptibilities occur while the coupling constants remain small, which justifies the one-loop RG analysis SM ().
The phase diagram shows that singlet SC is adjacent to CDW, while triplet SC is adjacent to SDW. For parameters near the phase boundary between singlet SC and CDW instability, Fig. 4(b) shows a nontrivial evolution of susceptibilities with the RG scale, or equivalently decreasing temperature. Since the particle-hole channel has a weaker nesting than the BCS channel, the CDW susceptibility is initially smaller than the SC susceptibilities at . As the RG scale increases, the CDW susceptibility grows steadily, and in an intermediate region, exceeds SC susceptibility. However, at a later stage, the SC susceptibility overtakes CDW and diverges first.
This intertwined behavior of singlet SC and CDW instabilities can be understood from the RG flow of the corresponding interaction strength [Eq. (8)]: for the spin-singlet SC and for CDW respectively. Both share , and hence a negative amplifies the two fluctuations simultaneously, leading to the intertwinement rather than competition. On the other hand, the two interactions differ in the sign of . A positive favors CDW instead of SC, hence makes CDW the leading instability at an intermediate RG scale Fig. 4(b). Nonetheless, the RG flow in Fig. 3 reveals that eventually flows from being positive to negative, making SC susceptibility to diverge first.
This intertwined behavior between spin-singlet SC and CDW is made possible by repulsive and attractive . The repulsive interaction naturally arises from the Coulomb interaction between electrons, while the attractive interaction can arise from the intervalley electron-phonon coupling. Note that typical phonon energy in graphene is much higher than the mini-band width, hence it is reasonable to integrate out the phonons and work directly with phonon-mediated electron-electron interaction Kivelson ().
With the simplification of , different patches related by three-fold rotations are decoupled. This leads to a degeneracy between - and -wave pairings, and between - and -wave. RG analysis with four coupling constants , , , and can be found in SM SM (). Moreover, depending on Fermi surface shape, the other two nesting parameters and can become important in determining possible instabilities. With three nesting parameters , we also find that interaction can drive instabilities associated with intervalley CDW and SDW at wave vector and PDW at wave vector SM ().
We now compare our main result—the intertwinement of superconductivity and density-wave instabilities—with the experiment of twisted bilayer graphene exp1 (). Near filling, resistivity measurement at zero magnetic field 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 CDW and spin-singlet SC instabilities, shown by the evolution of susceptibility with decreasing energy scale in Fig. 4(b). Finally, when superconductivity is destroyed by the magnetic field, resistivity becomes insulating at the lowest temperature. We interpret this insulating state as a CDW state where three ordering wave vectors related by symmetry are simultaneously present. At strong coupling, such three- CDW order is expected to gap the original Fermi surface completely, leading to an insulating state. We leave a further analysis of CDW states in twisted bilayer graphene for the future.
Acknowledgment. We thank Pablo Jarillo-Herrero, Yuan Cao, Valla Fatemi, Oskar Vafek, Leonid Levitov, Sankar Das Sarma, Allan MacDonald and especially Eva Andrei for helpful discussions. This work is supported by the DOE Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under award DE-SC0010526. LF is partly supported by the David and Lucile Packard Foundation.
In Supplemental Material, we explain the notion of Fermi surface nesting used in the main text, discuss singularities at the Van Hove energy, derive the formula for susceptibilities with corresponding interaction strengths, show several numerical solutions to the RG equations, analyze the RG equations for some simplified cases, and present the RG equations for a generalized model.
Appendix A Fermi surface nesting
a.1 Perfect and near nesting
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 . One finds a logarithmic divergence in a susceptibility with perfectly-nested Fermi surfaces in two dimensions, so that the susceptibility has 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 behavior at the energy scale . It allows to relax the nesting conditions at to be
When Fermi surfaces are perfectly nested, we have or , and dependence in the susceptibility holds down to the lowest energies. On the other hand, when Fermi surfaces are nearly nested with for , we see a logarithmic enhancement with in the range , and it is cut off by the lower bound .
Our RG analysis focus on the energy range , where Fermi surfaces are regarded as nearly nested, and hence assume the parameters are constant within the range. For energy scale below , can no longer be regarded as constant, and we need to consider the dependence of .
a.2 Inner and outer Fermi surfaces
When we consider Fermi surface nesting near a saddle point where multiple Fermi surfaces touch, we need to consider nesting of all those Fermi surfaces. For simplicity, we discuss here the Fermi surface at the Van Hove energy, but the discussion can be done in parallel even near the Van Hove energy. In the model of twisted bilayer graphene, focusing on the Fermi surface of one valley, we find one large hole Fermi surface encircling the point. At a saddle point (along the – line), it touches another Fermi surface 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.
In the main part, we discuss the nesting with wave vector , which involves the nesting of the inner Fermi surfaces from the different valleys. It leads to a singularity in and hence constant within the energy range for the nearly-nested Fermi surfaces. In addition, we can consider Fermi surface nesting involving both inner and outer Fermi surfaces. Figure S1 shows simplified Fermi surfaces, where the Fermi surfaces are approximated as corner-sharing triangles, to emphasize Fermi surface nesting. From the figure, we find four distinct Fermi surface nestings. First, the Fermi surfaces of patches and ( and 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: – and –. Second, as we have already mentioned, the inner Fermi surfaces are nested between and , leading to the logarithmic dependence in and constant . In addition, by looking at both inner and outer Fermi surfaces, we can find that the Fermi surfaces from and are nested both in the particle-hole and particle-particle channels. It gives logarithmic dependences in and , which results in constant and , respectively.
Appendix B Susceptibilities near the Van Hove energy
The density of states 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 , we obtain the density of states (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 eight susceptibilities in the particle-hole and particle-particle channels with four different wave vectors (see the main text). Among the eight susceptibilities, corresponds to the density of states; . The divergence in the limit is relaxed by finite chemical potential or away from the Van Hove energy:
The susceptibility in the BCS channel 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 is an intervalley susceptibility since the saddle points at and belong to different valleys. Calculating , we obtain
For (at the Van Hove energy), has a double log singularity. It is the leading divergence in the eight susceptibilities, and it sets the RG scale . It is important to note that the two logarithms have different origins: one logarithmic singularity is from the Van Hove singularity and the other from the Fermi surface nesting. The logarithmic divergence from Fermi surface nesting has already been discussed. The divergence of the density of states 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, ), 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. s_Gonzalez ().) Then, has the same singularity as , and the corresponding parameter becomes constant as we define the RG scale by . In contrast, Fermi surface nesting is absent for , the corresponding parameter decays as for s_Furukawa (); s_Nandkishore ().
Appendix C Susceptibilities for ordering instabilities
c.1 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 as
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. S2(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; see Fig. S2(b). Various interaction strengths corresponding to ordering with all 16 coupling constants are as follows:
In addition to instabilities for superconductivity and density waves, we also write down the interaction strengths for the charge compressibility and the uniform spin susceptibility .
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. S3). Each - and -wave pairing has two different pairings, but those two give the same interaction strength . We assume that the system respects the point group , which has the representations , , and . From this viewpoint, - and -wave pairings belong to the one-dimensional representations and , respectively, and each - and -wave pairing belong to the two-dimensional representation . For the other interaction strengths, we assume the representation.
c.2 Derivation of susceptibilities for ordering instabilities
where is the susceptibility calculated without interaction. are equal to with depending on : (PDW), (all SC), (CDW, SDW), (CDW, SDW), (CDW, SDW), (PDW), and (PDW).
Appendix D Various instabilities and validity of the one-loop RG
Depending on the bare values of the coupling constants, instabilities occurring as the RG scale grows (i.e., in low temperatures) varies. For all numerical results in this paper, we fix the parameters , which characterize the degree of Fermi surface nesting, to be , , and .
First, we show the results for the case where only two coupling constants and are nonzero. This case has already been considered in the main part, and the results are partly shown in Fig. 4. In SM, we show the RG flows of the coupling constants and susceptibilities for all four possible instabilities (Fig. S4), in order to check the validity of the one-loop RG calculation. For all four cases, we confirm that the divergences of susceptibilities happen within the region where the coupling constants remain small. This observation justifies the perturbative one-loop RG calculation.
Finite and in addition to and lift the degeneracies of two pairing symmetries each belonging to spin-singlet and spin-triplet superconductivities. Figure S5 displays the four different instabilities toward superconductivity with different pairing symmetries. As we see later from analytic expressions, the tendency to superconductivity is enhanced by including and , since those two scattering processes belong to the BCS interactions. We also confirm that the divergences of susceptibilities occur within the weak-coupling regime.
With nonzero , , and , it is also possible to find CDW, SDW (at wave vector ), and PDW (at wave vector ). Those can be realized by considering , , , and . We note that all those four coupling constants can be induced by an on-site density-density interaction. Figure S6 shows that CDW, SDW, and PDW instabilities are indeed possible depending on the values of the coupling constants. Those instabilities are again confirmed within the weak-coupling regime.
Appendix E Analysis of the RG equations for simplified cases
e.1 Nonzero and
The coupling constants and distinguish the pairing symmetries of superconductivity, but superconductivity itself depends largely on and , including whether it is spin-singlet or spin-triplet pairing. It is suggestive to point out that those two coupling constants can be regarded simultaneously as corrections to the Landau parameters (zero sound) and the BCS-type couplings.
If we neglect other coupling constants than and , the RG equations for the two coupling constants are
These two coupled RG equations can be rewritten as two decoupled RG equations:
The RG flows on the – plane are shown in Fig. 3 (in the main text) by changing the parameter . These decoupled first-order differential equations can be readily solved analytically, leading to the analytic expressions of the two coupling constants and susceptibilities.