Low-frequency Raman active modes of twisted bilayer MoS2

We study the low-frequency Raman active modes of twisted bilayer MoS2 for several twist angles using a force-field approach and a parametrized bond polarizability model. We show that twist angles near high-symmetry stacking configurations exhibit stacking frustration that leads to significant buckling of the moiré superlattice. We find that atomic relaxation due to the twist is of prime importance. The periodic displacement of the Mo atoms shows the realization of a soliton network, and in turn, leads to the emergence of a number of frequency modes not seen in the high-symmetry stacking systems. Some of the modes are only seen in the XZ Raman polarization setup while others are seen in the XY setup. The symmetry of the normal modes, and how this affects the Raman tensors is examined in detail.


Introduction
Transition metal dichalcogenides (TMDCs), represented by the formula MX 2 (where M represents a transition metal and X denotes a chalcogen), are well-known examples of van der Waals (vdW) materials.These materials exhibit a wide array of electronic characteristics.For example, some, such as Ta, Pt, Ti, Hf, Zr, Mo, and W -based dichalcogenides, display insulating or semiconducting properties.Other dichalcogenides, such as the V-, Nb-, and Ta-based dichalcogenides, demonstrate metallic or semimetallic behavior.This variance Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. in electronic properties is attributed to the gradual filling of the d-bands, which are non-bonding, by electrons from the transition metal element [1].A notable aspect of vdW materials is how even minimal interactions between layers can significantly influence their electronic properties.This effect is evident when two-layered compounds are compared with their monolayered or bulk counterparts.Bilayer graphene serves as a prime example of this phenomenon [2].In the emerging field of twisted bilayer structures, the possible configurations are limitless, which enables for a rich landscape of electronic properties.Xian et al adopted an ab initio-based approach to characterize the electronic properties of twisted bilayer MoS 2 (tBLM).[3] They studied the collective behavior of tBLMs in the presence of interactions, and characterize an array of different magnetic and orbitally-ordered correlated phases, which may be susceptible to quantum fluctuations giving rise to exotic, purely quantum, states of matter.Furthermore, Zhang et al studied tBLMs' properties by means of an ab initio-based tight-binding model [4].Twist angle governs the crystal symmetry and the development of moiré patterns, grain boundaries, and atomic relaxation due to the non-homogenous crystal field can lead to the emergence of number of physical behaviors, such as Hofstadter's spectra, unconventional superconductivity, moiré excitons, tunneling conductance, nonlinear optics, and structural super-lubricity [5].
The study of phonon properties in twisted bilayer MoS 2 is a fascinating and complex field that merges concepts from materials science, physics, and nanotechnology [6,7].MoS 2 , or molybdenum disulfide, is a semiconductor in the TMDC family, known for its unique electrical and mechanical properties.The manufacturing processes which yield thin film and bilayer crystalline MoS 2 can have defects which have notable impacts on the mechanical and electronic properties.While mechanical exfoliation can produce defect free single layers, there is large interest in scalable manufacturing techniques of MoS 2 thin films.To this end, radio frequency magnetron sputtering has been shown to exhibit tunable electronic confinement in Mo ions via modulating voltage bias [8].Chemical vapor deposition can be used with varying annealing to adjust the layer thickness of substrates [9].However, the higher the annealing temperature, the more likely defects are in the interfacial structure.Alternatively, spin-coated MoS 2 films have also been found to exhibit sulfur vacancy induced dielectric peaks [10].When two layers of MoS 2 are stacked with a twist, intriguing new properties emerge, notably in terms of phonons, which are quasiparticles representing the quantum mechanical description of vibrational motion in solids.This twist introduces a moiré patterna large-scale interference pattern that significantly alters the material's electronic and vibrational properties.The degree of twist-the angle between the two layers-is a crucial factor that determines these properties.One of the first studies to investigate the interlayer coupling in tBLM used photoluminescence spectroscopy [11], revealing a tunability of the interlayer coupling.The photoluminescence intensity ratio of trion and exciton was found to reach its maximum value for twist angles 0 • or 60 • , while for twist angles 30 • or 90 • the situation is the opposite [12].This technique has also been employed to investigate folded MoS 2 bilayers, uncovering the influence of interfacial coupling on lattice vibrations and photon transitions [13].Huang and cowokers also investigated the low-frequency Raman modes of interlayer shear and breathing Raman modes (located below 50 cm −1 ) in tBLMs by Raman spectroscopy and first-principles modeling [14].They demonstrated the sensitivity of the behavior of low-frequency Raman modes for probing the interfacial coupling and environment of tBLMs and potentially other two-dimensional materials and heterostructures [12].In the as-grown/transferred tBLMs, Lin et al showed that the periodic potentials determined by moiré patterns also modify the properties of phonons of its monolayer MoS 2 constituent to generate Raman modes related to moiré phonons [15].Zhu et al proposed a general moiré-templated nanoscale morphology engineering method based on bilayer TMDCs [16].Evidence was also provided on how the morphology can modulate band gap and optical absorption of an MoS 2 monolayer, envisioned as a potential constituent layer in a moiré-templated, strain-engineered vdW heterostructure of TMDCs [17].The authors demonstrated the systematic evolution of the interlayer coupling strength with twist angle in bilayer MoS 2 using a combination of Raman spectroscopy and classical simulations [17].Maity et al studied the effects of twisting on low-frequency shear modes and layer breathing modes in twisted bilayer TMDCs from the power spectra of mode-projected velocity autocorrelation function based on classical molecular dynamics simulations, which enabled the phonon calculations of small twist angles with large moiré superstructures [18].Finally, Quan et al developed a low-energy continuum model for phonons that overcomes the outstanding challenge of calculating the properties of large moiré supercells and successfully captures the essential experimental observations [19].That study, coupled with Raman measurements of tBLMs as a function of twist angles revealed phonon renormalization in this moiré superlattice.
Phonons are critical in understanding the thermal and mechanical properties of materials.In twisted bilayer MoS 2 , the interaction between layers significantly alters the phonon dispersion relations.The twist introduces a moiré pattern, a periodic spatial configuration, which influences the way phonons propagate through the material.This pattern can lead to changes in phonon band structures, potentially opening or closing phonon bandgaps, altering the material's heat conduction properties.The degree of twist between the two layers can modify the phonon spectrum in several ways.For small twist angles, the interlayer interaction is strong, leading to hybridized phonon modes.Conversely, at larger twist angles (closer to the 2D stacking configurations), the interaction weakens and the phonon properties start resembling those of individual monolayers.TBLMs exhibit a varied thermal conductivity as a result of the modified phonon dispersion.This is due to the fact that the moiré patterns can act as scattering centers for phonons, affecting the material's ability to conduct heat.This scattering can be tuned by changing the twist angle, offering a way to engineer the thermal properties of the material for specific applications, such as in nanoelectronics where heat management is of paramount importance [20].
In this article, we elucidate the Raman spectra of lowfrequency Raman active modes of twisted bilayer MoS 2 for several twist angles using an all-atom force-field and bond polarizability model (BPM) for MoS 2 .We show that twist angles near high-symmetry angles exhibit stacking frustration that leads to significant buckling of the moiré superlattice.We also demonstrate that this buckling realizes a rich soliton network and, in turn, permits a greater breadth of low-frequency Raman active modes.

Interatomic potential and phonons
We use classical molecular dynamics, assuming a nearequilibrium approximation for interatomic potentials to relax the atomic positions of the unit cells and to complete all the phonon calculations using a finite-difference scheme.Reactive empirical bond order (REBO) potentials [21] have proven to be an efficient and robust scheme for large atomistic simulations where chemical reactivity dominates interatomic interactions.Specifically, REBO potentials asssume a general potential of the form for atomic pairs {i, j}.Here, f(r ij ) is the cutoff function and V A (r ij ), V R (r ij ) are the attractive and repulsive parts of the potential,respectively.We consider a REBO potenital designed for MoS 2 [22] to describe Mo-S, S-S, and Mo-Mo interactions specifically.While inter-layer interactions are explicitly paramaterized in some schemes, such as the Kolmogorov-Crespi potential [23], we find that the MoS 2 layer thickness precludes an explicit need for normal vector considerations, since the innermost sulfur atoms dominate the interlayer interactions, rearranging in the local tangent plane to each layer to stabilize the stacking frustrations.
The coarse energy landscape of twisted bilayer MoS 2 complicates an accurate relaxation of the structure by traditional means, such as conjugate gradient descent.To address this complexity and guarantee the mechanical stability of the relaxed structure, we utilize the symmetries of the Γpoint phonons by imposing 3 translational (ω = 0) modes during relaxation.To calculate the modes and frequencies of the phonons, we construct the dynamical matrix, given for a general point q in the First Brillouin Zone by the mass-reduced Fourier transform of the force constant matrix and diagonalize at q = 0.If complex frequencies are found as a solution, the structure is not in mechanical equilibrium.Thus, the conjugate gradient descent is iterated until there are 3 translational phonon frequencies such that |ω| < 10 −5 cm −1 .In practice, the computational feasibility of this process is bounded by the stacking frustrations near high-symmetry angles, as well as the number of atoms in the unit cell of the moiré superlattice.We construct the Raman susceptibility tensors α ij using a BPM for MoS 2 , using polarizabilities obtained from density functional theory (DFT) data [24].In the BPM scheme, the derivatives of the polarizabilities with respect to the bond lengths ∂α ij /∂B βk are given in terms of bond-specific contributions whose parameters only depend on distance.In particular, where B is the unit vector between a pair of connected atoms in bond β, and α l and α p are the parallel and perpendicular components of the polarizability.From this, we see that the BPM model can be fully defined by the three coefficients ).We find that only considering Mo-S and Mo-Mo interactions for the construction of the polarizability tensor for bilayer MoS 2 is sufficient to match experimental data [11] and previous theoretical results [25,26].After energy minimization of the REBO potential, the structures exhibit an average length of 3.18 Å for Mo-Mo bonds and 2.45 Å for Mo-S bonds.
With this proviso, we can now obtain the Raman tensors , where n k is the Bose occupation factor and ω k is the frequency of the kth phonon mode.We note that the theory used here is a first-order theory, and only those modes at the Γ point of the Brillouin zone yields a Raman signal in the scheme employed.From the constructed Raman tensors, we are in a position to find the Raman intensities for each mode of tBLM for a pair of incoming and scattered polarization vectors e i , and e s respectively in the conventional method: [27] To more easily compare to experiment, we use a Gaussian convolution scheme across all of the intensities, assuming a standard deviation of the peak of 2 cm −1 to construct the continuous Raman spectrum, using Here, ω R , and R are a corresponding phonon frequency (measured in cm −1 ) and Raman tensor pair from the set of all phonon modes.Finally, we note that when the laser excitation energy corresponds to the energies of excitons such as the A, B, or C excitons in MoS 2 , Raman scattering process becomes resonant and Raman intensities are substantially enhanced in general.In other words, excitons enhance the electron-photon and electron-phonon couplings, which are both crucial for determining Raman intensities.[28,29] The BPM can capture the effect of excitons if the bond polarizability parameters are extracted based on the resonant Raman intensities of monolayer MoS 2 when the laser excitation energy corresponds to A, B, or C exciton [24].The magnitude of the polarizability parameters change accordingly based on resonant Raman intensities.If we use polarizability parameters corresponding to different laser excitation energies to generate Raman spectra of twisted bilayer MoS 2 , we should be able to capture the excitonic effects and show how Raman intensities evolve with the excitons.For interlayer excitons [30][31][32], we may need to fit the polarizability parameters based on resonant Raman data of bilayer MoS 2 when the laser excitation energy corresponds to the interlayer exciton.However, these considerations are beyond the scope of this work, and they will be studied in the future.

Structure definition
Our work here uses periodic boundary conditions and the construction of unit cells for moiré superlattices is based on the identification of commensurations between the top and bottom layers.In particular, such a commensuration can be defined by finding matrices obeying where T and B are 2 × 2 matrices whose rows are the basis vectors of the top and bottom layers respectively, and ∀z ∈ Z T , Z B , z ∈ Z. From this premise, one can show that any rational valued twist angle can be produced from coprime, non-negative integers p, q. [33] In particular, the angle of the commensuration is given by cos θ = 3p 2 − q 2 3p 2 + q 2 (6) and the number of lattice sites for the commensuration is given by Although there is a dense cover of possible twist angles in theory, many of these have more lattice sites than are computationally feasible.In this work, we examine 14 twist angles θ ∈ [0, 60 • ], which have fewer than 1305 atoms.We follow the procedure outlined in reference citemichael with the difference that atomic sites now alternate between a S and a Mo atom, and such that a twist angle of 0 • corresponds to the ideal AB stacking of the MoS 2 bilayer [26].

Results
For every twist angle considered here, we compute the Raman tensors of all phonon modes and classify a mode as active if the mode's peak intensity over all polarization angles considered is at least 5% of the maximum Raman intensity over all phonon modes and polarizations.If an E mode is included, its degenerate pair is also included as active even if the peak intensity across the polarization does not cross this threshold.The independence of the E mode activity is considered and discussed later.We consider polarization angles in the XY, and XZ planes, where e i = ⟨1, 0, 0⟩ and e s = ⟨cos (θ XY ) cos (θ XZ ) , sin (θ XY ) , sin (θ XZ )⟩.
We consider each of the polarization axes separately, so that the data for the XY polarization correspond to θ XZ = 0, and the data for the XZ polarization correspond to θ XY = 0.
Figure 1 shows the result of this procedure for the nontwisted bilayer MoS 2 in the AB configuration.The bandstructure obtained using the classical force-field potential is in good agreement with that obtained from DFT [25] and the Raman polarization plots indicate how the intensities change with the choice of polarization angle.Figures 2 and 3 show the result of the Raman calculation for four different angles.Those figures also show how the structural relaxation leads to signification displacement of the atoms, especially for small angles around the high-symmetry angles of 0 • and 60 • .In these cases, the periodic displacement of the Mo atoms clearly shows the realization of the soliton network which generates the lowfrequency modes not present in intermediate twist angles or at high-symmetry configurations.We see that twist angles near but not at high-symmetry angles exhibit unique Raman active modes.Data for intermediate angles listed in table 3, are similar to that found for θ = 49.0 • (not shown).

Space and point group analysis
Using the relaxed coordinates for each case, we can find the space and point group for each angle from which we can determine the general forms of the Raman tensors which are permitted by symmetry, as shown in table 2. We identify the Raman tensors allowed for point group 3 (C 3 ) as and for 3 m (C 3v ), −3 m (D 3d ), and 32 (D 3 ) as and For −3 m (D 3d ) the tensors are mapped as In our numerical approach, the E modes are found as a degenerate pair which are a linear combinations of the Frequencies at each q point are calculated by direct diagonalization of the dynamical matrix through the FBZ (right).The two polarization axes examined in this study are shown with respect to untwisted bilayer MoS 2 .We take impinging polarization e i = ⟨1, 0, 0⟩ and polarization of scattered light is varied from 0 • to 90 • in both schemes.At Γ, AB stacked MoS 2 include three high-frequency modes (one A 1g and two Eg modes) and three low-frequency modes (one single degenerate breathing mode and two degenerate shear modes), in addition to other modes that are not Raman active.Raman spectra and displacement maps of selected normal modes of representative low twist angles of tBLM, at 6.0 • and 3.9 • angles.(a) Shows the crystal structure of each of the bilayer systems viewed from the c axis and the displacement of Mo atoms after relaxation (symmetry guarantees that these are indeed identical for each of the layers).We find that the periodic displacement of the Mo atoms shows the realization of the soliton network which generates the low-frequency modes not present in intermediate twist angles or twist angles at high-symmetry angles.We conjecture that only near a high-symmetry angle can one achieve the structural periodicity of displacement that would influence the low-frequency Raman spectrum.(b) Raman spectra up to 120 cm −1 for XY (left) and XZ (right) polarization setups.(c) Shows the corresponding polar Raman plots, from which the A 1g symmetries become more apparent, vanishing at 90 • .There is an interestingly present E mode near 89 cm −1 that is only visible in the XZ polarization, due to the element e, which is the only non-zero component of the corresponding Raman tensor.
must generate the same Raman intensity as E α and E β do, denoted Table 3.All Raman active modes of twisted bilayer MoS 2 for the 14 twist angles considered in this study.Only three modes are seen in the high-symmetry structures (0 • and 60 • ); the breathing mode at ≈ 30 cm −1 , and the well-studied A and E modes in the higher frequency region of the spectra.Modes shown in blue indicate modes that are not present in the high-symmetry structures, indicating that the mode splitting at Γ is driven by the formation of the soliton network only seen in twist angles near high-symmetry angles [34,37].Note that many of the E modes only have one non-zero component, c.This implies that only one of these modes is active when in the XZ polarization, whereas the mode has a net isotropic intensity with respect to the polarization angle in the XY polarization.

Twist Angle
Stacking patches present in low (3.9 • ) and high (52.7 • ) twist angle tBLMs.We find that low-twist angles exhibit AA and AB stacking patches, while high-twist angles exhibit AA', A'B, and AB' stacking patches.These stacking patches correspond to buckling in the relaxed moiré superlattice as each of the stackings has slightly different equilibrium interlayer separations.
Noting that the E modes for 3 m (C 3v ), −3 m D 3d , and 32 (D 3 ) are a subset of the modes for 3 (C 3 ) with c = f = 0, we will now detail all Raman active modes and corresponding Raman tensors for the polarizations across XY and XZ.The result of this process, shown in table 3, is the main result of this work.

Discussion
The richness of the low-frequency modes near 0 • and 60 • can be understood by considering a superposition of the pure stacking modes of MoS 2 .Near high-symmetry twist angles, multiple stackings are present locally in large domains of the moiré superstructure's unit cell, thus driving a geometric frustration that leads to large periodic atomic displacements of comparable size to the unit cell.Near 0 • , both the AA and AB stacking regions are present and near 60 • , AA', AB', and A'B stackings are present as illustrated in figure 4, for twist angles 3.9 • , 6 • .The different equilibrium interlayer separations between the different stacking patterns cause this buckling.These modes which are only active in the presence of large stacking frustrations will be referred to as soliton-driven modes.To examine the behavior of the soliton-driven modes, we visualize the A modes for some of the near-high-symmetry angles, which can be seen to the right of figure 5. Indeed, these modes appear to rotate around sites of constant buckling, indicating that the soliton networks stabilize the quasi-circular trajectory of the modes.This indicates that soliton-driven E modes are only possible very close to the high-symmetry angles, such as seen in the 3.9 • , 6.0 • , and 52.7 • systems, whereas soliton-driven A 1g modes can be seen slightly further, such as in the 13.2 • and 49.0 • systems, though a more comprehensive analysis is required to determine exactly the threshold for this to occur and if it is a smooth or abrupt transition.We further note that the inversion symmetry, as determined by the point groups, of the 60 • and 0 • structures means that the Raman and infrared (IR) active modes are mutually exclusive and therefore Raman modes are not IR active at 60 • and 0 • .However, the point groups of all intermediate twist angle studied lack inversion symmetry, indicating that Raman modes at intermediate twist angles could be IR active as well.This is to our knowledge the only study of tBLM active modes that considers off-axis polarization.We find that this consideration is useful for highlighting certain low-frequency E modes, where the XZ elements of the Raman tensors may be appreciably larger than other off-diagonal terms.Visualized eigenstates of soliton driven phonons in tBLMs.A view normal to the x-axis and the z-axis is is shown for each mode, where atomic displacements were normalized to a maximum of 2Å per mode and shown in red along unit cells of molybdenum (white cycles) and sulfur (black circles).All twist angles that exhibit soliton-driven phonons have 2 A 1g modes, while only of the twist angles 3.9 • , 6 • , and 52.7 • exhibit soliton-driven E modes.The classification of E(x) and E(y) is based on the typical displacements generated by the normal modes.The visualizations were created with pyputil [38].

Figure 1 .
Figure 1.(left) The phonon band structure of AB stacked bilayer MoS 2 for a Γ-M-K-Γ path through the first Brillouin Zone (FBZ).Frequencies at each q point are calculated by direct diagonalization of the dynamical matrix through the FBZ (right).The two polarization axes examined in this study are shown with respect to untwisted bilayer MoS 2 .We take impinging polarization e i = ⟨1, 0, 0⟩ and polarization of scattered light is varied from 0 • to 90 • in both schemes.At Γ, AB stacked MoS 2 include three high-frequency modes (one A 1g and two Eg modes) and three low-frequency modes (one single degenerate breathing mode and two degenerate shear modes), in addition to other modes that are not Raman active.

Figure 2 .
Figure 2.Raman spectra and displacement maps of selected normal modes of representative low twist angles of tBLM, at 6.0 • and 3.9 • angles.(a) Shows the crystal structure of each of the bilayer systems viewed from the c axis and the displacement of Mo atoms after relaxation (symmetry guarantees that these are indeed identical for each of the layers).We find that the periodic displacement of the Mo atoms shows the realization of the soliton network which generates the low-frequency modes not present in intermediate twist angles or twist angles at high-symmetry angles.We conjecture that only near a high-symmetry angle can one achieve the structural periodicity of displacement that would influence the low-frequency Raman spectrum.(b) Raman spectra up to 120 cm −1 for XY (left) and XZ (right) polarization setups.(c) Shows the corresponding polar Raman plots, from which the A 1g symmetries become more apparent, vanishing at 90 • .There is an interestingly present E mode near 89 cm −1 that is only visible in the XZ polarization, due to the element e, which is the only non-zero component of the corresponding Raman tensor.

Figure 3 .
Figure 3. Raman spectra and displacement maps of representative high twist angles of tBLM, at 52.7 • and 49.0 • angles.The panel description is similar to that of figure 2. In contrast to the low-twist angle structures, these displacements indicate smaller periodic buckling.Interestingly, there is an E mode present near 107 cm −1 that is only visible in the XZ polarization due to the comparatively large e element of the Raman tensor.

Figure 5 .
Figure 5. Visualized eigenstates of soliton driven phonons in tBLMs.A view normal to the x-axis and the z-axis is is shown for each mode, where atomic displacements were normalized to a maximum of 2Å per mode and shown in red along unit cells of molybdenum (white cycles) and sulfur (black circles).All twist angles that exhibit soliton-driven phonons have 2 A 1g modes, while only of the twist angles 3.9 • , 6 • , and 52.7 • exhibit soliton-driven E modes.The classification of E(x) and E(y) is based on the typical displacements generated by the normal modes.The visualizations were created with pyputil[38].

Table 2 .
[35]e and point groups for all studied twist angles.Symmetry assignment, after relaxation, was performed using Phonopy with a 10 −5 threshold[35]CIF files for each structure are available at reference citedataset.