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

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 (ta) is not yet well understood. Here, we monitor the change of twisted bilayer MoS2 characteristics as function of ta. We identify distinct structural regimes, 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 0.2-59.6deg. Some systems include up to half a million atoms, making structure optimization and electronic property calculation challenging. For 13<ta<47, the structure is well-described by a moir\'e regime composed of two rigidly twisted monolayers. At small ta (ta<3 and 57<ta), 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 ta approaching 0 and 60deg. Our results show at which ta range the characteristic features of the reconstruction emerge and give rise to exciting electronics. We expect our findings also to be relevant for other twisted bilayer systems.


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, θ, 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 • and 60 • , 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 lowenergy 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], MoS 2 [2,15], or MoS 2 /WSe 2 [16]. Several recent works also give strong experimental evidence of such reconstruction effects. [10,17,18] Depending on the value of θ (0 • ≤ θ ≤ 60 • ), the characteristics of the systems change, and different structural regimes are obtained: A θ of exactly 30 • results, as in bilayer graphene, in a quasicrystalline arrangement. [19] For large twist angles (θ → 30 • ), in the so-called moiré regime [10], structural relaxation of individual layers is minor. [12,15] Here, the weak oop deformation exhibits a sinusoidal modulation and decreases further for θ → 30 • . [2,3,20] At small twist angles, θ → 0 • and θ → 60 • , 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 θ → 0 • (see Fig. S1 for comparison to θ → 60 • ). 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 • or 60 • . 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., ≈ 1 • in tBL graphene [10,12,20], ≈ 4 • in tBL WSe 2 [27], and up to 13 • in tBL MoS 2 [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 θ = 2 • , there are 5000 atoms, this number increases to 20,000 at 1 • and reaches more than 100,000 for twist angles below 0.44 • . 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 [22,28] approaches.
In this work, we elucidate how the structure and electronic properties (band structure, local density of states (LDOS), local gap, and electronic transport) of tBL MoS 2 change with the twist angle, θ. To provide a com-plete picture, we obtain our results for a very large set of systems, allowing us to describe structural changes with θ almost continuously. We derive a complete characterization of the structural regimes in tBL MoS 2 and the transition between these regimes. We analyze the strong atomic reconstruction for very small θ 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 θ. Our approach can be adapted to other twisted vdW structures where atomic relaxation effects play a role.

METHODS
Initial tBL MoS 2 structures with periodic boundary conditions were constructed from ideal monolayer (ML) equilibrium geometries placed on top of one another with twist angles θ between 0 • and 60 • . 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 The R-type corresponds to a 0 • rotation, the H-type to a 60 • 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, ∆E, relative to the lowest-energy stacking H h h . Energies are given in meV per unit cell (uc).
[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 θ in the range of 0.2 • − 59.6 • . All fully-relaxed tBL MoS 2 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). While the interlayer distances in (t)BL structures optimized with ReaxFF and PBE+TS are in close agreement, we observe a consistent overestimation by about 0.35Å in PBE+MBD. The moiré patterns introduced by rigidly twisted layers change with θ, and the resulting size of the periodic unit strongly increases for small twist angles (θ → 0 • and θ → 60 • ). While at θ = 1.05 • , there are already 17,862 atoms in the unit cell, this value reaches 493,026 atoms in the largest studied structure with θ = 0.2 • . 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 θ is available as part of the ZENODO repository. [34] 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 1 /Å −1 along the reciprocal lattice vectors, and non-self consistent spin-orbit coupling (SOC). [41] Using bilayer and twisted bilayer MoS 2 systems, we verified that band structures calculated using ReaxFF-and DFToptimized structures are consistent (see Figs. S6, S7, S8, and Table S3). Electronic properties were analyzed from band structure and density of state (DOS) for twist angles between 3.89 • and 56.11 • . 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 × 2 × 1 for large supercells with θ ≤ 7 • and to 5 × 5 × 1 for θ > 7 • . 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 ≈ 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 × 5 × 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.

Structural and energetic properties
The primary reference systems needed for the analysis of tBL MoS 2 are the high-symmetry stacking configurations in untwisted bilayer (BL) MoS 2 , as shown in Fig. 2. The calculated structural properties of the BL MoS 2 stackings are consistent with previous reports in the literature [47] with more details given in the SI (see Table S4 Structural and electronic properties of tBL MoS 2 gradually change as a function of θ. As most properties of tBL MoS 2 show a continuous change with θ, 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 θ particular structural transitions emerge. Here, we provide a comprehensive characterization of structural properties as a function of θ. 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 θ and can be modeled with rigidly twisted layers. When going towards smaller θ, the characteristics of the system change as the influence of relaxation becomes more prominent. The structural elements start to appear, indicating the beginning transi- tion 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 θ, 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 θ are: (i) The interlayer distance d with its minimum, d min , maximum, d max , and mean, d mean , values as shown in Fig. 4(a). Here, gaps in the curves around θ = 21.79 • , 38.21 • , and 13.17 • 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 θ. (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 θ 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, ϵ, gives the local strain relative to an unstrained ML, the anisotropic strain field, ⃗ 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 ≈ 13 • and ≈ 47 • , an almost constant value of d 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 ≈ 0.15 mJ/m² per 1 • . The anisotropic strain field remains non-zero, as the anisotropy becomes zero only for untwisted bilayers at 0 • and 60 • (see Fig. S19). The displacement field being close to zero confirms that the ip reconstruction is negligible (see figure S18(b)).
In the transition regime, 6 • < θ < 13 • (equivalently, 47 • < θ < 54 • ), the atomic reconstruction becomes more prominent. Here, the interlayer distance shows a slow decrease of d mean with decreasing θ. This is reflected in the exfoliation energy, which slowly increases by ≈ 1.4 mJ/m² per 1 • . A strongly increasing interval ∆d = d max − d min (from 0.16Å at 13 • to 0.41Å at 6 • ) 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 θ is observed as a strong increase in their height and width, as shown in Fig. 4(c, d). While at 13 • , the base width of the nodes (solitons) is 13Å (8Å), it increases to 30Å (15Å) at 6 • . The nodes are the dominating structural element of the transition regime as they appear already at larger values of θ near 13 • (47 • ). Solitons appear only for lower θ near 6 • (54 • ), 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 θ (see figure S19).
In the soliton regime, 3 • < θ ≤ 6 • (54 • ≤ θ < 57 • ), the atomic reconstruction continues to become more prominent, and all three structural elements are present. The d max value reaches its upper bound, equivalent to R h h (H M h ) stacking (S@S type), at θ = 6 • (θ = 54 • ). This correlates with fully pronounced nodes of local highenergy stacking type. The height of the structural elements increases and reaches its maximum value at θ = 3 • (θ = 57 • ). Similarly, the width of nodes and solitons also increases. For θ down to ≈ 3 • (≈ 57 • ), the width of the nodes coincides with a 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 min decreasing further (see figure 4(a)). Their small size is also visible from d mean and the exfoliation energy, which show a small change with decreasing θ. 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 • < θ ≤ 3 • (57 • ≤ θ < 60 • ). 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 ≈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 min reaches its lower bound, equivalent to R M h and H h h stackings, showing that the stacking domains are fully pronounced. As ∆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 θ ≤ 1 • (θ ≥ 59 • ): 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 ≈15Å at 3.48 • . This is in excellent agreement with our results, giving w =22.7Å for 3.48 • , 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 M , as shown in Fig. 4(d). Decreasing θ 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 • , only 1.5% of the moiré cell can be accounted for the stacking domains, compared to up to 80% at 0.2 • . The formation of extended stacking domains is reflected in d 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 ≈ 21 mJ/m² per 1 • (see Fig. S17 and Table S5). Strain effects account for only ≈0.1% of the exfoliation energy, the deformation of individual layers ≈1%, making the interlayer adhesion terms the dominating contributions (see Fig. S16). The values of the exfoliation energy (0.29 J/m² for small θ and 0.24 J/m² for large θ) are in good agreement with results by Emrem et al., [51] who reported 0.33 J/m² for BL MoS 2 , calculated at random-phase approximation (RPA) level. The structural properties described for the domain-soliton regime are reflected in the strain fields, as shown for examples with θ → 0 • and θ → 60 • in Fig. 5. The stacking domains exhibit a constant isotropic strain ϵ 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 θ → 0 • and θ → 60 • deviates strongly for θ < 1 • and θ > 59 • (see Fig. 5 and H X h 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 θ → 60 • being stronger bound by up to 1.6 mJ/m² compared to θ → 0 • (see Fig. S20), resulting from a stronger adhesion between the layers. We observe that even though there are such substantial differences between θ → 0 • and θ → 60 • , there is an almost ideal symmetry θ i = 60 • − θ ′ 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 M . They find 1.2 • 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 M and the width of the solitons. They find values from 1 • to 1.6 • for tBL graphene, similar to the value of 3 • in the more rigid MoS 2 systems. [53] Maity et al. [15] find a maximum deformation of the individual layers in tBL MoS 2 at 3 • and a transition angle of 13 • . 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 • to 5.1 • for tBL WSe 2 .
This agrees with the soliton regime, where the structural elements start to vanish. Enaldiev et al. [54] find values of 2.5 • and 59 • 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 • ≤ θ < 60 • 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.

Electronic properties
We now investigate the impact of structural relaxation on the global and local electronic properties of tBL MoS 2 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.
Changes between the regimes are also found in the atom-projected band structures, as shown in Fig. 6(b, c). All tBL MoS 2 systems are indirect band gap semiconductors, with the VBM at Γ, similar to the untwisted bilayer limits. The CBM is positioned at K, which changes to be along the path between M and Γ with decreasing θ 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 θ (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 θ in the moiré regime, the band structure differences between relaxed ( Fig. 6(b)) and rigidly twisted With decreasing θ, the band gap slowly decreases and approaches the value of the S@Mo BL system, showing a linear dependence on d 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 MoS 2 with θ = 3.89 • , 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.
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)). It is also described in several other publications. [25,26,[55][56][57] 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 θ approaches 0 • , as shown in Fig. S22(a). Their bandwidth strongly decreases with decreasing θ (see Fig. S22(b)) and is expected to result in flat bands for the smallest θ. 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 correspond-ing 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 θ = 3.89 • 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 θ 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.
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 MoS 2 with θ = 3.89 • is reduced compared to the high-symmetry BL stackings with deviations of ≈50 meV from both R M h 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 [6,50,55,[59][60][61] 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 in-sights 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 θ is highlighted.

CONCLUSIONS
In this work, we studied the structural, energetic, and electronic properties of tBL MoS 2 by changing the twist angle θ 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 θ values for which the transition between the regimes occurs. For small twist angles 3 • < θ ≤ 6 • (54 • ≤ θ < 57 • ) 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 θ ≤ 3 • (θ ≥ 57 • ), 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 ≈ 13 • and ≈ 47 • inside the moiré regime, the model of rigidly twisted layers is applicable. Between these regimes for 6 • < θ < 13 • (47 • < θ < 54 • ), 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 θ. In rigidly twisted systems, flat bands form at small θ. 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 θ → 0 • 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 ≈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 MoS 2 systems. Otherwise, the predicted properties might be wrong not only quantitatively but also qual-itatively. 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 MoS 2 systems.

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. The research data associated with this work is available as a ZENODO repository, URL:
• Movie files, showing the change in the interlayer distance landscape and the strain fields with θ.
• The Python script used to create the plots of the interlayer distance, which was also the base for plotting the local electronic properties.
• The inputs and outputs of the geometry optimization calculations with ReaxFF, using LAMMPS.
• The inputs and outputs of the geometry optimization calculations at the DFT PBE level, using FHI-aims, to verify the results from ReaxFF.
• The inputs and outputs of the electronic structure calculations at the DFT PBE level, using FHI-aims, and at the DFTB level with QN13 parameters, [2] using Quantu-mATK.

Structure generation
Rigidly twisted BL MoS 2 structures without inherent interlayer strain were obtained from ideal MoS 2 monolayer (ML) supercells using a set of indices m, n, p, q so that fulfill the condition ⃗ A = ⃗ B and, thus, m 2 + mn + n 2 = p 2 + pq + q 2 . The lattice vectors ⃗ a i and ⃗ b i are of the primitive cell of ML A and B. The interlayer twist angle, θ, is then given as: and the number of atoms Z as: It is interesting to note that the size of the simulation cell and the moiré cell only agree for a few structures. The number of moiré cells N in the simulation cell is mostly larger than one. Due to the approach used to construct the initial structures from the indices m, n, p, q, we find that for θ < 30 • the value of N is given as: The different branches in Fig. S2 originate from this effect, with each branch belonging to tBL structures with a given value of N . The reason is that the size of the moiré cell is determined by the continuous change in the atomic registry and, as a consequence, gives the smallest periodic unit of the local interlayer separation landscape. The simulation cell, however, also needs to fulfill periodic boundary conditions in its atomistic structure, making it necessary, in most cases, to use a moiré supercell to obtain a valid tBL structure.

Geometry optimization with ReaxFF
Geometry optimization was performed for all structures using ReaxFF as parameterized   The relative cell constant (see Fig. S3(a)) shows a consistent behavior between ReaxFF and PBE+vdW with an exception at 60 • (H h h BL). The mean interlayer distance (see Fig. S3(c)) shows an excellent agreement between ReaxFF and PBE+TS. PBE+MBD gives a larger d, analogous to the BL stackings. However, d mean relative to the 0 • low-energy stacking type R M h (see Fig. S3(d)) shows a very good agreement between all methods. The main difference is the curvature of the individual lines, indicating differences in the deforma-tion behavior of the layers. This is approximated by the variation of the z-coordinate of the Mo atoms inside the layers, ∆z (see Fig. S3(b)). ReaxFF shows the largest ∆z and, thus, the strongest waviness of the layers. This can be understood by the layers being described as less rigid in ReaxFF compared to PBE+vdW. Consequently, the layers start to deform at larger θ and, thus, the θ values for the transition between the different regimes can be expected to vary.

Analysis of structural properties
Our analysis is based on the moiré cell, which we define as the smallest repeating unit in the local interlayer separation landscape of the simulation cell. We construct it from the positions of the nodes as areas of maximum local interlayer distance, placed at the vertices.
The solitons are located on the edges and the short diagonal of the moiré cell. This allows for a description independent of the original simulation cell and for an easier analysis of the structural elements. The general procedure is as follows: 1. Import structure (assumption: oriented in x, y-plane).

Identify Mo atoms and sort them into upper and lower layers.
3. Calculate the local interlayer distance via interpolation: • Extend Mo atoms into supercell to avoid edge effects during interpolation (extending by 20% of lattice constant or 35Å in each direction is sufficient).
• Interpolate z-coordinates of Mo atoms at points of interpolation grid (x, y) with suitable step size for both layers (interpolate.griddata from SciPy with method=cubic [8] is used in this work).
• Obtain local interlayer distance d from d(x, y) = z UL − z LL .

Identify the moiré cell:
• Find position and number N of local maxima (=nodes) in local interlayer distance.
• If N = 1, use the position of nodes as origin of moiré cell, and apply lattice parameters of original cell.
• If N > 1, use global maximum of d as origin of moiré cell (reference node), and use its nearest neighbor nodes to construct the hexagonal moiré cell.
• Interpolate z-coordinates of Mo atoms at (x, y) for both layers.
• Calculate d along ⃗ p.
6. Analyse structure of nodes and solitons: • Determine cross-section of nodes and solitons from interpolation along ⃗ p.
• Calculate height h from maximum of d in cross-section and global minimum of d inside the simulation cell. • against cross-section along ⃗ p.
• Calculate base width w as 4σ.
The results for an interpolation along the path ⃗ p to obtain the cross-section of the nodes and solitons are shown for two examples in Fig. S4. Fitting a Gaussian to approximate their shape is based on the findings of Carmesin et al. [9], who describe the curvature of TMDC BLs as being Gaussian. The consistently high R 2 values (R 2 ≥ 0.97) for all structures with pronounced structural elements confirm this approach. We use the base width 4σ as a characteristic property of Gaussians to determine the spatial extension of the structural elements. This approach tends to overestimate the size of nodes and solitons slightly. Both binding energy, E binding , and exfoliation energy, E exfol , are commonly encountered terms. They differ in sign, but are both calculated from the energy of the fully relaxed tBL system, E relaxed tBL , and the free-standing, unstrained monolayers, E ideal ML , using: The separation into individual energy contributions is fully mathematical and is not possible in real systems. It assumes three processes taking place consecutively (see Fig. S5).
from the total energy of the flat ML used to construct the rigidly twisted BL with ideal and optimized (strained) cell constant. In the next step, the individual MLs, still separated, are deformed by the relaxation of the atomic coordinates, introducing an intralayer waviness.
The deformation energy, E def , is defined as: with E optim ML,i as the total energy of the individual layers in the fully relaxed tBL system. Finally, the monolayers are bound by the interlayer interactions to form the fully relaxed tBL system. We further separate this adhesion contribution into two parts: The first part originates from the adhesion between the flat MLs in the rigidly twisted system. This interlayer (IL) adhesion energy, E IL , is obtained from: using the total energy of the rigidly twisted tBL system, E rigid tBL , with an interlayer distance of 6.30Å(intermediate between the high-symmetry BL stackings). The second part is the change in the interlayer adhesion when the individual layers relax. This change in the IL adhesion energy, E ∆IL is obtained from: To allow for a comparison between different tBL MoS 2 systems, the energies were normalized as follows: with a tBL as the simulation cell constant.

Twisted bilayer MoS 2
In tBL systems, the band structures calculated using geometries optimized with ReaxFF and PBE+vdW also agree well. Identical to the BL systems, the band structures based on the PBE+TS geometry match best with the ReaxFF-based band structures. Thus, the ReaxFF-optimized geometry can be used to avoid a computationally expensive geometry optimization using DFT.  The quantum transmission was computed at zero bias using the Landauer-Büttiker and NEGF approach [10,11] T where G indicates the total Green's function of the scattering region coupled to the electrodes. Γ α = −2Im(Σ α ) is the broadening function. Σ α with α = L, R are self-energies, calculated through self-consistent iteration. [12] Transmission pathways [13] were calculated to allow for an analysis option that splits the transmission coefficient into local bond contributions, T i,j . If the system is divided into two parts (A, B), the pathways are across the boundary between A and B, summing up the total transmission coefficient The direction of the pathway (either from i to j, or j to i) shows the flow of the electrons in the device. The device configuration used in these calculations is shown for an exemplary system with θ = 3.89 • in Fig. S12. The transport direction is from left to right.

STRUCTURAL AND ENERGETIC PROPERTIES
MoS 2 BL stacking types

Cell constant
The cell constant of tBL MoS 2 as a function of θ is studied by projecting the cell constant of the simulation cell, A, to the cell constant a prim of a primitive MoS 2 ML using FIG. S13. Projected primitive cell constant a prim in fully optimized tBL MoS 2 as a function of θ. The cell constant of the S@S and S@Mo BL stacking types is indicated for comparison.
Compared to the freestanding ML (a ML = 3.186Å), the cell constant in tBL MoS 2 is compressed by a maximum of 0.05%. Even though these changes are small, we found that optimizing the cell constant was necessary to obtain sufficiently relaxed structures that satisfy the convergence conditions.

Moiré periodicity
The moiré periodicity a M is the cell constant of the periodic unit in the moiré pattern. It can be evaluated numerically from the local interlayer distance landscape or analytically as  This was also observed by van Wijk et al. [14] in tBL graphene, who describe a strong decrease of the formation energy for θ < 2 • . We performed a piecewise linear fit:

Displacement field
The displacement field ⃗ v(⃗ r) [15] is the change of the atomic coordinates ⃗ x i between the fully relaxed and the initial, rigidly twisted system: The compression of the cell constant is neglected as it gives only a small relative change in the atomic coordinates Thus, the already lattice-optimized initial structure is used. Furthermore, the out-of-plane (oop) component is neglected as the initial interlayer distance can be chosen arbitrarily, adding ambiguity to its value.

Strain fields
The isotropic strain field ϵ is a scalar field: with d ML = 3.186Å. The anisotropic strain field ⃗ u is a vector field: They use the vector ⃗ d ij from Mo atom i to its next neighbor, j, and N = 6 as the number of next neighbors. The isotropic strain field can be understood as the local variation of a prim relative to the unstrained ML. Inside the moiré regime, it shows a situation similar to a sinusoidal modulation as described by van Wijk et al. [14]. The anisotropic strain field ⃗ u shows the local variation and the direction in which the Mo-Mo distances change the strongest. Contributions from the solitons dominate the bands below, which become the kagome-type bands at smaller θ. This links the emergence of kagome-type bands to a kagome lattice formed by the center of the solitons (see Fig. S22(c)).

Band gap vs. interlayer distance
The band gap E gap shows a linear dependence on the mean interlayer distance d mean :