Quantum twisting microscopy of phonons in twisted bilayer graphene

Cryogenic (T = 4 K) QTM

The cryogenic QTM consists of two nanopositioner towers facing each other. One tower is equipped with three translational degrees of freedom (XYZ) and carries a flat sample. The opposing tower holds an AFM cantilever and has a rotational degree of freedom (θ), as well as two lateral degrees of freedom for positioning the AFM tip at the centre of rotation. The entire assembly is cooled down in vacuum to liquid helium temperatures. In the experiment, we bring into contact a vdW heterostructure on the tip with a flat vdW heterostructure on the bottom sample, creating a two-dimensional interface, typically several hundred nanometres across in both directions. vdW attraction between the two heterostructures self-cleans the contact area, leading to a pristine interface. The tip and sample remain in continuous contact throughout the experiment, including during any twisting or scanning operation. Self-sensing AFM cantilevers equipped with piezoresistive elements are used for monitoring and maintaining a constant force, ensuring that the area of the contact interface remains constant during the scans. The QTM junctions in this paper do not have a tunnelling barrier separating the two conducting sides (apart from Fig. 3b,c, in which we used defects in the WSe₂ barrier to image the tip shape). This means that the resistance of the tunnel junction, especially at low twist angles or at high bias, can be comparable with or lower than the resistance of the contacts or of the bulk of the two-dimensional layers leading from the contact to the junction. We obtain this contact resistance from measurement close to 0° twist angle and then directly calculate how much bias decreases on the contact resistance and how much is decreasing on the junction, V_b, which is the variable we use throughout the paper (Methods section ‘Removing contact resistance in the two-probe measurement’).

Fabrication of the vdW-devices-on-tip

The fabrication of the vdW-devices-on-tip follows the procedure we described previously¹. In brief, we use tipless AFM cantilevers that have a thin metallic line reaching its end, which we use for electrically connecting to the active vdW layer. We make a custom tip using focused-ion-beam-deposited platinum pyramid, with base dimensions of 2 × 2 μm and a height of 1–2 μm. The vdW layers are transferred onto the pyramid using the polymer membrane transfer technique.

Conductance measurements

Voltages are applied using a custom-built digital-to-analogue converter consisting of a.c. + d.c. sources. Current is measured using a FEMTO amplifier, followed by an NI sampler. We use a Zurich Instruments lock-in amplifier (MFLI) for force feedback of the piezoresistive AFM cantilevers.

High-resolution elastic momentum-resolved tunnelling near the commensurate angle of 21.8° and the upper bound it puts on strain in the tip and sample

In Extended Data Fig. 1, we show high-resolution momentum-resolved imaging experiments conducted near the commensurate angle of 21.8°. These cryogenic experiments, performed during the same cooldown and using the same sample and tip as in the experiments in the main text, reveal an exceptionally sharp momentum-resolved map that aligns well with theoretical predictions. Below we explain how this measurement allows us to put a tight upper bound on the strain of the graphene layers on the tip and flat sample and show that the strains are very small, less than 0.1%.

Our experiment involves a TBG interface that is formed in situ, by bringing into contact a graphene layer on a tip and one on a flat substrate. It is instructive to ask whether this interface has substantial strains, owing to a strain present in one or both participating layers. Although monolayer graphene typically exhibits minimal strain when placed on a flat hBN substrate, we may imagine that placing graphene on the rough topography of the tip could introduce substantial strain, particularly near its apex.

To estimate the strain, we use elastic momentum-conserving tunnelling experiments. In such experiments, performed in a QTM geometry, we map the electronic energy bands by rotating two samples relative to each other and measuring the tunnelling current as a function of the rotation angle and bias¹. Because current flows between the layers only when their states match in both momentum and energy, a two-dimensional θ–V_b scan effectively maps the energy bands of the system. For these experiments, we typically use a thin insulating barrier between the two layers, ensuring that the experiment remains in the tunnelling regime and allowing for the application of a well-defined bias across the junction. However, in this study, the two graphene layers are in direct contact, with no barrier material in between. Notably, when the experiment is conducted around a twist angle of 21.8°, we still observe elastic momentum-resolved current and the experiment is operating in the tunnelling regime even without a barrier. To understand this, we recall that, although at a twist angle of 21.8° the Fermi surfaces of the two graphene layers do not overlap in the first BZ, they do overlap in the third extended BZ^1,37,39, as shown in Extended Data Fig. 2. Namely, the momentum-resolved tunnelling is performed between very high, in-plane momentum components of the wavefunction, which correspondingly decay much more rapidly in the z direction. As a result, despite the direct contact between the layers, the experiment operates in the tunnelling regime.

Extended Data Fig. 1a–c shows the measured tunnelling current and its derivatives versus twist angle and bias. The dashed lines in these panels correspond to theoretically calculated crossing conditions between the Fermi surface in each layer and the empty energy bands in the other layer, as illustrated in Extended Data Fig. 1d. The agreement between the measurement and theory is notable, given the fact that the theory is performed without free parameters (we used the capacitance values measured independently in the main text).

Even a tiny strain in the sample or in the tip would make the measurement look completely different—without strain, the BZ of graphene is a perfect hexagon with its six Dirac cones exactly 60° apart from each other. With strain, this perfect hexagon deforms and the Dirac points are no longer equidistant in angle. For example, 1% strain will make the angle between adjacent Dirac points vary between 59.5° and 60.5°. Remembering that our experiment examines momentum-conserving tunnelling collected simultaneously from all six Dirac points in the tip matching their corresponding six Dirac points of the bottom sample, it is clear that, if the layers are strained, we will see several copies of the crossing conditions, displaced along the angle axis. Each copy will result from one of the six Dirac points crossing its corresponding Dirac point in the other layer. The same holds for the experiments performed near 21.8°, at which the six crossing Dirac points are at the corners of the third BZ (Extended Data Fig. 2). Indeed, occasionally, we obtain tips that have strains and then we can clearly resolve such multiple copies. The fact that in the measurement in Extended Data Fig. 1 we see only one copy means that all six Dirac points are crossing at approximately the same twist angle. In fact, if we assume that the angular width of the lines in the measurement comes from smearing owing to several copies crossings, then we can put an upper bound on the strain in the tip and sample. In the measurement, we can identify lines whose full width at half maximum is approximately 0.05°, which puts an upper bound of 0.1% on the strain in our experiments.

Flat samples and QTM tips

Extended Data Fig. 4 shows the optical image of vdW devices used in the experiment. Extended Data Fig. 4a shows a monolayer graphene QTM tip. Extended Data Fig. 4b,c shows the graphite sample used in Fig. 1 and multilayer graphene (MLG)/hBN/graphite sample used in Figs. 2 and 3.

Symmetry of phonon modes in voltage bias

Extended Data Fig. 5 shows the mirrored I–V curve (in voltage bias) of the graphite–graphite tunnelling junction (as shown in Fig. 1d). The phonon peak positions appear at the same energy for both positive and negative voltage biases. However, there is an asymmetry in phonon peak intensity, which may arise from the complex, electron–hole asymmetric band structure of graphite.

Imaging the tip contact area in situ and determining the pressure in the experiment

The QTM geometry creates a clean two-dimensional interface between the QTM tip and the sample side. We can image this interface by scanning the tip laterally across a region of the sample that has a transition metal dichalcogenide (TMD) layer with atomic defects (three-monolayer-thick WSe₂ in the current experiment, placed adjacent to the area in which the main experiment is performed). The defects provide an extra channel for the tunnelling current (Extended Data Fig. 6a). Thus, whenever the tip overlaps with a defect, we observe an increased current. As a result, in a measurement of the current as a function of lateral scanning (X and Y coordinates), each defect produces an image of the tip’s contact area. Extended Data Fig. 6b shows such imaging, in which we resolve several copies of the bulk and outline of the tip. Knowledge of the exact area of the tip is crucial for the quantitative analysis of the EPC strength. The left panel of Extended Data Fig. 6b shows the tip contact area of the hBN-backed monolayer graphene tip used in this work. The area of tip 1 is 37,602 nm² and that of tip 2 is 18,106 nm². We apply a d.c. bias close to the band edge (on the order of ±1 V for WSe₂); in this energy window, many defects are active for tunnelling.

In the experiment, we use the built-in AFM in our cryogenic QTM to measure and maintain a constant force throughout the experiment. We generally apply small forces F ≈ 100 nN. Combined with the area of tip 1 mentioned above, this amounts to a pressure on the interface of P ≈ 2.5 MPa, which is much smaller than the vdW adhesion pressure, which is on the order of 1 GPa.

Theory model for inelastic tunnelling through phonon emission and the two mechanisms of EPC

The phonon-assisted inelastic tunnelling is modelled using Fermi’s golden rule:

in which we assume the interlayer tunnelling current δI_inter is flowing from the tip (top layer) to the sample (bottom layer). δI_inter results from the tunnelling of an electronic state |k, s⟩ in the top layer to another electronic state |p′, s′⟩ in the bottom layer, with the emission of a phonon mode r and momentum Q. s and s′ are the band index of graphene layers. E_k,s and ω_r,Q are the corresponding electronic state energy and phonon mode energy. N_f = 4 accounts for the spin and valley degeneracy. f_T(E_k,s) and f_B(E_p′,s′) are the Fermi–Dirac distribution functions for the top and bottom layers. The chemical potential for the top and bottom layers and their relative energy shift are determined from the electrostatic model (Methods section ‘Electrostatic model for the QTM junction’). Here we denote Q for the phonons emitted in the top layer and Q′ for the phonons emitted in the bottom layer. For each layer, there are three possible Q vectors contributing to the tunnelling current in the leading order, which correspond to three Bragg scatterings within the first BZ.

The above formula is a general expression for phonon-assisted tunnelling and the matrix element ⟨k, s|T_in-el|p′, s′; r, Q⟩ contains two types of EPC⁴⁹:

H_inter is the interlayer EPC, for which the phonon field affects the interlayer tunnelling amplitude t_⟂ and gives rise to the coupling between the layers, as shown in Fig. 2h. The other two terms correspond to the in-layer mechanism shown in Fig. 2g. The microscopic expression includes: H_T, accounting for the pure electronic interlayer tunnelling; G₀, the Green’s function of the uncoupled layers accounting for the energy denominators of the virtual intermediate states; H_in-layer, the in-layer EPC originating from the changes of hopping amplitudes in the layer. The detailed derivation of the two types of EPC is shown in ref. ⁴⁹. Here we give the final expressions of the interlayer EPC:

These two coupling matrices are for phonons emitted in the top and bottom layers, respectively. Inside the square brackets in equations (4) and (5), the first term describes in-plane phonons and the second term describes out-of-plane phonons (both ZA mode and ZO mode). α and β are sublattice indices and d_z is the interlayer distance. ω_r,Q and \(\boldsymbol\varepsilon _r,\bfQ^\alpha \) are the phonon energy and polarization vector, respectively, for phonon mode r at momentum Q. M is the total mass in the unit cell. N is the number of unit cells. K (K′) and G_j=1,2 (\(\bfG_j=1,2^\prime \)) are the Dirac point and reciprocal lattice vectors of top (bottom) layer, respectively. G_j=0 = 0 corresponds to the no Bragg scattering case. d is the shift vector between the two layers. q_j is the vector connecting the Dirac points between the top and bottom layers, for which j = 0, 1, 2 accounts for three Bragg scatterings within the first BZ. They correspond to three emitted phonon vectors, as seen by the delta functions \(\delta _\bfQ,\bfq_j\) and \(\delta _\bfQ^\prime ,\bfq_j\) in equations (4) and (5). Note that q₀ = q_M in our notation in the main text. The interlayer electronic tunnelling matrix \(T_j^\alpha \beta \) is expressed as

in the sublattice space and \(\zeta =\frac2\pi 3\). For the interlayer mechanism, the theory predicts: (1) the dominant twist angle dependence in equations (4) and (5) is coming from the ZPM factor \(\sqrt{\frac\hbar 2M\omega _r,\bfQ}\), which results in an enhanced EPC as the twist angle decreases for the acoustic modes, whereas for optical modes, the coupling approximately remains constant. (2) There is a geometric factor between the K–G_j vector and phonon polarization vector, as seen by the term \((\bfK^\prime -\bfG_j^\prime )\cdot \boldsymbol\varepsilon _r,\bfQ^\alpha \) in equation (4) and \((\bfK-\bfG_j)\cdot \boldsymbol\varepsilon _r,\bfQ^\prime ^\beta \) in equation (5). For the transverse mode, this dot product gives \(\cos \left(\frac\theta 2\right)\) factor, in which θ is the twist angle, whereas for the longitudinal mode, it gives \(\sin \left(\frac\theta 2\right)\) factor, which suppresses the EPC. This explains the missing LA mode in the experiment. In the experiment, the measured EPC is averaged over the Fermi surface of the top and bottom layers, and we define it as:

in which θ_k and θ_p′ are the polar angle of the electron momentum relative to the top and bottom layer Dirac point, respectively. Here we can change the phonon field basis of the top (Q) and bottom (Q′) layers into layer-symmetric and antisymmetric mode and account \(| g_\rminter,r,\bfQ=\bfq_\rmM| \) as the EPC only for one phason mode, as the layer-symmetric mode does not couple to tunnelling electrons.

As for the second mechanism, the in-layer EPC is given by⁴⁹:

in which β₀ = ∂t_∥/∂r_∥ and \(\widehat\bfe_j\) is the unit vector along the in-plane carbon–carbon bond direction. The in-layer contribution to the tunnelling current is given by \(T_\textin-el^\textin-layer=H_\rmTG_H_\textin-layer+H_\textin-layerG_H_\rmT\). After averaging over the Fermi surface, the matrix element is approximated by⁴⁹:

in which the factor \(\alpha =\fract_\perp \hbar \nu _\rmD\bfq_\rmM\) is from the virtual tunnelling process and \( g_\textin-layer,r,\bfQ=\bfq_\rmM^2\) is the Fermi surface averaged in-layer EPC. As q_M gets smaller, \(g_\textin-layer,r,\bfQ=\bfq_\rmM\) for the acoustic mode tends to zero as \(\sqrtq_\rmM\), whereas the coupling to the optical mode remains constant. Compared with the phason mode, we have the extra N_layer = 2 factor coming from phonons emitted in the top and bottom layers.

In the numerical calculations for the interlayer and in-layer EPC in Fig. 4, we use the following parameters⁵⁰: t_⟂ = 0.1 eV, ν_D = 10⁶ m s⁻¹, β₀ = 4.22 eV Å⁻¹. The phonon dispersion and polarization vectors are solved by diagonalizing the dynamic matrix of graphene⁵¹. The numerical calculations are shown in Extended Data Fig. 7, for Fermi surface averaged interlayer \(| g_\textinterlayer,r,\bfq_\rmM| \) and intralayer \(\alpha \sqrtN_\rmlayer| g_\textin-layer,r,\bfq_\rmM| \) contributions as a function of rotation angle. To simplify notation, we denote them as |g_inter,r| and |g_in-layer,r| in Extended Data Fig. 7. For the TA mode, the interlayer mechanism dominates over the in-layer mechanism. For the LA mode, the in-layer process is stronger than the interlayer one, owing to the geometric factor that suppresses the interlayer coupling. The overall coupling strength of the LA mode at the momentum q_M that the QTM tip examines is weak compared with the TA mode. For TO and LO modes, we plot both the separate and averaged value in Extended Data Fig. 7c. The averaged contribution is dominated by in-layer coupling contributions, which become stronger as the twist angle decreases. The experiment measures the EPC to the TO and LO modes together, near the K point, along an arc in the momentum space. Along this trajectory, the TO and LO modes have an avoided crossing, as seen in Extended Data Fig. 7d. The corresponding bare in-layer EPC of TO and LO also differs, as shown in Extended Data Fig. 7e. Note that, near the K point, the TO phonon has larger bare EPC than the LO phonon. In Extended Data Fig. 7f, we further decompose the TO and LO contribution in terms of \(\alpha ^2N_\rmlayer g_\textin-layer,r,\bfq_\rmM^2\) and \(^2\) (r = TO, LO) and compare with the data: ΔG_optical/2βν_Tν_BA_tip from Fig. 4a. In the experiments, we cannot separate LO and TO phonon contributions, and the averaged value is consistent with the theory estimate, for which TO phonon has a stronger contribution. In short, both the interlayer and in-layer mechanism are, in general, the same order of magnitude. However, parametrically, the interlayer mechanism is much stronger for the acoustic TA mode, whereas the in-layer mechanism dominates for the optical modes.

Extracting EPC from the measured inelastic conductance steps

Equation (2) can be simplified (derivations are detailed in ref. ⁴⁹). Here we show the result:

in which Ω is the tunnelling area of the twisted junction and, in our experiment, Ω = A_tip, imaged directly by atomic defects (see Methods section ‘Imaging the tip contact area in situ and determining the pressure in the experiment’). ν_T and ν_B are the DOS of the top and bottom layers, respectively. θ(eV_b − ω_r,Q) is the step function, which onsets the conductance \(\frac\rmdI_\rmin-el\rmdV\) when the bias voltage reaches the phonon energy. We defined the angular averaged coupling:

in which N is the number of unit cells in the junction. Inserting into equation (10), we have:

in which \(A_=\fraca^2\sqrt32\) is the unit-cell area and ν_T and ν_B are the DOS on the top and bottom layers, respectively. Because the phonon wave vector Q is related by the C₃ symmetry (as seen by \(\delta _\bfQ,\bfq_j\) factor in equations (4) and (5)), we have the conductance step owing to the emission on a phonon mode r as:

in which \(\beta =\frac2\pi e^2\hbar N_\rmfN_\rmb\,A_\), such that N_b = 3 accounts for the Bragg scatterings and N_f = 4 accounts for the flavour degeneracy (spin/valley). Given the two types of microscopic EPC, \(^2\) includes:

For the in-layer mechanism, N_layer = 2 because phonons in both layers contribute. For the interlayer mechanism, only the layer-antisymmetric phason mode contributes. For extracting the in-layer coupling in Fig. 4, we use the parameters t_⟂ = 0.1 eV and ν_D = 10⁶ m s⁻¹ for estimating the value of the α factor.

Electrostatic model for the QTM junction

The QTM junction is modelled by a three-capacitor model¹. A graphene layer on the tip side and a graphene layer on the sample side are separated by a vacuum gap. The third layer is the graphite back gate, separated by a hBN spacer on the bottom sample side. We have the following electrostatic equations:

in which V_b is the voltage bias applied between the two graphene layers and V_bg is the back-gate voltage with respect to the bottom graphene layer. n_B (n_T) and μ_B (μ_T) are the carrier density and chemical potential, respectively, of the bottom (top) layer graphene. C_g is the geometric capacitance per unit area between two graphene layers and C_bg is that between the back-gate and the bottom graphene layer. We use the values C_g = 11 μF cm⁻² and C_bg = 38.6 nF cm⁻² to fit the data. From the electrostatic model, we plot the charge-neutral Dirac point lines (dashed red and blue lines) in Fig. 3d,e and Extended Data Figs. 9 and 11.

Bulk-graphite phonon dispersion model

In the main text, the phonon dispersion of graphite is modelled using a dynamic matrix through the force constant approaches. For the in-layer force constants, we use the parameters from ref. ⁵¹ (GGA calculation, column 3 in table 3 in that reference), which includes the fourth nearest-neighbour coupling that fits to DFT calculations. To obtain the bulk-graphite phonon spectrum, we further include the interlayer mechanical coupling in the dynamic matrix, as detailed in ref. ⁵². We include only the nearest-neighbour out-of-plane radial and tangential coupling, denoted by \(\widehat\phi _\rmr\) and \(\widehat\phi _\rmt\), and rescale these two parameters together to match the ZA phonon gap measured in experiments, giving the values \(\widehat\phi _\rmr=\mathrm3,296\,\rmdyn\,\rmcm^-1\) and \(\widehat\phi _\rmt=-\mathrm7,121\,\rmdyn\,\rmcm^-1\). For other phonon modes, introducing these interlayer coupling parameters generates an energy splitting of less than 0.5 meV. In Extended Data Fig. 8, we compare the spectrum obtained from ref. ⁵² (Extended Data Fig. 8b) to the model described above (Extended Data Fig. 8a) and we can see that the latter fits the experiments much better.

Further data: G as a function of V
_b and V
_bg measured at two more twist angles, θ = 22.7° and 9.4°

In Fig. 3, we show the junction conductance as a function of V_b and V_bg at twist angle θ = 16.8°. Extended Data Fig. 9 shows the same map at two more angles, θ = 22.7° and 9.4°. These figures show similar features as in Fig. 3d,e: horizontal lines owing to inelastic phonon emission processes and conductance dips tracing the neutrality points of the top and bottom graphene layers. The dashed blue and red lines correspond to the neutrality points calculated using the electrostatic model for the QTM junction.

Further data from different samples and QTM tips

Here we show further measurements performed with different bottom graphene samples and QTM tips to those shown in the main text. Our scanning set-up allows us to create the twisted interface at different locations along the bottom graphene, which are separated by several micrometres. Because the tip contact area is only several hundreds of nanometres, measuring at such spatially separated positions effectively means that we have switched to a completely different bottom sample. Also, within the QTM set-up, we have the capability to modify in situ the tip shape. One example is shown in Fig. 3b, in which the area of the tip is changed by approximately a factor of two (tips 1 and 2). Because the area added to the tip is substantial and can potentially have different strains, edges and so on, the measurements with these two tips effectively amount to experiments performed with different tips. In Fig. 2, we show the phonon spectrum measured using tip 1. Extended Data Fig. 10a,b shows further inelastic phonon spectroscopy data obtained using tip 2, which was also measured in a different sample location, overlaid with the theoretical phonon dispersion of bulk graphite. This scan was performed only down to 10°, because, at this angle, the tip caught some dirt. Performing the same analysis as in Fig. 4, we extract the electron–phason and electron–optical couplings and compare them with the data in Fig. 4 over the angular range that they overlap. Extended Data Fig. 10c shows this comparison for the extracted intervalley optical TO/LO EPC and Extended Data Fig. 10d shows the comparison for the electron–phason coupling. We can see that there is a good agreement both in terms of the magnitudes of the coupling as well as on their twist angle dependence.

Extended Data Fig. 11 shows further measurements performed with a completely different tip and sample in a different cooldown. Extended Data Fig. 11a,b shows the graphene–graphene junction conductance G and d²I/dV² and as a function of V_b and V_bg at twist angles θ = 30° and 20°. We can clearly resolve in both scans horizontal lines that correspond to inelastic processes with the different phonon modes. We also see that the step in conductance changes gradually with back-gate voltage. Performing a similar analysis as in the main text (Fig. 3g) and plotting the height of the conductance steps for the TA and ZA modes as a function of gate voltage (Extended Data Fig. 11c), we clearly resolve the predicted linear dependence. Moreover, unlike the sample in the main text, which could not be measured at large negative back-gate voltage owing to a gate leak, in the sample below, we also explored the hole side and we can clearly see how that, also there, the dependence is linear with a slope with inverted sign (dashed lines in Extended Data Fig. 11c are the model calculations).

Removing contact resistance in the two-probe measurement

In the experiments, we apply a voltage V_applied between the contacts connected to the vdW-devices-on-tip and to the flat sample side and record the current I flowing through the twistable junction and its derivative, \(\widetildeG=\frac\rmdI\rmdV_\rmapplied\). However, as these are two-terminal measurements, there will be a contact resistance R_c leading to/from the twistable junction area. Therefore, the voltage bias V_b falling on the junction will be smaller than V_applied. The difference between V_applied and V_b becomes notable when the junction resistance is low (either at high applied voltage or at low twist angles). To remove this effect, we rotate to zero twist angle (for which the resistance of the junction is the lowest) and record the resistance R_c (about 50 kΩ for the monolayer graphene–monolayer graphene experiment (Figs. 2–4) and about 4 kΩ for the graphite–graphite experiment (Fig. 1)). For each measurement, we remove the corresponding R_c to obtain:

Limitations at small twist angles

In the experiments, the QTM is capable of controlling the twist angle of TBG continuously down to 0°. However, the measurements show that phonon emissions steps in G start to disappear when the twist angle is less than 6° and, at low bias, we measure substantial junction conductance. There are several possible reasons for this. (1) The tail of elastic tunnelling processes becomes more and more dominant and overwhelms the inelastic phonon emission processes. (2) At smaller angles near the magic angle, the two layers become strongly hybridized such that Fermi’s golden rule approximation, which is used to interpret our measurements, breaks down, as electrons can tunnel back and forth several times across the junction coherently (Fermi’s golden rule assumes just a single pass). (3) Processes involving the absorption of thermal phonons can be notable at low angles and may dominate the phonon emission processes of interest to us. The energy of a phonon at q = q_M becomes smaller with decreasing twist angle. For θ = 6°, it is already ω_ph ≈ 10 meV. For lower twist angles, it will be even smaller. Because the same EPC element is relevant for both the emission and absorption of phonons, absorption of thermally activated phonons at T = 4 K may well be responsible for the large conductive background at low angles⁴⁴. (4) Contact resistance: at lower angles, the EPC gets stronger as shown in the main text and junction conductance becomes better. The tail of the onset of phonon emission already leads to a junction resistance that can be much smaller than the contact resistance. We can account for some of this effect with our theoretical model, but if the junction resistance becomes negligibly small, then it becomes impossible to resolve its behaviour.

The EPC for the ZA (ZO′) phonon mode

In the main text, we show the extracted EPC for the TO/LO and layer-antisymmetric TA phason mode. Here we show the ZA mode (the layer-antisymmetric ZA phonon is often also called ZO′) (Extended Data Fig. 12). This mode shows a gapped dispersion in the measurement. In the angle range we measure, ZA mode phonon energy decreases as the twist angle becomes smaller. The crossing between the ZA and TA modes is around θ = 6°. The theory discussions of the ZA mode are given in ref. ⁴⁹, which shows that, similar to the TA mode, the ZA mode is dominated by the interlayer mechanism. Its coupling strength increases with decreasing angle at approximately \(\frac1{\sqrt\omega _\rmZA}\), as seen by equations (4) and (5). Extended Data Fig. 12b shows the extracted g_ZA as a function of θ, which shows the coupling strength of around 25 meV at θ = 6°. At θ = 8°, at which this mode overlaps with the TA mode in bias, we divided the total conductance step according to the ratio of the conductance steps measured at larger angle. This assumption most likely overestimates the ZA mode coupling and underestimates the TA mode coupling at θ = 8°.

Implications for superconductivity in TBG

Given the EPC we measure, we can estimate the coupling constant λ in the Bardeen–Cooper–Schrieffer (BCS) theory. The attractive electron–electron interaction from EPC is given by:

in which g_q,r is the EPC for phonon mode r and momentum q. A_unit is the unit-cell area. Our experiments measure the EPC at q_M (\(g_q_\rmM\)). For the acoustic phonons and phasons, we assume the q dependence of the EPC in the phonon mini-BZ to follow \(\sqrtq\). Because for these modes ω_q ≈ q, we get that their U_eff,r(q) is independent of momentum. For optical phonons, we assume that the EPC is independent of q. Because for these modes ω_q is independent of q, we get that, also for them, U_eff,r(q) is independent of momentum. We can thus readily integrate equation (18) over the phonon mini-BZ to obtain a U_eff,r that depends only on \(g_q_\rmM\), measured in our experiments. Consequently, the coupling constant λ is given by:

We approximate the DOS \(\nu _=\frac1WA_\rmmoire\), in which W is the bandwidth and A_moire is the moiré cell area. We then get:

For the phason mode, we measured the \(g_q_\rmM\) and \(\omega _q_\rmM\) at θ = 6°. Assuming the \(\frac1\sqrtq_\rmM\) dependence of \(g_q_\rmM\) that we measured extends to θ = 1.1°, we get:

For the optical modes, we get:

Here we further discuss the limitations of the perturbative theory approach and possible breakdowns when extending it to the magic-angle regime:

1. a.

   Hybridization of the energy bands at small angles owing to the strong interlayer coupling: this strong coupling will cause the lowest-order perturbative expansion of the interlayer coupling, resulting from the phonon field, to break down.

2. b.

   Gapping of the TA layer-antisymmetric mode owing to interlayer mechanical coupling: in our theory, we assume that the interlayer symmetric and layer-antisymmetric phonons have the same energy. We may expect, however, that, at small twist angles, the layer-antisymmetric mode will be gapped as a result of mechanical coupling between the two layers. Such a mechanical gap was measured for the layer-antisymmetric mode in Bernal bilayer graphene using Raman²⁸ and already they turned out be rather small, approximately 4 meV. We note that, in TBG, we expect the mechanical gapping of the TA mode to be much smaller—this gapping is directly related to the mechanical friction between the layers. In a bilayer, this friction is high because of the perfect commensuration between the layers. However, once the two layers are twisted, their friction reduces substantially as a result of superlubrication. Superlubrication should be relevant at angles already as small as the magic angle.

3. c.

   Relaxation near the magic angle: lattice reconstruction effects are substantial below 2°. For twist angles near or above the first magic angle of TBG, linear response theory^53,54 provides an adequate description of the relaxation effects, as confirmed by more accurate computational methods^55,56. The applicability of linear response theory indicates that the displacements associated with relaxation do not greatly affect the phonon frequencies. In terms of the EPC, reconstruction will effectively increase the areas of AB stacking and decrease those of AA stacking. This will effectively decrease the tunnelling between the layers. Because the phason EPC is dominated by the interlayer coupling, namely the modulation of the tunnelling elements between the layers, scaling down this tunnelling is expected to scale down the EPC in proportion.

4. d.

   Phonon–phonon interactions. As shown in ref. ²³, the phason mode can be damped at low energies owing to the moiré superlattice disorder. This is very important for understanding phenomena such as the strange-metal behaviour, as studied in that paper, but as long as this does not reach the strong coupling phonon–phonon limit, this will not influence the EPC estimated by our model.

Possible future experiments for measuring neutral collective modes with the QTM

The inelastic momentum-resolved scanning technique demonstrated here should be applicable to a wide range of vdW materials, including both carbon-based systems and TMDs. Beyond phonons, this method could potentially measure other neutral collective modes as well. Two fundamental requirements make these measurements particularly straightforward:

1. 1.

   Materials with small Fermi pockets, preferably away from the Γ point.

2. 2.

   The Fermi velocity is larger than the velocity of the collective mode.

Requirement 1 is well met by many carbon-based materials, for which the band structure is derived from small Dirac points. More importantly, this condition is also satisfied by the broader class of TMD materials, for which small Fermi surfaces are typically found at the K or Q points in the BZ for both electron and hole carriers.

For requirement 2, the phonon sound velocity is much smaller than the Fermi velocity in all realistic materials within these families, including both carbon-based materials and TMDs. Other neutral collective modes, such as magnons and spinons, are also expected to have much lower velocities than electrons.

Even when these conditions are not met (for example, for plasmons, which travel faster than electrons), the technique can still provide valuable insights. In such cases, however, a more thorough analysis would be required to separate the electronic contributions from those of the collective modes.

Another important point is that the experiment can be performed in two different configurations. The first is when the ‘launching’ and ‘absorbing’ layers for the collective modes are in contact, as in the current experiment. In this configuration, we can examine the collective modes of each individual layer or that of the combined twisted interface. The second configuration involves placing an insulating material, whose collective modes are to be studied, between the ‘launching’ and ‘absorbing’ layers. In fact, momentum-conserving tunnelling experiments performed in devices with a fixed twist angle have already successfully detected the phonons of a hBN barrier between the layers^41,42,43. Two examples of proposed experiments consider the possibility of studying spinons in quantum spin liquids⁴ and magnons in magnetic moiré heterostructures³.