Strain-shear coupling in bilayer MoS2.

Layered materials such as graphite and transition metal dichalcogenides have extremely anisotropic mechanical properties owing to orders of magnitude difference between in-plane and out-of-plane interatomic interaction strengths. Although effects of mechanical perturbations on either intralayer or interlayer interactions have been extensively investigated, mutual correlations between them have rarely been addressed. Here, we show that layered materials have an inevitable coupling between in-plane uniaxial strain and interlayer shear. Because of this, the uniaxial in-plane strain induces an anomalous splitting of the degenerate interlayer shear phonon modes such that the split shear mode along the tensile strain is not softened but hardened contrary to the case of intralayer phonon modes. We confirm the effect by measuring Raman shifts of shear modes of bilayer MoS2 under strain. Moreover, by analyzing the splitting, we obtain an unexplored off-diagonal elastic constant, demonstrating that Raman spectroscopy can determine almost all mechanical constants of layered materials.

Introduction

Three-dimensional (3D) layered materials such as graphite and transition metal dichalcogenides (TMDs) are formed by weak van der Waals force between two-dimensional (2D) crystals whereas atoms within each 2D crystal are bound through a strong covalent bonding. A huge difference between in-plane and out-of-plane interatomic interaction strengths often result in extremely anisotropic mechanical properties of 3D-layered systems. Since electronic and mechanical properties of multilayered 2D crystals critically depend on their interlayer interactions[1-3], understanding the influence of the external mechanical perturbations on the interlayer as well as intralayer interaction is very important. Mechanical properties of a single-layer 2D crystal can be understood by analyzing variations of high-frequency optical phonon modes with external mechanical perturbations[4-12]. Similarly, the effects of interlayer interactions on mechanical properties can be analyzed by probing the low-frequency interlayer shear and breathing modes[13-17]. Raman spectroscopy, an invaluable diagnostic tool for 2D crystals, have played an important role in measuring these modes for various layered materials such as graphene and TMDs[13-28]. For example, intralayer[11,12] and interlayer[13-17] elastic moduli as well as their variations due to rotational and translational stacking faults[29,30] can be obtained by examining the ultralow-frequency Raman spectrum of layered 2D materials. As introduced here, the effects of mechanical perturbations on physical properties related with either intralayer or interlayer interactions separately have been studied extensively whereas mutual interplay between them have rarely been studied.

In this study, we attempt to elucidate the overall effect of external mechanical strain on the interlayer and intralayer interactions, as well as the coupling between them. We first analyze the effect of tensile strain on bilayer MoS2 using the analytic linear elastic model. Then, we perform first-principles calculations on elastic constants of bilayer MoS2, a prototypical TMD material. Both the analytic model and the first-principles calculations predict that a uniaxial external mechanical strain along in-plane direction indeed couples to interlayer shear, splitting the degenerate low-energy shear phonon modes. This indicates relative sliding of two layers with an external strain along in-plane direction. Contrary to usual phonon softening under tensile strain, our calculations show that the softened shear mode in strained bilayer MoS2 is perpendicular to the applied tensile strain direction and vice versa. Next, we confirm the predicted frequency shift and splitting of the shear modes of bilayer MoS2 by carrying out ultralow-frequency Raman scattering measurements under uniaxial strain. From a careful analysis on the polarized Raman spectroscopy data on split shear phonon modes, we can measure the off-diagonal elastic constant of the bilayer system that has not yet been explored. Our results could therefore provide a way to measure a whole set of elastic constants of layered systems that characterize their mechanical properties completely.

Results

Coupling between uniaxial strain and shear

The bilayer 2H-MoS2 has a stacked trigonal prismatic structure with inversion and three-fold rotational symmetry (Fig. 1a, b). So, the compliance tensor of the system with D 3d symmetry in a matrix form[31] is given by

Fig. 1

Lattice structure of bilayer MoS2 under uniaxial strain. a Definitions of angles. b Side view of crystal structure of bilayer MoS2. The atoms in the top layer are shown as solid circles, whereas those in the bottom layer are shown as hashed circles. c Top view of unstrained single-layer MoS2 (left panel) and distorted lattice structure under armchair strain (right panel). Solid lines are for Mo–S covalent bonds. The double lines denote a longer bond distance compared to the single line. d Top view of hexagonal network of interface sulfur atoms in unstrained bilayer MoS2 (left) and distorted lattice structure under zigzag strain (right). Sulfur atoms belonging to the upper (lower) layer are denoted by yellow filled (hashed) circles. The dotted lines indicate interlayer van der Waals bonds between interface sulfur atoms. Here, the double and single lines indicate the lengthened and shortened interfacial sulfur bond distances, respectively. Zigzag strain in single-layer MoS2 and armchair strain in bilayer MoS2 induce similar lattice distortions

Here, E i(o) is the in-plane (out-of-plane) Young’s modulus, ν i(o) is the in-plane (out-of-plane) Poisson’s ratio, and the vector form of strain and stress are defined as ε 1 = ε xx, ε 2 = ε yy, ε 3 = ε zz, ε 4 = 2ε yz, ε 5 = 2ε xz, ε 6 = 2ε xy, σ 1 = σ xx, σ 2 = σ yy, σ 3 = σ zz, σ 4 = σ yz, σ 5 = σ xz, and σ 6 = σ xy. The tensor form of strain is further defined as ε ij = 1/2(∂u i/∂x j+∂u j/∂x i) with u i being the local displacement vector component, and the corresponding stress tensor is σ ij = ∂F/∂ε ij, where F is the free energy of the system[31]. A constant uniaxial tensile stress σ along a direction at an angle φ from the x-axis in Fig. 1a corresponds to a stress tensor with components σ 1 = σ cos2 φ, σ 2 = σ sin2 φ, σ 6 = σ cosφ sinφ, and σ 3 = σ 4 = σ 5 = 0, using rotational tensor transformation of σ ij. From Eq. (1), we find that such a stress leads to strain components and . These correspond to lateral sliding of the upper layer with respect to the lower layer by Δr = (Δx, Δy) with and . This implies that the uniaxial stress causes lateral sliding of the upper layer with respect to the lower layer in the direction at an angle of −2φ from the positive y-axis if s 14 > 0. If s 14 < 0, the direction of sliding is reversed toward the angle −2φ from the negative y-axis, which is π−2φ from the positive y-axis.

For a simple case of a constant uniaxial tensile stress σ along the x-direction (φ = 0°) as shown in Fig. 1d, σ 1 = σ and σ 2 = … = σ 6 = 0, which lead to an off-diagonal strain component ε 4 = s 14  σ. This indicates that the stress can induce a lateral sliding between the two layers given by , where d int is the interlayer distance.

First-principles calculations of shear mode splitting

Applying a tensile strain typically softens the phonons along the direction of the strain because the effective spring constant along it becomes weaker owing to the elongation of interatomic distances[4,6-10,32]. Such a trend is found in our first-principles calculations and other studies for the intralayer optical phonon modes of strained MoS2 [7-10]. For the intralayer E g mode of bilayer MoS2, the calculated frequency shifts match our experiments very well and are found to be linear up to 2% strain (see Supplementary Fig. 1). However, opposite results are obtained for the interlayer shear modes. Figure 2 shows the calculated phonon frequency shift of the interlayer shear modes of the same system with respect to the uniaxial strain; the higher frequency mode (S +) corresponds to the vibrations along the strain axis, whereas the lower frequency mode (S −) corresponds to the vibrations perpendicular to it. Our first-principles calculations show that the shifts barely depend on the strain direction (see Fig. 2 and Supplementary Fig. 2). By inspecting the relaxed atomic structures for phonon calculation under strain, we have found that the two layers slide with respect to each other as uniaxial strain is applied. This indeed agrees with the aforementioned analysis of the compliance tensor of the system. The amount of strain-induced interlayer sliding based on the first-principles calculation is not negligible at all and is about a half percent of the lattice constant for the applied strain of 1%. It is also confirmed that the direction of sliding matches well with our model prediction.

Fig. 2

Splitting of interlayer shear mode under uniaxial strain. Shifts of split interlayer shear modes under uniaxial tensile strain for two different strain directions of φ = 15.7° (red) and 23.8° (blue). Calculation (open circle and square) and experimental (filled symbols) results are compared. For the experimental data, values obtained from the Stokes and anti-Stokes spectra are plotted. Experimental error bars in strain come from the uncertainty in the estimate of the strain using the equation ε = t 0 θ/L 0 as explained in Methods. The error bars in Δω are defined by the standard deviation of experimental values

We have found that the stress-induced layer sliding is the origin of the anomalous shear mode splitting. We first assume that the interlayer coupling for the shear modes is mainly determined by the interaction between adjacent sulfur atoms at the interface. These interface sulfur atoms also form a buckled hexagonal lattice, where each sublattice belongs to a different sulfur layer (Fig. 1d). Therefore, the splitting of shear modes would solely depend on the way of breaking the hexagonal symmetry in the hexagonal network between the nearest sulfur atoms at the interface between the top and bottom layers of bilayer MoS2. From the compliance tensor analysis explained earlier, the strain along the zigzag direction (σ 1 = σ and otherwise zero), for example, induces an interlayer shift along the armchair direction (Fig. 1d). Such a layer sliding along the armchair direction breaks the symmetry in the interface hexagonal lattice similar to the distorted intralayer hexagonal lattice under strain along the armchair direction (Fig. 1c). Thus, the effective bond distance between interface sulfur atoms parallel to the strain direction is shortened, while the perpendicular one is elongated as schematically illustrated in Fig. 1. This is indeed opposite to the in-plane bond elongation within a single layer under the same strain. Thus, one can expect that the polarization dependence of Raman scattering intensity of the hardened shear mode should be similar to that of softened intralayer optical phonon modes and vice versa. From our first-principles calculations, we have also checked that the frequency splitting for the case of uniaxial stretching without sliding. It has shown that the splitting becomes ten times smaller compared to the splitting when the sliding is properly considered (see Supplementary Fig. 2). We developed a detailed interface phonon model using effective interatomic interaction between the sulfur atoms at the interface and obtain a quantitative expression for the splitting of shear phonon modes under strain agreeing well with our first-principles calculation results (Supplementary Note 1 and Supplementary Fig. 3). From all theoretical considerations, we expect that the anomalous splitting of shear phonon modes under uniaxial strain should occur irrespective of direction of applied in-plane strain.

Raman spectroscopic determination of shear mode splitting

We have performed polarized Raman experiments to compare with our theoretical prediction. Bilayer MoS2 samples are prepared directly on acrylic substrates by mechanical exfoliation from MoS2 flakes (Supplementary Fig. 4). The number of layers is determined by combination of optical contrast, Raman and photoluminescence (PL) measurements. Experimental details of application of strain and polarized ultralow-frequency Raman measurements are presented in “Methods” section.

Figure 3a shows the Raman spectra of bilayer MoS2 for different uniaxial strains up to 1.88%. The A 1g and E g modes of bilayer MoS2 correspond to intralayer vibrations, the interlayer in-plane mode (shear mode, S) corresponding to E g and the interlayer out-of-plane A 1g mode (breathing mode, B) are also observed (see Supplementary Fig. 4a for schematics of vibration modes)[14,15]. First, we focus on the intralayer high-frequency modes under strain. The in-plane E g mode redshifts and splits into two peaks as the strain increases, whereas the A 1g mode does not vary much[7-10]. Since the A 1g mode is an out-of-plane mode, the effect of in-plane strain is minimal[7-10]. The E g mode splits owing to the symmetry breaking of the hexagonal lattice[4,7,10,32]. We label the split mode with lower (higher) frequency as E g − (E g +). We note that since the strain lowers symmetry of the system, the vibration modes of the strained MoS2 are no longer E g [8,28]. Nevertheless, we label the split modes as E g − and E g + in order to indicate their origin. The overall behavior of intralayer Raman modes of bilayer MoS2 is similar to that of the single-layer case[7,10]. The interlayer shear (S) and breathing (B) modes show negligible dependence on strain. These results are summarized in Fig. 3b. Assuming linear dependence of the frequencies on the uniaxial strain, the shift rates of the E g −, E g +, and A 1g modes for bilayer are obtained to be −4.0 ± 0.1, −0.9 ± 0.2, and −0.3 ± 0.1 cm−1%−1, respectively, agreeing well with previous studies[7]. We compared the experimental data with calculations with different in-plane Poisson’s ratio (see Supplementary Fig. 1). The experimental results seem to fit best with calculations with the Poisson’s ratio of intrinsic bilayer MoS2 (ν i = 0.22)[33], which we used for the following analysis. From the shift rates, the Grüneisen parameter (γ) and the shear deformation potential (β) for the E g modes are calculated as γ = 1.0 ± 0.2 and β = 0.8 ± 0.2, which are also similar to previous results[7] (Supplementary Note 3). The polarization dependences of the in-plane and out-of-plane intralayer modes with strain are shown in Fig. 4a and summarized in Fig. 4b, c, respectively. By inspecting the polarized Raman spectrum of in-plane intralayer modes, we can determine the direction of strain with respect to the crystal orientation following the well-established procedure[4] (Supplementary Note 4). Analyzing the data in Fig. 4b, we obtain the crystallographic orientation of the sample, φ = 15.7 ± 1.1°. For another sample, we obtained φ = 23.8 ± 2.1° (Supplementary Fig. 5).

Fig. 3

Raman spectra of bilayer MoS2 under uniaxial strain. a Raman spectra of bilayer MoS2 as a function of uniaxial strain. From left to right, the shear, breathing, split E g, and A 1g modes are denoted by the triangle, square, diamond, and circle symbols, respectively. b The peak positions for intralayer phonon modes (upper panel) and interlayer shear and breathing modes (lower panel) as a function of uniaxial tensile strain. Error bars in strain come from the uncertainty in the estimate of the strain using the equation ε = t 0 θ/L 0. The error bars in Raman shift are defined by spectral resolution of the spectrometer

Fig. 4

Polarization dependence of intralayer Raman modes under strain. a Polarized Raman spectra of strained bilayer MoS2 (ε = 1.28%) for main intralayer modes. The incident polarization (θ i) is fixed at 0° and the spectra are measured as a function of ψ defined in Fig. 1a. Empty and filled diamond and circle symbols are defined in Fig. 3. Normalized polar plots of b E g (E g − and E g +) and c A 1g modes

Now we turn our focus onto the interlayer shear modes. Figure 5a, b compare the polarization dependence of the low-frequency Raman modes without and with uniaxial strain, respectively. The intensity of the breathing mode (B) shows a dependence on the polarization similar to the case of the A 1g mode irrespective of strain. Under strain, the shear mode peak intensity does not seem to depend on polarization, and no apparent splitting is observed at any given polarization. Upon close inspection, however, one can notice that the shear mode peak at ~22.8 cm–1 seems to move periodically with polarization unlike the unstrained case as shown in Fig. 5a, b. Since the phonon mode frequency itself should not depend on polarization, the apparent polarization-dependent shift can only be interpreted as being due to polarization dependence of the relative intensities of closely spaced peaks due to strain-induced splitting. This kind of a small splitting can be distinguished with a help of polarization[34,35] because degeneracy-lifted phonon modes have different polarization dependences[4,6,10,32]. Since the peaks are not resolved well at any particular polarization angle, a simple double Lorentzian fitting of individual spectrum would not be able to determine the splitting reliably. Therefore, we measured a whole set of 37 spectra as a function of polarization at fixed strain of 1.28%. Since the frequency of the modes should not depend on the polarization, the whole set of 37 spectra were fitted as a whole by requiring that the peak positions are the same in all 37 spectra and only the relative intensities vary (Fig. 5c). By doing so, one can greatly reduce the experimental uncertainty in the positions of the two peaks as shown in Fig. 5c and find that the peak splitting is ~0.91 ± 0.05 cm–1 under 1.28% uniaxial strain. The validity of this procedure is ascertained by Fig. 5d, in which the intensities of the split peaks obtained through the fitting procedure follow the dependences expected from strain-split peaks of E g phonon modes (Supplementary Note 4).

Fig. 5

Polarization dependence of interlayer Raman modes. Polarized Raman spectra of bilayer MoS2 a without and b with strain (ε = 1.28%) in the low-frequency region. The incident polarization (θ i) is fixed at 0° and the spectra are measured as a function of ψ defined in Fig. 1a. The insets compare spectra taken for ψ = 40° (blue) and 130° (red) and show the splitting for the case with strain. c Representative fitting results for split shear modes. Normalized polar plots of d shear (S, S −, and S +) and e breathing (B) modes

For the interlayer shear modes, we denote the low-frequency modes and high-frequency modes as S − and S +, respectively, as mentioned earlier. The polarization dependence of the shear and breathing modes are shown in Fig. 5d, e. We also find almost similar shifts in other samples with different strain directions (Supplementary Fig. 5). The experimentally obtained peak positions under strain are plotted in Fig. 2, which show good agreement with our first-principles calculations. As expected from our theoretical analysis, the splitting behavior of the interlayer shear modes is indeed in sharp contrast to the intralayer case in that the split mode in the direction parallel to the uniaxial strain hardens while the other mode in the direction perpendicular to it softens with strain. The polarization dependence of S + (Fig. 5d) is the same as that of E g − (Fig. 4b), and vice versa, confirming our theoretical analysis. The corresponding Grüneisen parameter and shear deformation potentials for the interlayer shear phonon are obtained to be 0.9 ± 0.3 and 2.4 ± 0.4, respectively.

Discussion

An off-diagonal elastic modulus s 14 can be obtained by comparing the experimental data with the model based on the interface hexagonal sulfur networks. With a simple spring model between the interface sulfur atoms, the amount of frequency splitting can be estimated to be , where ω 0, d int, E i(=1/s 11), and ε are the phonon frequency without strain, the interlayer distance, the in-plane Young’s modulus, and the strain along the applied uniaxial stress, respectively; k is the lateral component of the effective spring constant between nearest-neighbor interface sulfur atoms from different layers; is the linear scaling factor of k to the change of interatomic distance from the equilibrium; and s 14 is the element of the compliance tensor between ε 1 and σ 4 (see Supplementary Note 1 for derivation of Δω). On the other hand, the average shift of the two split frequencies is given by (see Supplementary Note 1 for derivation), where a SS is the interlayer sulfur-to-sulfur lateral distance, and v i is the in-plane Poisson’s ratio. Combining the expressions for Δω and , we obtain . With a ss = 0.18 nm, d int = 0.62 nm, ν i = 0.22, E i = 0.33 TPa[33,36,37] and the experimental values of Δω/ε = 0.60 ± 0.13 cm−1%−1, and Δ/ε = −0.15 ± 0.04 cm–1%–1, we find that of MoS2 is 1.46 ± 0.34 TPa–1. Most of the uncertainty comes from a small Δ value that goes into the denominator. Within this uncertainty, our experimentally obtained value is consistent with our first-principles calculation result of 0.84 TPa–1, which is obtained following a procedure similar to that in ref. [33]. From these considerations, we tabulate all relevant elastic parameters in Table 1, and the complete matrix elements of compliance tensor in Eq. (1) for bilayer MoS2 can be written as:where the matrix elements are in the unit of TPa–1. Here, we used experimental values wherever available.

Table 1

Elastic constants of bilayer MoS2

Elastic constants	Value	References
s 14	−1.46 ± 0.34 TPa−1 (exp)/−0.84 TPa−1 (cal)	This work
s 44	111 ± 8 TPa−1 (exp)a/154 TPa−1 (cal)	This work
E i	0.33 TPa (exp)	Ref. [36]
E o	0.058 ± 0.002 TPa (exp)b	This work
ν i	0.22 (cal)	Ref. [33]
ν o	0.18 (cal)	Ref. [33]

Note: Summary of elastic constants estimated from experimental results

aExperimental value of s 44 was estimated from c 44, obtained from the position of the (unstrained) shear modes14,15, assuming that off-diagonal elements are much smaller than diagonal elements (s 44≈1/c 44)

bExperimental value of E o was obtained from the position of the (unstrained) breathing mode[14,15]

By analyzing the low-frequency Raman spectrum, we have shown that the hitherto unexplored off-diagonal elastic constant s 14 of MoS2 can be estimated. Our study here can be easily extended to other multilayered 2D crystals and could open a new way to determine almost all elastic constants of layered materials.

Methods

Computational methods

For the theoretical analysis, we have carried out first-principles calculations using Quantum ESPRESSO package[38] with plane wave basis and norm-conserving pseudopotentials[39]. To include the interlayer van der Waals interaction properly, we used the revised version[40] of the nonlocal correlation functional method developed by Vydrov and van Voorhis[41]. Phonon frequencies are calculated using density functional perturbation theory[38,42]. The energy cutoff for the basis set expansion is 110 Ry. Such a high-energy cutoff is adopted for resolving the low frequencies of interlayer shear modes which is order of a few meV. The k-point grid of 8 × 8 × 1 is used. In order to obtain s 14, we have alternatively calculated s 56 (=−2s 14) by inverting the matrix form of the shear part of stiffness tensor c ijkl, where i, j, k, l = x, y, z. The stiffness tensor is defined by c ijkl = ∂ 2 F/∂ε ij ∂ε kl. Here, F is the total free energy of the system from the first-principles calculations as a function of discrete values of strain ε ij up to 5%, and the differentiation is done numerically. The thickness of a bilayer MoS2 is set to be twice the interlayer distance. For more details, refer to Supplementary Note 2.

Experimental details of Raman measurements

The bilayer samples of MoS2 were prepared directly on acrylic substrates by mechanical exfoliation from MoS2 flakes (SPI supplies). The number of layers was determined by the combination of optical contrast, Raman and PL measurements (Supplementary Fig. 4). The uniaxial strain is applied by a four-probe bending stage and calculated by ε = t 0 θ/L 0, where L 0 and t 0 is unstrained length and thickness of acrylic substrate, θ is an angle in arc by bending substrate[6-8]. The laser beam was focused onto the sample by a 50× microscope objective lens (0.8 N.A.), and the scattered light was collected and collimated by the same objective. The scattered signal was dispersed with a Jobin-Yvon Triax 550 spectrometer (1800 grooves mm–1) and detected with a liquid-nitrogen-cooled back-illuminated charge-coupled device detector. To access the low-frequency range below 100 cm−1, reflective volume holographic filters (Ondax) were used to reject the Rayleigh-scattered light. The spectral resolution of our system is ~0.7 cm−1. The laser power was kept below 0.2 mW in order to avoid heating.

Data availability

The data that support the findings of this study are available from the corresponding author upon request.