Cascade of replica bands in flat band systems: predictions for twisted bilayer graphene

We investigate the effect of electron-phonon interactions (EPI) in systems exhibiting one or more flat electron bands close to the Fermi level and a comparatively large phonon energy scale. After solving the self-consistent full-bandwidth Eliashberg equations, we compute Angular Resolved Photo Emission Spectroscopy and Scanning Tunneling Spectroscopy/Microscopy spectra. We obtain a sequence of quasiparticle replica bands in both the normal and superconducting states that originate from frequency dependent features of the electron mass renormalization function. We show that these replica bands can be used to extract the relevant phonon energy scale from experiments. Focusing in particular on twisted bilayer graphene, we predict replica-band formation which, when observed, will shed light on the role of EPI in this archetypal flat-band system.

Introduction. Effects of electron-phonon interactions in metals and superconductors are most accurately modeled by the Migdal-Eliashberg formalism [1,2]. The connection between this theory and Scanning Tunneling Spectroscopy/Microscopy (STS/STM) has been well understood for many decades [3,4]. In more recent years, advances in Angular Resolved Photo Emission Spectroscopy (ARPES) have lead to the possibility of an even richer comparison between theory and experiment [5][6][7][8]. These techniques have been successfully applied to gain better understanding of many materials, such as the hightemperature superconducting cuprates [9,10] and monolayer FeSe on a SrTiO 3 substrate [11][12][13], to name only few examples.
In this work we focus on systems with one or more flat electron energy bands close to the Fermi level, where the term 'flatness' is to be understood in comparison to the phonon energy scale Ω. Calculating ARPES and STS/STM spectra using full-bandwidth Eliashberg theory, we reveal a sequence of quasiparticle replica bands outside the electron bandwidth W of the flat bands, occurring both at positive and negative frequencies. A closer analysis reveals that these spectral features are a direct manifestation of the electron mass renormalization function, and they are located at integer multiples of Ω along the frequency axis. Therefore they can serve as a means to extract the phonon frequency directly from the measured spectra. Although in our formalism we incorporate both the normal and superconducting state, we find that the intensity and location of the replica bands are invariant to the occurrence of superconductivity.
We apply our theory to twisted bilayer graphene (TBG) at the magic angle ∼ 1.1 • [14][15][16], for which we perform calculations of the quasiparticle spectrum using material specific input. This system exhibits two flat bands close to the Fermi energy, as reported by both theory [17][18][19] and experiment [20][21][22][23][24]. These bands are gaped out from the remaining electron energies, and, together with a reasonably large phonon frequency [25,26] that has been shown to be relevant for explaining superconductivity [27], constitute an ideal testing ground for the phenomenon discovered here. Although the predicted intensity of the replica bands is relatively weak, we propose here their detection in TBG in the foreseeable future.
Methodology. We consider a single-branch isotropic Einstein phonon spectrum with characteristic frequency Ω. For simplicity we also treat the electron-phonon scattering strength g 0 as momentum independent. Adopting the Eliashberg formalism in imaginary frequency space, with fermion frequencies ω m = πT (2m + 1) at temperature T , we decompose the inverse electron Green's function asĜ Above, ξ l (k) are electron energies at momentum k in band l, and we work in Nambu space [28] with Pauli matricesρ i . The mass enhancement Z(iω m ), superconducting order parameter φ(iω m ) and chemical potential χ(iω m ) do not acquire any momentum dependence due to the isotropic nature of the electron-phonon coupling. By using the electron self-energŷ with phonon propagator D(iq m−m ′ ) = D 0 (iq m−m ′ ) = −2Ω/(Ω 2 + q 2 l ) and boson frequencies q n = 2πT n, we derive a closed set of equations for Z(iω m ), χ(iω m ) and φ(iω m ) (see Appendix A).
The self-consistent results in Matsubara space are then analytically continued to the real-frequency axis, iω m → ω + iδ. This step is carried out in a formally exact and self-consistent manner, via the method first introduced by Marsiglio et al. [29] and extended for finite electron energy bandwidths as in Ref. [12], resulting in functions Z(ω), χ(ω) and φ(ω). With the real-frequency dependent electron Green's function at hand we can calculate the band and momentum resolved spectral function which can be compared to ARPES measurements when summed over band index l. By further summing up the momentum degree of freedom we have the means to compare our calculations to tunneling experiments, All our calculations are performed with the Uppsala Superconductivity (uppsc code) [12,[30][31][32][33][34][35]. For further details on the theory see also Refs. [12,27] and Appendix A. For sake of clarity, we note that the hereinvestigated phenomenon cannot be described within the Bardeen-Cooper-Schrieffer (BCS) model, which assumes a frequency-independent order parameter, absence of mass renomalization, and an electronic bandwidth W that is much larger than the phonon frequency, W ≫ Ω.
Results. We begin with a conceptually rather easy case of a nearest neighbor, one band tight-binding model on a 2D square lattice with bare electron energies ξ(k) = −2t[cos(k x )+ cos(k y )]− µ (we drop the band index). The hopping energies and chemical potential respectively are fixed at t = 0.425 meV and µ = −1 meV. Unless specified otherwise, we choose the electron-phonon scattering strength as g 0 = 2 meV and a relatively large phonon frequency Ω = 11 meV. With an electronic bandwidth of W = 3.4 meV the energies ξ(k) appear flat when compared to the phonon energy scale. We are not primarily interested in superconductivity, so, unless noted otherwise, we consider here T > T c corresponding to φ(ω) = 0. The replica bands under discussion occur outside the electron bandwidth of the flat band, hence not in a frequency regime where the superconducting energy gap alters the spectrum. Therefore superconductivity does not play any role for the current analysis, as we show in Appendix B. With this input we solve the full-bandwidth Eliashberg equations in Matsubara space and analytically continue the results to the real-frequency axis.
In Fig. 1(a) we show the self-consistent results for the mass renormalization Z(ω) = Z ′ (ω) + iZ ′′ (ω) in blue and chemical potential χ(ω) = χ ′ (ω) + iχ ′′ (ω) in red. Real and imaginary parts for both functions are drawn as solid and dashed lines, respectively. Within the electronic bandwidth, |ω| ≤ W/2, the real part of the mass renormalization takes on values close to unity, while χ ′ is in the range of µeV. This behavior is somewhat expected due to the relatively small coupling λ = 2g 2 0 N 0 /Ω ≃ 0.175 (N 0 : density of states at the Fermi level) and electron energy scale. In this frequency range Z ′′ (ω) and χ ′′ (ω) are negligible. Turning to |ω| > W/2, we observe a highly unexpected behavior of both Z(ω) and χ(ω). Apart from unusually large magnitudes in all four functions plotted in Fig. 1(a), we find large negative values for the mass renormalization at various frequencies. Although not straight-forward to physically interpret, it has been shown by Marsiglio and Carbotte that such values for Z(ω) can occur in the very strong coupling limit [36].
Next we look into the ARPES spectrum, which we compute from Eq. (3) with a smearing of δ = 0.01 meV. Our result for A(k, ω) is shown along high-symmetry lines of the Brillouin zone (BZ), and as function of frequency in Fig. 1(b). Here it is apparent that multiple nearly flat quasiparticle bands occur below and above the Fermi level. The frequencies corresponding to these rather coherent features seem to be separated by approximately the Einstein phonon frequency. To more reliably examine the energy positions of the observed replica bands, we show the logarithmic differential conductance in Fig. 1(c), where the frequency axis is normalized to Ω. The highest peak is observed at ω/Ω = 0 and corresponds to the one band electron dispersion. The next three peaks to both the left and right occur very accurately at multiples of Ω. All remaining signals can also be attributed to frequencies p Ω with p ∈ Z, but slightly shifted. Empirically we therefore conclude that a cascade of quasiparticle bands centered at p Ω exists, which represent replications of the original energy band. We provide a more rigorous proof of this argument in Appendix C, where we show that the positions of replica bands directly follow changes in Ω.
It is worthwhile investigating which part of A(k, ω) is responsible for producing the replica bands. For this purpose we use Eq. (1) to write the spectral function as: withξ(k, ω) = (ξ(k)+ χ(ω))/Z(ω) the renormalized electron energy dispersion. To make further progress, we analyze A(k, ω) for two different frequency regions: In one case the frequency lies within the bandwidth, |ω| ≤ W/2, and in the other we consider |ω| > W/2. Case |ω| ≤ W/2: From our numerical results we find that the imaginary partξ ′′ (k, ω) of the renormalized dispersion is to first order negligible in this frequency region. This is due to the fact that for |ω| ≤ W/2, The spectral function found from Eq. (5) in this case is and corresponds to the coherent part of the quasiparticle excitation spectrum. Since Z ′ (ω) is nearly constant for the frequencies under consideration, and its value is close to unity, we set for simplicity Z ′ (ω) = 1. Under this assumption Eq. (6) simplifies to A (1) (k, |ω| ≤ W/2) ≃ δ(ω −ξ ′ (k, ω)), so that the corresponding differential conductance can be approximated as To show that Eq. (7) provides the main contributions to the tunneling spectrum for |ω| ≤ W/2, we plot the result as red curve in Fig. 1(d). The delta-function is approximated as δ(x) ≃ exp(−x 2 /(2σ 2 ))/ √ 2πσ 2 with smearing σ = 0.01 meV. For comparison, dI/dV , as obtained by summing the non-simplified Eq. (5) over momenta, is drawn in blue. The boundaries of the renormalized electron dispersionξ ′ (k, ω) are indicated in yellow. We see that, despite the approximations made, the sum over delta-functions in Eq. (7) reproduces the full spectrum to very high accuracy for |ω| ≤ W/2. The spectral features for frequencies outside the electronic bandwidth, see Fig. 1(b), are due to the incoherent part of A(k, ω) as we will show in the following.
We can write the differential conductance, which results from A (2) (k, ω), as dI (2) The outcome of Eq. (9) is shown in Fig. 1(f) as dotted red curve. The full differential conductance as obtained from Eq. (4) is plotted in solid blue. We observe that the two curves fall precisely on top of each other, which shows explicitly that the replicas originate from A (2) only, and no contribution from A (1) enters for |ω| > W/2. As it turns out, we can reproduce the main features of the largefrequency spectrum by further simplifying Eq. (9). First, since |ω| > W/2, we can assume that ξ(k) has a comparatively minor influence on both, peak positions and amplitudes. Second, we might set the chemical potential renormalization to zero, which is a rather drastic simplification since neither χ ′ (ω) nor χ ′′ (ω) are negligibly small for |ω| > W/2, compare Fig. 1(a). However, assuming that these simplifications are valid we can write Note, that there is no longer a momentum dependence in Eq. (10) since we neglect ξ(k). The outcome of the above expression is shown in Fig. 1(f) in solid yellow. It is directly evident that the peak positions agree well with the full solution shown as blue curve. The heights do not precisely match the reference curve, which can be understood from the neglected '−χ ′′ (ω)' in the numerator of Eq. (9). Therefore we can conclude that the large-ω tunneling features are mainly mediated by the mass renormalization. For obtaining the correct intensities one needs to also include the chemical potential into the calculation. The bare electron dispersion ξ(k) plays a negligible role here. We note that our results show similarities to a study by Marsiglio and Carbotte [36], who investigated frequency dependent results of isotropic Eliashberg theory in the strong coupling limit. They showed the existence of quasiparticle-like excitations in the spectral function, that are located at ∆ 0 +p Ω, with ∆ 0 the superconducting gap edge. While we similarly find a sequence of replica bands, our results differ in that we are not depending on the limit λ → ∞, and the features detected in this work are independent of superconductivity, see Appendix B.
Twisted bilayer graphene. Let us now turn to TBG at a twist angle of ∼ 1.1 • , where we fix the phonon frequency at Ω = 11 meV [25,26] and use a faithful tenband tight-binding model for the electron energies [37]. This model has two flat bands near the Fermi level with a narrow bandwidth W ≈ 7 meV which are energetically separated from the rest of the bands by energy gaps over 20 meV. Further, we consider the normal state, T = 1.6 K > T c [15,16] and choose an electron-phonon scattering strength g 0 = 1.6 meV [27]. When perform-ing our analysis for the two flat bands only, the resulting spectra are very similar to our model calculations above and we present the outcomes in Appendix D. Next, we take into account a total of four energy bands with bandwidth W ≃ 127 meV > Ω. This includes the two flat bands close to the Fermi level, enclosing an energy window of around 7 meV < Ω, and an additional occupied and unoccupied band below and above. In Fig. 2 we show dI/dV for three different electron fillings n. Here, n (0) corresponds to half filling and for n (e) (n (h) ) the Fermi level lies exactly at the van Hove singularity of the unoccupied (occupied) band of the bare system. As apparent in panel (a), the replicas are superimposed with the additional non-flat energy bands, which are located at approximately |ω| 20 meV. However, the signals at p = ±1, ±2, ±3 are still clearly resolvable. The outcomes for the low-ω regime are drawn in Fig. 2(b). Comparing both panels reveals that the filling does not noticeably influence the intensity or location of the replicas, which is expected due to the comparatively small energy scale on which ξ l (k) is shifted.
Based on these findings, we are confident to predict the observation of signals at ω = p Ω in experiment. Since the replica bands do not represent coherent excitations of the system, their intensity is significantly lower than for actual poles of the Green's function (by about a factor O(50)). Consequently, the proposed tunneling features are not expected to be prominent in the experimental spectra however they should be detectable, given sufficient resolution. The observation of the here-predicted quasiparticle replica bands for TBG would provide strong support for the importance of the electron-phonon interaction for its low-temperature behavior.
Conclusions. To summarize, we predict a sequence of quasiparticle replicas in 2D systems that exhibit both, one or more flat bands around the Fermi level, and a comparatively large phonon energy scale. The conditions necessary to observe these features in ARPES or STS/STM measurements are a high experimental accuracy, and sufficient energy gaps between the flat band(s) and the remaining (un)occupied levels. The prime candidate for detecting such signals is TBG, where the flat bands close to the Fermi energy are isolated to a good approximation. Our calculations explicitly show that the replicas, occurring at multiples of the phonon frequency, are well distinguishable from the spectral signals of neighboring energy bands. Although TBG has been studied extensively, no such sequence of replicas has yet been discussed or observed experimentally, up to our knowledge. We argue that this is partially due to experimental resolution, since their rather low intensity might bury the replica signals in background noise. Furthermore, most experiments focus on a frequency range comparable to the bandwidth of the two flat bands, while for the observation of the here-discovered phenomenon an investigation of frequencies at least up to ∼ 50 meV is required. Hence, we predict that replica bands exist for TBG and should be observable under suitable conditions. This work has been supported by the Swedish Research Council (VR), the Röntgen-Ångström Cluster, the Knut and Alice Wallenberg Foundation (grant No. 2015.0060), and the Swedish National Infrastructure for Computing (SNIC). tors, with σ ∈ {↑, ↓} as spin label. The phonon frequency Ω and electron-phonon scattering elements g 0 are both assumed to be isotropic. The electron Green's function G l (k, iω m ) as defined in Eq. (1) obeys the Dyson equation withΣ(iω m ) given by Eq. (2) in the main text. The non-interacting Green's function in Eq. (A2) is (A4) As we describe in the main text, the phonon propagator is approximated by D(q, iq n ) = D 0 (q, iq n ), so that we get the electron-phonon interaction kernel (A5) The resulting Eliashberg equations for the mass renormalization Z(iω m ), chemical potential χ(iω m ) and superconducting order parameter φ(iω m ) in Matsubara space read The electron filling of the system is given by where L denotes the number of electronic bands.
Once we have solved the Eliashberg equations in Matsubara space, we can calculate ARPES and STM spectra to make direct contact with experiment [12]. For this purpose we analytically continue the solutions to Eqs. (A6-A8) self-consistently via with ζ(ω, z) = tanh ω−z 2T + coth z

2T
introduced for brevity (see [12,29]), and N l (0) the band-resolved density of states at the Fermi level. The solutions to Eqs. (A10-A12) are then used to compute the realfrequency matrix Green's function of the system, which in turn can be employed to find the momentum, band and frequency resolved spectral function Eq. (3), as well as the STS/STM spectrum via Eq. (4).
For solving Eqs. (A6-A8) and Eqs. (A10-A12) we do not make use of any further simplifications and keep the full complexity of the problem. Momentum and frequency grids have been checked for convergence. Our efficient implementation [30] makes use of fast Fourier transform (FFT) convolution schemes, and we exploit the known functional form of the interaction kernel to reach faster convergence in the number of Matsubara frequencies [38].
Appendix B: Results in the superconducting state As briefly mentioned in the main text, superconductivity does not play any significant role in the observation of replicated flat band(s). To explicitly prove this point we perform additional calculations for our model dispersion, using again the phonon frequency Ω = 11 meV. The electron-phonon scattering strength in this section is set to g 0 = 3 meV, such that we get a maximum superconducting gap of ∆ = φ(0)/Z(0) ≃ 215 µeV at T = 1.6 K. The corresponding tunneling spectrum for |ω| > W/2 is shown in Fig. 3 as blue curve. Once we increase the temperature to 2 K the gap closes and we obtain results for dI/dV as drawn in red. It is easily observed that the two spectra in Fig. 3 do neither differ in the locations of replica bands, nor are the corresponding intensities visibly deviating. This confirms that superconductivity has no observable effect on the tunneling features under investigation, since the gap opening occurs only inside the electron bandwidth of the flat band.

Appendix C: Influence of phonon frequency
In this section we want to examine the effect of phonon frequency on our results for the model system. In the main text we chose Ω significantly larger than the electronic bandwidth W = 3.4 meV, such that ξ(k) appears as flat band in comparison. Now we additionally consider the cases where Ω is 7, 3, or 0.5 meV, keeping T < T c and choosing the scattering strength g 0 , such that λ = 2N 0 g 2 0 /Ω = constant. Our results for the differential conductance are shown in Fig. 4(a), where the curves have been shifted vertically with respect to each other. We focus here on the large-frequency part of the tunneling spectrum and show the results only for ω > 0. The outcomes for Ω = 11 meV (blue), Ω = 7 meV (yellow) and Ω = 3 meV (red) rigorously prove that the replica bands occur at multiples of the respective phonon frequencies. In Fig. 4(b) we draw the complete tunneling spectrum for phonon frequency Ω = 0.5 meV and observe that the replica bands are absent. This behavior is to be expected, since for such small Ω the single electron band does no longer appear flat, i.e. the energy bandwidth is significantly larger than the phonon frequency.
Appendix D: Two-band case of TBG The replica bands in TBG are most easily observed when only the two flat bands close to the Fermi level are considered. In this section we therefore look into results obtained by this setup, additionally setting T > T c , Ω = 11 meV and g 0 = 1.6 meV. The electron bandwidth W ≃ 7 meV is then smaller than Ω, hence we expect to observe effects comparable to our model system in the main text and Appendix C. As before, we solve the Eliashberg equations in Matsubara space and analytically continue the self-consistent results to real frequencies. The ARPES and tunneling spectra are then obtained as function of ω.
In Fig. 5(a) we show the spectral function, summed over energy bands, at a filling n = n (e) along momenta and frequencies. The bare two-band dispersion corresponds to enhanced signals close to the Fermi level. All remaining features in this graph represent replications of the original energies. As is easily seen, and in agreement to our model calculations in the main text, these signals are almost constant in momentum space and occur approximately at integer multiples of the phonon frequency. In panels (b) and (c) of Fig. 5 we zoom into the frequency regions around ω/Ω ∼ 3 and ω/Ω ∼ −3, respectively. From these parts of the spectrum it is evident that the observed features are indeed direct replications of the original ξ l (k).
We show in Fig. 6 the corresponding differential conductance as it can be measured by STM experiments. For better visibility we plot in Fig. 6(a) the high-frequency spectra only, for fillings n (h) (blue), n (e) (red) and n (0) (yellow) as given in the main text. Multiples of Ω are indicated by dashed gray vertical lines. The region of small ω is shown in a magnified way in Fig. 6(b) for similar fillings and color code. We observe from this graph that the large-energy part of the spectrum contains a sequence of replicas similar to our model system. Further, the range of electron doping that is relevant for super- conductivity [27] does not affect the replica positions in a noticeable way, since the phonon energy scale is dominant. It should be noted that not all signals are equally well pronounced, so only some replicas could be be ob-servable due to experimental resolution.