Atomic reconstruction in twisted bilayers of transition metal dichalcogenides

Van der Waals heterostructures form a massive interdisciplinary research field, fueled by the rich material science opportunities presented by layer assembly of artificial solids with controlled composition, order and relative rotation of adjacent atomic planes. Here we use atomic resolution transmission electron microscopy and multiscale modeling to show that the lattice of MoS$_2$ and WS$_2$ bilayers twisted to a small angle, $\theta<3^{\circ}$, reconstructs into energetically favorable stacking domains separated by a network of stacking faults. For crystal alignments close to 3R stacking, a tessellated pattern of mirror reflected triangular 3R domains emerges, separated by a network of partial dislocations which persist to the smallest twist angles. Scanning tunneling measurements show that the electronic properties of those 3R domains appear qualitatively different from 2H TMDs, featuring layer-polarized conduction band states caused by lack of both inversion and mirror symmetry. In contrast, for alignments close to 2H stacking, stable 2H domains dominate, with nuclei of an earlier unnoticed metastable phase limited to $\sim$ 5$~$nm in size. This appears as a kagome-like pattern at $\theta\sim 1^{\circ}$, transitioning at $\theta\rightarrow 0$ to a hexagonal array of screw dislocations separating large-area 2H domains.


Abstract
Van der Waals heterostructures form a massive interdisciplinary research field, fueled by the rich material science opportunities presented by layer assembly of artificial solids with controlled composition, order and relative rotation of adjacent atomic planes. Here we use atomic resolution transmission electron microscopy and multiscale modeling to show that the lattice of MoS 2 and WS 2 bilayers twisted to a small angle, θ < 3 • , reconstructs into energetically favorable stacking domains separated by a network of stacking faults. For crystal alignments close to 3R stacking, a tessellated pattern of mirror reflected triangular 3R domains emerges, separated by a network of partial dislocations which persist to the smallest twist angles. Scanning tunneling measurements show that the electronic properties of those 3R domains appear qualitatively different from 2H TMDs, featuring layer-polarized conduction band states caused by lack of both inversion and mirror symmetry. In contrast, for alignments close to 2H stacking, stable 2H domains dominate, with nuclei of an earlier unnoticed metastable phase limited to ∼ 5 nm in size. This appears as a kagome-like pattern at θ ∼ 1 • , transitioning at θ → 0 to a hexagonal array of screw dislocations separating large-area 2H domains.
Moiré superlattices generated at the interfaces of two-dimensional crystals with a small relative twist or lattice mismatch have proved to be a powerful tool for controlling the electronic and optical properties of van der Waals heterostructures. In crystallographically aligned graphene/boron nitride (hBN) heterostructures, moiré superlattices lead to the formation of mini-bands for electrons and the Hofstadter butterfly effect [1,2]. Moiré superlattices in twisted bilayer graphene have already produced superconductivity [3], a Mott insulator-like state [4], and helical networks of topological edge states [5,6]. Significant breakthroughs have also been achieved in understanding the potential of moiré superlattices in transition metal dichalcogenide (TMD) twisted bilayers. Localisation of excitons by stacking-dependent modulation of the bandgap has been observed in WSe 2 /MoSe 2 [7,8] and WSe 2 /WS 2 [9] hetero-bilayer structures, while MoSe 2 /WS 2 heterostructures have demonstrated the formation of minibands of resonantly hybridized excitons, promoted by the close proximity of the conduction band edges in the two compounds [10].
The electronic structure of bilayers is known to depend critically on the local atomic stacking configuration. Despite its relative weakness, the coupling between neighboring van der Waals layers induces an atomic lattice reconstruction in the constitutent crystals, already ob-served in both graphene/hBN and twisted graphene bilayers [11,12]. For graphene bilayers with a small misalignment angle, θ ≤ 1 • , a monotonic variation of the interlayer lattice registry transforms into a pattern of triangular Bernal-stacked domains separated by a network of stacking faults observed from the diffraction contrast present in low magnification transmission electron microscopy (TEM) images [13,14,15]. The presence of lattice reconstruction has also been suggested [16,17,18] in TMD bilayers, but, despite enormous scientific interest in this system, an experimental study of lattice reconstruction in twisted TMDs is still missing. Here we use atomic resolution scanning transmission electron microscopy (STEM), accompanied by first-principle and multiscale modeling, to study how a network of domains with strong local commensuration develops at small twist angles in suspended WS 2 and MoS 2 homo-and heterobilayers. We find that, while lattice reconstruction in twisted bilayers close to 3R-type (parallel alignment) is morphologically similar to graphene, twisted homo-and hetero-bilayers of TMDs offer a broader diversity of physical properties prescribed by the lack of inversion symmetry in the constituent layers. This makes the resulting structure qualitatively different in twisted TMD bilayers not just morphologically, but also electronically. In particular, lattice reconstruction for P-oriented bilayers creates conditions for the formation of sizable domains of 3R stacking (which is rarely found in bulk TMD crystals [19,20]) and features intrinsic asymmetry of electronic wavefunctions, which has not been observed before in 2D materials. The latter observation was made using high-resolution conductive AFM tunneling experiments in agreement with Density Functional Theory (DFT) modeling. Moreover, strikingly different patterns emerge in 2H-type (anti-parallel) bilayers, where we find kagome-like dislocation network.
Bulk TMD crystals are composed of regularly stacked monolayers, each consisting of a triangular sublayer of metal atoms at the mirror-symmetry plane between two identical triangular chalcogen sublayers. In the 2H TMD polytype (space group P63/mmc), all metal (M) atoms are vertically aligned with the chalcogen (X) atoms in the nearest neighbor layers. In contrast, in the 3R polytype (space group R3m), chalcogens are aligned with the empty centers of the hexagons in the layers above and below. To recreate these polytypes in twisted homo-bilayers, we employ the tear-and-stamp transfer technique [21], where one half of an exfoliated monolayer flake is picked up and deposited onto the other half. In this case, perfect parallel (P) alignment of the crystalline lattices would result in 3R stacking, while an exact 180 • rotation (anti-parallel, AP) would produce 2H stacking order. A small deviation from a perfect P or AP alignment (to an angle θ ∼ 0.5 • − 3 • ) gives rise to a moiré pattern with a period, = a/2 sin θ 2 , much larger than the in-plane lattice constant, a, of the TMD. To fabricate hetero-bilayers with P or AP alignment we use a similar approach, where the parallel edges of the crystals are aligned during the transfer process [1]. The assembled stacks are then transferred onto custom-made support grids (see SI S1) to allow high-resolution STEM imaging. For scanning tunneling measurements, the assembled bilayers are transferred onto multilayer graphite pre-exfoliated onto an oxidized silicon substrate.
The presence of domains in the resulting moiré patterns is clearly seen from the diffrac-tion contrast in low-magnification STEM images, Fig.1  To understand the energetics of domain formation we computed interlayer binding energy densities, W (r 0 , d), for the various stacking configurations realized locally across the moiré supercell, encoded by a mutual shift r 0 (|r 0 | < a) of sulfur atom positions in the two layers ( Fig.2b,e) and interlayer distance d. The highest energy in both P and AP orientations is found for XX stacking (r 0 = 0). For P-type bilayers the most energetically favorable stacking is XM or MX (corresponding to r 0 = (0, a/ √ 3) and r 0 = (0, −a/ √ 3), respectively), which are symmetric with respect to the basal plane reflection and have the same energy. For the AP case, perfect 2H stacking has the lowest energy, with MM stacking being less favourable, which is reflected in the disparity of their sizes evident in Fig.2d.
The computed values of W (r 0 , d), shown in Fig.2b and 2e, are in a qualitative agreement with the earlier studies [16,17] and are well-described by an interpolation formula,  Table S1 in SI). We note that vdW-DFT overestimates the equilibrium layer spacing as compared to the experimentally measured bulk value, but the shape of W (r 0 , d) as a function of d is believed to be accurate [22]. Therefore, in the following we only analyze relative changes in the interlayer distance using a quadratic expansion for To achieve a quantitative description of the atomic reconstruction of the observed twisted TMD bilayers as a function of twist angle, we employ multiscale modeling, similar to the approach used in [18]. For two rigid monolayers, twisted to a small angle, θ ≤ 5 • , the interlayer shift, r 0 ≈ θẑ × r, periodically repeats the same local stacking configurations in consecutive moiré superlattice unit cells. Lattice relaxation introduces an additional shift such that, r 0 ≈ θẑ × r + u tb , where relative displacements are determined by the strain distribution, u We also take into account out-of-plane distortions of the layers by finding the optimal interlayer distance for the various local stacking configurations using the above-described expansion around d 0 , with a height profile illustrated in the insets in Fig.2b,e. At this point, we neglect the energy of bending deformation, estimated to be much less than other energy scales in the problem, especially for longer-period moire structures, as the energy cost of bending scales as inverse 4th order of the moiré pattern period. Note that AFM topography of a twisted bilayer (presented in SI S6 for a MoS 2 homo-bilayer) reveals height variation at the boundaries, e.g. between XM and MX domains, where the interlayer distance is expected to swell, according to eq.(1) and Fig.2b.
Note that such adjustment of the interlayer distance to the local stacking promotes the lattice reconstruction, as compared to the predictions made in [18].
Hence, we calculate the expected lattice reconstruction from the competition between elastic and adhesion energies, A 1 e −qd 0 cos (g n r + G n u tb ) + A 2 e −Gd 0 sin g n r + G n u tb + ϕ AP/P , d AP/P (r, u tb ) = 1 2ε [qA 1 e −qd 0 cos(g n r + G n u tb ) + GA 2 e −Gd 0 sin g n r + G n u tb + ϕ AP/P ].
Here, λ and µ are the first Lamé parameter and shear modulus of a TMD monolayer, g n = θẑ × G n are the supercell reciprocal vectors. To find the optimal u tb , we expand it in a Fourier series over g n and self-consistently solve a system of non-linear equations for their Fourier amplitudes (see S9 SI).
The results of such simulations for P and AP bilayers of WS 2 with θ = 1.29 • /1.09 • are shown in Fig.2c and f for comparison with the local lattice reconstruction determined using high-resolution STEM, Fig.2a  Notably, for θ < 1.5 • the MM domain seed does not grow larger than ∼5 nm, even as the moiré period increases, which leads to the kagome-like appearance of the 2H and MM stacked lattice reconstruction shown in Fig.1c,d. Furthermore, for θ < 1 • two adjacent 2H/MM boundaries merge together to form a 2H/2H boundary characterised by a screw dislocation with the Burgers vector parallel to the zigzag direction of the domain lattice (the deformation field calculated for such a single dislocation is described in SI, S10). As a result, for the smallest twist angle, the moiré pattern evolves into an array of dislocations running through a 2H stacked bilayer (such as in Fig.1d) and the energy gain from growing 2H areas fully overcomes the energy cost of forming the 2H/2H boundaries. This structure qualitatively differs from the simple triangular network found in P-bilayers, Fig.1b and Fig.S10, where MX and XM domains both represent identical lowest energy configurations (Fig.2b). The domain walls between MX and XM are partial screw dislocations with Burgers vector a/ √ 3 parallel to the armchair direction of the domain lattice, that converge to points of XX stacking (Fig.3g). The same trend was observed in hetero-bilayers (see Fig.1a and SI Fig.S8).
For both orientations the transition from almost rigid bilayers to fully developed domain structures can be traced using a "deformation parameter", determined as the mean square of the in-plane shift, δ, between the metal and sulfur atoms in the top and bottom layers (δ = 0 for perfect 2H and 3R stacking), averaged over one supercell for AP orientation and two triangular-shaped half supercells for P orientation. As δ = 0 inside perfectly stacked domains, the main contribution to δ 2 comes from domain boundaries (dislocations) so that it scales as the reciprocal of superlattice period, −1 . The theoretically computed δ 2 /a 2 , Fig.3h, shows a clear δ 2 /a 2 ∝ θ trend at small angles, signaling formation of ideal stacking domains, and a trend towards saturation starting at θ * P ≈ 2 • (θ * AP ≈ 0.9 • ) which we identify as the critical angles for formation of the dislocation network. Note that such dislocation networks may temporally evolve by dislocations reaching the sample edges, however, such recrystallization is slow, and, here, the twist angle (which determines the network period) is pinned by clamping the crystal to the grid or a substrate.
The presented analysis demonstrates that, below critical angles, twisted homo-and heterobilayers with a small lattice mismatch undergo a strong structural reconstruction into large area equilibrium domains separated by a network of dislocations. Such reconstruction necessarily changes electronic and optical properties of the twisted bilayers, as compared to the previously considered rigidly rotated crystalline lattices [7,8,9,10]. In particular, for θ < θ * P/AP the properties of a bilayer sheet would be dominated by intrinsic properties of 2H/3R domains. For Compared to twisted graphene structures, we show that twisted homo-and hetero-bilayers of TMDs offer a broader diversity of physical properties prescribed by the lack of inversion symmetry in the constituent layers. In particular, we find that lattice reconstruction for Poriented bilayers creates condition for the formation of sizeable domains of 3R stacking (which is rarely found in bulk MoS 2 and WS 2 crystals) and features intrinsic asymmetry of electronic wavefunctions which can be used for electric tuning of the optical properties of such twisted bilayers by a displacement field. For both P and AP-stacking, the formation of dislocation networks may also result in interesting properties due to the strong deformation fields produced at domain walls and their intersections (see SI S11 for pseudo-magnetic fields calculations). were performed with QSTEM [25], matching the experimental parameters for high resolution imaging using the FEI Titan microscope, using an effective source size of 1.0Å. Image filtering was performed using a patch based principal component analysis (PCA) denoising algorithm [26] implemented using the open source python package Hyperspy [27]. Briefly, small patches of the image, centred on local bright peaks, were compiled to form a 3D image stack. PCA denoising was then performed on the image stack, before the denoised patches were averaged at their original locations in the image to generate a denoised image. A high-pass filter was then applied to the image to remove long range signals associated with local surface contamination.

STEM imaging and analysis
See SI for details.

Calculations
For DFT calculations of the interlayer binding energy densities we used van der Waals density functional theory (vdW-DFT) with optB88 functional implemented in Quantum Espresso [28]. In these calculations we neglected spin-orbit coupling, used a plane-wave cutoff of 816.34 eV (60 Ry), and kept the monolayer structure rigid varying interlayer distances d and stacking configurations. Details concerning the multiscale modeling can be found in the Supplementary Materials. The density functional theory (DFT) bands were computed in the local density approximation, as implemented in the VASP code [29], with spin-orbit coupling taken into ac-count using projector augmented wave pseudopotentials. The cutoff energy for the plane-waves is set to 600 eV, and the in-plane Brillouin zone is sampled by a 12x12 grid. The bilayers are placed in a periodic three-dimensional box with a separation of 20 angstroms between repeated images to ensure no interaction would occur between them. The structural parameters were taken from experiments [30]. Gwyddion software package [31].