Pairing in magic-angle twisted bilayer graphene: role of phonon and plasmon umklapp


Identifying the microscopic mechanism for superconductivity in magic-angle twisted bilayer graphene (MATBG) is an outstanding open problem. While MATBG exhibits a rich phase-diagram, driven partly by the strong interactions relative to the electronic bandwidth, its single-particle properties are unique and likely play an important role in some of the phenomenological complexity. Some of the salient features include an electronic bandwidth smaller than the characteristic phonon bandwidth and a non-trivial structure of the underlying Bloch wavefunctions. We perform a systematic theoretical study of the cooperative effects due to phonons and plasmons on pairing in order to disentangle the distinct role played by these modes on superconductivity. We consider a variant of MATBG with an enlarged number of fermion flavors, \(N \gg 1\)​, where the study of pairing instabilities reduces to the conventional (weak-coupling) Eliashberg framework. In particular, we show that certain umklapp processes involving mini-optical phonon modes, which arise physically as a result of the folding of the original acoustic branch of graphene due to the moiré superlattice structure, contribute significantly towards enhancing pairing. We also investigate the role played by the dynamics of the screened Coulomb interaction on pairing, which leads to an enhancement in a narrow window of fillings, and study the effect of external screening due to a metallic gate on superconductivity. We propose a smoking-gun experiment to detect resonant features associated with the phonon-umklapp processes in the differential conductance and also discuss experimental implications of a pairing mechanism relying on plasmons.

1 Introduction

Since the discovery of superconductivity (SC) [1] near an interaction-induced insulator [2] in magic-angle twisted bilayer graphene (MATBG), the field has evolved dramatically [3]–[5]. One of the defining characteristic features of MATBG is the emergence of isolated, (nearly) flat bands separated by a large energy gap to higher dispersive bands [6]–[8]. A variety of symmetry broken states and phase transitions have been observed experimentally [4], [9]–[14] as a function of electron filling (\(\nu\in [-4,4]\)​) measured in units of electrons per moiré unit cell [15]. In spite of some similarities between the phenomenology in MATBG and copper-oxide based materials [16], MATBG bears its own set of unique properties, which likely play a central role in the underlying microscopic origins of superconductivity and other experimentally observed features. For example, MATBG has spatially extended Wannier wavefunctions [17], a non-trivial band topology [18]–[21], a complex phonon band structure [22], [23] that may extend beyond the electronic bandwidth, and plasmons that decouple from the particle-hole continuum associated with the narrow bands [24].

These unique features call for a detailed analysis of their influence on the low-energy electronic properties in MATBG. For example, while the enhancement of Coulomb interactions relative to the narrow bandwidth undoubtedly plays an important role in stabilizing the insulators at various commensurate fillings [2]–[4], the extent to which they are crucial for superconductivity remains unclear. A number of recent experiments have attempted to study the role of Coulomb interactions, either by varying the distance between the MATBG layer and a nearby metallic screening gate across different devices [25], [26],1 or, by independently varying the density of carriers in a layer of Bernal-stacked bilayer graphene [28], situated a fixed distance away from MATBG in the same device. Not surprisingly, the enhanced screening suppresses the various insulating phases, that develop originally at a sequence of commensurate fillings. On the other hand, SC is affected only weakly, if at all, by the gate. A number of recent theoretical studies of superconductivity in MATBG have relied on a variety of BCS mean-field and other weak-coupling based approaches [29]–[39]; see also two recent studies [40], [41] focusing on the interplay of insulating and superconducting states in MATBG, where the underlying band topology plays a crucial role.

Inspired by the growing number of such interesting experiments, we focus on a possible microscopic mechanism for the superconducting instabilities in MATBG due to an interplay of attraction generated by phonons and purely electronic collective modes, such as plasmons. Importantly, we take into account the retarded interaction mediated by the acoustic phonons and incorporate the dynamics of the screened Coulomb interaction within Eliashberg theory [42]. In the experimental regime of interest, MATBG is likely defined by an "intermediate-coupling" problem with no natural small parameter. In order to make theoretical progress, we introduce a controlled large\(-N\)​ expansion, such that our results are asymptotically exact in the weak coupling limit (defined by \(N\rightarrow\infty\)​). We argue that in spite of this approximation, the results shed new light on the interplay of different sources of attraction on SC in MATBG. Our theoretical approach is not suited to discuss the insulating states observed in MATBG at the integer fillings [4], [12], [13]. However, it is entirely possible that the interactions responsible for SC are distinct from those responsible for the insulating phases [43]. This work addresses only the possible role played by the interplay of phonons and screened Coulomb interactions on SC in MATBG in a systematic fashion. It is entirely conceivable that the strong-interaction induced insulating states, that develop at commensurate fillings, simply punctuate the SC phase-diagram that we discuss and present in this study.

One of the unique features of MATBG is the large size of the unit-cell, compared to the inter atomic distance. Consequently, higher order umklapp processes involve momenta which are still small compared to the original Brillouin zone (BZ) and thus, play an important role in phonon-mediated SC.2 We show below that the SC transition temperature, \(T_c\)​, due to acoustic phonons of graphene is enhanced significantly upon including the effects of the umklapp processes for a wide range of dopings. On the other hand, pairing due to plasmons is dominant only in a narrow range of fillings \(\nu \approx 2-3\)​ and is much less sensitive to the inclusion of umklapp processes.

Given that phonons can play a crucial role in pairing, we propose a tunneling-based experiment to investigate the spectral signatures of the umklapp processes at an energy that is related to an integer times the size of the BZ and the speed of sound in graphene. More generally, we propose here that probing the replica bands in MATBG spectroscopically is a promising avenue to study the coupling of electrons to the various bosonic modes (including phonons) in the system. Finally, we study the effect of a screening layer, coupled to MATBG by Coulomb interactions, on the density dependence of \(T_c\)​.

The remainder of this paper is organized as follows: in Sec. 2, we review the continuum model for MATBG and introduce the Coulomb and electron-phonon interactions, projected to the nearly flat bands. We also highlight some of the technical aspects associated with umklapp scattering and the large\(-N\)​ formulation in this section. In Sec. 3, we focus on setting up the (linearized) Eliashberg equations for an inter-valley, spin-singlet \(s-\)​wave pairing gap and highlight some salient features associated with the pairing kernel. Sec. 4 is devoted to our results, which includes a comprehensive study of the filling dependence of \(T_c\)​ due to the different sources of attraction, the role of external screening and a discussion of our proposed experimental setup for investigating the fingerprint of umklapp phonon processes. Sec. 5 concludes with a discussion and future outlook for interesting directions. We describe some of the technical aspects of our work in appendices 6 - 8.

Figure 1: (a) A small twist-angle between two graphene layers leads to formation of mini-Brillouin zones located in the vicinity of the \(\boldsymbol K\) and \(\boldsymbol K'\) points. We focus on intervalley (\(\xi \neq \xi'\)) pairing (dashed line) to form Cooper pairs at zero center of mass momentum. (b) The moiré interlayer coupling gives rise to formation of two (spin degenerate) nearly-flat bands for each valley \(\xi\). \(\bar{K},\bar{\Gamma},\bar{M},\bar{K}'\) correspond to high symmetry points of the MBZ. Positive (\(\nu>0\)) and negative (\(\nu<0\)) fillings correspond to the electron and hole band, respectively. (c) The energy landscape of an electron band for \(\xi=1\). The red arrows depict the umklapp processes (\(0\)G, \(1\)G, \(2\)G) between states originating from different MBZs (one shown as white hexagon) in the extended zone scheme. Here \(K=|\bar{K}-\bar{\Gamma}|\).

2 Model

The non-interacting part of our model, \(\mathcal{H}_0\)​, leads to the action,

\[\begin{equation}\mathcal{S}_{0}=\sum_{i=1}^N\sum_{\boldsymbol{k},{\omega},\{{\gamma}\}} \left(-i{\omega}+E_{\boldsymbol{k}\{{\gamma}\}}\right) c^\dagger_{{\omega}\boldsymbol{k}\{{\gamma}\};i} c_{{\omega}\boldsymbol{k}\{{\gamma}\};i},\end{equation}\]

where \(c_{{\omega}\boldsymbol{k}\{{\gamma}\};i}\)​ is the electron annihilation operator in the eight-dimensional space of \(\{{\gamma}\}\equiv\{\xi,{\sigma},n\}\)​, with \(\xi\)​ being the valley index of the original graphene layers, \(\sigma\)​ labeling the electron spin and \(n\)​ tracking the electronic band index of MATBG. We have introduced a fictitious index \(i=1,..,N\)​ which labels \(N\)​ identical copies of the non-interacting MATBG Hamiltonian, with an eye for carrying out a controlled \(1/N\)​ expansion.3 We denote \(\boldsymbol{k}\)​ as the crystal momentum in the original graphene layer, such that the two mini-Brillouin zones (MBZ) are located near \(\boldsymbol{K}\)​ and \(\boldsymbol{K'}\)​ points of the original graphene layers, and the \(\boldsymbol{k}\)​ summation ranges over the MBZ, see schematic in Fig.1a.

The energy bands, \(E_{\boldsymbol{k} \{{\gamma}\} }\)​, are computed from the effective continuum Hamiltonian introduced in Ref. [45] for MATBG at a twist angle \(\theta=1.05^\circ\)​.4 The band structure of the two narrow bands, which are relevant for superconductivity are shown in Fig.1b. In the analysis that follows, we introduce a chemical potential (\(\mu\)​) in the usual way by shifting the single-particle energies, \(E_{\boldsymbol{k} \{{\gamma}\} } \to E_{\boldsymbol{k} \{{\gamma}\} } - \mu\)​, and control the electron filling. Additional details associated with the bandstructure appears in Appendix 6.

Let us now turn our attention to the interaction terms in our model. As a result of the large momentum separation between \(\boldsymbol{K}\)​ and \(\boldsymbol{K'}\)​ in the original BZ (see Fig.1a), the states near \(\xi=\pm\)​ valleys in the two MBZs are effectively decoupled. The interaction is then given by, \[\begin{gather} \label{eq:ham-95int} &\mathcal S_{\textrm{int}} = \frac{1}{2}\sum_{\substack{\boldsymbol{q},{\omega}\\ \xi,\xi',i}}{\cal V}^{\mathrm{int}}_{\xi,\xi'}(\boldsymbol{q},i{\omega})\nonumber\rho_{\xi; i}(\boldsymbol{q},i{\omega})~\rho_{\xi';i}(-\boldsymbol{q},-i{\omega}),\\ &\rho_{\xi; i}(\boldsymbol{q},i{\omega}) = \sum_{\boldsymbol{k}, \nu} \Lambda_{{\gamma}{\gamma}'}(\boldsymbol{k}+\boldsymbol{q},\boldsymbol{k})~c^\dagger_{\nu+{\omega}\boldsymbol{k}+\boldsymbol{q}\{{\gamma}\};i} c_{\nu\boldsymbol{k}\{{\gamma}'\};i}\,, \end{gather}\]

where the density operators, \(\rho_{\xi; i}(\boldsymbol{q},i{\omega})\)​, include the form-factors, \(\Lambda_{{\gamma}{\gamma}'}(\boldsymbol{p},\boldsymbol{k})=\delta_{\xi\xi'}\delta_{{\sigma}{\sigma}'}\left\langle \boldsymbol{p},\{{\gamma}\}\left|e^{i (\boldsymbol{p-k})\cdot \boldsymbol{r}} \right|\boldsymbol{k},\{{\gamma}'\}\right \rangle\)​, which necessarily involve states located near the same valley \(\xi\)​ as the Hamiltonian is block-diagonal in the valley space. Here \(|\boldsymbol{k},\{{\gamma}'\} \rangle\)​ denotes a Bloch wavefunction of the non-interacting Hamiltonian (see Appendix 6). Note that we have assumed explicitly that the interaction term is diagonal in the fictitious \(i=1,\dots,N\)​ index. The interaction vertex above is a sum of two contributions — the dynamically screened Coulomb interaction (\({\cal V}^{\textrm{C}}_{\xi,\xi'}\)​) and the phonon mediated interaction (\({\cal V}^{\textrm{ph}}_{\xi,\xi'}\)​): \[\begin{equation}{\cal V}^{\mathrm{int}}_{\xi,\xi'}(\boldsymbol{q},i{\omega}) = {\cal V}^{\textrm{C}}_{\xi,\xi'}(\boldsymbol{q},i{\omega})+{\cal V}^{\textrm{ph}}_{\xi,\xi'}(\boldsymbol{q},i{\omega})\,.\label{eq:Int} \end{equation}\]

To stress the intervalley nature of the pairing in what follows, we make the labels \(\xi\)​ and \(\xi'\)​ explicit in both the Coulomb and phonon interactions. The microscopic form of the interaction however is the same for both inter-/intravalley interactions with the dependence on valley indices \(\xi,\xi'\)​ contained only in the density operators \(\rho_{\xi}(\boldsymbol{q},i\omega)\)​ introduced in Eq. \(\ref{eq:ham-95int}\).

The first term in Eq. \(\ref{eq:Int}\) is the dynamically screened Coulomb interaction, which is given by \[\begin{equation}{\cal V}^{\textrm{C}}_{\xi,\xi'}(\boldsymbol{q},i{\omega}) = {1\over N}{2\pi e^2\over {\varepsilon}_{\mathrm{RPA}}(\boldsymbol{q} ,i{\omega}) q}\,,\label{eq:V-95C} \end{equation}\]

where the RPA dielectric function is given by \({\varepsilon}_{\mathrm{RPA}}(\boldsymbol{q},i{\omega}) = \kappa - 2\pi e^2 \Pi_{ee}(\boldsymbol{q} ,i{\omega})/q\)​. Here \(\Pi_{ee}(\boldsymbol{q}, i{\omega})\)​ is the electronic polarization of MATBG [24] and \(\kappa\)​ is the background dielectric constant, which can in principle be frequency and momentum dependent as well.

The second term, the phonon-mediated interaction, arises from the electron-phonon coupling Hamiltonian, \[\begin{equation}{\cal H}_{\textrm{el-ph}} = -i\sqrt{\frac{gc_s}{2N}}\sum_{\xi, \boldsymbol{q}} \sqrt{ q}~ \rho_{\xi; i}(\boldsymbol{q}) (a_{\boldsymbol{q}}+a^\dagger_{-\boldsymbol{q}}),\label{eq:ham-95phonon} \end{equation}\]

where \(a_{\boldsymbol{q}},~a^\dagger_{\boldsymbol{q}}\)​ represent the phonon annihilation and creation operators, respectively5. The phonon coupling constant, \(g = D^2/\rho_m c_s^2\)​, is related to the deformation potential, which we set as \(D~(= 25~ \textrm{eV})\)​, and we also use \(c_s~ (= 12 000~\textrm{m/s})\)​ for the speed of sound in graphene and \(\rho_m ~(= 7.6\times 10^{-8} ~\mathrm{g}/\mathrm{cm}^2)\)​ for the atomic mass density [46], [47]. After integrating out the phonons, we obtain the density-density interaction term of the form Eq. \(\ref{eq:ham-95int}\) where, \[\begin{equation}{\cal V}_{\xi,\xi'}^{\textrm{ph}}(\boldsymbol{q},i{\omega}) = -{g\over N}{\omega^2_{\textrm{ph}}(q) \over {\omega}^2 + \omega^2_{\textrm{ph}}(q)},\label{eq:V-95ph} \end{equation}\]

and \(\omega_{\textrm{ph}}(q)=c_s q\)​ is the acoustic phonon dispersion for graphene. Note that we have chosen the individual \(N-\)​dependent normalizations such that both \({\cal V}^{\textrm{c}}_{\xi,\xi'}\)​ and \({\cal V}_{\xi,\xi'}^{\textrm{ph}}\)​ have the same \(1/N\)​ prefactor and can therefore, be compared in a meaningful fashion in the large\(-N\)​ limit.

Let us now clarify an important aspect of the underlying scattering processes, that play a crucial role in our subsequent analysis. The summation over the crystal momenta, \(\boldsymbol{k}\)​, in Eq. \(\ref{eq:ham-95int}\) is limited to the first MBZ. However, the momentum transfer, \(\boldsymbol{q} \equiv \boldsymbol{G} + \tilde {\boldsymbol{q}}\)​, can scatter to Bloch states belonging to different MBZs, where \(\tilde {\boldsymbol{q}}\)​ is now restricted to lie in the first MBZ. The lattice vector, \(\boldsymbol{G}=m_1 \boldsymbol{G}_{1}^{M}+ m_2 \boldsymbol{G}_{2}^{M}\)​ (\(m_1,m_2 \in \mathbb{Z}\)​), accounts for scattering by the multiplicity of reciprocal moiré lattice vectors, \(\boldsymbol{G}_{1}^M\)​, \(\boldsymbol{G}_{2}^M\)​. In order to evaluate the importance of the umklapp processes described above, from now on we restrict our analysis to the \(m^{\textrm{th}}\)​ MBZ, such that \(m_1,m_2 \in [-m,m]\)​. We refer to such an analysis as including \(m\)​G processes, see Fig.1c. Varying \(m\)​ allows us to assess the relative importance of the successive umklapp processes. As we shall highlight later on, while including a few umklapp processes beyond the simplest 0G process leads to qualitatively new features, there is also rapid convergence with increasing \(m\)​, such that we do not need to include processes with arbitrarily large values of \(m\)​.

It is worth noting that these umklapp processes can be equivalently viewed as processes where mini-optical phonons are exchanged. These optical-like modes are a natural consequence of folding of the original acoustic phonon branch of graphene due to the moiré potential [22], and have no relation to the optical phonons of the original (decoupled) graphene layers.6

3 Eliashberg Equations

We are interested in the linearized gap equation in the pairing channel, ignoring the self-energy corrections to the electron dispersion and quasiparticle weight (which are subleading in \(1/N\)​). The gap function in the spin-singlet, \(s\)​-wave pairing channel with zero center of mass momentum is defined as, \[\begin{gather} \Delta(i\omega,{k}) \equiv {\langle}{\varepsilon}_{{\sigma}{\sigma}'} {\varepsilon}_{\xi\xi'} c_{\omega,\boldsymbol{k}\{{\gamma}\}} c_{-\omega,-\boldsymbol{k}\{{\gamma}'\}}{\rangle},\end{gather}\]

where \({\varepsilon}_{\alpha\beta}\)​ is the fully antisymmetric tensor in the indices \(\alpha,~\beta\)​ and we have further assumed that the gap function has no explicit dependence on the angle of \(\boldsymbol{k}\)​. Note that the above requires an intervalley scattering, Fig.1a, which imposes the condition \(\xi \neq \xi'\)​ in the interaction Eq. \(\ref{eq:ham-95int}\). The gap equation thereby reduces to an eigenvalue problem \[\begin{gather} \Delta(i\omega,{k}) = -T \sum_{\nu} \sum_{{p}} K(i{\omega},{k};i\nu,{p})\Delta(i\nu,{p}) \label{eq:sc-95gap-95equation},\end{gather}\]

where the kernel \(K(...)\)​ is given by \[\begin{gather} \label{eq:Kernel} K(i{\omega},&{k};i\nu,{p}) \equiv \\&{1\over (2\pi)^2} \int d{\Omega}_{\boldsymbol{p}} {\cal V}^{\textrm{int}}_{-\xi,\xi}(\boldsymbol{k}-\boldsymbol{p},i\omega-i\nu) \frac{\Lambda( \boldsymbol{p} ,\boldsymbol{k})\Lambda(-\boldsymbol{p},\boldsymbol{k}) }{\nu^2+{E}^2_{\xi,\boldsymbol{p}}}\nonumber\end{gather}\]

and where \(\int d {\Omega}_{\boldsymbol{p}}\)​ denotes integration over the angle between vector \(\boldsymbol{k}\)​ and \(\boldsymbol{p}\)​ for a fixed direction of \(\boldsymbol{k}\)​.7 In the above expression, we have suppressed the indices \(\gamma,\gamma'\)​ on each of the two form-factors for clarity, but the two \(\Lambda(...)\)​ terms carry opposite time-reversed labels.

The Eliashberg equation for the pairing gap in Eq. \(\ref{eq:sc-95gap-95equation}\) ignores a number of possibly important contributions. We ignore an interaction induced momentum-dependent (“Hartree”) renormalization of the MATBG dispersion, which has been argued to lead to qualitative changes in the bandstructure [48], [49]. We have (artificially) introduced a control parameter that selects a subset of the relevant processes, making the Eliashberg equation asymptotically exact in the large\(-N\)​ limit. With the possible qualitative changes in mind, in the analysis that follows we clearly identify and distinguish the features that follow from generic properties associated with moiré narrow-band systems and the ones that are tied to the specific aspects of the model and approximations used. This allows us to shed interesting light on the importance of different interactions on the origin of pairing and various qualitative aspects thereof in the limit where the Eliashberg framework is applicable.


Figure 2: Solution to the gap-equation in Eq. \(\ref{eq:sc-95gap-95equation}\) due to (a) phonon and (b) plasmon for \(\nu\approx 1\). In both cases we consider only \(0\)G processes giving \(T_C=0.0118\) meV and \(T_C=0.0072\) meV respectively. The phonon-induced solution is always positive, indicative of a purely attractive pairing, whilst the plasmon-induced solution has a sign-change associated with the gap — a characteristic behavior for an overall repulsive interaction (see main text). The solutions are normalized such that \({\rm max}\{\Delta(i\pi T,k)\}=1\). We intentionally show two different matsubara grids: linearly spaced in (a), and, variably spaced (exact matsubara at low frequencies and then approximate) in (b); see Appendix 7. A wide range of matsubara frequencies is required to capture the plasmonic-gorge (black arrow) in (b). Note that momenta, \(k\), are measured with respect to the center of the MBZ, see Fig. 1b.

To solve for \(T_c\)​ we seek the temperature, at which the kernel has an eigenvector with corresponding eigenvalue \(-1/T_c\)​. To perform the analysis numerically we must first perform the angular average Eq. \(\ref{eq:Kernel}\). We then need to perform two additional approximations: The first is to introduce frequency and momentum cutoffs and choose an appropriate mesh discretization. Furthermore, we select only a subset of momentum and frequency points under the sum. For more details on the numerical procedure we refer the reader to Appendix 7.

Before proceeding with our analysis of the results, let us pause to discuss some of the properties of the eigenvector in Eq. \(\ref{eq:sc-95gap-95equation}\). The optimal solution has a large negative weight where the interaction is most repulsive. Therefore, the optimal solution for a phonon vs. plasmon mechanism is qualitatively different; see Fig. 2a,b for the solutions to Eq. \(\ref{eq:sc-95gap-95equation}\) due to a phonon and plasmon solution, respectively. While the phonon solution is almost featureless at frequencies below the characteristic pairing energy (see later), it has a suppression of the solution at small \(k\)​. As discussed later, this is connected to a suppression of the superconducting dome at large values of the filling. On the other hand, the plasmon solution presents a sharp resonance-like feature at a characteristic pairing momentum (denoted the plasmon gorge) [50] as well as an eigenstate solution that changes sign — a necessary requirement to satisfy the gap equation \(\ref{eq:sc-95gap-95equation}\) when the Coulomb interaction is repulsive at all frequencies. The observed frequency dependencies, of the phonon solution in particular, leaves behind a resonant feature in the real-time Green’s function of single-fermion excitations, which can be detected in the single particle density of states.

4 Results

Let us now present our results for the superconducting transition temperature, obtained by solving the eigenvalue problem in Eq. \(\ref{eq:sc-95gap-95equation}\) numerically. In order to disentangle the effects of a purely phonon mediated attraction, including the effects of umklapp scattering (highlighted above), vs. the combination of phonon and plasmon mediated superconductivity, we study their effects individually in the next two sections. We also investigate various spectral features associated with the electron-phonon coupling and the effect of a metallic screening gate in subsequent sections. In the following analysis, we set the parameter \(N=20\)​ to ensure that we are firmly in the weak-coupling regime.9

4.1 Phonon mediated superconductivity

We start our discussion by focusing exclusively on the effects of the phonon mediated interaction, Eq. \(\ref{eq:V-95ph}\). In order to investigate the importance of phonons in a controlled fashion, we first set the Coulomb interaction (\({\cal V}^{\textrm{C}}_{\xi,\xi'}\)​) in Eq. \(\ref{eq:V-95C}\) to zero and compute \(T_c\)​ due to the phonon mediated interaction (\({\cal V}_{\xi,\xi'}^{\textrm{ph}}\)​) in Eq. \(\ref{eq:ham-95phonon}\). Moreover, we include a sequence of umklapp processes, \(m\mathrm G\)​, with \(m=0,1,2\)​ and \(3\)​; the results for the transition temperature (\(T_c\)​) as a function of \(\nu\)​ are shown in Fig. 3a.

As is evident from our results, umklapp scattering processes up to \(m=2\)​ have a dramatic effect on \(T_c\)​, which saturates for \(m\geq3\)​. These umklapp phonons are essentially the lowest optical modes resulting from the folding of the original acoustic phonon branch into the MBZ. Thus, we find that the interaction with these lowest optical modes are crucial for understanding the electronic properties of MATBG. In order to clearly demonstrate the effect of umklapp processes on the pairing interaction, it is useful to study the phonon spectral function as a function of energy, defined as \[\begin{gather} {\alpha}^2 F({\omega}) &=& {N_0(0) g \over 2N}\times\\\nonumber &&\left\langle \sum_{\boldsymbol{q}} \omega_{\textrm{ph}}(q) \lvert\Lambda( \boldsymbol{k} + \boldsymbol{{q}} ,\boldsymbol{k})\rvert^2 ~{\delta}({\omega}- \omega_{\textrm{ph}}(q))\right\rangle_{\textrm{FS}},\end{gather}\]

where \(\boldsymbol{q}=\boldsymbol{\tilde{q}}+\boldsymbol{G}\)​, as defined previously (see paragraph following Eq. \(\ref{eq:V-95ph}\)). In the above, \(N_0(0)\)​ is the density of states at the Fermi surface (FS) and \(\langle \dots \rangle_{\textrm{FS}}\)​ denotes averaging \(\boldsymbol{k}\)​ over the FS. We stress the presence of the form-factors, \(\Lambda(\boldsymbol{p}, \boldsymbol{k})\)​, which encode the unusual dependence of the superconducting kernel on the large momentum umklapp processes. As can be seen in Fig. 3b, the inclusion of umklapp processes have a clear effect on the spectral function, most prominent of which include a significant enhancement in the electron-phonon coupling strength, combined with a drastic rearrangement of the spectral weight to higher energies. These same umklapp processes are responsible for increasing \(T_c\)​.

Figure 3: (a) Superconducting dome due to a phonon-only mechanism (\(N=20\)). Colored lines correspond to a different number of umklapp processes (\(m\)G with \(m=0,1,2,3\)) involved in the pairing interaction. The \(T_C\) values for \(\nu >0\) are higher than those for \(\nu<0\) as the conduction band, Fig. 1b, of the model in Ref. [45] is flatter than the valence band. The \(0\)G dome peaks near the van-Hove singularities of the bands. Saturation of the \(T_C\) increase with \(3\)G umklapp processes can be traced back to the nature of the Bloch wavefunctions (see main text and details of the bandstructure discussed in Appendix 6). (b) Electron-phonon spectral function for \(\nu\approx 1\) showing the importance of umklapp processes. (c) A plot of \(T_C\) obtained from a BCS-type formula, Eq. \(\ref{eq:BCS-95formula}\), with parameters obtained from the MATBG spectral function, Eq. ¿eq:spectral-95function?. As expected, the BCS dome peaks in the vicinity of the van-Hove singularity, which by virtue of the specific details of the model in Ref.  [45], occurs at \(\nu\approx 0.63\). Note the similar trend of \(T_C\) enhancement between panels (a) and (c) with umklapp processes.

The characteristic energy scale associated with the peak of the above spectral function, provided that the form-factors do not lead to a suppression, is set by the graphene sound velocity (\(c_s\)​) and the length-scales stemming from the moiré period. Using the magnitude of the moiré reciprocal vector \(G=4\pi/(\sqrt{3}L_M)\)​ as the characteristic size of the MBZ (\(L_M = \frac{a}{2}\sin(\frac{\theta}{2})\)​), we find the typical frequencies as

\[\begin{equation}\frac{1}{2} c_s (m+1) G = c_s (m+1) \frac{2 \pi}{\sqrt{3} L_M} \approx 2.1\,, 4.3\,, 6.4 ~{\rm meV}\end{equation}\]

for the first three \(m\)​G umklapp processes: \(0\)​G, \(1\)​G, \(2\)​G. These estimates are in reasonably good agreement with the location of the peaks in Fig. 3b. In MATBG, as mentioned in the introduction, the precise values of these characteristic frequency scales can be affected by the details of the phonon dispersion (including, for example, the presence of a gap [22], [23]).

The microscopic mechanism responsible for the significant contribution of umklapp processes to the spectral function, which leads to the enhancement of superconducting \(T_C\)​, is intimately tied to the origins of the flat moiré bands. To obtain the electronic dispersion, \(E_{\boldsymbol{k}}\)​, of carriers in these narrow-bands (see Fig. 1b) due to the slowly varying moiré interlayer potential, it is not sufficient to consider only plane wave Bloch states with crystal momentum \(\boldsymbol{k}\)​, but also those of nearby states that are coupled by multiple of the moiré reciprocal vectors, \(\boldsymbol{G}\)​ (see Appendix 6). As a result, the spectral weight of the resulting Bloch wavefunction is extended across several plane wave states. This results in a slowly vanishing Bloch wavefunction overlap, on the scale of the moiré reciprocal momentum scale \(G\)​, that enters into the form-factors, \(\Lambda_{\gamma\gamma'}(\boldsymbol{p},\boldsymbol{k})\)​. We stress that this property of narrow-band wavefunctions is independent of the finer details of the bandstructure and is intimately tied to many of the unusual properties of MATBG (e.g. it is partly responsible for the extended Wannier functions [45], [52]).

Let us now place the essence of the phonon umklapp-driven enhancement of the critical temperature, \(T_C\)​ in the context of BCS theory. We note upfront, however, that although it captures some of the trends seen in Fig. 4a, it is by no means a replacement for the full Eliashberg approach (see later). Within standard BCS theory with a coupling constant, \(g\)​, the expected transition temperature is given by [53],

\[\begin{equation}T_C \approx 1.14 \langle\omega\rangle ~e^{-1/\tilde{\lambda}}\label{eq:BCS-95formula}\end{equation}\]

where \(\langle\omega\rangle\)​ corresponds to a pairing energy range (usually the Debye frequency) and \(\tilde{\lambda}=g N_0(0)\)​. We can now use the spectral function defined in Eq. ¿eq:spectral-95function? to recast the effective coupling constant as [54]

\[\begin{equation}\tilde{\lambda} = 2 \int d \omega ~\frac{{\alpha}^2 F({\omega})}{\omega},\end{equation}\]

and similarly denote the characteristic pairing frequency \({\langle}\omega{\rangle}\)​ in terms of the moment of the same distribution as

\[\begin{equation}{\langle}\omega{\rangle}= \frac{2}{\tilde{\lambda}} \int d\omega ~{\alpha}^2 F({\omega})\,.\end{equation}\]

For the MATBG phonon spectral function plotted in Fig. 3b, corresponding to an electron filling near \(\nu\approx 1\)​, we list each of the two parameters, \(\tilde{\lambda}\)​ and \({\langle}\omega{\rangle}\)​, in the legend for all of the umklapp processes (\(m\)​G, with \(m=0-3\)​). A simple application of the BCS formula in Eq. \(\ref{eq:BCS-95formula}\) in terms of the two filling dependent parameters, \(\tilde{\lambda}\)​ and \({\langle}\omega{\rangle}\)​, captures many of the interesting trends that we saw previously, as shown in Fig. 3c, including the relevance of umklapp processes.

The simple BCS formula, parametrized by \(\tilde{\lambda}=g N_0(0)\)​, predicts the \(T_C\)​ dome to peak at the location of the van Hove singularity. To a large extent, this is also reproduced in the domes of Fig.3a; we note, however, that the peak of the dome shifts towards higher fillings upon including successive umklapp processes.10

We also note that the simple BCS estimate, where the characteristic pairing energy is naively given by \({\langle}\omega{\rangle}\)​, grossly overestimates the \(T_C\)​ as seen by comparing Fig. 3b and Fig. 3c. In fact the predicted \(T_C\)​ lies beyond the regime of validity of the above formula, e.g. \(T_C\gg |\mu|\)​. This is all to be expected given the overall strength of the coupling \(\tilde{\lambda}\)​, which can introduce further corrections in Eq. \(\ref{eq:BCS-95formula}\) and suppress \(T_C\)[54], [55]. The umklapp-driven enhancement of phonon-mediated superconductivity is not model specific and is a generic property of any moiré narrow-band system with a Bloch wavefunction overlap that is vanishing slowly, on a moiré momentum scale, as a function of the momentum exchange, \(\boldsymbol{q}\)​. In Sec. 4.4, we shall return to a discussion of the electron-phonon spectral function and propose a sharp experimental signature that can probe some of these features in a tunneling measurement.

4.2 Plasmon mediated superconductivity

We are now in a position to include the effects of Coulomb interactions on pairing. To that end, we begin by evaluating the dielectric function \({\varepsilon}_{\textrm{RPA}}(i{\omega},\boldsymbol{q})\)​ numerically [24] and reinstating \({\cal V}^{\textrm{C}}_{\xi,\xi'}\)​ to the pairing kernel in Eq. \(\ref{eq:Kernel}\).

Figure 4: Comparison of phonon and plasmon mediated pairing mechanisms as a function of an increasing number of umklapp processes. (a) \(0\)G (no umklapp), (b) \(1\)G and (c) \(2\)G. Blue curves correspond to a phonon only mechanism (Sec. 4.1), red curves correspond to a plasmon only mechanism (Sec. 4.2), and, yellow curves include effects of both phonon and plasmon on pairing (see Sec. 4.2). The phonon-driven attraction is strongly enhanced with the inclusion of umklapp scatterings; the plasmon mechanism is largely insensitive to umklapp. We choose \(N=20\).

Let us begin by studying the problem in the absence of the phonon mediated interaction, i.e. set \({\cal V}_{\xi,\xi'}^{\textrm{ph}}=0\)​, and include only the effects of the Coulomb interaction. In this case, the origin of pairing lies in the frequency dependence of the dielectric function close to the plasma resonance [50], [56]. Moreover, just as in the case of phonons, we include a sequence of \(m\)​G umklapp processes for the plasmons (we choose \(N=20\)​). The result for \(T_c\)​, due purely to the plasmonic mechanism, is shown in red in Fig. 4(a)-(c) as a function of filling. This analysis leads us to conclude that (i) the plasmonic mechanism of pairing leads to an enhancement of \(T_c\)​ for a narrow range of fillings near \(\nu \approx 2-3\)​, and, (ii) successive \(m>0\)​ umklapp processes have no appreciable effect on \(T_c\)​.

We now explain the microscopic origin for both of these observations. In contrast to a purely attractive phonon-mediated interaction, the dynamically screened Coulomb interaction is always repulsive. However, as argued in Sec. 3, the frequency dependent dynamically screened interaction becomes weak enough at certain momenta, such that the sign-changing gap function can minimize the overall effect of repulsive Coulomb interaction [57], giving rise to an effective attractive part.

In order to understand the structure of the frequency-momentum regions where the screened interaction becomes weak, it is helpful to plot \({\cal V}^{\textrm{C}}_{\xi,\xi'}(\boldsymbol{q},i{\omega})\)​ in Eq. \(\ref{eq:V-95C}\), as shown in Fig. 5a at a fixed filling (\(\nu\approx1\)​). The behavior differs from the one in a conventional 2D Fermi gas [58]. Most crucially, in the latter, given the form of the polarization function, one would expect the interaction to be most repulsive as \(q\to 0\)​ and then weaken with increasing \(q\)​. On the other hand, in MATBG, there is a local minimum of the interaction at a finite momentum \(q\)​. We can understand this behavior by focusing on the limit of \(\omega\to 0\)​, as discussed in Ref. [59], [60]. At the magic angle and at low fillings \(\nu\approx 0\)​, the static polarization function behaves as \(\Pi_{ee}(\boldsymbol{q},\omega\to 0) \propto q/v_F\)​. This form is reminiscent of the polarization function in a Dirac-like system [61], [62]. At momenta smaller than the moiré reciprocal momentum (\(G\)​), \(v_F\)​ corresponds to the renormalized MATBG Fermi-velocity near \(\nu\approx0\)​ with \(v_F\sim 10^4\)​ m/s. The polarization function is dominated by the inter flat-band transitions. On the other hand, at momenta comparable to and larger than the moiré scale (and at a similar filling), \(v_F\sim 10^6\)​ m/s, the bare graphene velocity. The polarization function is now dominated by the inter-band transitions between the flat and dispersive bands as the effect of the moiré interlayer potential becomes less relevant. As a result, the plasmons in MATBG [24] have interesting properties. As long as plasmons rise above the particle-hole continuum, its dispersion is controlled by the energy scale associated with inter-band transitions between flat-bands, and between the flat and the dispersive bands. As such, it therefore becomes weakly sensitive to the filling value.

The aforementioned local minimum of the pairing interaction occurs at momenta smaller that the moiré reciprocal lattice scale, \(G\)​. At momenta larger than \(G\)​, the interaction in Eq. \(\ref{eq:V-95C}\) reduces to the simple unscreened form, \(2\pi e^2/q\)​, suppressing any dynamical contribution to the superconducting gap. As a result, any contribution due to higher umklapp processes involving plasmons does not enhance \(T_C\)​ drastically; see Fig. 4a-b. In fact even at large enough momenta, on the scale of \(2\)​G, as the form-factors are still non-vanishing it can suppress it slightly; see Fig. 4c.

Figure 5: (a) Plot of \(-\log\left[N_0(0) {\cal V}^{\mathrm{C}}_{\xi,\xi'}(\boldsymbol{q},i{\omega})\right]\) for \(\nu\approx 1\). Note that a local minimum (black arrow) of the interaction occurs at a finite momentum (see discussion in the main text). (b) Fermi-surface averaged dependence of form-factors, Eq. \(\ref{eq:coherence-95factor-95overlap}\), with the same vertical axis as in panel (a). Note how small momentum (\(q < K\)) processes dominate over the large momentum terms (\(q > K\)). Near a filling, \(\nu\approx 2.45\), the form-factors are slightly larger for a wider range of momenta than at lower fillings. This behavior lies behind the narrow plasmon peak of the \(T_C\) dome, Fig. 4, which is exponentially sensitive to the coupling strength, c.f. Eq. \(\ref{eq:BCS-95formula}\). (c) Plasmon-mediated SC under different approximation schemes. For the plasmon-pole approximations (iii,iv) we use \(\omega_{\mathrm{pl}}\approx 6\) meV. There is hardly any difference between (i) and the original result. (d) Superconducting temperature in the plasmon pole approximation as a function of plasma frequency \(\omega_{\mathrm{pl}}\). Pairing occurs from frequencies close to \(\omega_{\mathrm{pl}}\), and leads to an increase in \(T_C\) as \(\omega_{\mathrm{pl}}\) approaches the chemical potential.

While the specific form of the screened interaction, Fig. 5a, affects the shape of the plasmon-induced superconducting dome, it is useful to disentangle it from the role played by other elements that appear in the pairing kernel, Eq.\(\ref{eq:Kernel}\): the form-factors \(\Lambda_{\gamma\gamma'}(...)\)​ and the energy denominator, \(1/(\nu^2+E^2_{\xi,\boldsymbol{p}})\)​. To that end, we focus on the \(0\)​G processes and selectively modify the different elements that enter into the kernel, Eq. \(\ref{eq:Kernel}\). In particular, we consider the following modifications: (i) we compute the dielectric function at one specific filling, \(\nu\)​, and use it for all fillings (thereby ignoring the \(\nu-\)​dependence of the dielectric function), (ii) we ignore the momentum dependence of the form-factors and replace them with unity, except for restricting the interaction to intervalley pairing, (iii) we use the “plasmon-pole” approximation instead of using the full RPA dielectric function, and, (iv) we invoke the same approximation as in (iii) above, but make the substitution for the form-factors as in (ii). A figure demonstrating all of these cases is shown in Fig. 5c. Notably, the results under the approximation in (i) above do not affect the results at all (not shown), thereby indicating that the density dependence of the dielectric function does not play a significant role.

Let us begin with a discussion of points (iii) and (iv) above, that rely on the plasmon-pole approximation. As described previously, the Coulomb interaction is always repulsive; the dynamical screening can, however, lead to a possible superconducting solution that overcomes the effect of the repulsion. To investigate this matter further, we explore here an idealized limit — the so-called “plasmon-pole” approximation— where the Coulomb interaction in Eq. \(\ref{eq:V-95C}\) is replaced by

\[\begin{equation}{\cal V}^{\textrm{C}}_{\xi,\xi'}(\boldsymbol{q},i{\omega}) \approx {1\over N}{2\pi e^2\over \kappa q} \left(1-\frac{\omega^2_{\textrm{pl}}}{\omega^2_{\textrm{pl}}+\omega^2}\right)\,.\end{equation}\]

Here, \(\omega_{\textrm{pl}}\)​ is the plasmon frequency. We note that the “plasmon-pole” approximation is strictly valid only in the \(\omega\gg v_Fq\)​ limit, which will be of interest to us below. For \(\omega \ll v_F q\)​, the interaction can be approximated by its purely repulsive and static Thomas-Fermi screened form. For a wide range of momenta, the plasmons of interest to us originate primarily from inter-band transitions between the nearly-flat and dispersive bands. As argued above, the resulting plasma frequency is independent of the filling and is instead set by the bandwidth, \(W\)​, of the nearly-flat band and the gap, \(\Delta_{\textrm{band}}\)​, between flat and dispersive bands as \(\omega_{\textrm{pl}}\approx\sqrt{W \Delta_{\textrm{band}}}\)​ (see Ref. [24] for details). We therefore choose a constant \(\omega_{\textrm{pl}}\)​ in our calculation. As is evident from Fig. 5b, the plasmon-pole result exceeds the RPA result. With increasing filling, \(T_C\)​ rises almost linearly as the Fermi energy approaches the plasma frequency, thereby enhancing the effect of the plasmon in driving pairing. However, for \(\nu\gtrsim 3\)​, \(T_C\)​ starts to drop rapidly — a behavior we attribute to the bands being more dispersive at these fillings, c.f. Fig. 1b.

We can verify our understanding of the interplay of these two results by varying \(\omega_{\textrm{pl}}\)​ as an external phenomenological parameter, as done in Fig. 5c. We notice immediately that the closer \(\omega_{\textrm{pl}}\)​ is to the relevant chemical potential, the higher is the \(T_C\)​. This is however not sufficient at large fillings, \(\nu > 3\)​. We note that the sharp fall-off at \(\nu>3\)​ is a property of the continuum model used. This fact is precisely the reason for the form of the gap solutions obtained in Fig. 2, where for \(k\lesssim0.25K\)​ the solution due to phonons vanishes (Fig. 2a), while it is positive in Fig. 2b, implying an overall repulsive contribution to the plasmon-induced gap equation. This analysis also suggests an interesting possible route towards enhancing \(T_C\)​ due to plasmons. If \(\omega_{\textrm{pl}}\)​ can be brought closer to the chemical potential, whilst maintaining the same strength of Coulomb interactions in a system, then it is possible to raise \(T_C\)​.

We conclude the analysis of the plasmon-pole approximation and its effect on pairing by pointing out that suppressing the momentum dependence of the form-factors leads to an enhancement of \(T_C\)​; see curves labeled (ii) and (iv) in Fig. 5a. This can be understood by realizing that the Bloch wavefunction overlap depends on the underlying “fidget-spinner” structure of the energy contours, c.f. Fig.1c, which can suppress certain scattering processes and thereby lower \(T_C\)​.

The plasmon-pole approximation captures most of the features we obtain within the full Eliashberg calculation (Fig. 4). However, it does not immediately lead to a simple explanation for the sharp peak associated with the (plasmon-)dome near \(\nu \approx 2-3\)​. It is natural to ask if this feature arises solely due to a change in the dielectric function as a function of \(\nu\)​. To explore this possibility, we compute \(T_C\)​ by fixing the dielectric function corresponding to \(\nu\approx 2.45\)​ (associated with the peak of the plasmon-dome in Fig. 4a) and not varying it as a function of \(\nu\)​; this corresponds to the approximation denoted as (i) above. We find that the overall shape of the resulting dome is completely identical to the full computation (result not shown). As explained previously, this behavior stems from the dielectric function of MATBG being dominated by inter-band transitions, which are largely insensitive to the filling.

Taking all of these observations into account, the peak of the dome at a filling of \(\nu \approx 2-3\)​ comes from an interplay of the form-factors along with the dielectric function. This conclusion stems from the analysis leading to Fig. 5b, which demonstrates that the plasmon-pole approximation, even with the appropriate form-factors, can not reproduce the sharp peak at \(\nu\approx 2-3\)​.

To assess this further, we now focus on the \(q\)​ dependence of the form-factors. To that end we plot in Fig. 5d, the following quantity

\[\begin{equation}F(q,\mu) = \int_{\textrm{MBZ}} d^2 \boldsymbol{k} \int d \theta_{\boldsymbol{k}'}~ k' \lvert \Lambda(\boldsymbol{k}',\boldsymbol{k}) \rvert^2 \delta(E_{\xi,\boldsymbol{k}}) \delta(E_{\xi,\boldsymbol{k'}}) \,,\label{eq:coherence-95factor-95overlap}\end{equation}\] where \(q=|\boldsymbol{k}-\boldsymbol{k'}|\)​ and \(\delta(E_{\xi,\boldsymbol{k}})\)​, \(\delta(E_{\xi,\boldsymbol{k'}})\)​ constrain the two states to lie on the Fermi-surface (within mesh resolution); \(\theta_{\boldsymbol{k}'}\)​ denotes angle of \(\boldsymbol{k}'\)​. The decrease in \(F(q,\mu)\)​ with \(q\)​ reflects the density dependence of the orbital hybridization of the Bloch wave-functions in the relevant bands. For comparison, \(F(q,\mu)\)​ would be momentum independent and equal to unity, when there is no orbital hybridization (i.e. when the band and orbital basis are identical). In contrast, it is known to diminish rapidly with the momentum exchange, \(q\)​, in a Dirac system [63]. Therefore, \(F(q,\mu)\)​ extends to higher momentum close to \(\nu \approx 2-3\)​, where orbital hybridization is minimized, taking a maximum at \(\nu\approx 2.45\)​. In turn, this maximizes the product with the screened interaction from, Fig.5a, leading to a higher \(T_C\)​.

Finally, we analyze the cooperative effect of both phonons and plasmons on pairing; the resulting \(T_C\)​ vs. \(\nu\)​ is shown for the combined (as well as individual) effect of phonons and plasmons in Figs. 4(a)-(c). As before, we include the effects of \(m\)​G umklapp processes with \(m=0-2\)​, which has a strong effect on the phonon-mediated contribution but barely affects the plasmonic mechanism. We find that with the inclusion of upto \(2\)​G processes, the umklapp-driven phonon-mediated mechanism clearly dominates over the plasmon mechanism.11 However, both mechanisms work in a cooperative fashion to give rise to pairing.

4.3 Role of external screening

Our analysis thus far has shed light on the distinct features associated with a purely electronic vs. a purely phonon-based mechanism, as well as their combined effect, on the emergence of SC in MATBG. In light of the recent experiments [25], [26], [28] that have studied the role of an external screening layer on the filling dependence of \(T_c\)​, let us explore the effect of similar screening within our setup.

We begin by studying the effects of a metallic gate, coupled to MATBG via Coulomb interactions, on pairing. The only role played by the gate is in a further renormalization of the Coulomb interaction, Eq. \(\ref{eq:V-95C}\), via the dielectric constant, \({\varepsilon}_{\textrm{RPA}}\)​. In the limit where the density of states associated with the metallic gate is higher than that of the screened substrate, the effect of the gate can be incorporated by modifying the bare Coulomb interaction as:

\[\begin{equation}\frac{2\pi e^2}{q}\to \frac{2\pi e^2}{q}\left(1-e^{-2q d}\right)\,,\label{eq:gate-95screening}\end{equation}\]

where \(d\)​ is the distance between MATBG and the gate. As a result, the dynamically screened interaction (Eq. \(\ref{eq:V-95C}\)) now changes to \[\begin{gather} {\cal V}^{\textrm{C}}_{\xi,\xi'}(\boldsymbol{q},i{\omega}) &= {1\over N}{2\pi e^2\over {\varepsilon}_{\mathrm{RPA}}(\boldsymbol{q} ,i{\omega}) q} \left(1-e^{-2q d}\right),\\ \nonumber \varepsilon_{\mathrm{RPA}}(\boldsymbol{q} ,i{\omega}) &= \kappa-\frac{2\pi e^2}{q} \Pi_{ee}(\boldsymbol{q},i\omega)\left(1-e^{-2q d}\right)\,.\end{gather}\]

The presence of the metallic gate results in an overall suppression of the Coulomb interaction, that scales exponentially with \(d\)​. Thus, processes involving \(q \lesssim 1/2d\)​ are not responsible for mediating SC.

Figure 6: (a) Effect of screening by a metallic gate for four different values of \(d\) (nm) on \(T_C\) with the inclusion of \(2\)G umklapp processes (\(N=20\)). Due to the cooperative interplay of phonons and plasmons (see main text) we find that the gate suppresses \(T_C\) for a range of fillings. (b) Effect of screening for few fillings over a wide range of gate distance. The dashed-dotted lines correspond to a purely phonon-driven mechanism, while the dashed lines include the cooperative effects due to both plasmons and phonons (see also Fig. 4c). Note the non-monotonic dependence of \(T_C\) on \(d\) as well as the wide range of \(d\) values over which \(T_C\) varies; see discussion in the main text for implications on plasmon-mediated pairing. (c) Normalized differential conductance from Eq.\(\ref{eq:diff-95cond}\) computed at \(T = T_C/10\) for each curve for phonon-mediated superconductivity. All umklapp processes result in conventional BCS gap behavior and exhibit resonances corresponding to relevant \(m\)G umklapp processes (indicated with arrows), c.f. Fig.3b. This modulation of \(\tilde{G}(\omega-\mu)\) allows to reconstruct the form of the superconducting pairing. Dashed line indicates the electronic bandwidth, which by virtue of the model, gives rise to a signal due to \(2\)G, \(3\)G phonons inside the gap between the flat and non-flat band, respectively (shaded area). As a result of mesh induced broadening, \(N_0(\omega-\mu)\) has a small finite value, and so does \(\tilde{G}\), inside the gap region.

In Fig. 6a, we plot \(T_C\)​ as a function of filling (\(\nu>0\)​) for a few different values of \(d\)​ (with the \(2\)​G umklapp processes included). As discussed previously in the context of Fig. 4, the dynamically screened Coulomb interaction assists the phonon-driven mechanism in a cooperative fashion, especially for momenta smaller than the characteristic scale of the MBZ (see Fig. 5a). A screening of the form in Eq. \(\ref{eq:gate-95screening}\) is expected to result in a strong suppression of \(T_C\)​ for a wide range of fillings, but especially so in the range of fillings where the plasmon contributes the most to pairing. This behavior is indeed observed in Fig. 6a, where the gate leads to a suppression of superconductivity over the entire range of fillings, with perhaps a slightly more pronounced effect near fillings \(\nu\approx 2-3\)​.

In addition to the direct experimental relevance [25], [26] for understanding the effect of screening from a metallic gate on the strength of Coulomb interaction in Eq.\(\ref{eq:gate-95screening}\), the above setup also helps elucidate the length-scales that participate in plasmon mediated pairing. Specifically, a metallic gate located a distance \(d\)​ away suppresses the interaction over length-scales of \(O(2d)\)​. Depending on how this length-scale compares with the other characteristic length-scales in the problem, which include e.g. the scale associated with twist-angle inhomogenities etc., will determine whether electron-electron interactions contribute towards pairing. This requires a careful analysis of the necessary length-scales associated with the plasmons that are required for plasmon-mediated pairing to be observable experimentally.

To address the question raised above, we study the effect of varying \(d\)​ over many orders of magnitude in Fig. 6b. For \(d\to \infty\)​ (i.e. significantly larger than any other relevant length-scale in the problem), we expect to reproduce the previously obtained values for \(T_C\)​ in Fig. 4c for both the phonon and plasmon mediated interaction. On the other hand as \(d\to 0\)​, we expect the contribution due to the plasmon to be suppressed significantly, such that \(T_C\)​ will almost entirely be determined by only the phonon-mediated interaction. We observe that while low-momentum processes corresponding to distances \(2d \approx 200~ {\rm nm} \approx 14.5 L_M\)​ are necessary for the plasmons to maximally assist in pairing, the plasmon substantially assists the phonon-driven mechanism already at lengthscales starting \(2d \approx 50~ {\rm nm} \approx 3.5 L_M\)​.

Interestingly, we also observe that \(T_C\)​ depends non-monotonically on the gate-distance, \(d\)​, as shown in Fig. 6c. As noted above, bringing the screening layer closer to MATBG starting from a very large separation, leads to a drop in the value of \(T_C\)​ as a result of the suppression of plasmon-induced pairing. However, once \(d\approx 10\)​ nm, the value for \(T_C\)​ starts to drop beyond what is expected even from a purely phonon-mediated mechanism (i.e. \(d\to 0\)​), eventually reaching a minimum value at \(d\approx 1\)​ nm. The additional reduction of \(T_C\)​ arises as a result of the metallic gate suppressing the attractive contribution to the pairing coming from momenta \(q \ll K\sim{G}/2\)​, c.f. Fig.5, leaving only the purely repulsive part. The minimum value for \(T_C\)​ is attained approximately at the same value of \(d\)​ at all fillings; this can be understood from the filling independence of the dielectric function for MATBG, which is dominated by the interband transitions (see Sec. 4.2). Beyond this value of \(d\)​, decreasing it further leads to an increase in \(T_C\)​, in agreement with the expectations in a purely phonon-mediated attraction [54].

As mentioned in the introduction, experimental setups that include a gate to control the strength of the Coulomb repulsion have already been realized in several groups. Namely, Refs. [25], [26], [28] find that the gate tends to suppress the correlated state significantly, while superconductivity remains largely unaffected. As such these results are broadly consistent with our conclusions regarding superconductivity within our model. In these experiments, the gate (or, the screening substrate) is closer to MATBG than \(30\)​ nm and thus our model would predict superconductivity to be primarily phonon-driven. We note however that the precise onset of the non-monotonic gate dependence seen in Fig. 6b is dependent on the finer details of the screening mechanism. For example, the extent to which the density of states for the metallic gate compares with that of MATBG, or, whether further modifications to the RPA dielectric function (such as local field effects) can become pronounced for \(q\)​ greater than moiré momentum lengthscale can affect these details. Finally, we also note that this filling independent onset of the non-monotonic mechanism is indicative of a geometric scale dictating properties of a dielectric function — this corresponds to the interband transitions that are set by the flat-band bandwidth, \(W\)​. To that end, if flat-bands were to become broader, e.g. due to the presence of Hartree corrections, we would then expect a dependence on filling to appear.

4.4 Spectroscopic signature of umklapp phonons

Let us finally discuss a sharp experimental fingerprint associated with the phonon umklapp processes and their role in pairing. In Sec. 4.1, we have demonstrated how the scattering of umklapp phonons contributes significantly to pairing. Consequently, we anticipate these modes to appear as resonant features in the single-particle density of states at the energies where the spectral weight of the phonons is pronounced [64] (see Fig. 3b). Consider the normalized differential conductance, defined as \[\begin{gather} \tilde{G}(\omega-\mu) = \frac{N_{\textrm{SC}}(\omega-\mu)}{N_0(\omega-\mu)},\label{eq:diff-95cond}\end{gather}\]

where \(N_{\textrm{SC}}(\omega-\mu)\)​ is the density of states in the superconducting state and \(N_0(\omega-\mu)\)​ is the corresponding density of states in the normal state; the result is shown in Fig. 6c. To obtain the superconducting density of states, we extended the gap equation, Eq. \(\ref{eq:sc-95gap-95equation}\), to its complete non-linear form. This allows us to calculate not only the critical temperature \(T_C\)​, but also the actual superconducting gap. For more details on the specific computational aspects associated with solving the non-linear gap equation as well as the details of the analytic continuation to real frequencies necessary for calculation of the spectral function, see Appendices 7 and 8. For all \(m\)​G processes, we see a well-defined \(s\)​-wave superconducting behavior with the usual coherence-peak at the onset of the pairing gap; the higher the \(T_C\)​, the larger is the gap and hence the onset of the peak. Therefore, we observe the coherence peak shifting in agreement with umklapp processes enhancing the \(T_C\)​. Away from the peak, \(\tilde{G}(\omega-\mu)\)​ exhibits the classic square root singularity until it reaches the \(0\)​G and \(1\)​G phonon-resonance features (indicated with the first two, left-most arrows). We discuss now the final signature seen in the differential conductance.

We observe sharp “Fano-shaped” resonance-like features at frequencies of the order of 5 meV, which corresponds precisely to the energy where the density of states deduced from \(\alpha^2F(\omega)\)​ for the \(2\)​G and \(3\)​G umklapp phonon processes is the largest, as seen earlier in Fig. 3b. The origin of these features lies in the inelastic tunneling process into the bottom of the Bogoliubov quasi-particle spectrum that also involves the production of a phonon excitation [64], [65]. We stress however that unlike the role of umklapp processes in enhancing superconductivity, which are model independent and stem purely from the underlying moiré physics (see Sec. 4.1), the presence of sharp Fano resonances located in the bandgap between nearly flat and dispersive bands relies on the narrow bandwidth of the former (\(W\approx 4\)​ meV). If instead the bandwidth associated with the nearly flat band were larger than the characteristic phonon pairing frequency, the corresponding resonance would appear inside the band.12 Thus, when the resonance occurs inside the band gap, it becomes more clearly visible in the proposed experimental setup.

5 Discussion and Outlook

In this work, we have focused on the pairing instabilities in MATBG within the framework of Eliashberg approximation. We have focused on two sources of interaction, mediated by acoustic phonons of the original graphene layers, and the dynamically screened Coulomb repulsion, respectively. Interestingly, we find that both sources contribute to pairing in a cooperative fashion with a relatively similar strength, Fig. 4. However, the screened Coulomb repulsion plays a role in a narrow range of density, while the phonons contribute to pairing over a wide range of fillings that nearly extends over the entire narrow band. One of the key findings of our work is the prominent role played by the umklapp processes in the phonon mediated interaction in enhancing the transition temperature, \(T_c\)​. Umklapp processes that scatter states up to the third MBZ (“\(m\)​G” with \(m=3\)​) have a marked effect on the enhancement of \(T_c\)​ and eventually saturate for any higher order processes (\(m>3\)​), c.f. Fig.3a. On the other hand, the dynamics associated with the screened Coulomb repulsion, while being relatively insensitive to the umklapp processes, plays a cooperative role in pairing for a wide range of fillings (i.e. by aiding the phonon-mediated mechanism), but in particular in the vicinity of \(\nu \approx \pm 2-3\)​. As a result, the superconducting \(T_C\)​ exhibits a non-monotonic dependence as a function of the distance to a nearby metallic (screening) gate. The last observation is dependent on both the properties of the screening material and the extent to which long-wavelength plasmons can participate in SC pairing.

It is natural to ask if there are sharp experimental signatures associated with any of the scattering processes considered above in MATBG. We have demonstrated that the electron-phonon interaction leaves behind an imprint in a spectroscopic tunneling experiment. In particular, the superconducting coherence peaks are replicated at frequencies that are resonant with the frequencies of the corresponding bosonic modes that are responsible for superconductivity. Thus, specifically for the umklapp phonon processes (which correspond to optical-like modes arising from a folding of the original acoustic branch in the MBZ), the resonant features are present at frequencies that are independent of the details associated with the electronic band structure. The replica features are likely to be resolved most clearly when the resonant frequencies occur at energies that lie inside the bandgap between the nearly flat and dispersive electronic bands associated with MATBG.

Whenever possible, we have explicitly pointed out the universal, model-independent aspects of our predictions, in contrast to the features that rely on the non-universal aspects of the model. We stress that the importance of phonon-umklapp processes in enhancing the pairing temperature is a model-independent feature. It relies solely on the geometric properties of the MBZ as well as on the presence of a slowly varying moiré interlayer potential that gives rise to slowly vanishing form-factors (Bloch wavefunction overlaps). As such, we therefore expect to see features in the differential conductance at the phonon resonance frequencies, which may or may not lie inside the bandgap — the latter property being model specific. Likewise, for the plasmon-mediated mechanism for superconductivity, we can identify various universal features. Firstly, plasmons that originate from the flat bands necessary have a characteristic frequency that is of scale similar to the flat-band bandwidth (or, chemical potential). Therefore, the dynamical Coulomb screening is bound to play an important role in describing superconductivity, i.e. the effect of Coulomb interaction cannot be simply reduced to a static repulsion. Secondly, as a result of the microscopic origin of the behavior of the polarization function (i.e. interband transitions between flat-bands at small momenta and between flat/dispersive bands at large momenta), it is evident that the plasmon-mediated mechanism is weakly affected by umklapp processes. The extent to which the plasmon dome however exhibits a narrow peak at high fillings \(\nu\approx 2-3\)​ is dependent on the extent to which the continuum model [45] accurately captures behavior of the form-factors in MATBG.

There are a number of interesting directions that remain to be explored, based on the formalism developed here. A natural extension would be to include the filling-dependent bandstructure renormalization arising from the interaction itself [48], [49] and then study the pairing instabilities as a result of the same interactions. While we have considered the effects of a single acoustic phonon mode on pairing in this study, MATBG hosts multiple in-plane as well as out-of-plane phonon modes [22], [23]; some of these modes also develop small gaps as a result of folding due to the moiré potential. Studying the combined effects of these modes in the presence of umklapp scattering on pairing and on the resonant tunneling spectra is clearly an interesting problem.

Recent theoretical works have indicated the possibility of low-energy goldstone modes associated with spontaneously broken continuous symmetries near some of the correlated insulating phases at commensurate fillings [17], [40], [66]. It would be interesting to explore and clarify to what extent the umklapp processes considered in this work in the context of phonons are also important for these goldstone modes of electronic origin in general — the recent analysis in Ref. [66] includes some such (“\(2\)​G”) processes for one specific example. It would also be interesting to study the effect of screening gates on the associated phenomenology. With regard to the relevance of the umklapp processes considered here on other aspects of the phenomenology, it is likely that they also play an important role in electrical charge transport. To what extent these might be relevant for some of the recent unconventional results reported in Refs. [67], [68] is an interesting open question.

Another interesting question concerns the unavoidable role of disorder, on the superconducting phenomenology in MATBG [69]–[71]. Disorder can significantly affect the ability of the flat band to screen the Coulomb interaction, which will modify the dynamics of the Coulomb interaction and the results we have obtained here. Disorder in MATBG comes in many avatars and includes long-ranged impurities which are poorly screened at low energy, as well as various forms of correlated disorder and perhaps, the most dominant source being the unconventional “twist-angle” disorder [69], [70].

Within our weak-coupling approach, we have focus on s-wave superconductivity with a uniform phase over the entire Fermi surface. While disorder might prompt such a state, recent experiments suggest that the superconducting state near \(\nu = -2\)​ breaks the rotational symmetry [14]. This implies that there is a strong indication that the SC order parameter is not purely s-wave and, moreover, belongs to a multi-component representation. Thus, it is important to understand under what conditions such a state is preferred and whether the mixture of Coulomb repulsion and umklapp phonons may favor such a state.

We end by noting that in two spatial dimensions, the true SC transition is described by the Berezinskii-Kosterlitz-Thouless (BKT) transition, with \(T_c=\pi D_s^-/2\)​, where \(D_s^-\)​ is the superfluid stiffness at \(T\rightarrow T_c^-\)​. In this paper, we identify \(T_c\)​ with the scale associated with the formation of Cooper pairs, instead of the phase-ordering scale, which serves as an upper bound on the transition temperature. The geometric properties associated with the non-exponentially localizable Wannier functions for nearly flat bands [72]–[76] leads to an additional contribution to \(D_s\)​. We leave a detailed study of such effects within our present framework for future work.

C.L. acknowledges support from the STC Center for Integrated Quantum Materials, NSF Grant No. DMR-1231319, and from the Gordon and Betty Moore Foundation through Grant GBMF8682. D.C. is supported by a faculty startup grant at Cornell University. J.R. acknowledges funding by the Israeli Science Foundation under grant number 994/19.

6 Continuum model

In this appendix, we provide additional details of the bandstructure model used in the calculations presented in the main text. We use the continuum model introduced in Ref.[45] and restated here for completeness. The Hamiltonian for a valley \(\xi=-1,1\)​ and spin \(\sigma=\downarrow,\uparrow=-1,1\)​ takes the form

\[\begin{equation}H^{(\xi,\sigma)}=\begin{pmatrix} H_1 & U^\dagger \\ U & H_2 \end{pmatrix}\label{eq:ham-95cont}\end{equation}\] in the basis of \((A_1, B_1, A_2, B_2)\)​ sites of the original two layers (\(l=1,2\)​). The matrices \(H_{l}\)​ correspond to the intralayer Hamiltonian of the layer \(l\)​ and are explicitly given as

\[\begin{equation}H_{1}= -\frac{\hbar v}{a}\begin{pmatrix} 0 &e^{-i\xi \theta/2} k_{-} a + \frac{4 \pi}{3} \\ e^{i\xi \theta/2} k_{+} a + \frac{4 \pi}{3} & 0 \end{pmatrix}\,,\quad H_{2}= -\frac{\hbar v}{a}\begin{pmatrix} 0 &e^{i\xi \theta/2} k_{-} a + \frac{4 \pi}{3} \\ e^{-i\xi \theta/2} k_{+} a + \frac{4 \pi}{3} & 0 \end{pmatrix}\end{equation}\] where \(k_{\pm}=k_x \pm i k_y\)​ and \(k_x\)​,\(k_y\)​ are crystal lattice momenta measured with respect to the original \(\Gamma\)​ points of the graphene layers. The \(4\pi/3\)​ terms in the matrix are remnants of the low-energy expansion of the graphene monolayer Hamiltonians around the \(\boldsymbol{K}\)​ and \(\boldsymbol{K}'\)​ points of the original layers. The MBZ of MATBG is defined as in the inset of Fig.1a with the two reciprocal lattice vectors being

\[\begin{equation}\boldsymbol{G}_1^{M} = -\frac{2\pi}{\sqrt{3} L_M} \begin{pmatrix} 1 \\ \sqrt{3} \end{pmatrix}\,,\quad \boldsymbol{G}_2^{M} = \frac{4\pi}{\sqrt{3} L_M} \begin{pmatrix} 1 \\ 0 \end{pmatrix}\,.\end{equation}\] Here, we use the moiré real space lattice constant, \(L_M = a/2\sin(\theta/2)\)​, and \(\theta=1.05^\circ\)​. The matrix, \(U\)​, is the effective moiré interlayer coupling given by: \[\begin{gather} U = \begin{pmatrix} u & u'\\ u' & u\\\end{pmatrix}+\begin{pmatrix} u & u'\nu^{-\xi}\\ u'\nu^{\xi} & u\\\end{pmatrix}e^{i\xi \boldsymbol{G}_{1}^M\cdot\boldsymbol{r}}+\nonumber &+\begin{pmatrix} u & u'\nu^{\xi}\\ u'\nu^{-\xi} & u\\\end{pmatrix}e^{i\xi \left(\boldsymbol{G}_{1}^M+\boldsymbol{G}_{2}^M\right)\cdot\boldsymbol{r}}\,,\end{gather}\]

where \(\nu=e^{i2\pi/3}\)​. We take the energy scale as \(\hbar v/a = 2.1354\)​ eV and the lattice constant \(a = 0.246\)​ nm. The interlayer coupling terms \(u\)​ and \(u'\)​ are taken as \(u=0.0797\)​ eV and \(u'=0.0975\)​ eV. For a detailed analysis of the origins of the Hamiltonian and discussion of the significance and numerical value of the coefficients, we refer the reader to Ref. [45] and references therein. In practice, the integers \(m_1\)​ and \(m_2\)​ in total cover at most only around \(\sim 60\)​ combinations, which stems from using the cutoff procedure explained in Ref. [45]. We stress that no qualitative change to the bandstructure would be observed if the cutoff were to be increased. This is because majority of the spectral weight is in fact present only in the \(m_1,m_2\in [-3,3]\)​ range covering \(49\)​ possible combinations - that fact lies also behind the saturation of the phonon \(T_C\)​ domes with the inclusion of \(3\)​G umklapp processes, Fig. 3a. The bandstructure for the two flat bands for valley \(\xi=1\)​ is shown in Fig.1b. We refer to the bands Fig. 1b as the flat bands each with a bandwidth \(W \approx 4\)​ meV. Accordingly, all other bands are called non-flat bands and are separated from the flat-bands by a bandgap, \(\Delta_{\rm band} \approx 12\)​ meV.

7 Numerical solution of the (non-)linear gap equations

Here we detail the procedure used to solve the linearized gap equation introduced in the main text. In the second half of this section we extend this calculation to the full non-linear gap problem, which is then employed in the process of analytic continuation, necessary for computation of the tunneling density of states.

As explained in the main text, the linearized gap equation is posed as an eigenvalue problem, Eq. \(\ref{eq:sc-95gap-95equation}\): \[\begin{gather} \Delta(i\omega,{k}) = -T \sum_{\nu} \sum_{{p}} K(i{\omega},{k};i\nu,{p})\Delta(i\nu,{p}) \label{eq:app-95sc-95gap-95equation},\end{gather}\]

where the kernel \(K(...)\)​ is given by \[\begin{gather} \label{eq:app-95Kernel} K(i{\omega},{k};i\nu,{p}) \equiv \frac{1}{(2\pi)^2} \int d{\Omega}_{\boldsymbol{p}} {\cal V}^{\textrm{int}}_{-\xi,\xi}(\boldsymbol{k}-\boldsymbol{p},i\omega-i\nu) \frac{\Lambda( \boldsymbol{p} ,\boldsymbol{k})\Lambda(-\boldsymbol{p},\boldsymbol{k}) }{\nu^2+{E}^2_{\xi,\boldsymbol{p}}}\end{gather}\]

and where \(\int d{\Omega}_{\boldsymbol{p}}\)​ denotes integration over the angle between vectors \(\boldsymbol{k}\)​ and \(\boldsymbol{p}\)​ for a fixed direction of \(\boldsymbol{k}\)​. To proceed with the analysis, we first make a simplifying assumption that the superconducting order parameter is a function of Matsubara frequency, \(i\omega\)​, and magnitude of the momentum, \(k\)​, but not on the angle, respectively. This allows us to solve the resulting integral equation on a two-dimensional grid of \(k\times\omega\)​ points. With this assumption, the eigenvalue problem simplifies then to \[\begin{gather} \Delta(i\omega,k) = -T \sum_{\nu} \sum_{k} K(i{\omega},k;i\nu,p)\Delta(i\nu,p) \label{eq:app-95sc-95gap-95equation-95mom-95mag}\end{gather}\]

with the momentum direction-averaged kernel, described in detail below, is given by \[\begin{gather} K(i{\omega},k;i\nu,p) \equiv {1\over (2\pi)^2} \int d{\Omega}_{\boldsymbol{p}} {\cal V}^{\textrm{int}}_{-\xi,\xi}(|\boldsymbol{k}-\boldsymbol{p}|,i\omega-i\nu) \frac{\Lambda( \boldsymbol{p} ,\boldsymbol{k})\Lambda(-\boldsymbol{p},\boldsymbol{k}) }{\nu^2+{E}^2_{\xi,\boldsymbol{p}}} \label{eq:app-95Kernel-95mom-95mag}\end{gather}\] In the above kernel we have made the assumption that the interaction \(V^{\textrm{int}}(...)\)​ has no explicit dependence on the angle associated with the momentum, and only depends on its magnitude. This simplification allows us to compute the dynamical dielectric function along just one momentum direction and then estimate it along the other directions by simply comparing magnitudes of \(|\boldsymbol{k}-\boldsymbol{p}|\)​.

We start the calculation by pre-computing a 2D MBZ mesh of points, with their associated Bloch wavefunctions and energies. In the calculations we use a mesh of \(\sim{8000}\)​ MBZ points. To carry out an angle-average of the kernel from Eq. \(\ref{eq:app-95Kernel-95mom-95mag}\), we first fix a particular direction of vector \(\boldsymbol{k}\)​ (upon verifying that the conclusions are not dependent on the specific direction) and then identify all \(\boldsymbol{p}\)​ points that are of magnitude \(p\)​ (within a resolution admitted by the mesh). We then estimate the angular average Eq. \(\ref{eq:app-95Kernel-95mom-95mag}\) by averaging over these \(\boldsymbol{p}\)​ points — in practice, \(\sim 50-100\)​ points are used for each pair of \(k\)​ and \(p\)​ momentum values.

To determine the critical temperature, we seek the temperature \(T\)​ for which Eq. \(\ref{eq:app-95sc-95gap-95equation-95mom-95mag}\) has an eigenvalue of unity. In practice, we carry out a bisection method search for a \(T\)​ giving an eigenvalue within \(\pm 0.001\)​ of the unity. In the calculations we use a linearly spaced grid of \(30\)\(k\)​ points ranging from the center of the MBZ \(\bar{\Gamma}\)​ to the \(\bar{K}\)​ point. For the Matsubara grid we employed both an exact Matsubara frequency summation as well as an approximate scheme, upon verifying that the approximate Matsubara grid agrees (in determining the value of the critical temperature) with the exact summation to 3 significant figures. The approximate Matsubara grid was chosen to consist of 10 first Matsubara frequencies followed by 20 linearly spaced frequencies starting from the 11th Matsubara frequency to the UV cutoff. Several UV cutoffs were tested: multiples of plasma frequency, multiples of the Debye frequency, or, multiples of the system bandwidth. All cutoffs were found consistent with each other provided that the UV cutoff exceeds, roughly, \(30\)​ meV for the bandstructure used.

To determine a self-consistent gap, we extend the calculation by modifying the kernel in the integral equation Eq. \(\ref{eq:app-95sc-95gap-95equation-95mom-95mag}\) as \[\begin{gather} K(i{\omega},k;i\nu,p)\to K_{\textrm{SC}}(i{\omega},k;i\nu,p) \equiv {1\over (2\pi)^2} \int d{\Omega}_{\boldsymbol{p}} {\cal V}^{\textrm{int}}_{-\xi,\xi}(|\boldsymbol{k}-\boldsymbol{p}|,i\omega-i\nu) \frac{\Lambda( \boldsymbol{p} ,\boldsymbol{k})\Lambda(-\boldsymbol{p},\boldsymbol{k}) }{\nu^2+{E}^2_{\xi,\boldsymbol{p}}+\Delta^2(i\nu,\boldsymbol{p})} \label{eq:app-95Kernel-95mom-95mag-95SC}\end{gather}\] where \(\Delta(i\nu,\boldsymbol{p})\)​ is now the self-consistent gap. Starting with the eigenvector of the linearized gap equation as an input, we then self-consistently solve the gap equation. For temperatures, \(T \ll T_C\)​, as used in Fig. 6c, the gap converges within \(10-20\)​ iterations to below \(0.1\%\)​ total relative error (difference between successive self-consistent steps). The resulting gap function follows the BCS result; see Fig.7(a).

Figure 7: (a) Superconducting gap as a function of temperature obtained from the self-consistent calculation. Plotted quantity (blue) corresponds to an average of the superconducting gap \(\Delta(i\nu,\boldsymbol{p})\) over the Fermi surface. Orange line corresponds to the BCS interpolating solution (accurate withing few percent). (b) Analytically continued superconducting gap \(\Delta(\omega)\) as a function of real frequency obtained from the self-consistent calculation. As in (a) the gap was averaged over the Fermi surface. In (a),(b) we consider phonon mediated superconducting with only 0G processes at a characteristic filling of \(\nu = 1\). Note how \({\rm Im}[\Delta(\omega)]\) is zero until gap opens, followed by a peak at the characteristic frequency for 0G processes, as expected from Fig. 3b.

8 Analytic continuation and single particle density of states

To determine the tunneling density of states it is necessary to obtain the self-consistent gap as a function of real frequency [65], [77]. To avoid this complex procedure, few schemes of carrying out an analytic continuation from the Matsubara frequency to real frequencies have been introduced in literature with the two most common relying on either Padé approximants [78], [79], or, an iterative solution to the Eliashberg equations [80]. The latter approach is not directly applicable in the case of MATBG as it relies on an assumption of infinite system bandwidth. We employ therefore a Padé approximation scheme, in particular, the approach detailed in Ref. [81].

As an input, we start with the self-consistent gap obtained following the procedure detailed in the previous section. This gap is then approximated with the help of Padé approximant and analytically continued to real frequencies. An example of real frequency dependence of a gap is shown in Fig. 7(b). This self-consistent gap defined at real-frequencies is then used to compute the spectral function [57],

\[\begin{equation}\mathcal{A}_{\{\gamma\}}(\omega,\boldsymbol{k}) = -\frac{1}{\pi}{\rm Im}\left[\frac{(\omega+i0)+E_{\boldsymbol{k},\{\gamma\}}}{(\omega+i0)^2-E^2_{\boldsymbol{k},\{\gamma\}}-\Delta^2(\omega)}\right]\,,\end{equation}\]

and the tunneling density of states \(N_{\textrm{SC}}(\omega-\mu)\)​,

\[\begin{equation}N_{SC}(\omega-\mu) = \sum_{\{\gamma\}} \int_{MBZ} d^2 \boldsymbol{k}~ \mathcal{A}_{\{\gamma\}}(\omega,\boldsymbol{k}),\end{equation}\]

shown in Fig. 6c.


[1] Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,”

[2] Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L. Tomarken, Jason Y. Luo, Javier D. Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, Ray C. Ashoori, and Pablo Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,”

[3] Matthew Yankowitz, Shaowen Chen, Hryhoriy Polshyn, Yuxuan Zhang, K. Watanabe, T. Taniguchi, David Graf, Andrea F. Young, and Cory R. Dean, “Tuning superconductivity in twisted bilayer graphene,”

[4] Xiaobo Lu, Petr Stepanov, Wei Yang, Ming Xie, Mohammed Ali Aamir, Ipsita Das, Carles Urgell, Kenji Watanabe, Takashi Taniguchi, Guangyu Zhang, Adrian Bachtold, Allan H. MacDonald, and Dmitri K. Efetov, “Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene,” 10.1038/s41586-019-1695-0.

[5] Leon Balents, Cory R. Dean, Dmitri K. Efetov, and Andrea F. Young, “Superconductivity and strong correlations in moiréflat bands,” 10.1038/s41567-020-0906-9.

[6] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, “Graphene bilayer with a twist: Electronic structure,”

[7] E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, “Flat bands in slightly twisted bilayer graphene: Tight-binding calculations,” 10.1103/PhysRevB.82.121407.

[8] Rafi Bistritzer and Allan H. MacDonald, “Moiré bands in twisted double-layer graphene,” 10.1073/pnas.1108174108.

[9] Alexander Kerelsky, Leo J. McGilly, Dante M. Kennes, Lede Xian, Matthew Yankowitz, Shaowen Chen, K. Watanabe, T. Taniguchi, James Hone, Cory Dean, Angel Rubio, and Abhay N. Pasupathy, “Maximized electron interactions at the magic angle in twisted bilayer graphene,”

[10] Youngjoon Choi, Jeannette Kemmer, Yang Peng, Alex Thomson, Harpreet Arora, Robert Polski, Yiran Zhang, Hechen Ren, Jason Alicea, Gil Refael, Felix von Oppen, Kenji Watanabe, Takashi Taniguchi, and Stevan Nadj-Perge, “Electronic correlations in twisted bilayer graphene near the magic angle,”

[11] Yuhang Jiang, Xinyuan Lai, Kenji Watanabe, Takashi Taniguchi, Kristjan Haule, Jinhai Mao, and Eva Y. Andrei, “Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene,” 10.1038/s41586-019-1460-4.

[12] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, Ady Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, “Cascade of phase transitions and dirac revivals in magic-angle graphene,”

[13] Dillon Wong, Kevin P. Nuckolls, Myungchul Oh, Biao Lian, Yonglong Xie, Sangjun Jeon, Kenji Watanabe, Takashi Taniguchi, B. Andrei Bernevig, and Ali Yazdani, “Cascade of electronic transitions in magic-angle twisted bilayer graphene,”

[14] Yuan Cao, Daniel Rodan-Legrain, Jeong Min Park, Fanqi Noah Yuan, Kenji Watanabe, Takashi Taniguchi, Rafael M. Fernandes, Liang Fu, and Pablo Jarillo-Herrero, “Nematicity and Competing Orders in Superconducting Magic-Angle Graphene,” arXiv e-prints , arXiv:2004.04148 (2020),

[15] Y. Cao, J. Y. Luo, V. Fatemi, S. Fang, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, “Superlattice-induced insulating states and valley-protected orbits in twisted bilayer graphene,” 10.1103/PhysRevLett.117.116804.

[16] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, “From quantum matter to high-temperature superconductivity in copper oxides,”

[17] Hoi Chun Po, Liujun Zou, Ashvin Vishwanath, and T. Senthil, “Origin of mott insulating behavior and superconductivity in twisted bilayer graphene,”

[18] Liujun Zou, Hoi Chun Po, Ashvin Vishwanath, and T. Senthil, “Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions,”

[19] Hoi Chun Po, Liujun Zou, T. Senthil, and Ashvin Vishwanath, “Faithful tight-binding models and fragile topology of magic-angle bilayer graphene,”

[20] Zhida Song, Zhijun Wang, Wujun Shi, Gang Li, Chen Fang, and B. Andrei Bernevig, “All magic angles in twisted bilayer graphene are topological,”

[21] Junyeong Ahn, Sungjoon Park, and Bohm-Jung Yang, “Failure of nielsen-ninomiya theorem and fragile topology in two-dimensional systems with space-time inversion symmetry: Application to twisted bilayer graphene at magic angle,”

[22] Mikito Koshino and Young-Woo Son, “Moiré phonons in twisted bilayer graphene,”

[23] Héctor Ochoa, “Moiré-pattern fluctuations and electron-phason coupling in twisted bilayer graphene,”

[24] Cyprian Lewandowski and Leonid Levitov, “Intrinsically undamped plasmon modes in narrow electron bands,” 10.1073/pnas.1909069116.

[25] Petr Stepanov, Ipsita Das, Xiaobo Lu, Ali Fahimniya, Kenji Watanabe, Takashi Taniguchi, Frank H. L. Koppens, Johannes Lischner, Leonid Levitov, and Dmitri K. Efetov, “The interplay of insulating and superconducting orders in magic-angle graphene bilayers,” (2019),

[26] Yu Saito, Jingyuan Ge, Kenji Watanabe, Takashi Taniguchi, and Andrea F. Young, “Independent superconductors and correlated insulators in twisted bilayer graphene,”

[27] T. Senthil, “What drives superconductivity in twisted bilayer graphene?” 10.36471/jccm_may_2020_03.

[28] Xiaoxue Liu, Zhi Wang, K. Watanabe, T. Taniguchi, Oskar Vafek, and J. I. A. Li, “Tuning electron correlation in magic-angle twisted bilayer graphene using coulomb screening,” (2020),

[29] Risto Ojajärvi, Timo Hyart, Mihail A. Silaev, and Tero T. Heikkilä, “Competition of electron-phonon mediated superconductivity and stoner magnetism on a flat band,”

[30] Hiroki Isobe, Noah F. Q. Yuan, and Liang Fu, “Unconventional superconductivity and density waves in twisted bilayer graphene,”

[31] Teemu J. Peltonen, Risto Ojajärvi, and Tero T. Heikkilä, “Mean-field theory for superconductivity in twisted bilayer graphene,”

[32] Fengcheng Wu, A. H. MacDonald, and Ivar Martin, “Theory of phonon-mediated superconductivity in twisted bilayer graphene,”

[33] Dante M. Kennes, Johannes Lischner, and Christoph Karrasch, “Strong correlations and \(d+\mathit{id}\) superconductivity in twisted bilayer graphene,”

[34] Young Woo Choi and Hyoung Joon Choi, “Strong electron-phonon coupling, electron-hole asymmetry, and nonadiabaticity in magic-angle twisted bilayer graphene,” 10.1103/PhysRevB.98.241412.

[35] Yi-Zhuang You and Ashvin Vishwanath, “Superconductivity from valley fluctuations and approximate so(4) symmetry in a weak coupling theory of twisted bilayer graphene,” 10.1038/s41535-019-0153-4.

[36] Biao Lian, Zhijun Wang, and B. Andrei Bernevig, “Twisted bilayer graphene: A phonon-driven superconductor,” 10.1103/PhysRevLett.122.257002.

[37] Yu-Ping Lin and Rahul M. Nandkishore, “Chiral twist on the high-\({T}_{c}\) phase diagram in moiré heterostructures,”

[38] Fabian Schrodi, Alex Aperis, and Peter M. Oppeneer, “Prominent cooper pairing away from the fermi level and its spectroscopic signature in twisted bilayer graphene,” 10.1103/PhysRevResearch.2.012066.

[39] Girish Sharma, Maxim Trushin, Oleg P. Sushkov, Giovanni Vignale, and Shaffique Adam, “Superconductivity from collective excitations in magic-angle twisted bilayer graphene,” 10.1103/PhysRevResearch.2.022040.

[40] Eslam Khalaf, Shubhayu Chatterjee, Nick Bultinck, Michael P. Zaletel, and Ashvin Vishwanath, “Charged skyrmions and topological origin of superconductivity in magic angle graphene,” (2020),

[41] Maine Christos, Subir Sachdev, and Mathias Scheurer, “Superconductivity, correlated insulators, and wess-zumino-witten terms in twisted bilayer graphene,” (2020),

[42] Douglas J Scalapino, “The electron-phonon interaction and strong-coupling superconductors,” in Superconductivity, edited by R.D. Parks (Marcel Dekker, Inc., New York, 1969) pp. 449–560.

[43] Yang-Zhi Chou, Yu-Ping Lin, Sankar Das Sarma, and Rahul M. Nandkishore, “Superconductor versus insulator in twisted bilayer graphene,”

[44] Fengcheng Wu, Euyheon Hwang, and Sankar Das Sarma, “Phonon-induced giant linear-in-\(t\) resistivity in magic angle twisted bilayer graphene: Ordinary strangeness and exotic superconductivity,”

[45] Mikito Koshino, Noah F. Q. Yuan, Takashi Koretsune, Masayuki Ochi, Kazuhiko Kuroki, and Liang Fu, “Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene,”

[46] Dmitri K. Efetov and Philip Kim, “Controlling electron-phonon interactions in graphene at ultrahigh carrier densities,”

[47] Jian-Hao Chen, Chaun Jang, Shudong Xiao, Masa Ishigami, and Michael S. Fuhrer, “Intrinsic and extrinsic performance limits of graphene devices on sio2,” 10.1038/nnano.2008.58.

[48] Francisco Guinea and Niels R. Walet, “Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers,”

[49] Zachary A. H. Goodwin, Valerio Vitale, Xia Liang, Arash A. Mostofi, and Johannes Lischner, “Hartree theory calculations of quasiparticle properties in twisted bilayer graphene,” arXiv e-prints , arXiv:2004.14784 (2020),

[50] Yasutami Takada, “Plasmon mechanism of superconductivity in two- and three-dimensional electron systems,”,

[51] Dan Phan and Andrey V. Chubukov, “Kohn-luttinger correction to \({T}_{c}\) in a phonon superconductor,” 10.1103/PhysRevB.101.024503.

[52] Jian Kang and Oskar Vafek, “Symmetry, maximally localized wannier states, and a low-energy model for twisted bilayer graphene narrow bands,”

[53] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of superconductivity,”

[54] W. L. McMillan, “Transition temperature of strong-coupled superconductors,” 10.1103/PhysRev.167.331.

[55] P. B. Allen and R. C. Dynes, “Transition temperature of strong-coupled superconductors reanalyzed,” 10.1103/PhysRevB.12.905.

[56] Jonathan Ruhman and Patrick A. Lee, “Pairing from dynamically screened coulomb repulsion in bismuth,”

[57] Piers Coleman, (Cambridge University Press, 2015).

[58] Frank Stern, “Polarizability of a two-dimensional electron gas,” 10.1103/PhysRevLett.18.546.

[59] J. M. Pizarro, M. Rösner, R. Thomale, R. Valentı́, and T. O. Wehling, “Internal screening and dielectric engineering in magic-angle twisted bilayer graphene,” 10.1103/PhysRevB.100.161102.

[60] Zachary A. H. Goodwin, Fabiano Corsetti, Arash A. Mostofi, and Johannes Lischner, “Attractive electron-electron interactions from internal screening in magic-angle twisted bilayer graphene,” 10.1103/PhysRevB.100.235424.

[61] E. H. Hwang and S. Das Sarma, “Dielectric function, screening, and plasmons in two-dimensional graphene,”

[62] B Wunsch, T Stauber, F Sols, and F Guinea, “Dynamical polarization of graphene at finite doping,”

[63] Jonathan Ruhman and Patrick A. Lee, “Pairing from dynamically screened coulomb repulsion in bismuth,”

[64] J. R. Schrieffer, D. J. Scalapino, and J. W. Wilkins, “Effective tunneling density of states in superconductors,” 10.1103/PhysRevLett.10.336.


[66] Vladyslav Kozii, Michael P. Zaletel, and Nick Bultinck, “Superconductivity in a doped valley coherent insulator in magic angle graphene: Goldstone-mediated pairing and Kohn-Luttinger mechanism,” arXiv e-prints , arXiv:2005.12961 (2020),

[67] Yuan Cao, Debanjan Chowdhury, Daniel Rodan-Legrain, Oriol Rubies-Bigorda, Kenji Watanabe, Takashi Taniguchi, T. Senthil, and Pablo Jarillo-Herrero, “Strange metal in magic-angle graphene with near planckian dissipation,”

[68] Hryhoriy Polshyn, Matthew Yankowitz, Shaowen Chen, Yuxuan Zhang, K. Watanabe, T. Taniguchi, Cory R. Dean, and Andrea F. Young, “Large linear-in-temperature resistivity in twisted bilayer graphene,”

[69] A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani, D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and E. Zeldov, “Mapping the twist-angle disorder and landau levels in magic-angle graphene,”

[70] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, Ady Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, “Cascade of phase transitions and dirac revivals in magic-angle graphene,”

[71] Hyobin Yoo, Rebecca Engelke, Stephen Carr, Shiang Fang, Kuan Zhang, Paul Cazeaux, Suk Hyun Sung, Robert Hovden, Adam W. Tsen, Takashi Taniguchi, Kenji Watanabe, Gyu-Chul Yi, Miyoung Kim, Mitchell Luskin, Ellad B. Tadmor, Efthimios Kaxiras, and Philip Kim, “Atomic and electronic reconstruction at the van der waals interface in twisted bilayer graphene,” 10.1038/s41563-019-0346-z.

[72] Sebastiano Peotta and Päivi Törmä, “Superfluidity in topologically nontrivial flat bands,”

[73] Fang Xie, Zhida Song, Biao Lian, and B. Andrei Bernevig, “Topology-bounded superfluid weight in twisted bilayer graphene,” 10.1103/PhysRevLett.124.167002.

[74] A. Julku, T. J. Peltonen, L. Liang, T. T. Heikkilä, and P. Törmä, “Superfluid weight and berezinskii-kosterlitz-thouless transition temperature of twisted bilayer graphene,”

[75] Xiang Hu, Timo Hyart, Dmitry I. Pikulin, and Enrico Rossi, “Geometric and conventional contribution to the superfluid weight in twisted bilayer graphene,” 10.1103/PhysRevLett.123.237002.

[76] Johannes S. Hofmann, Erez Berg, and Debanjan Chowdhury, “Superconductivity, pseudogap, and phase separation in topological flat bands: a quantum Monte Carlo study,” arXiv e-prints , arXiv:1912.08848 (2019),

[77] E. R. Margine and F. Giustino, “Anisotropic migdal-eliashberg theory using wannier functions,” 10.1103/PhysRevB.87.024505.

[78] H. J. Vidberg and J. W. Serene, “Solving the eliashberg equations by means ofn-point padéapproximants,”

[79] C.R. Leavens and D.S. Ritchie, “Extension of the n-point padé approximants solution of the eliashberg equations to t \(~\) tc,”

[80] F. Marsiglio, M. Schossmann, and J. P. Carbotte, “Iterative analytic continuation of the electron self-energy to the real axis,”

[81] K. S. D. Beach, R. J. Gooding, and F. Marsiglio, “Reliable padé analytical continuation method based on a high-accuracy symbolic computation algorithm,”

  1. The devices are nominally at different twist angles and have varying levels of disorder; see [27] for a further discussion.↩︎

  2. A detailed study of these umklapp processes and on their role in superconductivity has not been clarified in earlier studies [32], [36], [44].↩︎

  3. \(N = 1\)​ corresponds to the physical limit of MATBG.↩︎

  4. We note that the eigenstates and eigenvalues of the two valleys \(\xi=\pm\)​ are time-reversed partners and accompanied by a simultaneous complex conjugation and \(\boldsymbol{k}\to-\boldsymbol{k}\)​ transformation.↩︎

  5. Here we use the notation that \(\rho_{\xi; i}(\boldsymbol{q})\equiv\rho_{\xi; i}(\boldsymbol{q},i\omega=0)\)​ as per the definition in Eq. \(\ref{eq:ham-95int}\)↩︎

  6. Ref. [22] points out that certain in-plane phonon modes can develop small “gaps” as a result of the moiré potential, which we ignore in our study.↩︎

  7. Note that here we have neglected the extended \(s\)​-wave of the \(D_3\)​ symmetry group.↩︎

  8. On the other hand, the precise values of the transition temperature, \(T_c\)​, are not meant to be taken literally.↩︎

  9. We stress that if the resulting \(T_C\)​ for a given filling becomes of order of the chemical potential, \(\mu\)​, then it is necessary to determine the chemical potential self-consistently [51]. In our results, however, this limit is never reached, but the theory at very low fillings \(|\nu| \lesssim 0.2\)​ for higher-order (\(m \geq 2\)​) umklapp processes starts to approach, \(T_C\approx 0.5 \mu\)​.↩︎

  10. The exact location of the van-Hove singularity at \(\nu\approx 0.63\)​ is a consequence of the underlying details of the bandstructure in Ref. [45].↩︎

  11. The bare coupling constants for the Coulomb and electron-phonon interactions are chosen to be close to the experimentally reported values relevant for MATBG; we choose \(N=20\)​.↩︎

  12. When appearing inside the band we expect the resonance to be possibly weaker as it will get suppressed by a finite electronic density of states.↩︎