Mean-field theory for superconductivity in twisted bilayer graphene

Recent experiments show how a bilayer graphene twisted around a certain magic angle becomes superconducting as it is doped into a region with approximate flat bands. We investigate the mean-field $s$-wave superconducting state in such a system and show how the state evolves as the twist angle is tuned, and as a function of the doping level. We argue that part of the experimental findings could well be understood to result from an attractive electron--electron interaction mediated by electron--phonon coupling, but the flat-band nature of the excitation spectrum makes also superconductivity quite unusual. For example, as the flat-band states are highly localized around certain spots in the structure, also the superconducting order parameter becomes strongly inhomogeneous.

Experiments on strongly doped graphene [1][2][3][4] have shown that with proper preparations, graphene can be driven to the superconducting state. Such experiments indicate that the lack of superconductivity in undoped graphene is not necessarily due to a lack of an (effective) attractive electron-electron interaction with strength λ that would drive graphene superconducting, but rather the small density of states (DOS) close to the Dirac point. Technically, in contrast to the Cooper instability for metals taking place with arbitrarily small λ, superconductivity in an electron system with a massless Dirac dispersion 2 p = v 2 F p 2 and an energy cutoff c has a quantum critical point λ c = π 2 v 2 F /(2 c ) [5] such that for λ < λ c , meanfield superconductivity does not show up at any temperature. From this perspective, doping to a potential µ leads to an increased DOS, and thereby to a non-vanishing critical temperature T c ≈ |µ| exp[−(λ c /λ − 1) c /|µ| − 1]. An alternative approach would be to change the spectrum and increase the density of states close to the Dirac point. The extreme limit would be an approximatively flat band of size Ω FB , where the group velocity tends to zero. In such systems the critical temperature is a linear function of the coupling strength T c = λΩ FB /π 2 [6] and quite high T c can be expected even without extra doping [7][8][9][10][11][12][13].
Recent observations of superconductivity in twisted bilayer graphene [TBG, see Fig. 1(a)] [14] take place in systems where theoretical studies have predicted the occurrence of asymptotically flat bands [15][16][17][18][19][20][21][22][23][24][25]. Here we study the mean field theory of superconductivity in such systems, starting from the hypothesis [26] that the observations could be explained from the flat-band perspective. In particular, we use the model of Refs. [15,20] for the spectrum of the twisted bilayer, add an on-site attractive interaction of strength λ, and evaluate the mean-field order parameter profile.
We find that the order parameter, and along with it the mean-field critical temperature, have a similar nonmonotonous behavior with respect to the twist angle as in the experiments. We also predict the behavior of the density of states in the superconducting state, resulting from the peculiarities of the flat-band eigenstates and from the position dependence of the superconducting gap [ Fig. 1(b)]. Even if our pairing interaction is quite simple, the resulting energy dependent density of states is quite unusual. In addition, we show how doping away from the flat band eventually destroys superconductivity.
Normal state. We describe the normal state of TBG with the model of Refs. [15,20]. With this model, we can describe the twist angles θ at which the lattices L and L θ of the two graphene layers are commensurate, so that the system as a whole is periodic in the moiré superlattice SL. For small θ, the set of all twist angles is dense in commensurate angles, in which case we can approximate an incommensurate structure by a commensurate one [20]. Hence for simplicity we study only the structures for which the exact superlattice matches the approximate superlattice. These structures are characterized by a rotation parameter m ∈ N, for which the rotation angle is given by cos(θ) = 3m 2 + 3m + 1/2 3m 2 + 3m + 1 . (1) The primitive vectors of the superlattice SL are given by t 1 = ma 1 +(m+1)a 2 , t 2 = −(m+1)a 1 +(2m+1)a 2 and the primitive vectors of the reciprocal superlattice SL * are G 1 = 4π 3||t1|| 2 ((3m + 1)a 1 + a 2 ), G 2 = 4π 3||t1|| 2 (−(3m + 2)a 1 + (3m + 1)a 2 ), where the lattice constant of the superlattice is ||t 1 || = √ 3m 2 + 3m + 1 a and the graphene lattice primitive vectors are a 1 = (1, √ 3)a/2 and a 2 = (−1, √ 3)a/2 with a the lattice constant [15]. The first Brillouin zone of the TBG is denoted by R 2 /SL * . In the following, we assume that G ∈ SL * and k ∈ R 2 /SL * , and also that the sums and integrals are restricted to these sets.
In the normal state, TBG is described by a low energy effective Hamiltonian [15] H ρk (G, G ) = (2) [ where the matrix structure corresponds to the layer structure and ρ ∈ {+, −} is the valley index with + corresponding to K and − to K = −K. Furthermore, each entry is a 2 × 2 matrix due to the sublattice structure in graphene. The diagonal terms in Eq. (2) describe the Dirac dispersion in the two layers and are diagonal also in G. Here, σ ρ = (ρσ x , σ y ). For the second layer we include the rotation θ so that σ ρ θ = e +iθσz/2 σ ρ e −iθσz/2 . ∆K = K θ −K is the relative shift of the Dirac cones between the layers. The coordinates correspond to those of layer 1 as measured from the K-point, but shifted with a vector +∆K/2 for layer 1 and −∆K/2 for layer 2. With this choice, the relative momentum k on both layers corresponds to the same absolute momentum. Furthermore µ is the chemical potential describing the effect of doping, here taken to be equal in both layers.
The off-diagonal terms in the Hamiltonian describe the coupling between the two layers. The matrix element at valley ρ between a state in sublattice α in layer 1 and a state in sublattice β in layer 2 is where δ αβ (r) is the horizontal displacement vector between the site at r, sublattice α in layer 1 and the nearest neighbor at sublattice β in layer 2. δ 1 denotes one of the nearest neighbor vectors connecting the graphene A and B sublattices. The sum is over the graphene lattice sites in the superlattice unit cell, and N denotes the number of these sites. For the interlayer hopping energy t ⊥ (δ) we use the same Slater-Koster parameterization as in Ref. [15]. Furthermore we approximate the interlayer coupling by only including the matrix el- , since they are an order of magnitude larger than the rest. For θ ≈ 1 • , the electronic dispersion becomes almost flat [19] and the group velocity d p /dp tends towards zero. In Fig. 2 we plot the resulting normal-state dispersion (ac) and the (local and total) density of states (d-i) close to this "magic" angle. The exact value of this magic angle depends on the details of the hopping model. In our case it is around 0.96 • , i.e., somewhat lower than what was found in [19]. However, the qualitative behavior of the local density of states is rather similar to the previous models. In particular, there are two closely spaced DOS peaks signifying the flattening of the bands. The local density of states is plotted along the line shown in Fig. 1, including three high-symmetry points with AB, AA, and BA stacking. These correspond to r = −1/3, 0, and 1/3, respectively.
Superconducting state. We assume that there is a local attractive interaction with strength λ (see Appendix for further details). We do not consider the specific nature of the pairing interaction and for the purposes of this paper it can be mediated by phonons or other bosonic modes. This model disregards the retardation effects due to such a mechanism, but is a valid approximation to the more general Eliashberg approach for weak coupling [27,28]. However, we also disregard direct Coulomb interaction. As shown in [28], within mean field theory and for weak coupling, it leads to an effective attractive interaction strength λ eff , smaller than the bare λ. However, as long as λ eff > 0, there is a possibility for a superconducting state.
Within a mean-field theory in the Cooper channel we find a self-consistency equation for a local superconducting order parameter (see Appendix for further details). Assuming furthermore that this order parameter shares the periodicity of the moiré superlattice, we find the selfconsistency equation for the superconducting order parameter, where the band sum b is calculated over the positive energy bands, α ∈ {A, B} is the sublattice index and i ∈ {1, 2} is the layer index. u ρbk and v ρbk are the eigenvectors of the Bogoliubov-de Gennes equation, We solve this self-consistent order parameter with a few values of the interaction constant λ and for a few different twist angles close to the magic angle. We include in the sum the energy levels closest to zero energy. We have checked that the results do not change much if the higher-energy levels are also included, up to energy cutoffs of the order of 100 meV. For comparison between different angles, we choose the chemical potential corresponding to the charge neutrality, marked in Figs. 2(g-i) with a dashed line. The resulting total density of states is plotted in Fig. 2(j-k), to allow for a comparison to the normal state. The corresponding local density of states (not shown) has the same localized structure as in the normal state, but the energy-dependence is modified similarly as the total DOS. The maximum of the position dependent ∆, which according to numerics is equal in both layers and sublattices, is plotted in Fig. 3 the magic angle, the Fermi speed v F (θ) increases so that the chosen λ is below the critical value λ c . This is why ∆ vanishes for angles θ 1.1 • . We can analyze the resulting magnitude of ∆ based on a flat-band result (assuming a constant ∆ and E ρbk ≈ ∆ for an extreme flat band) according to which (see Appendix for further details) ∆ = λΩ FB /π 2 , where Ω FB ≈ Ω moiré = 8π 2 /( √ 3||t 1 || 2 ). This yields ∆ = 1.3×10 −3 λ/a 2 for m = 34 corresponding to the magic angle. For comparison a linear fit to the linear region in Fig. 3(b) gives max(∆) = −0.07 meV+0.96×10 −3 λ/a 2 . The magnitude hence agrees very well with this simple model.
In Fig. 4 we show the temperature dependence of ∆ for m = 34, from which we may infer the approximate linear relation k B T c ≈ 0.25 max(∆(T = 0)) for the critical temperature. The prefactor is somewhat lower than for an extreme flat band with a constant ∆, for which k B T c = ∆/2 (see Appendix for further de- tails). The difference is most likely explained by the nonvanishing bandwidth and the position-dependent ∆ of our model. The maximum critical temperatures for the models calculated in Fig. 3(a) range from 3 K to about 20 K. The lower end of these values, calculated with λ = 1 eVa 2 , is thus quite close to that found in [14]. We stress that this is the mean-field critical temperature; the observed resistance transition is most likely rather a Berezinskii-Kosterlitz-Thouless transition and hence somewhat lower than the mean-field one. Note that despite the flatness of the bands, the supercurrent can be non-vanishing in the case when the eigenstate Wannier functions overlap [29] as is the case for TBG. Dependence on electrostatic doping. Besides θdependence, we can check how doping away from the center of the two DOS peaks affects the superconducting state. Figure 5(a) shows the charge density vs. chemical potential in the normal state. In Fig. 5(b) we plot the order parameter max(∆(δµ)) for different values of the doping δµ as measured from the charge neutrality point. Close to the magic angle, the energy scale of superconductivity exceeds that of the normal-state dispersion, and hence the only effect of the doping is to move away from the flat-band regime, suppressing superconductivity [30]. For the lowest value of λ in the plot, max(∆) is smaller than the bandwidth, and hence doping to the DOS peaks from their center can enhance superconductivity, resulting in a non-monotonous doping dependence. Note that as the effect of direct Coulomb interaction neglected here depends on the amount of screening, it also most likely depends on the charge density. As a result, this may lead to a more complicated ∆(δµ). However, as this depends a lot on the electromagnetic environment of the system, we will study this effect in a further work.
Concluding, we find that a BCS-type mean field model with relatively weak attractive interaction constant possibly even due to electron-phonon coupling (see Appendix for further details) can explain the occurrence of superconductivity in twisted bilayer graphene. We also make a number of predictions concerning the mean-field superconducting state, in particular the density of states and doping dependence. Our results form hence a checkpoint for further studies that use a simplified picture of the TBG flat band states, or consider mechanisms beyond the one in this paper. Our results could also have relevance in explaining the observations of superconductivity in twisted interfaces of graphite [31][32][33]. Clearly our mean-field theory fails to explain the insulator state [34] found experimentally in TBG at n ≈ ±2e/A moiré as well as the lack of superconductivity for hole doping [14,35]. There have been many suggestions of an unconventional (or more unconventional than here) superconducting state both for regular graphene [36,37] and for TBG [24,[38][39][40][41][42][43][44][45][46], typically directly related with the Coulomb interaction, and in some cases related with non-local interactions. We point out that our simple BCS model disregards the strain effects in moiré bands, as well as the possible dependence of the interaction constant on the twist angle. Whereas such mechanisms may play a role in TBG, we believe that the simplest BCS-type mean field superconductivity should also be considered as a viable effective model of the observations. Nevertheless, even in this case superconductivity would be highly exceptional, for example because it can be so strongly controlled by electrostatic doping.
The Hamiltonian for the attractive interaction is (A1) where α ∈ {A, B} is the sublattice index and i ∈ {1, 2} is the layer index. Doing the mean field approximation in the Cooper channel for Eq. (A1), we find a selfconsistency equation for a local superconducting order parameter where u ρbk,αi is the (α, i)-component of the spinor u ρbk . The spinors u ρbk and v ρbk are determined by solving the Bogoliubov-de Gennes equation (A3) Substituting the Bloch wave expansion for the eigenstates into Eq. (A3) and assuming ∆(r) to be periodic in the superlattice, we find the Fourier space Bogoliubov-de Gennes equation [Eq. (5) in the main text] and the Fourier space version of the self-consistency equation [Eq. (4) in the main text].

Charge density
The non-coupled system of twisted bilayer graphene is charge neutral at the chemical potential µ = 0. The charge density due to the electrons at that potential is where e is the electron charge and the factor of 2 comes from the spin. We formulate the calculation so that the k-sum goes over the superlattice Brillouin zone L * BK /SL * , B is the set of bands and 0,bk is the noninteracting dispersion. f 0 is the Fermi-Dirac distribution function at zero temperature. In the second step we introduce a cutoff by dividing the sum over the bands into two terms; to a sum over a set of low-energy bands Ω and to a sum over high-energy bands B\Ω.
In the presence of interlayer coupling, (normal state) dispersion changes to bk . The number of bands stays constant and if the interactions, temperature and chemical potential are small compared to the energy of the lowest energy band (in absolute value) of B\Ω in the non-interacting case, the index set B can be chosen so that the bands in B\Ω that were full (empty) in the noninteracting case, are still full (empty) in the interacting case. The interacting charge density is (A6) where f is the Fermi-Dirac distribution at temperature T and n high has the the same value as in (A5). The above has been formulated in the non-linearized theory. To calculate the excess charge relative to the charge neutrality point in the linearized theory, we split the bands between the two valleys and find where n is the excess charge density and Ω is now the set of bands in one valley.
The charge neutrality point µ * is determined from the equation n(µ * ) = 0. It is shown for different twist angles in Figs. 2(g-i) of the main text, and is always located in the middle between the two DOS peaks.

Simplified model of the superconducting state
The notion of weak or absent electron-phonon mediated superconductivity in pristine graphene is widely known. Here we reconcile this notion with our results claiming that a quite simple BCS-style model could describe the observations of superconductivity in twisted bilayer graphene. These results are not new, but we follow especially the treatments in [5,30] and adopt to the notation of the main paper, along with some estimates. We start from the generic self-consistency equation for the mean-field order parameter ∆. If ∆ is position independent, the Bogoliubov-de Gennes equation can be solved to yield where the prefactor 4 comes from summation over the valley and band indices, where in the band sum we include only the doubly degenerate lowest positive energy band. The cutoff k c is specified more below. We moreover assume that E k = 2 k + ∆ 2 . Here and below, without loss of generality we assume ∆ = |∆| ≥ 0. Our idea is to solve the self-consistency equation in three cases: (i) at the Dirac point for a Dirac dispersion 2 k = 2 v 2 F k 2 , (ii) for a Dirac dispersion at non-zero doping µ, i.e., 2 k = ( v F k − µ) 2 , and (iii) for a flat band with and without doping, k ≈ µ for k ∈ Ω FB . In each case we have the normal-state solution ∆ = 0, which we exclude by dividing both sides in Eq. (A8) by ∆.
Note that Eq. (A8) does not represent the full selfconsistency equation solved in the main text. Rather, we use it here simply to provide estimates of the behavior of ∆ in various limits.

Linear dispersion, no doping
Far away from the magic angle, the twisted bilayer behaves as if the two graphene layers would be almost uncoupled. This means that the low-energy dispersion exhibits two separate copies of the graphene Dirac dispersion. Inserting an ultraviolet energy cutoff c = v F k c and performing the integral for the T = 0 gap function, the self-consistency equation goes to the form where λ c = π 2 v 2 F /(2 c ). Since ∆ ≥ 0, this solution makes sense only if λ > λ c , and otherwise the only possible solution is the normal state ∆ = 0.
In pristine graphene, the critical interaction strength can be written also in terms of the nearest-neighbour hopping term γ 0 ≈ 3 eV [47].
Namely, within a nearest-neighbour tight-binding model the Fermi speed of graphene is v F = √ 3γ 0 a/(2 ), where a is the graphene lattice constant. We hence get If the attractive interaction results from electron-phonon coupling, a typical cutoff energy could be of the order of the Debye energy 200 meV [48]. In this case λ c ≈ 50 eVa 2 , 5 to 50 times larger than the values of λ used in our work. 100 to 200 meV is also the range of the maximum cutoff energy that we have used in our numerical results when including the contribution from higher bands. Even if the cutoff c would be of the order of γ 0 , the resulting λ c would be one order of magnitude larger than the smallest λ used in our results.

Linear dispersion, with doping
Let us try to reconcile the observations of superconductivity in Li or Ca doped graphene with the above idea. These cases are more accurately described by [49] within the Eliashberg theory. Here we just show in which sense doping fits into the above picture. Assuming k = ± v F k − µ and ∆ < c , and cutting the integral at k = c the self-consistency equation at T = 0 becomes [5] π 2 v 2 F 2λ = ∆ 2 + 2 c − ∆ 2 + µ 2 +|µ| ln |µ| + ∆ 2 + µ 2 ∆ .
(A12) Let us assume that ∆ |µ|, c so that we can expand the right hand side in ∆. In this case we find an analytic solution for ∆, Let us assume a cutoff energy c = 200 meV and a coupling strength λ = λ c /22 (corresponding to about 5 eVa 2 with the above estimates). In this case, with µ = 0.7 eV we would get ∆ = 1.3 meV. This corresponds to a critical temperature of 9 K, in the same range as the one that was measured in Li or Ca doped graphene [1][2][3][4].

Flat band estimate
Let us now make similar estimates for the flat-band case of the moiré superlattice. In this case, we assume that ∆ is larger than the bandwidth of the lowest-energy band. Within that band, we can hence approximate E k ≈ ∆ in Eq. (A8) and at T = 0 the integral is over a constant function. As a result, we get where Ω FB = 8π 2 /[ √ 3(3m 2 + 3m + 1)a 2 ] is the area of the first Brillouin zone of the moiré superlattice. Within the model adapted in the main text, the magic angle is around m ≈ 34, in which case we would get ∆ FB = 1.3 × 10 −3 λ/a 2 . In Fig. 3(b) of the main text, the solid line has a slope of 1 × 10 −3 λ/a 2 , i.e., very close to this simple estimate.
The temperature dependent ∆ in the flat-band case is obtained by solving At the critical temperature, ∆ → 0, and we can hence expand the right hand side to the linear order in ∆/(2k B T c ). This directly yields T c = ∆ FB /(2k B T ).
In the case of a non-zero potential µ, we can use E k ≈ µ 2 + ∆ 2 in the self-consistency equation. It then becomes (for ∆ > 0) In this case superconductivity is hence suppressed when the absolute value of the chemical potential is larger than ∆ FB .