Mechanical Properties and Strain Transfer Behavior of Molybdenum Ditelluride (MoTe₂) Thin Films

The discovery of, and pioneering works on graphene introduced the remarkable potential of two-dimensional (2D) materials [1,2]. Transition metal dichalcogenides (TMDs) are a class of 2D semiconductors with varying bandgap depending on chemistry, unlike pristine graphene which has the undesirable zero bandgap for microelectronic applications [3]. Of these, Group VI TMDs (MX₂ where, M = Mo, W and X = S, Se, Te) are increasingly studied for their multifunctional applications ranging from highly tunable mechanical and electronic properties to exotic quantum phases and superconductivity [4]. For instance, TMDs with sulfur as the chalcogen atom (MoS₂, WS₂) have well-known tribological applications such as lubrication and wear resistance [5,6] in addition to their electrical properties. Molybdenum disulfide MoS₂ has been reported to sustain relatively large deformation (11% fracture strain) indicating flexible mechanical behavior appealing for electronic devices [7]. Coupling between such unique mechanical and electrical properties of layered TMDs highlights the need for further exploration [8,9]. Moreover, the existence of different phases with large variations in physical properties offers another degree-of-freedom for these materials [10].

Here, we focus on the 2H phase of molybdenum ditelluride, MoTe₂, which has an indirect bulk bandgap close to Si (∼1 eV), making it promising for on-chip integration [11]. In recent years, MoTe₂ has been explored for wide-ranging applications including as photodetectors [12], in energy storage [13], and as novel piezoelectrics [14]. A thickness-dependent transition from indirect to direct optical bandgap in its monolayer form provides additional tunability [15]. Compared to MoS₂, MoTe₂ has a smaller first-principles calculated energy difference between semiconducting 2H and metallic 1T′ phases [16], suggesting lower-energy switching between these two phases in MoTe₂. Many techniques have been used to induce phase transitions in MoTe₂ including laser irradiation [17], plasma treatment[18], defects [19,20], and mechanical loading [21]. Among these, mechanical straining has shown a room temperature phase transition at as low as 0.2% strain [22]. However, devices based on monolayer or few-layer TMDs are challenging to characterize in experiment due to the nanometer scale thicknesses involved. Additionally, MoTe₂ is highly reactive to oxygen creating further complexities in measurement of its mechanical properties [23]. Despite these challenges, some experimental characterizations of mechanical properties of TMDs have been reported [24,25]. There have been recent advances in novel approaches to strain engineer TMD thin films for non-volatile room-temperature phase change transistors. In particular, a MoTe₂ device has been demonstrated to surpass the performance of conventional field-effect transistors by eliminating static and dynamic power consumption problems [26]. Yet, many aspects of the deformation behavior and strain transfer in these devices are still unknown.

Computational approaches can offer valuable insights into the fundamentals of deformation in these thin-film structures. Within this context, while first-principles calculations such as density functional theory (DFT) inherently incorporate quantum mechanical information, feasible calculations are limited to a few thousand atoms that imposes limitations on simulating nanoscale experiments. To understand larger-scale features such as nanoscale deformation mechanisms, alternative techniques need to be employed. Molecular dynamics (MD) is an appropriate tool to extend modeling capabilities to nanometer and micron length scales. Furthermore, MD allows analysis of macroscopic factors such as temperature and the application of complex mechanical loads such as nanoindentation for the effective mapping of nanomechanical behavior [27]. The performance of MD models largely depends on using an interatomic potential that accurately describes the forces between atoms. Empirical potentials have been successfully used in MD to compute several properties of 2D materials [28,29]. The Stillinger–Weber (SW) is a suitable potential for covalently bonded monolayer TMD systems as it can describe the dominant bond-stretching and angle-bending interactions along with their nonlinear effects [30]. Additionally, this potential is more computationally efficient than bond-order potentials such as the reactive empirical bond order (REBO) and Tersoff [31,32]. In many cases, the SW potential can be an order of magnitude faster than bond order potentials [33], allowing its use in larger structures while maintaining reasonable accuracy.

This paper aims to describe the deformation behavior in MoTe₂ thin films. We first develop parameterization of an SW interatomic potential that can accurately replicate mechanical properties of MoTe₂ using MD. The parameters are fitted to equilibrium structures, elastic constants, and vibrational frequencies obtained from DFT calculations. Comparison with mechanical properties such as Young’s modulus at room temperature, uniaxial stress-strain behavior, and earlier micron-scale experimental nanoindentation provides validation of the potential. We then characterize mechanical behaviors by applying uniaxial deformation until failure using MD simulations. We consider macroscopic quantities such as crystal orientation and temperature and show their influence on the nanoscale deformation mechanisms of MoTe₂ monolayers. Analysis of stress-strain responses beyond the linear elastic regime reveals directional anisotropy in mechanical behavior and failure properties. Finally, we combine atomistic models and Raman spectroscopy experiments to quantify strain transfer behavior in MoTe₂ multilayered films stressed by a capping layer. Our findings help in further understanding the strain-induced properties of MoTe₂ thin films in strain-engineered electrical devices.

MoTe₂ films have a layered structure, but in contrast to graphene, each layer consists of three atomic sub-layers where the Mo atoms reside between two Te sub-layers forming a Te–Mo–Te sequence. Here, we focus on the most stable form of MoTe₂ having a 2H hexagonal stacking sequence, i.e., Te (Mo) atoms of one layer reside directly on top of the Mo (Te) atoms of the following layer as shown in Fig. 1(a). The intralayer atoms (Mo and Te atoms within a single layer) form a trigonal prismatic coordination where each Mo atom has six Te atoms as first neighbors and each Te atom is directly bonded with three Mo atoms [34]. Atoms within a single layer form strong covalent bonds and interlayer atoms are coupled by weak van der Waals interactions. The covalent bonds are predominantly governed by three bond-stretching (Mo–Te, Mo–Mo, and Te–Te) and three angle-bending motions (two Mo–Te–Te angles with Mo being the central atom and one Te–Mo–Mo). For the convenience of creating the structures in MD, we used an equivalent orthogonal unit cell (⁠$b=3a$⁠, and α = β = γ = 90 deg) instead of the primitive hexagonal unit cell.

Density functional theory (DFT) calculations were carried out using the Vienna ab initio Simulation Package (vasp)[35–38] with projector augmented wave (PAW) pseudopotentials [39,40] including Mo 4spd 5s and Te 5sp states as valence. The plane-wave energy cutoff was 500 eV, and a Gamma-centred 9 × 9 × 2 k-point grid was used. These provided energy convergence to 0.5 meV per formula unit and convergence of the elastic constants to within 1%. The convergence criterion for the electronic self-consistent loop was set to 10⁻⁸ eV. The optB86b-vdW [41] exchange-correlation functional was used to account for the disperse interlayer interactions, which gave lattice parameters within 1% of experiment. Monolayers were separated by a 20 Å vacuum region. Structures were optimized until the residual forces on the ions were less than 0.0001 eV Å⁻¹. Elastic constants were calculated with the finite difference method implemented in vasp. Phonon dispersion spectra were calculated using the finite displacement method implemented in phonopy[42] using supercell sizes 4 × 4 × 1 for the bulk structure and 5 × 5 × 1 for the monolayer structure, and non-analytical corrections were applied. For calculating stress-strain relationships, lattice parameters were fixed to a uniaxial strain value and ion positions were allowed to relax until the residual forces were less than 0.0001 eV Å⁻¹.

We used LAMMPS [43] open source MD program (August, 2019 version) to study several mechanical behaviors for both monolayer and multilayer single-crystal films. Before the application of deformation, structural parameters were optimized using energy minimization (convergence of energy norm below 10⁻¹⁴). A uniaxial tensile deformation was applied to monolayer MoTe₂ with respect to different crystal orientations defined by the chiral angle (θ) between loading direction and zigzag crystal direction as shown in Fig. 1(b). Boundary conditions were kept periodic in the planar directions (X − Y) and a vacuum region (∼80 Å) was inserted in the out-of-plane (Z) direction to avoid interaction with periodic images. Monolayer film sizes were approximately 30 nm × 30 nm in the armchair and zigzag directions, whereas the size was adjusted accordingly to ensure periodic boundary conditions for structures with other orientations. The MD time-step was 0.5 fs, and the systems were equilibrated using the NPT ensemble [44] at the desired temperatures. For uniaxial deformation, a constant engineering strain rate of 10⁹ s⁻¹ was applied.

For strain transfer length scale simulations, a molecular statics method was used which is inherently at 0 K temperature. An in-plane biaxial strain was applied only to the top layer of multilayered nanosheets. We constructed the loading condition so that a fully encapsulated sample with a thin film stressor was simulated, as was used in experiment. A biaxial strain of Δε = 0.14% (ε = $ϵx2+ϵy2$⁠) was incremented up to a ε = 0.85% strain magnitude. In each step, atoms in the top layer were kept stationary, and energy minimization (convergence of energy norm below 10⁻¹⁶) was performed between loading steps to allow the transfer of strain to subsequent layers. Using energy minimization ensures that all atomic displacements are solely caused by atomic interactions among layers and not due to inertial effects arising from high MD strain rates. For comparison, strain transfer in MoS₂ films with two to seven layer thicknesses was simulated using the same setup except with a REBO potential obtained from the literature [31,45]. Further details and original results on strain engineering of MoS₂ are described in a recent work by the authors [46]. OVITO open visualization tool was used to visualize the structures [47]. Percentage change in the length of each layer was used for calculating the layer-by-layer average strain magnitudes. Atomic strain was used for the local strain information along with its spatial distribution [48] through the thickness, which are not accessible from average quantities. However, atomic strain is a statistical measure and has more fluctuations near the free surface. Therefore, we use the atomic strain to show spatial distributions and length changes to report the average strain to avoid inhomogeneities associated with the statistical nature of atomic strain.

In order to compare the performance of the developed SW potential with experimental nanoindentation results, we conducted finite element analysis (FEA) nanoindentation simulations of circular MoTe₂ membranes in commercially available FEA program abaqus. FEA models used non-linear elastic constitutive behavior with the uniaxial stress versus strain data obtained from MD results (by the SW potential) defining the material property. We considered MoTe₂ membranes of at least 7 nm thickness (approximately ten layers) with a diameter of 1 μm, in order to be comparable with experiments. A structure of this size would contain too many atoms to be within the scope of MD. Suspended membranes were modeled with a total of 4912 four-node plane stress elements and were immobilized at the periphery to constrain displacements. A rigid spherical indenter was used with frictionless interaction with the membranes. The incremental displacement was 0.1 nm per load step and the fracture point was identified as the step beyond which equilibrium convergence could not be achieved.

MoTe₂ samples for Raman spectroscopy were mechanically exfoliated from the bulk crystal onto flat MgO single-crystal substrates (average roughness < 0.125 nm) inside a humidity-controlled environment (<1 ppm H₂O and O₂). The samples were then placed inside a vacuum chamber for thin film stressor deposition. Thin film stressors were e-beam evaporated at 5 × 10⁻⁶ torr, with evaporation rates kept between 1 and 2 Å/s. Raman spectra of the given samples were acquired with a WITec Alpha300R Confocal Raman Microscope. The power of the 532 nm Raman laser was carefully monitored to be 0.5 mW, in order to avoid heating-induced damage to the samples. A spectrometer grating of 1800 l/mm was chosen to ensure a spectral resolution of ±0.1 cm⁻¹. Finally, the peak positions of the given Raman spectra were determined by fitting to a Voigt spectral profile.

In order to check the accuracy of our newly developed potential, we evaluate its performance on known quantities starting with elastic properties such as the elastic constants and Young’s modulus of MoTe₂. Elastic constant calculations are performed at 0 K in LAMMPS to be comparable with DFT. The Young’s modulus is extracted from uniaxial stress-strain data up to 1% strain at 300 K for comparison with experimental results [54]. Based on the comparison demonstrated in Table 5, the developed empirical potential accurately predicts elastic properties with a maximum of 4% error.

The phonon dispersion spectrum encompasses a material’s vibrational properties, and the longitudinal and transverse acoustic phonon modes are also related to the elastic properties of a material. We generate the phonon dispersion curves of monolayer MoTe₂ using GULP along the high symmetry Γ − M direction of the Brillouin zone. The phonon dispersion spectrum calculated by the developed SW force-field is shown in Fig. 2, and for comparison, corresponding frequencies from DFT calculations are also presented. Figure 2 clearly illustrates that acoustic phonon frequencies from the proposed SW potential closely follow the DFT results, especially near the Γ-point where they are closely associated with elastic properties. Furthermore, the stress-strain data obtained by using our potential in MD is shown in Fig. 3. It is clear that even up to 4% strain, the stress-strain curves predicted by MD in both zigzag and armchair loading closely follow the DFT results, pointing to the reliability of our empirical potential away from the equilibrium configuration. At higher strains, the MD results start to deviate as the formulation of the SW potential is inherently different from DFT at that range.

As a final validation step of the potential, we implement the uniaxial response obtained from MD to model the material for micron-scale suspended nanoindentation in FEA. In Fig. 4, the resultant load-displacement data using a 30 nm indenter tip is presented. As a comparison, recent experimental data reported for atomic force microscopy (AFM) nanoindentation of a 2H MoTe₂ thin film are also plotted [54]. It can be seen that the FEA prediction aligns well with the experimental measurements even at large displacements. The breaking force predicted by FEA is 1.46 μN compared to 1.55 μN observed in experiment equaling to a <6% difference. The breaking strength (σ₀) of the film is 3.95 GPa from FEA whereas the experimental data reported a strength of 4.6 GPa for a film with 6.7 nm thickness. There is a similar match between FEA and experimental results when the indenter size is changed to 100 nm radius (Fig. S1 in the Supplemental Material available on the ASME Digital Collection). The close agreement of the nanoindentation force-displacement curve, as well as the fracture quantities with AFM experiments from literature, is a strong implication that the interatomic potential developed is reasonably accurate well beyond the linear elastic regime. With these sets of supporting evidence, one can be confident in using the potential for predicting general mechanical properties of MoTe₂, which is the aim for the remaining parts of the paper.

In this section, we use MD simulations to examine the mechanical behavior of MoTe₂ monolayers under uniaxial loading conditions, including the effects of loading direction and temperature.

To explore the influence of the loading direction on the stress-strain response, we apply uniaxial tension in orientations varying between the zigzag (θ = 0 deg) and armchair directions (θ = 30 deg), applying the deformations up to the point of failure. Due to the hexagonal symmetry of the crystal, orientations beyond 30 deg are symmetrically equivalent. In Fig. 5, we present the stress-strain curves at θ = 0 deg, 9.5 deg, 21 deg, and 30 deg orientations, simulated at 1 K. The four stress-strain responses are overlapping in the linear elastic region (up to 4% strain using the 0.2% offset method for the zigzag orientation), indicating that the material is elastically isotropic. However, the complete stress-strain curves up to failure are quite different revealing anisotropy in nonlinear mechanical behavior. These findings are in accordance with previous reports of first-principles calculations showing different stress-strain responses for MoTe₂ in armchair and zigzag loading [55,56].

Figure 5 also shows orientation-dependent behavior in the failure properties. Plotted in Fig. S2, there is an approximately linear increase in fracture strain with θ, and the fracture strain in the armchair direction (34.7%) is considerably higher than the fracture strain in the zigzag direction (22.6%). Similarly, the material has a higher ultimate tensile strength (UTS) under armchair loading than zigzag. Interestingly, there is no yield point in the stress-strain response under zigzag loading, whereas distinct yield points (marked by a sharp drop in stress) are evident in all other directions (Fig. 5). Moreover, there is a decreasing trend in yield stress and corresponding strain with increasing θ. Similar orientation-dependent yield behavior has been reported in other 2D materials [57].

The anisotropy in the mechanical response can be attributed to the alignment of crystal bonds relative to the loading direction. In the case of zigzag loading, 33% of the bonds (two out of six in the hexagonal motif) are perpendicular to the stretching direction and thus cannot contribute to load-bearing. As a result, a large proportion of bonds are left unstretched while the remaining bonds are subject to additional local strain, as illustrated in a contour plot in Fig. S3(a). In contrast, the armchair configuration has 33% of bonds aligned parallel to the applied load, and none are directly perpendicular. Consequently, all bonds can stretch and contribute to load-bearing (Fig. S3(b)). Unfavorable bond alignments reduce the capacity of the structure to sustain macroscopic deformation, resulting in fracture at lower strains. It can be inferred that as the orientation rotates through θ = 9.5 deg and θ = 21 deg toward the armchair direction, the bonds align more favorably with macroscopic loading, resulting in prolonged failure.

Besides crystal orientations, temperature is known to affect the mechanical behavior of TMDs. To understand the temperature dependency, we uniaxially stretch single-layer films at different temperatures. Six temperatures ranging from 1 K to 300 K are considered. In Fig. 6, the stress-strain responses of the armchair MoTe₂ monolayer at different temperatures are presented. With increasing temperature, a significant degree of deterioration in mechanical properties is observed. Fracture properties in terms of UTS, failure strain, and fracture toughness all display a decaying trend with temperature. In particular, failure strain drops from 34.7% to 19.5% when temperature is increased from 1 K to 300 K, underlining the adverse influence of temperature on mechanical properties. This behavior is expected, as there is an increased degree of disorder in the system at higher temperatures due to the rapid thermal fluctuation of atoms. This phenomenon is further compounded with applied deformation, weakening the covalent bonds and leading to premature bond breakage. Furthermore, the elastic properties (Young’s modulus) monotonically decrease with temperature which is consistent with theoretical predictions of graphene [58]. For other crystal orientations, the effect of temperature shows the same trend (Fig. S4). For practical applications, the interplay of temperature and orientation on the final mechanical properties highlighted here needs to be considered. Overall, our findings can be used as a guide to tune macroscopic parameters for achieving desired response when deploying the material.

Strain engineering of TMDs immediately arises as an opportunity for various technological applications as these materials exhibit impressive strain-tunable optical and electronic properties combined with high mechanical limits. As the community moves closer to using TMDs for such optoelectronic/electronic applications, there must be an understanding of how TMDs interact with thin films deposited during micro/nanofabrication. Thin films are crucial for functional uses in technological applications, ranging from protective optical coatings to contact electrodes. However, the development of a thin film’s microstructure during growth will inevitably create some amount of residual stress that transfers into whatever it is deposited onto. Without proper characterization, residual stress in a thin film can lead to failure modes such as delamination, or undesired alteration in device performance.

Alternatively, taking advantage of process induced stress of thin films, one could intentionally engineer residual stresses in thin films to manipulate strain in micro/nanofabricated devices to enhance device performance. Engineering thin film stress has been well-implemented in silicon-based transistors, where the first demonstration of depositing silicon nitride stressors allowed carrier mobilities to be enhanced while scaling down transistor dimensions [59]. This technique was quickly adopted in commercial silicon transistor manufacturing by 2004 and has been part of almost all electronics since that point. One can deposit tunable thin films with tensile or compressive stress and consequently allows for applied strain to be tailored locally to individual devices in a highly dense integrated circuit. From a broader perspective, in the pursuit of realizing TMD-based electronics or optoelectronics, TMDs will inevitably be engineered with thin films that contain process induced stress since any fabrication process will create residual stress, and it may be freely used to enhance device performance. For this purpose, significant effort still needs to be made in characterizing strain transfer in 2D-bonded TMDs and how it differs from the well understood 3D-bonded strain engineering techniques of the past.

As TMDs are van der Waals bonded in the out-of-plane direction, strain engineering multilayer TMD structures requires additional insight of strain transfer throughout the c-axis. When one layer is under an applied in-plane load, several works have observed incomplete transfer of shear traction between that layer and the subsequent weakly bonded adjacent layers. This incomplete shear traction has been quantified with an interfacial stiffness constant, unique to the interlayer adhesion of the given 2D van der Waals heterostructure [60]. The work of separation between 2L van der Waals systems varies with interlayer adhesion, 2L MoTe₂ requiring twice the work of separation than 2L MoS₂ [61]. Based off these previous findings, we suspect the strain transfer length scale of the c-axis in MoTe₂ to be twice as long as observed in MoS₂. We first seek to extract the strain transfer length scale of MoTe₂ in the c-axis via molecular statics (MS) simulations with our newly developed SW potential, to confirm a longer strain transfer length scale than that observed in MoS₂.

Our previous work has implemented the MS simulations for MoS₂ samples varying in thickness (2L-7L), where we only observe strain transferred into the top two layers before decreasing below 0.1% strain in the succeeding layers [46]. We next perform the same simulations with the described MoTe₂ SW potential, the strain transferred in MoS₂ and MoTe₂ is presented visually with contour plots of atomic strain magnitude (Figs. 7(a) and 7(b)). The simulations are conducted in MoTe₂ by applying 0.85% biaxial strain onto the top layer in samples varying in thickness (2L, 4L, 6L, 8L, and 10L samples); we note that 0.85% biaxial strain is still under the limit where interlayer slippage may occur. The bottom layer of these structures is set to be 75% fixed, which is computed using a weighted average of two separate sets of MS simulations with fixed and free bottom layer boundary conditions. In a 6L sample, strain through each subsequent layer after the top layer decreases by a factor of ∼1.75 from the top layer max value of 0.85% strain (Fig. 7(c)). We extract that strain penetrates into the top four layers before it is below the 0.1% level.

To further validate these findings, we attempt to extract this strain transfer length scale experimentally with Raman spectroscopy. We fully encapsulate various MoTe₂ flakes by evaporating a Al₂O₃ (10 nm)/MgF₂ (100 nm)/ Al₂O₃ (10 nm) trilayer thin film stressed in tension, therefore applying biaxial strain onto the top layer of the MoTe₂ samples. The Al₂O₃ layers are used to enhance adhesion to the 2D material and also to protect the film from strain relaxation from environmental interaction. The optical transparency of the thin film stressors allows for us to easily probe the strain transferred into the MoTe₂ samples via Raman scattering. The in-plane Raman mode (⁠$E2g1$⁠) has been theoretically and experimentally verified to shift in peak position with applied biaxial strain onto TMDs [62]. Therefore, we extract the peak shifts (⁠$ΔE2g1$⁠) by quantifying the differences in the $E2g1$ peak positions of stressed and control samples with respect to MoTe₂ flake thickness (2L-11L). To compare the strain transfer length scale, we present the measured $ΔE2g1$ versus thickness for MoS₂ and MoTe₂ (Fig. 7(d)). When fitting both $ΔE2g1$ versus thickness to that of an exponential decay function (⁠$ΔE2g1=C0exp−(x/λ)$⁠, where x = # of layers), we extract a ratio of $λMoTe2λMoS2∼1.78±0.39$⁠. This ratio confirms approximately double the strain transferred throughout the c-axis in MoTe₂ than MoS₂, as predicted from increased van der Waals coupling.

We next confirm what is being measured experimentally matches the strain predicted from our MS simulations. As there is heterostrain present within this system, the measured Raman signature is a superposition of responses of the given layers with varying strain. The final Raman signature being a superposition of Raman responses is most apparently observed in heterostructures [63]. This superposition of Raman responses manifests as the smooth exponential decay trend as seen in $ΔE2g1$ with thickness. Therefore, we calculate the suspected final Raman signature by creating a superposition of Lorentzian signatures from each layer within a given sample thickness, where each peak position is calculated to match the strain determined from our MS results. Other parameters such as intensities and full-width-half-maximums are extracted from our experimental Raman results. In order to convert the peak positions to represent the strain transferred, we utilize translation factors from other works that relate $ΔE2g1$ versus biaxial strain for MoS₂.

In-plane Raman mode (⁠$E2g1$⁠) shift per unit strain, referred to as “translation factors” for the rest of this text, have been observed to be dependent on both the material and its thickness. These translation values have been reported to decrease with increasing material thickness, suggesting the Raman modes become less responsive to strain with increasing thickness. While some studies suggested that a decrease in translation values originates from substrate adhesion issues, other works have debunked this possibility by extracting translation values of few-layer graphene samples optimized with an epoxy encapsulation [64]. This trend has repeatedly been verified to exist in MoS₂, confirming that these translation values are related to the intrinsic properties of the material and its given thickness with respect to the nature of the applied strain [65]. For our calculations, we utilize $ΔE2g1$ versus biaxial strain translation factors for MoS₂ of 5.2, 4.2, 3, and 2.5 (⁠$cm−1ϵ%$⁠) corresponding to 1L, 2L, 3L, and 20L sample thicknesses, respectively. We note that we utilize these specific translation values since they have been experimentally reported with biaxial strain [8,62]. The values for 4L–19L are interpolated using an exponential function.

Translation values have not been characterized as thoroughly in MoTe₂ as MoS₂, however, we suspect the values to be similar as these materials are similar in composition. The translation values of MoTe₂ being close to those of MoS₂ is supplemented by the fact we observe similar $ΔE2g1$ values in both 2L TMD samples, when applying the exact same thin film stressor (same magnitude of strain applied to the top TMD layer). Theoretical work has also observed ∼5.8 $cm−1ϵ%$ slope with biaxial strain onto 1L-MoTe₂; therefore, we choose to approximate the MoTe₂ translation values to be uniformly 1.12 times larger than those of MoS₂ to match the 5.8 $cm−1ϵ%$ 1L result (i.e., 5.8, 4.7, and 3.36 $cm−1ϵ%$ for 1L-3L, respectively) [66].

After taking all factors into account and calculating $ΔE2g1$ with the peak positions from MS simulation, we find the calculated exponential decay matches well with our experimental results for both MoS₂ and MoTe₂ (dashed lines presented in Fig. 7(d)). This confirms that our result is self-consistent and represents an experimental verification of the MS simulations as well as confirmation that the SW potential is valid. Since this experimental confirmation is not unique to MoTe₂, it can be used as a method to quantify out-of-plane strain transfer for various 2D materials.

The computational approach should be generally applicable to other types of interface rather than the equilibrium commensurate stacking considered here. Several models have been recently developed that describe interlayer interactions such as the Frenkel–Kontorova model [67], shear-lag model [68], and registry-dependent Kolmogorov–Crespi potential [69]. These models consider stacking order as well the energy barrier for relative shear between layers at the interface and can potentially be more robust for studying systems such as heterostructures, twisted bilayers, and the effect of substrate materials. Although we used a simple LJ potential for the interlayer van der Waals energy, it still can predict complex nonlinear phenomena. An example is the formation of strain solitons in the case of the bilayer structure when strain exceeds the slippage strain (calculated to be 1.95%) in the top layer (Fig. S5). The presence of such strain solitons has been reported from both computational works [60,67] using the above mentioned models and from experimental[70,71] works, similar to our results.

In this paper, we have presented a study of the deformation behavior of MoTe₂ TMD thin films, an important class of engineering material. We have implemented a computational framework to explore the mechanical response of monolayer and multilayer MoTe₂ at the nanoscale by atomistic simulations. We first developed an accurate empirical potential in the form of a Stillinger–Weber force field, which successfully predicted elastic and failure properties as compared to first-principles methods and experiments, respectively. Upon using the potential to study the mechanical behavior in MD simulations, our major findings can be summarized as follows:

• Uniaxial tensile loading reveals directional anisotropy in the nonlinear mechanical response.

• Macroscopic factors such as crystal orientation and temperature influence the elastic and failure properties.

• Out-of-plane strain transfer in multilayered films continues predominantly up to the top four layers, with the resultant calculated Raman peak shifts matching with experiments.

The interdependency of different macroscopic factors and their contributions to the mechanical response shown here, along with the calculation of strain transfer length scale, provide vital details needed for strain engineering this material. Our results will facilitate the development of future models for the study of more complex MoTe₂ structures such as heterostructures and Moiré patterns in twisted bilayers.

We wish to acknowledge support from the National Science Foundation (OMA-1936250) and National Science Foundation Graduate Research Fellowship Program (DGE-1939268). We performed Raman spectroscopy at the Cornell Center for Materials Research Shared Facilities, who are funded through NSF MRSEC programme (DMR1719875). KI and SMG (first principles calculations) were supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator (QSA). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant No. ACI-1548562. Work at the Molecular Foundry was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Hesam Askari remembers Prof. Hussein M. Zbib for his invaluable inspiration, friendly mentorship, and casual conversations. He is sorely missed.

There are no conflicts of interest.