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] ¹, 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, or by stabilizing superconductivity at angles far away from the magic angle where electronic correlations are suppressed[29]. 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 [30]–[38], whilst other works extend the treatment to an Eliashberg framework [39]–[41]. In particular, Ref.[41] applied Eliashberg theory to study pairing mediated by a single Einstein Einstein phonon, where the electron-phonon coupling strength was chosen phenomenologically. Ref.[40] considered plasmon-mediated pairing in a hexagonal lattice model for MATBG. See also two recent studies [42], [43] focusing on the interplay of insulating and superconducting states, 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 [44]. 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 large\(-N\) expansion, such that our results correspond to the weak-coupling limit (obtained by taking \(N\gg 1\)). We argue that in spite of this approximation, the results shed 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 [45]. 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 ² Below, we show that the SC transition temperature, \(T_c\), due to acoustic phonons of graphene is enhanced significantly upon including the effects of these 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.

To verify our predictions, we propose a tunneling experiment. As is well known when the coupling is strong enough “dip-hump" features in the tunneling spectra appear at energies that can be associated with the frequency of the bosonic modes that contribute most to pairing [47]. In the present case the umklapp phonons have a very distinct frequency scale associated with the speed of sound in graphene and an integer multiple of the MBZ size. Thus, this allows us to make a rather sharp prediction for the energy at which such features will appear. Note that this prediction does not depend on the model for the electronic band structure.

Finally, we also study the effect of a screening layer, coupled to MATBG by Coulomb interactions, on the density dependence of \(T_c\). The distinct effect of screening on the plasmonic modes in graphene is discussed and proposed as a method to distinguish their contribution to pairing.

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 - 9.

2 Model↩︎

The non-interacting part of our model, \(\mathcal{H}_0\), leads to the action, \[\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},\] 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 \(1/N\) expansion ³. 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. [48] for MATBG at a twist angle \(\theta=1.05^\circ\) ⁴. 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{align} \label{eq:ham95int} &\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{align}\tag{1}\] 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'}\)): \[\label{eq:Int} {\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})\,.\tag{2}\] 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 as well as make the Coulomb and phonon interaction manifestly intervalley only through the \(\delta_{\xi,-\xi'}\) (see below). 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. 1 .

The first term in Eq. 2 is the dynamically screened Coulomb interaction, which is given by \[\label{eq:V95C} {\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}\delta_{\xi,-\xi'}\,,\tag{3}\] 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, \[\label{eq:ham95phonon} {\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}}),\tag{4}\] where \(a_{\boldsymbol{q}},~a^\dagger_{\boldsymbol{q}}\) represent the phonon annihilation and creation operators, respectively⁵. 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 [49], [50]. After integrating out the phonons, we obtain the density-density interaction term of the form Eq. 1 where, \[\label{eq:V95ph} {\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)} \delta_{\xi,-\xi'},\tag{5}\] 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. 1 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. In what follows we will also relate our findings to previous works that focused on pairing mechanisms mediated by optical phonons[32], [34], [36], [41], [46] irrespective of their microscopic origins or frequency, i.e. whether they are mini-optical phonons arising from folding of the acoustic dispersion, actual optical phonons of graphene or just phenomenological Einstein modes. We also clarify that the phonons we consider correspond to symmetric motion of atoms in both graphene layers, but as Ref. [22] points out, there are also in-plane phonon modes which we ignore in our study corresponding to asymmetric motion of atoms in both layers that can develop small “gaps” as a result of the moiré potential.

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 (for details on how to include them we refer the reader to Ref. [51]). The gap function in the spin-singlet, \(s\)-wave pairing channel with zero center of mass momentum is defined as, \[\begin{align} \Delta(i\omega,{k}) \equiv {\langle}{\varepsilon}_{{\sigma}{\sigma}'} c_{\omega,\boldsymbol{k}\{{\gamma}\}} c_{-\omega,-\boldsymbol{k}\{{\gamma}'\}}{\rangle}, \end{align}\] 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}\). The latter assumption is consistent with our assumption of intervalley Cooper pairing from valleys that are time-reversed partners, as explained before. Note that the above requires an intervalley scattering, Fig.1a, which imposes the condition \(\xi \neq \xi'\) in the interaction Eq. 1 . The gap equation thereby reduces to an eigenvalue problem \[\begin{align} \Delta(i\omega,{k}) = -T \sum_{\nu} \sum_{{p}} K(i{\omega},{k};i\nu,{p})\Delta(i\nu,{p}) \label{eq:sc95gap95equation}, \end{align}\tag{6}\] where the kernel \(K(...)\) is given by \[\begin{align} \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{align}\tag{7}\] 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}\). ⁶ In the above expression, as we are focusing on an intervalley pairing, 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. Accordingly the valley indices are omitted henceforth.

The Eliashberg equation for the pairing gap in Eq. 6 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 [52], [53]. We have (artificially) introduced a control parameter that selects a subset of the relevant processes and weakens the strength of the pairing mechanism. This places our analysis in the weak-coupling regime allowing us to focus on solving one equation only for the superconducting order parameter. Namely, we have neglected the mass and dispersion renormalizations, \(Z\) and \(\chi\) (as these are subleading in \(1/N\)). For more details on their inclusion we refer the readers to Ref. [51]. We also note that the mass renormalization \(Z\) tends to reduce \(T_c\) by a factor [54], while the \(\chi\) becomes most important when \(T_c\) is of the order of the Fermi energy \({\varepsilon}_F\) [54], [55] (i.e. at strong coupling).

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. On the other hand, it is important to recall that the large-\(N\) model we consider here is somewhat different from the experimental situation. Therefore, the actual value of \(T_c\) takes meaning only by comparison and not by its absolute value.

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. 7 . 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 8.

Before proceeding with our analysis of the results, let us pause to discuss some of the properties of the eigenvector in Eq. 6 . 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. 6 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) [56] as well as an eigenstate solution that changes sign — a necessary requirement to satisfy the gap equation 6 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. 6 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 ⁷.

4.1 Phonon mediated superconductivity↩︎

We start our discussion by focusing exclusively on the effects of the phonon mediated interaction, Eq. 5 . 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. 3 to zero and compute \(T_c\) due to the phonon mediated interaction (\({\cal V}_{\xi,\xi'}^{\textrm{ph}}\)) in Eq. 4 . 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{align} {\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}},\label{eq:spectral95function} \end{align}\tag{8}\] where \(\boldsymbol{q}=\boldsymbol{\tilde{q}}+\boldsymbol{G}\), as defined previously (see paragraph following Eq. 5 ). 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\).

The enhancement of \(T_C\) with umklapp processes can be anticipated by looking at the form of Eq. 5 . ⁸ The effective electron-electron interaction at the lowest frequency, \(\omega=0\), is a constant independent of \(q\). Momentum dependence enters in Eq.@eq:eq:V95ph only for finite frequency. For the first such frequency (\(n=1\), \(\omega= 2\pi T\)), the interaction has a reversed lorentzian shape in \(q\)-space where it is zero for \(q\to0\) and saturates to the (\(\omega=0\)) value as momentum becomes large. With the above observations, one can predict that as \(q\) increases, the coupling strength and therefore the pairing tendency will be enhanced until some value where the effect of taking higher umklapp contributions into account will lead to a saturation. We highlight, however, that this behavior alone is not behind the sharp saturation seen in Fig.3a, but rather occurs in combination with the presence of form-factors that rapidly vanish past some critical \(mG\) (here \(3\)G).

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 \[\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}\] 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]). It is worthwhile to comment that similar multi-peak features in el-ph spectral functions were seen also in previous ab-initio calculations[34], however the precise connection to umklapp processes was not emphasized.

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 [48], [57]).

Let us now briefly comment on the overall similarity of the dome shape to that of other phonon-pairing based approaches. This is best seen by a comparison with, for example Ref. [41], where an Eliashberg treatment together with a pairing mechanism relying on an optical-like phonon is considered. This similarity simply stems from the fact that in a phonon-only mechanism that produces an attractive interaction over all frequency range, the overall shape of the dome is largely dictated by the electron bandstructure rather than the details of the pairing. As such, any analysis based on a bandstructure model, either tight-binding/continuum/ab-initio in nature, will produce similar-looking domes if the underlying bandstructures are qualitatively similar.

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 [58], \[T_C \approx 1.14 \langle\omega\rangle ~e^{-1/\tilde{\lambda}}\label{eq:BCS95formula}\tag{9}\] 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. 8 to recast the effective coupling constant as [59] \[\tilde{\lambda} = 2 \int d \omega ~\frac{{\alpha}^2 F({\omega})}{\omega},\] and similarly denote the characteristic pairing frequency \({\langle}\omega{\rangle}\) in terms of the moment of the same distribution as \[{\langle}\omega{\rangle}= \frac{2}{\tilde{\lambda}} \int d\omega ~{\alpha}^2 F({\omega})\,.\] 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. 9 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. ⁹

The simple BCS formula, where the characteristic pairing energy is naively given by \({\langle}\omega{\rangle}\), was used here only for a rough estimate and comparison with the results of the Eliashberg theory. It should be mentioned that in this case, where \(\langle {\omega}\rangle \gtrsim {\varepsilon}_F\), the frequency cutoff is sometimes substituted by \({\varepsilon}_F\) [60]–[62]. However, given that this substitution is also non-rigorous [54], we do not perform it here.

Finally, we note that 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 consider the experimental fingerprints of these modes in the tunneling density of states.

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] (see also Appendix 7) and reinstating \({\cal V}^{\textrm{C}}_{\xi,\xi'}\) to the pairing kernel in Eq. 7 .

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 [56], [63]. 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 [64], 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. 3 , as shown in Fig. 5a at a fixed filling (\(\nu\approx1\)). The behavior differs from the one in a conventional 2D Fermi gas [65]. 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. [66], [67]. 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 [68], [69]. 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. 3 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.

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.@eq: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. 7 . 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. 3 is replaced by \[{\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) \delta_{\xi,-\xi'}\,.\] 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] and App. 7 for details). We therefore choose a constant \(\omega_{\textrm{pl}}\) in our calculation. As is evident from Fig. 5c, 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. 5d. 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. 5c. 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. 5d, 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. 5b, the following quantity \[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:coherence95factor95overlap}\tag{10}\] 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\). According to our analysis the strong non-monotonic dependence of \(T_C\) on filling, when resulting from a plasmon mechanism, is a consequence of the form factors. This should be compared with the earlier work relying on a plasmon mechanism [40], where these form-factors are not included.

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 up to \(2\)G processes, the umklapp-driven phonon-mediated mechanism clearly dominates over the plasmon mechanism. ¹⁰ 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. 3 , 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: \[\frac{2\pi e^2}{q}\to \frac{2\pi e^2}{q}\left(1-e^{-2q d}\right)\,,\label{eq:gate95screening}\tag{11}\] where \(d\) is the distance between MATBG and the gate. As a result, the dynamically screened interaction (Eq. 3 ) now changes to \[\begin{align} {\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)\delta_{\xi,-\xi'},\\ \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{align}\] 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.

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. 11 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.@eq:eq:gate95screening , 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. 6b. 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 [59].

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, and in fact in majority of MATBG experiments to date, 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. As such therefore if the thickness of the substrate is increased, our analysis predicts an increase in critical temperature if plasmons participate strongly in mediating the pairing. 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 an 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 [47] (see Fig. 3b). Consider the normalized differential conductance, defined as \[\begin{align} \tilde{G}(\omega-\mu) = \frac{N_{\textrm{SC}}(\omega-\mu)}{N_0(\omega-\mu)},\label{eq:diff95cond} \end{align}\tag{12}\] 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. 6 , 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 8 and 9. 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 “dip-hump” 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. These features are spectral fingerprints of the phonon modes  [47], [70]–[72], which are most visible for modes that are strongly coupled. Therefore, our conclusion that the exchange of umklapp phonons is the strongest contributor to pairing in MATBG leads to the experimental prediction of such "dip-hump" features at these specific frequencies. This frequency of the "dip-hump" features does not rely on the specific choice of electronic band structure model, as it is only related to the sound velocity in graphene and the MBZ size. However, the visibility of the features depends strongly on the electronic density of states background, which is here much enhanced due to the band gap. We caution however that similar-looking resonances, refereed in the literature as replica bands, can also be present in the spectroscopic signal and bear little relation to superconductivity — for a detailed analysis we refer to Refs. [73], [74].

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. As was predicted by Ref. [47] processes where phonons are emitted in the tunneling process lead to these fingerprints, which thus appear at a frequency \({\Delta}(0)+{\omega}_0\), where \({\omega}_0\) is the frequency of the bosonic mode contributing to pairing. 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 “dip-hump” features in the tunneling density of states are most visible when they appear 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 [48] 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 [52], [53] 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], [42], [75]. 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. [75] 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. [76], [77] is an interesting open question.

Another interesting question concerns the unavoidable role of disorder, on the superconducting phenomenology in MATBG [12], [78], [79]. 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 [12], [78].

Within our weak-coupling approach, we have focused 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. Nonetheless, it should also be mentioned that the same umklapp enhancement mechanism for pairing considered here can also play a role in unconventional pairing states and we only considered the effect of these sources of attraction on the s-wave channel for simplicity.

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 [80]–[84] 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.

References↩︎