In 2004, Geim and Novoselov successfully exfoliated monolayer graphene, a true two-dimensional (2D) materials^[1], which stimulated intense interest on other 2D materials (2DMs)^{[2, 3]}, such as hexagon boron-nitrogen (hBN), transition metal dichalcogenides (TMDs) (Fig. 1). The physical properties of 2DMs exhibit a very broad range from superconductor, semimetal, semiconductor to insulator and show obvious difference from their bulk counterpart, such as thickness-dependent band structures^{[4, 5]} and lattice vibrations^[6]. 2DMs can be stacked layer by layer in precise sequence to form van der Waals heterostructures (vdWHs) without limitation of lattice match (Fig. 1), combining the desirable properties of individual constituents into a single structure. Two flakes with the same monolayer counterpart can be stacked with each other with a certain twist angle , forming the so-called twisted bilayer or multilayer 2DMs and generating moiré patterns. Physical properties of twisted bilayer graphenes and twisted multilayer graphenes (tMLG) had been intensely studied^[7–9]. The properties in vdWHs can be tunable by the choice of constituents, structure symmetry, interface and thickness, due to the modification of quasiparticles, including electron, exciton and phonon. This provides a plethora of opportunities for observation of numerous exciting physical phenomena and quantum design^[10-12]. For example, moiré phonons are found and modulated by moiré potential in twisted bilayer MoS ( 2LM)^[13], interlayer exciton in vdWHs exhibits feature of moiré exciton with modulation of interfacial moiré potential^[14], which are observed in several experiments recently^[15–17].

The intrinsic or modulated properties of 2DMs and vdWHs has been characterized by a large variety of spectroscopic methods^{[6, 18]}. Among them, Raman spectroscopy, one of the most important nondestructive, fast and relatively inexpensive characterization tools, gives the maximum structural and electronic information with high spatial resolution at both laboratory and mass-production scales^{[9, 19]}. Phonon, i.e. elementary excitation of lattice vibration, is involved in the Raman scattering process. Raman spectra can provide rich information, such as crystal quality, lattice structure, disorder, strain, doping and interlayer coupling of 2DMs^[20]. Due to the momentum and energy conversion rules, the involved phonon in first-order Raman process locates at the point of the Brillouin zone (BZ), and its frequency can be determined by the energy difference between incident and scattered photons. Lattice vibrations in 2DMs and vdWHs can be classified into two types, intralayer and interlayer modes, which mainly stem from the intralayer chemical bonds and interlayer var dan Waals forces, respectively^[21]. In general, due to wake nature layer-layer interaction, interlayer modes can be observed below 100 cm^–1, close to strong Rayleigh lines. The recent upgrades of low frequency Raman scattering technique improves detection of low frequency Raman signals, by employing notch filter based on the volume Bragg grating technique, in a single monochromator Raman system^[22]. Thanks to the improved technique, investigation on interlayer Raman modes recently becomes increasingly popular^[21].

The interlayer modes in 2DMs and vdWHs, corresponding to layer-layer vibrations, are sensitive to interlayer coupling, which can be well described by linear chain mode (LCM)^{[8, 19, 22, 23]}. Based on LCM, each layer of an -layer 2DM can be treated as a single ball so that the 2DM can be simplified as a linear chain with balls in which only nearest-neighbor interlayer interaction is considered. According to the relative motions of the planes themselves either perpendicular or parallel to their normal, the interlayer modes can be denoted as shear (S) and layer breathing (LB) modes, respectively. Here, for the isotropic 2DMs, there are –1 degenerate pairs of S modes and –1 LB modes, denoted as S and LB ( = 1, 2, , –1), where S (LB ) (i.e., = –1) has the highest frequency and S (LB ) (i.e., = 1) has the lowest frequency. Assuming that the nearest-neighbor force constants per unit area for the S and LB modes are and , respectively, the frequencies (in cm^–1) of the S ( (S)) or LB ( (LB)) modes can be calculated by solving the corresponding (tridiagonal) dynamical matrix as follows^[22]:

$ \left(\omega_{i}\right)^{2} { u}_{i} = \dfrac{1}{2 {\text{π}}^{2} c^{2} \mu} { D} { u}_{i}, $ (1)

where the is the phonon eigenvector with mode and frequency , is the monolayer mass per unit area, cm/s is the speed of the light, and D is the shear or LB part of the force constant matrix. Thus, the frequency of the S and LB modes can be given by^{[8, 19]}:

$ \begin{split} & \omega\left(\mathrm{S}_{N,N-i}\right) = \dfrac{1}{{\text{π}} \mathrm{c}} \sqrt{\alpha_{0}^{\scriptstyle\text{‖}} / \mu} \sin (i {\text{π}} / 2 N), \\ & \omega\left(\mathrm{LB}_{N,N-i}\right) = \dfrac{1}{{\text{π}} \mathrm{c}} \sqrt{\alpha_{0}^{\perp} / \mu} \sin (i {\text{π}} / 2 N). \end{split} $ (2)

And the th displacement eigenvector for th layer is

$ v_{j}^{(i)} = \cos [i(2 j-1) {\text{π}} / 2 N]. $ (3)

For bulk counterpart of 2DM flakes, , and . Thus, Eq. (2) can be simplified as follows^{[9, 19]}:

$\begin{split} & \omega({\rm S}_{N,N-i}) = \omega({\rm S}_{\rm{bulk}}) {\rm{sin}}(i {\text π}/ 2 N), \\ & \omega({\rm{LB}}_{N,N-i}) = \omega({\rm{LB}}_ {\rm{bulk}}) {\rm{sin}} (i {\text{π}}/ 2 N) . \end{split} $ (4)

Several branches of the S and LB modes have been observed in 2H-MoS₂ flakes^[23]. and show significant dependence on number of layers of MoS₂ flakes, which can be well explained by LCM in which only nearest-neighbor interlayer interaction is considered. LCM can be used to understand the thickness-dependent and of other 2DMs^{[19, 6, 9]}, and be extended to probe the interfacial coupling in vdWHs^[24]. Furthermore, the frequency of the intralayer modes in multilayer and twisted 2DMs can also be modified by the interlayer coupling, as reported for the Davydov splitting^{[25, 26]} and moiré phonons in 2LM^[13], respectively.

Here, we briefly give an up-to-date overview on the lattice vibrations in vdWHs studied by Raman spectroscopy. First, we discuss interlayer coupling in twisted multilayer graphene (tMLG). Then, twist-angle ( ) dependent moiré phonons in t2LM are discussed. Furthermore, LCM and interlayer bond polarizability model are extended to understand the frequencies and relative intensities of layer-breathing (LB) modes in MoS /graphene and hBN/WS heterostructures, respectively.

We first discuss vdWHs formed by 2DM flakes with the same monolayer counterpart, i.e., the so-called twisted multilayer 2DMs. We show that the twist in twisted multilayer 2DMs modify the lattice symmetry from their constitutes and make more interlayer vibration modes be observed. The different interfacial coupling between the S and LB modes in twisted multilayer 2DMs make the atomic displacement of the S modes be localized in the constitutes and that of the LB modes be extended over the whole twisted multilayer 2DMs. The moiré patterns can be formed in twisted multilayer 2DMs to create moiré superlattices stemming from periodic interlayer interaction potentials. The moiré superlattices can result in the zone folding effect for the phonon dispersion curve of each constituent, which can be used to probe the phonon dispersion curves of monolayer constituents.

In AB-stacked NLG (AB-NLG), only the branch of the shear modes with the highest frequency has been observed^[22]. The shear mode in NLG is also referred to as the C mode because it provides a direct measurement of the interlayer Coupling^[22]. The N-dependent can be well fitted by Eq. (2) and N/m³. However, the LB modes are absent in the Raman spectra of AB-NLG because of their Raman inactivity or weak electron-phonon coupling^[22]. Besides the AB and ABC stacking configurations, NLG can also be formed by assembling mLG (AB-mLG if ) and nLG (AB-nLG if ) flakes together with a twist angle ( ) at the twisted interface, which is denoted as tNLG or t(m+n)LG^{[7, 8]}. Due to the zone folding effect resulting from the moiré superlattices in tNLG, LO and TO phonons of mLG and nLG can be folded back to the point and become Raman active, making so-called R and R' modes be observed in Raman spectra, respectively. The frequencies of R and R' modes can be used to identify in tNLG. In tNLG, moiré superlattice induces periodically modulated electron interaction between different constituents, which results in -dependent van Hove singularities (VHSs) in tNLG. By changing , it is possible to tune the optical absorption of tNLG.

As presented in Fig. 2(a), optical contrast of t LG exhibits additional optical absorption features around 2.0 eV in comparison to that of AB-stacked NLG, which corresponds to the energies of VHSs in the optically-allowed joint density of states. When the excitation energy (E) is close to energies of the VHSs, all the Raman modes can be resonantly enhanced, and the C and LB modes can be observed in the Raman spectra. As shown in Fig. 2(b), the C modes of 3LG and the LB modes of 4LG are observed in t(1+3)LG, where the LB modes are absent in cross (XY) polarization configuration due to the Raman selection rule.

To reproduce the interlayer modes in t LG, the modified LCM is presented in Fig. 2(c), as an example for t(2+3)LG. The shear coupling strength (i.e., force constant) at the twisted interface ( ) is found to be about 20% of the bulk case ( N/m³), and the force constant between the layers adjacent to the interface is about 90% of . The twisted interface obstructs the interfacial shear coupling between the constituents of tNLG so that only the C modes of the AB-stacked constituents can be observed in Raman spectra of t LG. For the LB modes, the experimental results cannot be well reproduced if one only take into account in the LCM, and second-nearest-neighbor layer-breathing force constant ( ) is necessary to be included into the LCM. The experimental results indicate that interface twist exhibits little influence on the layer-breathing coupling strength in t LG, and (LB) is only determined by the total number ( ) of layers of tNLG. Thus, the interfacial LB coupling should be almost identical to the interlayer LB coupling in AB-NLG, and the interlayer LB force constants of AB-stacked NLG can be obtained from the experimental (LB) of tNLG, i.e., N/m³ and N/m³. Based on the above discussions and the deduced force constants in tNLG, one can conclude that, under the resonant excitation conditions, the layer number of tNLG and its AB-stacked constituents can be identified according to the observed (C) and (LB), respectively. This method had been applied to identify the stacking order and total number of layers of multilayer graphene grown by chemical vapor deposition^{[27, 28]}.

The lattice structure of twisted bilayer 2DMs can be rigorously periodic to form the crystallographic superlattice, where the superlattice vectors of the two monolayer constituents match with each other, giving a finite unit cell. The lattice structure of the crystallographic superlattice can be represented by a pair of mutual prime numbers ( ) ( ). In addition to the crystallographic superlattices, the periodic moiré patterns are also generated in twisted bilayer 2DMs and can be also denoted by ( ). The lattice constant of moiré pattern in twisted bilayer MoS (t2LM) or twisted bilayer graphene ( 2LG) can be given , and its corresponding value of moiré reciprocal lattice can be evaluated as , where is the in-plane lattice constant of monolayer constitute and . For a given twisted bilayer 2DM, and are only dependent on . However, those in crytallographic superlattice are also dependent on the index of ( , ), i.e., and . In general, the period of moiré superlattices can be equal to or smaller than that of the corresponding crytallographic superlattice, while both of them might induce phonon folding effect in the Raman spectra^[13]. The BZ of crystallographic (solid green lines) and moiré superlattices (dashed red lines) are presented in the left panel of Fig. 3(a), taking t2LM with = as an example. Raman measurement is applied to twisted bilayer MoS to clarify moiré superlattices or crystallographic superlattices is responsible for the zone folding effect in twisted bilayer 2DMs^[13].

Raman spectra in the low frequency region of t2LM, 2H- and 3R-stacked MoS are presented in Fig. 3(b), where the LB mode provides a direct evidence for the good interfacial coupling between the constituents. The 2H stacking at and 3R stacking at are the first and second most stable stacking for bilayer MoS (2LM), and thus t2LMs tend to exhibit 3R and 2H stacking when its twist angle is and , respectively. Fig. 3(b) shows that the S mode in the t2LM with exhibits weaker intensity relative to the corresponding LB mode but similar frequency with respect to that in the 2LM with , in agreement with theoretical prediction.

Due to the periodic potential in crystallographic and moiré superlattices, the off-center phonons in the monolayer constituent linked with the corresponding reciprocal lattice vectors can be folded back to the point of the superlattices. These folded phonons may become Raman active in twisted bilayer 2DMs. Indeed, as shown in Fig. 3(c), many new Raman modes are observed in t2LMs and whose peak positions vary monotonously in the ranges of – and – . This is consistent with -dependent . However, it is not such case for the crystallographic superlattices in which is also sensitive to . Therefore, phonon folding effect of 2LMs results from moiré potential and links with the basic vector of moiré reciprocal space. The folded phonons are referred as moiré phonons. In Fig. 3(c), the seven series of -dependent Raman modes are identified as Raman scattering from moiré phonons associated with out-of-plane acoustic (ZA), transverse acoustic (TA), longitudinal acoustic (LA), (TO₁), (LO₁), (TO ) and A phonon branches in 1LM, respectively, where and are respectively related to longitudinal optical (LO₁) and transverse optical (TO₁), LO₂ and TO₂ phonons branches. Fig. 3(d) presents the frequencies of moiré phonons as a function of . According to , the phonon branches of monolayer constituent in BZ of monolayer constituent can be obtained, as shown in Fig. 3(e).

Because the phonon folding effect related to moiré phonons is mastered by the basic vector of moiré reciprocal lattices, rather than by that in the reciprocal lattices of crystallographic superlattices of t2LM, this study can be extended to other twisted bilayer 2DMs and the related vdWHs.

The interfacial coupling in vdWHs between different 2DM constitutes with distinct properties is also important for their fundamental research and application. Here, we discuss two vdWHs formed by different 2DMs, such as semimetal graphene and semiconductor MoS₂, and insulator hBN and semiconductor WS₂. It shows that interfacial coupling in vdWHs plays significant role in their vibration properties.

Fig. 4(a) presents Raman spectra of vdWHs formed from bilayer MoS₂ (2LM) and nLG, denoted as 2LM/nLG, where three branches of the LB modes (LB , LB and LB ) are observed. The red-shift of their peak positions with increasing n of nLG indicates the existence of significant interfacial coupling between the two constituents although 2LM and nLG are quite different 2DMs, semiconductor and semimetal. In general, only nearest LB force constant of MoS₂, N/m³ can be considered in the LCM to well reproduce the thickness-dependent (LB) in multilayer MoS ^[23]. However, in nLG, besides the nearest LB force constant , the second nearest LB force constant is necessary to be considered. In vdWHs, graphene and MoS should be considered as an overall system, and the nearest interfacial LB coupling ( N/m³) between MoS₂ and graphene constituents are also necessary, as shown in Fig. 4(b). The modified LCM in Fig. 4(b) can well reproduce all the experimental (LB) as a function of layer number of MoS₂ and graphene constituents in vdWHs, as shown in Fig. 4(c). This result suggests that the interfacial layer-breathing couplings in the vdWHs formed by MoS₂ and graphene flakes are not sensitive to their stacking order and twist angle between the two constituents.

The S modes of 2LM/nLG vdWHs show almost identical peak position to that of 2LM, which indicates that the observed S mode in 2LM/nLG vdWHs is mainly localized in the 2LM constituent. The weaken interfacial shear coupling had also been observed in other mLM/nLG vdWHs and tMLG. The mismatched lattices between graphene and MoS layers can cause the local interfacial shear coupling either attractive or repulsive, which makes the total shear coupling weakened.

Interfacial layer-breathing coupling in vdWHs formed by semimetals and semiconductors induces new LB modes, which not only can be used to measure the interfacial interactions in vdWHs but also is beneficial to identify the total thickness of vdWHs. However, the S modes are localized in the constituents and thus they can be used to determine the numner of layers of the constituents in vdWHs.

The intralayer Raman modes in graphene/MoS vdWHs only consist of those of the two constituents. The peak position of the MoS constituent can be slightly modified by the interfacial coupling between the two constituents.

The LCM demonstrated by the examples of MLG and MoS /graphene vdWHs can also be generally applied to other vdWHs, such as that formed by n-layer hBN and m-layer WS₂ (nL-hBN/mLW) vdWHs. The LB modes from the collective vibrations of all the stacking layers in vdWHs were also reported^[24] even though the vdWHs are over tens and hundreds of layer thickness, as shown in Fig. 5(a). According to the interlayer LB coupling of WS₂ ( N/m^3[30] and the experimental frequencies of the LB modes in nL-hBN/mLW vdWHs, the interlayer LB coupling of hBN ( N/m³) and the interfacial coupling between nL-hBN and mLW ( N/m³) can be obtained. The stacking-order independent interfacial LB coupling is comparable to the interlayer coupling in WS₂ and hBN constituents, confirming that the LB phonons are from the collective vibrations of all the stacking layers in vdWHs with bulk-like features. Furthermore, these LB modes can be significantly resonantly enhanced when the excitation energy approaches the C exciton energy of standalone mLW flakes, suggesting these three-dimensional LB modes can strongly couple to the 2D electron states confined within the few-layer WS₂ constituents.

To understand this cross-dimensional electron–phonon coupling (EPC), a microscopic picture mediated by the interfacial coupling is proposed to calculate the relative Raman intensity of the LB modes in vdWHs^[24]. Due to the significant interfacial coupling, the LB vibrations in mLW constituent efficiently interact with those in hBN constituent, leading to the bulk-like LB vibrations. The Raman intensity of the LB modes is determined by the EPC strength, which is to the first order, a function of the vibration displacement (Fig. 5(a)). Because E is directly resonant with the C exciton of the mLW constituent, the EPC strength of the LB modes in vdWHs can be estimated by the sum of its weighting factor of interlayer displacements from all the LB modes in the standalone mLW flakes. The weighting factor can be given by the projection between the vdWH LB phonon wavefunction ( ) and that ( ) of LB ( = 1, 2, , m–1) phonon in a standalone mLW flake, i.e., . And the Raman intensity is proportional to , where is the relative Raman intensity of the corresponding LB in the standalone LW flake. The calculated relative Raman intensity of different LB modes in 39L-hBN/3LW is shown in Fig. 5(b), in line with the experimental results. The relative Raman intensity of different LB modes in nL-hBN/mLW vdWHs can also be understood by the interlayer bond polarizability model^[24]. In this case, each layer is simplified as a single object and the Raman intensity is proportional to the square of system's polarizability change ( ), which is related to interlayer bond polarizability and bond vector^[31], i.e., , where is the polarizability derivative of the entire layer i with respective to the direction normal to the basal plane ( ) and is the normal displacement of layer i from the LCM. In nL-hBN/mLW, only the top and bottom layer of hBN and WS constituents should be considered, as shown in Fig. 5(a). Thus, the polarizability change of L-hBN/3LW is , where (W), (BN) and (I) are fitting parameters, related to the properties of the interlayer bond in WS , hBN constituents and that at the interface, respectively. The relative Raman intensity of the LB modes in 39L-hBN/3LW can be well fitted with (I)/ (W) = 0.3 and (BN)/ (I) = 0.003, as shown in Fig. 5(c). Similar analysis can be applied to vdWHs involved different layers^[24].

The models presented above can be generalized and extended to other vdWHs. The cross-dimensional electron-phonon coupling in vdWHs provides a way to manipulate both the designable phonon excitations in vdWHs and their coupling to the electronic states by varying the constituents and engineering the interface.

This review presents application of Raman spectroscopy to investigate the lattice vibrations in twisted bilayer and multilayer 2DMs and vdWHs formed by different constituents. We introduce LCM to understand the interlayer vibration modes in 2DMs, which can also be extended to vdWHs and reproduce their interlayer Raman modes. From the experimental shear and LB modes, the interlayer shear and LB coupling strength in individual components and interface can be obtained, respectively. Furthermore, the interlayer coupling in twisted bilayer 2DMs induces periodic moiré potential, which modifies the lattice vibrations of monolayer constituent to form moiré phonons. By varying the twist angle, the moiré phonons of twisted bilayer 2DMs can be exploited to map the phonon dispersions of the monolayer constituent. Because of the significant interfacial layer-breathing couplings between the two constituents, many new LB modes with frequencies dependent on their layer numbers are observed in MoS₂/graphene and WS₂/hBN vdWHs. Because of the large lattice mismatch between two constituents, the interfacial LB couplings are not sensitive to their stacking order and twist angle in vdWHs. The flexible van der Waals stacking in vdWHs leads to multiple opportunities to engineer the interlayer phonon modes for cross-dimensional electron-phonon coupling.

We acknowledge support from the National Key Research and Development Program of China (Grant No. 2016YFA0301204), the National Natural Science Foundation of China (Grant Nos.11874350 and 11434010).