Superconductivity from collective excitations in magic angle twisted bilayer graphene

A purely electronic mechanism is proposed for the unconventional superconductivity recently observed in twisted bilayer graphene (tBG) close to the magic angle. Using the Migdal-Eliashberg framework on a one parameter effective lattice model for tBG we show that a superconducting state can be achieved by means of collective electronic modes in tBG. We posit robust features of the theory, including an asymmetrical superconducting dome and the magnitude of the critical temperature that are in agreement with experiments.


Introduction:
The remarkable experimental observations of superconducting and insulating phases in twisted bilayer graphene (tBG) [1][2][3] close to the magic angle ushers a new paradigm in attempts to study strongly correlated phases of matter by bandstructure engineering.On stacking two graphene monolayers at a small relative twist angle θ, emerges a large wavelength moiré superlattice potential.Further, band calculations reveal [4] that there is substantial renormalization of the Fermi velocity giving rise to flat bands around half filling once θ is sufficiently small.At some very specific 'magic angles' (θ M ∼ 1.1 • ), the Fermi velocity is predicted to vanish facilitating strong electronic correlations.It is therefore plausible that the superconducting instability observed in twisted bilayer graphene is a consequence of these correlations.
Given the rich underlying physics at play, a lot of theoretical effort has already been devoted towards understanding the Mott-like physics of strongly correlated electrons in tBG [5][6][7][8][9][10][11][12][13][14][15][16][17], and attributing the observed superconducting state to weak electron-phonon coupling [18][19][20][21][22].It is however not obvious that the latter is a conventional Bardeen-Cooper-Schriefer (BCS) state even though the phonon mediated superconducting pairing is amplified in tBG due to the enhanced density of states at the Fermi surface (see Fig. 1).Neither is there a universal consensus or understanding that Mott-like physics is in action.In contrast, Coulomb interactions and associated plasmonic effects are known to be very strong in graphene [23], and theoretically [24], may lead to a superconducting state.In this Letter we explore the possibility of unconventional superconductivity in tBG mediated by the purely collective electronic modes.
We first emphasize that the electronic mechanism is quite different from the standard BCS interaction, as no phonon modes are necessary.The collective oscillations can generate an effective dynamic attractive interaction, thus Cooper-pairing two electrons.Starting from a oneparameter effective lattice model for tBG [25] we first calculate the dynamical polarization function Π(q, iω) and thereby the dynamically screened Coulomb interaction V (q, iω).Following Grabowski and Sham [24], we average the interaction kernel to a dimensionless momentumindependent interaction parameter, λ nm , and obtain the gap equation as follows: where ω n = (2n+1)πT C is the Matsubara frequency, and Z m is a renormalizing function calculated below.Solving Eq. ( 1) we find the superconducting critical temperature T C for various twist angles θ and carrier con-arXiv:1909.02574v1[cond-mat.mes-hall]5 Sep 2019 centrations linked to the Fermi energy E F via density of states.The prominent predictions of our theory are: (i) the plasmon-assisted Cooper pairing is much stronger than the conventional phonon-related one in tBG, (ii) the superconducting state does not occur at low electron concentrations but is prominent at electron densities around n = 10 12 cm −2 , (iii) the obtained critical temperature is of O(10K), similar to that observed in experiments [1][2][3] and larger than those calculated by the phonon contribution [18], (iv) our theory is not limited by particular plasmon modes as our model generically captures the effect of plasmons (whenever well-defined), but also the other density-fluctuation excitations.Our work specifically points out to the mechanism for superconductivity that is likely to be at play in tBG.
Phonons or plasmons?The striking scenario which develops at small twist angles is that the relevant electronic energy scales are shrunk to the order of a few meV (similar to the acoustic phonon energy scale), while the density of states is amplified.This raises the intriguing possibility of observing collective electronic modes in the same energy window.For example, it was recently [25] pointed out that plasmon modes in tBG can remain intrinsically undamped protected by the bandgap between the flat bands and higher bands.It is also worth noting the possibility to have hybrid acoustic phonon and plasmon modes still uncommon in condensed matter systems.
Let us compare the phonon and plasmon mediated superconductivity mechanisms in tBG with θ close to θ M .Here, we assume a simplest two-dimensional massless Dirac fermion model with the renormalized Fermi velocity v F .The Fermi-surface is then just a circle of radius k F , and the averaged electron-phonon coupling , where ζ, ρ, and c ph are deformation potential, mass density and the sound velocity, respectively.To the first approximation we can assume the attractive pairing potential to be finite only below the Debye frequency ω D (see Fig. 1a).The McMillan formula then suggests T C ∼ ω D exp(−1/λ) ∼ O(1K), as also shown recently [18].In contrast, the plasmon-mediated mechanism suggests that the dynamic Coulomb interaction V (q, iω) is responsible for the superconducting pairing.For such pairing the dynamical nature of the coupling is crucial because the Coulomb interaction is screened by the dielectric function (q, iω) = 1 + e 2 E F q/(2κω 2 ) [26,27] providing stronger interactions in the high-frequency limit.Here, the renormalized dielectric constant κ ∼ 12 accounts for effects of interband polarization in tBG [25].The corresponding dimensionless coupling λ nm reduces to the bare Coulomb interaction characterized by r s = e 2 /(κ v F ) in the high frequency limit, while remaining zero below a certain energy determined by ω b [24,28].This behavior being quite opposite to the phonon-assisted pairing is schematically shown in Fig. 1b.
Equation ( 1) is still too complicated, as it in fact represents an infinite number of coupled equations.We design a simple analytically tractable model by reducing the number of equations to three considering only the terms with m = 0, ±M , where M 1. Neglecting the selfenergy corrections for now (these will be included later on) and setting the diagonal elements λ mm = 0 we arrive at the following equation for T C : (2) We use the single-mode Eq. 2 with M ∼ ω b /E F to estimate T C in Fig. 1c.In the low-T C limit (T C E F ) and strong coupling (r s 1), we find an elegant expression that can be seen as a plasmonic analogue of the McMillan formula for T C .Since the left-hand side of Eq. ( 3) cannot exceed π/2 the solution for T C exists if and only if ω b /(r 2 s E F ) < 1, i.e. the electron concentration and coupling strength values must be large enough.This is indeed the case in tBG.Estimating ω b ∼ 15 meV, v F ∼ 1.5 × 10 4 m/s, we find r s ∼ 12, E F = v F πn/2 ∼ 1.5 meV (for electron density n ∼ 1.5 × 10 12 cm −2 ), and M ∼ 10, resulting in a reasonable value T C ∼ 2.6 K matching the observations [1].In conventional graphene r s ∼ 1, ω b > E F [27], and solution of Eq. ( 3) never exists making the plasmon model tBG-specific.This is consistent with the fact that superconductivity has been observed in tBG but not in monolayer graphene.

Dynamical Coulomb interaction:
To go beyond the single-mode model, we consider the following one-parameter effective nearest-neighbor tight binding Hamiltonian defined on a hexagonal lattice, which mimics the low energy Hamiltonian of twisted bilayer graphene [25] where c k,η and c † k,η are the annihilation and creation operators for electrons with momentum k in the Brillouin zone (BZ) on the sublattice η.The effective hopping matrix element t eff = W/3 corresponds to the bandwidth (W ) of the nearly flat bands in tBG close to the magic angle.The summation δ j is over the nearest neighbors δ j = (cos(2πj/3), sin(2πj/3))ã √ 3, where ã = (a/2 sin(θ/2)) is the periodicity of the moiré superlattice.The constants a = 2.46 Å is the graphene lattice constant, while θ is the twist angle.For our calculations we obtain the bandwidth W from the tBG continuum model [4].The advantage of the above tightbinding model is that it reproduces the symmetry of the actual tBG and has a natural ultraviolet smooth cutoff scale W .The divergent density of states at the van Hove singularity is also manifested in this model, which will be important for our analysis.The eigenvalues are E n k = n|h k | and the four-fold degenerate eigenstates are where n = ±1 is the band index.
We are particularly interested in the dynamic polarization function Π(q, iω), which is given by where n, m are the band indices, f n k is the Fermi-Dirac distribution, and F nm k,k+q = |ψ † n,k+q ψ m,k | 2 is the graphene chirality factor.Since h k has a complicated momentum dependence, we evaluate the above function numerically expanding up to second order in q.As recently shown [25], this is actually sufficient to describe collective excitations such as plasmons over a wide range of energy scales 0 < ω < 2W .The dynamical polarizability is the used to evaluate the dielectric constant within the random phase approximation (RPA) as (q, iω) = 1 − V q Π(q, iω), where V q = 2πe 2 /κq is the bare Coulomb interaction.Finally, the dynamically screened Coulomb interaction is given by V (q, iω) = V (q)/ (q, iω).The dynamical dielectric constant (q, iω) is related to the collective propagator χ co (q, iω) as (q, iω) −1 = 1 + V (q)χ co (q, iω).
The gap equation: The Migdal-Eliashberg superconductivity theory suggests the following gap equation [29] ∆(p, where the product of the Green's functions B(k, iω m ) given by Here, Σ(k, iω m ) represents the normal-state self energy, which is given by In the above expressions the energy dispersion k is measured from the Fermi energy.Note that we will now restrict our attention to the conduction band intersecting the Fermi energy for electron doped system.We point out that one may include the phonon contribution in this framework by adding the phononic propagator D(q, ω) into the gap equations.Since, we are interested in evaluating the pure electronic contribution, we do not attempt this calculation here and reserve it for future studies.It has been shown [24] that superconductivity in a Coulomb system with no attractive interactions is essentially determined by the frequency dependence of the screened Coulomb interaction.Therefore, we can work with an effective Coulomb interaction, which is averaged over the momenta up to k c = 2k F [24], thus retaining only the crucial frequency dependence.The momentum averaged interaction is given by The dimensionless coupling is given by λ(iω) = V (iω) N (E F ), where N (E F ) is the density of states.
The limiting cases of the coupling are given by lim ω=0 V (iω) N (E F ) → 0 and lim ω→∞ V (iω) N (E F ) = µ, where we define µ to be the high frequency limit of the coupling.To proceed further, the dimensionless coupling is then mapped onto the following explicit expression: where ω b is a boson frequency that sets the scale of transition from the low to the high frequency limit.The mapping of the kernel onto the above Lorentzian form allows for an analytical treatment that is indeed found to resemble closely to the actual numerical solution for λ nm .The parameters µ and ω b are then extracted by fitting the actual coupling λ(iω) to Eq. ( 10).The propagating boson here (with frequency ω b ) is a collective density-fluctuation excitation such as electronhole excitation or a plasmon.For the momentum averaged interaction, Eq. ( 6) becomes momentum independent and is given by Eq. ( 1) with

accounting for the self-energy corrections
The gap equation is now of the form of an eigenvalue equation ∆ = Ĉ ∆.At T C the largest eigenvalue of Ĉ is exactly one.Equation (1) must be solved numerically to obtain a reliable value for critical temperature.The rate of convergence of the numerical solution depends on the ratio T C /E F .If T C E F , the dimensions of the matrix involved can become prohibitively large to allow a numerical solution.However, in this regime the pseudopotential method [24], which assumes T C E F as a premise, gives us a good estimate of the superconducting critical temperature.For our calculations we use a combination of numerical solution and pseudopotential method depending on the ratio T C /E F .For θ close to θ M , the numerical solution converges within a reasonable computational time.As θ increases the pseudopotential method is used.We also briefly comment on the nature of the superconducting gap function ∆(iω).Since we have purely repulsive interactions (unlike the case of phonon mediated BCS coupling), the gap function changes sign as a function of the Matsubara frequency.This is necessary for the gap equations to have a non-trivial solution.
Superconductivity: We will now discuss salient predictions of our theory.In Fig. 2 we plot the critical temperature T C as a function electron density in tBG close to the magic angle.The presence of an asymmetrical superconducting dome around n = 10 12 cm −2 and the calculated T C = O(10K) are the main predictions of our theory.The obtained T C closely resembles the measured order of magnitude in experiments [1][2][3], and also the calculated T C from the phonon contribution [18].As the carrier density is increased the first maximum of the superconducting dome occurs when the Fermi energy intersects the M of the Brillouin zone point, where the van Hove singularity occurs.As the density is increased further, the dome is found to be asymmetric around this point, which can be understood from the fact that the density of states is not symmetric around the M point (see inset of Fig. 3).This asymmetry combined with effect of strong electronic interactions results in a second neighbouring maximum in the dome, which is more prominent at lower twist angles.When θ is increased further the width of the dome shrinks and the second maximum becomes less prominent since the electronic interactions become comparatively weaker.We observe superconductivity even close to θ ∼ 2 • , but the corresponding density window is quite narrow.In the monolayer limit the superconducting dome is practically non-existent within our model.The appearance of these domes is a special feature which arises from chosen realistic lattice model as opposed to the simple Dirac approximation for tBG.By examining the dependence on dielectric constant, we notice that the magnitude of T C is set primarily by r s (or µ), while the shape of the domes is set by the density of states of the non-interacting bands.Also note that we specifically focus on superconductivity.Other competing states (such as density waves, magnetism etc) may obviously affect the shape of the dome, when taken into consideration.In Fig. 3 we plot the obtained critical temperature T C as a function of twist angle for different carrier densities.The dependence is observed to be non-monotonic, and for large angles T C is eventually suppressed, as expected.Fig. 3 also compares the numerical solution to the T C obtained from single-mode model in Eq. (2).
Concluding remarks: We first point out that the vertex corrections are neglected in the Migdal-Eliashberg formalism, which may have a quantitative impact on the calculated T C .For twist angles close to the magic angle, we find that the propagating boson frequency (ω b ) can be several times larger than the Fermi energy E F .However in this regime, the vertex corrections can be ignored, because they turn out to be insignificant for processes much larger than E F as pointed out earlier in the literature [30].Therefore the calculated O(T C ) is expected to be robust especially close to the interesting regime of the magic angle.Vertex corrections may become more important for large twist angles causing suppression of T C , but evaluating them is beyond the scope of the current manuscript.
To conclude, this work specifically points out to an important mechanism which is likely to be at play in superconducting tBG close to the magic angle.Collective excitations of strongly coupled electrons can mediate pairing, which may be dubbed as 'plasmonic superconductivity', although the requirement of undamped plasmons is not strict in our formalism.We predict features of the theory, namely the superconducting dome and the magnitude of T C , which are in good agreement with recent experimental data [1,3].The nature of the gap function in momentum space may be inferred by solving the Eliashberg equations without any momentum averaging.We reserve this technically harder problem for future studies.

Figure 1 .
Figure 1.The single-mode model for the dimensionless pairing potential N (EF ) V eff (ω) (N (EF ) being the density of states and V eff (ω) the effective interaction) in tBG for (a) moiré phonons (with Debye frequency ωD) within BCS theory [18], and (b) collective electronic excitation (with frequency ω b ) such as a plasmon within Eliashberg theory.In this system with purely repulsive interactions, the frequency dependence of the pairing potential can induce superconductivity.(c) TC evaluated within the single-mode approximation (Eq.2) for a Dirac model with renormalized vF ∼ 1.5 × 10 4 ms −1 , ω b ∼ 15 meV.Below a threshold density (indicated by the dotted lines), no solution for TC exists.

Figure 2 .
Figure 2. Critical temperature TC as a function of electron density in tBG (close to the magic angle).The presence of an asymmetrical superconducting dome around n = 10 12 cm −2 and TC = O(10K) are the main predictions of our theory.

Figure 3 .
Figure 3. Critical temperature TC as a function of twist angle in tBG for three different carrier densities.The solid lines show the TC obtained from Eq. 2 using the parameters ω b , EF , and µ(rs) from the numerical solution.The inset shows the normalized density of states (for any angle) at the Fermi energy as the Fermi surface moves up in energy and intersects various symmetry points in the Brillouin zone K → M → Γ.