Signature of lattice dynamics in twisted 2D homo/hetero-bilayers

Few layers of hBN are mechanically exfoliated from bulk material (epitaxial solidification technique grown, purchased from 2D semiconductors) and transferred onto a Si substrate with 300 nm SiO₂ via Nitto tape. Monolayer WSe₂ is exfoliated from bulk 2H phase WSe₂ (from HQ Graphene) onto polydimethylsiloxane (PDMS) by Nitto tape followed by a deterministic all-dry transfer technique [29]. A tear-and-stack method [30] is performed to fabricate tB WSe₂. A large area of a monolayer WSe₂ flake is first exfoliated onto PDMS. Bringing the monolayer WSe₂ flake partially in contact with hBN and then slowly peeling off leaves the monolayer WSe₂ partially on hBN while the rest part still sticks on the PDMS stamp. In the next step a homo-bilayer sample is prepared by rotating the substrate with a defined angle and bringing the remaining part of the WSe₂ flake to contact. Via this procedure, the twisting angle can be well defined. The procedure is shown in detail in figure S1.

Twisted hetero-bilayers are prepared by an advanced tear-and-stack method. After preparing a 30^∘ twisted homo-bilayer, the top layer of different TMDC materials is stamped via PDMS. The edge of the top layer is aligned to any of the edges of the bottom homo-bilayer with an angle of $30^{\circ}\,\times\,n$ (n is any integer number since monolayer TMDCs belong to the D$_{3h}$ point group). This results in two areas of heterostructures, one with a relative twisting angle of 30^∘ and the other with 0^∘ or 60^∘. The procedure is shown in detail in figure S2.

Temperature-dependent Raman and photoluminescence (PL) measurements are conducted in a Linkam THMSEL600 chamber under a 50×, 0.5 NA objective, while room temperature Raman spectra are measured in ambient condition under a 100×, 0.9 NA objective. Temperature-dependent PL and room temperature Raman measurements are performed using a Horiba Xplora Plus equipped with a spectrometer comprising 600 l mm⁻¹ (for PL) and 2400 l mm⁻¹ (for Raman, with 1.6 cm⁻¹ spectral resolution) gratings and an electron multiplying charge coupled device (CCD). A DPSS 532 nm continuous wave (CW) laser source was used to excite the samples with an excitation power of 100 µW measured under the objective. Temperature-dependent Raman spectra are acquired by a Horiba LabRAM HR spectrometer with a 2400 l mm⁻¹ grating and an excitation wavelength of 514.7 nm with a spectral resolution of 0.8 cm⁻¹.

Density functional theory (DFT) calculations are performed using the projector augmented wave approach as implemented in the Vienna ab initio Simulation Package [31–33]. The DFT-D3 method [34] is employed to capture the interlayer van der Waals forces accurately, which is crucial in the bilayer system. PBEsol [35] is used as the approximation functional for the exchange-correlation functional and the cutoff energy is set to be 450 eV. For 0^∘ and 60^∘ twisted systems, the unit cells consist of two triangular prisms. A 20 Å thick vacuum layer is set in the unit cell to simulate the boundary conditions. A k-grid of $16\,\times\,16\,\times\,1$ points is used for these two configurations. A set of commensurate twisted bilayer WSe₂ structures are also calculated with various numbers of atoms, including those with twisting angles of 21.8^∘, 27.8^∘, and 38.2^∘ [24]. A smaller k-grid of $4\,\times\,4\,\times\,1$ points is used for these configurations because of the much smaller reciprocal lattice. All structures are relaxed until all forces on individual atoms are below 0.1 meV Å⁻¹. Each individual monolayer undergoes a reconstruction during this process and this phenomenon is successfully modeled by our first-principles calculation. The phonon modes and their frequencies are determined using finite difference method. Density functional perturbation theory is also used to determine the ion-clamped dielectric constant [36, 37], which is used to determine the Raman tensor of selected phonon modes.

The Raman tensor $\skew{2.4}\tilde{R}$ of a given phonon mode can be calculated by the change of polarizability [38]:

$$\begin{align} \skew{2.4}\tilde{R}_{\alpha \beta} \propto \frac{\Delta \varepsilon^{\infty}_{\alpha\beta}}{\Delta Q(s)} \end{align} \tag{ 1 }$$

where ε stands for the ion-clamped dielectric constant, i.e. the dielectric constant at infinite high frequency, Q(s) stands for the amplitude of the normal mode s, α and β indicate the directions. The derivative can be calculated by freezing in the phonon mode s and calculate the change of the dielectric constant. The nonresonant Raman intensity is usually calculated by the Placzek approximation [39]: $I_\mathrm{Raman}\propto |\vec e_{\mathrm{i}}\cdot\skew{2.4}\tilde{R}\cdot\vec e_{\mathrm{s}}|$. The $\vec{e}_\mathrm{i}$ and $\vec{e}_\mathrm{s}$ are the polarization directions of incident and scattered light. In the experimental setup, the incident light is polarized in-plane, noted as x. Both directions of in-plane scattered light were collected, which results in a total Raman intensity $I\propto I_{xx}\,+\,I_{xy} \propto R_{xx}^2\,+\,R_{xy}^2$.

In bulk WSe₂, there are four Raman-active modes according to group theory analysis (frequency from low to high [40]): $E_{2g}^2$, $E_{1g}$, $E_{2g}^1$ and $A_{1g}$. Among these phonon modes, the inter-layer shearing mode $E_{2g}^2$ has a low frequency of 35 cm⁻¹, because it is dominated by van der Waals interaction [41]. The intensity of $E_{1g}$ is low and hard to measure [40]. Thus, in the following sections, we mainly focus on the last two modes and another interesting mode $B_{2g}$, which becomes Raman active with the reduction of layers. The schematics of these three modes can be found in the insets of figure 1(a): the $E_{2g}$ (short for $E_{2g}^1$) mode is an in-plane vibrational mode, where the W and Se atoms vibrate against each other. The $A_{1g}$ mode is known to be an out-of-plane vibrational mode, where the two Se atoms on the opposite sides of a W atom vibrate against each other while the W atom shows no relative motion. The $B_{2g}$ mode is also an out-of-plane one, for which the Se and W atoms within the same layer vibrate against each other with a 180^∘ phase difference with regard to the vibration in the adjacent layer.

Zoom In Zoom Out Reset image size

Bulk WSe₂ can be mechanically exfoliated into 60^∘ twisted bilayer structure [41], which is known as the intrinsic bilayer structure. During this process, the number of symmetry elements are reduced from 24 to 12, among which there is a mirror plane in the W atoms layer perpendicular to the c-axis [40]. The mirror plane is important because the $B_{2g}$ mode elongation is reversed when this c-axis mirror plane is applied but remains unchanged when the three-fold rotation is applied, neither symmetry operation transforms it like the basis function of yz or xz, which makes it Raman inactive. The reduction in symmetry elements makes the $B_{2g}$ mode Raman active. However, if the number of layers is reduced to one, the c-axis mirror plane symmetry makes the $B_{2g}$ mode inactive in the monolayer. It is noteworthy that the irreducible representations (irreps) for the same phonon mode might change in different structures, but we usually refer to these phonon modes by their irreps in the bulk structures to be consistent. A complete list of irreps of the phonon modes in different structures can be found in table 1.

Table 1. Phonon modes and their irreps, frequencies (in cm⁻¹) in different structures. Most of the phonon modes in this table are Raman-active, with a few exceptions indicated by the label 'inactive'. The space group numbers are noted in the parentheses next to each structure name [42, 43]. The frequencies of some models are given as a range. Details of these structures can be found in the supplementary materials.

	Out-of-plane shearing	Out-of-plane breathing	In-plane shearing
Bulk, 2H (#194)	$B_{2g}$: inactive	$A_{1g}$	$E_{2g}$
0∘-tB, 3R (#187)	A1	A1:	E
	311.2	252.3	248.2
21.8∘-tB (#150)	A1	A1	E
	311.4	251.6	248.2–249.3
27.8∘-tB (#149)	A1	A1	E
	311.4	251.8	248.8–249.8
38.2∘-tB (#149)	A1	A1	E
	311.4	251.5	248.1–249.2
60∘-tB, 2H (#164)	$A_{1g}$	$A_{1g}$	Eg
	310.8	251.1	248.6
Monolayer (#156)	$A_2^{^{\prime\prime}}$: inactive	E'	$A_1^{^{\prime}}$

When the bilayer structure is twisted, our group theory analysis that shows the $B_{2g}$ mode remains active. At most twisting angles besides 0^∘ and 60^∘, there are four symmetry elements: a three-fold axis and three two-fold axis perpendicular to each other. This three-fold axis can be found where two tungsten atoms from different layers overlap in the z-direction. This overlapping certainly happens in a moiré pattern. The $B_{2g}$ mode remains unchanged under these symmetry operations, which makes it behave exactly the same way as the basis functions of $x^2\,+\,y^2$ and z². At 0^∘, this statement holds although multiple other symmetry elements are present.

In table 1, we list the phonon modes by their irreps in different structures, along with their Raman-activity and frequencies calculated from density functional theory. These simulated frequencies match fairly well with our experimental values in the next section, which confirms our mode assignment. Our phonon calculation also shows that, at an arbitrary twisting angle, the phonon modes that are originally not at the Γ-point can fold back and couple with the ones at the Γ-point. This phenomenon is significant for the in-plane $E_{2g}$ mode: it can split into a bunch of peaks that are close to each other, however, this is very hard to be observed experimentally.

Figure 1(a) shows the Raman spectra of monolayer and intrinsic bilayer WSe₂. The most pronounced peak corresponds to the combination of two vibrational modes, $E_{2g}$ and $A_{1g}$ mode, almost degenerate at around 250 cm⁻¹ [40, 41, 44, 45]. The feature at around 260 cm⁻¹ is a second order peak caused by a double resonance effect involving the longitudinal acoustic phonon at the M point in the BZ assigned as 2LA (M) [28, 45]. The most interesting feature is the out-of-plane $B_{2g}$ mode located at around 309 cm⁻¹. According to the group theory analysis, the $B_{2g}$ mode is Raman inactive in both bulk and monolayer, but becomes active in few layers due to the reduction of symmetry elements [44–46].

Figure 1(b) displays Raman spectra for a full period ($0^{\circ}\,\leqslant\,\theta\,\leqslant\,60^{\circ}$) tB WSe₂. Raman intensity fluctuations during different Raman measurements are eliminated by normalizing all spectra with respect to the most intense $E_{2g}/A_{1g}$ peak. With increasing twisting angle, $I_{B_{2g}}$ decreases and almost vanishes at 30^∘, then increases again and reaches a second maximum at 60^∘. To visualize this effect better the intensity ratio $I_{B_{2g}}/I_{{E_{2g}}/{A_{1g}}}$ is plotted as a function of twisting angle in figure 1(c). All spectra are fitted with Voigt functions and the integrated areas are taken as a measure for the intensity (a detailed fitting example can be seen in figure S9 in the SI). The red solid line in figure 1(c) is the calculated moiré superlattice constant (moiré period) for $0^{\circ}\,\leqslant\,\theta\,\leqslant\,30^{\circ}$. Since TMDC monolayers belong to $D_{3h}$ point group, the moiré superlattice constant can be calculated using the formula: $a_{\mathrm{moir}\acute{e}}\,=\,0.5a/\sin(\theta/2)$, when $0^{\circ}\,\leqslant\,\theta\,\leqslant\,30^{\circ}$ and $a_{\mathrm{moir}\acute{e}}\,=\,0.5a/\sin(30^{\circ}\,-\,\theta/2)$, when $30^{\circ}\,\leqslant\,\theta\,\leqslant\,60^{\circ}$ [47], where a is the in-plane lattice constant and θ is the relative twisting angle between layers. As shown in figure 1(c), $I_{B_{2g}}/I_{{E_{2g}}/{A_{1g}}}$ follows the progression of the moiré period with twisting angle very nicely. To further confirm this correlation, Raman spectra of such tB systems are also simulated by theoretical calculations. Unlike the intensity, there are no significant changes in the frequencies of the $B_{2g}$, $A_{1g}$ and $E_{2g}$ modes observable with twisting angle in both experimental and theoretical results. The phonon modes in the twisted structures can be viewed as off-center phonons in the monolayer [48]. Therefore, the process of twisting is similar to mapping the phonon dispersion of the monolayer WSe₂ in the in-plane direction. All three phonon modes considered here are high-frequency ones, i.e. mainly strong and localized valence bonds are determinant for the frequencies. These phonon modes thus reveal weak dispersion in the in-plane direction and consequently the resulting phonon frequencies hardly change upon twisting.

In figure 1(c), the simulated intensity ratio of $I_{B_{2g}}/I_{{E_{2g}}/{A_{1g}}}$ shows a similar trend as the experimental one. This trend further confirms that the periodic change of the intensity ratio is closely related to the geometric arrangement of the atoms. However, compared to the experimental result, the intensity ratio is markedly weaker, which is likely due to the inability of taking the resonant Raman conditions into account theoretically. Experimental evidence indicates that the $B_{2g}$ mode has a strong resonance dependence [28]. However, our DFT simulations only take the polarizability change into account and that leads to a significantly lower intensity of the $B_{2g}$ mode. Therefore, also the ratio is expected to be markedly lower than the ratio observed experimentally. We chose to simulate only six representative twisting angles because a large and strict crystallographic superlattice is necessary to precisely model the tB structures, while for most twisting angles, the crystallographic superlattice is much larger than the moiré superlattice, the calculation of which is beyond the capability of modern density functional theory packages [48]. The complete simulation results of the Raman spectra as a function of different twisting angles can be found in the supporting information.

To explain the interesting behavior of the $I_{B_{2g}}/I_{{E_{2g}}/{A_{1g}}}$ intensity ratio, we propose the hypothesis that the change of intensity ratio is induced by the highly localized stacking structure, which results in different interlayer coupling. A schematic diagram of moiré superlattices with different twisting angles is shown in figure S5. One can clearly identify that with increasing twisting angle the area of the energetically favorable AB stacking decreases in a moiré superlattice as also confirmed in [5]. As seen in figure S5(d), no more AB stacking exists in such superlattice when θ = 30^∘. We thus predict that the lack of AB stacking causes reduced interlayer coupling and the two layers cannot maintain 180^∘ vibrational phase difference. In other words, in a 30^∘ tB WSe₂ the two contributing monolayers vibrate like two independent monolayers. If this is true, then the negligible interlayer coupling at 30^∘ does not only lead to the vanishing of the B$_{2g}$ mode, but also strongly modifies the optical response of the system. To further confirm this hypothesis, we performed a series of temperature-dependent Raman measurements, since the optical response of WSe₂ is highly influenced by the excitonic resonance condition [28]. To probe the excitonic resonance condition, one approach is to record Raman spectra with multiple wavelengths (fixed exciton energy while changing the excitation energy), another approach is to explore the temperature dependence of the Raman spectra (fixed excitation energy while tuning the exciton energy by temperature). Temperature-dependent Raman spectra for iB, monolayer, and 30^∘ tB WSe₂ are shown in figure 2 (for simplicity, selected Raman spectra are shown in the main text, more Raman spectra can be found in figure S6). With increasing temperature, Raman features for iB, monolayer and 30^∘ tB WSe₂ shift towards lower frequency together with a broadening of the full width at half maximum (FWHM), which is in good agreement with what was reported in [49, 50]. The most interesting observation in figure 2(a) is that the Raman intensity of iB WSe₂ decreases drastically with increasing temperature. The signal to noise ratio becomes so low that the Raman features are almost invisible when the temperature is higher than 473 K, whereas for the monolayer and 30^∘ tB WSe₂ in figures 2(b) and (c) the Raman spectra remain intense at high temperature. The $E_{2g}$, $A_{1g}$ vibrational modes in WSe₂ are sensitive to the resonance condition with the A' exciton, where the A' arises from the splitting of the ground and excited states of the A transition (K–K transition in the BZ) due to the intralayer perturbation of the d electron band by the Se p orbitals [27, 51, 52]. However, information about the A' exciton energy in tB WSe₂ systems is still lacking in the literature. We therefore run a temperature-dependent PL measurement to determine the A exciton energy, which is used to derive the A' exciton energy. It is noteworthy that the Stokes shift for WSe₂ is 3–5 meV [27, 53], which is negligible for the following calculation, thus the PL peak position is taken as the exciton transition energy. Temperature-dependent PL spectra of iB, monolayer, and 30^∘ tB WSe₂ are shown in figures 3(a)–(c). The peaks guided by red dashed line correspond to the A exciton PL emission. The peak positions of the A exciton in 30^∘ tB WSe₂ are plotted in figure 3(e) and are fitted to a modified Varshni's equation (for iB and monolayer see figures S7(d) and (e)), which describes the temperature dependence of a semiconductor bandgap [54–57]:

$$\begin{align} E_g(T)\,=\,E_g(0)\,-\,S\left \langle \hbar\omega \right \rangle \left[\mathrm{coth}\left(\frac{\left \langle \hbar\omega \right \rangle}{2k_\mathrm{B}T}\right)\,-\,1~\right] \end{align} \tag{ 2 }$$

where $E_g(0)$ is the transition energy at $T\,=\,0\,\mathrm{K}$, S is a dimensionless constant that describes the strength of electron–phonon coupling and $\left \langle \hbar\omega \right \rangle$ describes the average acoustic phonon energy involved in electron–phonon interactions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For iB WSe₂, the A exciton energy $E_{A,\mathrm{iB}} (298\,\mathrm{K})\,=\,1.64$ eV, $E_{A,\mathrm{iB}} (473\,\mathrm{K})\,=\,1.56$ eV, and the energy difference between 473 K and 298 K $\Delta E_{A,\mathrm{iB}}\,=\,0.08$ eV, the A' exciton energy $E_{A^{^{\prime}},\mathrm{iB}} (298\,\mathrm{K})\,=\,2.28$ eV [27, 28], $E_{A^{^{\prime}},\mathrm{iB}} (473\,\mathrm{K})\,=\,E_{A^{^{\prime}},\mathrm{iB}} (298\,\mathrm{K})\,-\,\Delta E_{A,\mathrm{iB}}\,=\,2.20$ eV. Since strong resonance enhancement only occurs in an excitation energy range approximately ±0.1 eV from the A' exciton energy [28], our excitation energy of $E_\mathrm{laser}\,=\,2.41$ eV is far from the resonance condition at 473 K and results in strongly decreasing intensity in the Raman spectra of figure 2(a). Following the same procedure for monolayer WSe₂, $E_{A^{^{\prime}},\mathrm{1L}} (473\,\mathrm{K})\,=\,2.34$ eV, which still fulfills the resonance condition. Therefore, the Raman spectra remain intense at 473 K in figure 2(b). For 30^∘ tB WSe₂, $E_{A,30^{\circ} \mathrm{tB}} (298\,\mathrm{K})\,=\,1.62$ eV, $E_{A,30^{\circ} \mathrm{tB}} (473\,\mathrm{K})\,=\,1.55$ eV, and $\Delta E_{A,30^{\circ} \mathrm{tB}}\,=\,0.07$ eV. However, the information about the A' exciton energy is still lacking. Therefore, we take the values of the energy difference between the A' and A excitons in both iB and monolayer WSe₂ for the calculation. When taking $\Delta E_{A^{^{\prime}}-A,30^{\circ} \mathrm{tB}}\,=\,\Delta E_{A^{^{\prime}}-A,\mathrm{iB}}\,=\,0.70$ eV, following the previous procedure, $E_{A^{^{\prime}},30^{\circ} \mathrm{tB}} (473\,\mathrm{K})\,=\,2.25$ eV, which is far away from the resonance condition and disagrees with our experimental observation. However, when taking $\Delta E_{A^{^{\prime}}-A,30^{\circ} \mathrm{tB}}\,=\,\Delta E_{A^{^{\prime}}-A,\mathrm{1L}}\,=\,0.74$ eV, which results in $E_{A^{^{\prime}},30^{\circ} \mathrm{tB}} (473\,\mathrm{K})\,=\,2.29$ eV, so that the excitation energy is a little further away from the resonance condition than in monolayer, but is still close enough to cause resonance enhancement. This is in good agreement with our experimental results in figure 2(c), namely that we can still observe the Raman features from 30^∘ tB WSe₂ at 473 K well, but the spectra have relative lower signal-to-noise ratio compared to the monolayer results in figure 2(b). It is therefore evident from the resonant Raman scattering point of view that the Raman spectra of 30^∘ tB WSe₂ resemble those of the monolayer indicating that the two layers behave like two independent monolayers. Moreover, figure 3(d) displays the PL spectra of tB WSe₂ with different twisting angle. One can clearly observe the interlayer transition in iB and 1^∘, 10^∘, 60^∘ tB WSe₂. Yet the interlayer transition seems to be absent in 30^∘ tB WSe₂, which also agrees with our hypothesis that 30^∘ tB WSe₂ has a similar optical response as two independent monolayers.

Furthermore, we carry out similar Raman investigations on twisted hetero-bilayer (thB) TDMCs. Inspired by the tear-and-stack method, we developed an advanced tear-and-stack method, which enables to prepare samples with different crystal orientation (details described in figure S2) on the same substrate. This method offers the possibility to carry out large area spectroscopic mapping and to spatially resolve and study the optical properties of 2D materials.

The optical image of thB MoSe₂/WSe₂ prepared by such method is shown in figure 4(a). The areas marked by red and white dashed lines have twisting angles of 30^∘ and 0^∘/60^∘, respectively. Figures 4(b) and (d) display the Raman intensity ratio map of WSe₂ $I_{B_{2g}}/I_{{E_{2g}}/{A_{1g}}}$ and corresponding Raman spectra normalized with respect to the WSe₂ ${E_{2g}}/{A_{1g}}$ peak. The WSe₂ $B_{2g}$ peak is marked by a dashed line in figure 4(d), which is absent in monolayer and relatively weak in the 30^∘ thB MoSe₂/WSe₂, while becomes very intense in 0^∘/60^∘ thB. The strong contrast can be visualized in the Raman intensity map in figure 4(b). The argument still holds in figure 4(c) when the Raman intensity ratio map of MoSe₂ $I_{B_{2g}}/I_{A_{1g}}$ is plotted and in figure 4(e) when the spectra are normalized with respect to the MoSe₂ $A_{1g}$ peak. Similar experiments were also performed on different combinations of materials such as WSe₂/MoS₂ and MoSe₂/MoS₂. As shown in figure S8, all experiments provide similar results. Our study thus provides quite universal results for all kinds of Mo- and W-based TMDC thB.

Zoom In Zoom Out Reset image size