Phonons in Twisted Transition Metal Dichalcogenide Bilayers (“Twistnonics”): Ultra-soft Phasons, and a transition from Superlubric to Pinned Phase

The tunability of the interlayer coupling by twisting one layer with respect to the another layer of two dimensional materials provides a unique way to manipulate the phonons and related properties. We refer to this engineering of phononic properties as “Twistnonics”. We study the effects of twisting on low frequency shear (SM) and layer breathing (LBM) modes in transition metal dichalcogenide (TMD) bilayer using atomistic classical simulations. We show that these low frequency modes are extremely sensitive to twist and can be used to infer the twist angle. We find unique “ultra-soft” phason modes (frequency . 1 cm−1, comparable to acoustic modes) for any non-zero twist, corresponding to an effective translation of the moiré lattice by relative displacement of the constituent layers in a non-trivial way. Unlike the acoustic modes, the velocity of the phason modes are quite sensitive to twist angle. As twist angle decreases, (θ . 3◦, & 57◦) the ultra-soft modes represent the acoustic modes of the “emergent” soft moiré scale lattice. Also, new high-frequency SMs appear, identical to those in stable bilayer TMD (θ = 0◦/60◦), due to the overwhelming growth of stable stacking regions in relaxed twisted structures. Furthermore, we find remarkably different structural relaxation as θ → 0◦, → 60◦ due to sub-lattice symmetry breaking. Our study reveals the possibility of an intriguing θ dependent superlubric to pinning behavior and of the existence of ultra-soft modes in all two-dimensional (2D) materials.


I. INTRODUCTION
Twisting one layer of a bilayer system with respect to another provides a unique degree of freedom for tuning the properties of two dimensional (2D) materials. For example, in case of bilayer graphene, twisting leads to (a) structural changes, such as the observation of topological point defects, domain walls and layer buckling 1-8 , (b) significant change in electronic properties including superconductivity at "magic" twist angles [9][10][11][12][13][14][15][16][17][18][19][20][21] , and (c) superlubricity, a state of ultra-low friction [22][23][24] . An important facet of twisting is the evolution of low frequency vibrational modes, which has largely remained unexplored. Since the low frequency modes are solely determined by interlayer coupling and are accessible in Raman measurements, they provide a direct, non-destructive probe of the interlayer interaction [25][26][27][28] . The existing theoretical reports on the evolution of vibrational modes in twisted structures are restricted to large twist angle and use the Lennard-Jones potential 29,30 to describe the interlayer interaction, which is insufficient for capturing the stacking dependent energetics 31 . Although existing experimental studies 27,28,32 have explored small twist angles, they can only probe Raman active modes with frequencies > 10 cm −1 .
In this work, we computationally investigate the effects of twisting on low frequency shear (SM) and layer breathing (LBM) modes in bilayer MoS 2 , a prototypical transition metal dichalcogenide (TMD). Relative inplane and out-of-plane displacements of the constituent layers give rise to SM and LBM, respectively. The coexistence of several stackings in the moiré superlattice (MSL) that results from the twist leads to inhomogeneous interlayer coupling. As a consequence, the low fre-quency modes mix and become quite sensitive to twist. Our calculations show the existence of ultra-soft phason modes, large variation in LBM frequencies, appearance of multiple LBMs and high frequency SM in twisted bilayer MoS 2 (tBLMoS 2 ). Moreover, we find the velocity of the phason modes are quite sensitive to twist angle. These observations are generic to TMDs and we confirm our results for MoSe 2 as well. The domain walls and point defects present in twisted structures that are inevitable consequences of structural relaxation, are likely to influence the electronic properties [33][34][35][36] .

II. SIMULATION DETAILS
We use the Twister code 33 to create the MSL of bilayer MoS 2 with several different commensurate twist angles 1 • < θ < 59 • . The rigidly twisted structures are relaxed using LAMMPS [37][38][39] with the Stillinger-Weber and Kolmogorov-Crespi potentials to capture the intralayer 40,41 and interlayer interactions of tBLMoS 2 31,42 , respectively. We use PHONOPY 43 to compute the zero temperature vibrational spectra of the relaxed tBLMoS 2 . The interlayer KC potential consists of two types of interaction : (i) nearest neighbor interlayer X − X coupling, and (ii) nearest neighbor interlayer M − X coupling 31 , where M denotes metal and X denotes chalcogen atoms.
Independently, we also compute the low frequency modes from the power spectra of mode projected velocity auto-correlation function (mVACF) 44,45 from classical molecular dynamics simulations with periodic boundary conditions in the canonical ensemble using the Nosé-Hoover thermostat in LAMMPS. The time averaged where j denotes the atom type in the unit cell, andê q,s denotes the eigenvector. The mass weighted momentum projected velocities are defined as, where r jk are the atomic coordinates, k denote atoms belonging to particle type j and m j denotes atomic mass. We equilibrate the system in a NVT ensemble for ∼ 150 − 300 ps at T = 300 K. The supercell contains 32 MSL in order to correctly capture all the anharmonic effects. The power spectra of the mVACF projects the full phonon spectra to a particular branch for any q. Typically, we collect velocities of all atoms in the production run (480 ps, in N V E ensemble) every 20 timesteps (1 timestep = 1 fs). For computational efficieny, the 480 ps tracjectory is divided into 6 parts, each containing 80 ps. To compute the mVACF, we use the definition of unbiased estimator 45 . Finally, we use FFT to compute the power spectra. Instead of computing the mVACF for each eigenmodes in MSL, we use BLMoS 2 SM and LBM eigenvectors to compute the mVACF of the superlattice.
For θ → 0 • , there are two unique high-symmetry stacking regions 33 , AA (Mo, S of top layer are directly above Mo, S of bottom layer, respectively) and B Mo/S (Bernal (B) stacking with Mo of top layer directly above S of bottom layer, equivalent to B S/Mo ). For θ → 60 • , there are three unique high-symmetry stacking regions 33 , AB (Mo, S of top layer are directly above S, Mo of bottom layer, respectively), B Mo/Mo and B S/S . Among these high-symmetry stackings B Mo/S (AB) is the most stable with SM frequency ∼ 21 cm −1 , whereas AA (B S/S ) is unstable with strong imaginary SM frequency. Due to the difference in binding energies of different stackings (consequently, stability and interlayer separation (ILS)), upon relaxing the MSL the more stable stacking regions increase in area. The signatures of the growth of the stable stacking regions with θ are inherently embedded in the ILS landscape.

B. Relaxation : Interlayer Separation Landscape
For the calculation of vibrational properties i.e. perturbation with respect to the ground state, the relaxation of the twisted structures are necessary to obtain a suitable ground state. The relaxation of the rigidly twisted structures involves straining and buckling of each constituent layers. In order to demonstrate the consequences of relaxation, in Fig.2 we show the θ dependence of the ILS landscape and ILS av . We can identify the high-symmetry regions in the ILS landscape for θ → 0 • (Fig. 2a) : alternate triangles with the least ILS (B Mo/S and B S/Mo , deep blue), and red circles with the maximum ILS (AA). Both B Mo/S and B S/Mo regions grow equally as θ → 0 since they are degenerate. Six domain walls (light blue lines) meet at the "centers", where AA stacking (topological point defect) regions are located, similar to what happens in twisted bilayer graphene (tBLG) 1,2,5,6 . The ILS landscape for θ → 60 • , on the other hand, (Fig.2b) is very different due to sub-lattice symmetry breaking. AB regions (Reuleaux triangle like) grow overwhelmingly due to their relatively higher binding energy. Six curved domain walls meet at the "centers", where B S/S (red circles, maximum ILS) stacking regions are located. Generalized stacking fault energy in conjunction with continuum theory is also shown to predict similar in-plane features 46 for tBLMoS 2 . It is interesting to note that both the length and the shape of domain walls can be tuned with twists as θ → 60 • . Figure 2c, 2d capture θ dependence of ILS av , ILS min and ILS max of the MSL. For large twists (13 • < θ < 47 • ), the absence of any extended ideal high-symmetry stacking regions leads to θ-independent behavior of ILS av , ILS min , and ILS max . ILS min and ILS max saturate for θ < 3 • and θ > 57 • . Remarkably, ILS av doesn't saturate in this limit due to the presence of AA/B S/S and domain walls.

C. Inhomogeneous Interlayer Coupling
As a whole, twisting affects the phonon band structure in two ways. First, the shrinking of the Brillouin zone gives rise to folded phonon modes. Second, as there are multiple stackings in the MSL, the interlayer coupling is inhomogeneous, which leads to mode mixing among in-plane and out-of-plane modes. Relaxation of the rigidly twisted structures further changes this mode mixing and stabilizes the structure (SM frequency of unrelaxed tBLMoS 2 is strongly imaginary, implying instability). In order to separate the effects due to only zone folding and the mixing of modes due to presence of multiple stacking, we compute A ij = | ψ i BL |ψ j N I | 2 , projections of bilayer eigenmodes (|ψ i BL ) onto individual layer modes (|ψ j N I ) at Γ for both untwisted and twisted structures (Fig.3a, 3b). The untwisted and twisted structures are of similar dimensions. The projections would have been identical if the zone folding effects were the only factor determining phonons in twisted structures. It is clear from the figure that inhomogeneity in the interlayer coupling in the twisted structure, leads to greater mode mixing. The comparison of mode mixing between untwisted and twisted structures are essential as the mixing strongly depends on twist angle and produces nontrivial effects as we illustrate below. This observation is also crucial, as it invalidates the usage of a simple linear chain model to compute SM and LBM frequencies.

D. Shear and Layer Breathing Modes
Next, we focus on the effects of twisting on the SM and LBM frequencies with (at T = 300 K) and without thermal fluctuations (T = 0 K). While computing the mVACF (for T = 300 K) we use the SM and LBM eigenvectors of BLMoS 2 . As discussed above, due to mode mixing the eigenvectors of the twisted structures can be composed of several normal modes of the MSL. Any nondegenerate eigenmode involving relative displacements of the layers of the MSL should appear as a distinct peak in the power spectra of the mVACF. In order to compare with finite T results we also project the eigenmodes of MSL onto the BLMoS 2 SM and LBM eigenvectors, p = | ê MSL |ê BL | 2 at T = 0 K (Fig.3c). It should also be noted that, multiple zone "folded" modes from SM and LBM branches of different q points of the BLMoS 2 unit-cell Brillouin zone will appear at the Γ point of the MSL Brillouin zone. However, the projections of these "folded" modes onto BLMoS 2 SM and LBM eigenvectors at Γ point will be significantly smaller than that of the modes not arising due to zone folding. This is due to orthogonality of the eigenvectors at different q points. Therefore, considering the evolution of the eigenmodes with the largest projections (in fact, we consider all modes with p SM , p LBM > 0.1 as SM and LBM), we infer the twist angle dependence of shear and layer breath-ing modes. We can categorize the θ dependence of the low frequency SM and LBM into three regions. Since the change of the low frequency modes with twist angle is a result of complicated evolution of the interlayer coupling and mode mixing, we present both the frequencies and the eigenvectors of the SM and LBMs. The twist angle dependence of the eigenvectors provide deeper insight to the evolution of the low frequency modes, as we illustrate below.
Region(I) (7 • θ 53 • ) : In this region, we find an averaged LBM and exceedingly small SM frequencies (0-2 cm −1 , "ultra-soft"). The LBM frequency decreases monotonically as θ → 7 • /53 • from larger twists (Fig.4). Furthermore, the projections p LBM > 0.9, p SM > 0.9 indicate that, the nature of vibrations of SM and LBMs remain similar to that of BLMoS 2 . As a representative of this region, we show the eigenvectors corresponding to SM and LBM in Fig.6g,6h,7g,7h,5c for θ = 23.48 • . Clearly, the SM (LBM) eigenvectors correspond to the relative horizontal (vertical) uniform displacement of the constituent layers. This is due to the absence of any extended high-symmetry stacking in the MSL, which leads to nearly uniform interlayer coupling. However, as we approach Region(II), the LBM starts mixing with inplane modes giving rise to the monotonic decrease in the LBM frequencies. The change in LBM frequencies (by ∼ 1.5 cm −1 ) can be used to reliably infer large θ. Moreover, the presence of LBM also indicates that the layers in the twisted structures are not completely decoupled (in the out-of-plane direction). in-plane modes. As an example, we show the LBM eigenvector with maximum p LBM for θ = 5 • (Fig.5b). The inplane modes corresponding to the LBM is of similar order in magnitude. In this region all the high-symmetry stackings, and domain walls occupy comparable area-fraction of the MSL. This enhances the mode mixing. The ultrasoft SMs are also present in this region. However, the corresponding SM eigenvectors are no longer completely uniform (Fig.6e,6f,7e,7f) and start mixing with out-ofplane modes. In short, this region represents the transition from completely mismatched stacking (Reg.I) to highly ordered stable stacking regions separated by domain walls (Reg.III). ther (θ → 0 • , or θ → 60 • ) we find one LBM (with significantly large p LBM ) with frequency similar to that of stable BLMoS 2 . This is due to the overwhelming growth of the stable stacking regions (Fig.2) in this region. The mixing of LBM with in-plane modes although exists, are smaller compared to Reg.II (Fig.5). In order to establish this, we plot the LBM eigenvector with largest p LBM for θ = 1.9 • (Fig.5a). It is evident from the figure that, the LBM primarily arises from the out-of-plane vibration of the B Mo/S stacked region. On the other hand, the SM frequencies split essentially in two branches : one ultra-soft in nature (similar to Reg.I, II) and one with high frequency (similar to BLMoS 2 , 22-28 cm −1 ). Moreover, as θ decreases from 3 • to 0 • (57 • to 60 • ) the highfrequency SM redshifts (Fig.3d,3e). The splitting of the SM has important consequences in frictional properties. The ultra-soft SMs are localized on domain walls and AA stacking (domain walls and B S/S near 60 • , Fig.6a, 6b,  7a, 7b). On the contrary, high-frequency SMs primarily originate from the relative displacement of stable stacking regions (Fig.6c, 6d, 7c, 7d). Interestingly, the saddle point nature of the domain walls is evident from the inplane displacement of the eigenvectors. For the ultra-soft modes the eigenvectors are parallel to the domain walls connecting AA stacking and for high frequency modes perpendicular to domain walls connecting B Mo/S , B S/Mo (Fig.6,7). This splitting of SM frequencies occurs due to significant growth of the stable stacking region in the

E. Ultra-soft Shear Modes : Phasons
The presence of the ultra-soft modes with twist angle is one of the major finding of our work. Thus, we investigate three important aspects of this finding : the origin of these modes, their twist angle dependence and consequences on frictional properties.
The ultra-soft modes represent an effective translation of the MSL by local relative displacements of the atoms in the constituent layers (see attached movies). Since the frequencies associated with these modes are also very small (almost acoustic mode like), this implies that under the relative displacements of two layers (following the eigenvectors shown in Fig.6,7) the energy of the tBLMoS 2 remains invariant (or almost invariant), which can be a consequence of continuous symmetry breaking. For example, the in-plane acoustic modes (LA, TA modes) which represent invariance of the total energy under global translation originate due to translational symmetry breaking. However, strictly speaking, we show that the ultra-soft SMs found in our calculation are optical modes and does not represent continuous symmetry breaking. The optical nature is clearly reflected in the dispersion relation (dω/dq ≈ 0, for small q very near to Γ point) unlike the acoustic modes (Fig.8, small negative values at Γ , −0.2 cm −1 , are within numerical accuracies of our calculation). We highlight both the acoustic and ultra-soft modes in order to show this difference clearly. However, the ultra-soft nature of these SMs can be understood from one-dimensional Frenkel-Kontorova (FK) model. In this model, a linear chain of atoms is subjected to an external periodic potential. Depending on the ratio of the periodicity of the potential and linear chain, the structures can be commensurate or incommensurate. For simplicity, we assume the incommensurate (commensurate) structure as a lattice with infinite (finite, small) periodic length. Without considering dissipative coupling, the incommensurate structure possesses a gapless Goldstone mode (phason) with linear dispersion due to invariance of the phases of two mass density waves under uniform relative displacement. When a commensurate state is approached, the phason becomes gapped 47 with dω/dq ≈ 0. This is exactly what happens in the case of ultra-soft SMs found in our calculations. The linear chain of atoms and the external periodic potential in the FK model are replaced by the constituent MoS 2 layers of tBLMoS 2 , and stacking dependent binding energy, respectively. Since we only simulate commensurate angles, the phasons are always gapped (although ultra-soft) optical modes with dω/dq ≈ 0. Hence, in our calculation the phase invariance associated with these modes are always approximate. In the case of incommensurate twist angles, however, the phase invariance becomes exact and we expect corresponding gapless phason modes with ω ∝ q.
The exact eigenvectors for the ultra-soft SMs can not be directly probed by experiments. Keeping in mind the possible experimental signatures of twist angle dependence, we compute the dω/dq (group velocity) of these ultra-soft modes. We compare them with the velocities of the acoustic modes ( Fig.9) away from Γ point using the linear dispersion of these modes. The velocity of the acoustic modes remains invariant with respect to the change of twist angle, whereas the velocity of the ultrasoft phason modes change significantly (by a factor of 2-3 at θ = 1.9 • , 58.1 • ). Physically, the in-plane acoustic modes represent the long-wavelength vibrations of the entire twisted bilayer MoS 2 , with all the atoms of two layers participating. At Γ point, these modes represent uniform translation of all the atoms in the MSL. Therefore, the acoustic mode velocities within the harmonic approximation can be written in terms of the Lamé coefficients, µ, λ in the following manner : v LA = (λ + 2µ)/ρ and v TA = µ/ρ 47 . The presence of the MSL only affects the motion involving relative displacement of the constituent layers. Since, the LA and TA modes only correspond to "in-phase" motion of the layers, they remain unaffected irrespective of twist angle. Moreover, the rigidity (governed by λ, µ) of the lattice for the in-plane acoustic modes are the same as that of single layer MoS 2 . This is evident from the twist angle independent behavior of the velocities of the LA, TA modes (Fig.9). On the other hand, the ultra-soft phason modes are localized on domain walls and AA stacking for θ → 0 • (domain walls and B S/S stacking). They represent "acoustic" modes of the "emergent" soft moirè scale lattice (contains only domain walls and AA near 0 • , whereas domain walls and B S/S near 60 • ). Strictly speaking, they become acoustic modes only for incommensurate structures, as pointed out earlier. The soft nature of the lattice is noticeable from the change of velocity of these modes (Fig.9). Similar to the LA, TA modes the velocities of the phason modes are functions of λ E , µ E , where λ E , µ E represents effective Lamé coefficients of the moirè scale lattice. The decrease in ultra-soft modes velocity imply reduction in strength of the effective Lamé coefficients. This is because ultra-soft modes originate from the "out-of-phase" motion of the two layers. Thus, the strength of the interlayer coupling gives rise to effective Lamé coefficients. We leave details of the consequences of this emergent soft lattice and investigation of their rigidity for a subsequent work.
Before discussing the consequence of the twist angle dependence of the ultra-soft SM frequencies on frictional properties, we briefly summarize previous experimental studies on this aspect. The introduction of a twist between two layers of two dimensional materials has been shown to modify the frictional properties dramatically [22][23][24][48][49][50][51][52][53][54] . Depending on the twist angle, the frictional force can be exceedingly small (structural superlubricity) or moderately large. For instance, the frictional force can drop by 2-3 orders of magnitude if the twist angle is > 5 • or < 55 • in the case of graphene on graphite, twisted bilayer MoS 2 22,48,50,55 . Also, for 5 • θ 55 • , the frictional properties remain almost constant. In all these examples, maximum friction is obtained when the two layers are untwisted. Although microscopic details of the interlayer sliding process giving rise to superlubric behavior can be quite involved, the primary reason behind superlubricity is often attributed to incommensurability between two surfaces. Here, we identify that the change of frictional properties with θ are intimately related to the existence and evolution of the ultra-soft modes. For θ = 0 • /60 • , the BLMoS 2 is unit cell commensurate. When trying to shear one layer with respect to another in BLMoS 2 , all unit cells have to cross the interlayer-sliding barrier simultaneously. This leads to large SM frequencies (∼ 21 cm −1 ) and high-friction. This can also be understood from the aforementioned one dimensional FK model. In the unit-cell commensurate case ( (Fig.10a)), while shearing the atoms globally with respect to the external periodic potential, every atom rises toward the hill. Thus, the cost of shearing is large. However, if the system is large-scale periodic (Fig.10b), then every atom has to cross variable sliding barrier. This reduces the shearing energy, due to cancellation of cost while rising up the hill and gain while falling from it. Similarly, tBLMoS 2 for large twist angles (θ 5 • , 55 • ) is periodic at larger scale. While shearing unit cells have to cross a variable interlayer-sliding barrier (due to coexistence of multiple stackings with different binding energies and stability). In effect, this drastically reduces SM frequencies and hence friction as well (superlubric). This can be easily confirmed from the existence of the ultrasoft phason modes and the absence of high-frequency SM. However, as θ decreases further not only the highfrequency SM similar to the untwisted bilayer appear, but also the ultra-soft modes start to localize on the domain walls (splitting of SM). This immediately implies that, to shear one layer with respect to the other globally a large barrier has to be overcome. Because, only a small fraction of atoms participate in the ultra-soft modes in this case unlike large twist angle. A significant fraction of atoms participate in the high-frequency SM implying pinning. The origin of the pinning lies in the significant growth of stable stacking when θ → 0 • /60 • (Fig.2) leading to the development of large interlayer sliding barrier against shearing (Fig.10d,10e). Similarly, if the potential is strong enough in the FK model, most atoms sit at the minima. Only a small number of atoms occupy the energetically unfavorable hills, known as discommensuration (Fig.10c). This leads to pinned phasons in the FK incommensurate structures. The nature of vibrations of ultra-soft phason modes in our calculations (uniform at large θ due to complete stacking mismatch and nonunifrom at small θ due to overwhelming growth of stable stacking) indicate the possibility of having pinned phasons (Aubry-like transition) in the small twist incommensurate structures 47,[56][57][58] . Interestingly, similar pinning behavior has also been realized in systems like colloidal monolayers in optical lattices 59,60 , and physisorbed sub-monolayers on crystal surfaces 61 .

IV. DISCUSSIONS
Other 2D materials: Twisting one layer with respect to another in bilayers of 2D materials leads to the presence of multiple types of stacking in the MSL with different binding energies and stability, irrespective of the materials electronic properties. The intra-layer interaction in 2D materials are far stronger than the interlayer coupling. The combination of strong in-plane stiffness and weak variable interlayer coupling of twisted structures should produce similar behavior of low frequency vibrational modes in any 2D material. Thus, the existence of ultra-soft phason modes, the twist angle dependence of the corresponding eigenvectors and velocity are expected to be generic to any 2D materials. The twist angle dependence of the interlayer coupling is also revealed in the relaxed twisted structures. In order to illustrate, we compute the interlayer separation landscape and low . frequency vibrational modes for twisted bilayer MoSe 2 using mVACF at T = 300 K (Fig11). The trends are quite similar as in the case of tBLMoS 2 justifying our conclusion.
Raman Spectroscopy : We have demonstrated the evolution of shear and layer breathing mode frequencies in twisted TMD bilayers. Irrespective of their Raman sensitivity, these modes are present in the twisted bilayer structures. For BLMoS 2 (untwisted, most stable stacking), both the SM and LBM are found to be Raman active 26 . The observation of these modes using Raman spectroscopy depends on the Raman scattering intensity 26 . Since both the SM and LBM in the case of BLMoS 2 are Raman active, we expect in the twisted structures SM and LBM with largest p SM , p LBM will also be Raman active. However, this is speculative and explicit calculations of Raman intensity in twisted structures are unfeasible at present. Hence, we compare our results directly to Raman measurements in the case of tBLMoSe 2 to show the usefulness of our calculations. The trends of the twist angle dependence of low frequency vibrational modes in tBLMoSe 2 are in excellent agreement with experiment 28 . Similar to tBLMoS 2 (Fig.3), we categorize the θ dependence of SM and LBM frequencies. In Reg.(III), the high-frequency SM appears and LBM frequencies are similar to that of BLMoSe 2 , as predicted by our calculations. Multiple LBM and significant variation in LBM frequencies are also present in Reg.(II). Both in Reg.(I), (II) the SM frequencies are absent (because they are ultra-soft in nature) in Raman spectroscopy. However, two prominent features seem to be missing in the Raman experiment 28 : (i) Monotonic decrement of LBM frequencies in Reg.I (as in Fig.4) and (ii) slight stiffening of the SM frequencies compared to BLMoSe 2 (Reg.III). This might be due to the fact that, the experiment is carried out with CVD grown samples, which is known to exhibit lesser mobility and more disorder. We expect samples of better quality, such as mechanically exfoliated structures can show these missing features. Brillouin-Mandelstam Spectroscopy : Due to ultrasoft nature ( 1 cm −1 ) of the phason modes Raman spectroscopy can not be used to probe them. Brillouin-Mandelstam Spectroscopy (BMS) which can probe small frequencies (typically, 0.1 − 6 cm −1 ), is often used in mineral physics and material science to probe acoustic modes 62 . The phonon dispersion can also be mapped out by changing the incident light angle in BMS 63 . Therefore, the sound speed along with elastic rigidity can be experimentally obtained. Using BMS spectroscopy our predictions of the existence of the ultra-soft modes can be verified. Also, by measuring the twist angle dependence of the dispersion of the acoustic and ultra-soft modes our results on velocity dependence of ultra-soft modes can be tested. Although, it is more likely to find incommensurate twisted structures in experiments and thus, the phason modes may become completely gapless. Still, we expect the phason velocity should be twist angle dependent as predicted here.
More Twistnonics : We have demonstrated the manipulation of low frequency vibrational modes with twist and related effects in frictional properties in transition metal dichalcogenide bilayers. Since phonon plays important role in determining many other material properties, the ability to control the phonon dispersion will facilitate engineering of those properties 64 . Here, we outline a few of them : (i) The existence of ultra-soft modes can strongly modify the low temperature specific heat 65,66 , one of the key variable that dictates thermodynamic properties of a material. In order to illustrate this, we calculate the difference in specific heat of twisted structures with respect to BLMoS 2 . Clearly, there are significant changes in δC v for T < 200 K. Corresponding single layer MoS 2 value is also shown. (ii) Acoustic phonons are known to be the dominant heat carriers in insulating materials. The pres- ence of ultra-soft phason modes can significantly modify the thermal conductivity of the twisted structures. Furthermore, as the group velocities of the ultra-soft phason modes can be tuned with twist angles, the thermal conductivity can also be engineered. (iii) The electronphonon coupling for the ultra-soft modes can play important role in determining electrical resistivity. For example, in the case of twisted bilayer graphene electron-"acoustic" phonon coupling has been identified as the dominant source of it's high T resistivity behavior 67 . The authors assumed the relative displacement of the two layers in twisted bilayer graphene as the additional "acoustic" modes of the system with identical acoustic modes (LA, TA) dispersion. Although, this model captures the features of the resistivity qualitatively, the D/v ph (D is deformation potential for the electron-"acoustic" phonon coupling, v ph is the phonon velocity) appearing in theoretical calculations seems to be off by a factor of 2-3 compared to experiment 67,68 . The assumption of the shear modes as the "acoustic" modes is reasonable in the twisted structures, as they are ultra-soft in nature. However, the assumption that the dispersion of the ultra-soft modes and acoustic modes are identical is not accurate. Our calculations for the tBLMoS 2 clearly shows the twist angle dependence of the velocity of the phason modes compared to that of the acoustic modes (by a factor of 2-3). This might provide an explanation for the missing factor in the case of tBLG.
Since the SM and LBM are interlayer coupling dependent, by modifying the interlayer interaction strength (with external pressure, for example) the phonon spectra can be further engineered. Another important degree of freedom to tune phononic properties is the choice of materials, for instance, van der Waals heterostructures (stacking of two dissimilar 2D materials).
Structural Relaxation: The ILS landscape of the twisted TMD bilayer shows distinct relaxation pattern for θ → 0 • and θ → 60 • . Other 2D materials with sub-lattice symmetry breaking (e.g. hBN 69 ) should exhibit similar distinct structural relaxation as well, which can be probed using electron microscopy, infrared nanoimaging 2,3 etc.

V. CONCLUSION
Here, we have shown that the low frequency modes are extremely sensitive to twist in twisted TMD bilayer and can be used as a probe to determine twist angle. We also make falsifiable predictions of the presence of ultra-soft phason modes and the twist angle dependence of their eigenvectors, velocities. Our results indicate a twist angle dependent transition from superlubric to pinnied state.
Our study provides the first step to twistnonics in 2D materials.

VI. ACKNOWLEDGEMENT
The authors thank Rahul Debnath, Shinjan Mandal, Sumilan Banerjee, Subroto Mukerjee, Arindam Ghosh and Rahul Pandit for useful discussions and thank the Supercomputer Education and Research Center (SERC) at IISc for providing computational resources. S.R. thanks SERB India and the Tata Education and Development Trust for support.
Note Added: While submitting the revised version we become aware of two recent studies, carried out within an elastic continuum approximation, reporting the existence of similar ultra-soft phonon modes in tBLG 70,71 .