Stacking sequence determines Raman intensities of observed interlayer shear modes in 2D layered materials – A general bond polarizability model

2D layered materials have recently attracted tremendous interest due to their fascinating properties and potential applications. The interlayer interactions are much weaker than the intralayer bonds, allowing the as-synthesized materials to exhibit different stacking sequences, leading to different physical properties. Here, we show that regardless of the space group of the 2D materials, the Raman frequencies of the interlayer shear modes observed under the typical configuration blue shift for AB stacked materials, and red shift for ABC stacked materials, as the number of layers increases. Our predictions are made using an intuitive bond polarizability model which shows that stacking sequence plays a key role in determining which interlayer shear modes lead to the largest change in polarizability (Raman intensity); the modes with the largest Raman intensity determining the frequency trends. We present direct evidence for these conclusions by studying the Raman modes in few layer graphene, MoS2, MoSe2, WSe2 and Bi2Se3, using both first principles calculations and Raman spectroscopy. This study sheds light on the influence of stacking sequence on the Raman intensities of intrinsic interlayer modes in 2D layered materials in general, and leads to a practical way of identifying the stacking sequence in these materials.


Results
Interlayer modes illustrated by few layer graphene. We begin by discussing the interlayer modes present in 2D layered materials. In a 2D material with N layers, there are N times the number of normal modes as that in the monolayer. Each normal mode in 1 layer (1L) evolves into N modes with slightly different frequencies in N layers (NL), keeping the same intralayer displacement while varying the phase difference between adjacent layers 12 . The difference in frequencies arises from the difference in relative phases in adjacent layers. In this way, the interlayer modes in NL correspond to the acoustic mode in 1L, in which all atoms within the single layer move together. Since the interlayer interactions are weak, the corresponding frequencies are also low.
As an illustration, we show in Fig. 1 the ultralow frequency Raman spectra of AB and ABC stacked FLG with varying thickness, computed using density functional theory (DFT) within the local density approximation (LDA) (see Methods for details). These LDA calculated frequencies are in excellent agreement with experiments (Supplementary Figure S1). We note that unless otherwise stated, the Raman intensities discussed in this work refer to the non-resonant Raman intensities. The Raman intensities for each line are normalized by the largest value in that system. Unless otherwise mentioned, we shall consider in this work the Raman intensities corresponding to those obtained in the most common ( ) z xx z polarization configuration. We use the notation "S" to label the interlayer shear modes and the subscripts, 0 to (N-1), denote in order the lowest to highest frequency shear modes in each system, with S 0 corresponding to the acoustic mode, S 1 the lowest frequency shear mode and S N-1 the highest frequency shear mode. Supplementary Tables S1 and S2 show the frequencies of all the interlayer modes in AB and ABC stacked FLG respectively -the stacking order has negligible impact on the frequency values, indicating that the interlayer force constants are similar in both systems (see also caption of Supplementary Figure  S1). However, distinct frequency trends are found in AB and ABC stacked FLG (Fig. 1b). The peaks with largest intensity correspond to S N-1 and S 1 in AB and ABC stacked FLG respectively. The other trends (S N-3 and S N-5 in AB stacked FLG, and S 3 and S 5 in ABC stacked FLG) are barely visible in Fig. 1b, but can be seen from the details in Supplementary Tables S1 to S2.
As we have seen, the subscript of the S mode determines its relative Raman intensity. Furthermore, this subscript, i.e. whether a mode is the lowest or highest frequency shear mode, or the third lowest or third highest frequency mode, dictates whether the frequencies red shift or blue shift with increasing thickness. In general we find that the lower-frequency shear modes red shift with increasing thickness, and the opposite is true for the higher-frequency shear modes. To understand this, we consider the Scientific RepoRts | 5:14565 | DOi: 10.1038/srep14565 atomic displacements of the S N-1 , S N-3 , S 1 and S 3 modes as shown in Fig. 2. First we note that the highest frequency mode S N-1 corresponds to maximally out-of-phase displacement between adjacent layers, while the lowest frequency mode S 1 corresponds to minimum out-of-phase displacement. Since the adjacent graphene layers vibrate out of phase (180 o ) in the S N-1 modes, the restoring forces accumulate with increasing thickness, resulting in the blue shift with thickness. In contrast, for the S 1 modes, the graphene layers are equally divided into two parts, with one part moving in one direction and the other part moving in the opposite direction, and the atomic displacements are gradually reduced towards the interface of the two parts. The relative displacement between adjacent layers therefore decreases as the number of layers increases, resulting in a red shift. Similar analysis can apply to the S N-3 and S 3 modes.
For completeness, we note that besides interlayer shear modes, interlayer breathing modes are also present; the atomic displacements being similar to those of the shear modes, but in the out-of-plane direction. In FLG, the intensities of such modes are much smaller than the shear modes. We note that unlike the shear modes, the stacking order has no influence on the relative Raman intensities of the breathing modes (Supplementary Figure S1, Tables S1-S2).
It is interesting to note that the stacking order dependent frequency trends for FLG are consistent with the trends for AB stacked TMD and ABC stacked Bi 2 Se 3 and Bi 2 Te 3 that are reported in the literature [12][13][14]21,23,24 -in other words, the frequency of the observed interlayer shear mode blue (red) shifts for AB (ABC) stacked materials as the number of layers increases. Just as in the case of FLG, these frequency trends are determined by which shear modes (lowest or highest frequency) have the largest Raman intensities [12][13][14]23,25 . The observation of these same trends in vastly different materials with similar stacking order begs the question of whether the stacking order (AB or ABC stacking) can in general determine the Raman intensities of the interlayer shear modes and therefore the observed frequency trends.
The Raman activity of a phonon mode is assigned using group theory analysis, which tells us about the symmetries present in the material. It is natural to question whether the space group (symmetries) of these 2D layered materials can also explain the stacking-order dependent Raman intensities. We summarize in Table 1 the space group (symmetries) of multi-layered graphene, MoS 2 (MoSe 2 and other 2H TMD) and Bi 2 Se 3 (Bi 2 Te 3 ) under AB and ABC stacking, with different number of layers. It is clear that the space group (symmetry) does not explain the Raman intensities for the interlayer shear modes. For example, both AB and ABC stacked graphene with an even number of layers have the same symmetry, even though the interlayer shear modes for these two systems have very different Raman intensities. Likewise, ABC stacked Bi 2 Se 3 and Bi 2 Te 3 have the same space group (D 3 3d ) as AB and ABC stacked graphene with an even number of layers, but the observed frequency trends follow only those of ABC stacked graphene. On the other hand, although ABC stacked graphene (D 3 3d ) and ABC stacked MoS 2 and MoSe 2 (C 1 3v ) have different symmetry, the Raman spectra for their interlayer shear modes share the same frequency trends. Therefore, it is the stacking order that determines the Raman intensities of the interlayer shear modes. The symmetry determines the Raman activity of each mode, but not the Raman intensities.
In what follows, we present a bond polarizability model that uses information about the relative atomic coordinates and atomic displacements of the vibration modes, without explicitly considering group symmetries. This model provides an intuitive explanation for the stacking-order dependent Raman intensities and frequency trends.  Bond Polarizability Model. In the first principles calculations, the nonresonant Raman intensity of a phonon mode k is computed in the Placzek approximation 26 : where η and η′ are the unit vectors for the polarization of the incident and scattered light,  R k is the second rank Raman tensor, ω k and ω 1 are the frequency and the Boltzmann distribution function of phonon mode k, respectively.
is the derivative of the electronic polarizability tensor αβ P with respect to the atomic displacement. γ u l is the displacement of lth atom in the γ direction for normal mode k and χ γ l k is the γth (x, y, or z) compononent of the eigenvector of phonon mode k. Here we consider the backscattering geometry with the polarization of η and η′ parallel to the in-plane a axis, i.e. ( ) z xx z in Porto notations 27 . From equation (1), one can deduce that for the same frequency and occupation numbers, the Raman intensity is proportional to the change in polarizability of the system, when the atoms are displaced from equilibrium in the direction of the phonon eigenvector. The polarizability is defined as the induced dipole moment relative to the applied electic field. For the ( ) z xx z configuration, the relavant polarizability component would be the xx component, i.e. the dipole moment in the x direction, induced by an applied electric field in the x direction.
As the atoms are being displaced and bonds stretched/compressed, it is reasonable to consider that the major change in polarizability will arise from changes in the relevant dipole moments of the bonds. This has been quantified in an empirical bond polarizability model 28 which can quantitatively predict the Raman intensities of fullerene 28 and graphene ribbons 29 . In this approach, the polarizability is written as a sum of individual bond polarizabilities, which are assumed to be roughly independent of the chemical environment 29 : is the bond vector connecting atom l to one of its nearest neighbor atoms l′ connected by bond B, the vector being normalized to unity. α and α ⊥ are the static longitudinal and perpendicular bond polarizability, respectively, which are further assumed to only depend on the bond length R. We note that equation (2) has further assumed cylindrical symmetry around the principal axis of each bondalthough this may not be exactly correct in general systems, in this work, we are only concerned with relative magnitudes of the Raman intensities, for which these details become unimportant (see later discussion). The derivative of the bond polarizability αβ, P k , which determines the Raman intensity, can be written as: is the bond vector at equilibrium configuration, normalized to unity, R 0β is the β component of ( , ) l B R 0 , while R 0 is the bond length at equilibrium. α′ ⊥ and α ′ are the radial derivatives of the bond polarizability with respect to the bond length. For the ( ) z xx z configuration, we have: For the sake of clarity in explanation, we note that the RHS of equation (4) has five terms, which we will analyse in detail later.
Here, we apply the above bond polarizability model to the interlayer vibrations in layered materials. The justification for applying the bond polarizability model to the interlayer modes is shown in Fig. 3, in which we find distinct charge accumulation in the interlayer regions along the axes connecting nearest neighbouring atoms, for bilayer graphene, MoS 2 and Bi 2 Se 3 . This indicates that although the interlayer interactions are widely believed to be of the van der Waals type, they also have some covalent character. The charge accumulated in these weak covalent bonds can produce dipole moments in the presence of an applied electric field imposed by incident radiation, resulting in Raman intensities which depend on the bond direction, i.e. the stacking order. It is because of this small amount of covalency that experimentally, different interlayer modes are observed for differently stacked materials -if there is no covalency at all, but rather completely delocalized interlayer charge distributions, the stacking order would not have such a significant impact on the Raman intensities.
Before applying the bond polarizability model in full detail, we note that we can understand quite simply the effect of stacking order on Raman intensities of the interlayer modes by using the concept of interlayer bond polarizabilities. We shall focus on the shear mode which is of interest here. For systems with in-plane isotropy, such as graphene, we can without loss of generality consider the shear modes with layers moving in the x direction. When two adjacent layers move against each other, we consider how the x component of the dipole moments along all nearest neighboring atoms would change. For AB and ABC stacked materials such as graphene, TMDs and Bi 2 Se 3 , each atom will have one bond B* with the largest x component, and this will determine the effective direction of the x component dipole moment induced by the field. In AB stacked materials, B* is pointing in opposite directions as we move from layer to layer. However, in ABC stacked materials, B* is pointing in the same direction moving from layer to layer. Therefore to maximize the change in dipole moment (for the largest Raman intensity), adjacent layers should be moving completely out-of-phase for AB stacked systems, while the whole system should be stretched out in the x direction (like a deck of cards) for ABC stacked systems. This suggests that S N-1 and S 1 should have the largest Raman intensities in AB and ABC stacked systems respectively, consistent with the above discussions.
We now illustrate the above qualitative picture more rigorously using equation (4) for AB and ABC stacked three layer graphene (3LG) (Fig. 4). We shall sum over all interlayer bonds with bond length not larger than R for the atoms in the unit cell (the intralayer bonds will not contribute to the change in polarization since the atoms within each layer move in phase). The largest x component of the interlayer bond is a = R*sinλ (see Fig. 4). We note that α ⊥ , α , α′ ⊥ and α ′ are all constants. Considering the interlayer shear mode with layers moving in the x direction, only the x component of χ l k is non-zero. We can now see that terms I, III and IV are proportional to ∑ R lB x 0 which is equal to zero from symmetry. On the other hand, terms II and V are proportional to ∑ R lB x 0 3 which is non-zero. For example, by summing the x component of the interlayer bonds connected to atom A in the AB stacked 3LG shown in Fig. 4a  3 will be zero, giving zero contribution for all terms in equation (4). We plot in Fig. 4e-f the layer by layer top view of the interlayer bonds that will give non-zero contributions to equation (4). For AB stacked 3LG, the non-zero ∑ R B x 0 3 term is equal to j, − 2j and j for atoms A, D and E respectively, while for ABC stacked by χ lx k for each atom l in the unit cell and sum over l. As shown in Fig. 5a, the atomic displacements of the interlayer shear modes are the same for AB and ABC stacked 3LG. Setting χ lx k to be Δ 1 , Δ 2 and Δ 3 for atoms in the 1 st , 2 nd and 3 rd layers respectively, we obtain: for ABstacked 3LG, and for ABC stacked 3LG. For the S 1 mode (lowest frequency) in 3LG, Δ 1 = − Δ 3 , and Δ 2 = 0. From equation (5), and therefore the Raman intensity is zero in AB stacked 3LG, while from equation (6), for ABC stacked 3LG. Thus, from this model, we can explain why, although the lowest frequency modes are Raman active in both AB and ABC stacked FLG, only the lowest frequency shear mode in ABC stacked FLG has non-zero Raman intensity. A similar analysis applies to the highest frequency S N-1 mode. The Raman intensities estimated from the bond polarizability model compare well with LDA calculated intensities for the interlayer shear modes of 2-7 L graphene and bulk graphite (Supplementary Figure S2). In AB stacked bulk graphite, there is a Raman peak located around 44 cm −1 corresponding to the E g 2 mode, but no Raman peaks are found in ABC stacked bulk graphite in the low frequency range. In the limit of ABC stacked bulk graphite, we find that the individual bond polarizabilities for the interlayer modes will always cancel out when the periodic boundary condition is considered, thus giving zero Raman intensity.
For the breathing modes, term IV in equation (4) is zero, but in general the other terms are non-zero. However, since differences in stacking order are reflected in differences in the in-plane components of ( , ) l B R 0 , while it is the z component of ( , ) l B R 0 that shows up in equation (4), the stacking order has no influence on the relative Raman intensities of the breathing modes, consistent with our DFT calculations in the previous section. Generalization to other 2D materials. Besides graphene multilayer systems, similar stacking dependent interlayer vibration modes are expected to exist in other 2D layered materials, such as MoSe 2 , MoS 2 , WSe 2 and Bi 2 Se 3 . The building block of TMD is made up of three atomic layers (trilayer), with the transition metal atom covalently bonded to the chalcogen atom within the trilayer, while the building block of Bi 2 Se 3 consists of "quintuple layers". In the interlayer vibration modes, the atoms within the building block are moving in phase with similar displacements, so the interlayer modes represent the relative displacements between the building blocks separated by the interlayer gap. Similar to graphite, the most stable structure for MoSe 2 , MoS 2 and WSe 2 is the 2H (AB) stacked sequence, with the rhombohedral (ABC) stacked system being a metastable structure. However, for Bi 2 Se 3 , the stable structure has the rhombohedral (ABC) stacked order. The symmetry of AB stacked TMD is the same as that for AB stacked FLG, while both ABC stacked multilayer Bi 2 Se 3 and ABC stacked FLG belong to the symmorphic space group D d 3 3 (P m 3 1) with inversion symmetry. In contrast, ABC stacked TMD belongs to the point group C v 3 without inversion symmetry, with the irreducible representations of zone center phonon modes given by ( + ) N E A 2 1 . We note that in both TMDs and Bi 2 Se 3 , the interlayer bonds are the same as those represented by the A atom in Fig. 4 (except for the vertical bond between A and the atom below it). For example, each chalcogen atom has three nearest neighbouring chalcogen atoms in the adjacent layer, and from the top view, it sits in the middle of the triangle formed by these three neighbours (Fig. 6a,b). Therefore the above bond polarizability analysis for FLG also applies to these systems. Indeed, according to our DFT LDA calculations for AB and ABC stacked TMDs as well as for Bi 2 Se 3 (ABC stacking), the shear modes with largest Raman intensity are S N-1 and S 1 for AB and ABC stacked systems respectively (Fig. 7). Detailed LDA results are shown in Supplementary Tables S3 to S10. Experimental evidence. Next, we turn to experimental evidence of the predicted trends. Recently, the group of Rui He has reported the Raman spectra of suspended AB and ABC stacked trilayer graphene (3LG) 30 . At the ultralow frequency range, they observed an interlayer breathing mode (59 cm −1 ) and highest frequency shear mode (37 cm −1 ) in AB stacked three-layer graphene (3LG), and the same breathing mode (59 cm −1 ) and the lowest frequency shear mode (22 cm −1 ) in ABC stacked 3LG. This observation is in good agreement with our analysis for the 3LG system. Our experiments on 5LG also found that the highest frequency shear mode (42 cm −1 ) only appears in AB stacked 5LG but disappears in ABC stacked 5LG (see Supplementary Figure S3). For the other 2D materials, besides the experimental results already published for AB stacked MoS 2 and WSe 2 , and for ABC stacked Bi 2 Se 3 (shown in Fig. 7), we have measured the ultralow Raman spectra for AB and ABC stacked 3L MoSe 2 (Fig. 8).
In Fig. 8, we show results from Raman scattering experiments on AB and ABC stacked MoSe 2 samples grown from chemical vapor deposition (CVD) 31 . Figure 8a-c shows the optical contrast and atomic force microscopy (AFM) image of a 3L MoSe 2 sample; the measured Raman spectra on different 3L MoSe 2 samples are shown in Fig. 8d. Interestingly, 35% of all 43 samples exhibit only a peak around 13.3 cm −1 , very close to the theoretical prediction of 14.4 cm −1 for the lowest frequency shear mode S 1 in 3L ABC stacked MoSe 2 samples, while 28% of all samples exhibit a peak at 23.1 cm −1 , close to our theoretically predicted value of 25.5 cm −1 for the highest frequency shear mode S N-1 in 3L AB stacked MoSe 2 . Since exfoliated 2H MoSe 2 has the AB stacking order, we compare our results with Raman spectra on the exfoliated sample. The Raman peak of the exfoliated AB stacked sample is around 23.3 cm −1 , confirming that the samples that show a peak around 23.1 cm −1 are AB stacked. It is noted that 37% of all samples detect both S 1 and S N-1 modes, because the spot size of our laser is at least 1μ m, and the 2H and 3R phases can coexist in the same sample with either a sharp boundary 32 , or a few-hundred nanometer wide transition area 31 . Furthermore, scanning tunneling electron microscopy (STEM) experiments performed on the CVD-grown samples 31,32 can be used to distinguish the AB and ABC 3L stacked samples by comparing the STEM images with simulated images 31 . Experimental data for MoSe 2 thin films with different number of layers is plotted in Fig. 7c; the results agree well with the theoretical predictions.

Discussion
The bond polarizability model is based on the formula for nonresonant Raman scattering intensities, which are typically computed within first principles density functional perturbation theory within the Placzek approximation. The computed non-resonant Raman intensities match reasonably well with intensities in experimental ultralow frequency Raman spectra obtained for ABA and ABC stacked materials with typical laser wavelengths of 532 nm 12,14,30,31,33 . Recently, resonant Raman scattering of interlayer vibration mode has been observed in the twisted few layer graphene system [34][35][36][37] . It is noted that the symmetry of the system is much lowered because of the twist, so that many more modes will become Raman active. On the other hand, the twist between two graphene layers results in an overlap of the two Dirac cones in a small region of k-space, leading to Van Hove singularities in the joint density of states, which allow for more optically allowed transitions 35,38,39 . Due to electron-phonon coupling, Raman intensities are enhanced when the laser excitation wavelength matches those of the optical transitions 38 . While  resonant Raman scattering reveals new insights into 2D layered materials 16 , our analysis and conclusions are still applicable in the majority of ultralow frequency non-resonant Raman experiments, and in particular for simply AB and ABC stacked layered materials [12][13][14]21,23,24,30,31,33,[40][41][42] , as we have discussed above.
The advantage of a bond polarizability model is that it can be efficiently applied to more complicated stacking sequences in thicker samples, such as the ABACA stacking in 5 layer MoSe 2 . We show these predictions in Fig. 5c. As the ABACA stacking order could not be distinguished in STEM, the predicted Raman intensities along with experimental Raman spectra were used to assign the stacking sequence in a few of our CVD-grown 5 layer MoSe 2 samples 31 . We also note that the bond polarizability predictions are consistent with our LDA results for these complicated stacking sequences. The above results provide evidence that although the space groups are not all the same for materials with the same stacking order, the frequency evolution trends as predicted by our bond polarizability model are general. Furthermore, these predictions can be used to determine the stacking orders in experimental samples, even when the stacking order cannot be distinguished by microscopy.
Finally, we note that in few layer black phosphorus (Fig. 6d), the in-plane shear modes cannot be detected under the ( ) z xx z polarization configuration. Considering both shear modes in the x and y directions, we can see that any change in bond polarizability will cancel out when we sum over all nearest neighbor bonds, because of the mirror planes along the y and x directions respectively. On the other hand, the breathing modes can have non-zero Raman intensity in the ( ) z xx z configuration. These predictions are consistent with first principles calculations and Raman spectroscopy experiments 40 .
In conclusion, we have shown that the interlayer bond polarizabilities can allow us to understand, both intuitively and semi-quantitatively, the Raman intensities of interlayer modes in general 2D layered materials. Specifically we find that the change in polarizability is maximized for the lowest frequency shear mode in ABC stacked materials, but for the highest frequency shear mode in AB stacked materials, regardless of the details of the space group. The resultant differences in Raman intensity result in clear and distinct frequency trends for AB and ABC stacked systems -as the number of layers increases, the interlayer shear mode red shifts for ABC stacked systems, and blue shifts for AB stacked systems. Because these trends are distinct, and furthermore, do not overlap, they provide a general way to distinguish AB and ABC stacking in 2D layered materials. This bond polarizability model also provides consistent results when applied to interlayer modes in other materials such as black phosphorus, and can be used as a tool to make quick predictions on Raman intensities for more complicated stacking orders. Our predictions are substantiated with extensive first principles calculations as well as Raman spectroscopy measurements on different 2D materials.

Methods
First principles calculations of vibrational Raman spectra are performed within density-functional perturbation theory (DFPT) as implemented in the plane-wave code QUANTUM-ESPRESSO 43 . The local density approximation (LDA) 44 to the exchange-correlation functional with projector-augmented wave potentials is employed for the calculation of phonon frequencies, while the non-resonant Raman intensity is calculated using the norm-conserving pseudopotential, within the Placzek approximation 26 . With the frequency and Raman intensity, a scale parameter of 1 is used in the Lorentzian broadening function to get the calculated Raman spectrum. To get the converged results, a plane-wave kinetic energy cutoff of 65 Ry is used for the wave functions, and the convergence threshold is set to 10 −9 eV and 10 −18 eV in the electron and phonon self-consistent calculation, respectively. The structures are considered as relaxed when the maximum component of the Hellmann-Feynman force acting on each atom is less than 0.003 eV/Å. A Monkhorst-Pack k-point mesh of 44 × 44 × 1, 17 × 17 × 1 and 11 × 11 × 1 are used to sample the Brillouin Zones for the FLG, TMD and Bi 2 Se 3 systems, respectively. We use a vacuum thickness of 16 Å in the direction perpendicular to the slabs to prevent interactions between periodic slab images. The spin-orbit coupling effect is included self-consistently by using fully relativistic pseudopotentials for the valence electrons in Bi 2 Se 3 .
Experimentally, the graphene layers and MoSe 2 are prepared by the mechanical exfoliation method 45 and chemical vapor deposition method 32 , respectively. The Raman spectroscopy measurements are conducted in a backscattering geometry, excited with a Helium-Neon laser with λ = 532 nm for MoSe 2 and graphene layers. The detection of ultralow frequency is achieved by filtering out the laser side bands through the adoption of a reflecting Bragg grating, the ultralow frequency of MoSe 2 is achieved by triple-grating setup (Horiba-JY T64000). The laser power is kept below 0.05mW on the sample surface to avoid laser-induced heating.