# Phonons in magic-angle twisted bilayer graphene

Xiaoqian Liu  School of Physical Science and Technology, ShanghaiTech
University, Shanghai 200031, China  Ran Peng  School of Physical Science and
Technology, ShanghaiTech University, Shanghai 200031, China  Zhaoru Sun
<EMAIL_ADDRESS>School of Physical Science and Technology,
ShanghaiTech University, Shanghai 200031, China  Jianpeng Liu
<EMAIL_ADDRESS>School of Physical Science and Technology,
ShanghaiTech University, Shanghai 200031, China

###### Abstract

Magic-angle twisted bilayer graphene (TBG) has attracted significant interest
recently due to the discoveries of diverse correlated and topological states
in this system. Despite the extensive research on the electron-electron
interaction effects and topological properties of the electrons, the phonons
of magic-angle TBG are relatively less explored. In this work, we study the
phonon properties in magic-angle TBG based on ab initio deep potential
molecular dynamics. We have calculated phonon band structures and 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é $\Gamma$ point, we have discovered a number of soft modes which
can exhibit dipolar-like, stripe-like, and octupolar-like vibrational patterns
within the moiré supercell, as well as some “vortical” modes with nonzero curl
in real space. At the moiré $K$/$K^{\prime}$ points, there are time-reversal
breaking chiral phonon modes with nonzero local phonon polarizations. We have
further studied the phonon effects on the electronic structures by freezing
certain soft phonon modes. We find that if a soft “stripe” phonon mode at
moiré $\Gamma$ point is assumed to be frozen, the system would exhibit a
charge order which naturally explains the recent observations from scanning
tunnelling microscopy. Moreover, there are also low-frequency
$C_{2z}$-breaking modes at moiré $\Gamma$ point, which would gap out the Dirac
points at the charge neutrality point once these modes get frozen. This
provides a new perspective to the origin of correlated insulator state at the
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, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14], quantum
anomalous Hall states [9, 15, 16, 17, 18, 19, 20, 13, 21, 22, 23, 24],
unconventional superconductivity[25, 26, 27, 4, 10, 11, 12, 28, 23], and
linear-in-temperature resistivity. Around magic angle 1.05∘[29], there are two
topologically nontrivial flat bands contributed by each valley and spin
degrees of freedom [30, 31, 32, 33, 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\rm{-}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, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63,
64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75].

However, some other phenomena observed in magic-angle TBG, such as
unconventional superconductivity [25, 26, 27, 4, 10, 11, 12, 28, 23] and
linear-in-temperature resistivity [76, 77], are relatively less understood.
One of the perspectives is that these unusual transport phenomena may result
from electron-phonon couplings [78, 79, 80, 81], and are thus intimately
related to phonon vibrational properties. In particular, recent experiments
show that the correlated insulator states that are driven by $e\rm{-}e$
interactions tend to compete with superconductivity, and the latter is getting
suppressed when the $e\rm{-}e$ interaction strength increases [12, 11]. 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 [82, 83, 78, 79, 84, 80, 85, 86, 87, 81, 88],
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. To be
specific, a machine-learning-based reactive potential method [89, 90] is
adopted, which allows for an accurate multi-body description for the
interatomic potentials of large-scale systems such as magic-angle TBG (with
$\sim 11000$ 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). 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é $\Gamma$ point, we have discovered a number of soft
phonon modes with frequencies $\sim 0.05\textrm{-}0.1\,$THz, which exhibit
various intriguing vibrational patterns on the moiré length scale. At the
moiré $K$/$K^{\prime}$ 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
stripe-type 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 $C_{2z}$-breaking modes at
$\Gamma$ 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.

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
(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 THz to 2.4 THz. (e) Phonon dispersions of magic-
angle TBG from 0 to 2.4 THz.

We first apply the DPMD method to the problem of structural relaxation of
magic-angle TBG. 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)[91] , with the twist angle
$\theta$ varying from 21.79∘ to 4.41∘. In particular, for the moiré primitive
cell at each twist angle ($21.79^{\circ}\leq\theta\leq\\!4.41^{\circ}$), we
have performed ab initio MD simulations for at least 100 fs to obtain plenty
enough training data. Then the 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 $\theta\\!\approx\\!1.08^{\circ}$) with
larger moiré primitive cells with the accuracy comparable to direct DFT
simulations. The details of this method are presented in Supplementary
Information [92]. The results of the structural relaxations are presented in
Fig. 1. To be specific, in Fig. 1(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 a 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 Fig. 1(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 [93, 94, 95, 96].

We have further calculated the atomic displacements of TBG with different
twist angles. In Fig. 1(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,
$\theta\gtrapprox 4^{\circ}$, the atomic displacements in the relaxed
structure are very weak $\sim 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.

Figure 2: Low frequency optical phonon modes at $\Gamma$. (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.

Based on the fully relaxed structure, we further study the phonon properties
of magic-angle TBG. In Fig. 2(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 [83]. 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 [83, 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
$\sim\\!0\textrm{-}2.4\,$THz, as marked by the blue lines in Fig. 1(d) (the
inset shows zoom-in phonon DOS in this low-frequency regime). In Fig. 1(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 $\sim\\!0.01\textrm{-}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.

We first focus at the phonon modes at $\Gamma$. In Fig. 2 we show the
vibrational patterns of several prototypical low-frequency optical phonon
modes at $\Gamma$. To be specific, in Fig. 2(a) we show the out-of-plane
vibrational pattern of the lowest optical phonon at $\Gamma$ point with
frequency $0.043\,$THz, which exhibit maximal amplitudes in the $AA$ region.
In Fig. 2(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 $C_{2z}$ and $C_{3z}$ symmetries. This
mode is doubly degenerate with another similar dipolar mode forming a two
dimensional representation of $C_{3z}$ operation. In Fig. 2(c) we show the
out-of-plane vibration amplitudes of the two degenerate “stripe-type” modes
with frequency $0.063\,$THz. As will be discussed later, if these two modes
were frozen, the system would generate a stripe charge order that is
consistent with the STM observation reported in Ref. 6. In Fig. 2(e) and (f)
we show two “octupolar-type” out-of-plane vibrational modes with the
frequencies being $0.073\,$THz and $0.163\,$THz respectively. These two modes
break $C_{2z}$ symmetry, and may open up a gap at the CNP. In Fig. 2(g) we
show the in-plane vibration amplitudes of another optical phonon mode with
frequency $\sim 1.218\,$THz. The atomic displacement vectors show an
interesting vortex pattern winding around the $AA$ point, which exhibits non-
vanishing curl as shown in Fig. 2(h).

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 non-degenerate “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
$\eta$, which is obtained by projecting a phonon eigenvector
$\mathbf{u}_{n}(\mathbf{q})$ ($n$ the phonon mode index and $\mathbf{q}$ is
the wavevector) onto the right-handed and left-handed basis vectors, then make
the subtraction [97]: $\eta=\sum_{i}\eta_{i}=\sum_{i}(|\langle
R_{i}|\mathbf{u}_{n}(\mathbf{q})\rangle|^{2}-|\langle
L_{i}|\mathbf{u}_{n}(\mathbf{q})\rangle|^{2}))$ where
$|R_{i}\rangle\\!=\\!(1/\sqrt{2})\times[1,i]^{\textrm{T}}$ and
$|L_{i}\rangle\\!=\\!(1/\sqrt{2})\times[1,-i]^{\textrm{T}}$ are the right-
handed and left-handed basis vectors for the displacement of the $i$th atom
(“T” stands for transpose conjugation). In Fig. 3(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
$\sqrt{3}\times\sqrt{3}$ 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
$\sim\\!0.38$. In Fig. 3(b) we show the local distribution a non-degenerate
“helical” phonon mode with frequency $\sim\\!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^{\prime}$) 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^{\prime}$ points may be coupled with the electronic degrees of freedom,
and may help in the formation of orbital magnetic states breaking moiré
translational symmetry.

Figure 3: The 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\,$THz and $1.074\,$THz respectively, where the primitive
cell and the $\sqrt{3}\times\sqrt{3}$ 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 $1.214\,$THz and $1.074\,$THz
respectively. The primitive cell and doubled moiré supercell are marked with
black and dashed lines.

In Fig. 3(c) and (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\,$THz 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 anti-phase 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 Fig. 3(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.

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 Fig. 2(f) ($0.163\,$THz) under frozen mode approximation, with the
average displacement amplitudes by the mode of 0.0430 ÅṪhe inset shows the gap
at moiré $K$ and $K^{\prime}$ points as a function of average displacement
amplitudes. (c)-(d): The local charge distribution with the two soft modes
shown in Fig. 1(b) (0.046 THz) and Fig. 1(c) (0.063 THz) being frozen.

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 Fig. 4(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 while the intervalley
scattering is significantly enhanced, generating a splitting of approximately
$0.19\,$meV between the two Dirac points (see inset of Fig. 4(a)). 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 Fig. 4(b) we show the band structure in a lattice
structure with an octupolar-type phonon mode around $0.163\,$THz (Fig. 3(f))
being frozen. A gap of $\sim 0.75\,$meV is created with an average
displacement amplitude $\sim 0.05\,\mbox{\AA}$, as shown in the inset of Fig.
4(b). The coupling between the soft phonon modes and the electrons can also
significantly change the local charge distributions, which may explain
previous STM observations [6]. In particular, we have calculated the local
charge distribution contributed by the valence flat bands at each atomic site,
with the “stripe-type” phonon mode shown in Fig. 2(c) ($\sim$0.063 THz) being
frozen, presented in Fig. 4(d), after subtracting the charge distribution of
the original relaxed structure without any phonon displacement. The difference
of charge distribution forms a quasi-1D stripe pattern, 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 ($\sim\\!0.046\,$THz, shown in Fig. 2(b)) being
frozen, and the local charge distribution (also with the original charge
distribution substracted) shows a diploar pattern as shown in Fig. 4(c). These
results indicate that the phonon modes and electron-phonon couplings play a
crucial role in the correlated states of magic-angle TBG.

To summarize, in this work we have comprehensively studied the phononic
properties of magic-angle TBG based on DPMD calculations. We have discovered a
number of low-frequency optical phonon modes at the high-symmetry points of
moiré Brillouin zone, which exhibit various dipolar, stripe-like, 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. We
find that if a soft stripe-type phonon mode were frozen, the system would
exhibit a charge order which naturally explains the recent STM observation.
Moreover, the freezing of certain $C_{2z}$-broken phonon mode would open a gap
at CNP, which may provides 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.

## Acknowledgement

R. P. and J. P. L. acknowledge support from National Key R & D program of
China (grant no. 2020YFA0309601) and the National Science Foundation of China
(grant no. 12174257). Computational resources are provided by the High-
Performance Computing (HPC) Platform at ShanghaiTech University and Shanghai
HPC center. We thank the useful discussion from Kexin Zhang and Shihao Zhang.

## References

  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Nature Physics 16, 725 (2020).
  * Andrei  _et al._ [2021] E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, Nature Reviews Materials 6, 201 (2021).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Nature 556, 80 (2018a).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature 574, 653 (2019).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, _et al._ , Nature 572, 95 (2019).
  * Jiang  _et al._ [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Nature (London) 573, 91 (2019).
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature (London) 572, 101 (2019).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, _et al._ , Nature Physics , 1 (2019).
  * Serlin  _et al._ [2019] M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science (2019).
  * Stepanov  _et al._ [2020a] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Nature 583, 375 (2020a).
  * Saito  _et al._ [2020] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Nature Physics 16, 926 (2020).
  * Liu  _et al._ [2021] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. Li, Science 371, 1261 (2021).
  * Xie  _et al._ [2021a] Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, (2021a), arXiv:2107.10854 [cond-mat.mes-hall] .
  * Das  _et al._ [2021a] I. Das, C. Shen, A. Jaoui, J. Herzog-Arbeitman, A. Chew, C.-W. Cho, K. Watanabe, T. Taniguchi, B. A. Piot, B. A. Bernevig, and D. K. Efetov, (2021a), arXiv:2111.11341 [cond-mat.str-el] .
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Stepanov  _et al._ [2020b] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, (2020b), arXiv:2012.15126 [cond-mat.mes-hall] .
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 588, 610 (2020).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Nature Materials (2021).
  * Das  _et al._ [2021b] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K. Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efetov, Nature Physics 17, 710 (2021b).
  * Pierce  _et al._ [2021] A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, (2021), arXiv:2101.04123 [cond-mat.mes-hall] .
  * Shen  _et al._ [2021] C. Shen, J. Ying, L. Liu, J. Liu, N. Li, S. Wang, J. Tang, Y. Zhao, Y. Chu, K. Watanabe, T. Taniguchi, R. Yang, D. Shi, F. Qu, L. Lu, W. Yang, and G. Zhang, Chinese Physics Letters 38, 047301 (2021).
  * Choi  _et al._ [2021] Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Nature 589, 536 (2021).
  * Stepanov  _et al._ [2021] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, Phys. Rev. Lett. 127, 197701 (2021).
  * Liu and Dai [2021a] J. Liu and X. Dai, Nature Reviews Physics 3, 367 (2021a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Codecido  _et al._ [2019] E. Codecido, Q. Wang, R. Koester, S. Che, H. Tian, R. Lv, S. Tran, K. Watanabe, T. Taniguchi, F. Zhang, M. Bockrath, and C. N. Lau, Science Advances 5 (2019), 10.1126/sciadv.aaw9770.
  * Cao  _et al._ [2021] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, science 372, 264 (2021).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Phys. Rev. Lett. 123, 036401 (2019).
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. Yang, Phys. Rev. X 9, 021013 (2019).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Phys. Rev. B 99, 195455 (2019).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Liu  _et al._ [2019] J. Liu, J. Liu, and X. Dai, Phys. Rev. B 99, 155415 (2019).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8, 031089 (2018).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).
  * Kang and Vafek [2018] J. Kang and O. Vafek, Phys. Rev. X 8, 031088 (2018).
  * Zou  _et al._ [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Phys. Rev. B 98, 085435 (2018).
  * Yuan and Fu [2018] N. F. Q. Yuan and L. Fu, Phys. Rev. B 98, 045103 (2018).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Phys. Rev. Lett. 122, 246401 (2019).
  * Seo  _et al._ [2019] K. Seo, V. N. Kotov, and B. Uchoa, Phys. Rev. Lett. 122, 246402 (2019).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Phys. Rev. Lett. 124, 097601 (2020).
  * Bultinck  _et al._ [2020a] N. Bultinck, S. Chatterjee, and M. P. Zaletel, Phys. Rev. Lett. 124, 166601 (2020a).
  * Zhang  _et al._ [2019] Y.-H. Zhang, D. Mao, and T. Senthil, Phys. Rev. Research 1, 033126 (2019).
  * Bultinck  _et al._ [2020b] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Phys. Rev. X 10, 031034 (2020b).
  * Liu and Dai [2021b] J. Liu and X. Dai, Phys. Rev. B 103, 035427 (2021b).
  * Zhang  _et al._ [2020a] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Phys. Rev. B 102, 035136 (2020a).
  * Hejazi  _et al._ [2021] K. Hejazi, X. Chen, and L. Balents, Phys. Rev. Research 3, 013242 (2021).
  * Kang and Vafek [2020] J. Kang and O. Vafek, Phys. Rev. B 102, 035161 (2020).
  * Wu and Das Sarma [2020] F. Wu and S. Das Sarma, Phys. Rev. Lett. 124, 046403 (2020).
  * Chen  _et al._ [2021] B.-B. Chen, Y. D. Liao, Z. Chen, O. Vafek, J. Kang, W. Li, and Z. Y. Meng, Nature Communications 12, 5480 (2021).
  * Lu  _et al._ [2020] C. Lu, Y. Zhang, Y. Zhang, M. Zhang, C.-C. Liu, Z.-C. Gu, W.-Q. Chen, and F. Yang, arXiv preprint arXiv:2003.09513 (2020).
  * Da Liao  _et al._ [2021] Y. Da Liao, J. Kang, C. N. Breiø, X. Y. Xu, H.-Q. Wu, B. M. Andersen, R. M. Fernandes, and Z. Y. Meng, Phys. Rev. X 11, 011014 (2021).
  * He  _et al._ [2020] W.-Y. He, D. Goldhaber-Gordon, and K. T. Law, Nature Communications 11, 1650 (2020).
  * Su and Lin [2020] Y. Su and S.-Z. Lin, Phys. Rev. Lett. 125, 226401 (2020).
  * Khalaf  _et al._ [2020] E. Khalaf, N. Bultinck, A. Vishwanath, and M. P. Zaletel, “Soft modes in magic angle twisted bilayer graphene,” (2020), arXiv:2009.14827 [cond-mat.str-el] .
  * Bernevig  _et al._ [2021a] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205413 (2021a).
  * Lian  _et al._ [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Phys. Rev. B 103, 205414 (2021).
  * Bernevig  _et al._ [2021b] B. A. Bernevig, B. Lian, A. Cowsik, F. Xie, N. Regnault, and Z.-D. Song, Phys. Rev. B 103, 205415 (2021b).
  * Xie  _et al._ [2021b] F. Xie, A. Cowsik, Z.-D. Song, B. Lian, B. A. Bernevig, and N. Regnault, Phys. Rev. B 103, 205416 (2021b).
  * Soejima  _et al._ [2020] T. Soejima, D. E. Parker, N. Bultinck, J. Hauschild, and M. P. Zaletel, Phys. Rev. B 102, 205111 (2020).
  * Huang  _et al._ [2021] C. Huang, N. Wei, and A. H. MacDonald, Phys. Rev. Lett. 126, 056801 (2021).
  * Potasz  _et al._ [2021] P. Potasz, M. Xie, and A. H. MacDonald, (2021), arXiv:2102.02256 [cond-mat.str-el] .
  * Zhang  _et al._ [2021a] X. Zhang, G. Pan, Y. Zhang, J. Kang, and Z. Y. Meng, Chinese Physics Letters 38, 077305 (2021a).
  * Hofmann  _et al._ [2021] J. S. Hofmann, E. Khalaf, A. Vishwanath, E. Berg, and J. Y. Lee, (2021), arXiv:2105.12112 [cond-mat.str-el] .
  * Parker  _et al._ [2021] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Phys. Rev. Lett. 127, 027601 (2021).
  * Xie and MacDonald [2021] M. Xie and A. H. MacDonald, Phys. Rev. Lett. 127, 196401 (2021).
  * Ying  _et al._ [2021] X. Ying, M. Ye, and L. Balents, Phys. Rev. B 103, 115436 (2021).
  * Wagner  _et al._ [2021] G. Wagner, Y. H. Kwan, N. Bultinck, S. H. Simon, and S. A. Parameswaran, (2021), arXiv:2109.09749 [cond-mat.str-el] .
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, (2021), arXiv:2105.05857 [cond-mat.str-el] .
  * Zhang  _et al._ [2021b] S. Zhang, X. Lu, and J. Liu, (2021b), arXiv:2109.11441 [cond-mat.str-el] .
  * Călugăru  _et al._ [2021] D. Călugăru, N. Regnault, M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, A. Yazdani, O. Vafek, and B. A. Bernevig, (2021), arXiv:2110.15300 [cond-mat.str-el] .
  * Herzog-Arbeitman  _et al._ [2021] J. Herzog-Arbeitman, A. Chew, D. K. Efetov, and B. A. Bernevig, “Reentrant correlated insulators in twisted bilayer graphene at 25t ($2\pi$ flux),” (2021), arXiv:2111.11434 [cond-mat.str-el] .
  * Song and Bernevig [2021] Z.-D. Song and B. A. Bernevig, “Matbg as topological heavy fermion: I. exact mapping and correlated insulators,” (2021), arXiv:2111.05865 [cond-mat.str-el] .
  * Lin  _et al._ [2021] X. Lin, B.-B. Chen, W. Li, Z. Y. Meng, and T. Shi, “Exciton proliferation and fate of the topological mott insulator in a twisted bilayer graphene lattice model,” (2021), arXiv:2110.00200 [cond-mat.str-el] .
  * Polshyn  _et al._ [2019] H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Nature Physics 15, 1011 (2019).
  * Cao  _et al._ [2020] Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020).
  * Wu  _et al._ [2018] F. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. 121, 257001 (2018).
  * Lian  _et al._ [2019] B. Lian, Z. Wang, and B. A. Bernevig, Phys. Rev. Lett. 122, 257002 (2019).
  * Wu  _et al._ [2019] F. Wu, E. Hwang, and S. Das Sarma, Phys. Rev. B 99, 165112 (2019).
  * Sharma  _et al._ [2021] G. Sharma, I. Yudhistira, N. Chakraborty, D. Y. H. Ho, M. M. A. Ezzi, M. S. Fuhrer, G. Vignale, and S. Adam, Nature Communications 12, 5737 (2021).
  * Cocemasov  _et al._ [2013] A. I. Cocemasov, D. L. Nika, and A. A. Balandin, Phys. Rev. B 88, 035428 (2013).
  * Choi and Choi [2018] Y. W. Choi and H. J. Choi, Phys. Rev. B 98, 241412 (2018).
  * Eliel  _et al._ [2018] G. S. N. Eliel, M. V. O. Moutinho, A. C. Gadelha, A. Righi, L. C. Campos, H. B. Ribeiro, P.-W. Chiu, K. Watanabe, T. Taniguchi, P. Puech, M. Paillet, T. Michel, P. Venezuela, and M. A. Pimenta, Nature Communications 9, 1221 (2018).
  * Angeli  _et al._ [2019] M. Angeli, E. Tosatti, and M. Fabrizio, Phys. Rev. X 9, 041010 (2019).
  * Koshino and Son [2019] M. Koshino and Y.-W. Son, Phys. Rev. B 100, 075416 (2019).
  * Koshino and Nam [2020] M. Koshino and N. N. T. Nam, Phys. Rev. B 101, 195425 (2020).
  * Choi and Choi [2021] Y. W. Choi and H. J. Choi, Phys. Rev. Lett. 127, 167001 (2021).
  * Wang  _et al._ [2018] H. Wang, L. Zhang, J. Han, and W. E, Computer Physics Communications 228, 178 (2018).
  * Zhang  _et al._ [2020b] Y. Zhang, H. Wang, W. Chen, J. Zeng, L. Zhang, H. Wang, and W. E, Computer Physics Communications 253, 107206 (2020b).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 169 (1996).
  * [92] See Supplemental Information for: (a) details about the relaxed structure of magic-angle TBG, (b) details about the DPMD calculations, (c) more results about the optical phonon modes at high symmetry points, and (d) more results about the frozen-phonon band structures.
  * Uchida  _et al._ [2014] K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama, Phys. Rev. B 90, 155451 (2014).
  * Angeli  _et al._ [2018] M. Angeli, D. Mandelli, A. Valli, A. Amaricci, M. Capone, E. Tosatti, and M. Fabrizio, Phys. Rev. B 98, 235137 (2018).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Phys. Rev. B 96, 075311 (2017).
  * Cantele  _et al._ [2020] G. Cantele, D. Alfè, F. Conte, V. Cataudella, D. Ninno, and P. Lucignano, Phys. Rev. Research 2, 043127 (2020).
  * Zhang and Niu [2015] L. Zhang and Q. Niu, Physical Review Letters 115, 1 (2015).
  * Lucignano  _et al._ [2019] P. Lucignano, D. Alfè, V. Cataudella, D. Ninno, and G. Cantele, Phys. Rev. B 99, 195419 (2019).
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
  * Hamada [2014] I. Hamada, Phys. Rev. B 89, 121103 (2014).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comp. Phys. Comm. 271, 108171 (2022).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).

Supplementary Information of “Phonons in magic-angle twisted bilayer graphene”

I Relaxed structure of twisted bilayer graphene

Figure 5: Atomic displacement vectors of magic-angle TBG after structural
relaxation.

The atomic displacement vectors in the relaxed lattice structure of magic-
angle TBG have been presented in Fig. 5, where the color coding and the vector
fields indicate the amplitudes and directions of the atomic displacements.
Inspection shows that after structural relaxation the interlayer distance in
the $AA$ region tends to be increased, leading to the out-of-plane
corrugation. In particular, the interlayer distance at the $AA$ point is as
large as 3.62 Å, while at the $AB/BA$ point it is 3.36 Å, consistent with
previous report [93, 98]. This is clearly shown in Fig. 6, where the color
coding denotes the in-plane variation of the interlayer distance in the moiré
supercell. Moreover, the $AA$ ($AB/BA$) region is contracted (expanded) after
structural relaxation. The in-plane atomic displacement vectors wind around
the $AA$ point, forming a vortex-like pattern.

Figure 6: The interlayer distance as a function of real-space coordinates
$(x,y)$ in the moiré superlattice.

II Deep potential molecular dynamics

Figure 7: Validation of the neural network (NN) potential: (a) Comparison of
total energies calculated from NN potential and DFT.deep-potential model
prediction; (b-d) Comparison of forces calculated from NN potential and DFT in
$x,y,z$ directions, respectively.

Primitive moiré supercells with different twist angles (from 21.79°to 4.41°)
have been constructed to get complex enough local atomic environments [99].
The rotation center is chosen to be the center of the hexagon formed by the
carbon atoms. With such a choice of rotation center, the $D_{6}$ symmetry of
graphene is preserved for the TBG system. In order to describe the tiny
lattice distortions of TBG during structural relaxations, we constructed
several structures with small perturbative atomic displacements from ideal
moiré pattern as the initial configurations for ab initio molecular dynamics
(AIMD). For each configurations, we performed AIMD calculation for at least
100 fs to obtain plenty enough training data. The the AIMD simulation have
been carried out using the Vienna Ab initio Simulation Package (VASP) [91]
with vdW-DF2 exchange-correlation functional [100]. $D_{6}$ symmetry is
enforced in the AIMD simulation. Under periodic conditions, a vacuum layer of
20 Å in the $z$ direction has been set up. Then the reliable interatomic
potential based on neural networks can be trained from DFT energies and
forces, based on which we further studied the structural and phononic
properties of TBG at smaller twist angles (including the magic angle
$\theta\\!\approx\\!1.08^{\circ}$) with larger moiré primitive cells with the
accuracy comparable to DFT calculations. Otherwise it would be extremely
demanding to directly tackle with the structural properties of magic-angle TBG
with DFT accuracy. The structural relaxation and force-constant matrix
calculations were calculated by Large-scale Atomic/Molecular Massively
Parallel Simulator (LAMMPS) [101] with such NN potential[89]. Finally, based
on the force-constant matrix, phonons were calculated by phonopy [102] package
with a $2\times 2$ moiré supercell. The accuracy of our NN potential has been
validated by comparing the total energies and forces of various structures
with DFT calculations, as shown in Fig. 7.

III Soft phonon modes and their symmetry characters

Figure 8: The first 30 soft phonon modes at $\Gamma$

The real-space out-of-plane vibrational patterns of the first 30 soft optical
phonon modes of magic-angle TBG at the high symmetry points $\Gamma$ and $M$
have been shown in Fig. 8 and Fig. 9 respectively. We further present the
irreducible representations of the first 300 optical phonon modes at $\Gamma$,
$M$, and $K$ in 1.

Figure 9: The first 30 soft phonon modes at $M$ Table 1: The first 300 phonon modes in $\Gamma$, $M$, and $K$. Such representation is the point group in Schoenflies notation. The first 30 phonon modes have been shown in here, and the whole table could be seen in https://github.com/Simon-lxq/Phonon-symmetry eigen No. |  $\Gamma$ |  M |  K
---|---|---|---
1 |  $D_{6}$ |  $C_{2}$ |  $C_{6}$
2 |  $D_{6}$ |  $C_{2}$ |  $D_{6}$
3 |  $D_{6}$ |  $D_{2}$ |  $D_{6}$
4 |  $C_{6}$ |  $C_{2}$ |  $C_{2}$
5 |  $C_{2}$ |  $C_{2}$ |  $D_{6}$
6 |  $C_{2}$ |  $C_{2}$ |  $C_{6}$
7 |  $D_{2}$ |  $C_{2}$ |  $D_{6}$
8 |  $D_{2}$ |  $D_{2}$ |  $C_{2}$
9 |  $D_{3}$ |  $C_{2}$ |  $C_{6}$
10 |  $C_{2}$ |  $C_{2}$ |  $D_{2}$
11 |  $C_{2}$ |  $D_{2}$ |  $C_{2}$
12 |  $C_{6}$ |  $C_{2}$ |  $D_{3}$
13 |  $D_{2}$ |  $C_{2}$ |  $C_{6}$
14 |  $D_{2}$ |  $C_{2}$ |  $C_{2}$
15 |  $D_{3}$ |  $C_{2}$ |  $C_{2}$
16 |  $D_{2}$ |  $C_{2}$ |  $C_{2}$
17 |  $D_{2}$ |  $D_{2}$ |  $C_{1}$
18 |  $D_{3}$ |  $C_{2}$ |  $D_{6}$
19 |  $C_{2}$ |  $C_{2}$ |  $C_{6}$
20 |  $C_{2}$ |  $C_{2}$ |  $D_{6}$
21 |  $C_{6}$ |  $C_{2}$ |  $D_{6}$
22 |  $C_{6}$ |  $D_{2}$ |  $D_{6}$
23 |  $C_{2}$ |  $C_{2}$ |  $D_{3}$
24 |  $C_{2}$ |  $C_{2}$ |  $C_{2}$
25 |  $D_{2}$ |  $D_{2}$ |  $D_{3}$
26 |  $D_{2}$ |  $C_{2}$ |  $C_{1}$
27 |  $D_{3}$ |  $C_{2}$ |  $C_{6}$
28 |  $D_{2}$ |  $C_{2}$ |  $C_{2}$
29 |  $D_{2}$ |  $D_{2}$ |  $C_{2}$
30 |  $C_{2}$ |  $C_{2}$ |  $C_{6}$

IV Band structures under frozen phonon approximation

The energy bands with some typical soft phonon modes being brozen are shown in
Fig. 10. In particular, in Fig. 10(a)-(c) we show the band structures with
out-of-plane “octupolar” phonons assumed to be frozen, while in Fig. 10(d)-(f)
we show the band structures with some “vortical” phonon modes being frozen
(the phonon modes are shown in the insets). Our calculations indicate that the
octupolar phonon modes tend to open up a gap at the charge neutraility point
due to the broken $C_{2z}$ symmetry, while the “vortical” phonon modes tend to
enhance intervalley scattering, which splits two Dirac points.

Figure 10: The energy bands of typical phonon modes in high symmetry point
$\Gamma$. (a)-(c) The band structure of the soft modes with the frequencies
$0.0726\,$THz,$0.1634\,$THz,$0.2387\,$THz, respectively. The insets are the
general patterns of out-of-plane vibrational amplitudes, where the average
amplitudes is 0.05 Å(̇d)-(f) The band structure of the soft modes with the
frequencies $1.2184\,$THz,$2.1259\,$THz,$3.6348\,$THz, respectively. The
insets are the general patterns of in-plane vibrational amplitudes.

V Lattice vectors and reciprocal lattice vectors of a commensurate supercell

Figure 11: Real space base vector and inverted space base vector

We consider studying a phonon mode at moiré wavevector $\mathbf{q}$, which can
be written as

$\mathbf{q}=\frac{q_{1}}{p_{1}}\mathbf{g_{1}}+\frac{q_{2}}{p_{2}}\mathbf{g_{2}}$
(1)

where $(q_{1},p_{1})$ and $(q_{2},p_{2})$ are two pairs of co-prime integers.
In order to plot the real-space vibrational pattern of the phonon mode at
wavevector $\mathbf{q}$, we need to construct a minimal real-space supercell
with the superlattice vectors:

$\begin{split}\mathbf{R_{1}}=n_{11}\mathbf{t_{1}}+n_{12}\mathbf{t_{2}}\\\
\mathbf{R_{2}}=n_{21}\mathbf{t_{1}}+n_{22}\mathbf{t_{2}}\end{split}$ (2)

We require
$\mathbf{R_{1}}\cdot\mathbf{q}=2\pi,\mathbf{R_{1}}\cdot\mathbf{q}=0$, such
that

$\begin{split}\frac{n_{11}q_{1}}{p_{1}}+\frac{n_{12}q_{2}}{p_{2}}=1\\\
\frac{n_{21}q_{1}}{p_{1}}+\frac{n_{22}q_{2}}{p_{2}}=0\end{split}$ (3)

And from Eq.(3) we obtain:

$\begin{split}\mathbf{G_{1}}=\frac{n_{22}\mathbf{g_{1}}-n_{21}\mathbf{g_{2}}}{n_{11}n_{22}-n_{12}n_{21}}\\\
\mathbf{G_{2}}=\frac{n_{11}\mathbf{g_{2}}-n_{12}\mathbf{g_{1}}}{n_{11}n_{22}-n_{12}n_{21}}\end{split}$
(4)

We find integer solutions of Eq.(3), and construct a supercell with minimal
area. The area of the primitive supercell $\Omega_{S}$ is expressed as:

$\Omega_{S}=(\mathbf{R_{1}}\times\mathbf{R_{2}})\cdot\hat{z}=(n_{11}n_{22}-n_{12}n_{21})\Omega_{0}$
(5)