Strong electron–phonon coupling in magic-angle twisted bilayer graphene

Sample fabrication

MATBG devices A, B and C are devices B, F and C from ref. ¹¹ (the full description of the fabrication procedure can be found in this reference). In brief, devices were fabricated using a ‘tear-and-stack’ method⁵¹ in which a single graphene sheet is torn in half by van der Waals interaction with hBN. The two halves were rotated relative to each other and stacked to form MATBG. Graphene, graphite (in the aligned device) and hBN were picked up with polyvinyl alcohol. Then, to flip the heterostructure upside down, the heterostructure was pressed against an intermediate structure consisting of polymethyl methacrylate/transparent tape/Sylgard 184, and the polyvinyl alcohol was dissolved by water injection. The heterostructure was then transferred to a SiO₂/Si chip with pre-patterned titanium/gold electrodes. Residual polymer was dissolved in dichloromethane, water, acetone and isopropyl alcohol. These chips were annealed in ultrahigh vacuum at 170 °C overnight and at 400 °C for 2 h in previous STM and STS measurements¹¹. The other devices shared a similar fabrication and preparation method.

Spatial- and angle-resolved photoemission spectroscopy

Synchrotron-based μ-ARPES measurements were performed at Beamline 7.0.2 (MAESTRO) of the Advanced Light Source, USA and Beamline 07U of the Shanghai Synchrotron Radiation Facility (SSRF), China. The samples were annealed in ultrahigh vacuum at 300 °C for 3 h and measured under ultrahigh vacuum below 3 × 10⁻¹¹ torr. The photon energy of the incident beam was 95 eV and the measurement was performed at room temperature. Data were collected using an R4000 analyser upgraded with deflectors (Advanced Light Source) and a DA30 analyser (SSRF). The incoming photon beam was focused down to a 2-μm spot size by using a capillary mirror²⁸ and a 1-μm spot size using a Fresnel zone plate (SSRF). The total energy and angle resolutions were 20 meV and 0.1°, respectively.

Laplacian (sum of second partial derivatives) plots were used to highlight the non-dispersive bands, as presented in Figs. 2 and 3. These plots were generated by summing the second partial derivatives of the original ARPES spectra along both axes in pixel units, then transforming them into energy and momentum space.

Model calculation

Estimation of EPC

We adopted the following non-interacting-electron tight-binding Hamiltonian to calculate the electronic structure and EPC of MATBG

where \({c}_{{I}i\alpha }^{\dagger }\) and \({c}_{{J}j\beta }\) are creation and annihilation operators for the p_z orbital of the iα/jβ carbon atom in the I/Jth moiré superlattice (α and β are joint indices for sublattices and layers), and R is the corresponding atomic coordinate. The impact of atomic coordinates on the hopping parameter t is approximated by the Slater–Koster formula, whose dependence on atomic coordinates sets the stage for the EPC estimation. The Slater–Koster parameters are specified in Supplementary Information section VII.

To evaluate the coupling in the mini-Brillouin zone, we further projected the coupling matrix using the truncated atomic plane wave (TAPW) method^43,52. Thus, a moiré phonon near the Γ_M point can be approximated by in-plane transverse optical/in-plane longitudinal acoustic/in-plane longitudinal optical phonons at graphene K/K′ points (see Supplementary Information section VIII for details). We set the moiré phonon energy ω₀ = 150 meV, the mass of carbon atoms m_c = 2.0 × 10⁻²⁶ kg, and the characteristic phonon length \({l}_{{\rm{p}}}=\sqrt{\hbar /(2{m}_{{\rm{c}}}{\omega }_{0})}=34.0\,{\rm{m}}{\text{\AA }}.\) The TAPW electrons and projected EPC constitute the model Hamiltonian for the TBG system

where \({{\bf{c}}}_{\bar{{\bf{k}}}\sigma }^{({\rm{TAPW}})}\) is the column vector of electron annihilation operators with moiré momentum \(\bar{{\bf{k}}}\) and spin σ in the TAPW basis and \({h}_{\bar{{\boldsymbol{k}}}}\) is the electronic hopping matrix for a specific moiré momentum \(\bar{{\bf{k}}}\). The size of its moiré mini-Brillouin zone is denoted as N_m to distinguish from the system size N. For folded phonon branches, \({a}_{\bar{{\bf{q}}}\nu }\) is the annihilation operator with moiré momentum \(\bar{{\bf{q}}}\) and ν labels the index of branches. The average distance between phonons and electrons is small compared with the moiré length scale. Therefore, the coupling matrix \({M}_{\bar{{\bf{q}}}\nu }\) can be estimated by the momentum-independent M₀, which can be evaluated by the frozen-phonon calculations.

As the experimentally relevant electrons lie in low-energy flat bands, we further project the TAPW orbitals onto the four flat-band orbitals, which can be expressed as the projection operator \({{\bf{c}}}_{\bar{{\bf{k}}}\sigma }^{({\rm{flat}})}={{P}_{\bar{{\bf{k}}}}^{\dagger }{\bf{c}}}_{\bar{{\bf{k}}}\sigma }^{({\rm{TAPW}})}\), where P is the flat-band projection operator. In the projected Hamiltonian, the EPC matrix becomes a 4 × 4-dimensional \({\widetilde{{M}_{\bar{{\bf{k}}}}}(\bar{{\bf{q}}})=P}_{\bar{{\bf{k}}}{\boldsymbol{+}}\bar{{\bf{q}}}}^{\dagger }\,{M}_{\bar{{\bf{q}}}}\,{P}_{\bar{{\boldsymbol{k}}}}\). On the basis of its numerical distribution (Supplementary Information section VIII), we find that \(\widetilde{{M}_{\bar{{\bf{k}}}}}(\bar{{\bf{q}}})\) can be approximated by \({Q}_{\bar{{\bf{k}}}{\boldsymbol{+}}\bar{{\bf{q}}}}\,{g}_{\bar{{\bf{q}}}}\eta {Q}_{\bar{{\bf{k}}}}^{\dagger }\), where Q is the similarity transformation, g is the coupling strength and η is a constant diagonal matrix. Thus, the projected Hamiltonian can be written in an EPC-diagonal basis

where \({{\bf{d}}}_{\bar{{\bf{k}}}\sigma }={Q}_{\bar{{\bf{k}}}}^{\dagger }{{\bf{c}}}_{\bar{{\bf{k}}}\sigma }^{({\rm{flat}})}\) and the diagonal matrix \({\varepsilon }_{\bar{{\boldsymbol{k}}}}\) represents the energy of the flat bands.

ARPES spectral simulation for TBG

The strong EPC for the flat-band electrons leads to non-perturbative polaronic dressing effects. We consider the Lang–Firsov transformation for the coupled Hamiltonian

where \({{\bf{d}}}_{\bar{{\bf{R}}}\sigma }\) is the real-space annihilation operator of electron in the EPC-diagonal basis. Owing to the separation of energy scales for electrons and phonons, we can employ the polaron ansatz for the ground-state wavefunction \(| {\varPsi }_{{\rm{G}}}\rangle ={U}_{{\rm{LF}}}^{\dagger }| {\psi }_{{\rm{e}}}\rangle \bigotimes | {\psi }_{{\rm{ph}}}\rangle \), where |ψ_e⟩ and |ψ_ph⟩ are the electronic and phononic wavefunctions. Moreover, as both the transformed coupling strength \({g}_{\bar{{\bf{q}}}}\) and temperature are much less than the phonon energy ω₀ = 150 meV, we further assume that the \(|{\psi }_{{\rm{ph}}}\rangle \) can be approximated by a vacuum state \(|{0}_{{\rm{ph}}}\rangle \). Thus, the electronic part \(| {\psi }_{{\rm{e}}}\rangle \) is determined by an effective Hamiltonian

The \({h}_{\bar{{\bf{R}}}\bar{{{\bf{R}}}^{{\prime} }}}^{* }\) is the phonon-dressed electronic hopping matrix, whose specific form is shown in Supplementary Information section X. In the flat-band limit of TBG, this matrix is close to identity and, therefore, is irrelevant in determining the polaronic dressing. The ground state of the above equation determines the \(| {\psi }_{{\rm{e}}}\rangle \) in \(| {\varPsi }_{{\rm{G}}}\rangle \).

The photoemission spectrum also involves excited states (denoted as \(| \varPhi \rangle \))

where Γ is the Lorentzian broadening. With the aforementioned ground-state ansatz, the intensity of the Mth replica peak is explicitly determined as

It follows a Poisson distribution with the factor \(\sum _{\bar{{\bf{q}}}}\frac{{N}_{v}\,{| {g}_{\bar{{\bf{q}}}}| }^{2}}{{N}_{{\rm{m}}}\,{\omega }_{0}^{2}}\approx 0.11\), according to the frozen-phonon simulations. Focusing on the relative intensity of replica features and ignoring the interacting nature of electrons inside each replica, we produce the spectral simulation in Fig. 4d.

Variational non-Gaussian ansatz for more general models

The Lang–Firsov transformation is suitable for the flat-band TBG model. The experiments presented in this paper also include TBG under more complicated conditions, including with hBN substrates and finite bandwidth away from the magic angle. These generalized cases can be modelled in the form of

where \({h}_{\bar{{\bf{k}}}}\) describes the band structure in general. The derivation and specific forms of Hamiltonian with hBN and bandwidth are detailed in the Supplementary Information sections X and XI.

It is important to note that the first term in equation (8) does not commute with the second term in this generalized case, making the Lang–Firsov transformation unsuitable for this case. To simulate the polaronic dressing in this generalized Hamiltonian, we employ a variational ansatz of the ground state \(| {\varPsi }_{{\rm{G}}}\rangle ={U}_{{\rm{NGS}}}^{\dagger }(\lambda )| {\psi }_{{\rm{e}}}\rangle | {0}_{{\rm{ph}}}\rangle \) with

where λ is the variational parameter in contrast to the fixed g/ω₀ of the Lang–Firsov transformation for TBG. Owing to the high phonon energy compared with any energy scales in equation (8), we still assume the post-transformation phonon state is vacuum. Therefore, the variational ansatz gives the total energy as a function of λ:

with the normalized kinetic energy E_kin follows the expression in Supplementary Information sections X and XI. Unlike the analytically solvable Lang–Firsov transformation, the variational parameter λ is obtained by numerical optimization of E_tot(λ). The self-consistent equations for different situations of the generalized Hamiltonian are derived in Supplementary Information sections X and XI. This numerically determined λ describes the relative strength of the polaronic dressing, where λ reproduces g/ω₀ when H takes the flat-band TBG form in equation (3). The Poisson factor for the relative replica intensity is then obtained by p = N_vλ².