# Low-energy moiré phonons in twisted bilayer van der Waals heterostructures

Jonathan Z. Lu<EMAIL_ADDRESS>Department of Physics, Harvard
University, Cambridge, Massachusetts 02138, USA  Ziyan Zhu
<EMAIL_ADDRESS>Department of Physics, Harvard University, Cambridge,
Massachusetts 02138, USA  Mattia Angeli  John A. Paulson School of Engineering
and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA
Daniel T. Larson  Department of Physics, Harvard University, Cambridge,
Massachusetts 02138, USA  Efthimios Kaxiras  Department of Physics, Harvard
University, Cambridge, Massachusetts 02138, USA  John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts
02138, USA

###### Abstract

We develop a low-energy continuum model for phonons in twisted moiré bilayers,
based on a configuration-space approach. In this approach, interatomic force
constants are obtained from density functional theory (DFT) calculations of
untwisted bilayers with various in-plane shifts. This allows for efficient
computation of phonon properties for any small twist angle, while maintaining
DFT-level accuracy. Based on this framework, we show how the low-energy phonon
modes, including interlayer shearing and layer-breathing modes, vary with the
twist angle. As the twist angle decreases, the frequencies of the low-energy
modes are reordered and the atomic displacement fields corresponding to phonon
eigenmodes break translational symmetry, developing periodicity on the moiré
length scale. We demonstrate the capabilities of our model by calculating the
phonon properties of three specific structures: bilayer graphene, bilayer
molybdenum disulfide (MoS2), and molybdenum diselenide-tungsten diselenide
(MoSe2-WSe2).

††preprint: APS/123-QED

##  I Introduction

Two-dimensional (2D) van der Waals (vdW) multi-layer structures can be
constructed by stacking 2D materials in a layer-by-layer fashion. If the
stacking involves a twist of one layer by an angle $\theta$ with respect to
the next layer, the relative misorientation between atomic lattices produces a
moiré interference pattern with a length scale larger than the lattice
constants of the individual layers. Such moiré bilayers can give rise to
localized electronic states and tunable electronic properties [1, 2, 3, 4, 5,
6, 7]. In twisted bilayer graphene, superconductivity and other unconventional
correlated states have been observed when the twist angle $\theta$ takes
values near the “magic angle” [8, 9, 10, 11]. Similar behavior has been shown
in twisted bilayers composed of transition metal dichalcogenides (TMDCs),
albeit often with decreased sensitivity to the twist angle [12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24]. Determination of the role of phonons is a
crucial step towards an understanding of these strongly correlated electronic
states, particularly superconductivity [25, 26, 27, 28]. The collective
vibration of a twisted bilayer can be strongly affected by the moiré pattern,
forming the so-called moiré phonons [29, 28]. Moiré phonons are also of
general experimental interest, and have recently been observed in graphene
[30]. The variation of phonon properties with the twist angle allow Raman
spectroscopy to be used to help characterize moiré bilayers [31, 32, 33, 34,
35].

The interplay of two distinct length scales and the large number of atoms
within a moiré cell make calculations of the phonon properties extremely
challenging. Due to the rapid scaling in computational cost with respect to
the moiré length, direct calculations of phonons using density functional
theory (DFT) are currently not feasible for twist angles less than a few
degrees. Moreover, DFT calculations using periodic supercells can only be
carried out for a specific set of commensurate angles. Instead, a continuum
approximation has been successfully used to model the electronic structure and
atomic relaxation of twisted bilayers and trilayers [36, 37, 38, 39, 40, 41,
42]. Continuum models have also been considered for moiré phonon analysis. The
most efficient of these are empirical continuum models, but they leave out
important details such as the out-of-plane degree of freedom and phonon modes
that require coupling, for instance the shearing and layer-breathing modes
[29].

We propose a first-principles-based continuum model of moiré phonons, similar
to the one developed by Quan  _et al._ [35] to help interpret the Raman
spectra of twisted MoS2 bilayers, but we have extended it to momenta beyond
the $\Gamma$-point and have applied the model to several different bilayer
systems. To bypass the need to compute the interatomic forces for all atoms in
a large moiré supercell, we adopt a configuration space formalism to describe
the local atomic environment. This approach has been used as the basis for a
continuum model of atomic relaxations in twisted vdW heterostructures [43, 44,
45, 46, 41]. Using this model we obtain the low-energy phonon properties of
bilayer graphene, bilayer MoS2, and a MoSe2-WSe2 heterostructure, and show how
a relative twist between the layers modifies the phonon frequencies and
eigenmodes.

The paper is organized as follows. In Section II we derive the moiré dynamical
matrix in configuration space, including twist-angle-dependent relaxation
[41]. We then compare our model predictions to the results of three other
approaches in Section III. In Section IV we apply our model to three
representative moiré bilayers and analyze their band structure and real space
atomic displacements. Finally, we summarize our results and their implications
in Section V.

##  II Configuration Space Continuum Model

The workflow of our configuration space continuum (CSC) model is summarized in
Fig. 1. We first construct the moiré dynamical matrix by computation of the
dynamical matrices of many rigidly shifted bilayers—one for each configuration
of atoms in layer 1 with respect to layer 2—and the individual monolayers.
From the monolayer we obtain the optimized lattice constant and the elastic
constants that are used for relaxation. We relax the ionic positions of each
shifted bilayer as well, keeping the horizontal positions of one atom per
layer fixed to maintain the registry between the layers. We then compute
atomic forces using the frozen phonon approach, which under the CSC formalism
produces the unrelaxed (rigid) moiré dynamical matrix. Finally, we include
twist-dependent relaxation which yields the relaxed moiré form of the
dynamical matrix. We provide a derivation of the model in this section, and
details of the calculations can be found in Appendices A and B. Throughout the
derivation, we use tildes to indicate a quantity in the moiré supercell, and
when referring to the dynamical matrix we use a bar to denote the matrix in
real space (the absence of the bar means the matrix is expressed in Fourier
space).

Figure 1: Schematic representation of the numerical implementation of the CSC
model. A collection of monolayer (blue) and bilayer (pink) calculations are
combined into a rigid (unrelaxed) moiré dynamical matrix and used to determine
the bilayer relaxation pattern (both purple). All of these elements are
combined to form the relaxed moiré dynamical matrix. The square labeled
“relaxed moiré dynamical matrix” formally represents the dynamical matrix
construction: the diagonal blocks are from the monolayer calculations whereas
the off-diagonal blocks are constructed from the stacking dependent force
field calculations.

###  II.1 Moiré Dynamical Matrix

The phonon spectrum of a crystal can be obtained by diagonalizing the
dynamical matrix with the basis expanded in momentum space. The moiré
dynamical matrix at a point $\widetilde{\bm{k}}$ in the moiré Brillouin zone,
which we call the center site, is represented generically as a layer-block
matrix encoding intralayer (diagonal) and interlayer (off-diagonal) couplings:

$\displaystyle
D_{\text{moir\'{e}}}(\bm{\widetilde{k}})=\begin{bmatrix}D^{1}(\bm{\widetilde{k}})&D^{12}\\\
D^{12*}&D^{2}(\bm{\widetilde{k}}).\end{bmatrix},$ (1)

Each block of the moiré dynamical matrix depends on two reciprocal-space
degrees of freedom, $\bm{k}^{(\ell)}$ and $\bm{k}^{(\ell^{\prime})}$, which
are momenta from the monolayer Brillouin zones of layers $\ell$ and
$\ell^{\prime}$. This is because in a twisted system, the interatomic forces
depend not only on the pair-wise distance of two atoms, but also the position
of the atom due to the change of the local environment.

###  II.2 Moiré Scattering Selection Rule

We first derive how the degrees of freedom in the basis of Eq. (1) are
connected to each other through a selection rule. We begin with a monolayer,
which is the intralayer term that has no stacking dependence. Considering only
one layer $\ell$, the monolayer phonon equation of motion can be expressed as

$\displaystyle\sum_{j\nu\beta}\overline{D}_{\mu\nu\alpha\beta}^{\ell}(\bm{R}_{i}-\bm{R}_{j})\delta\widehat{u}_{j\ell\nu\beta}=\omega^{2}\delta\widehat{u}_{i\ell\mu\alpha},$
(2)

where $\overline{D}_{\mu\nu\alpha\beta}^{\ell}(\bm{R}_{i}-\bm{R}_{j})$ is the
dynamical matrix element, $\bm{R}_{i}$ is a monolayer lattice vector,
$\mu,\nu$ are Cartesian degrees of freedom, $\alpha,\beta$ are sublattice
degrees of freedom, and $\delta\widehat{u}_{i\ell\mu\alpha}$ is the phonon
displacement component that corresponds to atom $\alpha$ in direction $\mu$.
We can define a phonon displacement in the Bloch basis as follows,

$\displaystyle\begin{aligned}
\delta\widehat{u}_{i\ell\mu\alpha}(\bm{k}^{(\ell)})&=\frac{1}{\sqrt{|\Gamma^{*}|}}e^{i\bm{k}^{(\ell)}\cdot\bm{R}_{i}}\delta
u_{\ell\mu\alpha}(\bm{k}^{(\ell)}),\end{aligned}$ (3)

where $\bm{k}^{(\ell)}$ is a momentum degree of freedom in $\Gamma^{(\ell)*}$,
the monolayer Brillouin zone of layer $\ell$ with area $|\Gamma^{*}|$. We
first note the orthogonality relation

$\displaystyle\begin{aligned}
\sum_{\mu\alpha}\delta\widehat{u}_{i\ell\mu\alpha}^{*}\delta\widehat{u}_{j\ell\mu\alpha}=\delta_{ij}.\end{aligned}$
(4)

To expand the equation of motion, Eq. (2), in the Fourier basis, we multiply
both sides by $\delta\widehat{u}_{k\ell\mu\alpha}^{*}$, then sum over $k$,
$\mu$, and $\alpha$.

$\displaystyle\sum_{jk\mu\nu\alpha\beta}\overline{D}_{\mu\nu\alpha\beta}^{\ell}(\bm{R}_{i}-\bm{R}_{j})\delta\widehat{u}_{k\ell\mu\alpha}^{*}\delta\widehat{u}_{j\ell\nu\beta}$
$\displaystyle=\omega^{2}\sum_{k\mu\alpha}\delta\widehat{u}_{k\ell\mu\alpha}^{*}\delta\widehat{u}_{i\ell\mu\alpha}.$
(5)

Using Eq. (4):

$\displaystyle\sum_{ij\mu\nu\alpha\beta}\overline{D}_{\mu\nu\alpha\beta}^{\ell}(\bm{R}_{i}-\bm{R}_{j})\delta\widehat{u}_{i\ell\mu\alpha}^{*}\delta\widehat{u}_{j\ell\nu\beta}$
$\displaystyle=\omega^{2}.$ (6)

Relabeling $\bm{R}=\bm{R}_{i},\Delta\bm{R}=\bm{R}_{i}-\bm{R}_{j}$ and changing
to the Bloch basis using Eq. (3) yields

$\displaystyle\frac{1}{|\Gamma^{*}|}\sum_{\bm{R}}\left[e^{-i(\bm{k}^{(\ell)}-\bm{k}^{\prime(\ell)})\cdot\bm{R}}\sum_{\mu\nu\alpha\beta}\sum_{\Delta\bm{R}}\overline{D}_{\mu\nu\alpha\beta}(\Delta\bm{R})e^{-i\bm{k}^{\prime(\ell)}\cdot\Delta\bm{R}}\delta
u_{\ell\mu\alpha}^{*}(\bm{k}^{(\ell)})\delta
u_{\ell\nu\beta}(\bm{k}^{\prime(\ell)})\right]$ $\displaystyle=\omega^{2}.$
(7)

Note the Poisson resummation rule:

$\displaystyle\sum_{\bm{R}}e^{-i(\bm{k}^{(\ell)}-\bm{k}^{\prime(\ell)})\cdot\bm{R}}=|\Gamma^{*}|\sum_{\bm{G}^{(\ell)}}\delta(\bm{k}^{(\ell)}-\bm{k}^{\prime(\ell)}-\bm{G}^{(\ell)}),$
(8)

where $\bm{G}^{(\ell)}$ is a reciprocal lattice vector in layer $\ell$.
Defining $D_{\mu\nu\alpha\beta}(\bm{k}^{\prime(\ell)})$ as:

$\displaystyle
D_{\mu\nu\alpha\beta}(\bm{k}^{\prime(\ell)})=\sum_{\Delta\bm{R}}\overline{D}_{\mu\nu\alpha\beta}(\Delta\bm{R})e^{-i\bm{k}^{\prime(\ell)}\cdot\Delta\bm{R}},$
(9)

the phonon equation of motion becomes

$\displaystyle\begin{aligned}
\sum_{\mu\nu\alpha\beta}D_{\mu\nu\alpha\beta}(\bm{k}^{(\ell)})\delta
u_{\ell\mu\alpha}^{*}(\bm{k}^{(\ell)})\delta
u_{\ell\nu\beta}(\bm{k}^{(\ell)})=\omega^{2},\end{aligned}$ (10)

for each $\bm{G}^{(\ell)}$, since the phonon displacement $\delta
u_{i\ell\mu\alpha}$ is the same at $\bm{k}^{(\ell)}$ and
$\bm{k}^{(\ell)}+\bm{G}^{(\ell)}$. Equation (10) is equivalent to:

$\displaystyle\begin{aligned}
D(\bm{k}^{(\ell)})\cdot\delta\bm{u}(\bm{k}^{(\ell)})=\omega^{2}\delta\bm{u}(\bm{k}^{(\ell)}),\end{aligned}$
(11)

which is the eigenvalue equation that describes the intralayer phonon at
$\bm{k}^{(\ell)}$, and the dynamical matrix $D(\bm{k}^{(\ell)})$ can be
computed using either DFPT or frozen phonon approaches.

In a moiré cell the dynamical matrix gains an additional degree of freedom
because the forces depend not only on the difference in lattice vectors, but
also the real space position itself due to the twist. The general phonon
equation of motion is given as the following:

$\displaystyle\sum_{j\nu\beta\ell^{\prime}}\overline{D}_{\mu\nu\alpha\beta}(\bm{R}_{i},\bm{R}_{j})\delta\widehat{u}_{j\ell^{\prime}\nu\beta}=\omega^{2}\delta\widehat{u}_{i\ell\mu\alpha},$
(12)

where $\bm{R}_{i}$ and $\bm{R}_{j}$ correspond to the atomic positions on
layers $\ell$ and $\ell^{\prime}$ respectively. Similar to the intralayer
term, we multiply both sides of the equation by
$\delta\widehat{u}^{*}_{i\ell\mu\alpha}$ and sum over $\mu$, $\alpha$, and
$\ell$. Plugging in the definition in Eq. (3):

$\displaystyle\sum_{\mu\nu\alpha\beta\ell\ell^{\prime}}D_{\widetilde{\bm{k}}\mu\nu\alpha\beta}(\bm{k}^{(\ell)},\bm{k}^{(\ell^{\prime})})\delta
u^{*}_{\ell\mu\alpha}(\bm{K}^{(\ell)})\delta
u_{\ell^{\prime}\nu\beta}(\bm{K}^{(\ell^{\prime})})=\omega^{2},$ (13)

where

$\displaystyle\begin{aligned}
D_{\widetilde{\bm{k}}\mu\nu\alpha\beta}(\bm{k}^{(\ell)},\bm{k}^{(\ell^{\prime})})=\frac{1}{|\Gamma^{*}|}\sum_{ij}\overline{D}_{\mu\nu\alpha\beta}(\bm{R}_{i},\bm{R}_{j})e^{-i\bm{K}^{(\ell)}\cdot\bm{R}_{i}}e^{i\bm{K}^{(\ell^{\prime})}\cdot\bm{R}_{j}}.\end{aligned}$
(14)

and $\bm{K}^{(\ell)}=\widetilde{\bm{k}}+\bm{k}^{(\ell)}$ is an expansion for
small momentum $\bm{k}^{(\ell)}$ about a given point $\widetilde{\bm{k}}$ in
the moiré Brillouin zone. Making a two-center approximation such that
$\overline{D}_{\mu\nu\alpha\beta}(\bm{R}_{i},\bm{R}_{j})=\overline{D}_{\mu\nu\alpha\beta}(\bm{R}_{i}-\bm{R}_{j})$,
we can then perform a Fourier expansion on
$\overline{D}_{\mu\nu\alpha\beta}(\bm{R}_{i}-\bm{R}_{j})$:

$\displaystyle\overline{D}_{\mu\nu\alpha\beta}(\bm{R}_{i}-\bm{R}_{j})=\int\frac{\differential^{2}{\bm{p}}}{(2\pi)^{2}}\,e^{i\bm{p}\cdot(\bm{R}_{i}-\bm{R}_{j})}D_{\mu\nu\alpha\beta}(\bm{p}),$
(15)

Note that despite the two-center approximation, the momentum-space dynamical
matrix $D_{\tilde{\bm{k}}\mu\nu\alpha\beta}$ still depends on two different
momenta from the opposite layers $\ell$ and $\ell^{\prime}$. Using Eq. (15),
as well as Eq. (8) in Eq. (14) yields

$\displaystyle\begin{aligned}
D_{\widetilde{\bm{k}}\mu\nu\alpha\beta}(\bm{k}^{(\ell)},\bm{k}^{(\ell^{\prime})})&=\frac{1}{|\Gamma^{*}|}\int\frac{\differential^{2}{\bm{p}}}{(2\pi)^{2}}\;D_{\mu\nu\alpha\beta}(\bm{p})\sum_{ij}e^{-i(\bm{K}^{(\ell)}+\bm{p})\cdot\bm{R}_{i}}e^{i(\bm{K}^{(\ell^{\prime})}+\bm{p})\cdot\bm{R}_{j}}\\\
&=\frac{|\Gamma^{*}|}{(2\pi)^{2}}\sum_{\bm{G}^{(\ell)}\bm{G}^{(\ell^{\prime})}}\int\differential^{2}{\bm{p}}\,D_{\mu\nu\alpha\beta}(\bm{p})\,\delta(\bm{K}^{(\ell)}+\bm{p}-\bm{G}^{(\ell)})\,\delta(\bm{K}^{(\ell^{\prime})}+\bm{p}-\bm{G}^{(\ell^{\prime})})\\\
&=\frac{1}{|\Gamma|}\sum_{\bm{G}^{(\ell)}\bm{G}^{(\ell^{\prime})}}D_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}+\bm{k}^{(\ell)}-\bm{G}^{(\ell)})\,\delta_{\bm{k}^{(\ell^{\prime})}-\bm{k}^{(\ell)},\bm{G}^{(\ell^{\prime})}-\bm{G}^{(\ell)}},\\\
\end{aligned}$ (16)

where $|\Gamma|$ is the area of the monolayer unit cell and we note that
$|\Gamma|=(2\pi)^{2}/|\Gamma^{*}|$. The $\delta$-function imposes a selection
rule that constrains the allowed values of $\bm{k}^{(\ell)}$ and
$\bm{k}^{(\ell^{\prime})}$. For concreteness, we define each reciprocal
lattice vector as

$\displaystyle\bm{G}^{(\ell)}_{mn}=m\bm{g}^{(\ell)}_{1}+n\bm{g}_{2}^{(\ell)},$
(17)

where $\bm{g}_{1}^{(\ell)}$ and $\bm{g}_{2}^{(\ell)}$ are a reciprocal lattice
basis for layer $\ell$. Neglecting higher-order coupling, we assume that only
terms with equal $m,n$ in $\bm{G}^{(\ell)}$ and $\bm{G}^{(\ell^{\prime})}$ are
nonzero, so the above simplifies to

$\displaystyle
D_{\widetilde{\bm{k}}\mu\nu\alpha\beta}(\bm{k}^{(\ell)},\bm{k}^{(\ell^{\prime})})=\frac{1}{|\Gamma|}\sum_{mn}$
$\displaystyle
D_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}+\bm{k}^{(\ell)}-\bm{G}^{(\ell)}_{mn})\,\delta_{\bm{k}^{(\ell)}-\bm{k}^{(\ell^{\prime})},\widetilde{\bm{G}}_{mn}}.$
(18)

The difference $\bm{G}^{(\ell)}_{mn}-\bm{G}^{(\ell^{\prime})}_{mn}$ between
corresponding reciprocal lattice vectors from each layer is a moiré reciprocal
lattice vector $\widetilde{\bm{G}}_{mn}$. Without loss of generality, taking
$\ell=1$ and $\ell^{\prime}=2$, Eq. (18) imposes a constraint on the momenta
given by

$\displaystyle\bm{k}^{(1)}-\bm{k}^{(2)}=\bm{G}^{(2)}_{mn}-\bm{G}^{(1)}_{mn}=-\widetilde{\bm{G}}_{mn}.$
(19)

We henceforth drop the subscripts on $\widetilde{\bm{G}}$ for clarity.
Equation (18) is the fundamental equation governing the moiré dynamical
matrix. The stacking-dependent part of the moiré phonon equation of motion
becomes

$\displaystyle\begin{aligned}
\frac{1}{|\Gamma|}\sum_{\nu\beta\ell^{\prime}}\sum_{\widetilde{\bm{G}}}D_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}+\bm{k}^{(\ell)}-\bm{G}^{(\ell)})\delta
u_{\ell^{\prime}\nu\beta}(\bm{K}^{(\ell^{\prime})}-\widetilde{\bm{G}})&=\omega^{2}\delta
u_{\ell\mu\alpha}(\bm{K}^{(\ell)}),\end{aligned}$ (20)

where we note that $\widetilde{\bm{G}}$ and $\bm{G}^{(\ell)}$ are in one-to-
one correspondence via Eq. (19). To simplify the above we choose a gauge
$\bm{k}^{(\ell)}=\bm{G}^{(\ell^{\prime})}$ (though other gauges will work, and
we will choose a different one to justify the construction of the moiré
dynamical matrix) with $\ell\neq\ell^{\prime}$. Without loss of generality, we
let $\ell=1$:

$\displaystyle\begin{aligned}
D_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}+\bm{k}^{(1)}-\bm{G}^{(1)})&=D_{\mu\nu\alpha\beta}(\bm{k}+\bm{G}^{(2)}-\bm{G}^{(1)})\\\
&=D_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}-\widetilde{\bm{G}}).\end{aligned}$
(21)

Note that reciprocal space has sufficient symmetry that any choice of a pair
of signs for $\bm{k}$ and $\widetilde{\bm{G}}$ produce the same term above.

Combining the monolayer term (Eq. (11)) and the twist-dependent term (Eq.
(20)), we obtain the moiré dynamical matrix. We verified numerically that the
magnitude of the dynamical matrix decays rapidly with distance from the
$\Gamma$-point, so it suffices to expand up to the first shell of the moiré
reciprocal lattice, which contains six vectors
$\widetilde{\bm{G}}_{1},\ldots,\widetilde{\bm{G}}_{6}$. For the off-diagonal
terms, we ignore $\bm{\widetilde{k}}$-dependence. The interlayer term is
explicitly given as follows:

$\displaystyle
D^{12}=\begin{bmatrix}D_{\mu\nu\alpha\beta}^{12}(\widetilde{\bm{G}}_{0})&D_{\mu\nu\alpha\beta}^{12}(\widetilde{\bm{G}}_{1})&\cdots&D_{\mu\nu\alpha\beta}^{12}(\widetilde{\bm{G}}_{6})\\\
D_{\mu\nu\alpha\beta}^{12}(-\widetilde{\bm{G}}_{1})&D_{\mu\nu\alpha\beta}^{12}(\widetilde{\bm{G}}_{0})&&\\\
\vdots&&\ddots&\\\
D_{\mu\nu\alpha\beta}^{12}(-\widetilde{\bm{G}}_{6})&&&D_{\mu\nu\alpha\beta}^{12}(\widetilde{\bm{G}}_{0})\end{bmatrix},$
(22)

where $D^{ij}_{\mu\nu\alpha\beta}(\widetilde{\bm{G}})$ is a submatrix of the
dynamical matrix at $\widetilde{\bm{G}}$ with rows (columns) corresponding to
degrees of freedom due to layer $i$ ($j$). The $(i,j)$ block thus represents
phonon hoppings between reciprocal lattice sites $\bm{\widetilde{G}}_{i}$ and
$\bm{\widetilde{G}}_{j}$. We have disregarded a small number of first-shell
terms above that do not live along the first row, first column, or the
diagonal, after having verified that their inclusion negligibly affects the
numerical results. That is, we ignore hoppings that do not include
$\bm{\widetilde{G}}_{0}$. Similarly for the intralayer term,

$\displaystyle
D^{i}(\bm{\widetilde{k}})=\begin{bmatrix}D_{P}^{i}(\bm{\widetilde{k}}+\widetilde{\bm{G}}_{0})&D_{\mu\nu\alpha\beta}^{ii}(\widetilde{\bm{G}}_{1})&\cdots&D_{\mu\nu\alpha\beta}^{ii}(\widetilde{\bm{G}}_{6})\\\
D_{\mu\nu\alpha\beta}^{ii}(-\widetilde{\bm{G}}_{1})&D_{P}^{i}(\bm{\widetilde{k}}+\widetilde{\bm{G}}_{1})&&\\\
\vdots&&\ddots&\\\
D_{\mu\nu\alpha\beta}^{ii}(-\widetilde{\bm{G}}_{6})&&&D_{P}^{i}(\bm{\widetilde{k}}+\widetilde{\bm{G}}_{6})\end{bmatrix},$
(23)

where the diagonal terms are the dynamical matrix of the pristine monolayer
$D^{i}_{P}(\bm{\widetilde{k}}+\bm{\widetilde{G}})$ and the off-diagonal terms
being the twist-dependent intralayer term constructed from Eq. (20) when
$\ell=\ell^{\prime}$.

###  II.3  Dynamical Matrix in Configuration Space

We are particularly interested in structures with a small twist angle, which
means the moiré cell contains a large number of primitive cells and
consequently many atoms. In order to make the model computationally feasible,
we express $D_{\mu\nu\alpha\beta}$ in terms of a collection of local atomic
environments, known as configuration space [41]. Within the continuum
approximation we replace the primitive cell lattice vectors $\mathbf{R}$ with
a continuous variable $\bm{r}$ defined throughout the moiré cell, and allow
quantities of interest to be defined at all $\bm{r}$. At each location
$\bm{r}$ the local environment can be approximated by an untwisted bilayer
with a horizontal shift of one layer relative to the other. We specify the
interlayer shift by $\bm{b}(\bm{r})$, and the space formed by all $\bm{b}$ is
configuration space. It is in one-to-one correspondence with real space by a
linear transformation

$\displaystyle\bm{b}(\bm{r})=(\mathbb{1}-A_{2}A_{1}^{-1})\bm{r},$ (24)

where $A_{i}$ is a matrix whose columns are the primitive lattice vectors of
layer $i$ and $\mathbb{1}$ is the $2\times 2$ identity matrix. This
correspondence also gives a useful relationship between the the reciprocal
lattice vectors of the moiré cell, $\bm{\widetilde{G}}$, and those of the
primitive unit cell, $\bm{G}$:

$\displaystyle\bm{\widetilde{G}}\cdot\bm{r}=\bm{G}\cdot\bm{b}(\bm{r}).$ (25)

Using these concepts, we can express the momentum-space dynamical matrix as
follows:

$\displaystyle\begin{aligned}
D^{\text{rig}}_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}+\widetilde{\bm{G}})=&\frac{1}{\sqrt{|\Gamma^{*}|}}\int
d^{2}\bm{r}\,\overline{D}_{\mu\nu\alpha\beta}(\bm{r})\,e^{i(\tilde{\bm{k}}+\tilde{\bm{G}})\cdot\bm{r}}\\\
=&\frac{1}{\sqrt{|\Gamma^{*}|}}\int
d^{2}\bm{b}\,\overline{D}_{\mu\nu\alpha\beta}(\bm{b})\,e^{i(\bm{k}+\bm{G})\cdot\bm{b}}\\\
\approx&\frac{1}{\sqrt{|\Gamma^{*}|}}\frac{|\Gamma|}{N^{2}}\sum_{\bm{b}\in\mathcal{S}(N)}\overline{D}_{\mu\nu\alpha\beta}^{(\bm{b})}\,e^{i(\bm{k}+\bm{G})\cdot\bm{b}},\end{aligned}$
(26)

The factor of $|\Gamma|$, arising from discretization, cancels out the
normalization from Eq. (18). Configuration space allows us to replace
$\overline{D}(\bm{r})$, which we cannot easily determine, with
$\overline{D}(\bm{b})$, which we can readily compute with DFT for a primitive
bilayer cell with any specified shift $\bm{b}$. In the small $\theta$ limit
the local environment varies slowly and smoothly with $\bm{b}$, so we can
approximate the integral over $\bm{b}$ with the sum over a $N\times N$ uniform
mesh $\mathcal{S}(N)$ of $\bm{b}$ values. Since $\mathcal{S}(N)$ is uniform,
$D^{\text{rig}}_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}+\widetilde{\bm{G}})$
describes a moiré bilayer with a uniform distribution of local configurations.
In other words, it represents a rigid moiré structure with no twist-induced
continuum relaxation, which is a good approximation for $\theta$ sufficiently
large. We describe the procedure for including relaxation effects in the next
subsection.

###  II.4 Continuum Relaxation

For moderate to small twist angles, twist-induced atomic relaxation becomes
important. The range varies by material, but relaxation typically begins to
have a significant impact below about $2^{\circ}$. In real space, atomic
relaxation results in enlarged areas of the lowest energy stacking
configurations (AB and BA in graphene) at the expense of high energy stackings
(AA), with domain walls forming in between. In configuration space relaxation
can be viewed as a higher density of configurations in the lowest-energy
stackings [41, 42]. This can be accounted for by shifting each $\bm{b}$ in the
uniform mesh $\mathcal{S}(N)$ by an angle-dependent relaxation displacement
vector $\bm{u}_{\theta}(\bm{b})$, yielding a non-uniform mesh
$\bm{B}=\bm{b}+\bm{u}_{\theta}(\bm{b})$ denoted $\mathcal{S}_{\theta}(N)$. We
use linear elasticity theory to describe the intralayer energy, combined with
the generalized stacking fault energy (GSFE) functional for the interlayer
energy. (The GSFE parameters are given in Appendix A.) By minimizing the total
energy we can compute $\bm{u}_{\theta}(\bm{b})$ in configuration space [41].
The dynamical matrix can then be updated by Fourier interpolation, using an
inverse transform

$\displaystyle\overline{D}_{\mu\nu\alpha\beta}^{(\bm{B})}(\bm{\widetilde{k}})=\frac{1}{\sqrt{|\Gamma^{*}|}}\sum_{\bm{\widetilde{G}}}D^{\text{rig}}_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}+\widetilde{\bm{G}})e^{-i(\bm{k}+\bm{G})\cdot\bm{B}},$
(27)

followed by a second forward transform, which arrives at the following relaxed
moiré dynamical matrix:

$\displaystyle
D_{\mu\nu\alpha\beta}(\widetilde{\bm{k}}+\widetilde{\bm{G}})=\frac{1}{\sqrt{|\Gamma^{*}|}}\frac{|\Gamma|}{N^{2}}\sum_{\bm{b}\in\mathcal{S}(N)}\overline{D}_{\mu\nu\alpha\beta}^{(\bm{B})}(\widetilde{\bm{k}})e^{i(\bm{k}+\bm{G})\cdot\bm{b}}.$
(28)

To obtain the real space phonon vibration pattern, we perform an inverse
Fourier transform by summing over components of the eigenvectors as follows,

$\displaystyle\delta
u_{n\ell\mu\alpha\widetilde{\bm{k}}}(\bm{r})=\sum_{\bm{\widetilde{G}}}\delta
u_{n\ell\mu\alpha\widetilde{\bm{G}}}(\widetilde{\bm{k}})e^{i(\bm{\widetilde{G}}+\bm{\widetilde{k}})\cdot\bm{r}},$
(29)

where $n$ is the band index and $\delta u_{n\ell\mu\alpha\widetilde{\bm{G}}}$
is the eigenvector component that corresponds to $\widetilde{\bm{G}}$.

###  II.5 Acoustic Sum Rules

Even though the total forces are fixed to be zero in the force field
calculations of individual stackings, the combination of all configurations in
the moiré dynamical matrix construction modifies the total force of the
combined system. Therefore, an acoustic sum of the moiré dynamical matrix is
required to ensure the total force stays zero in addition to the sum rules on
the unrotated bilayers and monolayers. We use a correction matrix (chosen
without loss of generality to be diagonal in the latter pair of indices)

$\displaystyle
C_{\mu\nu\alpha\alpha}=\sum_{\beta}\sqrt{\frac{M_{\beta}}{M_{\alpha}}}D_{\text{moir\'{e}},\mu\nu\beta\alpha}^{*}(\Gamma),$
(30)

where $D_{\text{moir\'{e}},\mu\nu\beta\alpha}^{*}(\Gamma)$ is the moiré
dynamical matrix by replacing the pristine monolayer blocks
$D_{P}(\bm{k}+\bm{\widetilde{G}})$ with $D_{P}(\Gamma)$, as we enforce
translational invariance only at the $\Gamma$ point and enforce continuity
everywhere else. We then adjust
$D_{\text{moir\'{e}}}(\bm{\widetilde{k}})\leftarrow
D_{\text{moir\'{e}}}(\bm{\widetilde{k}})-C$, observing numerically that the
above approximation satisfies adequately the translational invariance
constraint. The sum rule is adapted from Quan  _et al._ [35].

##  III Model Validation

To validate our CSC model, we compare results obtained from it with two first-
principles models, namely direct DFT on the moiré supercell for a small system
where such calculations are feasible, and with molecular dynamics (MD)
simulations that employ an interatomic potential [28, 47]. We also compare our
CSC results to another empirical continuum (EC) model adapted from Koshino and
Son [29].

###  III.1 Comparison with Molecular Dynamics and DFT

Direct DFT phonon calculations on the moiré cell are the most accurate but are
restricted by the size of the cell. MD simulations allow for larger supercell
sizes, but have inherent limitations on accuracy due to the choice of
interatomic potential. Both direct DFT and MD are computationally expensive
and require a commensurate supercell. Since the CSC model uses DFT inputs in a
well-defined approximation, it should reflect the direct DFT and MD
calculations reasonably well and does not suffer constraints based on moiré
supercell size or commensurability. We choose a twist angle,
$\theta=7.34^{\circ}$, that gives a commensurate pattern in bilayer graphene
for these comparisons.

Figure 2: Comparison of phonon band at $\theta=7.34^{\circ}$ between: DFT
(left panel), CSC (middle panel), and MD (right panel) calculations for
graphene. The LB mode is highlighted in red.

Figure 2 compares the phonon bands of the CSC model, MD simulations, and
direct DFT phonon calculations on the moiré supercell at
$\theta=7.34^{\circ}$. The layer breathing (LB) mode, which represents layer 1
(2) moving in the $-z$ ($+z$) direction and a much weaker in-plane motion, is
highlighted in red in each case. The results from the three models are
qualitatively similar, and to a reasonable extent quantitatively close to each
other, across the entire Brillouin Zone.

Figure 3: Comparison of the LB atomic displacements between direct DFT, CSC,
and MD. Each model has normalized its displacements $\delta\bm{u}$ separately,
in arbitrary units. The insets show the frequency and average displacement
magnitude in-plane ($\delta\overline{u}_{xy}$) and out-of-plane
($\delta\overline{u}_{z}$).

In addition to the phonon dispersion relations, Fig. 3 shows a comparison of
the atomic displacements of the LB mode in graphene, corresponding to the red
band in Fig. 2 at the $\Gamma$-point. We focus on a qualitative comparison of
the geometrical structure. All models exhibit a swirling in-plane motion
around the AA regions at the corner of the moiré supercell. The primary
difference between the models is the variation in the out-of-plane motion. In
the DFT calculation the maximum vertical displacement occurs in the AB
regions, whereas for the CSC model the maximum is in the AA regions. While the
MD calculation shows some variation in the vertical displacement between AA
and AB spots and the maximum vertical displacement also occurs in the AB
regions as in DFT, it is mostly uniform throughout the moiré cell. The reason
for the quantitative difference between the CSC and the first-principles
models is likely due to the breakdown of the continuum approximation at a
large angle. When the twist angle is large, as is the case for
$\theta=7.34^{\circ}$ in graphene, the local stacking order varies drastically
from one real space position to a neighboring one. However, the CSC still
assumes a smooth variation of the force fields in real space, which causes the
real space LB mode to be similar for all twist angles and thus it fails to
capture the change in the large angle. Despite this slight disagreement, both
the LB frequency and the band structure from the CSC show an excellent
agreement with first-principles calculations (Fig. 3). There agreement in the
LB mode out-of-plane displacement field is better at smaller twist angles (not
shown in Fig. 3), that is, all three models exhibit maxima in the displacement
field around the AA regions.

###  III.2 Comparison with Empirical Continuum Model

The EC model, based on empirical interlayer interactions, is the most
efficient of the models discussed in this section [29]. In Fig. 4 we compare
the phonon dispersion relations and density of states for MoS2, obtained with
the EC model and with our CSC model. To obtain the in-plane phonon modes, we
perturb around the relaxed equilibrium positions by adding a kinetic energy
contribution to the interlayer and intralayer terms of the dynamical matrix.
We then solve the equation of motion to obtain the phonon modes. We use the
GSFE coefficients and elastic constants provided in Tables 1 and 2 for the CSC
model. A detailed derivation is given in Appendix D.

Figure 4: Comparison of the phonon bands and DOS for bilayer MoS2 at (a)
$0.75^{\circ}$ and (b) $8.0^{\circ}$ as calculated using the EC model (red
lines) and the CSC model (black lines). For $\theta=0.75^{\circ}$, the EC DOS
has been scaled by a factor of 1000 for visualization.

The EC model has a free empirical parameter that we determine by scaling the
bands such that the lowest-frequency in-plane translational modes at K
coincide in the two models; the scaling factor is 2.3 for
$\theta=0.75^{\circ}$ and 22.5 for $8.0^{\circ}$.

There is little agreement between the two models for small $\theta$, as shown
in Fig. 4(a). The EC model predicts flatter bands and exhibits a gap between 2
and 3 cm-1, whereas no such features occur in the CSC model. Moreover, the
slopes of the in-plane acoustic modes cannot be matched well. Intrinsically,
the EC model also does not include any out-of-plane motion, such as the
$z$-translational mode, nor does it include interlayer modes such as the
shearing (S) mode. The LB mode is absent in the EC model for these two
reasons. Our continuum model does not suffer from such limitations.

In the large twist-angle limit, Fig. 4(b), the results from the two models
agree for low frequencies near $\Gamma$; they both exhibit linear DOS
relations and the position of the first folded band cluster is near 35 cm-1.

Figure 5: Phonon bands of graphene as obtained from our CSC model at four
representative angles of (a) $8^{\circ}$, (b) $2^{\circ}$, (c) $1.1^{\circ}$,
and (d) $0.4^{\circ}$. The two shearing modes are highlighted in red, and the
first (lowest energy) LB mode in blue.

##  IV Moiré Phonon Properties

Using our CSC model, we calculate the phonon properties of three different
moiré bilayers: graphene, MoS2, and a MoSe2-WSe2 heterostructure. Near the
$\Gamma$-point the lowest energy modes include two shearing modes (S) and one
layer breathing (LB) mode. These represent in-plane and out-of-plane motions,
respectively, with each layer moving in opposite direction.

###  IV.1 Low-energy mode frequencies

The phonon dispersion relations are qualitatively similar for twisted bilayer
graphene and TMDCs. Here, we focus on graphene. Figure 5 shows the graphene
phonon bands at different twist angles $\theta$. As $\theta$ decreases, the
frequency where folded clusters of phonon bands appear also decreases, due to
the shrinking of the moiré Brillouin zone area. Depending on the amount of
atomic relaxation introduced by the twist, the values of $\theta$ may be
partitioned into three regimes: a “decoupled” regime (large $\theta$), a
“soft” regime (moderate $\theta$) and a “relaxed” regime (small $\theta$). The
exact values of $\theta$ separating these regimes are specific to the
material. For graphene they are approximately $4.2^{\circ}$ and $0.5^{\circ}$,
respectively; for MoS2 they are $7.6^{\circ}$ and $1^{\circ}$; and for
MoSe2-WSe2 they are $5.5^{\circ}$ and $1.0^{\circ}$. In the decoupled regime
(Fig. 5(a)) the two layers are twisted sufficiently so that interlayer
coupling is very weak. Here, the S mode is slightly positive at $\Gamma$ with
twofold degeneracy (when the distinction between the two S modes becomes
important, we shall refer to the lower-energy curve as S1 and the higher-
energy curve as S2), while the LB mode is essentially dispersionless. The
zero-frequency modes at $\Gamma$ are the three translational modes. In the
soft regime (Fig. 5(b)), the S mode begins to exhibit imaginary frequencies at
the $\Gamma$ point (plotted as negative values for continuity in the phonon
bands), with the negativity increasing until the relaxed regime is reached
(Fig. 5(d)), at which point the S modes return to positive values. The
negative shearing mode frequencies are most likely numerical artifacts due to
the finite discretization of the local configuration for the following
reasons: at large twist angles, the interlayer interaction is negligible, and
thus the phonon band structure should be nearly identical to folded AB bilayer
graphene bands (see Fig. 5(a)). At these intermediate to large angles where
the shearing mode frequency is negative, the local configuration varies
relatively quickly on the unit cell scale within a moiré supercell. However,
force fields can have complex structure as a function of the local stacking
and some off-diagonal components may not respect the crystal symmetry (unlike
the GSFE, which can be captured by fitting the first few Fourier components).
Our $9\times 9$ sampling grid may not be dense enough to capture the full
structure of the force fields since phonon frequencies depend very sensitively
on the accuracy of the forces. In the relaxed regime, the effect of relaxation
becomes strong, and the local stacking varies smoothly on the moiré scale and
is dominated by the equilibrium AB/BA stackings. As a result, the shearing
mode frequeny becomes positive again.

In the soft regime, the LB mode becomes dispersive, as seen in Fig. 5(b)-(c),
but progressively flattens out again as $\theta$ approaches the relaxed
regime, Fig. 5(d). At $1.1^{\circ}$, the LB mode has a crossing with a cluster
of folded bands whose frequency decreases as the Brillouin zone shrinks
(5(c)), and similarly for $\theta=0.4^{\circ}$ (5(d)).

Figure 6: Frequency of the shearing (red line) and layer breathing (blue line)
modes at $\Gamma$ as a function of twist angle $\theta$ for various bilayers.
The MoSe2-WSe2 heterostructure has an additional LB-like mode that breaks
symmetry, labeled SBLB (black line). The background color shades illustrate
the transition from the decoupled region (large $\theta$, S frequencies
positive or nearly 0) in blue shade, to the soft region (from initial dip to
negative S frequencies to angle of maximally negative frequency) in white
shade, to the relaxed region (upwards curve to positive S frequencies) in
brown shade.

We may examine the frequency dependence on twist angle more thoroughly by
focusing on the $\Gamma$ point. Figure 6 compares the S and LB frequencies for
graphene, $\mathrm{MoS_{2}}$, and MoSe2-WSe2. The decoupled, soft, and relaxed
regimes are shaded blue, white, and brown, respectively. The soft regime
begins when the S mode takes negative frequency values, and ends where the
slope of $\omega$ becomes negative at small $\theta$. Aside from the location
and width of the soft regime, however, the frequency dependence on $\theta$ of
the S mode is qualitatively similar for all three materials. For MoS2 and
MoSe2-WSe2, the LB varies smoothly and slowly with $\theta$. By contrast, the
LB mode in graphene has a comparatively higher frequency in the decoupled
regime, which drops sharply near the end of the soft regime. Since the
MoSe2-WSe2 heterostructure breaks the layer symmetry, it admits an additional
LB-like mode that does not respect the symmetry under layer inversion
($\bm{r}\rightarrow-\bm{r}$). We refer to this mode as the symmetry-breaking
LB (SBLB) mode, and note that it is the only mode that varies linearly in
frequency with respect to $\theta$.

###  IV.2 Real Space Atomic Displacements at the $\Gamma$ Point

By performing the inverse Fourier transform of eigenmodes in reciprocal space
(Eq. (29)), we obtain the atomic displacement field of each phonon mode in
real space. The eigenmodes in reciprocal space are unit vectors, but
application of the Fourier transform gives the magnitude of the displacement
vectors a physical interpretation. Specifically, the overall magnitudes across
various twist angles serve as an indicator for the magnitude of each eigenmode
component $\delta\bm{u}(\widetilde{\bm{G}})$ in reciprocal space. For if most
of the magnitude in reciprocal space is concentrated in the
$\widetilde{\bm{G}}_{0}$ component, the real space magnitude will be large;
but if most of the magnitude is concentrated in the first shell
$\widetilde{\bm{G}}_{1},\ldots,\widetilde{\bm{G}}_{6}$, the phases in the
inverse Fourier transform reduce the real space magnitude to a much smaller
quantity. In the absence of the moiré interlayer couplings, the magnitude in
Fourier space should be concentrated in the $\widetilde{\bm{G}}_{0}$
component, but at smaller values of $\theta$, emergence of periodicity of the
moiré length scale imply that the magnitude becomes instead concentrated in
higher shells. Thus, we should expect the magnitude in real space to fall
substantially near the angle at which the magnitudes of the
$\widetilde{\bm{G}}_{0}$ component and the
$\widetilde{\bm{G}}_{1},\ldots,\widetilde{\bm{G}}_{6}$ components intersect.

Figure 7: Average atomic displacement magnitudes in real space
($\delta\overline{u}$, red), of the $\widetilde{\bm{G}}_{0}$ component in
reciprocal space (blue), and of the average of the first-shell components
$\widetilde{\bm{G}}_{1},\ldots,\widetilde{\bm{G}}_{6}$ of the LB mode at
$\Gamma$ for MoS2, graphene, and MoSe2-WSe2. The shading of various regimes
has the same meaning as in Fig. 6. The vertical black lines indicate the
largest value of the twist angle at which $\delta\overline{u}\to 0$.

Figure 7 demonstrates these features in the LB mode by approximating the
average real space magnitude through a $13\times 13$ mesh discretization of
the moiré supercell. For brevity, we have averaged the magnitudes of the six
first-shell reciprocal lattice vectors. Note that the magnitude of the
$\widetilde{\bm{G}}_{0}$ component can be measured via Raman spectroscopy
[35], and thus the real space magnitude across $\theta$ serves to indicate at
which twist angle the mode will be difficult to detect experimentally.

The magnitudes of $\delta\bm{u}(\bm{r})$ also allow for comparison of the
relative strength between in-plane ($\delta u_{xy}$) and out-of-plane ($\delta
u_{z}$) motion; for example, the S modes satisfy $\delta u_{z}\ll\delta
u_{xy}$. This is seen explicitly by visualization of the S and LB atomic
displacement fields, from which we may also deduce the geometric structure of
the atomic displacements as a function of $\theta$. In all figures showing
atomic displacement field results from the CSC model, we normalize the
magnitudes that correspond to each atomic degree of freedom by the atomic mass
$\sqrt{M_{\alpha}}$ for each layer, and subsequently average the vectors over
the atomic degrees of freedom. We verified numerically that for all modes
analyzed in the present work, the vectors for each atomic degree of freedom
are identical up to mass normalization, which justifies the averaging. We then
renormalize by multiplying the atom-averaged field by
$\sqrt{\sum_{\alpha}M_{\alpha}}$, where $\alpha$ is the atomic index.

Figure 8: Real space atomic displacements (in arbitrary units) for the twofold
degenerate shearing modes, S1 and S2, of the first layer of graphene at
various twist angles. In-plane ($x,y$) displacements are proportional to the
arrows and vertical ($z$) displacements are shown in color. The arrows are
rescaled for visualization purposes, with the average values of in-plane
($\delta\overline{u}_{xy}$) and out-of-plane ($\delta\overline{u}_{z}$)
magnitudes, rounded to 2 decimal places, shown in the lower insets. The
displacements of the layer 2 are exactly opposite to those of layer 1. The
grey parallelogram outlines a single moiré cell.

Figure 8 illustrates the evolution of the atomic displacement field of the two
S modes, S1 and S2, as a function of $\theta$. In all cases, the S mode is
dominated by in-plane displacements and is opposite in the two layers. In the
decoupled regime ($7.34^{\circ}$) the displacement is uniform. As $\theta$
decreases, translational symmetry on the graphene unit cell scale is broken
and we observe structure at the moiré scale.

Figure 9: LB atomic displacements of graphene at twist angles (a)
$1.00^{\circ}$ and (b) $7.34^{\circ}$. Color scale and insets are the same as
in Fig. 8. For each angle, the left (right) panel shows layer 1 (2).

In contrast to the S modes, the geometrical structure of the LB mode remains
invariant as $\theta$ changes. Figure 9 shows the LB modes at two different
angles in which the two layers have equal and opposite displacements along the
$z$-direction, similar to the LB mode in AB bilayer graphene. However, while
the LB mode in an untwisted bilayer system corresponds to a negligible in-
plane motion, the moiré LB modes correspond to swirling motion around the
corners of the moiré supercell in opposing directions between layers, which
resembles the relaxation pattern in twisted bilayer graphene [41]. This is
because the $\Gamma$ point LB mode changes the interlayer separation on the
moiré scale, which modifies the local interlayer energy. To compensate for
such an energy cost, the mode adopts an in-plane rotation that changes the
total AAarea. Consistent with Fig. 7, the relative magnitude of the LB mode
displacement field varies substantially with respect to twist angle, as shown
in the insets of Fig. 9.

Figure 10: MoSe2-WSe2 SBLB mode real space atomic displacements at
$\theta=1.00^{\circ}$. The insets and color scale are the same as in Fig. 9.

Finally, we examine the real space structure of SBLB mode in the MoSe2-WSe2
heterostructure. Unlike the LB mode where the in-plane motions are opposite at
every real space position, the SBLB mode in the heterostructure breaks layer
symmetry. As shown in Fig. 10, both layers share the same in-plane rotation
direction. Note that there is a preferential in-plane rotation direction in
the heterostructure; we did not observe a degenerate SBLB mode with the other
rotation direction. The out-of-plane motions are small but opposite between
layers.

###  IV.3 Higher energy LB (LB2 and LB3) Modes

In addition to the lowest-energy LB mode analyzed above (which we now call
LB1), we observe two LB modes of higher frequency. The key characteristics of
moiré LB modes are the oppositely oriented swirling in-plane motion and the
opposite out-of-plane motion between each layer at every point $\bm{r}$; these
are properties that all three LB modes exhibit. However, in LB1, as shown in
Fig. 9, the out-of-plane displacement has a uniform direction within each
layer. This does not hold for the LB2 and LB3 modes for all $\theta$. To
quantify the uniformity in out-of-plane displacement direction, we define the
following quantity:

$\mathcal{U}_{\text{LB}}^{12}=\frac{1}{N^{2}}\sum_{i=1}^{N^{2}}\delta_{\operatorname{sgn}[\delta\bm{u}_{z}^{(1)}(\bm{r}_{i})],1}\delta_{\operatorname{sgn}[\delta\bm{u}_{z}^{(2)}(\bm{r}_{i})],-1},$
(31)

where $N$ is the sampling grid size in the moiré supercell and
$\delta\bm{u}_{z}^{(\ell)}(\bm{r}_{i})$ is the $z$ component of the phonon at
point $\bm{r}_{i}$ in layer $\ell$. Equation (31) essentially estimates the
fraction of area where $\delta u_{z}^{(1)}$ points in the positive direction
while $\delta u_{z}^{(2)}$ points in the negative direction. Since it does not
matter which layer breathes in which out-of-plane direction, we define the
breathing uniformity value
$\mathcal{U}_{\text{LB}}=\max(\mathcal{U}_{\text{LB}}^{12},\mathcal{U}_{\text{LB}}^{21})$,
which estimates the fraction of the moiré cell that breathes in the same out-
of-plane direction. A translational mode with $z$ invariance has
$\mathcal{U}_{\text{LB}}=0$, and a perfectly uniform breathing motion (such as
that of LB1) has $\mathcal{U}_{\text{LB}}=1$. A maximally nonuniform breathing
motion has $\mathcal{U}_{\text{LB}}=0.5$, which means that at each point, the
field is equally likely to point in either direction.

Figure 11: The frequencies of the LB2 and LB3 modes color-coded by their
uniformity value $\mathcal{U}_{\text{LB}}$ (common scale on the right) for (a)
MoS2, (b) graphene, and (c) MoSe2-WSe2. LB1 in grey and SBLB in black are
shown for comparison.

Figure 11 shows the two higher-frequency LB modes for each material, color-
coded by $\mathcal{U}_{\text{LB}}$. We observe that the uniformity of the
breathing motion is dependent of $\theta$. In fact, at any value of $\theta$,
one mode is maximally uniform while the other is not. Hence, we define the LB2
mode as the maximally uniform mode and the LB3 mode as the other mode. This
identification is not well-defined in degenerate transition regions, such as
that near $\theta=6.0^{\circ}$ in MoS2 and $\theta\lesssim 1.0^{\circ}$ in all
three materials. In MoS2 only, the LB2 and LB3 mode change order in the
frequency scale.

In Fig. 12, we provide an example of a uniform versus a non-uniform higher-
frequency LB mode in MoS2 at $\theta=3.0^{\circ}$. The non-uniformity
manifests in the difference in direction of out-of-plane motion near the AA
regions versus the AB/BA regions.

Figure 12: Atomic displacement fields of LB2 and LB3 in MoS2 at
$\theta=3.0^{\circ}$. The insets and color scale are the same as in Fig. 8.

##  V Summary

We developed an accurate and efficient model for calculating phonon properties
in bilayer moiré vdW heterostructures, which is based on the configuration
space continuum (CSC) formalism. Within this formalism, we showed that phonons
related to moiré patterns in twisted bilayers can be calculated directly from
ensembles of phonons related to rigidly shifted bilayers at various
arrangementss, with first-principles accuracy.

By analyzing the low-frequency shearing (S) and layer breathing (LB) modes in
three representative materials (twisted bilayer graphene, bilayer MoS2, and a
MoSe2-WSe2 heterostructure), we showed that the physics of moiré-pattern
phonons depend on the twist angle $\theta$. Phonon frequencies vary slowly and
smoothly with $\theta$, while the atomic displacement magnitude is more
sensitive to small changes in twist angle near certain $\theta$ values. Our
CSC model can be applied more generally to classes of bilayers beyond the
three materials on which we focus here, such as bilayers that are oppositely
oriented (the pristine bilayer has one layer rotated by $180^{\circ}$ relative
to the other) or Janus materials.

Despite its efficiency and accuracy, the CSC model is limited in its present
form to lowest-frequency modes and negative shearing frequencies due to the
approximations inherent in the model. In the work presented here, we only
expand the moiré dynamical matrix to the first shell of the moiré Brillouin
zone. Increasing the cutoff radius would allow for the extension to higher
frequency modes, although a new derivation of the acoustic sum rule is
required. Future work to broaden the scope of the CSC model may include
introducing $\widetilde{\bm{k}}$ dependence on the off-diagonal components of
the moiré dynamical matrix in order to remove the low-frequency approximation,
and conducting a denser sampling in configuration space to better capture the
force field structure and thus produce non-negative shearing bands.

Building upon previous experimental studies with Raman spectroscopy [31, 32,
33, 34, 35], the phonon features derived by our CSC model may be used to guide
further experimental comparison for other materials. The most direct
application of the CSC model is as a phonon framework for a model of electron-
phonon coupling in moiré materials. Ultimately, such studies will provide a
theoretical foundation for the origin of strongly correlated states in moiré
systems, and the extent to which phonons play a role in such phenomena.

###### Acknowledgements.

We thank Lukas Linhart, Florian Libisch, and Stephen Carr for their helpful
advice. The calculations described in this paper were performed on the FASRC
Cannon cluster supported by the FAS Division of Science Research Computing
Group at Harvard University. JL acknowledges funding support from the Harvard
Herchel-Smith and PRISE fellowships. ZZ, DTL, and EK acknowledge funding from
the STC Center for Integrated Quantum Materials, NSF Grant No. DMR-1231319;
NSF DMREF Award No. 1922172; and the Army Research Office under Cooperative
Agreement Number W911NF-21-2-0147. DTL also acknowledges funding from the U.S.
Department of Energy, Office of Science, under Award Number DE-SC0019300. MA
acknowledges funding from the NSF under Award No. DMR-1922172 and the Army
Research Office under Grant Number W911NF-14-1-0247.

##  Appendix A Generalized Stacking Fault Energy (GSFE)

The GSFE functional measures the energy as a function of stacking
configuration. We evaluate the GSFE landscape $V^{\text{GSFE}}(\bm{b})$
according to a first shell expansion from Carr  _et al._ [41]. Since the
Fourier expansion is written in coefficients of $e^{i\bm{G}\cdot\bm{b}}$, with
$\bm{G}$ a reciprocal lattice vector of the bilayer at a configuration
$\bm{b}$, denote for a Bravais lattice basis
$\bm{a}_{1}=a_{0}(1,0),\bm{a}_{2}=a_{0}(1/2,\sqrt{3}/2)$:

$\displaystyle\begin{pmatrix}\mu\\\
\nu\end{pmatrix}=\frac{2\pi}{a}\begin{bmatrix}1&-1/\sqrt{3}\\\
0&2/\sqrt{3}\end{bmatrix}\bm{b}.$ (32)

Using $120^{\circ}$ rotational and $(\mu,\nu)\to(\nu,\mu)$ symmetries, the
GSFE can then be expanded to first shell as

$\displaystyle\begin{aligned}
V^{\text{GSFE}}=&\,c_{0}+c_{1}[\cos(\mu)+\cos(\nu)+\cos(\mu+\nu)]\\\
&+c_{2}[\cos(\mu+2\nu)+\cos(\mu-\nu)+\cos(2\mu+\nu)]\\\
&+c_{3}[\cos(2\mu)+\cos(2\nu)+\cos(2\mu+2\nu)]\\\
&+c_{4}[\sin(\mu)+\sin(\nu)-\sin(\mu+\nu)]\\\
&+c_{5}[\sin(2\mu+2\nu)-\sin(2\mu)-\sin(2\nu)]\end{aligned}$ (33)

with $\mathbf{c}=(c_{0},\cdots,c_{5})\in\mathbb{R}^{6}$. We fit
$V^{\text{GSFE}}$ with least-squares linear regression in the above basis. For
quantities with inversion reciprocal-space symmetry, or equivalently
$V^{\text{GSFE}}(\mu,\nu)=V^{\text{GSFE}}(-\mu,-\nu)$, $c_{4}=c_{5}=0$. The
parameters in $\mathbf{c}$ are important for computing the relaxation in
configuration space [41, 42]. Figure 13 shows the fits to the GSFE for each of
the 3 bilayers we studied. The coefficients we obtain for each material are
given in Table 1.

Figure 13: GSFE of MoS2, graphene, and MoSe2-WSe2 fitted with first-order Fourier basis regression. | MoS2 | Graphene | MoSe2-WSe2
---|---|---|---
$c_{0}$ | $25.91$ | $6.985$ | $27.07$
$c_{1}$ | $12.79$ | $4.111$ | $13.31$
$c_{2}$ | $-2.339$ | $-0.3104$ | $-2.335$
$c_{3}$ | $-8.0262$ | $-0.1023$ | $-0.8243$
$c_{4}$ | $0$ | $0$ | $0.1312$
$c_{5}$ | $0$ | $0$ | $0.2338$
Table 1: GSFE coefficients of MoS2, graphene, and MoSe2-WSe2 heterostructure.
All in units of meV per unit cell.

##  Appendix B Computational Details

We study the moiré phonons of graphene, MoS2, and the MoSe2-WSe2
heterostructure, with lattice constants $2.457,\,3.178,\,3.306$ Å,
respectively. Following the implementation in Fig. 1, we sample a $9\times 9$
uniform mesh in configuration space. For each configuration, we perform DFT
calculations using the Vienna Ab initio Simulation Package (VASP) [48, 49, 50,
51] with the r2SCAN-rVV10 van der Waals functional [52], an 800 eV plane-wave
energy cutoff, $17\times 17$ $k$-point grid, and $10^{-6}$ eV electronic
convergence threshold. The horizontal positions of the carbon atoms (graphene)
and metal atoms (TMDCs) are held fixed, but the remaining ionic coordinates
are allowed to relax until the forces are below $10^{-6}$ eV/Å. We then
compute the interatomic force constants using the frozen phonon approach with
a $3\times 3$ supercell, converging electronic self-consistent loops to
$10^{-8}$ eV. We repeat the frozen phonon calculations for the monolayers, but
use a $6\times 6$ supercell for improved accuracy.

We then construct the rigid moiré dynamical matrix in the manner described in
Section II.3 of the main text, enforcing pristine monolayer and bilayer
symmetries numerically with the phonopy and hiPhive libraries as well as the
acoustic sum rule from Appendix II.5 [53, 54]. Finally, continuum relaxation
is implemented by the procedure outlined in Section II.4 from Carr  _et al._
[41]. Continuum relaxation requires input parameters of the GSFE coefficients
from Table 1 and the monolayer bulk ($\mathcal{K}$) and shear ($\mathcal{G}$)
moduli, given in Table 2.

| MoS2 | Graphene | MoSe2 | WSe2
---|---|---|---|---
$\mathcal{K}$ | 73546 | 99265 | 81194 | 73745
$\mathcal{G}$ | 59028 | 99275 | 71562 | 63681
Table 2: Elastic constants $\mathcal{K}$ (bulk modulus) and $\mathcal{G}$
(shear modulus) in units of meV per unit cell.

Direct DFT calculations of moiré phonons were performed for a commensurate
7.34∘ graphene supercell containing 244 carbon atoms constructed from a
primitive cell with lattice constant $a=2.467$ Å. These calculations used VASP
and the PBE exchange correlation functional [55] with zero damping DFT-D3 van
der Waals corrections [56], $3\times 3\times 1$ $k$-point sampling, and a
plane wave energy cutoff of 400 eV. The atomic displacements, force constants,
and phonon band structure were computed with the phonopy package.

MD calculations of the moirè phonons are performed using the LAMMPS package
[57]. Phonons are computed by direct diagonalization of the system’s dynamical
matrix [28] obtained after relaxing the twisted bilayer’s atomic positions
[47]. The carbon-carbon intralayer interactions are modelled via the second
generation REBO potential [58]. The interlayer interactions are instead
modelled via the Kolmogorov-Crespi (KC) potential [59], using the
parametrization of Ref. [60]. The starting intralayer carbon-carbon distance
is set equal to $a_{0}=1.3978$ Å , corresponding to the equilibrium bond
length of the adopted REBO potential, giving a lattice parameter of $a\approx
2.42$ Å. Geometric optimizations are performed using the FIRE algorithm [61].
The atomic positions are relaxed toward equilibrium until the total forces
acting on each atom become less than $10^{-6}$ eV/atom.

##  Appendix C Real Space Displacement at Other High-Symmetry Points

Figure 14: Displacement fields of the first folded band in MoS2 at $M$. For
each of (a) and (b), the left (right) side displays layer 1 (2). The
periodicity is of a $2\times 2$ moiré supercell. Figure 15: Displacement
fields of the first folded band in MoS2 at $K$. For each of (a) and (b), the
left (right) side displays layer 1 (2). The maximum periodicity is of a
$3\times 3$ moiré supercell.

The real space analysis in this paper has focused on the low-energy modes at
$\Gamma$. However, the CSC model permits analysis of fields at any $k$-point,
computed in the same way but with an extra phase given by the $k$-point to the
inverse Fourier transform. Every field at $\Gamma$ has the period of one moiré
cell $\lambda$. However, at $M$, some (but not all) bands have the period of 2
moiré lengths, while at $K$ some have the period of 3 moiré lengths. This is
to be expected, as for a Bravais hexagonal lattice with basis vectors oriented
$120^{\circ}$ apart,
$M=\frac{1}{2}\widetilde{\bm{g}}_{1}+\frac{1}{2}\widetilde{\bm{g}}_{2}$ and
$K=\frac{2}{3}\widetilde{\bm{g}}_{1}+\frac{1}{3}\widetilde{\bm{g}}_{2}$ where
$\\{\widetilde{\bm{g}}_{i}\\}$ is the reciprocal lattice basis. Examples of
these periodicities are shown respectively in Figs. 14 and 15 on MoS2 on the
first folded band, which is the lowest-energy band at each $k$-point stemming
from the first folded cluster. For example, the analogous band for graphene is
the cluster between the S and LB modes of Fig. 5(a).

##  Appendix D Empirical Continuum Model

One approach to obtaining moiré phonons is through an empirical continuum
model based on the relaxation model introduced in Carr  _et al._ [41]. The
approach is introduced by Koshino and Son [29]; we briefly review and adapt it
here. To obtain the in-plane phonon modes, we perturb around the equilibrium
positions after relaxation. Here, we assume a homobilayer, for which the
relaxation in layers 1 and 2 have symmetry
$\bm{u}^{(1)}(\bm{r})=-\bm{u}^{(2)}(\bm{b})\equiv\bm{u}.$ The equilibrium
position is calculated by minimizing the total energy, which is the sum of the
interlayer energy that is the integral of the GSFE in configuration space (Eq.
(33)),

$\displaystyle E^{\mathrm{inter}}(\bm{u})=\int
V^{\mathrm{GSFE}}(\bm{u}+\bm{b})\,\differential^{2}\bm{b},$ (34)

and the intralayer elastic energy:

$\displaystyle E^{\mathrm{intra}}(\bm{u})$
$\displaystyle=2\int\Phi^{\mathrm{intra}}(\bm{u})\,\differential^{2}\bm{b},$
(35)

with

$\displaystyle\Phi^{\mathrm{intra}}(\bm{u})=\frac{1}{2}\nabla_{\bm{r}}\bm{u}(\bm{b}):C:\nabla_{\bm{r}}\bm{u}(\bm{b}),$
(36)

where $\nabla_{\bm{r}}$ is the real space gradient and $C$ is the linear
elasticity tensor. For in-plane deformations, $C$ is a rank-4 tensor with its
components defined as follows:

$\displaystyle C_{11ij}=\begin{pmatrix}\mathcal{K}+\mathcal{G}&0\\\
0&\mathcal{G}\end{pmatrix}\quad
C_{12ij}=\begin{pmatrix}0&\mathcal{K}-\mathcal{G}\\\
\mathcal{G}&0\end{pmatrix}\quad$ (37) $\displaystyle
C_{21ij}=\begin{pmatrix}0&\mathcal{G}\\\
\mathcal{K}-\mathcal{G}&0\end{pmatrix}\quad
C_{22ij}=\begin{pmatrix}\mathcal{G}&0\\\
0&\mathcal{K}+\mathcal{G}\end{pmatrix}.$

A lattice vibration perturbs the equilibrium positions, leading to an
additional kinetic energy contribution:

$T=\sum_{\ell=1}^{2}\int\frac{\rho}{2}\left[\left(\dot{U}_{x}^{(\ell)}\right)^{2}+\left(\dot{U}_{y}^{(\ell)}\right)^{2}\right]\differential^{2}\bm{r},$
(38)

where $\rho$ is the area density of a single layer ($\rho=7.61\times
10^{-7}\,\mathrm{kg/m^{2}}$ for graphene), $\ell$ labels the layer, and
$\dot{U}^{(\ell)}_{\alpha}$ is the time derivative of the cartesian component
$\alpha=x,y$ of $U$ in layer $\ell$.

The Lagrangian of the system is given by
$\mathcal{L}=T-(E^{\mathrm{intra}}+E^{\mathrm{inter}})$ as a function of
$\bm{U}^{(\ell)}(\bm{r})$. We define
$\bm{U}^{\pm}=\bm{U}^{(2)}\pm\bm{U}^{(1)}$ and rewrite $\mathcal{L}$ as a
function of $\bm{U}^{\pm}$. Let us rewrite the GSFE from Eq. (33) in a more
compact form:

$\displaystyle
V^{\mathrm{GSFE}}(\bm{b})=c_{0}+\sum_{j=1}^{3}\sum_{n=1}^{3}c_{j}\cos(\bm{G}_{j}^{n}\cdot\bm{b}),$
(39)

where $\bm{G}_{j}^{n}$ are the monolayer reciprocal lattice vectors defined in
Fig. 16. Here, the subscript $j$ labels the shell (vectors that are the same
distance from the origin) and the superscript $n$ labels the element within
the shell. For simplicity we assume symmetry between AB and BA stacking and
only consider the cosine components of the GSFE here ($c_{1},c_{2},c_{3}$),
but the following results can be easily generalized to include the sine
components ($c_{4},c_{5},c_{6}$).

The equations of motion for $\bm{U}^{-}$ are given as follows:

$\displaystyle\begin{aligned}
\frac{1}{2}\rho\ddot{U}_{x}^{-}=\frac{1}{2}\left(\mathcal{K}+\frac{\mathcal{G}}{3}\right)\left(\frac{\partial^{2}U_{x}^{-}}{\partial
x^{2}}+\frac{\partial^{2}{U_{y}^{-}}}{\partial{x}\partial{y}}\right)+\frac{\mathcal{G}}{2}\left(\frac{\partial^{2}{U_{x}^{-}}}{\partial{x}^{2}}+\frac{\partial^{2}{U_{x}^{-}}}{\partial{y}^{2}}\right)+\sum_{j=1}^{3}\sum_{n=1}^{3}G_{j,x}^{n}c_{j}\sin(\bm{G}_{j}^{n}\cdot(\bm{b}+\bm{u^{-}}))\\\
\frac{1}{2}\rho\ddot{U}_{y}^{-}=\frac{1}{2}\left(\mathcal{K}+\frac{\mathcal{G}}{3}\right)\left(\frac{\partial^{2}U_{y}^{-}}{\partial
y^{2}}+\frac{\partial^{2}{U_{x}^{-}}}{\partial{x}\partial{y}}\right)+\frac{\mathcal{G}}{2}\left(\frac{\partial^{2}{U_{y}^{-}}}{\partial{x}^{2}}+\frac{\partial^{2}{U_{y}^{-}}}{\partial{y}^{2}}\right)+\sum_{j=1}^{3}\sum_{n=1}^{3}G_{j,y}^{n}c_{j}\sin(\bm{G}_{j}^{n}\cdot(\bm{b}+\bm{u^{-}})),\end{aligned}$
(40)

where $G_{j,\alpha}^{n}$ is the $\alpha=x,y$ component of $\bm{G}_{j}^{n}$.
Note that $\bm{U}(\bm{r})$ is the real space displacement and $\bm{u}(\bm{b})$
is the configuration space displacement and they are related by the linear
transformation in Eq. (24). The terms on the 2nd lines of Eq. (40) that are
proportional to $c_{j}$ come from the GSFE, which depends only on
$\bm{U}^{-}$, the difference between relaxation in layer 1 and layer 2, and
not on $\bm{U}^{+}$. Thus the equations of motion for $\bm{U}^{+}$ are given
by replacing $\bm{U}^{-}$ with $\bm{U}^{+}$ and dropping the terms
proportional to $c_{j}$.

Figure 16: Definition of reciprocal lattice vectors in Eq. (39). where the
subscript labels different shells (different colors) and the superscript
labels different components within the same shell (same color). The gray
scatter points are reciprocal lattice vectors.

Because $\bm{U}^{+}$ represents equal motion in both layers, its fluctuations
yield the acoustic phonons of a monolayer.

To determine the moiré phonons we study $\bm{U}^{-}$ and consider a small
vibration around the static equilibrium relaxation:

$\displaystyle\bm{U}^{-}(\bm{r},t)=\bm{U}_{0}^{-}(\bm{r})+\delta\bm{U}^{-}(\bm{r},t),$
(41)

where $\bm{U}_{0}^{-}$ is the static relaxation that minimizes the total
energy, $E^{\mathrm{intra}}+E^{\mathrm{inter}}$, and
$\delta\bm{U}^{-}(\bm{r},t)$ is the perturbation around $\bm{U}_{0}^{-}$,
which constitutes the phonon modes. We define the following Fourier
coefficients to express the displacement and equations of motion in Fourier
space:

$\displaystyle\bm{U}_{0}^{-}(\bm{r})=\sum_{\bm{G}}\bm{U}_{0,\widetilde{\bm{G}}}e^{i\widetilde{\bm{G}}\cdot\bm{r}},$
$\displaystyle\delta\bm{U}^{-}(\bm{r},t)=e^{-i\omega
t}\sum_{\bm{q}}\delta\bm{U}^{-}_{\bm{q}}e^{i\bm{q}\cdot\bm{r}},$
$\displaystyle\sin(\bm{G}_{j}^{n}\cdot(\bm{b}+\bm{u^{-}}))=\sum_{\widetilde{\bm{G}}}f_{j,\bm{G}}^{n}e^{i\widetilde{\bm{G}}\cdot\bm{r}},$
$\displaystyle\cos(\bm{G}_{j}^{n}\cdot(\bm{b}+\bm{u^{-}}))=\sum_{\widetilde{\bm{G}}}h_{j,\bm{G}}^{n}e^{i\widetilde{\bm{G}}\cdot\bm{r}},$
(42)

where $\widetilde{\bm{G}}$ is the moiré reciprocal lattice vector. Note that
$\widetilde{\bm{G}}\cdot\bm{r}=\bm{G}\cdot\bm{b}$, where $\bm{G}$ is the
monolayer reciprocal lattice vector. We expand Eq. (40) around the static
solution $\bm{U}_{0}$ to the first order in $\delta\bm{U}$ and perform a
Fourier expansion to obtain the equation of motion for the perturbation part:

$\rho\omega^{2}\delta\bm{U}^{-}_{\widetilde{\bm{G}}+\bm{q}}=\widehat{K}_{\widetilde{\bm{G}}+\bm{q}}\delta\bm{U}^{-}_{\widetilde{\bm{G}}+\bm{q}}-4\sum_{\widetilde{\bm{G}^{\prime}}}\widehat{V}_{\widetilde{\bm{G}}-\widetilde{\bm{G}}^{\prime}}\delta\bm{U}^{-}_{\widetilde{\bm{G}^{\prime}}+\bm{q}},$
(43)

where

$\widehat{K}_{\widetilde{\bm{G}}}=\begin{pmatrix}(\mathcal{K}+\frac{4}{3}\mathcal{G})G_{x}^{2}+\mathcal{G}G_{y}^{2}&(\mathcal{K}+\frac{\mathcal{G}}{3})G_{x}G_{y}\\\
(\mathcal{K}+\frac{\mathcal{G}}{3})G_{x}G_{y}&(\mathcal{K}+\frac{4}{3}\mathcal{G})G_{y}^{2}+\mathcal{G}G_{x}^{2}\end{pmatrix}$
(44)

and

$\widehat{V}_{\widetilde{\bm{G}}}=\sum_{j=1}^{3}c_{j}\sum_{n=1}^{3}h_{j,\widetilde{\bm{G}}}^{n}\begin{pmatrix}G_{j,x}^{n}G_{j,x}^{n}&G_{j,x}^{n}G_{j,y}^{n}\\\
G_{j,x}^{n}G_{j,y}^{n}&G_{j,y}^{n}G_{j,y}^{n}\end{pmatrix}.$ (45)

Here, the $\widehat{V}_{\widetilde{\bm{G}}-\widetilde{\bm{G}}^{\prime}}$ is
essentially an interlayer scattering selection rule that couples
$\widetilde{\bm{G}}$ and $\widetilde{\bm{G}}^{\prime}$ After obtaining the
equilibrium relaxation displacement vectors $\bm{U}_{0}^{-}$ self-consistently
[41], we solve the eigenvalue problem in Eq. (43) to obtain the eigenmodes
$\omega^{2}$.

## References

  * Carr  _et al._ [2017] S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle, Phys. Rev. B 95, 075420 (2017), arXiv:1611.00649 [cond-mat.mes-hall] .
  * Yankowitz  _et al._ [2018] M. Yankowitz, J. Jung, E. Laksono, N. Leconte, B. L. Chittari, K. Watanabe, T. Taniguchi, S. Adam, D. Graf, and C. R. Dean, Dynamic band-structure tuning of graphene moiré superlattices with pressure, Nature (London) 557, 404 (2018), arXiv:1707.09054 [cond-mat.mes-hall] .
  * Liu  _et al._ [2019] X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Spin-polarized Correlated Insulator and Superconductor in Twisted Double Bilayer Graphene, arXiv e-prints , arXiv:1903.08130 (2019), arXiv:1903.08130 [cond-mat.mes-hall] .
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, All Magic Angles in Twisted Bilayer Graphene are Topological, Phys. Rev. Lett. 123, 036401 (2019), arXiv:1807.10676 [cond-mat.mes-hall] .
  * Burg  _et al._ [2019] G. W. Burg, J. Zhu, T. Taniguchi, K. Watanabe, A. H. MacDonald, and E. Tutuc, Correlated Insulating States in Twisted Double Bilayer Graphene, Phys. Rev. Lett. 123, 197702 (2019), arXiv:1907.10106 [cond-mat.mes-hall] .
  * Cao  _et al._ [2020] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable correlated states and spin-polarized phases in twisted bilayer-bilayer graphene, Nature (London) 583, 215 (2020).
  * Moriyama  _et al._ [2019] S. Moriyama, Y. Morita, K. Komatsu, K. Endo, T. Iwasaki, S. Nakaharai, Y. Noguchi, Y. Wakayama, E. Watanabe, D. Tsuya, K. Watanabe, and T. Taniguchi, Observation of superconductivity in bilayer graphene/hexagonal boron nitride superlattices, arXiv e-prints , arXiv:1901.09356 (2019), arXiv:1901.09356 [cond-mat.supr-con] .
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature (London) 556, 43 (2018a), arXiv:1803.02342 [cond-mat.mes-hall] .
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature (London) 556, 80 (2018b), arXiv:1802.00553 [cond-mat.mes-hall] .
  * 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, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature (London) 574, 653 (2019), arXiv:1903.06513 [cond-mat.str-el] .
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the Correlated Insulator States in Twisted Bilayer Graphene, Phys. Rev. Lett. 124, 097601 (2020), arXiv:2109.13528 [cond-mat.mes-hall] .
  * Wang  _et al._ [2019] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. Pasupathy, and C. R. Dean, Magic continuum in twisted bilayer WSe2, arXiv e-prints , arXiv:1910.12147 (2019), arXiv:1910.12147 [cond-mat.mes-hall] .
  * Li  _et al._ [2019] F. Li, W. Wei, H. Wang, B. Huang, Y. Dai, and T. Jacob, Intrinsic electric field-induced properties in janus mosse van der waals structures, The Journal of Physical Chemistry Letters 10, 559 (2019), https://doi.org/10.1021/acs.jpclett.8b03463 .
  * Wu  _et al._ [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard Model Physics in Transition Metal Dichalcogenide Moiré Bands, Phys. Rev. Lett. 121, 026402 (2018), arXiv:1804.03151 [cond-mat.mes-hall] .
  * Liu  _et al._ [2014] K. Liu, L. Zhang, T. Cao, C. Jin, D. Qiu, Q. Zhou, A. Zettl, P. Yang, S. G. Louie, and F. Wang, Evolution of interlayer coupling in twisted molybdenum disulfide bilayers, Nature Communications 5, 4966 (2014), arXiv:1406.6487 [cond-mat.mes-hall] .
  * Wu  _et al._ [2019] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Topological Insulators in Twisted Transition Metal Dichalcogenide Homobilayers, Phys. Rev. Lett. 122, 086402 (2019), arXiv:1807.03311 [cond-mat.mes-hall] .
  * Zhang  _et al._ [2021] K. Zhang, Y. Guo, D. T. Larson, Z. Zhu, S. Fang, E. Kaxiras, J. Kong, and S. Huang, Spectroscopic Signatures of Interlayer Coupling in Janus MoSSe/MoS2 Heterostructures, arXiv e-prints , arXiv:2108.10995 (2021), arXiv:2108.10995 [cond-mat.mtrl-sci] .
  * Tran  _et al._ [2019] K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, J. Embley, A. Zepeda, M. Campbell, T. Autry, T. Taniguchi, K. Watanabe, N. Lu, S. K. Banerjee, K. L. Silverman, S. Kim, E. Tutuc, L. Yang, A. H. MacDonald, and X. Li, Evidence for moiré excitons in van der Waals heterostructures, Nature (London) 567, 71 (2019), arXiv:1807.03771 [cond-mat.mes-hall] .
  * Jin  _et al._ [2019] C. Jin, E. C. Regan, A. Yan, M. Iqbal Bakti Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng, S. Shi, K. Watanabe, T. Taniguchi, S. Tongay, A. Zettl, and F. Wang, Observation of moiré excitons in WSe2/WS2 heterostructure superlattices, Nature (London) 567, 76 (2019), arXiv:1812.09815 [cond-mat.mes-hall] .
  * Alexeev  _et al._ [2019] E. M. Alexeev, D. A. Ruiz-Tijerina, M. Danovich, M. J. Hamer, D. J. Terry, P. K. Nayak, S. Ahn, S. Pak, J. Lee, J. I. Sohn, M. R. Molas, M. Koperski, K. Watanabe, T. Taniguchi, K. S. Novoselov, R. V. Gorbachev, H. S. Shin, V. I. Fal’ko, and A. I. Tartakovskii, Resonantly hybridized excitons in moiré superlattices in van der Waals heterostructures, Nature (London) 567, 81 (2019), arXiv:1904.06214 [cond-mat.mes-hall] .
  * Seyler  _et al._ [2019] K. L. Seyler, P. Rivera, H. Yu, N. P. Wilson, E. L. Ray, D. G. Mandrus, J. Yan, W. Yao, and X. Xu, Signatures of moiré-trapped valley excitons in MoSe2/WSe2 heterobilayers, Nature (London) 567, 66 (2019), arXiv:1809.04562 [cond-mat.mes-hall] .
  * Rivera  _et al._ [2015] P. Rivera, J. R. Schaibley, A. M. Jones, J. S. Ross, S. Wu, G. Aivazian, P. Klement, K. Seyler, G. Clark, N. J. Ghimire, J. Yan, D. G. Mandrus, W. Yao, and X. Xu, Observation of long-lived interlayer excitons in monolayer MoSe2-WSe2 heterostructures, Nature Communications 6, 6242 (2015), arXiv:1403.4985 [cond-mat.mes-hall] .
  * van der Zande  _et al._ [2014] A. M. van der Zande, J. Kunstmann, A. Chernikov, D. A. Chenet, Y. You, X. Zhang, P. Y. Huang, T. C. Berkelbach, L. Wang, F. Zhang, M. S. Hybertsen, D. A. Muller, D. R. Reichman, T. F. Heinz, and J. C. Hone, Tailoring the Electronic Structure in Bilayer Molybdenum Disulfide via Interlayer Twist, Nano Letters 14, 3869 (2014).
  * Angeli and MacDonald [2021] M. Angeli and A. H. MacDonald, $\Gamma$ valley transition metal dichalcogenide moiré bands, Proceedings of the National Academy of Science 118, e2021826118 (2021).
  * Wu  _et al._ [2018] F. Wu, A. H. MacDonald, and I. Martin, Theory of phonon-mediated superconductivity in twisted bilayer graphene, Phys. Rev. Lett. 121, 257001 (2018).
  * Lian  _et al._ [2019] B. Lian, Z. Wang, and B. A. Bernevig, Twisted bilayer graphene: A phonon-driven superconductor, Phys. Rev. Lett. 122, 257002 (2019).
  * Cea and Guinea [2021] T. Cea and F. Guinea, Coulomb interaction, phonons, and superconductivity in twisted bilayer graphene, Proceedings of the National Academy of Sciences 118, e2107874118 (2021), https://www.pnas.org/doi/pdf/10.1073/pnas.2107874118 .
  * Angeli  _et al._ [2019] M. Angeli, E. Tosatti, and M. Fabrizio, Valley jahn-teller effect in twisted bilayer graphene, Phys. Rev. X 9, 041010 (2019).
  * Koshino and Son [2019] M. Koshino and Y.-W. Son, Moiré phonons in twisted bilayer graphene, Phys. Rev. B 100, 075416 (2019), arXiv:1905.09660 [cond-mat.mes-hall] .
  * Gadelha  _et al._ [2021] A. C. Gadelha, D. A. A. Ohlberg, C. Rabelo, E. G. S. Neto, T. L. Vasconcelos, J. L. Campos, J. S. Lemos, V. Ornelas, D. Miranda, R. Nadas, F. C. Santana, K. Watanabe, T. Taniguchi, B. van Troeye, M. Lamparski, V. Meunier, V.-H. Nguyen, D. Paszko, J.-C. Charlier, L. C. Campos, L. G. Cançado, G. Medeiros-Ribeiro, and A. Jorio, Localization of lattice dynamics in low-angle twisted bilayer graphene, Nature (London) 590, 405 (2021).
  * Campos-Delgado  _et al._ [2013] J. Campos-Delgado, L. G. Cançado, C. A. Achete, A. Jorio, and J. P. Raskin, Raman-scattering study of the phonon dispersion in twisted bi-layer graphene, arXiv e-prints , arXiv:1301.3795 (2013), arXiv:1301.3795 [cond-mat.mes-hall] .
  * Wu  _et al._ [2014] J.-B. Wu, X. Zhang, M. Ijäs, W.-P. Han, X.-F. Qiao, X.-L. Li, D.-S. Jiang, A. C. Ferrari, and P.-H. Tan, Resonant Raman spectroscopy of twisted multilayer graphene, Nature Communications 5, 5309 (2014), arXiv:1408.6017 [cond-mat.mtrl-sci] .
  * Huang  _et al._ [2016] S. Huang, L. Liang, X. Ling, A. A. Puretzky, D. B. Geohegan, B. G. Sumpter, J. Kong, V. Meunier, and M. S. Dresselhaus, Low-Frequency Interlayer Raman Modes to Probe Interface of Twisted Bilayer MoS2, Nano Letters 16, 1435 (2016).
  * Lin  _et al._ [2018] M.-L. Lin, Q.-H. Tan, J.-B. Wu, X.-S. Chen, J.-H. Wang, Y.-H. Pan, X. Zhang, X. Cong, J. Zhang, W. Ji, P.-A. Hu, K.-H. Liu, and P.-H. Tan, Moiré phonons in twisted bilayer MoS2, ACS Nano 12, 8770 (2018), pMID: 30086224, https://doi.org/10.1021/acsnano.8b05006 .
  * Quan  _et al._ [2021] J. Quan, L. Linhart, M.-L. Lin, D. Lee, J. Zhu, C.-Y. Wang, W.-T. Hsu, J. Choi, J. Embley, C. Young, T. Taniguchi, K. Watanabe, C.-K. Shih, K. Lai, A. H. MacDonald, P.-H. Tan, F. Libisch, and X. Li, Phonon renormalization in reconstructed MoS2 moiré superlattices, Nature Materials 20, 1100 (2021), arXiv:2009.10650 [cond-mat.mes-hall] .
  * Dai  _et al._ [2016] S. Dai, Y. Xiang, and D. J. Srolovitz, Twisted Bilayer Graphene: Moiré with a Twist, Nano Letters 16, 5923 (2016).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Lattice relaxation and energy band modulation in twisted bilayer graphene, Phys. Rev. B 96, 075311 (2017), arXiv:1706.03908 [cond-mat.mtrl-sci] .
  * Español  _et al._ [2017] M. I. Español, D. Golovaty, and J. P. Wilber, Discrete-to-continuum modeling of weakly interacting incommensurate chains, Phys. Rev. E 96, 033003 (2017), arXiv:1705.00072 [cond-mat.mtrl-sci] .
  * Español  _et al._ [2018] M. I. Español, D. Golovaty, and J. P. Wilber, Discrete-to-continuum modelling of weakly interacting incommensurate two-dimensional lattices, Proceedings of the Royal Society of London Series A 474, 20170612 (2018), arXiv:1708.00974 [cond-mat.mtrl-sci] .
  * Leconte  _et al._ [2019] N. Leconte, S. Javvaji, J. An, A. Samudrala, and J. Jung, Relaxation Effects in Twisted Bilayer Graphene: a Multi-Scale Approach, arXiv e-prints , arXiv:1910.12805 (2019), arXiv:1910.12805 [cond-mat.mes-hall] .
  * Carr  _et al._ [2018] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras, Relaxation and domain formation in incommensurate two-dimensional heterostructures, Phys. Rev. B 98, 224102 (2018), arXiv:1805.06972 [cond-mat.mes-hall] .
  * Zhu  _et al._ [2020] Z. Zhu, P. Cazeaux, M. Luskin, and E. Kaxiras, Modeling mechanical relaxation in incommensurate trilayer van der Waals heterostructures, Phys. Rev. B 101, 224107 (2020), arXiv:1911.05324 [cond-mat.mes-hall] .
  * Massatt  _et al._ [2016] D. Massatt, M. Luskin, and C. Ortner, Electronic Density of States for Incommensurate Layers, arXiv e-prints , arXiv:1608.01968 (2016), arXiv:1608.01968 [math-ph] .
  * Cazeaux  _et al._ [2016] P. Cazeaux, M. Luskin, and E. B. Tadmor, Analysis of rippling in incommensurate one-dimensional coupled chains, arXiv e-prints , arXiv:1606.02619 (2016), arXiv:1606.02619 [math.NA] .
  * Cancès  _et al._ [2017] E. Cancès, P. Cazeaux, and M. Luskin, Generalized Kubo formulas for the transport properties of incommensurate 2D atomic heterostructures, Journal of Mathematical Physics 58, 063502 (2017), arXiv:1611.08043 [math-ph] .
  * Cazeaux  _et al._ [2019] P. Cazeaux, M. Luskin, and D. Massatt, Energy Minimization of Two Dimensional Incommensurate Heterostructures, Archive for Rational Mechanics and Analysis 235, 1289 (2019), arXiv:1806.10395 [physics.comp-ph] .
  * Angeli  _et al._ [2018] M. Angeli, D. Mandelli, A. Valli, A. Amaricci, M. Capone, E. Tosatti, and M. Fabrizio, Emergent ${D}_{6}$ symmetry in fully relaxed magic-angle twisted bilayer graphene, Phys. Rev. B 98, 235137 (2018).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal-amorphous-semiconductor transition in germanium, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Furness  _et al._ [2020] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, Accurate and numerically efficient r2SCAN meta-generalized gradient approximation, arXiv e-prints , arXiv:2008.03374 (2020), arXiv:2008.03374 [cond-mat.mtrl-sci] .
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, arXiv e-prints , arXiv:1506.08498 (2015), arXiv:1506.08498 [cond-mat.mtrl-sci] .
  * Eriksson  _et al._ [2018] F. Eriksson, E. Fransson, and P. Erhart, The hiphive package for the extraction of high-order force constants by machine learning, arXiv e-prints , arXiv:1811.09267 (2018), arXiv:1811.09267 [cond-mat.mtrl-sci] .
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, The Journal of Chemical Physics 132, 154104 (2010), https://doi.org/10.1063/1.3382344 .
  * 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, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications 271, 108171 (2022).
  * Brenner  _et al._ [2002] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons, Journal of Physics Condensed Matter 14, 783 (2002).
  * Kolmogorov and Crespi [2005] A. N. Kolmogorov and V. H. Crespi, Registry-dependent interlayer potential for graphitic systems, Phys. Rev. B 71, 235415 (2005).
  * Ouyang  _et al._ [2018] W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod, Nanoserpents: Graphene nanoribbon motion on two-dimensional hexagonal materials, Nano Letters 18, 6009 (2018), pMID: 30109806, https://doi.org/10.1021/acs.nanolett.8b02848 .
  * Bitzek  _et al._ [2006] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Structural relaxation made simple, Phys. Rev. Lett. 97, 170201 (2006).