# Moiré of Moiré: Modeling Mechanical Relaxation in Incommensurate Trilayer
van der Waals Heterostructures

Ziyan Zhu  Department of Physics, Harvard University, Cambridge, Massachusetts
02138, USA  Paul Cazeaux  Department of Mathematics, University of Kansas,
Lawrence, Kansas 66045, USA  Mitchell Luskin  School of Mathematics,
University of Minnesota, Minneapolis, Minnesota 55455, 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

The incommensurate stacking of multi-layered two-dimensional materials is a
challenging problem from a theoretical perspective and an intriguing avenue
for manipulating their physical properties. Here we present a multi-scale
model to obtain the mechanical relaxation pattern of twisted trilayer van der
Waals (vdW) heterostructures with two independent twist angles, a generally
incommensurate system without a supercell description. We adopt the
configuration space as a natural description of such incommensurate layered
materials, based on the local environment of atomic positions, bypassing the
need for commensurate approximations. To obtain the relaxation pattern, we
perform energy minimization with respect to the relaxation displacement
vectors. We use a continuum model in combination with the Generalized Stacking
Fault energy to describe the interlayer coupling, obtained from first-
principles calculations based on Density Functional Theory. We show that the
relaxation patterns of twisted trilayer graphene and $\mathrm{WSe_{2}}$ are
“moiré of moiré”, as a result of the incommensurate coupling two bilayer moiré
patterns. We also show that, in contrast to the symmetry-preserving in-plane
relaxation in twisted bilayers, trilayer relaxation can break the two-fold
rotational symmetry about the xy-plane when the two twist angles are equal.

††preprint: APS/123-QED

##  I Introduction

Assemblies of multilayers of two-dimensional (2D) materials, referred to as
van der Waals (vdW) heterostructures, have become a favorite platform for
exploring strongly correlated states, because they offer a large parameter
space and great tunability. Superconductivity and strongly-correlated states
have been discovered in these heterostructures with one small twist angle
(e.g., twisted bilayer graphene, twisted double-bilayer graphene) Cao  _et
al._ (2018a, b); Yankowitz  _et al._ (2018); Liu  _et al._ (2019); Cao  _et
al._ (2019); Burg  _et al._ (2019); Shen  _et al._ (2019); Lu  _et al._ (2019)
and systems with a small misalignment (e.g., hBN + graphene) Chen  _et al._
(2019); Moriyama  _et al._ (2019). Due to the twist angle or the small lattice
misalignment, the system forms moiré patterns and the periodicity becomes much
larger than the original unit cell size. To understand the electronic and
mechanical properties of the multilayered structures, we view them as a series
of conventional crystals with a weak perturbative interaction between layers
Tritsaris  _et al._ (2016).

Crystal relaxation in vdW heterostructures has been studied in the continuum
limit by energy minimization San-Jose  _et al._ (2014); Dai  _et al._ (2016);
Nam and Koshino (2017); Español  _et al._ (2017, 2018); Carr  _et al._ (2018);
Leconte  _et al._ (2019). Relaxation in incommensurate stacked 2D materials
has been known to form domain walls separating large commensurate areas
Coppersmith  _et al._ (1982). For example, in twisted bilayer graphene,
relaxation enlarges the AB/BA stacking regions (the equilibrium stacking) and
forms a thin domain line in between adjacent AB/BA stacking regions, while the
AA stacking regions (high energy stacking order) shrink to localized spots.
The domain formation in twisted bilayer systems has been observed
experimentally Alden  _et al._ (2013); Woods  _et al._ (2014); van Wijk  _et
al._ (2015); Yankowitz  _et al._ (2016); Kim  _et al._ (2017); Yoo  _et al._
(2019). The relaxation pattern has a dependence on the twist angle: the
magnitude of relaxation displacement vectors decreases as the twist angle
increases. At a twist angle $\theta>\theta^{*}$ ($\theta^{*}\sim 1^{\circ}$ in
twisted bilayer graphene), the relaxation effects are weak. As the twist angle
becomes smaller, the domain wall thickness stays roughly constant and only the
AB/BA stacking regions get enlarged.

A natural extension is to add another layer, which introduces a new degree of
freedom — a second twist angle — to tune the system. We emphasize that the
system we discuss here is different than multi-layered systems such as double-
bilayer graphene, in which there is only one independent parameter or one
twist angle despite the presence of four layers. The experimental technique to
realize twisted trilayer graphene is readily available and unconventional
correlated states have been observed Tsai  _et al._ (2019). There have been
preliminary studies of electronic states of twisted trilayer graphene using
$k\cdot p$ perturbation theory Mora  _et al._ (2019); Ma  _et al._ (2019);
Tsai  _et al._ (2019). In these models, mechanical relaxation has not been
taken into consideration. However, atomic relaxation will modify the band
structure by opening up a single-particle gap in twisted bilayer graphene Nam
and Koshino (2017); Carr  _et al._ (2019a); Guinea and Walet (2019) and thus
is an important effect in obtaining realistic description of electronic
behaviors.

The extension from twisted bilayers to twisted trilayers is challenging for
the following reasons. (1) The system size becomes much larger. In twisted
bilayer systems the relevant length scale is on the order of $1/\theta$. In
twisted trilayers, assuming the two twist angles are identical, the length
scale becomes $1/\theta^{2}$, where $\theta$ is the relevant twist angle,
typically a small quantity for systems of interest. (2) Because of the
independent second twist angle, a supercell description is no longer valid,
except for specific cases. The trilayer system is truly incommensurate, while
there always exists a moiré supercell in the continuum limit or the small
angle limit of twisted bilayers.

To address the challenge of the large system size and the lack of periodicity,
we introduce a multiscale method based on the concept of configuration space
to minimize the energy directly in the large system limit without having to
use a supercell approximation. Specifically, we perform energy minimization in
the continuum approximation over a collection of local atomic environments
with respect to the other layers, which we call the configuration space
Massatt  _et al._ (2016); Cazeaux  _et al._ (2016); Cancès  _et al._ (2017);
Cazeaux  _et al._ (2019). This approach has recently been applied to
describing the crystal relaxation in stacked bilayer materials Carr  _et al._
(2018); Cazeaux  _et al._ (2019).

The rest of the paper is organized as follows: in Section II, we introduce the
configuration space to describe the trilayer system and formulate the energy
minimization problem in configuration space using a continuum approximation.
In Section III, we present the results for twisted trilayer graphene and
$\mathrm{WSe_{2}}$, as representative heterostructures of a 2D semi-metal and
a 2D semiconductor. Finally, we summarize the results and discuss potential
applications in Section IV.

Figure 1: (a) An example of a twisted trilayer honeycomb lattice in real space
with $\theta_{12}=5.3^{\circ}$ and $\theta_{23}=7.7^{\circ}$. (b) Magnified
view at the black box marked in (a). The twist angle between $L_{1}$ and
$L_{2}$, $\theta_{12}$, and the twist angle between $L_{2}$ and $L_{3}$,
$\theta_{23},$ are marked by black arrows.

##  II Continuum Approximation and Energy Minimization in Configuration Space

###  II.1 Setup and Notation

We first define the Bravais lattice basis vectors of a monolayer, $\bm{a}_{1}$
and $\bm{a}_{2}$ respectively, in a $2\times 2$ matrix form as the columns of
the following matrix:

$\displaystyle A=a_{0}\begin{pmatrix}\sqrt{3}/2&\sqrt{3}/2\\\
-1/2&1/2\end{pmatrix}=(\bm{a}_{1}\quad\bm{a}_{2}),$ (1)

where $a_{0}$ is the lattice constant. We denote the unit cell of layer $i$ by
$\Gamma^{(i)}.$

We consider the twisted trilayer systems with two independent twist angles; an
example of such a system is given in Fig. 1a. Suppose the second layer $L_{2}$
is fixed at the origin. The first layer ($L_{1}$) is twisted counterclockwise
by $\theta_{12}$ with respect to $L_{2}$, and $L_{3}$ is twisted clockwise by
$\theta_{23}$ with respect to $L_{2}$ (see Fig. 1b). The primitive Bravais
lattice basis vectors in the matrix form of $L_{i}$ is defined as $A_{i}$.
With $A_{2}=A$, $A_{1}$ and $A_{3}$ are given as follows,

$\displaystyle A_{1}=$
$\displaystyle\begin{pmatrix}\cos{\theta_{12}}&-\sin{\theta_{12}}\\\
\sin{\theta_{12}}&\cos{\theta_{12}}\end{pmatrix}A=(\bm{a}_{1}^{(1)}\quad\bm{a}_{2}^{(1)}),$
$\displaystyle A_{3}=$
$\displaystyle\begin{pmatrix}\cos{\theta_{23}}&\sin{\theta_{23}}\\\
-\sin{\theta_{23}}&\cos{\theta_{23}}\end{pmatrix}A=(\bm{a}_{1}^{(3)}\quad\bm{a}_{2}^{(3)}).$
(2)

The primitive reciprocal lattice vectors are given by the columns of $2\pi
A_{i}^{-T}.$

Figure 2: The definition of the local shift or the disregistry $\bm{b}^{(j)}$.
Here, we take $j=2$ as an example, defining $\bm{b}^{(2)}$ between $L_{1}$ and
$L_{2}$.

Figure 3: Mapping between real and configuration space. Left: Real space
atomic positions of a trilayer triangular lattice. The first layer (blue) is
twisted by $2^{\circ}$ with respect to the second layer (red), and the third
layer (green) is twisted by $3^{\circ}$ with respect to the second. Middle:
magnified view of the three boxed area on the left. Red parallelograms are the
unit cell of the layer 2. Arrows are pointing from the atom of interest (the
red atom on the left corner) to the neighboring atoms on the other two layers.
Right: the corresponding configuration space to the atom of interest. Blue,
red, and green parallelograms are the unit cells of $L_{i}$ for $i=1,2,3$
respectively. Figure 4: A demonstration of the bilinear interpolation of the
local shift vector $\bm{b}^{2\rightarrow 3}$. This corresponds to the same
real space position in the green (the second) box in Fig. 3 but with a
different origin. Here, the origin is chosen to be a lattice site of $L_{1}$
with $\bm{b}^{(1)}=\bm{0}$, which is marked by a blue cross in the figure,
while in Fig. 3, the origin is the red atom on the left corner.
$\bm{b}_{mn}^{2\rightarrow 3}$ are the relative shifts between $L_{2}$ and
$L_{3}$ at the four nearest lattice sites to the origin.

###  II.2 Relaxation in Configuration Space

The issue with applying the real space continuum model to obtain the
relaxation pattern in trilayer systems is the lack of periodicity. While
bilayer systems always have a moiré supercell Carr  _et al._ (2018); Cazeaux
_et al._ (2019), trilayer systems are generally incommensurate and thus lack a
periodic supercell even in the continuum limit. Therefore, a more general
description beyond the supercell approximation is required. Here, we introduce
the configuration space to describe the local environment of every position in
the continuum and parametrize the trilayer system in configuration space
Massatt  _et al._ (2016); Cazeaux  _et al._ (2016); Cancès  _et al._ (2017)
based on the formalism introduced by Cazeaux et al. Cazeaux  _et al._ (2019).
We also reformulate the energy minimization problem in configuration space.

The local configuration or the environment for an arbitrary position $\bm{r}$
in layer $L_{i}$, which may be a lattice position on layer $L_{i}$, can be
uniquely determined by three relative shift vectors or disregistries
$\bm{b}^{(j)}$ describing the relative position from this position in $L_{i}$
with respect to the lattice of layer $L_{j}$ (see Fig. 2 for an example),
which is defined as the lattice position with respect to $\bm{r}$ modulo the
unit cell $\Gamma^{(j)}$. More explicitly,

$\bm{b}^{(j)}(\bm{r})=\bm{R}^{(j)}-\bm{r}\in\Gamma^{(j)},$ (3)

where $\bm{R}^{(j)}$ is any lattice site in layer $L_{j}$. The shift vectors
$\bm{b}^{(j)}$ take values in the periodic unit cell $\Gamma^{(j)}$ of layer
$L_{j}$ by the translational invariance of the lattice. Fig. 3 illustrates
this mapping between real space and the configuration space of the three
relative shift vectors for $\bm{r}$ sampled at lattice points of the middle
layer $L_{2}$ of a triangular lattice trilayer for a representative atomic
structure. The collection of these shift vectors forms a six-dimensional
configuration space:

$\Omega=\\{\omega\mid\omega=(\bm{b}^{(1)},\bm{b}^{(2)},\bm{b}^{(3)})\\}.$ (4)

The shifts $\bm{b}^{(j)}(\bm{r})$ vary quickly on the scale of the lattice
size, and thus do not correspond to the intuitive notion of disregistry which
varies on the scale of the moiré pattern. Instead, we construct functions
$\bm{b}^{i\rightarrow j}(\omega)$ measuring smoothly the disregistry between
layer $L_{i}$ and $L_{j}$ as a function of the local configuration $\omega$,
by using bilinear interpolation between lattice points Carr  _et al._ (2018);
Cazeaux  _et al._ (2019)

$\bm{b}^{i\rightarrow
j}(\omega)=\bm{b}^{(j)}-A_{j}A_{i}^{-1}\bm{b}^{(i)}\in\Gamma^{(j)}.$ (5)

This expression is periodic as a function of $\bm{b}^{(i)}$ and $\bm{b}^{(j)}$
respectively in the unit cells of $L_{i}$ and $L_{j}$, and the disregistry
$\bm{b}^{i\rightarrow j}(\omega(\bm{r}))$ varies slowly on the atomic scale.

Fig. 4 demonstrates this construction. We first choose a reference point
$\bm{r}$; here, we use a lattice site of $L_{1}$ as an example. At the
reference point, the configuration is
$\omega=(\bm{b}^{(1)},\bm{b}^{(2)},\bm{b}^{(3)})$. At each one of the four
lattice sites in layer $L_{i}$ (respectively $L_{j}$) denoted by
$\bm{R}_{mn}^{(i)}$ (respectively $\bm{R}_{mn}^{(j)}$) neighboring the
reference point $\bm{r}$, the shift vector $\bm{b}^{(j)}(\bm{R}_{mn}^{(i)})$
coincides with the local disregistry of layer $L_{j}$ with respect to layer
$L_{i}$:

$\bm{b}^{i\rightarrow
j}_{mn}=\bm{R}_{mn}^{(j)}-\bm{R}_{mn}^{(i)}\equiv\bm{b}^{(j)}(\bm{R}_{mn}^{(i)}),$
(6)

where we use the definition in Eq. (3). Because the lattices are almost
aligned, these four disregistries only differ slightly from each other. We can
then define $\bm{b}^{i\rightarrow j}(\omega)$ by bilinear interpolation of the
four relative shifts $\bm{b}_{mn}^{2\rightarrow 3}$, yielding Eq. (5) Cazeaux
_et al._ (2019). Note that the interpolated disregistry given in Eq. (5) is
consistent with Eq. (6) at any lattice site of $L_{i}$.

In practice, we need to only consider reference points that coincide with
lattice points in layer $L_{i}$, $\bm{r}=\bm{R}^{(i)}$, in order to
parameterize atomic displacements. This leads us to introduce the reduced,
four-dimensional atomic configuration spaces:

$\Omega_{i}=\\{\omega\mid\omega=(\bm{b}^{(1)},\bm{b}^{(2)},\bm{b}^{(3)})\ |\
\bm{b}^{(i)}=\bm{0}\\}.$ (7)

The relaxation displacement $\bm{u}^{(i)}(\omega)$ at each lattice point is a
function of its local environment of shifts $\omega\in\Omega_{i}.$ The above
construction allows us to define an interpolated displacement
$\widehat{\bm{u}}^{(j)}(\omega)$ of another layer $L_{j}$ from the viewpoint
of a lattice site in $L_{i}$ by

$\displaystyle\widehat{\bm{u}}^{(j)}(\omega)=\bm{u}^{(j)}(\bm{b}^{j\rightarrow
1}(\omega),\bm{b}^{j\rightarrow 2}(\omega),\bm{b}^{j\rightarrow 3}(\omega)),$
(8)

where we simply use the definition given in Eq. (5). For example, the
interpolated displacement $\widehat{\bm{u}}^{(2)}(\omega)$ for
$\omega\in\Omega_{1}$ is

$\displaystyle\widehat{\bm{u}}^{(2)}(\omega)$
$\displaystyle=\bm{u}^{(2)}(\bm{b}^{2\rightarrow
1}(\omega),\bm{b}^{2\rightarrow 3}(\omega))$ (9)
$\displaystyle=\bm{u}^{(2)}(-A_{1}A_{2}^{-1}\bm{b}^{(2)},\bm{b}^{(3)}-A_{3}A_{2}^{-1}\bm{b}^{(2)}),$

where we omit $\bm{b}^{j\rightarrow j}(\omega)=0$ from the argument of
$\bm{u}^{(j)}.$ We can then describe the relaxation-modulated local shift
$\bm{B}^{i\rightarrow j}(\omega)$ of $L_{i}$ with respect to $L_{j}$ by

$\bm{B}^{i\rightarrow j}(\omega)=\bm{b}^{i\rightarrow
j}(\omega)+\widehat{\bm{u}}^{(j)}(\omega)-\bm{u}^{(i)}(\omega)\ \mathrm{for}\
\omega\in\Omega_{i},$ (10)

where $\widehat{\bm{u}}^{(j)}(\omega)$ is the interpolated displacement of
$L_{j}.$

The relaxed site energy has two contributions, interlayer and intralayer site
energies, as a function of the relaxation displacement vectors,
$\bm{u}=(\bm{u}^{(1)},\bm{u}^{(2)},\bm{u}^{(3)})$. For a given layer $L_{i}$,

$\Phi_{i,\omega}(\bm{u})=\Phi^{\mathrm{intra}}_{i,\omega}(\bm{u})+\Phi^{\mathrm{inter}}_{i,\omega}(\bm{u}).$
(11)

By the uniform sampling of the continuum of local configurations due to the
ergodicity of incommensurate lattices Cancès  _et al._ (2017); Massatt  _et
al._ (2016), we can derive the total energy by integrating the site energies
over all configurations Cazeaux  _et al._ (2019)

$\displaystyle
E^{\mathrm{tot}}(\bm{u})=E^{\mathrm{intra}}(\bm{u})+E^{\mathrm{inter}}(\bm{u})=\sum_{i=1}^{3}\int_{\Omega_{i}}\Phi_{i,\omega}(\bm{u})\,d\omega$
$\displaystyle\quad=\sum_{i=1}^{3}\int_{\Omega_{i}}\Phi^{\mathrm{intra}}_{i,\omega}(\bm{u})\,d\omega+\sum_{i=1}^{3}\int_{\Omega_{i}}\Phi^{\mathrm{inter}}_{i,\omega}(\bm{u})\,d\omega.$
(12)

We then minimize $E^{\mathrm{tot}}(\bm{u})$ to obtain the relaxation vector in
configuration space $\bm{u}$. We will provide the form of $E^{\mathrm{intra}}$
and $E^{\mathrm{inter}}$ in Sections. II.3 and II.4 respectively.

###  II.3 Intralayer Energy

We proceed to describe the intralayer energy contribution in configuration
space. From linear elasticity theory, the intralayer energy due to lattice
relaxation depends on the real space gradient of displacement vectors,

$\displaystyle E^{\mathrm{intra}}(\bm{u})$
$\displaystyle=\sum_{i=1}^{3}\int_{\Omega_{i}}\Phi^{\mathrm{intra}}_{i,\omega}(\bm{u})d\,\omega,$
(13)

where $\omega=(\bm{b}^{(j)},\bm{b}^{(k)})\in\Omega_{i}$, and

$\displaystyle\Phi^{\mathrm{intra}}_{i,\omega}(\bm{u})=\frac{1}{2}\gradient_{\bm{r}}\widehat{\bm{u}}^{(i)}(\omega):C:\gradient_{\bm{r}}\widehat{\bm{u}}^{(i)}(\omega).$
(14)

In Eq. (14), $C$ is the linear elasticity tensor. For in-plane deformation,
$C$ is a rank 4 tensor with its components defined as follows:

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

where $G$ and $K$ are the shear and bulk modulus of a monolayer (anisotropic
models can be used for anisotropic 2D materials such as black phosphorous).
These values can be obtained from first-principles total-energy calculations
based on density functional theory (DFT) by isotropically straining and
compressing the monolayer and performing a quadratic fit of the total energy
as a function of the applied strain or shear. In order to incorporate out-of-
plane relaxation, the elastic tensor would be a $3\times 3\times 3\times 3$
tensor Dai  _et al._ (2016), which is beyond the scope of this work. We can
compute the real space gradient of the interpolated displacement
$\widehat{\bm{u}}^{(i)}$ at an arbitrary configuration $\omega\in\Omega_{i}$
by the chain rule applied to Eqs. (3),(5), and (8):

$\gradient_{\bm{r}}\widehat{\bm{u}}^{(i)}(\omega)=\sum_{i\neq
j}\gradient_{\bm{b}^{(j)}}\bm{u}^{(i)}\cdot(A_{j}A_{i}^{-1}-I).$ (16)

Defining the following matrix,

$M_{ji}=A_{j}A_{i}^{-1}-I,$ (17)

we then have the Fourier representation of $\widehat{\bm{u}}^{(i)}:$

$\displaystyle\gradient_{\bm{r}}$
$\displaystyle\widehat{\bm{u}}^{(i)}(\bm{b}^{(j)},\bm{b}^{(k)})=i\sum_{\bm{G}^{(j)}}\sum_{\bm{G}^{(k)}}\operatorname{e}^{i\bm{G}^{(j)}\cdot\bm{b}^{(j)}}\operatorname{e}^{i\bm{G}^{(k)}\cdot\bm{b}^{(k)}}$
$\displaystyle\times\tilde{\bm{u}}^{(i)}(\bm{G}^{(j)},\bm{G}^{(k)})\otimes\left[M_{ji}^{T}\bm{G}^{(j)}+M_{ki}^{T}\bm{G}^{(k)}\right],$
(18)

where $\bm{G}^{(i)}$ are reciprocal lattice vectors for $L_{i}$ and the
Fourier coefficients $\tilde{\bm{u}}^{(i)}$ are defined according to

$\displaystyle\bm{u}^{(i)}(\bm{b}^{(j)},\bm{b}^{(k)})$
$\displaystyle=\sum_{\bm{G}^{(j)},\,\bm{G}^{(k)}}\tilde{\bm{u}}^{(i)}(\bm{G}^{(j)},\,\bm{G}^{(k)})$
$\displaystyle\quad\times\operatorname{e}^{i(\bm{G}^{(j)}\cdot\bm{b}^{(j)}+\bm{G}^{(k)}\cdot\bm{b}^{(k)})},$
(19)

for $j,k\neq i$ and $j<k$. The contribution of the mode
$(\bm{G}^{(j)},\bm{G}^{(k)})$ to the variational elastic energy is

$\displaystyle\mathcal{E}_{i}(\bm{G}^{(j)},\bm{G}^{(k)})$
$\displaystyle=\frac{1}{2}\left[\tilde{\bm{u}}^{(i)}(\bm{G}^{(j)},\bm{G}^{(k)})\otimes\left(M_{ji}^{T}\bm{G}^{(j)}+M_{jk}^{T}\bm{G}^{(k)}\right)\right]:C:$
$\displaystyle\qquad\left[\tilde{\bm{u}}^{(i)}(\bm{G}^{(j)},\bm{G}^{(k)})\otimes\left(M_{ji}^{T}\bm{G}^{(j)}+M_{jk}^{T}\bm{G}^{(k)}\right)\right].$
(20)

In terms of the $\mathcal{E}_{i}$’s, the total intralayer energy of the
trilayer system is

$\displaystyle
E^{\mathrm{intra}}(\bm{u})=\sum_{i=1}^{3}\sum_{\bm{G}^{(j)},\bm{G}^{(k)}}\mathcal{E}_{i}(\bm{G}^{(j)},\bm{G}^{(k)}).$
(21)

Each layer in the trilayer system has a four-dimensional configuration space
$\Omega_{i},$ but the in-plane interactions only span a two-dimensional
submanifold. The Euler-Lagrange PDE for the intralayer strain energy is thus
not properly elliptic, but it is nonetheless nonsingular since
$\mathcal{E}_{1}(\bm{G}^{(2)},\bm{G}^{(3)})=0$ if and only if
$\tilde{\bm{u}}_{1}(\bm{G}^{(2)},\bm{G}^{(3)})=0$ or
$\bm{G}^{(2)}=\bm{G}^{(3)}=0$, and similarly for
$\mathcal{E}_{2}(\bm{G}^{(1)},\bm{G}^{(3)})$ and
$\mathcal{E}_{3}(\bm{G}^{(1)},\bm{G}^{(2)})$. To see this, we take $i=1$ in
Eq. (21) and
$\mathcal{E}_{1}(\bm{G}^{(2)},\bm{G}^{(3)})\propto|M_{21}^{T}\bm{G}^{(2)}+M_{31}^{T}\bm{G}^{(3)}|^{2}$.
We note that

$M_{21}^{T}\bm{G}^{(2)}+M_{31}^{T}\bm{G}^{(3)}={\bm{G}}^{\prime(1)}-\bm{G}^{(2)}-\bm{G}^{(3)},$
(22)

where
${\bm{G}^{\prime}}^{(1)}=A_{1}^{-T}A_{2}^{T}\bm{G}^{(2)}+A_{1}^{-T}A_{3}^{T}\bm{G}^{(3)}$
is in the reciprocal lattice of $L_{1}$, and
${\bm{G}^{\prime}}^{(1)}-\bm{G}^{(2)}-\bm{G}^{(3)}=0$ if and only if
${\bm{G}^{\prime}}^{(1)}=\bm{G}^{(2)}=\bm{G}^{(3)}=0$ by the
incommensurability of the trilayer Cazeaux  _et al._ (2019). Observe that
$|{\bm{G}^{\prime}}^{(1)}-\bm{G}^{(2)}-\bm{G}^{(3)}|^{2}$ and hence
$\mathcal{E}_{1}(\bm{G}^{(2)},\bm{G}^{(3)})$ can be small even though
$\bm{G}^{(2)}$ and $\bm{G}^{(3)}$ are large, which is contrary to the
ellipticity condition that $\mathcal{E}_{1}(\bm{G}^{(2)},\bm{G}^{(3)})$ is a
uniformly positive definite quadratic form in $(\bm{G}^{(2)},\bm{G}^{(3)})$
for all Fourier coefficients of the relaxation vectors
$\tilde{\bm{u}}^{(1)}(\bm{G}^{(2)},\bm{G}^{(3)}).$ Physically, nonellipticity
implies the instability due to the resonance of long wavelength modes in an
infinite system. However, note that a finite physical system does not support
such an instability if the elastic energy Eq. (20) is simply positive. In
contrast, the PDEs for the intralayer energy of bilayer vdW heterostructures
are elliptic. For bilayers, the variational elastic energy of $L_{i}$ takes
the following form:

$\displaystyle\mathcal{E}_{i}(\bm{G}^{(j)})$
$\displaystyle=\frac{1}{2}\left(\tilde{\bm{u}}^{(i)}(\bm{G}^{(j)})\otimes
M_{ji}^{T}\bm{G}^{(j)}\right):C:$ (23)
$\displaystyle\qquad\left(\tilde{\bm{u}}^{(i)}(\bm{G}^{(j)})\otimes
M_{ji}^{T}\bm{G}^{(j)}\right).$

Similarly, taking $i=1$, we have
$\mathcal{E}_{1}(\bm{G}^{(2)})\propto|M_{21}^{T}\bm{G}^{(2)}|^{2}=|\bm{G}^{\prime(1)}-\bm{G}^{(2)}|^{2}$,
where $\bm{G}^{\prime(1)}=A_{1}^{-T}A_{2}^{T}\bm{G}^{(2)}$. In this case,
$|\bm{G}^{\prime(1)}-\bm{G}^{(2)}|$ is bounded from zero and takes its minimum
values when $\bm{G}^{\prime(1)}-\bm{G}^{(2)}$ are one of the columns of
$2\pi(A_{1}^{-T}-A_{2}^{-T})$, which can be understood as the primitive
vectors of the bilayer moiré supercell reciprocal lattice.

###  II.4 Interlayer Energy

Table 1: Equilibrium lattice constant of the primitive unit cell $a_{0}$ (in Å), unit cell area $\cal{A}$ (in Å2), elastic constants $K$ and $G$, and GSFE Fourier components for graphene, parallel and anti-parallel $\mathrm{WSe_{2}}$. Units of elastic constants and GSFE components are in meV per unit cell.  System |  $a_{0}\,(\mathrm{\AA}$) |  $\cal{A}\,(\mathrm{\AA^{2}}$) | $K$ | $G$ | $c_{0}$ | $c_{1}$ | $c_{2}$ | $c_{3}$ | $c_{4}$ | $c_{5}$
---|---|---|---|---|---|---|---|---|---|---
Graphene | 2.47 | 5.28 | 69518 | 47352 | 6.832 | 4.064 | -0.374 | -0.095 | 0.000 | 0.000
Parallel $\mathrm{WSe_{2}}$ | 3.28 | 9.32 | 47123 | 31153 | 27.85 | 14.08 | -2.442 | -0.773 | 0.000 | 0.000
Anti-parallel $\mathrm{WSe_{2}}$ | 3.28 | 9.32 | 47123 | 31153 | 35.03 | 12.17 | -1.899 | -0.1501 | 3.865 | 0.542

We now describe the interlayer energy due to the misfit from the twist of the
relaxation displacement. The origin of the interlayer energy is the slight
misalignment between the two layers. This interaction is represented by the
Generalized Stacking Fault Energy (GSFE). The GSFE is obtained with DFT by
calculating the difference in total energy between configurations that involve
relative rigid shifts between the two layers of a bilayer, with the zero
corresponding to the reference ground-state-energy configuration. The GSFE can
be expressed as a Fourier sum of the first few terms Kaxiras and Duesbery
(1993), and its functional form is given in Appendix A. The computational
details of the GSFE are provided in Appendix B and the Fourier coefficients
for graphene and $\mathrm{WSe_{2}}$ are provided in Table 1. The GSFE for
these two materials are shown in Fig. 5. We note that due to a lower symmetry,
for $\rm{WSe_{2}}$, there are two possible stackings of a bilayer, referred to
as parallel and anti-parallel, which denote the relative orientation of the
two layers in the reference configuration. In graphene and parallel
$\rm{WSe_{2}}$, the origin (AA stacking) is when the two layers are exactly
aligned, and AB/BA stacking is when one sublattice is aligned with the other.
In anti-parallel $\rm{WSe_{2}}$, we refer to the origin (BB stacking) as the
stacking configuration when the chalcogen atoms (Se) are aligned, the AA
stacking order as when the metal atoms (W) are aligned, and the AB stacking
order is when the metal atom is aligned with the chalcogen atom. The stacking
configuration is consistent with the notation in Phillips and Hellberg (2019).

Unrelaxed localized misfit site energies $\Phi^{\mathrm{misfit}}_{i\pm}$:
$\Gamma^{i\pm 1}\mapsto\mathbb{R}$ are defined by

$\Phi^{\mathrm{misfit}}_{i+}(\bm{b}^{i\rightarrow
i+1}(\omega))=\Phi^{\mathrm{misfit}}_{(i+1)-}(\bm{b}^{i+1\rightarrow
i}(\omega)).$ (24)

These site energies thus satisfy the following symmetry relation for $\bm{s}\
\in\ A_{i}^{-1}\Gamma_{i}=A_{i+1}^{-1}\Gamma_{i+1}\equiv[0,1)^{2}:$

$\Phi^{\mathrm{misfit}}_{(i+1)-}\left(A_{i}\bm{s}\right)=\Phi^{\mathrm{misfit}}_{i+}\left(-A_{i+1}\bm{s}\right),$
(25)

and we see that the misfit energies $\Phi^{\mathrm{misfit}}_{(i+1)-}$ and
$\Phi^{\mathrm{misfit}}_{i+}$ can be computed from the same misfit energy
between $L_{i}$ and $L_{i+1}$. As mentioned above, the simulations in this
paper approximate the misfit energy $\Phi^{\mathrm{misfit}}_{i\pm}$ by the
GSFE, $V^{\mathrm{GSFE}}_{i\pm}$, defined in Eq. (34) in Appendix A.

When the layers are relaxed, the argument of $V^{\mathrm{GSFE}}$ is modified
to be the relaxation-modulated local shift $\bm{B}^{i\rightarrow j}$ defined
in Eq. (10). The interlayer contribution of each layer can be given by

$\Phi^{\mathrm{inter}}_{i,\omega}(\bm{u})=\begin{cases}\frac{1}{2}V^{\mathrm{GSFE}}_{1+}(\bm{B}^{1\rightarrow
2}(\omega)),\hfill\mathrm{for}\ \omega\in\Omega_{1},\\\
\frac{1}{2}\left[V^{\mathrm{GSFE}}_{2+}(\bm{B}^{2\rightarrow
3}(\omega))+V^{\mathrm{GSFE}}_{2-}(\bm{B}^{2\rightarrow
1}(\omega))\right],\hfill\\\ \hfill\mathrm{for}\ \omega\in\Omega_{2},\\\
\frac{1}{2}V^{\mathrm{GSFE}}_{3-}(\bm{B}^{3\rightarrow
2}(\omega)),\hfill\mathrm{for}\ \omega\in\Omega_{3}.\end{cases}$ (26)

The total interlayer energy can therefore be obtained by integrating over all
local configurations $\omega$:

$\displaystyle
E^{\mathrm{inter}}=\frac{1}{2}\int_{\Omega_{1}}V^{\mathrm{GSFE}}_{1+}(\bm{B}^{1\rightarrow
2}(\omega))\,\mathrm{d}\omega$ (27)
$\displaystyle\quad+\frac{1}{2}\int_{\Omega_{2}}\left[V^{\mathrm{GSFE}}_{2-}(\bm{B}^{2\rightarrow
1}(\omega))+V^{\mathrm{GSFE}}_{2+}(\bm{B}^{2\rightarrow
3}(\omega))\right]\,\mathrm{d}\omega$
$\displaystyle\quad+\frac{1}{2}\int_{\Omega_{3}}V^{\mathrm{GSFE}}_{3-}(\bm{B}^{3\rightarrow
2}(\omega))\,\mathrm{d}\omega.$

Note that for hetero-trilayer materials, $V^{\mathrm{GSFE}}_{2-}$ and
$V^{\mathrm{GSFE}}_{2+}$ have different functional forms.

Figure 5: Left: the GSFE for bilayer graphene, parallel $\mathrm{WSe_{2}}$,
and anti-parallel $\mathrm{WSe_{2}}$. Right: side views of the AA stacking
orientation for each bilayer corresponding to the three cases on the left.

Figure 6: Two-dimensional line cuts of the relaxation pattern and relaxed GSFE
in configuration space for (a)-(d) trilayer graphene and (e)-(h) trilayer
anti-parallel $\mathrm{WSe_{2}}$ at
$\theta_{12}=0.5^{\circ},\theta_{23}=0.2^{\circ}$. (a) and (b) show the change
in the local shift vectors, $(\bm{B}^{i\rightarrow 2}-\bm{b}^{i\rightarrow
2})$. (a) shows $(\bm{B}^{1\rightarrow 2}-\bm{b}^{1\rightarrow 2})$, averaged
over the $L_{3}$ configuration, and (b) shows $(\bm{B}^{3\rightarrow
2}-\bm{b}^{3\rightarrow 2})$, averaged over the $L_{1}$ configuration. (c) and
(d) show the relaxed GSFE. (c) shows the sum of the GSFE between the
$L_{1},L_{2}$ bilayer pair, averaged over the $L_{3}$ configuration, and (d)
shows the sum of the GSFE between the $L_{2},L_{3}$ bilayer pair, averaged
over the $L_{1}$ configuration. (e)-(h) are the corresponding plots for anti-
parallel ${\mathrm{W}Se_{2}}$.

###  II.5 Relation between Real Space and Configuration Space

Finally, we describe the mapping between real space and configuration space.
Each adjacent pair of layers $i$ and $j$ forms a moiré supercell,
$\Gamma^{ij},$ given by the columns of the matrix
$(I-A_{i}A_{j}^{-1})^{-1}A_{i}=\left(A_{i}^{-1}-A_{j}^{-1}\right)^{-1}$. At
the real space position at a lattice site in $L_{i}$, $\bm{R}^{(i)}_{mn}$, we
can calculate the local configuration relative to $L_{j}$ by calculating the
local shift vector (see Eq. (6)):

$\displaystyle\bm{b}^{(j)}(\bm{R}_{mn}^{(i)})$
$\displaystyle=\bm{R}_{mn}^{(j)}-\bm{R}_{mn}^{(i)}=(A_{j}-A_{i})\begin{pmatrix}m\\\
n\end{pmatrix}$
$\displaystyle=(A_{j}A_{i}^{-1}-I)\bm{R}_{mn}^{(i)}\in\Gamma^{(j)}.$ (28)

Note that this construction agrees with evaluating the interpolated
disregistry $\bm{b}^{(j)}$ given by Eq. (3) at a lattice site
$\bm{R}^{(i)}_{mn}$ as the reference point. This result can be generalized to
any real space position $\bm{r}$:

$\displaystyle\bm{b}^{i\rightarrow
j}(\omega(\bm{r}))=(A_{j}A_{i}^{-1}-I)\bm{r}\in\Gamma^{(j)},$ (29)

where we notice that $\bm{b}^{i\rightarrow j}(\bm{r})$ is periodic in the
bilayer moiré supercell $\Gamma^{ij}$ in the continuum approximation, meaning
it sends a position in the bilayer supercell to the unit cell of $L_{j}$. This
generalization is equivalent to evaluating the interpolated disregistry given
in Eq. (5) at a configuration $\omega$.

The relation between the displacement defined in real space,
$\bm{U}^{(i)}(\bm{r})$, and in configuration space, $\bm{u}^{(i)}(\bm{r})$,
can then be given by evaluating $\widehat{\bm{u}}^{(i)}(\omega)$ at
$\bm{b}^{i\to j}(\bm{r})$ and $\bm{b}^{i\to k}(\bm{r})$ with Eq. (II.5) above
at a position $\bm{r}$ to obtain

$\bm{U}^{(i)}(\bm{r})=\bm{u}^{(i)}(\bm{b}^{i\rightarrow
j}(\omega(\bm{r})),\bm{b}^{i\rightarrow k}(\omega(\bm{r}))),$ (30)

where $j,k\neq i$ and $j<k.$ We see that $\bm{U}^{(i)}(\bm{r})$ is not
generally periodic since $\bm{b}^{(j)}(\bm{r})$ and $\bm{b}^{(k)}(\bm{r})$ are
each periodic on the incommensurate moiré supercells $\Gamma^{ij}$ and
$\Gamma^{ik}.$ In comparison, for incommensurate bilayers
$\bm{U}^{(i)}(\bm{r})=\widehat{\bm{u}}^{(i)}(\omega)$ is periodic on
$\Gamma^{ij}.$

##  III Results

###  III.1 General Twist Angles, $\theta_{12}\neq\theta_{23}$

Figure 7: Real space relaxation pattern for twisted trilayer graphene with
$\theta_{12}=1.73^{\circ}$, $\theta_{23}=2.24^{\circ}$, using discretization
$N=81$. Insets on the upper right corner zoom in on the black boxed areas near
the origin. (a) the relaxation displacement $\bm{U}^{(i)}(\bm{r})$. (b) Color
plot shows the curl of the displacement vectors
$-\gradient_{\bm{r}}\times\bm{U}^{(i)}(\bm{r})$ indicating the rotation
direction: blue is clockwise rotation and red is counterclockwise rotation,
and white means no rotation. Black parallelograms on $L_{1}$ and $L_{3}$ show
the approximate supercell vectors that correspond to $\theta_{12}$ and
$\theta_{23}$ respectively, calculated using Eq. (32). The small insets on the
upper right corner zoom in on the small purple boxed areas in the black
inserts. In the purple insets, blue and red dots show A and B sublattices
respectively, and the green parallelogram shows the monolayer cell.  Figure 8:
Same as Fig. 7, for trilayer $\mathrm{WSe_{2}}$, with parallel-stack between
$L_{1}$ and $L_{2}$ at $\theta_{12}=1.0^{\circ}$, and anti-parallel-stack
between $L_{2}$ and $L_{3}$ at $\theta_{23}=1.1^{\circ}$.

We perform energy minimization in configuration space with the discretized
grid $N\times N\times N\times N$ using the JULIA OPTIM package K Mogensen and
N Riseth (2018). We present our results in Fig. 6 for
$\theta_{12}=0.5^{\circ},\theta_{23}=0.2^{\circ}$ by taking two-dimensional
line cuts of the configuration space, using $N=36$. We show the results for
twisted trilayer graphene in Fig. 6a-(d) and for trilayer anti-parallel
$\mathrm{WSe_{2}}$ in Fig. 6(e)-(h). For the $L_{1},L_{2}$ bilayer pair, we
take the average over the configurations of $L_{3}$ (left), and similarly for
the $L_{2},L_{3}$ bilayer pair, we take the average over the configurations of
$L_{1}$ (right). For the relaxation displacement vectors, we show the change
in modulated local shift vector compared to the local shift vector without
relaxation between the bilayer pair, i.e., $(\bm{B}^{i\rightarrow
2}-\bm{b}^{i\rightarrow 2})$ for $i=1,3$, where $\bm{b}^{i\rightarrow 2}$ and
$\bm{B}^{i\rightarrow 2}$ are the local shift vector and the modulated local
shift vector respectively, given in Eqs. (5) and (10) respectively. For
example, for the relative relaxation between the $L_{1}$ and $L_{2}$ pair, we
plot the following quantity:

$\displaystyle\bm{B}^{1\rightarrow
2}(\bm{b}^{(2)},\bm{b}^{(3)})-\bm{b}^{1\rightarrow
2}(\bm{b}^{(2)},\bm{b}^{(3)})$ (31)
$\displaystyle=\widehat{\bm{u}}^{(2)}(A_{1}A_{2}^{-1}\bm{b}^{(2)},\bm{b}^{(3)}-A_{3}A_{2}^{-1}\bm{b}^{(2)})-\bm{u}^{(1)}(\bm{b}^{(2)},\bm{b}^{(3)}),$

as a function of $\bm{b}^{(2)}$, averaged over the $\bm{b}^{(3)}$. Physically,
this means displaying the relative relaxation vectors from $L_{2}$ to $L_{1}$,
averaged over all the $L_{3}$ shifts. Similarly, for $L_{2}$ and $L_{3}$, we
show $(\bm{B}^{3\rightarrow 2}-\bm{b}^{3\rightarrow 2})$ as a function of
$\bm{b}^{(2)}$, averaged over all the $L_{1}$ configurations. The relaxation
between $L_{2}$ and $L_{3}$ for both graphene and $\mathrm{WSe_{2}}$ is
stronger because of the smaller corresponding twist angle $\theta_{23}$. The
relaxation displacement vectors point towards the AB/BA stacking in graphene
and AB stacking in anti-parallel $\mathrm{WSe_{2}}$, which are the equilibrium
stacking orders. For both graphene and $\mathrm{WSe_{2}}$, the relaxed GSFE
exhibits domain wall formation between the adjacent bilayer pairs with
different thicknesses. A smaller corresponding twist angle gives rise to
stronger relaxation and thinner domain walls, similar to twisted bilayer
graphene. The shape of the domains are different depending on the lattice
symmetry: graphene has triangular domains while anti-parallel
$\mathrm{WSe_{2}}$ has hexagonal domains due to its broken inversion symmetry.

Figure 9: Comparison between unrelaxed and relaxed atomic positions in twisted
trilayer $\mathrm{WSe_{2}}$ with (a)-(b) parallel and (c)-(d) anti-parallel
stackings at $\theta_{12}=2.6^{\circ}$ and $\theta_{23}=2.3^{\circ}$. In all
panels, top left: atomic positions of $L_{1}$ and $L_{2}$ only (removing
$L_{3}$); top right: atomic positions of $L_{2}$ and $L_{3}$ only (removing
$L_{1}$); bottom: atomic positions of the three layers. Red and blue arrows
are guides to eyes, showing the bilayer moiré lengths that correspond to
$\theta_{12}=2.6^{\circ}$ and $\theta_{23}=2.3^{\circ}$.

Unlike in the case of twisted bilayer heterostructures, there does not exist
in trilayers a linear mapping at the continuous level between configuration
space and a periodic moiré supercell in real space. From the energy
minimization in configuration space, we obtain a set of displacement vectors,
$\bm{u}^{(i)}(\bm{b}^{(j)},\bm{b}^{(k)}),$ for each layer. The real space
displacement $\bm{U}^{(i)}(\bm{R}^{(i)})$ at lattice points $\bm{R}^{(i)}$ in
$L_{i}$ is then given by evaluating
$\bm{U}^{(i)}(\bm{R}^{(i)})={\bm{u}}^{(i)}(\bm{b}^{(j)}(\bm{R}^{(i)}),\bm{b}^{(k)}(\bm{R}^{(i)})).$
In contrast, the configuration space for the bilayer case is two-dimensional
and the displacement is more simply given by
$\bm{U}^{(i)}(\bm{R}^{(i)})={\bm{u}}^{(i)}(\bm{b}^{(j)}(\bm{R}^{(i)})).$ Fig.
7 shows the real space relaxation displacement vectors $\bm{U}^{(i)}(\bm{r})$
and the curl of $-\bm{U}^{(i)}(\bm{r})$,
$-(\gradient_{\bm{r}}\times\bm{U}^{(i)}(\bm{r}))$, which measures the rotation
component of the modulation displacement vector field, for twisted trilayer
graphene Zhang and Tadmor (2018). Physically, this quantity corresponds to the
change of the local twist angle, which can be measured using the SQUID-on-tip
technique Uri  _et al._ (2019). In Fig. 7a, we approximate the supercell
lattice constant that corresponds to the twist angle $\theta_{i,i+1}$
according to

$a_{\mathrm{sc}}(\theta_{i,i+1})=\frac{a_{0}}{2\sin(\theta_{i,i+1}/2)},$ (32)

where $a_{0}$ is the lattice constant of the graphene/$\mathrm{WSe_{2}}$ unit
cell given in Table 1. The length scale $a_{\mathrm{sc}}(\theta_{i,i+1})$ can
be understood as the periodicity of a twisted bilayer with the twist angle
$\theta_{i,i+1}$ in the continuum limit. Fig. 7a shows the relaxation
displacement vectors. Near the origin, the relaxation tends to follow the
local rotation to increase the local twist angle to reduce the size of the AA
spots, which agrees with the bilayer results Zhang and Tadmor (2018).

To better understand the relaxation pattern, Fig. 7b shows
$-(\gradient_{\bm{r}}\times\bm{U}^{(i)}(\bm{r}))$. On the bilayer moiré scales
(see insets on the upper right corners), the relaxation pattern has two length
scales: in $L_{1}$, the length scale of the relaxation is dominated by the
continuum supercell lattice constant of $\theta_{12}$,
$a_{\mathrm{sc}}(\theta_{12})$; in $L_{3}$, the relaxation has a length scale
$a_{\mathrm{sc}}(\theta_{23})$ instead. In the middle layer ($L_{2}$), the
relaxation pattern is a result of the interference between the two bilayer
moiré lengths, as a result of the coupling to both $L_{1}$ and $L_{3}$.
Zooming out, we observe structures at a larger length scale than both of the
bilayer moiré lengths, which is the “moiré of moiré,” due to the interference
between the two bilayer moiré cells. Note that these larger structures are not
identical to each other (most visible in Fig. 7b, Layer 2) due to the
incommensurability of the system. The emergence of the larger length scales
can also be explained by the nonellipticity of the system (see Section II.3).
The reciprocal space variational elastic energy components given in Eq. (20),
$\mathcal{E}_{i}(\bm{G}^{(j)},\bm{G}^{(k)}),$ have multipliers
$|\bm{G}^{\prime(i)}-\bm{G}^{(j)}-\bm{G}^{(k)}|^{2}$ which can take an
arbitrarily small value, which means the real space relaxation pattern can
have arbitrarily large structures. The quasi-periodic large structures we
observe in Fig. 7 manifests the nonellipticity of the Euler-Lagrangian PDE,
and since these structures are aperiodic, zooming out further will lead to the
appearance of even longer quasi-periodic structures. This shows our model’s
capability of capturing the incommensurability of the trilayer relaxation
pattern.

To show the generality of our model, we can examine the results for other
materials. For example, Fig. 8 shows the relaxation displacement vectors and
the corresponding curl in trilayer $\mathrm{WSe_{2}}$ heterostructures, with
parallel stacking between $L_{1}$ and $L_{2}$ and anti-parallel stacking
between $L_{2}$ and $L_{3}$. Similarly to the graphene case, we observe the
coupling between two length scales that corresponds to $\theta_{12}$ and
$\theta_{23}$. On the bilayer moiré length scales, the relaxation forms
triangular domains between $L_{1}$ and $L_{2}$ and hexagonal domains between
$L_{2}$ and $L_{3}$. On the moiré of moiré length scale, the overall domains
are also hexagonal, since the overall inversion symmetry of the system is
broken by anti-parallel $\mathrm{WSe_{2}}$. The shape can also be understood
as the interference pattern between the triangular domains and hexagonal
domains.

We note that when $\theta_{12}=-\theta_{23}$, $L_{1}$ and $L_{3}$ are aligned
with a twisted middle layer. The system can also be modeled as a bilayer
system, the interlayer energy of $L_{2}$ doubles and the strength of
relaxation in $L_{2}$ doubles. We used our trilayer model to calculate the
relaxation pattern for such systems and obtain a similar result as the
relaxation pattern shown in Carr  _et al._ (2019b).

We examine the real space atomic positions after relaxation. Fig. 9 compares
the unrelaxed and relaxed atomic positions for parallel and anti-parallel
$\mathrm{WSe_{2}}$ at $\theta_{12}=3^{\circ}$, $\theta_{23}=2.3^{\circ}$. Note
that the unrelaxed atomic positions are different in the parallel and anti-
parallel stackings (Fig. 9a, c). In both cases, the relaxation forms domain
between $L_{1}$ and $L_{2}$, as well as between $L_{2}$ and $L_{3}$, but at
different length scales, to enlarge equilibrium stacking orders. Fig. 9b shows
that in parallel $\mathrm{WSe_{2}}$, AB/BA stacking regions are enlarged and
the AA stacking regions shrink. Fig. 9d shows the domain formation in anti-
parallel $\mathrm{WSe_{2}}$, in which AB stacking regions are enlarged.

###  III.2 Equal Twist Angles, $\theta_{12}=\theta_{23}$

In the special case where $\theta_{12}=\theta_{23}$, we found that there exist
two types of solutions: symmetric and asymmetric. An example is shown in Fig.
10. The symmetric solution (Fig. 10a-b) preserves the two-fold rotational
symmetry along the $z$-axis (out-of-plane direction) of the unrelaxed
trilayers, whereas the asymmetric solution (Fig. 10c, d) exhibits a pinwheel-
shaped relaxation pattern, which remarkably breaks this symmetry. Near the
origin (local AAA stacking) of $L_{2}$, the relaxation has zero magnitude for
the symmetric solution, whereas the relaxation is finite for the asymmetric
solution. The asymmetric solution is found to have a lower energy in this
case: at $\theta_{12}=\theta_{23}=3^{\circ}$ for parallel $\rm{WSe_{2}}$, the
total energies per monolayer unit cell after relaxation are 39.34$\,\rm{meV}$
and 39.30$\,\rm{meV}$ for the symmetric and asymmetric solutions respectively.
Note that there exists two asymmetric solutions with equal energy, in which
the pinwheel rotates in the opposite directions.

Physically, the asymmetric solution redistributes the local AA stacking
density between the neighboring bilayer pairs to reduce the total AA areas and
consequently the total interlayer energy, at the cost of additional elastic
energy due to the pinwheel relaxation pattern. This can be seen from the
insets in Fig. 10. In the insets of Fig. 10b for $L_{1}$ and $L_{3}$, the
black parallelograms (bilayer moiré cells) agrees with the length scale of the
relaxation pattern. In $L_{2}$, relaxation is canceled out at the AAA spots
(white color near the origin) to preserve the two-fold rotational symmetry. In
Fig. 10d, however, the insets show that the relaxation near the origin in
$L_{1}$ and $L_{3}$ now has different length scales than the bilayer moiré
length corresponding to $3^{\circ}$. We also observe that the relaxation in
$L_{2}$ obtains a preferred net rotation, which breaks the symmetry between
the top and bottom layers. These results suggest that relaxation can
spontaneously break the symmetry of the system, unlike in twisted bilayers in
which the in-plane relaxation always preserves the symmetry.

Figure 10: Relaxation pattern for parallel $\mathrm{WSe_{2}}$ at
$\theta_{12}=\theta_{23}=3^{\circ}$, same as Fig. 7, showing the symmetric
solution in (a)-(b) and the asymmetric solution in (c)-(d).

A bifurcation point occurs at a critical twist angle, $\theta_{c}\approx
4.55^{\circ}$. Above $\theta_{c}$, the only solution is the symmetric one;
below $\theta_{c}$, the symmetric solution becomes an unstable saddle point
and the two counter-rotating asymmetric solutions become the global minima.
The existence of a critical twist angle can be understood from the competition
between reducing the misfit energy and increasing elastic energy. As the twist
angle increases, the moiré of moiré length scale becomes smaller, and
therefore elastic energy cost for the pinwheel pattern increases. Eventually,
the reduction in the misfit energy is not enough to compensate for the elastic
energy cost of the pinwheel pattern and the system no longer supports the
asymmetric solution.

Figure 11: (a) Relaxed total energy as a function of the twist angle shows the
bifurcation critical point is around $\theta=4.58^{\circ}$. (b)-(d) The curl
of relaxation displacement vectors in $L_{2}$ at the global minimum, (b)
below, (c) at, and (d) above the bifurcation point.

Fig. 11a shows the total energy after relaxation as a function of twist angle
in twisted trilayer parallel $\rm{WSe_{2}}$ for the two types of solution,
confirming the existence of a critical bifurcation point at $\theta_{c}\approx
4.58^{\circ}$. Above this $\theta_{c}$, the symmetric solution is the only
solution. Below $\theta_{c}$, the asymmetric solution is the global minimum
and the symmetric solution becomes an unstable saddle point and has slightly
higher energy. Fig. 11b-d shows the three types of local minima below, at, and
above the critical twist angle. In particular, Fig. 11c shows that at
$\theta_{c}$, the asymmetric solution has very weak pinwheel-like rotation and
it almost becomes a symmetric one. While we use twisted parallel
$\mathrm{WSe_{2}}$ as an example, the bifurcation also occurs in graphene or
any trilayers with inversion symmetry between AB and BA stacking orders.

##  IV Conclusion

We have presented a general formalism to obtain relaxation patterns in twisted
trilayer systems, which cannot be achieved with the traditional continuum
model in real space due to the aperiodicity of the system. We give results for
twisted trilayer graphene and $\mathrm{WSe_{2}}$. We introduce configuration
space as a natural description of the twisted trilayer system. We show that
the relaxation pattern of trilayer systems is the “moiré of moiré,” as a
result of the incommensurate coupling between two bilayer moiré lengths that
correspond to the two twist angles. The relaxation pattern manifests the
nonellipticity of the energy minimization problem in the sense that there are
infinitely many large length scales emerging in the system, and our model is
capable of capturing these large length scale modes. Similar to the bilayer
relaxation, in trilayers, the relaxation forms domains at the bilayer moiré
lengths to maximize the equilibrium stacking orders and minimize the high
energy stacking order (AA stack). In contrast to the in-plane relaxation
pattern in twisted bilayer graphene where the symmetry of the system is
preserved, we found that the twisted trilayers can have symmetry-breaking
relaxation patterns. When $\theta_{12}=\theta_{23}$, the trilayer relaxation
has two types of solutions: one symmetric solution that preserves the two-fold
rotational symmetry along the $z$-axis and the other asymmetric solution that
breaks this rotational symmetry. This symmetry breaking might be
experimentally observable effect, for example, through second harmonic
generation or scanning transmission electron microscopy.

Although we have only presented results for homo-trilayers in this work, our
approach can be also applied to hetero-trilayer materials. For materials with
similar lattice constants (such as graphene and hBN), we can simply apply a
small strain to the monolayer by assuming all three layers have the same
lattice constant. In this way, we can calculate the GSFE by performing rigid
shifts in a unit cell using DFT total energies as the interlayer coupling,
which effectively takes care of the small strain. For materials with very
different lattice constants, for example, graphene on top of transition metal
dichalcogenides, we need to modify the definition of $A_{i}$ matrices and be
careful about mapping from configuration space to real space. The
configuration space approach can also be generalized to four or more layers,
and each additional layer adds two dimensions to configuration space. However,
the extension would be very computationally expensive - the computational cost
scales as $N^{2(d-1)}$, where $d$ is the dimension and $N$ is the
discretization. Therefore, more sophisticated minimization methods should be
introduced for computational efficiency. Moreover, we would expect a symmetry-
breaking solution to be present for heterostructures with an odd number of
layers.

The mechanical relaxation effect is an important correction to be taken into
account for in modeling the electronic structure of twisted trilayers. The
relaxation pattern in twisted trilayer may be used as a new platform to
engineer strain in vdW heterostructures. The second twist angle is an
additional degree of freedom to manipulate the electronic states and explore
unconventional strongly correlated states Tsai  _et al._ (2019).

###### Acknowledgements.

We thank Stephen Carr, Shiang Fang, Steven Torrisi, and Emine Kucukbenli,
Philip Kim, Liang Fu, Ke Wang for insightful discussions. This work was
supported by the STC Center for Integrated Quantum Materials, NSF Grant No.
DMR-1231319, ARO MURI Award W911NF-14-0247, NSF DMREF Award 1922165, and NSF
Award DMS-1819220. Calculations were performed on the Odyssey cluster
supported by the FAS Division of Science, Research Computing Group at Harvard
University.

##  Appendix A Generalized Stacking Fault Energy

The GSFE is the total energy at different stacking configurations without
rotation between a bilayer. Letting $(b_{x},b_{y})$ be the relative
displacement or the local stacking order between two adjacent layers, we
define the following vector $\bm{v}=(v,w)\in[0,2\pi]^{2}$:

$\begin{pmatrix}v\\\
w\end{pmatrix}=\frac{2\pi}{\alpha}\matrixquantity[\sqrt{3}/2&-1/2\\\
\sqrt{3}/2&1/2]\begin{pmatrix}b_{x}\\\ b_{y}\end{pmatrix}.$ (33)

In terms of $\bm{v}$, the GSFE of a hexagonal lattice can be expressed as a
Fourier sum as Kaxiras and Duesbery (1993); Zhou  _et al._ (2015); Carr  _et
al._ (2018),

$\displaystyle V^{\mathrm{GSFE}}_{j\pm}=c_{0}+$ $\displaystyle c_{1}(\cos
v+\cos w+\cos(v+w))$ $\displaystyle+$ $\displaystyle
c_{2}(\cos(v+2w)+\cos(v-w)+\cos(2v+w))$ $\displaystyle+$ $\displaystyle
c_{3}(\cos(2v)+\cos(2w)+\cos(2v+2w))$ $\displaystyle+$ $\displaystyle
c_{4}(\sin v+\sin w-\sin(v+w))$ (34) $\displaystyle+$ $\displaystyle
c_{5}(\sin(2v+2w)-\sin(2v)-\sin(2w)),$

where we use the subscript $j\pm$ to denote the GSFE between the $L_{j}$ and
$L_{j\pm 1}$. We note that we are only including the first few terms in the
Fourier series because of the lattice symmetry and higher order terms are
three or more orders of magnitude smaller than the leading order terms. For
other lattice geometries, it may be necessary to include higher order terms in
the Fourier series.

##  Appendix B Density Functional Theory Calculations

The coefficients for the interlayer and intralayer energies are obtained from
DFT calculations with the Vienna Ab initio Simulation Package (VASP) Kresse
and Hafner (1993); Kresse and Furthmüller (1996). The van der Waals force is
impletemented through the vdW-DFT method using the SCAN+rVV10 functional Dion
_et al._ (2004); Vydrov and Van Voorhis (2010); Klimeš  _et al._ (2011); Peng
_et al._ (2016). The values are presented in Table 1. The elasticity constants
and interlayer coupling coefficients of graphene and $\mathrm{MoS_{2}}$ are
also given in S.Carr et alCarr  _et al._ (2018).

To obtain the GSFE, we performed a rigid shift between a bilayer using a
$9\times 9$ grid in the unit cell. We extract the ground state energy from
DFT, by allowing relaxation in the out-of-plane direction but fixing the in-
plane direction. Fig. 5 shows the GSFE of graphene and $\mathrm{WSe_{2}}$,
with high-symmetry stacking orders, AA stacking and AB/BA stacking, marked. In
graphene, the AB/BA stacking (when the A sublattice of one layer is aligned
with the B sublattice of the second layer) is the equilibrium stacking
configuration and the AA stacking (the origin) has the highest energy. In
TMDCs, $\mathrm{MX_{2}}$, with M being a transition metal atom (e.g., Mo, W)
and X being a chalcogen atom (S, Se, Te), the bilayer can have two
orientations when performing the rigid shifts because the two sublattices have
different atomic species. In parallel $\mathrm{WSe_{2}}$, the two bilayers are
aligned. The GSFE energy landscape is similar to that of graphene. In anti-
parallel $\mathrm{WSe_{2}}$, however, the inversion symmetry is broken. In
this case, the equilibrium stacking order is when the metal atoms and the
chalcogen atoms are aligned (AB stacking).

The GSFE is a more accurate description of interlayer coupling compared to
empirical interlayer coupling potentials when they exist, such as for graphene
and hBN, and it enables the description of interlayer coupling in TMDCs for
which empirical interlayer coupling potentials do not exist.

## References

  * Cao  _et al._ (2018a) Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature (London) 556, 43–50 (2018a), arXiv:1803.02342 [cond-mat.mes-hall] .
  * Cao  _et al._ (2018b) Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L. Tomarken, Jason Y. Luo, Javier D. Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, Ray C. Ashoori, and Pablo Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature (London) 556, 80–84 (2018b), arXiv:1802.00553 [cond-mat.mes-hall] .
  * Yankowitz  _et al._ (2018) Matthew Yankowitz, Jeil Jung, Evan Laksono, Nicolas Leconte, Bheema L. Chittari, K. Watanabe, T. Taniguchi, Shaffique Adam, David Graf, and Cory R. Dean, “Dynamic band-structure tuning of graphene moiré superlattices with pressure,” Nature (London) 557, 404–408 (2018), arXiv:1707.09054 [cond-mat.mes-hall] .
  * Liu  _et al._ (2019) Xiaomeng Liu, Zeyu Hao, Eslam Khalaf, Jong Yeon Lee, Kenji Watanabe, Takashi Taniguchi, Ashvin Vishwanath, and Philip 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] .
  * Cao  _et al._ (2019) Yuan Cao, Daniel Rodan-Legrain, Oriol Rubies-Bigorda, Jeong Min Park, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, “Electric Field Tunable Correlated States and Magnetic Phase Transitions in Twisted Bilayer-Bilayer Graphene,” arXiv e-prints , arXiv:1903.08596 (2019), arXiv:1903.08596 [cond-mat.str-el] .
  * Burg  _et al._ (2019) G. William Burg, Jihang Zhu, Takashi Taniguchi, Kenji Watanabe, Allan H. MacDonald, and Emanuel Tutuc, “Correlated Insulating States in Twisted Double Bilayer Graphene,” arXiv e-prints , arXiv:1907.10106 (2019), arXiv:1907.10106 [cond-mat.mes-hall] .
  * Shen  _et al._ (2019) Cheng Shen, Na Li, Shuopei Wang, Yanchong Zhao, Jian Tang, Jieying Liu, Jinpeng Tian, Yanbang Chu, Kenji Watanabe, Takashi Taniguchi, Rong Yang, Zi Yang Meng, Dongxia Shi, and Guangyu Zhang, “Observation of superconductivity with Tc onset at 12K in electrically tunable twisted double bilayer graphene,” arXiv e-prints , arXiv:1903.06952 (2019), arXiv:1903.06952 [cond-mat.supr-con] .
  * Lu  _et al._ (2019) Xiaobo Lu, Petr Stepanov, Wei Yang, Ming Xie, Mohammed Ali Aamir, Ipsita Das, Carles Urgell, Kenji Watanabe, Takashi Taniguchi, Guangyu Zhang, Adrian Bachtold, Allan H. MacDonald, and Dmitri K. Efetov, “Superconductors, Orbital Magnets, and Correlated States in Magic Angle Bilayer Graphene,” arXiv e-prints , arXiv:1903.06513 (2019), arXiv:1903.06513 [cond-mat.str-el] .
  * Chen  _et al._ (2019) Guorui Chen, Lili Jiang, Shuang Wu, Bosai Lyu, Hongyuan Li, Bheema Lingam Chittari, Kenji Watanabe, Takashi Taniguchi, Zhiwen Shi, Jeil Jung, Yuanbo Zhang, and Feng Wang, “Evidence of a gate-tunable Mott insulator in a trilayer graphene moiré superlattice,” Nature Physics 15, 237–241 (2019), arXiv:1803.01985 [cond-mat.mes-hall] .
  * Moriyama  _et al._ (2019) Satoshi Moriyama, Yoshifumi Morita, Katsuyoshi Komatsu, Kosuke Endo, Takuya Iwasaki, Shu Nakaharai, Yutaka Noguchi, Yutaka Wakayama, Eiichiro Watanabe, Daiju Tsuya, Kenji Watanabe, and Takashi 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] .
  * Tritsaris  _et al._ (2016) Georgios A. Tritsaris, Sharmila N. Shirodkar, Efthimios Kaxiras, Paul Cazeaux, Mitchell Luskin, Petr Plecháč, and Eric Cancès, “Perturbation theory for weakly coupled two-dimensional layers,” Journal of Materials Research 31, 959–966 (2016).
  * San-Jose  _et al._ (2014) Pablo San-Jose, A. Gutiérrez-Rubio, Mauricio Sturla, and Francisco Guinea, “Spontaneous strains and gap in graphene on boron nitride,” Phys. Rev. B 90, 075428 (2014), arXiv:1404.7777 [cond-mat.mes-hall] .
  * Dai  _et al._ (2016) Shuyang Dai, Yang Xiang, and David J. Srolovitz, “Twisted Bilayer Graphene: Moiré with a Twist,” Nano Letters 16, 5923–5927 (2016).
  * Nam and Koshino (2017) Nguyen N. T. Nam and Mikito 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) Malena I. Español, Dmitry Golovaty, and J. Patrick 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) Malena I. Español, Dmitry Golovaty, and J. Patrick 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] .
  * Carr  _et al._ (2018) Stephen Carr, Daniel Massatt, Steven B. Torrisi, Paul Cazeaux, Mitchell Luskin, and Efthimios Kaxiras, “Relaxation and domain formation in incommensurate two-dimensional heterostructures,” Phys. Rev. B 98, 224102 (2018), arXiv:1805.06972 [cond-mat.mes-hall] .
  * Leconte  _et al._ (2019) Nicolas Leconte, Srivani Javvaji, Jiaqi An, and Jeil 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] .
  * Coppersmith  _et al._ (1982) SN Coppersmith, Daniel S Fisher, BI Halperin, PA Lee, and William F Brinkman, “Dislocations and the commensurate-incommensurate transition in two dimensions,” Physical Review B 25, 349 (1982).
  * Alden  _et al._ (2013) Jonathan S. Alden, Adam W. Tsen, Pinshane Y. Huang, Robert Hovden, Lola Brown, Jiwoong Park, David A. Muller, and Paul L. McEuen, “Strain solitons and topological defects in bilayer graphene,” Proceedings of the National Academy of Science 110, 11256–11260 (2013), arXiv:1304.7549 [cond-mat.mes-hall] .
  * Woods  _et al._ (2014) C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C. Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gorbachev, A. V. Kretinin, J. Park, L. A. Ponomarenko, M. I. Katsnelson, Yu. N. Gornostyrev, K. Watanabe, T. Taniguchi, C. Casiraghi, H. J. Gao, A. K. Geim, and K. S. Novoselov, “Commensurate-incommensurate transition in graphene on hexagonal boron nitride,” Nature Physics 10, 451–456 (2014), arXiv:1401.2637 [cond-mat.mes-hall] .
  * van Wijk  _et al._ (2015) M. M. van Wijk, A. Schuring, M. I. Katsnelson, and A. Fasolino, “Relaxation of moiré patterns for slightly misaligned identical lattices: graphene on graphite,” 2D Materials 2, 034010 (2015), arXiv:1503.02540 [cond-mat.mes-hall] .
  * Yankowitz  _et al._ (2016) Matthew Yankowitz, K. Watanabe, T. Taniguchi, Pablo San-Jose, and Brian J. Leroy, “Pressure-induced commensurate stacking of graphene on boron nitride,” Nature Communications 7, 13168 (2016), arXiv:1603.03244 [cond-mat.mes-hall] .
  * Kim  _et al._ (2017) Na Yeon Kim, Hu Young Jeong, Jung Hwa Kim, Gwangwoo Kim, Hyeon Suk Shin, and Zonghoon Lee, “Evidence of local commensurate state with lattice match of graphene on hexagonal boron nitride,” ACS Nano 11, 7084–7090 (2017).
  * Yoo  _et al._ (2019) Hyobin Yoo, Rebecca Engelke, Stephen Carr, Shiang Fang, Kuan Zhang, Paul Cazeaux, Suk Hyun Sung, Robert Hovden, Adam W. Tsen, Takashi Taniguchi, Kenji Watanabe, Gyu-Chul Yi, Miyoung Kim, Mitchell Luskin, Ellad B. Tadmor, Efthimios Kaxiras, and Philip Kim, “Atomic and electronic reconstruction at the van der Waals interface in twisted bilayer graphene,” Nature Materials 18, 448–453 (2019), arXiv:1804.03806 [cond-mat.mtrl-sci] .
  * Tsai  _et al._ (2019) Kan-Ting Tsai, Xi Zhang, Ziyan Zhu, Yujie Luo, Stephen Carr, Mitchell Luskin, Efthimios Kaxiras, and Ke Wang, “Correlated Superconducting and Insulating States in Twisted Trilayer Graphene Moire of Moire Superlattices,” arXiv e-prints , arXiv:1912.03375 (2019), arXiv:1912.03375 [cond-mat.mes-hall] .
  * Mora  _et al._ (2019) Christophe Mora, Nicolas Regnault, and B. Andrei Bernevig, “Flatbands and Perfect Metal in Trilayer Moiré Graphene,” Phys. Rev. Lett. 123, 026402 (2019), arXiv:1901.05469 [cond-mat.mes-hall] .
  * Ma  _et al._ (2019) Zhen Ma, Shuai Li, Ya-Wen Zheng, Meng-Meng Xiao, Hua Jiang, Jin-Hua Gao, and X. C. Xie, “Topological flat bands in twisted trilayer graphene,” arXiv e-prints , arXiv:1905.00622 (2019), arXiv:1905.00622 [cond-mat.mes-hall] .
  * Carr  _et al._ (2019a) Stephen Carr, Shiang Fang, Ziyan Zhu, and Efthimios Kaxiras, “An exact continuum model for low-energy electronic states of twisted bilayer graphene,” arXiv e-prints , arXiv:1901.03420 (2019a), arXiv:1901.03420 [cond-mat.mes-hall] .
  * Guinea and Walet (2019) Francisco Guinea and Niels R. Walet, “Continuum models for twisted bilayer graphene: Effect of lattice deformation and hopping parameters,” Phys. Rev. B 99, 205134 (2019), arXiv:1903.00364 [cond-mat.str-el] .
  * Massatt  _et al._ (2016) Daniel Massatt, Mitchell Luskin, and Christoph Ortner, “Electronic Density of States for Incommensurate Layers,” arXiv e-prints , arXiv:1608.01968 (2016), arXiv:1608.01968 [math-ph] .
  * Cazeaux  _et al._ (2016) Paul Cazeaux, Mitchell Luskin, and Ellad 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) Eric Cancès, Paul Cazeaux, and Mitchell 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) Paul Cazeaux, Mitchell Luskin, and Daniel Massatt, “Energy Minimization of Two Dimensional Incommensurate Heterostructures,” Archive for Rational Mechanics and Analysis , 122 (2019), arXiv:1806.10395 [physics.comp-ph] .
  * Kaxiras and Duesbery (1993) Efthimios Kaxiras and Michael S. Duesbery, “Free energies of generalized stacking faults in Si and implications for the brittle-ductile transition,” Phys. Rev. Lett. 70, 3752–3755 (1993).
  * Phillips and Hellberg (2019) Madeleine Phillips and C. Stephen Hellberg, “Commensurate structures in twisted transition metal dichalcogenide heterobilayers,” arXiv e-prints , arXiv:1909.02495 (2019), arXiv:1909.02495 [cond-mat.mtrl-sci] .
  * K Mogensen and N Riseth (2018) Patrick K Mogensen and Asbjørn N Riseth, “Optim: A mathematical optimization package for Julia,” The Journal of Open Source Software 3, 615 (2018).
  * Zhang and Tadmor (2018) Kuan Zhang and Ellad B. Tadmor, “Structural and electron diffraction scaling of twisted graphene bilayers,” Journal of Mechanics Physics of Solids 112, 225–238 (2018).
  * Uri  _et al._ (2019) Aviram Uri, Sameer Grover, Yuan Cao, J. A. Crosse, Kousik Bagani, Daniel Rodan-Legrain, Yuri Myasoedov, Kenji Watanabe, Takashi Taniguchi, Pilkyung Moon, Mikito Koshino, Pablo Jarillo-Herrero, and Eli Zeldov, “Mapping the twist angle and unconventional Landau levels in magic angle graphene,” arXiv e-prints , arXiv:1908.04595 (2019), arXiv:1908.04595 [cond-mat.mes-hall] .
  * Carr  _et al._ (2019b) Stephen Carr, Chenyuan Li, Ziyan Zhu, Efthimios Kaxiras, Subir Sachdev, and Alex Kruchkov, “Coexistence of ultraheavy and ultrarelativistic Dirac quasiparticles in sandwiched trilayer graphene,” arXiv e-prints , arXiv:1907.00952 (2019b), arXiv:1907.00952 [cond-mat.str-el] .
  * Zhou  _et al._ (2015) Songsong Zhou, Jian Han, Shuyang Dai, Jianwei Sun, and David J. Srolovitz, “van der Waals bilayer energetics: Generalized stacking-fault energy of graphene, boron nitride, and graphene/boron nitride bilayers,” Phys. Rev. B 92, 155438 (2015).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, “Ab initio molecular dynamics for liquid metals,” Phys. Rev. B 47, 558–561 (1993).
  * 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–11186 (1996).
  * Dion  _et al._ (2004) M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, “Van der Waals Density Functional for General Geometries,” Phys. Rev. Lett. 92, 246401 (2004), arXiv:cond-mat/0402105 [cond-mat.mtrl-sci] .
  * Vydrov and Van Voorhis (2010) Oleg A. Vydrov and Troy Van Voorhis, “Nonlocal van der Waals density functional: The simpler the better,” J. Chem. Phys. 133, 244103–244103 (2010), arXiv:1009.1421 [physics.chem-ph] .
  * Klimeš  _et al._ (2011) Jiří Klimeš, David R. Bowler, and Angelos Michaelides, “Van der Waals density functionals applied to solids,” Phys. Rev. B 83, 195131 (2011), arXiv:1102.1358 [cond-mat.mtrl-sci] .
  * Peng  _et al._ (2016) Haowei Peng, Zeng-Hui Yang, John P. Perdew, and Jianwei Sun, “Versatile van der Waals Density Functional Based on a Meta-Generalized Gradient Approximation,” Physical Review X 6, 041005 (2016).