## Abstract

Transition metal dichalcogenides (TMDs) offer superior properties over conventional materials in many areas such as in electronic devices. In recent years, TMDs have been shown to display a phase switching mechanism under the application of external mechanical strain, making them exciting candidates for phase change transistors. Molybdenum ditelluride (MoTe2) is one such material that has been engineered as a strain-based phase change transistor. In this work, we explore various aspects of the mechanical properties of this material by a suite of computational and experimental approaches. First, we present parameterization of an interatomic potential for modeling monolayer as well as multilayered MoTe2 films. For generating the empirical potential parameter set, we fit results from density functional theory calculations using a random search algorithm known as particle swarm optimization. The potential closely predicts structural properties, elastic constants, and vibrational frequencies of MoTe2 indicating a reliable fit. Our simulated mechanical response matches earlier larger scale experimental nanoindentation results with excellent prediction of fracture points. Simulation of uniaxial tensile deformation by molecular dynamics shows the complete non-linear stress-strain response up to failure. Mechanical behavior, including failure properties, exhibits directional anisotropy due to the variation of bond alignments with crystal orientation. Furthermore, we show the deterioration of mechanical properties with increasing temperature. Finally, we present computational and experimental evidence of an extended c-axis strain transfer length in MoTe2 compared to TMDs with smaller chalcogen atoms.

## 1 Introduction

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 (MX2 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 (MoS2, WS2) have well-known tribological applications such as lubrication and wear resistance [5,6] in addition to their electrical properties. Molybdenum disulfide MoS2 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, MoTe2, which has an indirect bulk bandgap close to Si (∼1 eV), making it promising for on-chip integration [11]. In recent years, MoTe2 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 MoS2, MoTe2 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 MoTe2. Many techniques have been used to induce phase transitions in MoTe2 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, MoTe2 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 MoTe2 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 MoTe2 thin films. We first develop parameterization of an SW interatomic potential that can accurately replicate mechanical properties of MoTe2 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 MoTe2 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 MoTe2 multilayered films stressed by a capping layer. Our findings help in further understanding the strain-induced properties of MoTe2 thin films in strain-engineered electrical devices.

## 2 Methodology

### 2.1 Structural Description.

MoTe2 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 MoTe2 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.

Fig. 1
Fig. 1
Close modal

### 2.2 Computational Details.

Density functional theory (DFT) calculations were carried out using the Vienna ab initio Simulation Package (vasp)[3538] 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−8 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 Å−1. 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 Å−1.

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−14). A uniaxial tensile deformation was applied to monolayer MoTe2 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 (XY) 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 109 s−1 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−16) 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 MoS2 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 MoS2 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 MoTe2 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 MoTe2 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.

### 2.3 Experimental Characterization.

MoTe2 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 H2O and O2). The samples were then placed inside a vacuum chamber for thin film stressor deposition. Thin film stressors were e-beam evaporated at 5 × 10−6 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−1. Finally, the peak positions of the given Raman spectra were determined by fitting to a Voigt spectral profile.

### 2.4 Interatomic Potential.

The covalent bonding of intralayer atoms can be described by a SW potential that includes two-body bond-stretching and three-body angle-bending terms in the following form [30]:
$E=∑i∑j>iϕ2(rij)+∑i∑j≠i∑k>jϕ3(rij,rik,θijk)$
(1)
Here, ϕ2 and ϕ3 denote two-body and three-body interactions, respectively, in the following form:
$ϕ2=Aij(Bijrij4−1)exp[ρijrij−r12max]$
(2)
and
$ϕ3=Kijkexp[ρijrij−r12max+ρikrik−r13max](cosθ−cosθ0)2$
(3)
where, rij is the instantaneous bond length between atoms i and j, and θ is the angle formed by two bonds with θ0 being their equilibrium angle. Cutoff parameters (rmax) control the distance up to which the interactions between neighboring atoms are included in the potential. The SW potential has been previously used to successfully capture nonlinear mechanical behavior of TMD materials [28,49], and it is already implemented in many molecular dynamics programs such as GULP [50] and LAMMPS. For empirical parameters fitting, GULP is particularly useful as it offers lattice dynamics calculations. It should be noted that to capture only the correct angle-bending terms, GULP uses an additional cutoff parameter $r23max$. In this way, large Mo–Te–Te angles such as θu shown in Fig. 1(a) are excluded from calculations as the Te–Te distance in these angles are beyond the cutoff value. However, LAMMPS does not contain this additional three-body cutoff parameter in its formulations. As a solution, we modify the SW potential source code in LAMMPS to avoid large non-physical angles. An additional angle cutoff condition is implemented that exclude angles outside a specific range (|cos θ − cos θ0| > 0.35). Finally, the weak interlayer interactions are modeled using a 12-6 Lennard–Jones potential (LJ) which has the following equation for energy:
$ELJ=4ϵ[(σrij)12−(σrij)6],rij≤rcut$
(4)
where rcut, the cutoff distance for long-range interactions, was chosen to be 10 Å.

### 2.5 Parameterization of the Potential.

Global optimization methods based on random search techniques offer several advantages for generating empirical potential parameters fitted to physical properties. Starting with a large domain of random parameter sets, they allow populating a wider region of the parameter search space. Therefore, they are more likely to avoid pitfalls such as local optima in the potential energy surface. Additionally, as there is no need for carefully estimating the initial guess, they are less sensitive to initial conditions. Within the global optimization framework, particle swarm optimization (PSO) [51,52] is a popular method due to its robustness, simplicity in implementation, and fast convergence. We use PSO to obtain the best parameter set of the SW potential that minimizes a single global cost function in the form of error measured with respect to physical quantities obtained by DFT. In the context of PSO, each initial set of parameters is called a particle ($x→$) in a swarm occupying the parameter space and is allowed to communicate with other particles in the swarm [53]. Starting from a randomly assigned position (in the parameter search space), every particle evaluates a cost function (fitness value) which is defined as the measurement error. At the end of each iteration, every particle updates its velocity and position as a combination of three vector quantities: (1) an inertia term consisting of the particle’s current velocity, (2) a cognitive term connecting current position and the position which resulted in best fitness value ($Pi→$) for the specific particle up to the current iteration, and finally, (3) a social component measuring the relative position from the best candidate particle (best fitness value) in the swarm ($g→$). It is important to emphasize that position and velocity within the context of PSO only refers to an arbitrary set of parameters and a perturbation to change them iteratively and have no association with physical positions or velocity of atoms. The following sets of equations are then used to update the position and velocity of individual particles.
$vi→(t+1)=wvi→(t)+r1c1[Pi→(t)−xi→(t)]+r2c2[g→(t)−xi→(t)]$
(5)
$xi→(t+1)=xi→(t)+vi→(t+1)$
(6)
Here, w, c1, and c2 are constants designed to avoid instability in velocity evolution, and r1 and r2 are random numbers drawn from [0, 1]. The iteration process continues until the cost function is minimized and the best parameter set is found. For the cost function, we use a weighted sum of the squared difference between reference quantities and their evaluated values using a specific parameter set. The cost function can be expressed as follows:
$cost(x→)=∑i=1Nwi[firef−fieval(x→)]2$
(7)
Here, $x→$ denotes a set of force-field parameters, N runs over the number of reference observables, and the goal of the algorithm is to find the best candidate that minimizes the defined sum of squared error cost function. For a monolayer MoTe2 system, the required parameter set is 13-dimensional considering the two-body and three-body interactions between different atoms in the SW potential. We use the lattice dynamics program GULP to evaluate the required properties using a generated SW potential parameter set in each PSO iteration. The reference dataset for cost function calculations included structural parameters in lattice constants, mechanical properties in elastic constants, and vibrational properties in the form of phonon frequencies (Γ − M direction), all computed by DFT. SW parameters that minimized the cost function do not rely on initial guess and are presented in Tables 13 in both GULP and LAMMPS formats. Only directly bonded interactions are included with non-bonded interactions set to zero. Mo–Mo–Mo and Te–Te–Te terms only have two-body contributions. It is to be noted that for the SW potential, quantities such as equilibrium angle θ0 are taken directly from DFT calculations and cutoff distances are set to ∼1.3 − 1.4 times the equilibrium bond lengths. For the LJ potential, the distance parameter σ is derived from the equilibrium interlayer atomic distance and the energy parameter is fitted to the out-of-plane elastic constant and phonon frequencies of bulk MoTe2. LJ parameters used in this study are shown in Table 4 for multilayered structures.
Table 1

Fitted two-body parameters of SW potential in GULP format

ParametersAρBrminrmax
Mo–Te3.10261.481637.95590.03.55
Te–Te0.32410.905119.81590.04.9
Mo–Mo2.17060.567520.52410.04.9
ParametersAρBrminrmax
Mo–Te3.10261.481637.95590.03.55
Te–Te0.32410.905119.81590.04.9
Mo–Mo2.17060.567520.52410.04.9
Table 2

Fitted three-body parameters of SW potential in GULP format

ParametersKθ0ρ1ρ2$rmax12$$rmax13$$rmax23$
Mo–Te–Te26.270680.5540.82160.82163.553.554.9
Te–Mo–Mo12.152680.5540.860560.860563.553.554.9
ParametersKθ0ρ1ρ2$rmax12$$rmax13$$rmax23$
Mo–Te–Te26.270680.5540.82160.82163.553.554.9
Te–Mo–Mo12.152680.5540.860560.860563.553.554.9
Table 3

Fitted SW potential parameters in LAMMPS format

Parametersε (eV)σ (Å)aλγcos θ0ABpqtol
Mo–Te–Te1.0001.48162.395926.27060.55450.16413.10267.87594.0000.0000.0
Te–Mo–Mo1.0001.48162.395912.15260.58080.16413.10267.87594.0000.0000.0
Te–Te–Te1.0000.90515.41350.00.00.00.324129.52264.0000.0000.0
Mo–Mo–Mo1.0000.56758.63360.00.00.02.1706197.81214.0000.0000.0
Parametersε (eV)σ (Å)aλγcos θ0ABpqtol
Mo–Te–Te1.0001.48162.395926.27060.55450.16413.10267.87594.0000.0000.0
Te–Mo–Mo1.0001.48162.395912.15260.58080.16413.10267.87594.0000.0000.0
Te–Te–Te1.0000.90515.41350.00.00.00.324129.52264.0000.0000.0
Mo–Mo–Mo1.0000.56758.63360.00.00.02.1706197.81214.0000.0000.0

Note: Units are in LAMMPS metal style.

Table 4

LJ potential parameters used for modeling interlayer interactions in MoTe2

Parametersσ (Å)ε (eV)
Te–Te3.5380.0304
Parametersσ (Å)ε (eV)
Te–Te3.5380.0304

### 2.6 Validation of the Fitted Parameters.

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 MoTe2. 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.

Table 5

Comparison of elastic constants Cij and Young’s modulus E using the fitted SW potential in LAMMPS with DFT calculations and experiment

MDDFT/experimentMagnitude difference (%)
C11127.8125.91.5
C22127.7125.91.5
C1230.731.52.5
C3345.545.70.4
C6649.147.24
E110112.5 [54]2.2
MDDFT/experimentMagnitude difference (%)
C11127.8125.91.5
C22127.7125.91.5
C1230.731.52.5
C3345.545.70.4
C6649.147.24
E110112.5 [54]2.2

Note: Units are in GPa.

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 MoTe2 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.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal

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 MoTe2 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 (σ0) 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 MoTe2, which is the aim for the remaining parts of the paper.

Fig. 4
Fig. 4
Close modal

## 3 Results and Discussions

### 3.1 Uniaxial Deformation.

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

#### 3.1.1 Orientation Dependence.

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 MoTe2 in armchair and zigzag loading [55,56].

Fig. 5
Fig. 5
Close modal

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.

#### 3.1.2 Temperature Effect.

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 MoTe2 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.

Fig. 6
Fig. 6
Close modal

### 3.2 Static Strain Transfer.

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 MoTe2 requiring twice the work of separation than 2L MoS2 [61]. Based off these previous findings, we suspect the strain transfer length scale of the c-axis in MoTe2 to be twice as long as observed in MoS2. We first seek to extract the strain transfer length scale of MoTe2 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 MoS2.

Our previous work has implemented the MS simulations for MoS2 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 MoTe2 SW potential, the strain transferred in MoS2 and MoTe2 is presented visually with contour plots of atomic strain magnitude (Figs. 7(a) and 7(b)). The simulations are conducted in MoTe2 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.

Fig. 7
Fig. 7
Close modal

To further validate these findings, we attempt to extract this strain transfer length scale experimentally with Raman spectroscopy. We fully encapsulate various MoTe2 flakes by evaporating a Al2O3 (10 nm)/MgF2 (100 nm)/ Al2O3 (10 nm) trilayer thin film stressed in tension, therefore applying biaxial strain onto the top layer of the MoTe2 samples. The Al2O3 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 MoTe2 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 MoTe2 flake thickness (2L-11L). To compare the strain transfer length scale, we present the measured $ΔE2g1$ versus thickness for MoS2 and MoTe2 (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 MoTe2 than MoS2, 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 MoS2.

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 MoS2, 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 MoS2 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 MoTe2 as MoS2, however, we suspect the values to be similar as these materials are similar in composition. The translation values of MoTe2 being close to those of MoS2 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-MoTe2; therefore, we choose to approximate the MoTe2 translation values to be uniformly 1.12 times larger than those of MoS2 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 MoS2 and MoTe2 (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 MoTe2, 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.

## 4 Conclusions

In this paper, we have presented a study of the deformation behavior of MoTe2 TMD thin films, an important class of engineering material. We have implemented a computational framework to explore the mechanical response of monolayer and multilayer MoTe2 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 MoTe2 structures such as heterostructures and Moiré patterns in twisted bilayers.

## Acknowledgment

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.

## Conflict of Interest

There are no conflicts of interest.

## References

1.
Novoselov
,
K. S.
,
Geim
,
A. K.
,
Morozov
,
S. V.
,
Jiang
,
D.
,
Zhang
,
Y.
,
Dubonos
,
S. V.
,
Grigorieva
,
I. V.
, and
Firsov
,
A. A.
,
2004
, “
Electric Field Effect in Atomically Thin Carbon Films
,”
Science
,
306
(
5696
), pp.
666
669
.
2.
Stoller
,
M. D.
,
Park
,
S.
,
Zhu
,
Y.
,
An
,
J.
, and
Ruoff
,
R. S.
,
2008
, “
Graphene-Based Ultracapacitors
,”
Nano. Lett.
,
8
(
10
), pp.
3498
3502
.
3.
Wang
,
Q. H.
,
Kalantar-Zadeh
,
K.
,
Kis
,
A.
,
Coleman
,
J. N.
, and
Strano
,
M. S.
,
2012
, “
Electronics and Optoelectronics of Two-Dimensional Transition Metal Dichalcogenides
,”
Nat. Nanotechnol.
,
7
(
11
), pp.
699
712
.
4.
Yang
,
H.
,
Kim
,
S. W.
,
Chhowalla
,
M.
, and
Lee
,
Y. H.
,
2017
, “
Structural and Quantum-State Phase Transitions in Van Der Waals Layered Materials
,”
Nat. Phys.
,
13
(
10
), pp.
931
937
.
5.
Savan
,
A.
,
Pflüger
,
E.
,
Voumard
,
P.
,
Schröer
,
A.
, and
Simmonds
,
M.
,
2000
, “
Modern Solid Lubrication: Recent Developments and Applications of MoS2
,”
Lubrication Sci.
,
12
(
2
), pp.
185
203
.
6.
Rapoport
,
L.
,
Leshchinsky
,
V.
,
Lvovsky
,
M.
,
Nepomnyashchy
,
O.
,
Volovik
,
Y.
, and
Tenne
,
R.
,
2002
, “
Friction and Wear of Powdered Composites Impregnated With WS2 Inorganic Fullerene-like Nanoparticles
,”
Wear
,
252
(
5-6
), pp.
518
527
.
7.
Bertolazzi
,
S.
,
Brivio
,
J.
, and
Kis
,
A.
,
2011
, “
Stretching and Breaking of Ultrathin MoS2
,”
ACS. Nano.
,
5
(
12
), pp.
9703
9709
.
8.
Nayak
,
A. P.
,
Bhattacharyya
,
S.
,
Zhu
,
J.
,
Liu
,
J.
,
Wu
,
X.
,
Pandey
,
T.
,
Jin
,
C.
,
Singh
,
A. K.
,
Akinwande
,
D.
, and
Lin
,
J. -F.
,
2014
, “
Pressure-Induced Semiconducting to Metallic Transition in Multilayered Molybdenum Disulphide
,”
Nat. Commun.
,
5
(
1
), pp.
1
9
.
9.
Desai
,
S. B.
,
Seol
,
G.
,
Kang
,
J. S.
,
Fang
,
H.
,
Battaglia
,
C.
,
Kapadia
,
R.
,
Ager
,
J. W.
,
Guo
,
J.
, and
Javey
,
A.
,
2014
, “
Strain-Induced Indirect to Direct Bandgap Transition in Multilayer WSe2
,”
Nano. Lett.
,
14
(
8
), pp.
4592
4597
.
10.
Tan
,
Y.
,
Luo
,
F.
,
Zhu
,
M.
,
Xu
,
X.
,
Ye
,
Y.
,
Li
,
B.
,
Wang
,
G.
,
Luo
,
W.
,
Zheng
,
X.
,
Wu
,
N.
,
Yu
,
Y.
,
Qin
,
S.
, and
Zhang
,
X.-A.
,
2018
, “
Controllable 2H-to-1T Phase Transition in Few-Layer MoTe2
,”
Nanoscale
,
10
(
42
), pp.
19964
19971
.
11.
Lezama
,
I. G.
,
Ubaldini
,
A.
,
Longobardi
,
M.
,
Giannini
,
E.
,
Renner
,
C.
,
Kuzmenko
,
A. B.
, and
Morpurgo
,
A. F.
,
2014
, “
Surface Transport and Band Gap Structure of Exfoliated 2H-MoTe2 Crystals
,”
2D. Mater.
,
1
(
2
), p.
021002
.
12.
Bie
,
Y.-Q.
,
Grosso
,
G.
,
Heuck
,
M.
,
Furchi
,
M. M.
,
Cao
,
Y.
,
Zheng
,
J.
,
Bunandar
,
D.
,
Navarro-Moratalla
,
E.
,
Zhou
,
L.
,
Efetov
,
D. K.
,
Taniguchi
,
T.
,
Watanabe
,
K.
,
Kong
,
J.
,
Englund
,
D.
, and
Jarillo-Herrero
,
P.
,
2017
, “
A MoTe 2-Based Light-Emitting Diode and Photodetector for Silicon Photonic Integrated Circuits
,”
Nat. Nanotechnol.
,
12
(
12
), pp.
1124
1129
.
13.
Panda
,
M. R.
,
Gangwar
,
R.
,
Muthuraj
,
D.
,
Sau
,
S.
,
Pandey
,
D.
,
Banerjee
,
A.
,
Chakrabarti
,
A.
,
Sagdeo
,
A.
,
Weyland
,
M.
,
Majumder
,
M.
,
Bao
,
Q.
, and
Mitra
,
S.
,
2020
, “
High Performance Lithium-Ion Batteries Using Layered 2H-MoTe2 As Anode
,”
Small
,
16
(
38
), p.
2002669
.
14.
Duerloo
,
K.-A. N.
,
Ong
,
M. T.
, and
Reed
,
E. J.
,
2012
, “
Intrinsic Piezoelectricity in Two-Dimensional Materials
,”
J. Phys. Chem. Lett.
,
3
(
19
), pp.
2871
2876
.
15.
Ruppert
,
C.
,
Aslan
,
O. B.
, and
Heinz
,
T. F.
,
2014
, “
Optical Properties and Band Gap of Single-and Few-Layer MoTe2 Crystals
,”
Nano. Lett.
,
14
(
11
), pp.
6231
6236
.
16.
Duerloo
,
K.-A. N.
,
Li
,
Y.
, and
Reed
,
E. J.
,
2014
, “
Structural Phase Transitions in Two-Dimensional Mo-and W-Dichalcogenide Monolayers
,”
Nat. Commun.
,
5
(
1
), pp.
1
9
.
17.
Cho
,
S.
,
Kim
,
S.
,
Kim
,
J. H.
,
Zhao
,
J.
,
Seok
,
J.
,
Keum
,
D. H.
,
Baik
,
J.
,
Choe
,
D.-H.
,
Chang
,
K. J.
,
Suenaga
,
K.
,
Kim
,
S. W.
,
Lee
,
Y. H.
, and
Yang
,
H.
,
2015
, “
Phase Patterning for Ohmic Homojunction Contact in MoTe2
,”
Science
,
349
(
6248
), pp.
625
628
.
18.
Nan
,
H.
,
Jiang
,
J.
,
Xiao
,
S.
,
Chen
,
Z.
,
Luo
,
Z.
,
Zhang
,
L.
,
Zhang
,
X.
,
Qi
,
H.
,
Gu
,
X.
,
Wang
,
X.
, and
Ni
,
Z.
,
2018
, “
Soft Hydrogen Plasma Induced Phase Transition in Monolayer and Few-Layer MoTe2
,”
Nanotechnology
,
30
(
3
), p.
034004
.
19.
Wang
,
Y.
,
Xiao
,
J.
,
Zhu
,
H.
,
Li
,
Y.
,
Alsaid
,
Y.
,
Fong
,
K. Y.
,
Zhou
,
Y.
,
Wang
,
S.
,
Shi
,
W.
,
Wang
,
Y.
,
Zettl
,
A.
,
Reed
,
E. J.
, and
Zhang
,
X.
,
2017
, “
Structural Phase Transition in Monolayer MoTe2 Driven by Electrostatic Doping
,”
Nature
,
550
(
7677
), pp.
487
491
.
20.
Rhodes
,
D.
,
Chenet
,
D. A.
,
Janicek
,
B. E.
,
Nyby
,
C.
,
Lin
,
Y.
,
Jin
,
W.
,
Edelberg
,
D.
,
Mannebach
,
E.
,
Finney
,
N.
,
Antony
,
A.
,
Schiros
,
T.
,
Klarr
,
T.
,
Mazzoni
,
A.
,
Chin
,
M.
,
Chiu
,
Y.-c.
,
Zheng
,
W.
,
Zhang
,
Q. R.
,
Ernst
,
F.
,
Dadap
,
J. I.
,
Tong
,
X.
,
Ma
,
J.
,
Lou
,
R.
,
Wang
,
S.
,
Qian
,
T.
,
Ding
,
H.
,
Osgood
,
R. M.
,
Paley
,
D. W.
,
Lindenberg
,
A. M.
,
Huang
,
P. Y.
,
Pasupathy
,
A. N.
,
Dubey
,
M.
,
Hone
,
J.
, and
Balicas
,
L.
,
2017
, “
Engineering the Structural and Electronic Phases of MoTe2 Through W Substitution
,”
Nano. Lett.
,
17
(
3
), pp.
1616
1622
.
21.
Ghasemi
,
A.
, and
Gao
,
W.
,
2020
, “
Atomistic Mechanism of Stress Modulated Phase Transition in Monolayer MoTe2
,”
Extreme Mechanics Letters
,
40
(Special Issue in Honor of Horacio D. Espinosa, recipient of the 2019 Prager Medal), p.
100946
.
22.
Song
,
S.
,
Keum
,
D. H.
,
Cho
,
S.
,
Perello
,
D.
,
Kim
,
Y.
, and
Lee
,
Y. H.
,
2016
, “
Room Temperature Semiconductor–Metal Transition of MoTe2 Thin Films Engineered by Strain
,”
Nano. Lett.
,
16
(
1
), pp.
188
193
.
23.
Diaz
,
H. C.
,
Chaghi
,
R.
,
Ma
,
Y.
, and
Batzill
,
M.
,
2015
, “
Molecular Beam Epitaxy of the Van Der Waals Heterostructure MoTe2 on MoS2: Phase, Thermal, and Chemical Stability
,”
2D. Mater.
,
2
(
4
), p.
044010
.
24.
Liu
,
K.
,
Yan
,
Q.
,
Chen
,
M.
,
Fan
,
W.
,
Sun
,
Y.
,
Suh
,
J.
,
Fu
,
D.
,
Lee
,
S.
,
Zhou
,
J.
,
Tongay
,
S.
,
Ji
,
J.
,
Neaton
,
J. B.
, and
Wu
,
J.
,
2014
, “
Elastic Properties of Chemical-Vapor-Deposited Monolayer MoS2, WS2, and Their Bilayer Heterostructures
,”
Nano. Lett.
,
14
(
9
), pp.
5097
5103
.
25.
Apte
,
A.
,
Kochat
,
V.
,
Rajak
,
P.
,
Krishnamoorthy
,
A.
,
Manimunda
,
P.
,
Hachtel
,
J. A.
,
Idrobo
,
J. C.
,
Syed Amanulla
,
S. A.
,
Vashishta
,
P.
,
Nakano
,
A.
,
Kalia
,
R. K.
,
Tiwary
,
C. S.
, and
Ajayan
,
P. M.
,
2018
, “
Structural Phase Transformation in Strained Monolayer MoWSe2 Alloy
,”
ACS. Nano.
,
12
(
4
), pp.
3468
3476
.
26.
Hou
,
W.
,
Azizimanesh
,
A.
,
Sewaket
,
A.
,
Peña
,
T.
,
Watson
,
C.
,
Liu
,
M.
,
Askari
,
H.
, and
Wu
,
S. M.
,
2019
, “
Strain-Based Room-Temperature Non-Volatile MoTe 2 Ferroelectric Phase Change Transistor
,”
Nat. Nanotechnol.
,
14
(
7
), pp.
668
673
.
27.
Shao
,
S.
,
Zbib
,
H.
,
Mastorakos
,
I.
, and
Bahr
,
D.
,
2013
, “
Effect of Interfaces in the Work Hardening of Nanoscale Multilayer Metallic Composites During Nanoindentation: A Molecular Dynamics Investigation
,”
ASME J. Eng. Mater. Technol.
,
135
(
2
), p.
021001
.
28.
Jiang
,
J.-W.
, and
Park
,
H. S.
,
2014
, “
Mechanical Properties of MoS2/graphene Heterostructures
,”
Appl. Phys. Lett.
,
105
(
3
), p.
033108
.
29.
Ma
,
F.
,
Sun
,
Y.
,
Ma
,
D.
,
Xu
,
K.
, and
Chu
,
P. K.
,
2011
, “
Reversible Phase Transformation in Graphene Nano-Ribbons: Lattice Shearing Based Mechanism
,”
Acta. Mater.
,
59
(
17
), pp.
6783
6789
.
30.
Stillinger
,
F. H.
, and
Weber
,
T. A.
,
1985
, “
Computer Simulation of Local Order in Condensed Phases of Silicon
,”
Phys. Rev. B
,
31
(
8
), p.
5262
.
31.
Liang
,
T.
,
Phillpot
,
S. R.
, and
Sinnott
,
S. B.
,
2009
, “
Parametrization of a Reactive Many-Body Potential for Mo–S Systems
,”
Phys. Rev. B
,
79
(
24
), p.
245110
.
32.
Tersoff
,
J.
,
1989
, “
Modeling Solid-state Chemistry: Interatomic Potentials for Multicomponent Systems
,”
Phys. Rev. B
,
39
(
8
), p.
5566
.
33.
Hossain
,
M.
,
Hao
,
T.
, and
Silverman
,
B.
,
2018
, “
Stillinger–Weber Potential for Elastic and Fracture Properties in Graphene and Carbon Nanotubes
,”
J. Phys.: Condens. Matter.
,
30
(
5
), p.
055901
.
34.
Dawson
,
W.
, and
Bullett
,
D.
,
1987
, “
Electronic Structure and Crystallography of MoTe2 and WTe2
,”
J. Phys. C: Solid State Phys.
,
20
(
36
), p.
6159
.
35.
Kresse
,
G.
, and
Hafner
,
J.
,
1993
, “
Ab Initio Molecular Dynamics for Liquid Metals
,”
Phys. Rev. B
,
47
(
1
), pp.
558
561
.
36.
Kresse
,
G.
, and
Hafner
,
J.
,
1994
, “
Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal–Amorphous-Semiconductor Transition in Germanium
,”
Phys. Rev. B
,
49
(
20
), pp.
14251
14269
.
37.
Kresse
,
G.
, and
Furthmüller
,
J.
,
1996
, “
Efficiency of Ab Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set
,”
Comput. Mater. Sci.
,
6
, pp.
15
50
.
38.
Kresse
,
G.
, and
Furthmüller
,
J.
,
1996
, “
Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set
,”
Phys. Rev. B
,
54
, pp.
11169
11186
.
39.
Blöchl
,
P.
,
1994
, “
Projector Augmented-Wave Method
,”
Phys. Rev. B
,
50
(
24
), pp.
17953
17979
.
40.
Kresse
,
G.
, and
Joubert
,
D.
,
1999
, “
From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method
,”
Phys. Rev. B
,
59
, pp.
1758
1775
.
41.
Klime
,
J.
,
Bowler
,
D. R.
, and
Michaelides
,
A.
,
2011
, “
Van Der Waals Density Functionals Applied to Solids
,”
Phys. Rev. B - Condens. Matter Mater. Phys.
,
83
(
19
), pp.
1
13
.
42.
Togo
,
A.
, and
Tanaka
,
I.
,
2015
, “
First Principles Phonon Calculations in Materials Science
,”
Scr. Mater.
,
108
, pp.
1
5
.
43.
Plimpton
,
S.
,
1995
, “
Fast Parallel Algorithms for Short-Range Molecular Dynamics
,”
J. Comput. Phys.
,
117
(
1
), pp.
1
19
.
44.
Hoover
,
W. G.
,
1985
, “
Canonical Dynamics: Equilibrium Phase-Space Distributions
,”
Phys. Rev. A.
,
31
(
3
), p.
1695
.
45.
Stewart
,
J. A.
, and
Spearot
,
D.
,
2013
, “
Atomistic Simulations of Nanoindentation on the Basal Plane of Crystalline Molybdenum Disulfide (MoS2)
,”
Modell. Simul. Mater. Sci. Eng.
,
21
(
4
), p.
045003
.
46.
Peña
,
T.
,
Chowdhury
,
S. A.
,
Azizimanesh
,
A.
,
Sewaket
,
A.
,
Askari
,
H.
, and
Wu
,
S. M.
,
2020
, “
Strain Engineering 2D MoS2 With Thin Film Stress Capping Layers
”. arXiv preprint arXiv:2009.10626.
47.
Stukowski
,
A.
,
2009
, “
Visualization and Analysis of Atomistic Simulation Data With Ovito–the Open Visualization Tool
,”
Model. Simul. Mater. Sci. Eng
,
18
, p.
015012
.
48.
Shargh
,
A. K.
, and
Abdolrahim
,
N.
,
2019
, “
Molecular Dynamics Simulation of Structural Changes in Single Crystalline Silicon Nitride Nanomembrane
,”
Ceram. Int.
,
45
, pp.
23070
23077
.
49.
Xiong
,
S.
, and
Cao
,
G.
,
2015
, “
Molecular Dynamics Simulations of Mechanical Properties of Monolayer MoS2
,”
Nanotechnology
,
26
(
18
), p.
185705
.
50.
Gale
,
J. D.
,
1997
, “
GULP: A Computer Program for the Symmetry-Adapted Simulation of Solids
,”
J. Chem. Soc., Faraday. Trans.
,
93
(
4
), pp.
629
637
.
51.
Kandemir
,
A.
,
Yapicioglu
,
H.
,
Kinaci
,
A.
,
Çağın
,
T.
, and
Sevik
,
C.
,
2016
, “
Thermal Transport Properties of MoS2 and MoSe2 Monolayers
,”
Nanotechnology
,
27
(
5
), p.
055703
.
52.
Mobaraki
,
A.
,
Kandemir
,
A.
,
Yapicioglu
,
H.
,
Gülseren
,
O.
, and
Sevik
,
C.
,
2018
, “
Validation of Inter-Atomic Potential for WS2 and WSe2 Crystals Through Assessment of Thermal Transport Properties
,”
Comput. Mater. Sci.
,
144
, pp.
92
98
.
53.
Kennedy
,
J.
, and
Eberhart
,
R.
,
1995
, “
Particle Swarm Optimization
,”
Proceedings of ICNN’95-International Conference on Neural Networks
, Vol.
4
,
IEEE
, pp.
1942
1948
.
54.
Sun
,
Y.
,
Pan
,
J.
,
Zhang
,
Z.
,
Zhang
,
K.
,
Liang
,
J.
,
Wang
,
W.
,
Yuan
,
Z.
,
Hao
,
Y.
,
Wang
,
B.
,
Wang
,
J.
,
Wu
,
Y.
,
Zheng
,
J.
,
Jiao
,
L.
,
Zhou
,
S.
,
Liu
,
K.
,
Cheng
,
C.
,
Duan
,
W.
,
Xu
,
Y.
,
Yan
,
Q.
, and
Liu
,
K.
,
2019
, “
Elastic Properties and Fracture Behaviors of Biaxially Deformed, Polymorphic MoTe2
,”
Nano. Lett.
,
19
(
2
), pp.
761
769
.
55.
Li
,
J.
,
Medhekar
,
N. V.
, and
Shenoy
,
V. B.
,
2013
, “
Bonding Charge Density and Ultimate Strength of Monolayer Transition Metal Dichalcogenides
,”
J. Phys. Chem. C
,
117
(
30
), pp.
15842
15848
.
56.
Mortazavi
,
B.
,
Berdiyorov
,
G. R.
,
Makaremi
,
M.
, and
Rabczuk
,
T.
,
2018
, “
Mechanical Responses of Two-Dimensional MoTe2; Pristine 2H, 1T and 1T and 1T/2H Heterostructure
,”
Extreme Mech. Lett.
,
20
, pp.
65
72
.
57.
Xiong
,
Q.-l.
,
Kitamura
,
T.
, and
Li
,
Z.-h.
,
2017
, “
Crystal Orientation-Dependent Mechanical Property and Structural Phase Transition of Monolayer Molybdenum Disulfide
,”
J. Appl. Phys.
,
122
(
13
), p.
135105
.
58.
Jiang
,
H.
,
Huang
,
Y.
, and
Hwang
,
K.
,
2005
, “
A Finite-Temperature Continuum Theory Based on Interatomic Potentials: Nanomaterials and Nanomechanics
,”
ASME J. Eng. Mater. Technol.
,
127
(
4
), pp.
408
416
.
59.
Pidin
,
S.
,
Mori
,
T.
,
Inoue
,
K.
,
Fukuta
,
S.
,
Itoh
,
N.
,
Mutoh
,
E.
,
Ohkoshi
,
K.
,
Nakamura
,
R.
,
Kobayashi
,
K.
,
Kawamura
,
K.
,
Saiki
,
T.
,
Fukuyama
,
S.
,
Satoh
,
S.
,
Kase
,
M.
, and
Hashimoto
,
K.
,
2004
, “
A Novel Strain Enhanced CMOS Architecture Using Selectively Deposited High Tensile and High Compressive Silicon Nitride Films
,”
IEDM Technical Digest. IEEE International Electron Devices Meeting
,
IEEE
, pp.
213
216
.
60.
Kumar
,
H.
,
Dong
,
L.
, and
Shenoy
,
V. B.
,
2016
, “
Limits of Coherency and Strain Transfer in Flexible 2D Van Der Waals Heterostructures: Formation of Strain Solitons and Interlayer Debonding
,”
Sci. Rep.
,
6
(
1
), pp.
1
8
.
61.
Levita
,
G.
,
Molinari
,
E.
,
Polcar
,
T.
, and
Righi
,
M. C.
,
2015
, “
First-Principles Comparative Study on the Interlayer Adhesion and Shear Strength of Transition-Metal Dichalcogenides and Graphene
,”
Phys. Rev. B
,
92
(
8
), p.
085434
.
62.
Lloyd
,
D.
,
Liu
,
X.
,
Christopher
,
J. W.
,
Cantley
,
L.
,
Wadehra
,
A.
,
Kim
,
B. L.
,
Goldberg
,
B. B.
,
Swan
,
A. K.
, and
Bunch
,
J. S.
,
2016
, “
Band Gap Engineering With Ultralarge Biaxial Strains in Suspended Monolayer MoS2
,”
Nano. Lett.
,
16
(
9
), pp.
5836
5841
.
63.
Nikam
,
R. D.
,
Sonawane
,
P. A.
,
Sankar
,
R.
, and
Chen
,
Y. -T.
,
2017
, “
Epitaxial Growth of Vertically Stacked P-MoS2/n-MoS2 Heterostructures by Chemical Vapor Deposition for Light Emitting Devices
,”
Nano Energy
,
32
, pp.
454
462
.
64.
Gong
,
L.
,
Young
,
R. J.
,
Kinloch
,
I. A.
,
Riaz
,
I.
,
Jalil
,
R.
, and
Novoselov
,
K. S.
,
2012
, “
Optimizing the Reinforcement of Polymer-Based Nanocomposites by Graphene
,”
ACS.Nano.
,
6
(
3
), pp.
2086
2095
.
65.
Rice
,
C.
,
Young
,
R.
,
Zan
,
R.
,
Bangert
,
U.
,
Wolverson
,
D.
,
Georgiou
,
T.
,
Jalil
,
R.
, and
Novoselov
,
K.
,
2013
, “
Raman-Scattering Measurements and First-Principles Calculations of Strain-Induced Phonon Shifts in Monolayer MoS2
,”
Phys. Rev. B
,
87
(
8
), p.
081307
.
66.
Shafique
,
A.
, and
Shin
,
Y.-H.
,
2017
, “
Strain Engineering of Phonon Thermal Transport Properties in Monolayer 2H-MoTe2
,”
Phys. Chem. Chem. Phys.
,
19
(
47
), pp.
32072
32078
.
67.
Woods
,
C. R.
,
2016
, “
Investigations Into the Interfacial Interaction of Graphene With Hexagonal Boron Nitride
,”
The University of Manchester
,
United Kingdom
.
68.
Kumar
,
H.
,
Er
,
D.
,
Dong
,
L.
,
Li
,
J.
, and
Shenoy
,
V. B.
,
2015
, “
Elastic Deformations in 2D Van Der Waals Heterostructures and Their Impact on Optoelectronic Properties: Predictions From a Multiscale Computational Approach
,”
Sci. Rep.
,
5
(
1
), pp.
1
11
.
69.
Naik
,
M. H.
,
Maity
,
I.
,
Maiti
,
P. K.
, and
Jain
,
M.
,
2019
, “
Kolmogorov–crespi Potential for Multilayer Transition-Metal Dichalcogenides: Capturing Structural Transformations in Moiré Superlattices
,”
J. Phys. Chem. C
,
123
(
15
), pp.
9770
9778
.
70.
Alden
,
J. S.
,
Tsen
,
A. W.
,
Huang
,
P. Y.
,
Hovden
,
R.
,
Brown
,
L.
,
Park
,
J.
,
Muller
,
D. A.
, and
McEuen
,
P. L.
,
2013
, “
Strain Solitons and Topological Defects in Bilayer Graphene
,”
Proc. Natl. Acad. Sci. U. S. A.
,
110
(
28
), pp.
11256
11260
.
71.
Woods
,
C.
,
Britnell
,
L.
,
Eckmann
,
A.
,
Ma
,
R.
,
Lu
,
J.
,
Guo
,
H.
,
Lin
,
X.
,
Yu
,
G.
,
Cao
,
Y.
,
Gorbachev
,
R.
,
Kretinin
,
A.
,
Park
,
J.
,
Ponomarenko
,
L.
,
Katsnelson
,
M.
,
Gornostyrev
,
Y. N.
,
Watanabe
,
K.
,
Taniguchi
,
T.
,
Casiraghi
,
C.
,
Gao
,
H.-J.
,
Geim
,
A.
, and
Novoselov
,
K.
,
2014
, “
Commensurate–incommensurate Transition in Graphene on Hexagonal Boron Nitride
,”
Nat. Phys.
,
10
(
6
), pp.
451
456
.