# Relaxation effects in twisted bilayer molybdenum disulfide: structure,
stability, and electronic properties

Florian M. Arnold  TU Dresden, Theoretical Chemistry, Bergstr. 66c, 01062
Dresden, Germany  Alireza Ghasemifard  TU Dresden, Theoretical Chemistry,
Bergstr. 66c, 01062 Dresden, Germany  Agnieszka Kuc<EMAIL_ADDRESS>Helmholtz-
Zentrum Dresden-Rossendorf, Abteilung Ressourcenökologie,Forschungsstelle
Leipzig, Leipzig, Germany  Jens Kunstmann  TU Dresden, Theoretical Chemistry,
Bergstr. 66c, 01062 Dresden, Germany  Thomas Heine<EMAIL_ADDRESS>TU Dresden, Theoretical Chemistry, Bergstr. 66c, 01062 Dresden, Germany
Helmholtz-Zentrum Dresden-Rossendorf, Abteilung
Ressourcenökologie,Forschungsstelle Leipzig, Leipzig, Germany  Yonsei
University and ibs-cnm, Seodaemun-gu, Seoul 120-749, Republic of Korea

###### Abstract

Manipulating the interlayer twist angle is a powerful tool to tailor the
properties of layered two-dimensional crystals. The twist angle has a
determinant impact on these systems’ atomistic structure and electronic
properties. This includes the corrugation of individual layers, formation of
stacking domains and other structural elements, and electronic structure
changes due to the atomic reconstruction and superlattice effects. However,
how these properties change with the twist angle, $\theta$, is not yet well
understood. Here, we monitor the change of twisted bilayer MoS2
characteristics as a function of $\theta$. We identify distinct structural
regimes, each with particular structural and electronic properties. We employ
a hierarchical approach ranging from a reactive force field through the
density-functional-based tight-binding approach and density-functional theory.
To obtain a comprehensive overview, we analyzed a large number of twisted
bilayers with twist angles in the range of $\theta=0.2^{\circ}\dots
59.6^{\circ}$. Some systems include up to half a million atoms, making
structure optimization and electronic property calculation challenging. For
$13^{\circ}\lessapprox\theta\lessapprox 47^{\circ}$, the structure is well-
described by a moiré regime composed of two rigidly twisted monolayers. At
small twist angles ($\theta\leq 3^{\circ}$ and $57^{\circ}\leq\theta$), a
domain-soliton regime evolves, where the structure contains large triangular
stacking domains, separated by a network of strain solitons and short-ranged
high-energy nodes. The corrugation of the layers and the emerging superlattice
of solitons and stacking domains affects the electronic structure. Emerging
predominant characteristic features are Dirac cones at $K$ and kagome bands.
These features flatten for $\theta$ approaching $0^{\circ}$ and $60^{\circ}$.
Our results show at which range of $\theta$ the characteristic features of the
reconstruction, namely extended stacking domains, the soliton network, and
superlattice, emerge and give rise to exciting electronics. We expect our
findings also to be relevant for other twisted bilayer systems.

##  I Introduction

The idea of combining two-dimensional (2D) crystals by stacking opened up an
intense research focus in the emergent field of van der Waals (vdW)
[hetero]structures.[1] The ability to control the stacking order and to
introduce an interlayer twist angle, $\theta$, gives seemingly endless
possibilities to tune the properties of these materials.[2, 3, 4, 5, 6] A
consequence of this is the formation of moiré patterns, [7, 8, 9, 10] which
exhibit a continuous change of the local stacking as visualized in Fig. 1(a).
In particular for small twist angles near $0^{\circ}$ and $60^{\circ}$,
domains of low-energy high-symmetry stackings are formed. They act as seeds
for a relaxation process, where energy minimization yields a maximization of
the size of low-energy domains, as visualized in Fig. 1(b). Relaxation affects
the in-plane (ip) structure, forming the stacking domains as dominating
structural elements with local high-symmetry stacking (see Fig. 1(c)). By
corrugation, also the out-of-plane (oop) structure is affected. So-called
strain solitons [11, 12, 13, 14], with the nodes at their intersections,
dominate the corrugated areas. Similar effects have been observed in
simulations of a range of different twisted bilayer (tBL) materials, such as
graphene [3, 12], MoS2 [2, 15], or MoS2/WSe2 [16]. Several recent works also
give strong experimental evidence of such reconstruction effects.[10, 17, 18]

Figure 1: Atomic structure of a twisted bilayer (tBL) MoS2 system with a twist
angle $\theta=1.05^{\circ}$ in the case of (a) rigidly twisted and (b) fully
relaxed layers (Mo - blue, S - red). Areas of local high-symmetry stacking are
indicated. (c) Structural elements and their respective local stacking for
$\theta\rightarrow 0^{\circ}$. The stacking domains (green triangles) exhibit
local $R_{h}^{M}$ stacking, either seen from the top (S-rich side) or the
bottom (Mo-rich side). The solitons (blue rectangles) show a transition
stacking type. The nodes (red circles) have local $R_{h}^{h}$ stacking, with
eclipsed S atoms. The local stackings for $\theta\rightarrow 60^{\circ}$ are
shown in the SI (Fig. S1).

Depending on the value of $\theta$ ($0^{\circ}\leq\theta\leq 60^{\circ}$), the
characteristics of the systems change, and different structural regimes are
obtained: A $\theta$ of exactly $30^{\circ}$ results, as in bilayer graphene,
in a quasicrystalline arrangement.[19] For large twist angles
($\theta\rightarrow 30^{\circ}$), in the so-called moiré regime [10],
structural relaxation of individual layers is minor.[15, 12] Here, the weak
oop deformation exhibits a sinusoidal modulation and decreases further for
$\theta\rightarrow 30^{\circ}$.[3, 20, 2] At small twist angles,
$\theta\rightarrow 0^{\circ}$ and $\theta\rightarrow 60^{\circ}$, the
structure is commonly described by a (strain) soliton regime, where an
extended atomic reconstruction both ip and oop takes place.[10, 21, 22, 23]
This results in a regular pattern in the interlayer distance landscape which
can be understood in terms of structural elements.[24] The first structural
element is the so-called stacking domain – large areas of minimum interlayer
distance and local low-energy stacking, marked with green triangles in Fig.
1(b) for $\theta\rightarrow 0^{\circ}$ (see Fig. S1 for comparison to
$\theta\rightarrow 60^{\circ}$). They are separated by domain walls of
transition stacking, usually called solitons or strain solitons, as a second
structural element.[11, 12, 13, 14] The solitons are marked by blue rectangles
in Fig. 1(b) and act as a network of domain walls, separating adjacent
stacking domains and vanishing only at exactly $0^{\circ}$ or $60^{\circ}$.
The stacking domains form a superlattice with honeycomb (hcb) symmetry, while
the centers of the solitons form a kagome (kgm) lattice.[25, 26] As the final
structural element, the so-called nodes are formed at the intersections of the
solitons, marked with red circles in Fig. 1(b). These are areas with maximum
interlayer distance and local high-energy stacking. The separation of the
nodes correlates with the moiré periodicity. The twist angles at which the
transition between structural regimes take place depend on which system is
studied, e.g., $\approx 1^{\circ}$ in tBL graphene [12, 10, 20], $\approx
4^{\circ}$ in tBL WSe2 [27], and up to $13^{\circ}$ in tBL MoS2 [15]. However,
the criteria to define these transition angles are not uniquely defined in the
literature. Different reports use deformation energy, moiré periodicity, or
the size of structural elements, amongst others.

Introducing a small interlayer twist into a homobilayer system results in
several hundreds or thousands of atoms in the moiré supercell. While at
$\theta=2^{\circ}$, there are 5000 atoms, this number increases to 20,000 at
$1^{\circ}$ and reaches more than 100,000 for twist angles below
$0.44^{\circ}$. This makes simulations challenging and calls for approximate
methods. Thus, to date, the literature is restricted to results from continuum
mechanics [11] and force field [28, 22] approaches.

In this work, we elucidate how the structure and electronic properties (band
structure, local gap, and electronic transport) of tBL MoS2 change with the
twist angle, $\theta$. To provide a complete picture, we obtain our results
for a very large set of systems, allowing us to describe structural changes
with $\theta$ almost continuously. We derive a complete characterization of
the structural regimes in tBL MoS2 and the transition between these regimes.
We analyze the strong atomic reconstruction for very small $\theta$ and the
emergence of Dirac cones and kagome bands in the electronic structure
associated with developing superlattices of the structural elements. These
features can only be visible for fully relaxed structures and are missing in
band structures of rigidly twisted layers. This calls for atomic relaxation,
especially for very small $\theta$. Our approach can be adapted to other
twisted vdW structures where atomic relaxation effects play a role.

##  II Methods

Figure 2: High-symmetry stackings of MoS2 bilayers in the top and side view:
(a) $R_{h}^{M}$ (corresponds to the 3$R$ stacking), (b) $R_{h}^{h}$, (c)
$H_{h}^{h}$ (corresponds to the 2$H$ stacking), (d) $H_{h}^{X}$, and (e)
$H_{h}^{M}$. The $R$-type corresponds to a $0^{\circ}$ rotation, the $H$-type
to a $60^{\circ}$ rotation between the layers. [48] The stacking type and the
interlayer distance of the fully relaxed structures are given together with
the total energy differences, $\Delta E$, relative to the lowest-energy
stacking $H_{h}^{h}$. Energies are given in meV per unit cell (uc).

Initial tBL MoS2 structures with periodic boundary conditions were constructed
from ideal monolayer (ML) equilibrium geometries placed on top of one another
with twist angles $\theta$ between $0^{\circ}$ and $60^{\circ}$. Details on
the structure generation are given in the Supporting Information (SI) [29]
(see Fig. S2). Full geometry optimization was performed employing the LAMMPS
code [30] using a reactive force field (ReaxFF) [31, 32, 33] parametrized for
MoS2.[28] Geometry convergence was considered to be reached at a maximum force
component of 0.1 eV/Å and the remaining pressure of the atomic system below 10
kbar. We obtained 1219 fully optimized tBL structures with $\theta$ in the
range of $0.2^{\circ}-59.6^{\circ}$. All fully-relaxed tBL MoS2 structures
used in this work and the calculation files are available as ZENODO
repository.[34] The results of the force field calculations were verified and
are in good agreement with density-functional theory (DFT), as described in
the SI (see Tables S1, S2, and Fig. S3). There is, however, a consistent
overestimation of the interlayer distance in (t)BLs by about 0.35 Å. The moiré
patterns introduced by rigidly twisted layers change with $\theta$, and the
resulting size of the periodic unit strongly increases for small twist angles
($\theta\rightarrow 0^{\circ}$ and $\theta\rightarrow 60^{\circ}$). While at
$\theta=1.05^{\circ}$, there are already 17,862 atoms in the unit cell, this
value reaches 493,026 atoms in the largest studied structure with
$\theta=0.2^{\circ}$.

The primary descriptor for analysis of the structural properties is the local
interlayer distance $d(x,y)$, defined as the vertical distance between the
planes spanned by the Mo atoms of the two layers. As the atoms in the two
layers are not necessarily eclipsed, multivariate interpolation was applied to
map the atomic coordinates to a quasi-continuous description of $d(x,y)$. The
size of the structural elements was determined using the so-called moiré cell
– the smallest repeating unit in the local interlayer distance landscape.
Additional information on the analysis of structural and energetic properties
is given in the SI (see Figs. S4 and S5). A movie file showing how the
interlayer distance landscape changes with $\theta$ is available as part of
the ZENODO repository.[34]

Figure 3: Atomistic structure of tBL MoS2 systems with representative $\theta$
corresponding to the different structural regimes: moiré regime, transition
regime, soliton regime, and domain-soliton regime. (a) Top views of the fully
relaxed atomistic structures with Mo atoms in blue and S atoms in red. (b)
Interlayer separation landscape, calculated from the positions of the Mo
atoms.

Electronic structure calculations were carried out using density functional
theory (DFT) as implemented in FHI-aims.[35, 36, 37] The Perdew-
Burke–Ernzerhof (PBE) exchange-correlation functional [38] was used together
with a light tier 1 numeric atomic orbital basis set [39], a converged
Monkhorst-Pack $k$-point grid [40] with a density of 5
$\nicefrac{{1}}{{\text{\AA}^{-1}}}$ along the reciprocal lattice vectors, and
non-self consistent spin-orbit coupling (SOC).[41] Using bilayer and twisted
bilayer MoS2 systems, we verified that band structures calculated using
ReaxFF- and DFT-optimized structures are consistent (see Figs. S6, S7, S8, and
Table S3). Electronic properties were analyzed from band structure and density
of state (DOS). Atom-projected DOS was calculated with a Gaussian broadening
of 50 meV. From this, the local gap of the Mo atoms was approximated using the
first maximum below the Fermi level as the valence band maximum (VBM) and the
onset above the Fermi level as the conduction band minimum (CBM). To further
show which parts of the structure contribute to the VBM, the atom-projected
DOS of the Mo atoms at the VBM energy of the whole structure was evaluated.

For exemplary systems, the electronic band structure and DOS were calculated
using the lowest-order density-functional based tight-binding (DFTB) method,
DFTB0, [42] as implemented in QuantumATK,[43] employing the parameterization
of Wahiduzzaman et al. [44] The $k$-grid was set to $2\times 2\times 1$ for
large supercells with $\theta\leq 7^{\circ}$ and to $5\times 5\times 1$ for
$\theta>7^{\circ}$. The electronic structure calculated using DFTB was
validated against DFT and found to be in fairly good agreement (see Figs. S9,
S10, and S11). We observe, however, a consistently larger band gap by $\approx
0.4$ eV.

Ballistic electron transport calculations within DFTB were performed,
employing the Landauer-Büttiker formalism and the non-equilibrium Green’s
function (NEGF) approach.[45, 46] The device simulated with NEGF consists of
three parts: the left and right (semi-periodic) leads attached to the finite
central scattering region (see Fig. S12 for an exemplary device
configuration). A dense $k$-point sampling of $1\times 5\times 99$ was used.
The whole device is periodic normal to the transport direction in the xy
plane. A vacuum of 20 Å normal to the bilayer plane was used to avoid any
interaction caused by periodic images.

Details for all methods used in this work are given in the SI.

##  III Results and Discussion

###  III.1 Structural and energetic properties

Figure 4: Structural, energetic, and electronic properties of tBL MoS2 as
function of $\theta$. (a) Maximum, minimum, and mean interlayer distance, $d$,
calculated from the Mo atoms in the neighbouring layers. The reference lines
correspond to the interlayer distance in the high-symmetry BL stackings. (b)
Cross-section through the interlayer distance landscape for
$\theta=1.05^{\circ}$. A Gaussian is fitted to the node (red) and soliton
(blue) to approximate their width, $w$. (c) Numerically determined height,
$h$, of solitons (blue) and nodes (red), obtained from the maximum interlayer
distance inside the structural element relative to the minimum value of the
structure. (d) Width, $w$, of solitons (blue) and nodes (red), obtained from
the base width $4\sigma$ of the Gaussian fitted to the cross-section as shown
in (b). The numerically determined moiré periodicity, $a_{\mathrm{M}}$, is
shown in grey. (e) Exfoliation energy of tBL MoS2. The $R_{h}^{M}$ and
$H_{h}^{h}$ stackings with the highest value of 0.30 J/m² are indicated with a
dashed blue line. (f) Contributions to the exfoliation energy originating from
the deformation of the layers and change of the interlayer adhesion due to the
corrugation. (g) Band gaps for rigidly twisted and fully relaxed structures.
As a reference, the band gap of the relaxed BL systems is shown. The relaxed
ML has a band gap of 1.65 eV. (h) Absolute energy of VBM and CBM when
including vacuum corrections. In all subplots, the structural regimes are
indicated with the (domain-)soliton regime in (dark) light red, the transition
regime in green, and the moiré regime in yellow.

The primary reference systems needed for the analysis of tBL MoS2 are the
high-symmetry stacking configurations in untwisted bilayer (BL) MoS2, as shown
in Fig. 2. The calculated structural properties of the BL MoS2 stackings are
consistent with previous reports in the literature [47] with more details
given in the SI (see Table S4). In contrast to heterobilayer systems, there
are only five possible stackings because $R_{h}^{M}$ and $R_{h}^{X}$ are
equivalent in homobilayer systems. $R_{h}^{M}$ and $H_{h}^{h}$ belong to the
S@Mo (read ”S on Mo”) stacking type – low-energy stackings with the smallest
interlayer distance of 6.03 Å, where the S atoms in the top layer eclipse the
Mo atoms in the bottom layer. $R_{h}^{h}$ and $H_{h}^{X}$ belong to the S@S
stacking type – high-energy stackings which are less stable by 54 meV/uc and
have a large interlayer distance of 6.50 Å. $H_{h}^{X}$ belongs to the Mo@Mo
stacking type with energy and interlayer distance comparable to that of
$H_{h}^{h}$. With this, three high-symmetry stacking configurations, namely
$R_{h}^{M}$, $H_{h}^{h}$, and $H_{h}^{X}$, with low energy and short
interlayer distances exist. Their relative energies are within 2 meV per
primitive BL unit cell (uc), consisting of 2 MoS2 formula units. This results
in the stacking domains for $\theta\rightarrow 0^{\circ}$ showing local
$R_{h}^{M}$ stacking, while for $\theta\rightarrow 60^{\circ}$ they alternate
between $H_{h}^{h}$ and $H_{h}^{X}$. This makes the $R$-type stacking domains
energetically equivalent, while the $H$-type stacking domains have different
stacking energy, resulting in a size difference between adjacent stacking
domains for twist angles near $60^{\circ}$ (see Fig. S1). [26] Local
$R_{h}^{h}$ ($H_{h}^{X}$) stacking arrangements are present within the nodes
for $\theta\rightarrow 0^{\circ}$ ($\theta\rightarrow 60^{\circ}$).

Figure 5: Selected structural properties of exemplary tBL MoS2: (a) Local
interlayer distance $d$. (b) Isotropic strain field, $\epsilon$, relative to
the unstrained MoS2 ML. Areas of compressive strain are shown in red, tensile
strain in blue. (c) Absolute value of the in-plane component of the
anisotropic strain field, $|\vec{u}_{\mathrm{ip}}|$, with the qualitative
directional component of the vector field indicated for $\theta=0.915^{\circ}$
using arrows. (d) Out-of-plane component of the anisotropic strain field,
$|\vec{u}_{\mathrm{oop}}|$. The properties in (b-d) are plotted only for the
upper layer and interpolated between the data points.

Structural and electronic properties of tBL MoS2 gradually change as a
function of $\theta$. As most properties of tBL MoS2 show a continuous change
with $\theta$, there is a continuous transition between the soliton and moiré
regimes rather than a single ”magic” twist angle at which this transition
occurs. Nevertheless, for certain $\theta$ particular structural transitions
emerge. Here, we provide a comprehensive characterization of structural
properties as a function of $\theta$.

Exemplary structures for each regime are given in Fig. 3 and details on the
criteria used for the in-depth description of the regimes are given in the SI.
The moiré regime is present for large $\theta$ and can be modeled with rigidly
twisted layers. When going towards smaller $\theta$, the characteristics of
the system change as the influence of relaxation becomes more prominent. The
structural elements start to appear, indicating the beginning transition
towards the soliton regime. These intermediate twist angles can be classified
by a transition regime. This accounts for the continuous transition towards
the (strain) soliton regime.[10, 21, 22] In this work, it is divided into two
parts due to structural and electronic property differences: For small twist
angles, we classify the soliton regime. It exhibits all characteristic
structural elements: The stacking domains form a honeycomb lattice but are
still relatively small. The solitons are present as bridges between the nearby
nodes. When going towards the smallest $\theta$, a more extensive atomic
reconstruction is observed: Extended stacking domains dominate the structure.
The solitons form a hexagonal network with the nodes positioned at the
intersection points. To account for these structural differences, We classify
the smallest twist angles as the domain-soliton regime.

The structural and energetic properties studied as a function of $\theta$ are:
(i) The interlayer distance $d$ with its minimum, $d_{\mathrm{min}}$, maximum,
$d_{\mathrm{max}}$, and mean, $d_{\mathrm{mean}}$, values as shown in Fig.
4(a). Here, gaps in the curves around $\theta=21.79^{\circ}$, $38.21^{\circ}$,
and $13.17^{\circ}$ are visible. They originate from ”magic angles” with
respect to the size of the moiré cell. A small cell can be realized at these
values while diverging to unfeasibly large structures at nearby $\theta$. (ii)
The spatial extension of nodes and solitons, analyzed by fitting a Gaussian
function to their cross-section in the interlayer distance landscape as
demonstrated in Fig. 4(b). Their height, $h$, as the maximum corrugation of
the interlayer distance inside a structural element, is visualized in Fig.
4(c). Their width, $w$, is shown in Fig. 4(d), using the diameter of the
circular nodes and the width at the center of the solitons to approximate
their spatial extension. This property is also compared to the moiré
periodicity (see Figs. S13 and S14 for more details) as the periodic unit of
the moiré pattern and used to estimate the area of the stacking domains (see
Fig. S15). (iii) The exfoliation energy is given in Fig. 4(e). It is further
separated into individual contributions from the deformation of the individual
layers and the change of the interlayer adhesion (see Fig. 4(f)) due to the
relaxation of the layers (for more details see Fig. S16). The exfoliation
energy as a function of $\theta$ can be described by piecewise linear parts
(see Fig. S17 and Table S5).

Additionally, the displacement and strain fields were analyzed. The
displacement field (see Fig. S18) allows to gain insight into the atomic
reconstruction during relaxation.[49] The intralayer strain distribution was
investigated using the isotropic and anisotropic strain fields, as shown in
Fig. 5, both giving access to different structural information.[3] While the
scalar isotropic strain field, $\epsilon$, gives the local strain relative to
an unstrained ML, the anisotropic strain field, $\vec{u}$, shows the
directional change of the strain as a vector field. Equations and additional
information are also given in the SI (see Fig. S19).

Inside the moiré regime, between $\approx 13^{\circ}$ and $\approx
47^{\circ}$, an almost constant value of $d_{\mathrm{mean}}$ with only minor
interlayer distance modulation is found (see Fig. 4(a)). This indicates that
the rigidly twisted model can be applied in this regime, and the structural
elements are not pronounced. Similarly to the interlayer distance, the
exfoliation energy is almost constant with a minor change of $\approx 0.15$
mJ/m² per $1^{\circ}$. The anisotropic strain field remains non-zero, as the
anisotropy becomes zero only for untwisted bilayers at $0^{\circ}$ and
$60^{\circ}$ (see Fig. S19). The displacement field being close to zero
confirms that the ip reconstruction is negligible.

In the transition regime, $6^{\circ}<\theta<13^{\circ}$ (equivalently,
$47^{\circ}<\theta<54^{\circ}$), the atomic reconstruction becomes more
prominent. Here, the interlayer distance shows a slow decrease of
$d_{\mathrm{mean}}$ with decreasing $\theta$. This is reflected in the
exfoliation energy, which slowly increases by $\approx 1.4$ mJ/m² per
$1^{\circ}$. A strongly increasing interval $\Delta
d=d_{\mathrm{max}}-d_{\mathrm{min}}$ (from 0.16 Å at $13^{\circ}$ to 0.41 Å at
$6^{\circ}$) is found, showing an intensifying oop deformation of the
individual layers (see Fig. 4(a)). Contrary to this, the ip displacement field
still shows a negligible magnitude (see Fig. S18(b)) with mean values below
0.03 Å. This confirms oop relaxation as the dominating effect of the atomic
reconstruction in this regime, associated with emerging structural elements.
An increasing spatial extension of nodes and solitons with decreasing $\theta$
is observed as a strong increase in their height and width, as shown in Fig.
4(c, d). While at $13^{\circ}$, the base width of the nodes (solitons) is 13 Å
(8 Å), it increases to 30 Å (15 Å) at $6^{\circ}$. The nodes are the
dominating structural element of the transition regime as they appear already
at larger values of $\theta$ near $13^{\circ}$ ($47^{\circ}$). Solitons appear
only for lower $\theta$ near $6^{\circ}$ ($54^{\circ}$), indicating the
beginning of the soliton regimes. At the same time, the anisotropy of the
strain field is becoming large and shows an emerging order with decreasing
$\theta$.

In the soliton regime, $3^{\circ}<\theta\leq 6^{\circ}$
($54^{\circ}\leq\theta<57^{\circ}$), the atomic reconstruction continues to
become more prominent, and all three structural elements are present. The
$d_{\mathrm{max}}$ value reaches its upper bound, equivalent to $R_{h}^{M}$
($H_{h}^{h}$) stacking (S@S type), at $\theta=6^{\circ}$
($\theta=54^{\circ}$). This correlates with fully pronounced nodes of local
high-energy stacking type. The height of the structural elements increases and
reaches its maximum value at $\theta=3^{\circ}$ ($\theta=57^{\circ}$).
Similarly, the width of nodes and solitons also increases. For $\theta$ of
$\approx 3^{\circ}$ ($\approx 57^{\circ}$), the width of the nodes coincides
with $a_{\mathrm{M}}$, showing that the nodes are the dominating structural
element of the moiré cell. Their large size also results in the solitons being
present only as short bridges between adjacent nodes, as shown in Fig. 3(b).
The stacking domains are still small but already visible in the interlayer
distance landscape, forming a honeycomb lattice. They exhibit their
characteristic local S@Mo stacking, but the interlayer distance of the S@Mo BL
is not yet reached, as visible from $d_{\mathrm{min}}$ decreasing further.
Their small size is also visible from $d_{\mathrm{mean}}$ and the exfoliation
energy, which show a small change with decreasing $\theta$. The dominating
factor in the atomic reconstruction is still the oop relaxation. The
displacement field becomes more sizeable in this regime with values up to 0.08
Å (see Fig. S18), showing the beginning influence of the ip reconstruction.

The final regime is the domain-soliton regime, $0^{\circ}<\theta\leq
3^{\circ}$ ($57^{\circ}\leq\theta<60^{\circ}$). Additionally to the oop
relaxation, a significant ip reconstruction is present in this regime. This
can be deduced from the displacement field, which shows strongly increasing
values of up to $\approx$0.8 Å (see Fig. S18(b)). Its largest values are at
the boundary between the structural elements. The center of nodes and stacking
domains form high-symmetry stackings in the rigidly twisted structure and,
thus, can be understood as seeds for the atomic reconstruction.

In the domain-soliton regime, $d_{\mathrm{min}}$ reaches its lower bound,
equivalent to $R_{h}^{h}$ and $H_{h}^{M}$ stackings, showing that the stacking
domains are fully pronounced. As $\Delta d$ becomes maximal, the height of
nodes and solitons reaches a plateau, with the nodes being at least three
times larger than the solitons (see Fig. 4(c)). A plateau is also observed for
their width at $\theta\leq 1^{\circ}$ ($\theta\geq 59^{\circ}$): the solitons
reach a maximum of 57 Å, and the nodes a value of 77 Å (see Fig. 4(d)). In the
literature, the appearance of such a plateau is also described.[8] Naik et
al.[50] find the full width at half maximum (FWHM) of the solitons to be
$\approx$15 Å at $3.48^{\circ}$. This is in excellent agreement with our
results, giving $w=$22.7 Å for $3.48^{\circ}$, corresponding to a FWHM of 15.7
Å. As the spatial extension of the nodes becomes maximal, their width no
longer coincides with the moiré periodicity $a_{\mathrm{M}}$, as shown in Fig.
4(d). Decreasing $\theta$ increases the separation between the nodes, creating
a pronounced network of solitons and extended stacking domains in the domain-
soliton regime (see Fig. 3(b)), giving the regime its name. The stacking
domains become the dominating structural element (see Fig. S15). At
$3^{\circ}$, only 1.5% of the moiré cell can be accounted for the stacking
domains, compared to up to 80% at $0.2^{\circ}$. The formation of extended
stacking domains is reflected in $d_{\mathrm{mean}}$ strongly decreasing
towards the value of the S@Mo stacking type (see Fig. 4(a)). It is also
visible as a strong increase of the exfoliation energy by $\approx 21$ mJ/m²
per $1^{\circ}$ (see Fig. S17 and Table S5). Strain effects account for only
$\approx$0.1% of the exfoliation energy, the deformation of individual layers
$\approx$1%, making the interlayer adhesion terms the dominating contributions
(see Fig. S16). The values of the exfoliation energy (0.29 J/m² for small
$\theta$ and 0.24 J/m² for large $\theta$) are in good agreement with results
by Emrem et al.,[51] who reported 0.33 J/m² for BL MoS2, calculated at RPA
level. The structural properties described for the domain-soliton regime are
reflected in the strain fields, as shown for examples with $\theta\rightarrow
0^{\circ}$ and $\theta\rightarrow 60^{\circ}$ in Fig. 5. The stacking domains
exhibit a constant isotropic strain $\epsilon$ with negligible anisotropy both
ip and oop, confirming that these are large, uniform S@Mo areas. The nodes and
boundary region between solitons and stacking domains, on the other hand, show
strong anisotropy. As indicated in Fig. 5(c), the directional component has a
circular arrangement around the nodes and, as reported by Alden et al., shows
the presence of sheer strain at the solitons.[52] Due to this strain
accumulation, the solitons are often called strain solitons. The oop
component, qualitatively comparable to the curvature of the individual layers,
points in opposite directions in the individual layers. It is largest in nodes
and solitons, as the main curvature of the layers is localized in these
structural elements.

It is important to note that the interlayer separation landscape in structures
with $\theta\rightarrow 0^{\circ}$ and $\theta\rightarrow 60^{\circ}$ deviates
strongly for $\theta<1^{\circ}$ and $\theta>59^{\circ}$ (see Fig. 5(a)). For
$\theta\rightarrow 0^{\circ}$, the stacking domains are equivalent, and the
solitons are linear. This changes for $\theta\rightarrow 60^{\circ}$, where
small and large stacking domains alternate, separated by bent solitons.
Different symmetries cause this: in $\theta\rightarrow 0^{\circ}$, all
stacking domains exhibit a local $R_{h}^{M}$ stacking, while in
$\theta\rightarrow 60^{\circ}$, the adjacent domains alternate between
$H_{h}^{h}$ and $H_{h}^{X}$ stackings, as shown in Fig. S1.[8, 15, 26] The
$H_{h}^{h}$ and $H_{h}^{X}$ stackings differ in energy (2 meV/uc, see Fig. 2).
Thus, the size of stacking domains also differs, with the regions of more
stable stacking becoming larger. This results in the formation of bent
solitons and a lowering of the overall symmetry. It is accompanied by
structures with $\theta\rightarrow 60^{\circ}$ being stronger bound by up to
1.6 mJ/m² compared to $\theta\rightarrow 0^{\circ}$ (see Fig. S20), resulting
from a stronger adhesion between the layers. We observe that even though there
are such substantial differences between $\theta\rightarrow 0^{\circ}$ and
$\theta\rightarrow 60^{\circ}$, there is an almost ideal symmetry
$\theta_{i}=60^{\circ}-\theta^{\prime}_{i}$ for the twist angles at which the
transition between the regimes takes place. It is interesting to note that
introducing an inherent interlayer strain also results in a bending of the
solitons, as shown in the SI (see Fig. S21).

Our characterization approach is fully consistent with the literature results
and shares common criteria other authors used to asses the transition between
the regimes. Garguilo et al.[12] analyze the twist energy – a property
comparable to the deformation energy used here – of tBL graphene as a function
of $a_{\mathrm{M}}$. They find $1.2^{\circ}$ as the crossover angle between a
”solitonic regime” and a ”rigid regime.” Dai et al.[20] describe the
transition between the regimes to happen at the coincidence point between
$a_{\mathrm{M}}$ and the width of the solitons. They find values from
$1^{\circ}$ to $1.6^{\circ}$ for tBL graphene, similar to the value of
$3^{\circ}$ in the more rigid MoS2 systems.[53] Maity et al.[15] find a
maximum deformation of the individual layers in tBL MoS2 at $3^{\circ}$ and a
transition angle of $13^{\circ}$. This is in direct agreement with the begin
of the domain-soliton regime and the moiré regime as defined in this work.
Wang et al.[27] describe the transition between the structural regimes to take
place in an interval from $4^{\circ}$ to $5.1^{\circ}$ for tBL WSe2. This
agrees with the soliton regime, where the structural elements start to vanish.
Enaldiev et al.[54] find values of $2.5^{\circ}$ and $59^{\circ}$ in tBL TMDC
systems as characteristic twist angles for forming extended stacking domains.
This is the only other work known by the authors to explicitly consider the
range $30^{\circ}\leq\theta<60^{\circ}$ and agrees with the newly defined
domain-soliton regime. Overall, we find an excellent agreement of the
characterization derived in this work with previously published results. Our
analyses are based on a much larger set of twisted bilayer structures,
allowing us to find a smooth transition between the regimes. The approach
proposed in this work is expected to be transferable to different twisted
homobilayers and, after small adaptations, to other systems, including
heterobilayers.

###  III.2 Electronic properties

We now investigate the impact of structural relaxation on the global and local
electronic properties of tBL MoS2 and how they change in the different
regimes. As a starting point, the electron density at the VBM is shown in Fig.
6 for examples of the different regimes. We observe that it follows
qualitatively the behavior of the interlayer distance landscape. Only small
modulations are found within the moiré regime, which increase in the
transition regime. In the soliton regime, the nodes become pronounced as areas
of minimum electron density. In contrast, extended stacking domains with
uniform maximum electron density and a network of solitons are found in the
domain-soliton regime.

Figure 6: (a) Electron density at the VBM, calculated at DFTB level, for
representative structures in the structural regimes, corresponding to Fig.3.
(b, c) Atom-projected band structures of the (b) fully relaxed and (c) rigidly
twisted systems, including the $R_{h}^{M}$ bilayer as a reference.

Changes between the regimes are also found in the atom-projected band
structures, as shown in Fig. 6(b, c). All tBL MoS2 systems are indirect band
gap semiconductors, with the VBM at $\Gamma$, similar to the untwisted bilayer
limits. The CBM is positioned at $K$, which changes to be along the path
between $M$ and $\Gamma$ with decreasing $\theta$ in the transition regime. An
exception are those systems with a direct gap due to backfolding (see SI,
Table S3; unfolding these band structures is beyond this work’s scope). The
atomic projection reveals that the frontier bands have a strongly mixed
character (contributions from both the Mo and S atoms). These bands show the
most substantial changes between rigid and relaxed models. Occupied bands with
an increased Mo character are away from the frontier bands. They show only
minor relaxation effects but are more sensitive against the method used for
geometry optimization.

The shape of bands and the band gap values change with $\theta$ (see Fig.
4(g)). The band gap difference between fully relaxed and rigidly twisted
systems is very small, with a maximum of 44 meV, showing that there is just a
minor influence of relaxation effects on this property. At large $\theta$ in
the moiré regime, the band structure differences between relaxed (Fig. 6(b))
and rigidly twisted (Fig. 6(c)) structures are also minor. The band gap does
not change with $\theta$ around $30^{\circ}$. This further supports the
validity of the model of rigidly twisted layers in this regime. With
decreasing $\theta$, the band gap slowly decreases and approaches the value of
the S@Mo BL system, showing a linear dependence on $d_{\mathrm{mean}}$ (see SI
for details). This change originates from an increase of the VBM and a
decrease of the CBM (see Fig. 4(h)), with the former being the dominating
factor. Furthermore, both valence and conduction bands show a decreasing
dispersion. This trend continues in the soliton regime.

Fig. 7 shows an exemplary band structure of tBL MoS2 with
$\theta=3.89^{\circ}$, calculated for relaxed and rigidly twisted cases. The
band gap values of both systems are about 1.25 eV. However, the electronic
structures differ significantly.

Figure 7: Electronic band structure, the corresponding DOS and conductance for
both the rigidly twisted (blue) and fully relaxed (red) structure for tBL MoS2
with $\theta=3.89^{\circ}$, computed at the DFTB level. The device
configuration consists of left and right semi-infinite leads and a finite
central region of three unit cells (see Fig. S12) to ensure a smooth
transition between the potential of the central region and the leads.

The top VB becomes flat in the rigidly twisted system, giving rise to van Hove
singularity in the DOS. Below that flat band, there is a gap of about 60 meV.
In the fully relaxed system, the corrugation of the layers results in the
formation of a Dirac cone at $K$ close to the Fermi level. [26] This can also
be seen as a parabola-like shape of the DOS, similar to graphene. Below the
Dirac bands, characteristic kagome bands are observed. Compared to the rigidly
twisted system, the stronger band dispersion flattens the DOS at the VBM and
results in a higher conductance in the relaxed system.

The emergence of Dirac bands in the fully relaxed system is related to the
formation of a honeycomb lattice of strongly reconstructed stacking domains in
the (domain-)soliton regime (see Fig. 8(a)).

Figure 8: (a) Interlayer separation landscape with a schematic visualization
of (left overlay) the structural elements and (right overlay) the
superlattices formed by the structural elements in fully relaxed tBL MoS2 with
$\theta=2.28^{\circ}$. The stacking domains form a honeycomb lattice (orange)
with two sites per unit cell, while the solitons form a kagome lattice (blue)
with three sites per unit cell. (b) Electronic band structure of fully relaxed
tBL MoS2 with $\theta=2.28^{\circ}$, calculated at the DFTB level. The bands
corresponding to the honeycomb and kagome lattice are highlighted.

It is also described in several other publications.[55, 25, 56, 57, 26] It can
further be seen that the three VBs below the Dirac bands form distorted
kagome-type dispersion [58] (blue in Fig. 8(b)), originating from a kagome
lattice formed by the centers of the solitons (see Fig. 8(a)). The formation
of Dirac bands due to the honeycomb superlattice of the stacking domains can
be observed both in the soliton and the domain-soliton regime. It becomes more
pronounced as $\theta$ approaches $0^{\circ}$, as shown in Fig. S22(a). Their
bandwidth strongly decreases with decreasing $\theta$ (see Fig. S22(b)) and is
expected to result in flat bands for the smallest $\theta$. We discuss these
superlattice effects in more detail in another publication.[26]

These results show that the atomic reconstruction results in changes in the
electronic properties corresponding to lattice types not present in the
initial atomic structure. To better understand the origin of these changes and
their consequences in the soliton regimes, we analyzed the local electronic
properties of an exemplary structure. We used $\theta=3.89^{\circ}$ as the
smallest twist angle which is reliably accessible both in DFTB and DFT. The
properties studied are: (i) The local gap, a property similar to the transport
gap, as shown in Fig. 9(a). It can further be linked to the moiré potential as
it is affected mainly by the energy of the VBM in different parts of the
structure. Additional information is given in the SI (see Fig. S23 and Tables
S6 and S7). (ii) The LDOS at the VBM, as shown in Fig. 9(b), with additional
information on changes with $\theta$ given in the SI (see Table S8). Changes
in this property are connected to the variation of the local gap within the
structure, influencing the localization of charge carriers. (iii) The electron
density at the VBM, calculated only at the DFTB level, as shown in Fig. 9(c).
This property reveals the most probable region where charge carriers are
located. (iv) The transport pathways, as shown in Fig. 9(d). It illustrates
the most viable path along which charge carriers can be transported within the
device configuration at the energy of the VBM.

Figure 9: Local electronic properties in (left) rigidly twisted and (right)
fully relaxed structures: (a) local gap, (b) atom-projected DOS (LDOS) at the
VBM, (c) electron density, $\rho$, at the VBM, and (d) transmission pathways
in tBL MoS2 with $\theta=3.89^{\circ}$. The moiré cell is indicated. Plots (a)
and (b) are interpolated between the data points. The atomic structure is
added to subplot (a) to show the placement of structural elements. The results
shown in (a) and (b) are calculated at the DFT PBE level, while (c) and (d) at
the DFTB level.

Studying the rigidly twisted systems (left part of Fig. 9) demonstrates the
electronic structure obtained when structural relaxation is neglected. The
smallest value of the local gap is found in the S@S areas, while the largest
ones are found in the S@Mo areas and areas of intermediate stacking.
Consequently, a hexagonal pattern of areas with increased LDOS and electron
density at the VBM is formed. The main contributions are localized in the S@S
areas, corresponding to areas with the smallest local gap. Lower contributions
are found in the areas with high local gap. This has the consequence that
rigidly twisted layers would transport mainly through the S@S areas, resulting
in a rather low transmission through the device. The localized states at the
VBM are reflected in the flat VB, as shown in Fig. 7.

The observed patterns change completely when the system is allowed to relax.
The corrugated system shows the nodes as areas of the largest local gap, while
the stacking domains have the smallest and the solitons intermediate values.
The range of the local gap in tBL MoS2 with $\theta=3.89^{\circ}$ is reduced
compared to the high-symmetry BL stackings with deviations of $\approx$50 meV
from both $R_{h}^{M}$ and $R_{h}^{h}$ reference systems. Analogous to the
rigidly twisted system, areas of high (low) local gap correspond to a low
(high) LDOS and reduced (increased) electron density at the VBM. Due to this,
the stacking domains and partially the solitons show increased contributions
to the VBM, while the nodes do not contribute. The fully relaxed structure in
the soliton regime, thus, exhibits a honeycomb-shaped network of areas with
high DOS and electron density at the VBM, corresponding to a delocalization of
the holes. This is also reported elsewhere[50, 55, 59, 60, 61, 6] and
confirmed in experiments.[62] Towards the transition and moiré regime, the
modulation of the LDOS contributions vanishes. A pronounced transmission
pathway through the device can be observed. This is in good agreement with the
conductance calculations (see Fig. 7), which show an increased transport in
the relaxed compared to the rigidly twisted system. The most promising
direction for charge carriers to travel are the S@Mo areas and the solitons
with a constant, high transmission coefficient. The observed change in the
local electronic properties throughout the structure can further be used to
gain insights into the moiré potential. As this term is usually not well-
defined, the energy of the VBM (see Tables S6 and S7) can be used as an
approximation. The stacking domains have the highest moiré potential, while it
is lowered inside the nodes.

The observed honeycomb lattice of high LDOS at the VBM can be directly linked
to the observed effect of the highest two VBs forming graphene $p_{z}$-like
bands with a Dirac cone. The analysis of the LDOS contributions of the
different structural elements at the frontier bands (see Fig. S23(b)) supports
this. The two Dirac bands show majority contributions from the stacking
domains, forming a honeycomb lattice. The bands below show the majority
contributions originating from the solitons. This links the emergence of the
kagome-type bands to a kagome lattice formed by the center of the solitons
(see Fig. 8). With this, the clear link between atomic reconstruction, the
local electronic structure, and the emergent features in the band structure at
small $\theta$ is highlighted.

##  IV Conclusions

In this work, we studied the structural, energetic, and electronic properties
of tBL MoS2 by changing the twist angle $\theta$ quasi-continuously. Using a
large number of fully relaxed structures and many different properties, we
derive a new, comprehensive characterization of how the structure can be
described in different twist angle regimes. We further identified the $\theta$
values for which the transition between the regimes occurs. For small twist
angles $3^{\circ}<\theta\leq 6^{\circ}$ ($54^{\circ}\leq\theta<57^{\circ}$) in
the soliton regime, there is a strong influence of relaxation effects and an
extensive atomic reconstruction, leading to a corrugation of the individual
layers. For the smallest twist angles $\theta\leq 3^{\circ}$ ($\theta\geq
57^{\circ}$), we describe the domain-soliton regime, where large,
reconstructed stacking domains with a network of sharp domain boundaries
(solitons) are found. For large twist angles between $\approx 13^{\circ}$ and
$\approx 47^{\circ}$ inside the moiré regime, the model of rigidly twisted
layers is applicable. Between these regimes for $6^{\circ}<\theta<13^{\circ}$
($47^{\circ}<\theta<54^{\circ}$), the structure changes continuously, which is
described by the transition regime.

We also show how the comprehensive knowledge of the structural regimes and
their structural elements can be directly linked to the electronic properties
and their change with $\theta$. In rigidly twisted systems, flat bands form at
small $\theta$. When the system is allowed to relax, graphene $p_{z}$-like
bands with Dirac cones at $K$ are found. This originates from the emergence of
extended stacking domains, forming a honeycomb superlattice of regions with
low-energy stacking. We further observe that these bands become flat for
$\theta\rightarrow 0^{\circ}$ and that lower valence bands resemble those of a
kagome lattice originating from the solitons. These differences between the
rigidly twisted and the relaxed system are also visible in the ballistic
transport properties. Analysis of the local electronic properties using the
atom-projected DOS reveals a variation of the local gap of $\approx$100 meV
throughout the relaxed structure. In the corrugated system, this originates
from the states at the VBM being predominantly localized in the areas of
stacking domains and solitons. This also has consequences on the electron
density and the transmission pathways.

We demonstrate that including relaxation effects is crucial when calculating
the electronic structure of tBL MoS2 systems. Otherwise, the predicted
properties might be wrong not only quantitatively but also qualitatively. The
model of rigidly twisted layers is only valid inside the moiré regime. We
further show that force field-based calculations allow for a computationally
cheap approximation of the relaxed geometry, avoiding the computationally
expensive step of geometry optimization with, e.g., DFT. The force-field
relaxed structures are, thus, reliable to use for electronic structure
calculations in DFT. We expect the findings presented here to also be
important for other tBL systems beyond the studied MoS2 systems.

##  V acknowledgments

This project has been funded by Deutsche Forschungsgemeinschaft (DFG) within
the Priority Program SPP 2244 ”2DMP” and CRC 1415, and the European Union’s
Horizon 2020 research and innovation program under agreement No 956813.
Computational resources were provided by ZIH Dresden (project nano-10) and
Forschungszentrum Jülich (project opt2d). The authors would like to thank
Roman Kempt for his technical support with the calculation of electronic
properties and Louis Stuber for IT support at Forschungszentrum Jülich.

## References

  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Van Der Zande  _et al._ [2014] A. M. Van Der Zande, J. Kunstmann, A. Chernikov, D. A. Chenet, Y. You, X. Zhang, P. Y. Huang, T. C. Berkelbach, L. Wang, F. Zhang, _et al._ , Nano letters 14, 3869 (2014).
  * Van Wijk  _et al._ [2015] M. Van Wijk, A. Schuring, M. Katsnelson, and A. Fasolino, 2D Materials 2, 034010 (2015).
  * Lu  _et al._ [2017] N. Lu, H. Guo, Z. Zhuo, L. Wang, X. Wu, and X. C. Zeng, Nanoscale 9, 19131 (2017).
  * Dong  _et al._ [2021] R. Dong, A. Jacob, S. Bourdais, and S. Sanvito, npj 2D Materials and Applications 5, 1 (2021).
  * Vitale  _et al._ [2021] V. Vitale, K. Atalar, A. A. Mostofi, and J. Lischner, 2D Materials 8, 045010 (2021).
  * Constantinescu and Hine [2015] G. C. Constantinescu and N. D. Hine, Physical Review B 91, 195416 (2015).
  * Carr  _et al._ [2018] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras, Physical Review B 98, 224102 (2018).
  * Ribeiro-Palau  _et al._ [2018] R. Ribeiro-Palau, C. Zhang, K. Watanabe, T. Taniguchi, J. Hone, and C. R. Dean, Science 361, 690 (2018).
  * Yoo  _et al._ [2019] H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, P. Cazeaux, S. H. Sung, R. Hovden, A. W. Tsen, T. Taniguchi, _et al._ , Nature materials 18, 448 (2019).
  * Kumar  _et al._ [2016] H. Kumar, L. Dong, and V. B. Shenoy, Scientific reports 6, 1 (2016).
  * Gargiulo and Yazyev [2017] F. Gargiulo and O. V. Yazyev, 2D Materials 5, 015019 (2017).
  * Gadelha  _et al._ [2021] A. C. Gadelha, D. A. Ohlberg, C. Rabelo, E. G. Neto, T. L. Vasconcelos, J. L. Campos, J. S. Lemos, V. Ornelas, D. Miranda, R. Nadas, _et al._ , Nature 590, 405 (2021).
  * Lamparski  _et al._ [2020] M. Lamparski, B. Van Troeye, and V. Meunier, 2D Materials 7, 025050 (2020).
  * Maity  _et al._ [2020] I. Maity, M. H. Naik, P. K. Maiti, H. Krishnamurthy, and M. Jain, Physical Review Research 2, 013335 (2020).
  * Kunstmann  _et al._ [2018] J. Kunstmann, F. Mooshammer, P. Nagler, A. Chaves, F. Stein, N. Paradiso, G. Plechinger, C. Strunk, C. Schüller, G. Seifert, _et al._ , Nature Physics 14, 801 (2018).
  * 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, _et al._ , Nature Nanotechnology 15, 592 (2020).
  * Zhao  _et al._ [2023] S. Zhao, Z. Li, X. Huang, A. Rupp, J. Göser, I. A. Vovk, S. Y. Kruchinin, K. Watanabe, T. Taniguchi, I. Bilgin, _et al._ , Nature Nanotechnology , 1 (2023).
  * Ahn  _et al._ [2018] S. J. Ahn, P. Moon, T.-H. Kim, H.-W. Kim, H.-C. Shin, E. H. Kim, H. W. Cha, S.-J. Kahng, P. Kim, M. Koshino, _et al._ , Science 361, 782 (2018).
  * Dai  _et al._ [2016] S. Dai, Y. Xiang, and D. J. Srolovitz, Nano letters 16, 5923 (2016).
  * Kumar  _et al._ [2015] H. Kumar, D. Er, L. Dong, J. Li, and V. B. Shenoy, Scientific reports 5, 1 (2015).
  * Naik  _et al._ [2019] M. H. Naik, I. Maity, P. K. Maiti, and M. Jain, The Journal of Physical Chemistry C 123, 9770 (2019).
  * 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).
  * Zhu and Johnson [2018] S. Zhu and H. T. Johnson, Nanoscale 10, 20689 (2018).
  * Angeli and MacDonald [2021] M. Angeli and A. H. MacDonald, Proceedings of the National Academy of Sciences 118 (2021).
  * Arnold  _et al._ [2023] F. M. Arnold, A. Ghasemifard, A. Kuc, and T. Heine, ChemRxiv 10.26434/chemrxiv-2023-rx2fz-v2 (2023).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, _et al._ , Nature materials 19, 861 (2020).
  * Ostadhossein  _et al._ [2017] A. Ostadhossein, A. Rahnamoun, Y. Wang, P. Zhao, S. Zhang, V. H. Crespi, and A. C. Van Duin, The journal of physical chemistry letters 8, 631 (2017).
  * [29] See Supplemental Information (SI) at [URL will be inserted by publisher] for additional details on computational methods, a benchmark of ReaxFF structures using DFT, and information supporting the results shown in the main text.
  * Plimpton [1995] S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Van Duin  _et al._ [2001] A. C. Van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, The Journal of Physical Chemistry A 105, 9396 (2001).
  * Chenoweth  _et al._ [2008] K. Chenoweth, A. C. Van Duin, and W. A. Goddard, The Journal of Physical Chemistry A 112, 1040 (2008).
  * Aktulga  _et al._ [2012] H. M. Aktulga, J. C. Fogarty, S. A. Pandit, and A. Y. Grama, Parallel Computing 38, 245 (2012).
  * [34] The calculation files, the rigidly twisted and fully ReaxFF-optimized structures, and movies showing the change of the interlayer distance landscape and the strain fields with the twist angle are available at ZENODO. URL: doi.org/10.5281/zenodo.7243735.
  * Blum  _et al._ [2009a] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Computer Physics Communications 180, 2175 (2009a).
  * Havu  _et al._ [2009] V. Havu, V. Blum, P. Havu, and M. Scheffler, Journal of Computational Physics 228, 8367 (2009).
  * Marek  _et al._ [2014] A. Marek, V. Blum, R. Johanni, V. Havu, B. Lang, T. Auckenthaler, A. Heinecke, H.-J. Bungartz, and H. Lederer, Journal of Physics: Condensed Matter 26, 213201 (2014).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Blum  _et al._ [2009b] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Computer Physics Communications 180, 2175 (2009b).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Physical review B 13, 5188 (1976).
  * Huhn and Blum [2017] W. P. Huhn and V. Blum, Physical Review Materials 1, 033803 (2017).
  * Porezag  _et al._ [1995] D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, and R. Kaschner, Physical Review B 51, 12947 (1995).
  * Smidstrup  _et al._ [2020] S. Smidstrup, T. Markussen, P. Vancraeyveld, J. Wellendorff, J. Schneider, T. Gunst, B. Verstichel, D. Stradi, P. A. Khomyakov, U. G. Vej-Hansen, _et al._ , J. Phys: Condens. Matter 32, 015901 (2020).
  * Wahiduzzaman  _et al._ [2013] M. Wahiduzzaman, A. F. Oliveira, P. Philipsen, L. Zhechkov, E. Van Lenthe, H. A. Witek, and T. Heine, Journal of chemical theory and computation 9, 4006 (2013).
  * Stradi  _et al._ [2016] D. Stradi, U. Martinez, A. Blom, M. Brandbyge, and K. Stokbro, Physical Review B 93, 155302 (2016).
  * Büttiker  _et al._ [1985] M. Büttiker, Y. Imry, R. Landauer, and S. Pinhas, Phys. Rev. B 31, 6207 (1985).
  * He  _et al._ [2014] J. He, K. Hummer, and C. Franchini, Physical Review B 89, 075409 (2014).
  * Yu  _et al._ [2018] H. Yu, G.-B. Liu, and W. Yao, 2D Materials 5, 035021 (2018).
  * Nam and Koshino [2017] N. N. Nam and M. Koshino, Physical Review B 96, 075311 (2017).
  * Naik and Jain [2018] M. H. Naik and M. Jain, Physical review letters 121, 266401 (2018).
  * Emrem  _et al._ [2022] B. Emrem, R. Kempt, K. Finzel, and T. Heine, Advanced Theory and Simulations , 2200055 (2022).
  * Alden  _et al._ [2013] J. S. Alden, A. W. Tsen, P. Y. Huang, R. Hovden, L. Brown, J. Park, D. A. Muller, and P. L. McEuen, Proceedings of the National Academy of Sciences 110, 11256 (2013).
  * Elder  _et al._ [2015] R. M. Elder, M. R. Neupane, and T. L. Chantawansri, Applied Physics Letters 107, 073101 (2015).
  * Enaldiev  _et al._ [2020] V. Enaldiev, V. Zólyomi, C. Yelgel, S. Magorrian, and V. Fal’ko, Physical Review Letters 124, 206101 (2020).
  * Naik  _et al._ [2020] M. H. Naik, S. Kundu, I. Maity, and M. Jain, Physical Review B 102, 075413 (2020).
  * Xian  _et al._ [2021] L. Xian, M. Claassen, D. Kiese, M. M. Scherer, S. Trebst, D. M. Kennes, and A. Rubio, Nature communications 12, 1 (2021).
  * Zhang  _et al._ [2021] Y. Zhang, T. Liu, and L. Fu, Physical Review B 103, 155142 (2021).
  * Springer  _et al._ [2020] M. A. Springer, T.-J. Liu, A. Kuc, and T. Heine, Chemical Society Reviews 49, 2007 (2020).
  * Zhang  _et al._ [2020a] Y. Zhang, Z. Zhan, F. Guinea, J. Á. Silva-Guillén, and S. Yuan, Physical Review B 102, 235418 (2020a).
  * Venkateswarlu  _et al._ [2020] S. Venkateswarlu, A. Honecker, and G. T. de Laissardiere, Physical Review B 102, 081103 (2020).
  * Zhan  _et al._ [2020] Z. Zhan, Y. Zhang, P. Lv, H. Zhong, G. Yu, F. Guinea, J. Á. Silva-Guillén, and S. Yuan, Physical Review B 102, 241106 (2020).
  * Zhang  _et al._ [2020b] Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi, K. Ueno, E. Tutuc, and B. J. LeRoy, Nature Physics 16, 1093 (2020b).

See pages 1 of SI.pdf

See pages 2 of SI.pdf

See pages 3 of SI.pdf

See pages 4 of SI.pdf

See pages 5 of SI.pdf

See pages 6 of SI.pdf

See pages 7 of SI.pdf

See pages 8 of SI.pdf

See pages 9 of SI.pdf

See pages 10 of SI.pdf

See pages 11 of SI.pdf

See pages 12 of SI.pdf

See pages 13 of SI.pdf

See pages 14 of SI.pdf

See pages 15 of SI.pdf

See pages 16 of SI.pdf

See pages 17 of SI.pdf

See pages 18 of SI.pdf

See pages 19 of SI.pdf

See pages 20 of SI.pdf

See pages 21 of SI.pdf

See pages 22 of SI.pdf

See pages 23 of SI.pdf

See pages 24 of SI.pdf

See pages 25 of SI.pdf

See pages 26 of SI.pdf

See pages 27 of SI.pdf

See pages 28 of SI.pdf

See pages 29 of SI.pdf

See pages 30 of SI.pdf

See pages 31 of SI.pdf

See pages 32 of SI.pdf

See pages 33 of SI.pdf

See pages 34 of SI.pdf

See pages 35 of SI.pdf

See pages 36 of SI.pdf

See pages 37 of SI.pdf

See pages 38 of SI.pdf

See pages 39 of SI.pdf

See pages 40 of SI.pdf

See pages 41 of SI.pdf