Complex Strain Scapes in Reconstructed Transition-Metal Dichalcogenide Moiré Superlattices

We investigate the intrinsic strain associated with the coupling of twisted MoS2/MoSe2 heterobilayers by combining experiments and molecular dynamics simulations. Our study reveals that small twist angles (between 0 and 2°) give rise to considerable atomic reconstructions, large moiré periodicities, and high levels of local strain (with an average value of ∼1%). Moreover, the formation of moiré superlattices is assisted by specific reconstructions of stacking domains. This process leads to a complex strain distribution characterized by a combined deformation state of uniaxial, biaxial, and shear components. Lattice reconstruction is hindered with larger twist angles (>10°) that produce moiré patterns of small periodicity and negligible strains. Polarization-dependent Raman experiments also evidence the presence of an intricate strain distribution in heterobilayers with near-0° twist angles through the splitting of the E2g1 mode of the top (MoS2) layer due to atomic reconstruction. Detailed analyses of moiré patterns measured by AFM unveil varying degrees of anisotropy in the moiré superlattices due to the heterostrain induced during the stacking of monolayers.


INTRODUCTION
Two-dimensional (2D) materials have attracted much attention in recent years due to their tunable optoelectronic properties that can be tailored by different external perturbations 1 and to the large amount of new and exotic features they exhibit when stacked in different combinations of multilayer 2D materials. 2,3 Transition-metal dichalcogenides (TMDCs) are room-temperature semiconductors that can be exfoliated down to the monolayer thickness, exhibiting strong photoluminescence (PL). 4 The extremely strong PL signal facilitates the identification of monolayers using standard optical characterization techniques. When combining two TMDC monolayers in vertical heterostructures, a charge transfer process occurs through the formation of a type-II band alignment. 5 The twist angle, θ, between the two monolayers determines the optical and electronic properties of the resulting material. 6,7 The recent emergence of twistronics based on the observation of superconductivity in twisted bilayer graphene (t-BLG) 8 and moiréexcitons in TMDC heterostructures 9−11 has motivated the study of other 2D homo-and heterostructures. For example, the quantum anomalous Hall effect was observed in t-BLG aligned to hexagonal boron nitride (h-BN), 12 and an incompressible Mott-like state of electrons was demonstrated in a long-period TMDC/h-BN moirésuperlattice. 13 In addition, ferroelectricity was observed in twisted h-BN 14 and in TMDC homo-and heterobilayers. 15,16 Small-angle twisted TMDC heterobilayers produce longrange-ordered moirésuperlattices where the interlayer exciton exhibits a relatively strong room-temperature PL 17 as well as more extravagant properties, such as sensitivity to circularly polarized light, 9 the presence of multiple excitonic resonances, 10 and long-lived excitons confined in the moirépotential. 18 Moreover, atomic reconstructions of the hexagonal lattice were linked to the formation of flat bands in the electronic structure of TMDC bilayers. 19−21 Moirépatterns can be visualized using high-resolution microscopy techniques, such as scanning transmission electron microscopy or scanning probe microscopy. 22 −25 In particular, the moirépatterns of TMDC heterobilayers were observed via piezoresponse force microscopy, Kelvin probe microscopy, and conductive atomic force microscopy. 26,27 While several studies have successfully observed lattice reconstruction in moireṕ atterns of TMDC homo-and heterostructures, 22,26 a fundamental understanding of the heterostrain occurring in the moirésuperlattices is still missing in the literature. This knowledge is crucial for the evolution of the field.
For TMDCs, recent studies have reported on the effect of varying the twist angle on the strain levels manifested in MoS 2 / MoS 2 homobilayers 28 and in MoS 2 /WSe 2 heterobilayers. 29 Experimentally, changes in shape and shift of the Raman inplane E mode were reported to be one of the indicators of the observed variations in strain. Through first-principles calculations, three regimes of atomic reconstruction were found in MoS 2 bilayers, where small twist angles (between 2 and 5°) exhibited the highest atomic-level strains. 28 However, the highest strains were observed with a twist angle of ∼0°in MoS 2 /WSe 2 heterostructures. 29 Such a difference in strains as a function of the twist angle is a consequence of the lattice mismatch concomitant to TMDC heterobilayers. Thus, the formation of moirépatterns is observed in heterobilayers at a 0°twist angle. In contrast, a small twist angle is required to form moirésuperlattices in bilayers without a lattice mismatch, such as in homobilayers; see Supplementary Figure S1.
Heterostrain is defined as the relative strain of the top layer with respect to the bottom one, which is deemed to be introduced during the fabrication process. 30 In t-BLG, it was found that the presence of small intralayer strains (between 0.1 and 0.7%) 30 might lead to the emergence of flat bands and correlated states. 31−33 In twisted MoS 2 bilayers, the strain was considered to be intrinsic to the bilayer formation due to lattice reconstruction at certain twist angles. 28,34 Furthermore, large uniaxial heterostrain levels were also shown to form onedimensional moirépatterns in both twisted bilayer graphene and TMDC heterobilayers. 35−37 Controlling and tuning the heterostrain is a viable option toward applications where anisotropy is required. For instance, interlayer PL of onedimensional TMDC moirésuperlattices can exhibit linearpolarization dependence as shown by Bai et al. 37 and could be implemented in polarization-sensitive optoelectronic devices. It is still unclear to which extent these strains arise as a result of the strong interaction between the two layers or whether heterostrain states are also introduced during the stacking process.
The fabrication of perfectly aligned heterobilayers composed of mechanically exfoliated monolayers is challenging. A straightforward method to precisely control the twist angle between two different TMDC monolayers is currently lacking. Although second-harmonic generation spectroscopy is the most widespread technique to experimentally determine the twist angle, this approach has some uncertainties 38 and is not accessible with standard optical characterization setups. Another characterization method that permits the quantitative measurement of a twist angle is highly desirable. Raman spectroscopy has traditionally been used to extract a plethora of useful information from a range of 2D materials. For instance, Raman measurements enable the determination of the number of layers in TMDC crystals, especially when analyzing the shear and layer-breathing modes in the low-frequency region. 39−41 A similar approach was utilized for quantitatively assigning the twist angle in twisted WSe 2 bilayers. 42 Parzefall et al. 38 reported on the correlation between the interlayer shear mode observed in the low-frequency Raman spectrum of MoSe 2 /WSe 2 heterostructures and the atomic reconstruction taking place in the system, which essentially indicates a near-0°alignment. The dependence of the interlayer exciton emission (both energy and intensity) with the twist angle has also been observed for different TMDC heterobilayers. 43,44 In the present work, we focus on MoS 2 /MoSe 2 heterobilayers. The ensuing intralayer strain states in reconstructed MoS 2 /MoSe 2 heterobilayers are examined via polarized Raman spectroscopy, AFM measurements, and molecular dynamics (MD) simulations. We show that the Raman in-plane E mode of MoS 2 can be used as a fingerprint of lattice reconstruction and, hence, of the presence of a near-0°twist angle in MoS 2 / MoSe 2 systems. Moreover, we extensively study the strain as a function of the twist angle by means of large-scale MD simulations. The analysis demonstrates the complexity of the strain distributions across the moirépatterns in MoS 2 /MoSe 2 heterostructures and evidences that lattice reconstruction occurs at small twist angles. Finally, we prove that upon applying uniaxial prestraining to the MoS 2 layer, the heterostrain that results from the formation of moirépatterns increases.

RESULTS AND DISCUSSION
Optical and Surface Characterization of MoS 2 /MoSe 2 Heterobilayers. We aim to study the strains occurring in reconstructed MoS 2 /MoSe 2 heterobilayers, that is, when two monolayers are stacked at θ ≈ 0°producing a large-periodicity moirésuperlattice. The samples were prepared by sequential transfer of mechanically exfoliated monolayers onto SiO 2 /Si substrates using a transfer stage with a controllable target rotation; see Experimental Methodology. Figure 1(a) shows the Raman spectra of a coupled MoS 2 / MoSe 2 bilayer stacked at θ ≈ 0°(in red) and a bilayer formed by uncoupled monolayers (in blue). The MoS 2 is on top in all cases. The coupling is induced by heating the sample above 80°C , as evidenced by the reduction of the interlayer spacing (Supplementary Figure S2) and the strong intralayer PL quenching. 45 The Raman spectrum of the uncoupled layers is just the combination of the noninteracting MoS 2 and MoSe 2 monolayers, whereas the spectrum of the coupled bilayer system (with θ ≈ 0°) exhibits a new mode emerging at ∼354 cm −1 that can be assigned to the A 2u mode of MoSe 2 . This mode becomes Raman active in the bilayer as predicted by group theory due to the reduction of symmetry in twisted heterostructures. 46 Also, it can be considered a fingerprint of the bilayer coupling in twisted MoS 2 /MoSe 2 heterobilayers, resembling the MoS 2 /WSe 2 case from Chiu et al. 47 We note that because the heterobilayer belongs to the C 3v point group, the correct notation for such a Raman mode should be A 1 . 46 (See Supplementary Figure S6 in Section S3 for the extended assignment of the Raman modes.) In addition, the Raman spectrum of the coupled bilayer unveils another interesting peculiarity related to the E mode (E 2g 1 for homobilayers), The PL spectra of the same samples are shown in Figure  1(b). The PL emission of the bilayer stacked at θ ≈ 0°u ndergoes a significant quenching of the intensity of both MoS 2 and MoSe 2 intralayer excitons due to the interlayer charge transfer. 50 When carefully comparing the energy of the MoS 2 excitons, a redshift of both A and B excitons can be observed. This is consistent with the observed behavior of the Raman features, as the intralayer tensile strains are expected to decrease the band gap energy, thereby redshifting the energy of the excitons in MoS 2 . 51 Another recognizable difference between the two samples is the appearance of a strong PL peak at ∼1.3 eV in the coupled bilayer. In agreement with previous reports, 52,53 we assign this peak to the interlayer exciton emission. These optical measurements provide a firm signature of the strong interaction between the two layers in the coupled systems, as also evidenced in the AFM topography image shown in Figure 1(c). A moirépattern with a periodicity of ∼8 nm is observed, which further corroborates the presence of a twist angle near 0°. We find such moirépatterns and Raman features in a set of five aligned MoS 2 /MoSe 2 samples; see Supplementary Figure S3 (Section S3). Additional AFM images of larger areas are also provided in Supplementary Figure S7 (Section S3). Note that in Figure 1(d), the moireṕ atterns are absent in the uncoupled samples.
The Raman spectra of a set of samples with large twist angles (θ > 10°) exhibit no E-mode splitting; see Supplementary Figure S8. This indicates that the observed Raman splitting is an exclusive feature of the aligned samples (θ ≈ 0°). We note that the activation of the A 2u mode of MoSe 2 remains unaffected by the twist angle and can be thus considered a fingerprint of the good coupling between the two monolayers. We verify that the E-mode splitting observed in MoS 2 /MoSe 2 heterostructures is universal for other heterosystems with similar lattice mismatches; see Supplementary Figures S12 and S13 (Section S4).
To shed more light on the type and levels of deformation in heterostructures stacked at θ ≈ 0°, we map out the dependence of the Raman spectra on the angle of linearly polarized scattered light (ϕ); see Figure 2(a). As expected, the intensity of the out-of-plane A 1 modes is highly dependent on the polarization angle, including the coupling-activated A 1 mode at 354 cm −1 . 54 We find that the E mode is a combination of (at least) two modes, similar to the peak splitting observed in uniaxially strained MoS 2 . 48,49,54 The intensities of the two peak components�centered at 381.5 and 383.5 cm −1 �are displayed in the polar plot of Figure  2(b). The two components can be assigned to E (−) and E (+) ; they arise due to the symmetry breaking of the doubly degenerate E phonon under anisotropic deformation. 55 The splitting and shift of the E mode�which is not accompanied by an apparent shift of the out-of-plane A mode�is characteristic of in-plane uniaxial strain in MoS 2 . 48 However, in pure uniaxial strain experiments, the intensity of the E-mode components can be fully controlled by modifying the polarization angle. 54 In our case, none of the peaks completely vanish at any polarization, thus evidencing that the strain is not purely uniaxial and that a more complex strain scenario develops in strongly interacting MoS 2 /MoSe 2 heterosystems.
Atomic Reconstruction of MoS 2 /MoSe 2 MoiréSuperlattices. We perform MD simulations to analyze the Here, δ is the lattice-parameter mismatch between the MoSe 2 and MoS 2 crystals (δ = d MoSed 2 /d MoSd 2 = 1.0375, in reference to the d MoSed 2 and d MoSd 2 values from the empirical potentials used in the MD simulations). However, for θ < 2°, we observe that the moirélattice constants obtained from FFT slightly deviate from the model (eq 1) due to the pronounced lattice reconstruction occurring in MoS 2 /MoSe 2 heterostructures. With θ ≤ 1°(i.e., near-0°conditions), the relaxed hexagonal patterns exhibit the largest anisotropy, as indicated by the marked fluctuations in the lattice constants extracted for the three directions of the hexagonal moirépatterns ( Figure  3(f1)).
Our MD results provide a mechanistic rationale for the atomic reconstruction occurring in MoS 2 /MoSe 2 heterobilayers. Figures 3(c,d) shows that the expansion of the AB and BA stacking domains plays a key role in the formation of moireś uperlattices with low twist angles. Furthermore, such a domain expansion leads to the shrinkage of the AA stacking sites in the relaxed configuration. A scheme of the stacking configurations formed in MoS 2 /MoSe 2 heterostructures is given in Figure 3(e). Similar structural reconstructions have been recently observed in WS 2 −WSe 2 moiréheterostructures. 57 With increasing θ values (2 < θ < 10°), the moiréperiodicity in our relaxed MoS 2 membranes drops as compared to the near-0°cases; see Figures 3(b,f2−f4). Lattice reconstruction, thus, is hindered with larger twist angles (θ > 10°), which produce moirépatterns characterized by a small periodicity (D < 2 nm; cf. Figure 3(b)) and mild out-of-plane deformations (Figure 3(f5)).
Twist-Angle-Dependent Strains in MoS 2 /MoSe 2 Heterostructures. The intralayer local strains, ε local , in our relaxed MD membranes are extracted from the components of the strain tensors defining the relaxation of the MoS 2 /MoSe 2 heterostructure; see Supplementary Section S7 for further details on the procedure. By means of this analysis, we gauge the evolution of the local strain distributions with increasing θ (from 0 to 30°) in twisted MoS 2 /MoSe 2 /Au heterosystems. To eliminate potential edge effects, the local strain levels are systematically evaluated in a central 30 × 30 nm 2 area of the relaxed MoS 2 membranes (marked with a blue square in the inset to Figure 3(a)).
The large-periodicity moirésuperlattices formed with very small twist angles (0 ≤ θ < 2°) lead to an intricate distribution of local strains in our relaxed MoS 2 membranes, where the largest ε local levels are predominantly located at AA regions in the reconstructed heterostructures; see Supplementary Figure S16 (Section S7). As a result, complex heterostrain states are obtained (Supplementary Figure S17). Our results indicate that the attainment of such complex strain patterns is a fingerprint of lattice reconstruction in TMDC heterosystems with θ < 10°; cf. Figure 4. Although the largest moireṕ eriodicity, D, in twisted MoS 2 /MoSe 2 heterosystems is anticipated for θ = 0°, 56 our ε local calculations indicate that, at first, the strain levels gradually increase in MoS 2 membranes with a twist angle that varies from θ = 0°to θ = 1°. We find that with θ = 1°, the MoS 2 /MoSe 2 heterosystem produces the largest average ε local value (∼1.15%). With increasing θ values (2 < θ < 10°), the average strain levels drop in our relaxed MoS 2 membranes, as compared to the near-0°cases. The small-periodicity moirépatterns associated with larger twist angles (θ > 10°) produce low ε local levels (with average values of ≈0.15%; cf. Figure 4).
While the maps of local strain values from Figure 4 exhibit only minor anisotropy and heterogeneity in the (nonprestrained) MoS 2 layers, a more revealing insight is provided by With θ = 0°, most of the Mo displacement vectors aim radially into the center of the AA sites. Only the displacement of the Mo atoms in between AA sites runs along the circumferential direction. However, with θ = 1°, the Mo displacements take on a direction that spirals toward the center of the AA site, resembling the "whirlpools" recently observed by Kim et al. 29 in MoS 2 /WSe 2 heterobilayers. Similar spiraling displacement patterns are observed with twist angles varying from 0.5 to 2°; see Supplementary Figure S18 (Section S8). These patterns gradually vanish with increasing θ. Under high twist angles (θ > 10°)�which hinder lattice reconstruction�an absence of such displacement fields in MoS 2 is observed; see Supplementary Figure S19 (Section S8).
The distributions of the strain components in the x and y directions (ε x and ε y , respectively) and their ratio (ε x /ε y ) are shown in Figure 5(c). In the ε x /ε y distributions, the green color corresponds to ε x /ε y = 1 ± 0.1, that is, to regions with mostly biaxial strain. Biaxial states are predominantly located at the AA sites, but their contribution vanishes with increasing θ; compare the ε x /ε y maps for θ = 0°and θ = 1°in Figure 5(c). Also, a localized strain distribution can be observed for approximately uniaxial strains (red hues) with − (ε x /ε y ) = ν ± 0.1, where ν is the Poisson's ratio set to 0.2 to cover the span of ν values (0.1−0.3) found in the literature. 58−60 All other colors (yellows, grays) belong to areas where ε x /ε y describes neither biaxial nor uniaxial straining states (i.e., shear strains).

Effect of Prestraining in MoS 2 /MoSe 2 MoiréSuperlattices.
By MD, we also evaluate the role of uniaxial prestraining in the resulting topography and strain maps of MoS 2 /MoSe 2 heterostructures with small twist angles, ranging from θ = 0°to θ = 2°. In our MD heterosystems, prestraining is introduced in the MoS 2 along the armchair (x) direction. Figure 6(a−d) displays the evolution of the moirépatterns together with the ε local maps with increasing prestrain levels under fixed θ = 0°. As expected, the ε local distribution and moirélattices become progressively more anisotropic with increasing prestraining. For comparison, the hexagonal moireć ells extracted from our relaxed MoS 2 membranes with different levels of prestraining are brought together in Figure  6(h). Figure 6(e) shows the experimental AFM-measured stiffness map of the coupled sample (displayed in Figure 1(c)), where the triangular-shaped AB (and BA) domains are clearly visible. For comparison, the resulting stacking domains found in the relaxed MD MoS 2 /MoSe 2 heterosystem with a uniaxial prestraining of 0.5% and θ = 0°are shown in Figure 6(f). Note that the mechanisms of stacking-domain shrinkage/ expansion in MoS 2 also lead to the complex atomic displacement fields observed in our MD heterostructures

ACS Nano www.acsnano.org
Article with near-0°twist angles ( Figure 5(b)). In this regard, the formation of large AB (and BA) domains is then limited to small twist angles. The AFM stiffness and MD topography profiles obtained from the linescans marked in Figure 6(e, f), respectively, are plotted together in Figure 6(i). Overall, these results reveal excellent agreement between the experiments and the MD simulations. Our FFT analysis of the MD MoS 2 topographies with twist angles θ = 0°and θ = 1°under different prestrain levels suggests the presence of three distinct values in the moireĺ attice constants (λ 1−3 ), thus revealing the onset of anisotropy in the moirécells as a result of prestraining (comparison of the λ 1−3 values from nonprestrained vs prestrained MoS 2 is shown in Figures 3(b) and 6(j), respectively). The same FFT analysis was performed for the AFM topography images, where the λ 1−3 were extracted under different scanning directions for several experimental samples; see Supplementary Figures S3 and S4 (Section S3). Averaging the corresponding λ 1−3 values for a given sample area allows us to reduce the effect of measurement fluctuations and piezo distortions (indicated as error bars in Figure 6(j)). This analysis reveals that the experimental moirépatterns also provide three unequal superlattice constants ranging from ∼7 nm to ∼8.5 nm. Such a large anisotropy provides another indication that our exfoliated MoS 2 membranes exhibit mild levels of prestraining (∼0.5%), which slightly vary between different samples and also between different areas in one sample.
To exclude possible AFM experimental artifacts, Supplementary Figures S10 and S11 (Section S3) display the outcome of measuring with parallel vs perpendicular scans and measuring with top-down vs bottom-up scans in the same sample area. The only negligible differences between the various scanning conditions confirm that the moirélattice distortion and the λ 1−3 values are real.
Simulated Raman Shift Distribution of Reconstructed MoS 2 . The strain distributions calculated from MD simulations are converted into Raman shifts of the MoS 2 E (+) (ω E + ) and E (−) (ω E − ) modes from the components of the strain tensor ε x , ε y , and ε xy (extracted from a central area of 30 × 30 nm 2 in the relaxed MoS 2 membranes) according to 61 where ω E 0 is the E-mode frequency without strain (taken as 385.5 cm −1 as a relative reference to the 1L MoS 2 on a SiO 2 /Si substrate; see the spectrum in Supplementary Figure S5 (Section S3)), γ E is the Gruneisen parameter, and β E is the shear-strain phonon deformation potential. We adopt γ E = 1.1 and β E = 0.78, which were obtained in a four-point bending (i.e., uniaxial) experiment conducted by Conley et al. 62 The "Raman spectrum" is then obtained by calculating the frequencies of the E (+) (ω E + ) and E (−) (ω E − ) modes for every atom contained within the simulated area (of 30 × 30 nm 2 ) and summing up all the contributions with the full-width-athalf-maximum (fwhm) set to 2 cm −1 for each calculated individual peak. It must be emphasized that the simplified simulations of the Raman shift merely capture the distribution of the strain-induced E-mode shifts. However, this allows for a comparative analysis of moirésuperlattices with large periodicity values that would not be accessible using common ab initio calculations.
The overall evolution of the simulated spectra of nonprestrained MoS 2 is plotted in Figure 7(a) for twist angles ranging from 0 to 30°. For twist angles >10°, the E mode is symmetrical. The fwhm reaches down 2.2 cm −1 , and the average value of ω E is 384.8 ± 0.2 cm −1 (indicated by the green dashed line in Figure 7(a)), which exactly matches the measured ω E of the noninteracting MoS 2 depicted in Figure  1(a). For comparison, the red dashed line in Figure 7(a) displays the ω E of the pristine MoS 2 (at 385.5 cm −1 ). For θ < 10°, the strain distribution in the Raman spectra takes on a markedly asymmetric shape, with two (or more) components clearly present.
In the nonprestrained MoS 2 , the higher-frequency component is always more intense than the lower-frequency one (Figure 7(b)), thereby strongly resembling the experimental Raman spectra. For the prestrained MoS 2 , the E-mode components evolve in a more complicated manner. However, they still maintain the asymmetric shape for low twist angles;  Figure 7(b). For comparison, the experimental spectra (obtained from the E-mode region in Figure 2(a)) are included in Figure 7(c) where the magenta-filled area outlines the range in which the spectral intensity reaches the given pixel for all analyzer angles and the gray dots represent their average. The positions of the two experimental E-mode maxima are also displayed in Figure 7(a) with magenta dashed lines. The net biaxial and uniaxial contributions to the simulated Raman spectra are marked in Figures 7(b, c) by green and red colors, respectively. The remaining areas under the curves (white) come from Raman shifts where ε x /ε y corresponds to various shear strains. More importantly, the spectrum areas filled by uniaxial states in the simulated Raman correlate with the changes in the experimental spectrum observed upon changing the analyzer angle; compare in Figure 7(c) the red areas (simulated spectra) with the magenta area (experimental spectrum). The agreement found between the simulated and measured strain distributions allows us to roughly estimate that the interrogated bilayer system has indeed a low twist angle (near 0°), as indicated by the E-mode splitting from Figure  7(c). Moreover, the resulting Raman shift indicates that our exfoliated MoS 2 membranes could be slightly prestrained.
We note, however, that the accuracy of the Raman-spectrum simulation may be hampered by the uncertainty in the parameters entering eq 2. There is a large span of γ E values in the literature from biaxial deformation tests, reaching down 0.60−0.64. 63,64 Even though we adopted the values reported by Conley et al. 62 for our calculations of the Raman shifts to maintain a single source for both parameters, the use of a significantly lower value for the Gruneisen parameter (i.e., γ E = 0.60) can change the outcome by shifting the E mode (∼1.5 cm −1 ) to higher wavenumbers and by narrowing the distribution (∼1 cm −1 ); see Supplementary Figure S20 (Section S9). The comparison between the experiments and the simulations is thus only tentative but it corroborates the observations made in the previous sections.

CONCLUSIONS
In this work, we evaluate local strains in reconstructed MoS 2 / MoSe 2 moirésuperlattices. Our study combines AFM and Raman measurements with MD simulations. We find that MoS 2 /MoSe 2 heterostructures with twist angles between 0 and 2°exhibit the highest levels of local strain due to specific structural reconstructions of the stacking domains that occur in the MoS 2 . Lattice reconstruction is hindered for heterostructures with larger twist angles (>10°). The MD simulations indicate that the formation of moirépatterns in TMDC heterostructures leads to a complex distribution of local strains with a combined deformation state of uniaxial, biaxial, and shear components. These findings are supported by our polarization-dependent Raman spectroscopy measurements. The twist-angle-dependent lattice reconstruction and highly localized strain distribution are in line with the previously reported observations of θ-dependent emission energy of the interlayer excitons 43,44 and confirm the potential for manipulating moiréquantum wells simply by proper oriented stacking. Furthermore, anisotropy of the moirésuperlattice may result from the introduction of heterostrain during the stacking process. The experiments and simulations suggest that moiréanisotropy can be tuned by controlling the uniaxial prestraining of the top layer, thereby enabling the development of TMDC moiréoptoelectronic devices with tailored directional properties.

EXPERIMENTAL METHODOLOGY
Sample Fabrication. Our MoS 2 /MoSe 2 heterostructures were prepared by mechanical exfoliation of bulk crystals (HQ graphene). The MoS 2 and MoSe 2 monolayers were separately exfoliated on polydimethylsiloxane stamps (Gelfilm by Gelpak) and sequentially transferred to SiO 2 /Si substrates (300 nm SiO 2 thickness) using a drytransfer technique. 65 Exfoliated monolayers were identified using PL spectroscopy.
Optical and Surface Characterization. Raman and PL measurements were performed in a LabRAM HR Evolution spectrometer (Horiba Scientific) with 532 nm laser excitation using a 100× objective (0.8 NA). Diffraction gratings of 1800 and 150 l/ mm, giving pixel-to-pixel resolutions of 0.5 and 7.0 cm −1 , were used for the Raman and PL characterizations, respectively. For all measurements, the laser power was set between 40 and 80 μW. The angle of the linear polarized incident light was fixed and the angle of the scattered light was selected by rotating the analyzer.
AFM topography and stiffness images were taken with a Bruker Dimension ICON in PeakForce Tapping mode using Scanasyst Air probes (Bruker Corp.).

COMPUTATIONAL METHODS
We perform large-scale MD simulations of MoS 2 /MoSe 2 /Au heterosystems employing the LAMMPS code. 66 The details of the construction of the all-atom heterosystems are given in Supplementary Section S5 and Figure S14. The Stillinger-Weber (SW) potential is employed to describe the intralayer atomic interactions in the MoS 2 and MoSe 2 membranes. We adopt the SW parameters obtained by Kandemir et al., 67 which generate consistent results with density functional theory. The constituent layers (the top MoS 2 membrane, the intermediate MoSe 2 membrane, and the bottom Au substrate) are mutually adhered through weak-bond van der Waals (vdW) forces. To model these interactions, the archetypal 12−6 Lennard-Jones (LJ) potential�Supplementary eq (S4) (Section S10)�is used to describe the vdW-type adhesion between the MoSe 2 layer and the substrate. The values of the LJ parameters for the Mo−Au and Se−Au interactions are given in Supplementary Table S2 (Section S10). The LJ cutoff distance is set to 8.5 Å.
In addition, to model MoS 2 −MoSe 2 interlayer forces, we employ the Kolmogorov-Crespi (KC) potential 68 parametrized for TMDCs by Naik et al. 69 This model incorporates a stacking-dependent term that enables interlayer sliding processes to generate different stacking configurations with different binding energies. The development of moirépatterns in 2D materials is fundamentally affected by this stacking dependence of the binding energy, which in turn cannot be captured by the simpler LJ model. 70 We adopt the z-normals simplification (KC-z), 69 and we set a KC cutoff radius of 14 Å. An assessment on the formation of stacking domains vis-a-vis the employed interlayer (LJ vs KC) potentials is provided in Supplementary Sections S11 and S12 and Figures