# Accurate force-field methodology capturing atomic reconstructions in
transition metal dichalcogenide moiré systems

Carl Emil Mørch Nielsen  Miguel da Cruz  Abderrezak Torche  Gabriel Bester
Institute of Physical Chemistry, University of Hamburg, 22607 Hamburg, Germany

###### Abstract

In this work, a generalized force-field methodology for the relaxation of
large moiré heterostructures is proposed. The force-field parameters are
optimized to accurately reproduce the structural degrees of freedom of some
computationally manageable cells relaxed using density functional theory. The
parameters can then be used to handle large moiré systems. We specialize to
the case of 2H-phased twisted transition-metal dichalcogenide homo- and
heterobilayers using a combination of the Stillinger-Weber intralayer- and the
Kolmogorov-Crespi interlayer-potential. Force-field parameters are developed
for all combinations of MX2 for $\text{M}\in\\{\text{Mo},\text{W}\\}$ and
$\text{X}\in\\{\text{S},\text{Se},\text{Te}\\}$. The results show agreement
within 20 meV in terms of band structure between density functional theory and
force-field relaxation. Using the relaxed structures, a simplified and
systematic scheme for the extraction of the interlayer moiré potential is
presented for both R- and H-stacked systems. We show that in-plane and out-of-
plane relaxation effects on the moiré potential, which is made both deeper and
wider after relaxation, are essential. An interpolation based methodology for
the calculation of the interlayer binding energy is also proposed. Finally, we
show that atomic reconstruction, which is captured by the force-field method,
becomes especially prominent for angles below 4-5∘, when there is no mismatch
in lattice constant between layers.

††preprint: APS/123-QED

##  I Introduction

Two dimensional (2D) moiré systems are currently an especially attractive
playground for new technological applications [1, 2, 3, 4]. Lattice mismatch
combined with the twist angle between the constituent layers allows for an
ingenious way of external mechanical control of the moiré period and thus the
resulting electronic properties. Without a doubt, the pioneering discovery of
twisted bilayer graphene and its magic angle of 1.05∘ [5] was the major
driving force toward the study of 2D heterostructures and constituted the
basis for the field of twistronics. An interesting and widely studied class of
moiré systems is the 2D family of transition metal dichalcogenides (TMDs),
featuring strong light-matter interaction and large spin-orbit coupling with a
sizable bandgap [6]. A fundamental advantage of TMDs is that flat minibands
are not only realised at specific angles, but exist in a continuum of small
angles [7]. An example of a moiré-structured TMD-system can be seen in Fig. 1.
Moreover, experimental and theoretical findings of the excited states in type-
II aligned heterostructured TMDs show evidence of spatially indirect excitons
localized within certain registries of the moiré structure [8, 9, 7, 10].
Moiré structured TMDs provide a platform for studying correlated quantum
phenomena [11] including hole Mott insulator states at integer and fractional
fillings with generalized Wigner crystallization, essentially creating a
Fermi-Hubbard system [12, 13, 14, 15, 16]. Moiré structured TMD bilayer
systems also allows for realization of Bose-Hubbard physics with excitons
trapped in a periodic triangular potential and subject to strong Coulomb
interactions [17].

Figure 1:  WS2 on MoS2 twisted at an angle of $6.0^{\circ}$. Sulphur atoms are
shown in orange, molybdenum in blue, and tungsten in grey. The moiré unit cell
is shown with black solid lines, and has a moiré period, $m_{0}$, of 30.1 Å
encompassing 546 atoms. The long diagonal, $r_{diag}$, is marked with a dashed
black line.

Moiré physics in TMDs are largely determined by the shape of the twist-induced
moiré potential, which arises from local stacking configurations, lattice
corrugation and, for small angles, atomic reconstruction [18, 19, 20, 7, 21,
22, 23, 24, 25]. As a consequence, relaxation effects are important for
numerical simulations that involve moiré structured TMD systems prone to
atomic reconstruction, and/or structures with a moiré period large enough to
corrugate the individual layers [26, 27]. From an ab initio standpoint, this
presents a large challenge owing to the fact that relaxation is a
computational bottleneck in such calculations. In an excellent paper by Naik
et al. [28], a method to overcome this problem is suggested by using a force-
field model based on a combination of the Stillinger-Weber (SW) [29, 30] and
Kolmogorov-Crespi (KC) [31, 32] potentials. The SW force-field accurately
describes the intralayer forces, while the KC potential captures van der Waals
(vdW) interaction between layers and includes a stacking-dependent term.
Previously, this had been parametrized and applied to graphene and hexagonal
boron nitride [33, 34, 35, 36], but is now also available for MX2
homobilayers, where $\text{M}\in\\{\text{Mo},\text{W}\\}$ and
$\text{X}\in\\{\text{S},\text{Se}\\}$ [28, 37].

However, the parameters presented in Ref. [28] are somewhat inaccurate when
comparing to density functional theory (DFT) calculated results, e.g. for some
structures, the bandgap is inaccurate by up to 100 meV. Even more importantly,
the band curvature and energetic position of e.g. the lowest conduction band
and highest valence band are skewed on similar scales. In Ref. [28], the
parameters are developed by fitting to DFT binding energies which will not
guaranty the force-field model to reproduce the DFT relaxed structure. In this
work, the structural parameters of the DFT optimized structures (i.e. atomic
positions and unit cell size) are used directly as target values for the
optimization of the force-field parameters. Furthermore, the KC-
parametrization of Ref. [28] is presented on a per interaction basis, meaning
that atom-atom interactions are considered the same for different systems,
e.g. S-S parameters for MoS2\- and WS2-bilayers are the same. However, from a
fundamental point of view, vdW interaction, being of long-range nature, is
known to be sensitive to the surrounding environment. As such, we
reparametrize the KC-potential on a per system basis, which yields more
accurate band structures. Furthermore, we expand the set of parameters to
include heterobilayers with and without lattice mismatch, essentially covering
all bilayer combinations of 2H-phased MX2 for
$\text{M}\in\\{\text{Mo},\text{W}\\}$ and
$\text{X}\in\\{\text{S},\text{Se},\text{Te}\\}$. However, the method presented
here is, in principle, extendable to any 2D moiré structure and not limited to
TMDs. Our force-field parameters, along with a variety of relaxed structures
can be found via Ref. [38].

Lastly, we present two interpolation-based schemes to describe the interlayer
exciton moiré potential of lattice-matched heterostrutures with type-II band
alignment by using a combination of the force-field method and DFT, which
provides easy access to the potential for almost any angle. We extend this
analogy to the binding energy, which allows for visualization of atomic
reconstruction and the rate at which the reconstructed domains form with
decreasing twist angle. Specifically, we see that atomic reconstruction
becomes significant for angles below 4-5∘ for the TMD heterostructures studied
here.

##  II Methodology

The first step is to develop the SW-parameters, which is done by considering
the constituting monolayers one at a time. For 2H-phased TMD monolayers, the
hexagonal symmetry reduces the structural degrees of freedom into two (target)
parameters only, namely the lattice constant, $a_{0}$, and the intralayer
distance, $d_{intra}$, i.e. the out-of-plane X-X distance. Therefore, the SW-
parametrization is carried out using $a_{0}$ and $d_{intra}$ as targets and
reproduces them extremely well. The force-field relaxations are performed
using the LAMMPS package [39], and the optimization of parameters is carried
out with use of the Dakota package [40].

For the optimization of the KC parameters we are following two strategies,
depending on whether the constituting layers are lattice matched or not.

###  II.1 Lattice matched bilayers

Bilayers that have the same chalcogen atom have a lattice constant mismatch
$\delta\sim 0.1\%$ and are treated as lattice matched. In this case, only one
additional structural parameter is considered, namely the interlayer spacing,
$d_{inter}$ (M-M distance). The KC-parameters are obtained by fitting to
$a_{0}$, $d_{intra}$ and $d_{inter}$ for the six high-symmetry stacking
configurations (HSSCs), while keeping the SW-parameters fixed. The HSSCs are
depicted in Fig. 2 and are divided in two groups, namely R- and H-stacking,
which differ by a rotation of one of the layers by 60∘.

Figure 2: The six high-symmetry stacking configurations of a bilayer with no
lattice-mismatch. (a)-(c) ((d)-(f)) belong to the R(H)-stacking group. The
stacking $\mathrm{R_{X}^{X}}/\mathrm{R_{M}^{M}}$ is also referred to as AA-
stacking, and $\mathrm{H_{X}^{M}}/\mathrm{H_{M}^{X}}$ as AB-stacking. The
dotted black lines indicate atoms that coincide along $z$, justifying the
naming convention.

This procedure follows the idea that the mechanical properties of the single
layer is well described by the SW potential and is not altered by the
interlayer interaction (KC potential). It is crucial to derive a force-field
that is transferable between the different stackings since the twisted
bilayers correspond to combinations of three different stackings, as will be
demonstrated subsequently.

Moreover, as we will indirectly show in Sec. IV.1, every subcell of a lattice-
matched moiré unit-cell is, to a certain extent, well described by a
superposition of the HSSCs. Note, that this is not the case for lattice-
mismatched systems where no local HSSCs can be identified. Justification of
our methodology becomes trivial for smaller angles, where domains of the HSSCs
make up a large fraction of the moiré unit cell. Finally, the small unit cells
constructed with merely six atoms, makes both DFT calculations and the
optimization schemes of Dakota and LAMMPS relatively fast.

###  II.2 Lattice mismatched bilayers

In Table 1 we show the lattice mismatch $\delta$ for the different
combinations of chalcogen atoms (the metal atom is nearly irrelevant for the
lattice constant). The lattice-mismatch of the systems investigated here
(X=S,Se,Te) is so large that the construction of small six atom unit cells as
done in the lattice matched case is not meaningful. The in-plane strain will
radically change the electronic properties [41].

Table 1:  Angles chosen for fitting lattice-mismatched structures accompanied by the lattice-mismatch ($\delta$, found using DFT), number of atoms ($n_{atom}$), and the moiré lattice constant $m_{0}$ (moiré period). X1 | X2 |  $\theta$ (∘) |  $\delta$ (%) | $n_{atom}$ |  $m_{0}$ (nm)
---|---|---|---|---|---
S | Se | $5.68$ | $4.1\pm 0.1$ | 525 | 3.0
Se | Te | $5.07$ | $7.0\pm 0.1$ | 471 | 3.0
S | Te | $0.00$ | $11.3\pm 0.1$ | 543 | 3.2

To circumvent this problem, we use relatively small (about 500 atoms) moiré
structures as targets for lattice-mismatched systems (See Table 1). Due to the
reduced symmetry of lattice-mismatched systems, the only valid targets are the
coordinates of all atoms of the moiré unit cell combined with the lattice
constant. However, using all atomic coordinates, i.e. three spatial dimensions
for each atom, renders the mesh adaptive search scheme for optimizing the KC-
parameters infeasible, as the number of target values greatly exceeds the
number of fitting parameters (Fig. 3, dashdotted green curve). As such, it is
necessary to reduce the number of target values. However, considering only the
three spatial coordinates of the metal atoms, thus reducing the target space
by one third, also yields sub-optimal KC-parameters (Fig. 3, dotted blue
curve). Lastly, optimizing only for the $z$-coordinates of the metal atoms,
which further reduces the target space by one third, results in a much better
fit (Fig. 3, dashed red curve). As such, we ultimately choose the
$z$-coordinates of the metal atoms and the lattice constant as target values
for lattice-mismatched systems, which yields satisfactory KC-parameters, as
discussed in Sec. III.

Figure 3:  $z$-coordinates of the metal atom in a 5.68∘ twisted WS2-MoSe2
bilayer along the long diagonal of the moiré unit cell, $r_{diag}$. The bottom
curvy layer corresponds to WS2, and the top more rigid layer corresponds to
MoSe2. The solid black curve is the DFT relaxed structure. The dashed red,
dotted blue, and dashdotted green curves are force-field relaxed with KC-
parameters (see text). (a) and (b) show R- and H-stackings, respectively.

###  II.3 The Kolmogorov-Crespi Potential

As mentioned previously, the KC potential, $V_{ij}$, is intended to model
interlayer effects between atom $i$ in one layer and atom $j$ in another, and
is given by

$\displaystyle V_{ij}$
$\displaystyle=\text{e}^{-\lambda(r_{ij}-z_{0})}\left[C+f(\rho_{ij})+f(\rho_{ji})\right]-A\left(\frac{r_{ij}}{z_{0}}\right)^{-6},$
$\displaystyle\rho_{ij}^{2}$
$\displaystyle=r_{ij}^{2}-(\bm{n}_{i}\bm{r}_{ij})^{2},$
$\displaystyle\rho_{ji}^{2}$
$\displaystyle=r_{ij}^{2}-(\bm{n}_{j}\bm{r}_{ij})^{2},$ $\displaystyle
f(\rho)$
$\displaystyle=\text{e}^{-(\rho/\delta)^{2}}\sum_{n=0}^{2}C_{2n}(\rho/\delta)^{2n}.$
(1)

$\bm{n}_{i}$ and $\bm{n}_{j}$ are the surface normals of the atom site $i$ and
$j$ in each layer. The choice of neighbors used to determine the surface
normals are the six nearest atoms in the respective strata (sublayer of the
monolayer). The last term of $V_{ij}$ contains the $r^{-6}$ vdW dependence,
and the first term has an exponentially decaying repulsion reflecting
interlayer wave-function overlap. The square bracket functions contain a
stacking dependent term, in contrast to e.g. the Lennard-Jones potential [31].
As seen, $V_{ij}$ leaves in total eight parameters to be fitted. As mentioned
in Ref. [28], it is possible to approximate $\bm{n}_{i,j}=\hat{z}$
corresponding to completely rigid layers, however, we do not make use of this
approximation in order to capture more accurately the corrugation caused by
the relaxation.

###  II.4 Computational details

We parametrize the potentials with different combinations of exchange
correlation plus vdW correction. We find that using PBE [42] from PseudoDojo
[43, 44] with Grimme’s DFT-D3 vdW correction [45] plus Becke-Johnson damping
[46] is best suited for parametrization. The structures are relaxed with
QuantumEspresso [47, 48] using a $k$-space density of $15\times 15$ ($1\times
1$) for high-symmetry (moiré) unit cells. DFT computations of moiré systems
are performed without spin-orbit coupling (SOC) to save computational
resources, since they are only used for comparing DFT to SW+KC relaxed
structures. We find that the lattice constant only converges at a cut-off
energy of 40 Ha in all cases. More importantly, the chosen cut-off energy
should be consistent between monolayers, untwisted bilayers and moiré
structured bilayers, when comparing DFT to SW+KC. We use the modified SW
implementation in LAMMPS for ease of use. For optimization in Dakota, we apply
a mesh adaptive direct search algorithm starting from the parameters presented
in Ref. [28].

##  III Results

For lattice-matched systems, i.e. homobilayers and heterobilayers having
identical chalgogen sites in both layers, which are developed by use of the
HSSCs, it is of high importance that the resulting structures can accurately
reproduce the electronic properties. In Fig. 4, a comparison between purely
DFT calculated parameters and SW+KC can be seen.

Figure 4: Comparison of $E_{g}$ at $\bm{K}_{\pm}$ in (a) and (d), interlayer
spacing ($d_{inter}$) in (b) and (e), and lattice constant $a_{0}$ in (c) and
(f) for the six HSSCs of a WS2 homobilayer and a MoS2-WS2 heterobilayer in
(a)-(c) (top panels) and (d)-(f) (bottom panels), respectively. DFT is marked
with black and SW+KC with red.

Note, that $E_{g}$ shown in Fig. 4(a,d) is the energetically lowest momentum-
conserving transition between the highest-lying valence band and the lowest-
lying conduction band, which occurs at the $\bm{K}_{\pm}$ points for all
stacking configurations and materials considered here. LAMMPS does not provide
$E_{g}$, instead this is calculated using DFT with the relaxed structures
generated by our SW+KC force-field method. Note, that for the purpose of
consistency, we adopt the notation that MoS2-WS2 implies that WS2 lies above
MoS2 with respect to $z$.

In the case of homobilayers, the maximum deviation of $E_{g}$ is 22 meV, and
occurs in the $\mathrm{H_{M}^{M}}$-stacking configuration. A similar maximum
deviation of 25 meV is seen for the heterobilayer, which occurs in the
$\mathrm{H_{M}^{M}}$-stacking configuration as well. For the remaining
lattice-matched structures, the deviations are of similar magnitude. Fig. 4
also demonstrates the high sensitivity of the bandgap with respect to changes
in the structural degrees of freedom.

Figure 5: Comparison of bandstructures and interlayer spacing profiles between
DFT in solid black and SW+KC in dotted red. WS2-WS2 homobilayer with
$\theta=6.0^{\circ}$ ($n_{atom}=546$) in (a),(b) and (c),(d) for R- and
H-stacking, respectively. MoS2-WS2 heterobilayer with $\theta=6.0^{\circ}$
($n_{atom}=546$) in (e),(f) and (g),(h) for R- and H-stacking, respectively.
WS2-MoSe2 lattice-mismatched heterobilayer with $\theta=5.1^{\circ}$
($n_{atom}=642$) in (i),(j) and (k),(l) for R- and H-stacking, respectively.
The bandstructures have the valence band maximum shifted to 0 in all cases,
and the Greek indices ($\gamma$, $\mu$ and $\kappa$) denote the high-symmetry
points of the moiré (mini) BZ (usually denoted $\Gamma$, M and K in the BZ of
the monolayer/untwisted bilayer). The interlayer spacing is interpolated and
plotted along the long diagonal of the unit cell.

Having established the SW+KC parameters of lattice-matched systems using the
HSSCs, we now tackle some larger moiré structures. As such, we use some
medium-scale moiré structures as benchmarks. Fig. 5 shows comparisons between
DFT- and SW+KC-relaxed band structures and interlayer spacing profiles for
different material cases. Greek indices denote the corners of the mini
Brillouin zone (BZ) associated with a moiré structured bilayer. The interlayer
spacing is plotted along the long diagonal of the unit cell (see Fig. 1),
which has a length of $\sqrt{3}m_{0}$, where $m_{0}$ denotes the moiré cell
lattice constant. Fig. 5(a)-(d) displays the case of a WS2-WS2 homobilayer
twisted at 6.0∘. For R-stacking ((a),(b)), the bands are well represented
using our SW+KC relaxed structure with only a 13 meV decrease of the bandgap,
likely due to the slight interlayer spacing profile discrepancy. In the case
of H-stacking ((c)-(d)), the same applies except the bandgap is a mere 5 meV
larger compared to the DFT relaxed structure results. In Fig. 5(e)-(h), the
case of a MoS2-WS2 heterobilayer with a twist angle of 6.0∘ is shown. For both
R- and H-stacking ((e),(f) and (g),(h) respectively), an excellent agreement
is obtained between DFT- and SW+KC-relaxed structures in terms of band
character. For the higher lying conduction bands around the $\bm{\gamma}$
point, there is only a 10 meV discrepancy. We again attribute this to the
slightly decreased interlayer spacing profiles of SW+KC in both cases, as seen
in Fig. 5(f),(h).

For lattice-mismatched systems, the optimization of the KC parameters was
performed for all possible combinations of metal and chalcogen atoms, as
explained in the methodology section (see Table 1). A good agreement is
obtained between the DFT- and SW+KC-relaxed structures for all lattice-
mismatched cases. For the sake of brevity, only the case of a WS2-MoSe2
heterobilayer rotated at 5.1∘ is shown in Fig. 5(i)-(l). For R-stacking
((i),(j)), the highest lying valence band is only 7 meV higher than the DFT
value at the $\bm{\kappa}$-point. The lowest-lying conduction band is only 6
meV above the DFT one. In general, we see small discrepancies between the
valence and conduction bands for the DFT and SW+KC-relaxed structures below 20
meV. For H-stacking ((k),(l)), the valence bands are well described except for
a 4 meV discrepancy of the highest-lying valence band near the $\bm{\kappa}$
point.

In general, we note that the slight difference in bandgap and band curvature
between DFT and our SW+KC-relaxed moiré structures arise from small
inaccuracies in the interlayer spacing profiles. Note, that this is not always
the case with the KC-parameters presented by Ref. [28], where the binding
energy was the target property. We also find that the accuracy of our lattice-
matched SW+KC parameters reduce with growing twist angle. This is expected,
since we fit to the untwisted HSSCs, which are not well represented in moiré
structures with such low periodicity. Conversely, the parameters are expected
to have better accuracy with decreasing twist angle. For angles below
$~{}3^{\circ}$, where large-scale atomic reconstruction starts to appear, the
accuracy of methodology is still maintained and most properties are well
captured, including the atomic reconstructions, as discussed in Sec. V.

##  IV Approximating Moiré Potentials

A defining feature of two-dimensional lattice-matched moiré structures is the
spatial variation of local stacking order across the structure, leading to
variation of local properties. Many combinations of TMDs possess type-II band
alignment [49, 50], and as such, the variation of the local bandgap at
$\bm{K}_{\pm}$ across the structure will, for many purposes, describe the
interlayer moiré potential [27, 51, 17]. However, it is worth mentioning, that
in the case of a large lattice-mismatch between the constituting layers,
developing such a potential becomes non-trivial.

We propose two interpolation-based methods for calculating the interlayer
moiré potential of lattice-matched systems. Moreover, any electronic property
that can be identified locally, can be accessed in the moiré structure
directly with these two methods, e.g. variation of the VBM, CBM etc. In both
methods, the moiré supercell is subdivided into small units the size of the
monolayer unit cell, for which local properties can be calculated. The first
method, which we call the high-symmetry interpolation method (HSIM), is based
on the local high-symmetry stacking character \- a geometrical quantity that
measures the similarity between the local stacking configuration within the
moiré cell and the HSSCs. Being based only on the six HSSCs, computing the
DFT-properties is fast and allows for high-throughput computations. It also
allows for easy visualization of reconstructed domains. The second method,
which we call the grid based interpolation method (GBIM), relies on computing
the local properties using DFT not only for the HSSCs, but also every local
stacking configuration in between, which can then be interpolated over the
moiré supercell. In principle, this scheme is more precise, since it relies
less on interpolation and more on ab initio calculations. However, it is time
consuming, as many DFT computations using different in-plane displacements and
interlayer spacing are needed. In what follows, both methods are explained in
detail and case studies are shown.

###  IV.1 High-symmetry interpolation method (HSIM)

For every metal site in one layer,
$\bm{\rho}_{\mathrm{M},i}=(x_{\mathrm{M},i},y_{\mathrm{M},i})$, we find the
transverse distance to the closest metal site in the adjacent layer, e.g.
$d_{\mathrm{M},i}^{\mathrm{M}}=\mathrm{min}(|\bm{\rho}_{\mathrm{M},i}-\bm{\rho}_{\mathrm{M},j}|)$,
where $j$ runs through every metal site in the adjacent layer (see Fig. 6).
The largest distance possible is $a_{0}/\sqrt{3}$.

Figure 6:  Close-up of an R-stacked lattice-matched moiré structure for
$\theta=6.0^{\circ}$. For a metal site $i$ in one layer, the nearest
transverse metal site, $\mathrm{M},j$, chalcogen site, $\mathrm{X},j$, and
hexagonal center, $\mathrm{H},j$, in the adjacent layer is seen.

As such, we can define the parameter
$c_{\mathrm{M},i}^{\mathrm{M}}=1-\sqrt{3}d_{\mathrm{M},i}^{\mathrm{M}}/a_{0}$,
which is unity for perfectly aligned metal atoms, e.g. $\mathrm{R_{X}^{X}}$\-
and $\mathrm{H_{M}^{M}}$-stacking, and zero for the remaining HSSCs. Eight
analogous parameters can be developed, e.g.

$\\{c_{\mathrm{S}_{1},i}^{\mathrm{S}_{2}}(\bm{\rho}_{\mathrm{S}_{1},i})\quad\mathrm{for}\quad\mathrm{S}_{1},\mathrm{S}_{2}\in\\{\mathrm{M},\mathrm{X},\mathrm{H}\\}\\},$

where $\mathrm{X}$ and $\mathrm{H}$ denote chalcogen sites and hexagonal
centers, respectively. For the purpose of consistency, it is assumed that
$\mathrm{S}_{2}$ lies above $\mathrm{S}_{1}$ with respect to $z$.
$\\{c_{\mathrm{S}_{1},i}^{\mathrm{S}_{2}}\\}$ is then interpolated on a skewed
grid that spans the moiré unit cell. Stacking coefficients are now found as

$\displaystyle
C_{\mathrm{R_{X}^{X}}}=c_{\mathrm{M}}^{\mathrm{M}}c_{\mathrm{X}}^{\mathrm{X}}c_{\mathrm{H}}^{\mathrm{H}}$
$\displaystyle,\quad
C_{\mathrm{H_{X}^{M}}}=c_{\mathrm{X}}^{\mathrm{M}}c_{\mathrm{M}}^{\mathrm{X}}c_{\mathrm{H}}^{\mathrm{H}},$
$\displaystyle
C_{\mathrm{R_{X}^{M}}}=c_{\mathrm{X}}^{\mathrm{M}}c_{\mathrm{H}}^{\mathrm{X}}c_{\mathrm{M}}^{\mathrm{H}}$
$\displaystyle,\quad
C_{\mathrm{H_{X}^{X}}}=c_{\mathrm{H}}^{\mathrm{M}}c_{\mathrm{X}}^{\mathrm{X}}c_{\mathrm{M}}^{\mathrm{H}},$
$\displaystyle
C_{\mathrm{R_{M}^{X}}}=c_{\mathrm{H}}^{\mathrm{M}}c_{\mathrm{M}}^{\mathrm{X}}c_{\mathrm{X}}^{\mathrm{H}}$
$\displaystyle,\quad
C_{\mathrm{H_{M}^{M}}}=c_{\mathrm{M}}^{\mathrm{M}}c_{\mathrm{H}}^{\mathrm{X}}c_{\mathrm{X}}^{\mathrm{H}}.$

Finally, the stacking coefficients are normalized such that
$\sum_{n}C_{n}({\bm{\rho}})=1$, where $n$ spans the HSSCs. $C_{n}$ is seen in
Fig. 7 for R-stacking. The $C_{n}$ with
$n\in\\{\mathrm{H_{X}^{M}},\mathrm{H_{X}^{X}},\mathrm{H_{M}^{M}}\\}$ are all 0
in this case.

Figure 7: Variation of $C_{\mathrm{R_{X}^{X}}}$, $C_{\mathrm{R_{X}^{M}}}$ and
$C_{\mathrm{R_{M}^{X}}}$ across an R-stacked MoS2-WS2 bilayer with
$\theta=4.41^{\circ}$ in (a), (b) and (c), respectively. The remaining
coefficients are 0. The structure was relaxed using LAMMPS.

The next step is finding the interlayer spacing profile,
$d_{inter}(\bm{\rho})$, where $\bm{\rho}=(x,y)$. Using the variation of
$E_{g}$, $E_{g}(\bm{\rho})$, as an example, it can be seen that

$E_{g}(\bm{\rho})=\sum_{n}C_{n}(\bm{\rho})E_{g}(n,d_{inter}(\bm{\rho})),$ (2)

assuming the variation of $E_{g}$ with $d_{inter}$ is known for all HSSCs.
Assuming that every subcell of the moiré structure can be described by a
superposition of HSSCs is an approximation, but has the benefit of easy
visualization of domains, as seen in Fig. 7. It shows great accuracy and
$E_{g}(n,d_{inter})$ can be extracted within few calculations, making it quite
fast to implement for all lattice-matched systems.

###  IV.2 Grid based interpolation method (GBIM)

A more general implementation can be developed by using the untwisted bilayer
with a transverse shift $\bm{\rho}_{s}=(x_{s},y_{s})$ between the layers,
where $\bm{\rho}_{s}=0$ corresponds to either $\mathrm{R_{X}^{X}}$\- or
$\mathrm{H_{M}^{M}}$-stacking. We calculate
$E_{g}^{\textrm{BL}}(\bm{\rho}_{s},d_{inter})$, where $\bm{\rho}_{s}$ is the
transverse distance between metal sites in each layer. Then, for a given
lattice-matched moiré system, for metal site $i$ in one layer, we can find the
vector $\bm{\rho}_{i}=\bm{\rho}_{\mathrm{M},j}-\bm{\rho}_{\mathrm{M},i}$,
where $j$ denotes the index of the nearest metal site in the adjacent layer.
Then, the value of $E_{g}$ at metal site $i$ is simply

$E_{g}(\bm{\rho}_{\mathrm{M},i})=E_{g}^{\textrm{BL}}(\bm{\rho}_{i},d_{inter}(\bm{\rho}_{\mathrm{M},i})).$
(3)

Note, that $\bm{\rho}_{i}$ should be adjusted relative to the rotation of the
individual layers, since the layers will likely be slightly angled compared to
the systems used in computing $E_{g}^{\textrm{BL}}(\bm{\rho}_{s},d_{inter})$.
Finally, $E_{g}$ is interpolated over the entire moiré unit cell.

In principle, the GBIM should be more accurate than the HSIM, but is also
computationally more expensive. We use twelve steps for $x_{s}$ and $y_{s}$
combined with sixteen increments for $d_{inter}$ when tabulating
$E_{g}^{\textrm{BL}}(\bm{\rho}_{s},d_{inter})$. This translates to 4608
separate DFT calculations to cover R- and H-stacking for one material, whereas
the HSIM needs only 96. In Fig. 8, a comparison between the HSIM and the GBIM
can be seen for the variation of $E_{g}$ in MoS2-WS2 with
$\theta=4.41^{\circ}$. At the high-symmetry points, both methods yield the
same value as expected, but the HSIM is slightly inaccurate in between.

Figure 8: Variation of $E_{g}$ along the long diagonal of an R- and H-stacked
MoS2-WS2 bilayer with $\theta=4.41^{\circ}$ in (a) and (b), respectively. The
solid black and dashed red curves represent bandgap variation found using the
GBIM and HSIM, respectively.

##  V Atomic Reconstruction and
Energetic Landscape

As mentioned, the energetic landscape of 2D moiré structures is constituted by
three codependent factors: the local stacking arrangement, the associated
interlayer spacing, and the atomic reconstruction. Often, the latter two,
being relaxation effects, are not considered in simulations [52, 51, 17, 53,
54, 55, 56, 57, 58, 59, 15, 60, 61, 62, 11], but can be managed with SW+KC
force-field relaxation.

For MoS2-WS2, which possesses type-II band alignment [49, 50], the interlayer
moiré potential is often described as the spatial variation of the local
bandgap at $\bm{K}_{\pm}$. In Fig. 9, the variation of
$E_{g}-\overline{E}_{g}$ across an R- and H-stacked MoS2-WS2 bilayer with
$\theta=1.01^{\circ}$ is seen, where $\overline{E}_{g}$ is the mean value
across the unit cell. In the rigidly twisted case, the average interlayer
spacing of the three R- or H-stacked HSSCs are used as interlayer spacing for
Fig. 9(a) and Fig. 9(d), respectively.

Figure 9: Variation of $E_{g}$ at $\bm{K}_{\pm}$ across a 1.01∘ twisted
R-stacked MoS2-WS2 bilayer without and with relaxation effects in the (a),(d)
and (b),(e), respectively. Comparison between the two cases along the long
diagonal of the moiré unit cells in (c),(f). (a)-(c) and (d)-(e) represent R-
and H-stackings, respectively.

The discrepancy between modeling the potential with- and without relaxation
effects is apparent. In the case of R-stacking, which has larger potential
depth than H-stacking, the depth with- and without relaxation effects are here
estimated to be 80 meV and 135 meV, respectively. For H-stacking, these
numbers are much lower, namely 10 meV and 28 meV for rigid and relaxed,
respectively. As a consequence, phenomena such as exciton trapping may be
realized more easily in R-stacked systems than H-stacked analogs.
Interestingly, for H-stacking, the minimum of the potential resides in the
$\mathrm{H^{X}_{X}}$-domain post-relaxation as opposed to the
$\mathrm{H^{M}_{M}}$-domain pre-relaxation. Finally, the effect of atomic
reconstruction also greatly changes the relative widths of the potential
wells, resulting in a sharper and more well-defined potential. We conclude
that atomic reconstruction significantly alters the range of $\theta$ in which
exciton trapping occurs.

In Table 2, the interlayer binding energy per atom, $E_{b}$, found as
$E_{b}=(E_{\mathrm{MoS_{2}}\text{-}\mathrm{WS_{2}}}-E_{\mathrm{MoS_{2}}}-E_{\mathrm{WS_{2}}})/6$
is shown, where $E_{\mathrm{MoS_{2}}\text{-}\mathrm{WS_{2}}}$ is the total
energy of the untwisted bilayer system, and $E_{\mathrm{MoS_{2}}}$ and
$E_{\mathrm{WS_{2}}}$ denote the total energies of constituting monolayers
found separately. As mentioned, the discrepancy in $E_{b}$ between DFT and
SW+KC is expected, since this was not the target property during development
of our KC parameters.

Table 2:  Binding energy of MoS2-WS2 in the six high-symmetry stacking configurations from DFT and from SW+KC. $E_{b}$ (meV) | $\mathrm{R^{X}_{X}}$ | $\mathrm{R^{X}_{M}}$ | $\mathrm{R^{M}_{X}}$ | $\mathrm{H^{M}_{X}}$ | $\mathrm{H^{X}_{X}}$ | $\mathrm{H^{M}_{M}}$
---|---|---|---|---|---|---
DFT | -21.8 | -34.3 | -34.5 | -34.6 | -22.4 | -31.6
SW+KC | -25.4 | -44.4 | -44.6 | -44.7 | -28.9 | -38.4

For R-stacking, the nearly identical $E_{b}$ of the $\mathrm{R^{X}_{M}}$\- and
$\mathrm{R^{M}_{X}}$-configurations facilitates a simultaneous growth of these
domains with decreasing $\theta$ (i.e. large moiré unit cells), while the
opposite is true for $\mathrm{R^{X}_{X}}$, explaining the formation of a mesh
of triangular domains, as seen in Fig. 9(b). For H-stacking, the
$\mathrm{H^{M}_{X}}$-configuration is energetically favorable, resulting in
hexagonal domains with decreasing $\theta$. The $\mathrm{H^{M}_{M}}$-like
domains shrink slower than those associated with $\mathrm{H^{X}_{X}}$, as seen
from the associated $E_{b}$ (see also [63]).

Lastly, $E_{b}(\theta)$ can be considered in order to access the formation
rate of domains. Using the HSIM, the variation of the local $E_{b}$ across a
moiré unit cell can be approximated, and the mean can be used to approximate
$E_{b}$ of the moiré unit cell, albeit neglecting the effects of strain
imposed by atomic reconstruction and corrugation from the varying interlayer
spacing. In the case of pure DFT, $E_{b}$, is found directly as

$E_{b}=E_{\text{moiré}}-(E_{\mathrm{MoS_{2}}}+E_{\mathrm{WS_{2}}})/2,$ (4)

where all energies are divided by the number of atoms, and $E_{\text{moiré}}$
denotes the total energy per atom of the moiré structure. However, $E_{b}$ has
contributions from the strain imposed by layer corrugation and atomic
reconstruction. The energy associated with these effects is denoted
$E_{\text{corr}}$ and is not captured by the HSIM. Instead, the $E_{b}$ found
by the HSIM should be compared to

$E_{b}-E_{\text{corr}}=E_{\text{moiré}}-(E_{\mathrm{MoS_{2},\text{moiré}}}+E_{\mathrm{WS_{2},\text{moiré}}})/2,$
(5)

where $E_{\mathrm{MoS_{2},\text{moiré}}}$ and
$E_{\mathrm{WS_{2},\text{moiré}}}$ denote the total energy per atom for the
corrugated and reconstructed constituting monolayers. This is computed in
separate DFT calculations having half the number of atoms as the moiré
structure they constitute.

Figure 10: $E_{b}(\theta)$ for MoS2-WS2. (a) and (b) are for R- and
H-stacking, respectively. Black points are found using the HSIM on the SW+KC-
relaxed structures, but with DFT-based parametrization of the HSIM as seen in
Eq. (2).

With SW+KC, $E_{b}$ is found analogously to Eq. (4), but $E_{\text{corr}}$ is
found directly by comparing the energy of the SW-potential in the two layers
to that of the constituting rigid monolayers. The variation of these
quantities with $\theta$ is seen in Fig. 10.

A common feature for all energy scales in Fig. 10 is the tendency towards the
value of the stable configurations for $\theta\rightarrow 0$. For vanishing
$\theta$, the relative size of the domain walls becomes negligible. As such,
$E_{\text{corr}}$ should vanish in the limit of vanishing $\theta$. The faster
convergence towards the $E_{b}$ of $\mathrm{R_{X}^{M}}$/$\mathrm{R_{M}^{X}}$
for R-stacking indicates that the triangular domains form more rapidly with
decreasing $\theta$ compared to the hexagonal $\mathrm{H_{X}^{M}}$-domains of
H-stacking. Although the values of $E_{b}^{\mathrm{SW+KC}}$ and
$E_{\text{corr}}^{\mathrm{SW+KC}}$ may appear off scale, they illustrate the
tendencies faithfully. Additionally, the graph of
$E_{b}^{\mathrm{DFT}}-E_{\text{corr}}^{\mathrm{DFT}}$ serves as a benchmark,
showing that the HSIM has accuracy within the 0.5 meV range, and further that
$E_{b}$ of SW+KC relaxed structures can be recovered to agree with DFT.

Fig. 11 shows the mean of the stacking coefficients $C_{n}$ over the moiré
unit cell of MoS2-WS2, which can be computed using the HSIM as described in
Sec. IV.1. $C_{n}$ represents the normalized contributions of the different
stacking configurations to the fully relaxed (reconstructed) moiré structure.

Figure 11: The mean of $C_{n}$ over the unit cell with respect to the twist
angle, $\theta$, of MoS2-WS2 found using SW+KC and the HSIM. (a) correspond to
R-stacking such that only $\mathrm{R_{X}^{X}}$, $\mathrm{R_{M}^{X}}$ and
$\mathrm{R_{X}^{M}}$ are non-zero and vice versa for (b). For each stacking
configuration, the corresponding $E_{b}$ of untwisted MoS2-WS2 is shown in
meV.

For R-stacking (H-stacking), the three possible domains are:
$\mathrm{R_{M}^{X}}$ (green), $\mathrm{R_{X}^{X}}$ (red), $\mathrm{R_{X}^{M}}$
(blue) ($\mathrm{H_{M}^{M}}$ (green), $\mathrm{H_{X}^{X}}$ (red),
$\mathrm{H_{X}^{M}}$ (blue)). For larger angles, the fraction of the unit cell
area occupied by each of the three domains is about $1/3$ for both R- and
H-stacking. At an angle of 1∘ the structure for R-stacking (Fig. 11(a)) is
already reconstructed in such a way that the energetically less favorable
$\mathrm{R_{X}^{X}}$ (red) domains represent only 2,5% of the structure. Both
$\mathrm{R_{M}^{X}}$ (green) and $\mathrm{R_{X}^{M}}$ (blue) domains are
energetically equivalent, and hence, occupy roughly 50% of the structure in
the limit of small $\theta$. For H-stacking (Fig. 11(b)), at the same angle of
1∘, the less favorable $\mathrm{H_{X}^{X}}$ and $\mathrm{H_{M}^{M}}$ domains
have significantly reduced contributions compared to the favorable
$\mathrm{H_{X}^{M}}$ region, but $\mathrm{H_{M}^{M}}$ still represents 20% of
overall structure.

Fig. 11 allows us to draw quantitative conclusions on the angle dependence of
the reconstruction effect. Indeed, neglecting reconstructions would lead to a
constant equal proportion of all three coexisting stackings (dotted lines in
Fig. 11). In the case of R-stacking the reconstruction is nearly complete at
an angle of 1∘, i.e., the moire structure is made of basically two type of low
energy domains ($\mathrm{R_{M}^{X}}$ (green), $\mathrm{R_{X}^{M}}$ (blue))
separated by a very narrow $\mathrm{R_{X}^{X}}$(red) energetically unfavorable
domain. For H-stacking at 1∘, the less favorable $\mathrm{H_{M}^{M}}$ (green)
domain still covers 15-20% of the area.

Fig. 12 shows the same graph as Fig. 11 for the remaining eight lattice-
matched materials. Generally, all R-stacked materials (left panels of Fig. 12)
display a simultaneous growth of $\mathrm{R_{M}^{X}}$ and $\mathrm{R_{X}^{M}}$
with decreasing $\theta$ except for MoTe2-WTe2, which can be attributed to the
discrepancy in $E_{b}$ for these stacking configurations. We conclude that for
both stackings and all materials considered here, except for MoSe2-MoSe2, that
atomic reconstruction becomes especially prominent below an angle of 4-5∘. For
MoSe2-MoSe2, atomic reconstruction occurs for angles below 6-7∘.

Figure 12: The mean of $C_{n}$ over the unit cell with respect to the twist
angle, $\theta$, of eight lattice-matched bilayers found using SW+KC and the
HSIM. Left panels correspond to R-stacking such that only
$\mathrm{R_{X}^{X}}$, $\mathrm{R_{M}^{X}}$ and $\mathrm{R_{X}^{M}}$ are non-
zero and vice versa for right panels. For each material and stacking
configuration, the corresponding $E_{b}$ of the untwisted bilayer is shown in
meV.

##  VI Conclusion

In conclusion, we have shown the dramatic consequences of incorporating
relaxation effects on the interlayer moiré potential of MoS$2$-WS2. For
R-stacking, this becomes about twice as deep at about 135 meV, and, for small
angles, much wider. For H-stacking, the potential depth is nearly tripled,
however, the width of the potential minima is still narrow, since it
corresponds to the energetically unfavorable
$\mathrm{H^{X}_{X}}$-configuration. Moreover, we have quantified the formation
rate of domains due to atomic reconstruction for nine lattice-matched TMD
moiré systems, and conclude that, in general, atomic reconstruction becomes
prominent for $\theta$ smaller than 4-5∘, but does so in a continuous manner.

Furthermore, we have presented a methodology for developing KC-parameters for
lattice-matched and -mismatched systems, and have developed such parameters
for TMD moiré heterostructures. The method shows excellent agreement between
DFT-calculated structural parameters and SW+KC-relaxed ones, which is further
reflected in the bandstructure and the interlayer binding energy with twist
angle dependence. The force-field parameters along with a variety of relaxed
structures can be found via Ref. [38]. We have further shown two methods for
capturing moiré induced fluctuations of local properties in lattice-matched
systems that do not require extensive ab initio treatment. These methods allow
for visualization of the importance of relaxation effects and further serve as
a first step in developing accurate moiré potentials. However, further
investigation is required to develop analogous tools for lattice-mismatched
moiré structures.

In summary, starting from the force-field model, it is now possible to tackle
excited state physics incorporating relaxation effects i.e. layer corrugation
and atomic reconstruction. For models such as tight-binding, this was not
possible before, and for ab initio studies, the cumbersome first step of
relaxation can be skipped, thus saving computational resources and time.
Furthermore, a thorough dissection of the formation rate of domains with
decreasing angle is required to gain quantitative insight into the mechanisms
behind it.

###### Acknowledgements.

The project is supported by the Deutsche Forschungsgemeinschaft (DFG) within
the Priority Program SPP2244 2DMP and by the Cluster of Excellence “Advanced
Imaging of Matter” of the Deutsche Forschungsgemeinschaft (DFG) – EXC 2056 –
project ID 390715994.

## References

  * Liu  _et al._ [2016] Y. Liu, N. O. Weiss, X. Duan, H.-C. Cheng, Y. Huang, and X. Duan, Nature Reviews Materials 1, 16042 (2016).
  * Novoselov  _et al._ [2016] K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, Science 353, aac9439 (2016), https://www.science.org/doi/pdf/10.1126/science.aac9439 .
  * Gibertini  _et al._ [2019] M. Gibertini, M. Koperski, A. F. Morpurgo, and K. S. Novoselov, Nature Nanotechnology 14, 408 (2019).
  * Yao  _et al._ [2021] X. Yao, Y. Wang, and S. Dong, International Journal of Modern Physics B 35, 2130004 (2021), https://doi.org/10.1142/S0217979221300048 .
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011), https://www.pnas.org/doi/pdf/10.1073/pnas.1108174108 .
  * Al-Ani  _et al._ [2022] I. A. M. Al-Ani, K. As’ham, O. Klochan, H. T. Hattori, L. Huang, and A. E. Miroshnichenko, Journal of Optics 24, 053001 (2022).
  * Naik and Jain [2018] M. H. Naik and M. Jain, Phys. Rev. Lett. 121, 266401 (2018).
  * Choi  _et al._ [2020] J. Choi, W.-T. Hsu, L.-S. Lu, L. Sun, H.-Y. Cheng, M.-H. Lee, J. Quan, K. Tran, C.-Y. Wang, M. Staab, K. Jones, T. Taniguchi, K. Watanabe, M.-W. Chu, S. Gwo, S. Kim, C.-K. Shih, X. Li, and W.-H. Chang, Science Advances 6, eaba8866 (2020), https://www.science.org/doi/pdf/10.1126/sciadv.aba8866 .
  * Baek  _et al._ [2020] H. Baek, M. Brotons-Gisbert, Z. X. Koong, A. Campbell, M. Rambach, K. Watanabe, T. Taniguchi, and B. D. Gerardot, Science Advances 6, eaba8526 (2020), https://www.science.org/doi/pdf/10.1126/sciadv.aba8526 .
  * Guo  _et al._ [2020] H. Guo, X. Zhang, and G. Lu, Science Advances 6, eabc5638 (2020), https://www.science.org/doi/pdf/10.1126/sciadv.abc5638 .
  * Wu  _et al._ [2018a] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. LDowett. 121, 026402 (2018a).
  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Nature 579, 353 (2020).
  * Pan  _et al._ [2020a] H. Pan, F. Wu, and S. Das Sarma, Phys. Rev. B 102, 201104 (2020a).
  * Slagle and Fu [2020] K. Slagle and L. Fu, Phys. Rev. B 102, 235423 (2020).
  * Pan  _et al._ [2020b] H. Pan, F. Wu, and S. Das Sarma, Phys. Rev. Research 2, 033087 (2020b).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl, and F. Wang, Nature 579, 359 (2020).
  * Götting  _et al._ [2022] N. Götting, F. Lohof, and C. Gies, Phys. Rev. B 105, 165419 (2022).
  * Rosenberger  _et al._ [2020] M. R. Rosenberger, H.-J. Chuang, M. Phillips, V. P. Oleshko, K. M. McCreary, S. V. Sivaram, C. S. Hellberg, and B. T. Jonker, ACS Nano 14, 4550 (2020), https://doi.org/10.1021/acsnano.0c00088 .
  * Carr  _et al._ [2018] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras, Phys. Rev. B 98, 224102 (2018).
  * Maity  _et al._ [2021] I. Maity, P. K. Maiti, H. R. Krishnamurthy, and M. Jain, Phys. Rev. B 103, L121102 (2021).
  * Enaldiev  _et al._ [2020] V. V. Enaldiev, V. Zólyomi, C. Yelgel, S. J. Magorrian, and V. I. Fal’ko, Phys. Rev. Lett. 124, 206101 (2020).
  * Li  _et al._ [2021a] H. Li, S. Li, M. H. Naik, J. Xie, X. Li, J. Wang, E. Regan, D. Wang, W. Zhao, S. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, A. Zettl, S. G. Louie, F. Wang, and M. F. Crommie, Nature Materials 20, 945 (2021a).
  * Weston  _et al._ [2020] A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, J. Zultak, D. Hopkinson, A. Barinov, T. H. Bointon, A. Kretinin, N. R. Wilson, P. H. Beton, V. I. Fal’ko, S. J. Haigh, and R. Gorbachev, Nature nanotechnology 15, 592—597 (2020).
  * Shabani  _et al._ [2021] S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy, Nature Physics 17, 720 (2021).
  * Sung  _et al._ [2020] J. Sung, Y. Zhou, G. Scuri, V. Zólyomi, T. I. Andersen, H. Yoo, D. S. Wild, A. Y. Joe, R. J. Gelly, H. Heo, S. J. Magorrian, D. Bérubé, A. M. M. Valdivia, T. Taniguchi, K. Watanabe, M. D. Lukin, P. Kim, V. I. Fal’ko, and H. Park, Nature Nanotechnology 15, 750 (2020).
  * Liu  _et al._ [2021] Y. Liu, C. Zeng, J. Yu, J. Zhong, B. Li, Z. Zhang, Z. Liu, Z. M. Wang, A. Pan, and X. Duan, Chem. Soc. Rev. 50, 6401 (2021).
  * Geng  _et al._ [2020] W. T. Geng, V. Wang, Y. C. Liu, T. Ohno, and J. Nara, The Journal of Physical Chemistry Letters 11, 2637 (2020), https://doi.org/10.1021/acs.jpclett.0c00605 .
  * Naik  _et al._ [2019] M. H. Naik, I. Maity, P. K. Maiti, and M. Jain, The Journal of Physical Chemistry C 123, 9770 (2019).
  * Stillinger and Weber [1985] F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).
  * Jiang and Zhou [2017] J.-W. Jiang and Y.-P. Zhou,  _Handbook of Stillinger-Weber Potential Parameters for Two-Dimensional Atomic Crystals_ (IntechOpen, Rijeka, 2017) chapter 1.
  * Kolmogorov and Crespi [2005] A. N. Kolmogorov and V. H. Crespi, Phys. Rev. B 71, 235415 (2005).
  * Kolmogorov and Crespi [2000] A. N. Kolmogorov and V. H. Crespi, Phys. Rev. Lett. 85, 4727 (2000).
  * Maaravi  _et al._ [2017] T. Maaravi, I. Leven, I. Azuri, L. Kronik, and O. Hod, The Journal of Physical Chemistry C 121, 22826 (2017), https://doi.org/10.1021/acs.jpcc.7b07091 .
  * Leven  _et al._ [2016] I. Leven, T. Maaravi, I. Azuri, L. Kronik, and O. Hod, Journal of Chemical Theory and Computation 12, 2896 (2016), https://doi.org/10.1021/acs.jctc.6b00147 .
  * Leven  _et al._ [2014] I. Leven, I. Azuri, L. Kronik, and O. Hod, The Journal of Chemical Physics 140, 104106 (2014), https://doi.org/10.1063/1.4867272 .
  * Ouyang  _et al._ [2018] W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod, Nano Letters 18, 6009 (2018).
  * Ouyang  _et al._ [2021] W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. Tkatchenko, L. Kronik, M. Urbakh, and O. Hod, Journal of Chemical Theory and Computation 17, 7237 (2021), https://doi.org/10.1021/acs.jctc.1c00782 .
  * [38] Force-field parameters, https://www.chemie.uni-hamburg.de/en/institute/pc/arbeitsgruppen/bester/downloads.html.
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comp. Phys. Comm. 271, 108171 (2022).
  * Eldred  _et al._ [2009] M. Eldred, K. Dalbey, W. Bohnhoff, B. Adams, L. Swiler, P. Hough, D. Gay, J. Eddy, and K. Haskell, N/A (2009).
  * Shi  _et al._ [2013] H. Shi, H. Pan, Y.-W. Zhang, and B. I. Yakobson, Phys. Rev. B 87, 155304 (2013).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * van Setten  _et al._ [2018] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, Computer Physics Communications 226, 39 (2018).
  * Hamann [2013] D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, The Journal of Chemical Physics 132, 154104 (2010), https://doi.org/10.1063/1.3382344 .
  * Becke and Johnson [2006] A. D. Becke and E. R. Johnson, The Journal of Chemical Physics 124, 221101 (2006), https://doi.org/10.1063/1.2213970 .
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Özçelik  _et al._ [2016] V. O. Özçelik, J. G. Azadani, C. Yang, S. J. Koester, and T. Low, Phys. Rev. B 94, 035125 (2016).
  * Zhang  _et al._ [2016] C. Zhang, C. Gong, Y. Nie, K.-A. Min, C. Liang, Y. J. Oh, H. Zhang, W. Wang, S. Hong, L. Colombo, R. M. Wallace, and K. Cho, 2D Materials 4, 015026 (2016).
  * Wu  _et al._ [2018b] F. Wu, T. Lovorn, and A. H. MacDonald, Phys. Rev. B 97, 035306 (2018b).
  * Wu  _et al._ [2017] F. Wu, T. Lovorn, and A. H. MacDonald, Phys. Rev. Lett. 118, 147401 (2017).
  * Rademaker [2022] L. Rademaker, Phys. Rev. B 105, 195428 (2022).
  * Yu  _et al._ [2017] H. Yu, G.-B. Liu, J. Tang, X. Xu, and W. Yao, Science Advances 3, e1701696 (2017), https://www.science.org/doi/pdf/10.1126/sciadv.1701696 .
  * Hichri  _et al._ [2021] A. Hichri, T. Amand, and S. Jaziri, Phys. Rev. Materials 5, 114002 (2021).
  * Lagoin and Dubin [2021] C. Lagoin and F. m. c. Dubin, Phys. Rev. B 103, L041406 (2021).
  * 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, Nature 567, 71 (2019).
  * Li  _et al._ [2021b] Z. Li, X. Lu, D. F. Cordovilla Leon, Z. Lyu, H. Xie, J. Hou, Y. Lu, X. Guo, A. Kaczmarek, T. Taniguchi, K. Watanabe, L. Zhao, L. Yang, and P. B. Deotare, ACS Nano 15, 1539 (2021b), https://doi.org/10.1021/acsnano.0c08981 .
  * Zheng  _et al._ [2021] H. Zheng, D. Zhai, and W. Yao, 2D Materials 8, 044016 (2021).
  * Zhai and Yao [2020] D. Zhai and W. Yao, Phys. Rev. Materials 4, 094002 (2020).
  * Brem  _et al._ [2020] S. Brem, C. Linderälv, P. Erhart, and E. Malic, Nano Letters 20, 8534 (2020), https://doi.org/10.1021/acs.nanolett.0c03019 .
  * Wu  _et al._ [2019] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Phys. Rev. Lett. 122, 086402 (2019).
  * Halbertal  _et al._ [2021] D. Halbertal, N. R. Finney, S. S. Sunku, A. Kerelsky, C. Rubio-Verdú, S. Shabani, L. Xian, S. Carr, S. Chen, C. Zhang, L. Wang, D. Gonzalez-Acevedo, A. S. McLeod, D. Rhodes, K. Watanabe, T. Taniguchi, E. Kaxiras, C. R. Dean, J. C. Hone, A. N. Pasupathy, D. M. Kennes, A. Rubio, and D. N. Basov, Nature Communications 12, 242 (2021).