Stacking domains and dislocation networks in marginally twisted bilayers of transition metal dichalcogenides

We apply a multiscale modeling approach to study lattice reconstruction in marginally twisted bilayers of transition metal dichalcogenides (TMD). For this, we develop DFT-parametrized interpolation formulae for interlayer adhesion energies of MoSe$_2$, WSe$_2$, MoS$_2$, WS$_2$, combine those with elasticity theory, and analize the bilayer lattice relaxation into mesoscale domain structures. Paying particular attention to the inversion asymmetry of TMD monolayers, we show that 3R and 2H stacking domains and a network of dislocations develop for twist angles $\theta^{\circ}<\theta^{\circ}_P\sim 1^{\circ}$ and $\theta^{\circ}<\theta^{\circ}_P\sim 2.5^{\circ}$ for, respectively, bilayers with parallel (P) and antiparallel (AP) orientation of the TMD unit cells and suggest how the 3R domain structure would manifest itself in tunneling characteristics of marginally twisted P-bilayers.

We apply a multiscale modeling approach to study lattice reconstruction in marginally twisted bilayers of transition metal dichalcogenides (TMD). For this, we develop DFT-parametrized interpolation formulae for interlayer adhesion energies of MoSe2, WSe2, MoS2, WS2, combine those with elasticity theory, and analize the bilayer lattice relaxation into mesoscale domain structures. Paying particular attention to the inversion asymmetry of TMD monolayers, we show that 3R and 2H stacking domains and a network of dislocations develop for twist angles θ • < θ • P ∼ 1 • and θ • < θ • P ∼ 2.5 • for, respectively, bilayers with parallel (P) and antiparallel (AP) orientation of the TMD unit cells and suggest how the 3R domain structure would manifest itself in tunneling characteristics of marginally twisted P-bilayers.
Beyond graphene and hBN, moiré superlattice effects have been observed in MX 2 /M X 2 heterostructures of transition metal dichalcogenides (TMD) [19][20][21][22][23][24][25], and it has been suggested theoretically [26][27][28] that twisted bilayers of TMD can undergo lattice reconstruction. Here, we determine parametric conditions for the formation of and the types of domain structures in twisted TMD homo-and heterobilayers, enriched by the lack of inversion symmetry of the individual 2D crystals. Different domain wall networks form for parallel-(P-bilayers) and antiparallel-(AP-bilayers) orientations of unit cells in the (a) Sketch of lattice relaxation across moiré supercell of marginally twisted AP-bilayers (θ < θAP ). Insets show local stacking configurations in 2H domain and MM seed. X(X ) and M(M ) label chalcogen and metal atoms, respectively, in the bottom (top) layer; (b) the left/right panel reveals side-view of 3R/2H-stacked TMD bilayers with the size of dots reflecting the layer-asymmetry of the electronic states at the conduction band edge. Central panel: the first reciprocal lattice star and the first Brillouin zone of the TMD bilayers with marked conduction (CBQ, CBK ) and valence band (VBΓ, and VBK ) extrema.
two layers, and we find crossover angles, θ • AP ∼ 1.0 • and θ • P ∼ 2.5 • , for the marginality of the twist and, then, dis-arXiv:1911.12804v1 [cond-mat.mes-hall] 28 Nov 2019 cuss how the resulting domain structures can be observed in scanning tunneling experiments.
In this study we employ the multiscale modeling: a combination of density functional theory (DFT) leading to interpolation formulae for adhesion energy, W P/AP (r 0 , d), between the layers at a distance d from each other and lateral offset r 0 , and elasticity theory for the lattice relaxation [27]. We perform this analysis for small misalignment angles, θ 1 (i.e θ • < 5 • ), and lattice mismatch δ 1. In this case energetics of local stacking can be described in terms of a lateral offset, r 0 (r) = θẑ × r + δr + u (t) − u (b) , between two aligned commensurate lattices, which varies across the moiré supercell and includes the effect of local deformations, u (b/t) (r), in the bottom/top layers. This multiscale approach enables us to overcome the system-size limitations of molecular dynamics ions [26,28].
For adhesion energy, we use the form, W P/AP (r 0 , d) = n f (P/AP ) n (d)e iGnr0 , where G n are the reciprocal lattice vectors of TMD. We truncate this sum at the first star of reciprocal lattice vectors, ±G 1,2,3 , (|G 1,2,3 | = G, Fig. 1) and set r 0 = 0 at the XX stacking for both P-and AP-bilayers. This choice -together with the D 3h lattice symmetry of TMD monolayers -suggests [29] that . Then, we inspect the adhesion energies for various H2DCs, computed using vdW-DFT with the optB88 functional [30] implemented in Quantum Espresso [31] for stacking configurations shown in Fig. 2. For P-bilayers, the most energetically favorable are configurations MX (r 0 = (0, − a / √ 3)) and XM (r 0 = (0, a / √ 3)), which have equal energies [32] and correspond to twins of a 3R bulk phase of a TMD. For AP-bilayers, 2H-stacking (r 0 = (0, − a / √ 3)) has the lowest energy (in agreement with the 2H bulk phase of these TMDs), rather than MM -stacking (configuration 5, r 0 = (0, a / √ 3)) suggested in Ref. [27]. Configurations 6, r 0 = − a /3(1, √ 3) and r 0 = − a /3(1, 0), are such that W P = f and W AP = f . The remaining two (2 and 4) have offsets and for AP(P)-bilayers. Using the DFT data for H2DCs shown in Fig. 2 and S1 in SI [30], we find the d-dependence of the factors f, f , f 1 , f 1 , and g 1 , and plot them in Fig. S2. By inspection, we find that f ≈ f and f 1 ≈ f 1 + g 1 , over the broad interval of interlayer distances that covers the minima of W P/AP (r 0 , d), and, then, consider the following factors to make an informed choice of functions, f , f 1 , and g 1 . (i) Coulomb potential of a lattice of ±q ions, whose potential decays exponentially, n α n e −Gnd e iGnr0 , with the distance from their plane (α 0 = 0 is due to the electroneutrality of each layer). This suggests a choice of g 1 = A 2 e −Gd for the first star of reciprocal lattice vectors  Table S1 in SI [30].
for chalcogen atoms in the outer (top/bottom) sublattices in each layer is determined by the exponential decay of atomic wave functions away from the plane, |ψ(z)| 2 ∝ e −|z|/ρ . Also, tunneling between the layers is suppressed for electrons with a larger in-plane wave number [33], n β n e −d √ G 2 n +ρ −2 cos (G n r 0 ), so that we (iii) Finally, following the earlier studies [34] of vdW interaction of H2DCs, which has a long-distance asymptotic, ∝ −C/d 4 , we combine the short range repulsion and long-range vdW attraction into f = − 3 n=1 C n /d 4n . Then, we use with ϕ P = π/2, ϕ AP = 0, and fit the values of parameters C 1,2,3 , A 1,2 , ρ to the DFT data listed in Table S1 in SI [30]. Lattice reconstruction in H2DCs happens when energy gain from the formation of favourable stacking overcomes elastic energy cost of strain produced by the local mutual adjustment of the two lattices, U = , are the first Lamé coefficient, shear modulus, and strain tensors related to the local in-plane displacements in top/bottom layer. We neglect the energy cost of flexural deformations allowing for the out-of-plane bending of the layers towards the optimal interlayer distance, d P/AP (r 0 ), for each offset r 0 . We describe the latter by expanding W P/AP into Taylor series about the minimum point, d 0 , of the zeroth Fourier Table S1 in SI [30]), and, then, find u (t/b) (r) that minimize energy, .
Here, g n = δG n − θẑ × G n are reciprocal lattice vectors of moiré superlattice, which we also use to expand u (t/b) (r) in Fourier harmonics up to the fourteenth reciprocal space star. Then, we minimize E with respect to those Fourier amplitudes numerically and obtain the displacements in each layer of the reconstructed bilayer. This method enables us to study moiré structure with Representative examples of lattice reconstructions in P-and AP-bilayers are shown in Fig. 3. For a larger angle, θ = 4 • (Fig. 3a,b) adhesion-induced displacements are small and the two layers behave as rigid lattices. For θ • < 1 • , twisted bilayers reconstruct into domain structures. For P-bilayers, Fig. 3d, each reconstructed moiré supercell comprises two equal area triangular domains of 3R(XM /MX ) stacking, separated by partial dislocations, with XX regions squeezed to the nodes on that partial dislocations network, Fig. 3g. For AP-bilayers, Fig. 3c, the reconstructed lattice features a honeycomb array of 2H domains separated by a dislocation network, where one half of the nodes are the seeds of the quasi-equilibrium MM phase (stacking configuration AP-5). As a quantitative measure for domain formation, we use a lateral distance, ℵ, between the closest metal and sulfur atoms in top and bottom layers (for ideal 2H and 3R domains, ℵ = 0) and compute its mean square, ℵ 2 /a 2 , over the supercell normalized by the TMD lattice constant, a. For rigidly rotated monolayers, ℵ 2 AP ≈ 0.14a 2 and ℵ 2 P ≈ 0.08a 2 . Upon the formation of 2H/3R domains, main contribution to ℵ 2 comes from domain boundaries so that the asymptotic behaviour, ℵ 2 ∝ 1/ ∝ θ, in Fig. 3h signals the formation of a domain structure at θ < θ P ≈ 0.044 (2.5 • ) and θ < θ AP ≈ 0.017 (1 • ). These quantitative estimates also explain why the molecular dynamics simulations performed in Refs. [26,28] for θ • > 3 • failed to establish the full picture of the lattice reconstruction in H2DCs, having interpreted MM areas in almost rigid AP-bilayers, Fig.  3a, as fully developed domains. Note that superlattice pattern -perfect domains and the dislocation network - For θ > θ AP/P the two layers behave as rigid, for marginally twisted bilayers (θ → 0) 2H for AP and 3R for P stacking domains emerge, separated by dislocations. (e-g) Intersections of dislocations, drawn over on the colour maps of the varied interlayer distance. (h) Crossover from rigid twist to the fully developed domain structure in P-and AP-bilayers, quantified using average value of parameter ℵ 2 described in the text.
The formation of domain structures takes place when the energy gain from the expanded 2H (for AP) or 3R (for P) areas overcomes the energy cost of domain walls. The latter are nothing but dislocations: screw dislocations in  2H and partial screw dislocations in 3R bilayers. The properties of such linear defects, analysed using energy functional (2), are shown in Fig. 4. Here, we set θ = 0 in Eq. (2) and replace u t − u b by ℵ, such that ℵ(−∞) = 0 and ℵ(+∞) = b AP/P , where b AP = a(−1, 0) is a Burgers vector of a dislocation in 2H-TMD (b AP = a), and b P = a(0, 1/ √ 3) is a partial dislocation Burgers vector in 3R bilayer (b P = a/ √ 3). A vector n = (cos α, sin α) determines the orientation of the dislocation line with respect to zigzag axis in the crystal. We find (see Fig.  4) that the energy of a dislocation in 2H-homobilayers is the lowest and the cross-sectional width narrowest for a zigzag orientation of the defect line (screw dislocation); for 3R-homobilayers, the most favourable orientation of the partial dislocation axis is along the armchair direction (partial screw dislocation). These choices coincide with the orientations of domain boundaries shown in Fig.  3. In addition, in Fig. 4, we show the profile of edge dislocations that would form in perfectly aligned (θ < δ) heterobilayers MS 2 /M S 2 and MSe 2 /M Se 2 .
Another feature of marginally twisted bilayers is a concentration of strain at the sites of their dislocation networks, highlighted in Fig. 3e-g. A generic feature of multivalley semiconductors (that is, both at CB Q and CB K ) consists in that inhomogeneous strain generates pseudomagnetic fields for electrons [63] with the alternating signs in time-reversed valleys, yy )]. For example, using parameters proposed [64] for electrons in ±K valleys of MoS 2 , we estimate that B * can reach values 20 − 30 Tesla over a 20 nm 2 area near the intersections of dislocations separating 2H-stacking domains (Fig. 5). These areas may represent interesting types of semiconductor quantum dots in 2D materials.
Finally, we point out that the described domain structures and dislocation networks result from the fabrication-induced twist in a bilayer. As such, they may slowly relax upon annealing, by dislocations escaping through the sample edges, resulting in larger-size domains of perfect 3R (for P) and 2H (for AP) stackings. Such heterostructure assembly may give an alternative (to the growth) method for fabrication of lattice-matched 2H and 3R MX 2 /M X 2 heterostructures. Acknowledgements.
Supplemental Material for "Stacking domains and dislocation networks in marginally twisted bilayers of transition metal dichalcogenides"

DENSITY FUNCTIONAL THEORY (DFT) FOR ADHESION ENERGIES IN TMD BILAYERS AND DATA ANALYSIS
In van der Waals-DFT (vdW-DFT) calculations of adhesion energies of TMD bilayers we neglected spin-orbit coupling, used a plane-wave cutoff of 816.34 eV (60 Ry), and kept the monolayer structure rigid, varying only interlayer distances and stacking (r 0 ). In modelling of heterobilayers (MoSe 2 /WSe 2 and MoS 2 /WS 2 ) we fixed the lattice constants to be the same in both layers (to ensure local commensurability), while keeping chalcogen-chalcogen distance in each monolayer equal to their individual monolayer values. For consistency, we performed these calculations using both (i) lattice constant of MoSe 2 (MoS 2 ) and (ii) lattice constant WSe 2 (WS 2 ) for both layers in MoSe 2 /WSe 2 (MoS 2 /WS 2 ) heterobilayer, comparing the differences. Such a comparison, Fig. S1e,f for sulphides (Fig. 2e,f of the main text for selenides), shows that these results are hardly distinguishable and their fitting using Eq. (1) produces almost identical parameters (the only noticeable difference is for the value of A 2 , as shown for lateral lattice constant a = 0.329/0.328 nm (a = 0.316/0.3153 nm) in Table S1). We use W P/AP (Eq. (1) in the main text) and DFT results for 2H-stacking to compute the frequencies of layerbreathing modes (LBM) for 2H and 3R TMD bilayers, and compare in Table S2 with those measured in Raman spectroscopy.

DFT BAND STRUCTURE ANALYSIS FOR 3R AND 2H BILAYERS
The DFT calculations of 2H and 3R bilayer band structures were carried out in the local density approximation (LDA), as implemented in the VASP code [68], with spin-orbit coupling taken into account using projector augmented wave pseudopotentials. The cutoff energy for the plane-waves was set to 600eV, and the Brillouin zone sampled by a 12x12x1 grid. The bilayers are placed in a periodic three-dimensional box with a separation of 20Å between repeated images to ensure no interaction between them. The structural parameters were taken from experimental measurements of bulk TMDC crystals [69,70].