Stacking-dependent interlayer phonons in 3R and 2H MoS₂

We can describe the frequencies of the breathing (shear) modes by a simple linear chain model in which adjacent layers are coupled harmonically with the same force constant [17]. An N-layer system hosts N  −  1 breathing modes and N  −  1 doubly-degenerate shear modes. Their mode frequencies are:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega _{N,B}^{(\,j)}={{\omega }_{B}}\cos \left( \frac{j\pi }{2N} \right)\nonumber \end{align} \tag{ 1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega _{N,S}^{(\,j)}={{\omega }_{S}}~\cos \left( \frac{j\pi }{2N} \right).\nonumber \end{align} \tag{ 2 }$$

Here j   =  1, 2, ... N  −  1 is the branch index, counting from the highest-frequency to the lowest-frequency branch; ${{\omega }_{B}}$ (${{\omega }_{S}}$) is the frequency of the highest breathing (shear) branch in the bulk limit. These two bulk frequencies are the only fitting parameters in the linear chain model.

Figure 4 compares the observed breathing and shear mode frequencies (symbols) with the prediction of the linear chain model (lines). The observed breathing modes in 3R and 2H MoS₂ can be fit well with a bulk breathing-mode frequency of ${{\omega }_{B}}$  =  53.85 cm⁻¹ and 57.40 cm⁻¹, respectively. The strongest breathing Raman mode corresponds to the lowest breathing branch (j   =  N  −  1), and other weaker breathing Raman modes correspond to higher branches with j   =  N  −  3, N  −  5, N  −  7 ... (solid blue lines in figure 4). Other breathing branches with j   =  N  −  2, N  −  4, N  −  6 ... are Raman inactive (dashed blue lines in figure 4).

Zoom In Zoom Out Reset image size

For the shear modes, the observed Raman frequencies in 3R and 2H MoS₂ can be fit well with a bulk frequency ${{\omega }_{S}}$  =  31.85 cm⁻¹ and 32.85 cm⁻¹, respectively. Both values match closely the measured shear frequency (32.9 cm⁻¹) of bulk 2H MoS₂ (the top spectrum in figures 3(c) and (d)). The linear chain model shows that the observed shear modes in 2H MoS₂ correspond to high-frequency shear branches with j   =  1, 3, 5 (red open dots and solid lines in figure 4). The branches with j   =  2, 4, 6 ... are Raman inactive (red dashed lines). In contrast, the observed shear modes in 3R MoS₂ correspond to low-frequency shear branches with j   =  N  −  1, N  −  3, N  −  5 (red solid dots). The branches with j   =  N  −  2, N  −  4, N  −  6 ... are Raman inactive (red dashed lines). The observed 2H and 3R shear modes therefore cover the high and low shear branches, respectively. For MoS₂ with odd layer numbers, the observable 2H and 3R shear modes complement each other. By combining the 2H and 3R data, we can reveal all the shear branches (j   =  1 to N  −  1) in MoS₂ with odd layer numbers N  =  3, 5, 7.

Our measured interlayer mode frequencies allow us to calculate the interlayer force constants and elastic moduli of 3R and 2H MoS₂ (table 1). The interlayer force constant (K) is related to the bulk interlayer vibration frequency ($\omega$) as $K={{\left( \omega \pi c \right)}^{2}}\mu$, where c is the speed of light and $\mu$  =  3.068  ×  10⁻⁶ kg m⁻² is the mass per unit area of MoS₂ monolayer [17, 31]. By using the extracted 3R and 2H bulk shear mode frequencies, we obtain the interlayer shear force constants: K_x  =  2.76  ×  10¹⁹ Nm⁻³ for 3R MoS₂ and K_x  =  2.94  ×  10¹⁹ Nm⁻³ (K_x and K_y are the same because MoS₂ is isotropic in the xy plane). Similarly, by using the extracted bulk breathing mode frequencies, we obtain the interlayer compressive force constants: K_z  =  7.89  ×  10¹⁹ Nm⁻³ for 3R MoS₂ and K_z  =  8.98  ×  10¹⁹ Nm⁻³ for 2H MoS₂. The elastic modulus (C) is related to the force constant as C  =  Kt, where t is the interlayer separation. 2H and 3R MoS₂ have slightly different interlayer distance: t  =  6.148 Å for 2H MoS₂ and t  =  6.123 Å for 3R MoS₂. Correspondingly, we obtain the shear modulus C₄₄  =  K_xt  =  16.90 GPa (18.08 GPa) for 3R (2H) MoS₂ and the compressive modulus C₃₃  =  K_zt  =  48.32 GPa (55.19 GPa) for 3R (2H) MoS₂ (table 1). Our 2H results agree with prior works [30, 31]. Notably, the 2H compressive modulus is 14% higher than the 3R modulus; this ratio matches roughly the prediction of recent first-principles calculations [69].

Table 1. Elastic parameters of 2H and 3R MoS₂ calculated from the interlayer phonon frequencies.

	Elastic parameters	3R MoS2	2H MoS2
Shear properties	${{\omega }_{S}}$ (cm−1)	31.85	32.85
	${{K}_{x}}$ (1019 Nm−3)	2.76	2.94
	${{C}_{44}}$ (GPa)	16.90	18.08
Compressive properties	${{\omega }_{B}}$ (cm−1)	53.85	57.40
	${{K}_{z}}$ (1019 Nm−3)	7.89	8.98
	${{C}_{33}}$ (GPa)	48.32	55.19

After calculating the elastic parameters, we examine the Raman properties of interlayer phonons. The interlayer phonon modes in 2H MoS₂ have been well studied in the literature [28, 31, 70, 71]. Their Raman activity behavior can be understood from their point group symmetry. 2H MoS₂ with even (odd) layer numbers have the D_3d (D_3h) point group. From the representations of the shear (breathing) modes, we can easily deduce that they are Raman active in every other branch, counting from the highest (lowest) branch, while being Raman inactive in other branches (see the supplementary material for a summary of group-theory analysis for 2H MoS₂). However, the point-group analysis is not applicable in the 3R structure, which has lower symmetry than the 2H structure. 3R MoS₂ has the C_3v point group for all the layer numbers, with no inversion nor mirror symmetry. All breathing modes have the A₁ representation and all shear modes have the E representation. They are all Raman active in the group-theory analysis. Therefore, a simple point-group symmetry analysis cannot explain the observed alternating appearance of the interlayer branches in the 3R Raman spectra.

To reveal the underlying physics of interlayer Raman modes in the 3R structure, we directly examine the influence of stacking order on the Raman tensor of the phonon modes. For a vibration mode k, its Raman intensity (I_k) is proportional to [72, 73]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{I}_{k}}\propto \frac{({{n}_{k}}+1)}{{{\omega }_{k}}}{{\left| {{\widehat{e}}_{i}}\cdot {{\widetilde{R}}_{k}}\cdot \widehat{e}_{s}^{T} \right|}^{2}}.\nonumber \end{align} \tag{ 3 }$$

Here ${{\omega }_{k}}$ is the phonon frequency; ${{n}_{k}}$ is the phonon occupation according to the Bose–Einstein distribution; ${{\widehat{e}}_{i}}$ and ${{\widehat{e}}_{s}}$ are the polarization unit vectors of the incident and scattered light, respectively; ${{\widetilde{R}}_{k}}$ is the Raman tensor. In a classical picture, ${{\widetilde{R}}_{k}}$ can be expressed as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\widetilde{R}}_{k}}=\frac{\partial \widetilde{\alpha }}{\partial {{Q}_{k}}} \Delta {{Q}_{k}}.\nonumber \end{align} \tag{ 4 }$$

Here $\widetilde{\alpha }$ is the polarizability tensor; ${{Q}_{k}}$ is the normal coordinate of the vibration normal mode; $\partial \widetilde{\alpha }/\partial {{Q}_{k}}$ is the derivative of the polarizability tensor at the equilibrium lattice position (${{Q}_{k}}$  =  0).

To determine the Raman activity of a phonon mode, we can directly examine the atomic configuration during the vibration. As an illustration, figure 5(a) shows the interlayer atomic configurations of the two shear modes $S_{3}^{1}$ and $S_{3}^{2}$ in 3L MoS₂ at positive (Q  >  0) and negative (Q  <  0) normal displacement. The 3R and 2H structures have distinct atomic configurations during the vibration. For the 3R $S_{3}^{1}$ mode, the interlayer atomic configurations are the same for Q  >  0 and Q  <  0. The differential polarizability $\partial \alpha /\partial Q$ is thus zero; this mode is Raman inactive. But the $S_{3}^{1}$ mode in the 2H structure has distinct atomic configurations for Q  >  0 and Q  <  0. This mode thus has non-zero $\partial \alpha /\partial Q$ and is Raman active. In contrast, the lower $S_{3}^{2}$ mode has opposite stacking dependence. The 3R $S_{3}^{2}$ mode exhibits different interlayer atomic configurations for Q  >  0 and Q  <  0, leading to finite $\partial \alpha /\partial Q$ and Raman response. However, the 2H $S_{3}^{2}$ mode has the same atomic configurations for Q  >  0 and Q  <  0, leading to zero $\partial \alpha /\partial Q$ and no Raman response. This simple analysis illustrates why the high and low shear branches appear exclusively in the Raman spectra from 3L MoS₂ of 2H and 3R phases, respectively (figure 5(b)). The above pictorial analysis can be applied to explain the Raman activity of shear and breathing modes in 2H and 3R structures with any number of layers (e.g. see the analysis of 5L MoS₂ in the Supplementary Material). Similar stacking-dependent shear modes are also observed in few-layer graphene with ABA and ABC stacking order [18, 26].

Zoom In Zoom Out Reset image size

Our qualitative analysis above can be quantified in an effective bond polarizability model [72, 73], which predicts the relative Raman intensity of different shear modes. This model considers the distortion of interlayer bonds by the interlayer vibration. The Raman tensor $\widetilde{R}$ is the sum of contributions from each interlayer bond. The diagonal element ${{R}_{xx}}$ is responsible for parallel-polarized Raman response; the off-diagonal element ${{R}_{xy}}$ is responsible for cross-polarized Raman response. For a shear mode k with x-direction layer displacement in an N-layer structure, the tensor element ${{R}_{\mu \nu }}$ ($\mu ,\nu =x,y,z$) may be expressed as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{\mu \nu }}(k)=\underset{l=1}{\overset{N-1}{\mathop \sum }}\,\alpha _{l}^{\prime}({{x}_{l}}-{{x}_{l+1}}).\nonumber \end{align} \tag{ 5 }$$

Here $l$ denotes the layer index; ${{x}_{l}}$ is the lateral displacement of layer l from the equilibrium position; ${{\alpha }_{l}}$ is the $\mu \nu$-element of the polarizability tensor of the interlayer bond between layers l and l  +  1. $\alpha _{l}^{\prime}=\partial {{\alpha }_{l}}/\partial {{x}_{l}}$ is the derivative of ${{\alpha }_{l}}$ with respect to ${{x}_{l}}$. For a shear mode $S_{N}^{\,j}$ with branch index j  in an N-layer structure, the displacement of layer l is: [17]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{x}_{l}}\propto \cos \left[ \frac{\left( N-j \right)\left( 2l-1 \right)}{2N}\pi \right].\nonumber \end{align} \tag{ 6 }$$

The 3R and 2H structure have the same layer displacement pattern in vibration, but different configurations of interlayer bonds. In 3R MoS₂, all the layers have the same orientation and identical interlayer bonds, and hence the same $\alpha _{l}^{\prime}$ (figure 6(a)). We may thus remove the l index and name it as $\alpha _{l}^{\prime}=\alpha _{R}^{\prime}$. By equation (5), the 3R Raman tensor element is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{\mu \nu }}\left( k \right)=&~\alpha _{R}^{\prime}\left[ \left( {{x}_{1}}-{{x}_{2}} \right)+\left( {{x}_{2}}-{{x}_{3}} \right)+\ldots \right.\nonumber \\&\left.+\left( {{x}_{N-1}}-{{x}_{N}} \right) \right]=\alpha _{R}^{\prime}\left( {{x}_{1}}-{{x}_{N}} \right).\nonumber \end{align} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

It depends only on the displacement of the top and bottom layers.

In contrast, the adjacent layers in 2H MoS₂ have opposite orientation. The interlayer bonds thus exhibit alternate lateral orientation for successive layer pairs, leading to opposite signs of $\alpha _{l}^{\prime}$ (figure 6(b)). By setting $\alpha _{l=1}^{\prime}=\alpha _{H}^{\prime}$, the other $\alpha _{l}^{\prime}$ would be ${{(-1)}^{l+1}}\alpha _{H}^{\prime}$. By equation (5), the 2H Raman tensor element is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{\mu \nu }}\left( k \right)&=\underset{l=1}{\overset{N-1}{\mathop \sum }}\,{{(-1)}^{l+1}}\alpha _{H}^{\prime}\left( {{x}_{l}}-{{x}_{l+1}} \right) \nonumber \\ & =\alpha _{H}^{\prime}\left[ \left( {{x}_{1}}-{{x}_{2}} \right)-\left( {{x}_{2}}-{{x}_{3}} \right)+\ldots \left( {{x}_{N-1}}-{{x}_{N}} \right) \right]. \nonumber \end{align} \tag{ 8 }$$

We can combine equation (3), (6)–(8) to obtain the Raman intensity. After some algebraic derivation (see supplementary material), the Raman intensities of the 3R and 2H shear modes $S_{N}^{\,j}$ show the following analytic forms:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{I}_{3R}}\left( k \right)\propto \frac{{{n}_{k}}+1}{{{\omega }_{k}}}{{\left| \alpha _{R}^{\prime} \right|}^{2}}{{\sin }^{2}}\left[ \frac{\left( N-j \right)\left( N-1 \right)}{2N}\pi \right]\left[ 1-{{\left( -1 \right)}^{N-j}} \right]\nonumber \end{align} \tag{ 9 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{I}_{2H}}\left( k \right)\propto \frac{{{n}_{k}}+1}{{{\omega }_{k}}}{{\left| \alpha _{H}^{\prime} \right|}^{2}}{{\sin }^{2}}\left( \frac{N-j}{2N}\pi \right)\nonumber {{\tan }^{2}}\left( \frac{N-j}{2N}\pi \right)\left[ 1-{{\left( -1 \right)}^{\,j}} \right].\nonumber \end{align} \tag{ 10 }$$

In equation (9), the factor $1-{{\left( -1 \right)}^{N-j}}$ is zero when $N-j$ is an even number; the 3R shear mode is non-zero only when the branch index j   =  N  −  1, N  −  3, N  −  5, .... In equation (10), the factor $1-{{\left( -1 \right)}^{\,j}}$ is zero when j  is an even number; the 2H shear mode is non-zero only when j   =  1, 3, 5, .... Such stacking-dependent Raman activity exactly matches our experimental results (figures 3 and 4).

In addition to accounting for the Raman activity, the above analytical theory enables us to quantitatively simulate the Raman spectra of shear modes for each layer number and stacking order. In the simulation, we use the linear chain model (equations (1) and (2)) to determine the mode frequencies and the bond polarizability model (equations (9) and (10)) to determine the Raman intensity ratio of different branches. The calculation gives the same results for parallel and cross polarizations. For comparison with experiment, we plot the theoretical Raman modes as Gaussian functions with the same width as the experimental peaks. Figures 5 and 7 compare our simulations with experiment for 3L, 5L and 9L MoS₂, which exhibit one, two and three shear modes, respectively, for both the 3R and 2H phases in our observed Raman spectra. Our calculated spectra show reasonably accurate intensity ratio between different shear branches (see supplementary material for simulations for other layer numbers).

Zoom In Zoom Out Reset image size

Our model can be readily extended to explain the breathing modes by redefining $\alpha _{l}^{\prime}=\partial {{\alpha }_{l}}/\partial {{z}_{l}}$, where ${{z}_{l}}$ is the vertical displacement of layer l. In this case, $\alpha _{l}^{\prime}$ does not depend on the lateral orientation of the interlayer bonds, so $\alpha _{l}^{\prime}$ is the same for all layer indices l in either stacking order. As a result, for both 2H and 3R breathing modes, their Raman tensor and intensity follow the same forms as those of 3R shear modes in equations (7) and (9). For instance, figure 8 displays our calculated breathing-mode Raman spectrum for 10L MoS₂ (the same for 2H and 3R phases), in comparison with the measured parallel-polarized (VV) spectrum of 10L 3R MoS₂. Our calculation quantitatively reproduces the experimental spectrum, which exhibits four breathing modes ($B_{10}^{9}$, $B_{10}^{7}$, $B_{10}^{5}$, $B_{10}^{3}$) from low to high frequency and strong to weak intensity.

Zoom In Zoom Out Reset image size