Moiré Phonons in Magic-Angle Twisted Bilayer Graphene.

Magic-angle twisted bilayer graphene (TBG) has attracted significant interest recently due to the discoveries of diverse correlated and topological states. In this work, we study the phonon properties in magic-angle TBG based on many-body classical potential and interatomic forces generated by a deep neural network trained with data from ab initio calculations. We have discovered a number of soft modes which can exhibit dipolar, quadrupolar, and octupolar vibrational patterns in real space, as well as some time-reversal breaking chiral phonon modes. We have further studied the phonon effects on the electronic structures by freezing certain soft phonon modes. We find that if a soft quadrupolar phonon mode is assumed to be frozen, the system would exhibit a charge order which is perfectly consistent with recent experiments. Moreover, once some low-frequency C2z-breaking modes get frozen, the Dirac points at the charge neutrality point would be gapped out, which provides an alternative perspective to the origin of correlated insulator state at charge neutrality point.

Twisted bilayer graphene (TBG) system around the magic angle is an ideal platform to realize various intriguing quantum phases[1,2] such as the correlated insulators,[3−14] quantum anomalous Hall states, and[9,13,15−24] unconventional superconductivity.[4,10−12,23,25−28] Around magic angle 1.05°,[29] there are two topologically nontrivial flat bands contributed by each valley and spin degrees of freedom.[30−34] A lot of the unusual phenomena, including correlated insulators and quantum anomalous Hall effects, can be attributed to the presence of such topologically nontrivial flat bands in the electronic degrees of freedom. The electron–electron (e–e) Coulomb interactions dominates over the kinetic energy near magic angle, and the interplay between the strong Coulomb correlations and the nontrivial topology of the flat bands give rise to diverse correlated insulator states and topological states, which have been extensively studied from the theoretical point of view over the past few years.[35−74]

However, some other phenomena observed in magic-angle TBG, such as unconventional superconductivity[4,10−12,23,25−28] and linear-in-temperature resistivity,[75,76] are relatively less understood. One of the perspectives is that these unusual transport phenomena may result from electron–phonon couplings,[77−80] and are thus intimately related to phonon vibrational properties. In particular, recent experiments show that the correlated insulator states that are driven by e–e interactions tend to compete with superconductivity, and the latter is getting suppressed when the e–e interaction strength increases.[11,12] This seems to imply that electron–phonon interactions are crucial in order to understand the origin of superconductivity in this system. On the other hand, despite a few pioneering works[77−87] the phononic properties and electron–phonon couplings in magic-angle TBG are much less explored compared to the comprehensive research on the electronic degrees of freedom.

In such a context, in this work we study the phonon properties of magic-angle TBG based on ab initio deep potential molecular dynamics (DPMD) method.[88,89] To be specific, a machine-learning-based reactive potential method[88,89] is adopted, which allows for an accurate multibody description for the interatomic potentials of large-scale systems such as magic-angle TBG (with ∼11 000 atoms in each moiré primitive cell). The accuracy of the calculated total energies and forces based on such DPMD method is comparable to that from first-principles calculations based on density functional theory (DFT), which has been verified in various previous studies.[90−98] Using this method, we have calculated phonon band structures and phonon density of states at the magic angle, and have systematically analyzed the phonon eigenmodes at high-symmetry points in the moiré Brillouin zone. In particular, at the moiré Γ point, we have discovered a number of soft phonon modes with frequencies ∼0.05–0.1 THz, which exhibit various intriguing vibrational patterns on the moiré length scale. At the moiré K/K′ points, there are time-reversal breaking chiral phonon modes with nonzero local phonon polarizations, which may be coupled with the orbital motions of electrons and boost the formation of orbital magnetic ground states. We have further studied the phonon effects on the electronic structures by freezing certain soft phonon modes. We find that if a soft quadrupolar phonon mode were frozen, the system would exhibit a charge order which naturally explains the recent observations from scanning tunnelling microscopy (STM).[6] Moreover, if some of the soft C2-breaking modes at Γ point were frozen, the electronic band structure becomes gapped at the charge neutrality point (CNP), which may provide a new perspective to the origin of correlated insulator state at CNP observed in experiments.

We first apply the DPMD method to the problem of structural relaxation of magic-angle TBG.[99] In order to obtain a reliable interatomic potential, we first performed first-principles calculations for TBG based on DFT as implemented by the Vienna Ab initio Simulation Package (VASP),[100] with the twist angle θ varying from 21.79° to 4.41°. Then, a reliable interatomic potential based on neural networks can be trained from DFT energies and forces, based on which we further study the structural and phononic properties of TBG at smaller twist angles (including the magic angle θ ≈ 1.08°) with larger moiré primitive cells with the accuracy comparable to direct DFT simulations. The details of this method are presented in Supporting Information. The results of the structural relaxations are presented in Figure . To be specific, in Figure a we show an overview of the relaxed structure of magic-angle TBG, where the color coding indicates the amplitudes of the atomic displacements deviating from the ideal moiré lattice positions obtained by a simple twist of the two graphene monolayers centered at the AA point. We see that the atomic displacements have largest amplitudes in the region between the AA and AB/BA points, forming a ring encircling the AA region. Further analysis reveal that these atomic displacement vectors of the relaxed structure are “winding” around the AA point, generating a vortex-like in-plane vector field, as clearly shown by the arrows (representing the directions of the in-plane displacements) in Figure b. Moreover, as AB/BA stacked bilayer graphene is energetically more stable than the AA stacked one, carbon atoms near the AA point tend to misalign with each other, resulting in a squeeze of the AA region. In addition to the in-plane displacements, our calculations indicate that there are also out-of-plane corrugations, i.e., the periodic modulations of interlayer distance in different regions of the moiré pattern, with the distance at the AA (AB/BA) point being 3.62 Å (3.36 Å), which is consistent with previous results.[101−107]

Figure 1

Structures of TBG after geometric optimization. (a) An overview of the relaxed structure and (b) a zoomed in view of the relaxed structure with the color coding and arrows denoting the amplitudes and directions of the in-plane atomic displacement vectors. (c) Amplitudes of atomic displacement vectors plotted along a line connecting AA-bridge-AB points (marked in panel b), where the twist angle is at 7.43°, 4.41°, 3.15°, 1.54°, and 1.08°. (d) The phonon density of states (DOS) of magic-angle TBG, where the inset shows the low-frequency DOS from 0 to 2.4 THz. (e) Phonon dispersions of magic-angle TBG from 0 to 2.4 THz.

We have further calculated the atomic displacements of TBG with different twist angles. In Figure c, we show the amplitudes of the atomic displacements in the relaxed TBG structures at different twist angles along a line connecting the AA-bridge-AB points. At large twist angles, θ > 4°, the atomic displacements in the relaxed structure are very weak ∼0.005 Å; but for smaller angles, the displacement amplitudes are greatly enhanced. Around the magic angle, the maximal displacement reaches 0.12 Å. Such significant structural relaxations are crucial in determining the low-energy electronic structures of magic-angle TBG.

Based on the fully relaxed structure, we further study the phonon properties of magic-angle TBG. In Figure a, we show the calculated phonon density of states for magic-angle TBG which exhibit a few peaks around 11.35, 19.36, 40.68, and 47.88 THz, consistent with previous reports.[82] Despite slight shift of DOS peaks, the phonon DOS of magic-angle TBG at the first glance are quite similar to those of bilayer graphene.[81,82] However, if one closely looks at the low-frequency phonon modes with long-wavelength vibrational patterns, it would be completely different. In magic-angle TBG, one would obtain numerous low-frequency moiré phonon modes which cannot be simply interpreted as folding atomic phonon modes of bilayer graphene into the moiré Brillouin zone. In this work, we specifically focus at the low-frequency, long-wavelength moiré phonon properties with frequency ∼0–2.4 THz, as marked by the blue lines in Figure d (the inset shows zoom-in phonon DOS in this low-frequency regime). In Figure e, we further show the phonon dispersions of magic-angle TBG with the frequency ranging from 0 to 2.4 THz. There are already hundreds of phonon modes within such a small frequency regime, and some of the optical phonon modes are extremely soft ∼0.01–0.1 Thz, implying that the system is likely to undergo structural transitions, which may significantly change the electronic structure. In what follows we will comprehensively analyze the low-frequency optical phonon modes of magic-angle TBG at the high-symmetry points in the moiré Brillouin zone.

Figure 2

Low frequency optical phonon modes at Γ. (a–f) The out-of-plane (normalized) vibrational amplitudes for the soft phonon modes with frequencies at 0.043, 0.046, 0.063, 0.063, 0.073, and 0.163 THz, respectively, where the black hexagon marks the moiré primitive cell. (g) The inplane vortex pattern of vortical mode with frequency 1.218 THz. (h) The curl field of the vortical mode.

We first focus at the phonon modes at Γ. In Figure we show the vibrational patterns of several prototypical low-frequency optical phonon modes at Γ. To be specific, in Figure a we show the out-of-plane vibrational pattern of the lowest optical phonon at Γ point with frequency 0.043 THz, which exhibit maximal amplitudes in the AA region. In Figure b, we show the out-of-plane vibrational amplitudes of the second optical phonon mode with frequency 0.046 THz, which exhibits a dipolar vibration pattern that clearly breaks C2 and C3 symmetries. This mode is doubly degenerate with another similar dipolar mode forming a two-dimensional representation of C3 operation. In Figure c, we show the out-of-plane vibration amplitudes of the two degenerate quadrupolar-type modes with frequency 0.063 THz. As will be discussed later, if these two modes were frozen, the system would generate a quadrupolar charge order that is consistent with the STM observation reported in ref (6). In Figure e,f we show two octupolar-type out-of-plane vibrational modes with the frequencies being 0.073 and 0.163 THz, respectively. These two modes break C2 symmetry, which may help opening up a gap at the CNP if these modes are coupled with electrons. In Figure g, we show the in-plane vibration amplitudes of another optical phonon mode with frequency ∼1.218 THz. The atomic displacement vectors show an interesting vortex pattern winding around the AA point, which exhibits nonvanishing curl as shown in Figure h. Further analysis reveal that these low-frequency moire phonon modes originate from the linear superpositions of the flexural modes of the two graphene monolayers. However, the interlayer couplings also play an important role, which would couple the flexural modes from the two layers at different moire reciprocal vectors, such that the moire phonon patterns of TBG are significantly reconstructed compared to those of two decoupled graphene monolayers. We refer the readers to Supporting Information for more details.

We continue to discuss the phonon properties at other high symmetry points in the moiré Brillouin zone. Interestingly, at moiré K point we find that there are doubly degenerate chiral phonon modes of opposite chiralities, as well as some nondegenerate “helical” phonon modes which have vanishing net chirality but nontrivial local distribution of phonon polarizations. To be specific, the chirality of a phonon mode can be characterized by its polarization denoted by η, which is obtained by projecting a phonon eigenvector u(q) (n the phonon mode index and q is the wavevector) onto the right-handed and left-handed basis vectors, then make the subtraction:[105] η = ∑η = ∑(|⟨R|u(q)⟩|2 – |⟨L|u(q)⟩|2)), where and are the right-handed and left-handed basis vectors for the displacement of the ith atom (“T” stands for transpose conjugation). In Figure a we show the calculated local distribution of phonon polarization of one of the doubly degenerate chiral modes at K with frequency 1.214 THz, where the dashed line marks a moiré supercell, and the solid black hexagon marks the moiré primitive cell. We see the phonon polarization reaches the maximal value in the AA region, and the net polarization of each moiré supercell ∼0.38. In Figure b we show the local distribution a nondegenerate “helical” phonon mode with frequency ∼1.074 THz, which has an octupolar-type distribution of local polarization with vanishing net chirality. It is worthwhile to note that although the atomic displacement vectors of a phonon mode at K (or K′) point breaks the primitive moiré translational symmetry, the phonon polarization still preserves the original translational symmetry as can be seen from the definition of phonon polarization. These intriguing time-reversal breaking phonon modes at K/K′ points may be coupled with the electronic degrees of freedom, and may help in the formation of orbital magnetic states breaking moiré translational symmetry. It is worthwhile to note that these moiré chiral phonon modes are qualitatively different from the chiral optical phonons at K and K′ points of monolayer graphene or Bernal bilayer graphene: the phonon polarization of the latter varies on atomic length scale with much higher vibrational frequency ≳16.4 THz, while the low-frequency moiré chiral phonons in TBG are speculated to originate from the complex combinations of the long-wavelength in-plane acoustic phonon modes of the two graphene monolayers.

Figure 3

Phonon polarizations at K point and soft modes at M point. (a,b) The local phonon polarizations of the optical modes at K point with frequencies 1.214 and 1.074 THz respectively, where the primitive cell and the moiré supercell are marked with solid line and dashed lines. (c,d) The out-of-plane vibrational amplitudes of the soft modes at M point with frequencies at 0.021 and 0.082 THz, respectively. The primitive cell and doubled moiré supercell are marked with black and dashed lines.

In Figure c,d, we show the real-space distributions of the out-of-plane vibrational amplitudes of the first and seventh optical phonon modes with extremely low frequencies 0.021 and 0.082 THz, respectively. The out-of-plane atomic displacement vectors of the first phonon mode at M point form a stripe pattern with antiphase modulation along a direction perpendicular to the lattice vector. Such an extremely soft phonon mode may be relevant with the unusual zero-Chern-number insulator states observed at filling factors 1 and 3 of the flat bands.[10] Once this soft mode is frozen, the system has to break moiré translational symmetry in such a way to double the moiré primitive cell as indicated by the black dashed rectangle, and within the doubled moiré supercell there would be two holes (six holes) at filling factor 1 (3) for each doubled supercell, which has the chance to realize a zero-Chern-number insulator state without the necessity to completely lift the valley-spin degeneracy of the topologically nontrivial flat bands. In Figure d, we show the real-space distribution of the out-of-plane vibrational amplitudes of the third phonon mode at M point, which exhibit an interesting quadrupolar pattern.

Before discussing the effects of the phonon modes, the structural relaxation effects on the electronic band structure of magic-angle TBG should be first understood. In Figure a we show the band structures of magic-angle TBG before and after relaxation as marked by the dashed and solid lines, respectively. Clearly the lattice relaxation enhances both the gaps between remote bands and flat bands and the bandwidth of the flat bands. In the meantime, the intervalley scattering is significantly enhanced, generating a splitting of approximately 0.19 meV between the two Dirac points (see inset of Figure a).

Figure 4

(a) Band structure of relaxed (solid lines) and unrelaxed (dashed lines) magic-angle TBG, where the inset indicates the gap opening at the Dirac point. (b) Flat bands of magic-angle TBG with the octupolar-type phonon modes shown in Figure f (0.163 THz) under frozen mode approximation, with the average displacement amplitudes by the mode of 0.0430 Å. The inset shows the gap at moiré K and K′ points as a function of average displacement amplitudes. (c,d) The local charge distribution with the two soft modes shown in Figure b (0.046 THz) and Figure c (0.063 THz) being frozen.

We continue to study the electron–phonon coupling effects on the electronic degrees of freedom by calculating the electronic band structures under the “frozen-phonon” approximation. In Figure b we show the band structure in a lattice structure with an octupolar-type phonon mode around 0.163 THz (Figure f) being frozen. A gap of ∼0.75 meV is created at an average displacement amplitude ∼0.05 Å, as shown in the inset of Figure b. Moreover, recently Jiang et al. reported that when the flat bands of magic-angle TBG are partially filled, a “pseudo-gap” phase accompanied by a global stripe charge order breaking the 3-fold rotational symmetry shows up.[6] Here we show that the coupling between the soft phonon modes and the electrons can also significantly change the local charge distributions, which can explain the origin of the peculiar stripe charge order observed in previous STM measurements.[6] In particular, we have calculated the local charge distribution contributed by the valence flat bands at each atomic site, assuming the quadrupolar phonon mode shown in Figure c (∼0.063 THz) is frozen, which is presented in Figure d after subtracting the charge distribution of the original relaxed structure without any phonon displacement. The difference of the charge distribution forms a quadrupolar order, which is in perfect agreement with previous STM observation.[6] Similarly, we have also calculated the local charge distribution contributed by the valence flat bands with the “dipolar” phonon mode (∼0.046 THz, shown in Figure b) being frozen, and the local charge distribution (also with the original charge distribution subtracted) shows a diploar pattern as shown in Figure c. These results indicate that the phonon modes and electron–phonon couplings play a crucial role in understanding the charge ordering and correlated states of magic-angle TBG. To evaluate the electron–phonon coupling in a more quantitative manner, we have explicitly calculated the electron–phonon coupling matrix elements within the flat-band subspace of magic-angle TBG for 12 low-frequency moiré phonon modes at Γ point, from which the effective electron–phonon coupling constant associated with a given phonon mode can be extracted. The details of the electron–phonon coupling calculations are provided in Supporting Information.

To summarize, in this work we have comprehensively studied the phononic properties and electron–phonon couplings of magic-angle TBG based on DPMD calculations and atomistic tight-binding modeling. We have discovered a number of low-frequency optical phonon modes which exhibit various dipolar, quadrupolar, octupolar, vortical, and chiral vibrational patterns on the moiré length scale. We have further studied the electron–phonon couplings by freezing certain types of soft phonon modes and study their effects on electronic band structures based on atomistic tight-binding models. We find that if a soft stripe-type phonon mode is frozen, the system would exhibit a charge order which naturally explains the recent STM observation. Moreover, the freezing of certain C2-broken phonon mode would open a gap at CNP, which may provide a new perspective to the origin of correlated insulator at CNP observed in experiment. Our work is a significant step forward in understanding the phononic and electronic properties of magic-angle TBG, and will provide useful guidelines for future experimental and theoretical studies.