Ferroelectric switching in bilayer 3R MoS2 via interlayer shear mode driven by nonlinear phononics

We theoretically investigate the mechanism of ferroelectric switching via interlayer shear in 3R MoS2 using first principles and lattice dynamics calculations. First principle calculations show the prominent anharmonic coupling of the infrared inactive interlayer shear and the infrared active phonons. The nonlinear coupling terms generates an effective anharmonic force which drives the interlayer shear mode and lowers the ferroelectric switching barrier depending on the amplitude and polarization of infrared mode. Lattice dynamics simulations show that the interlayer shear mode can be coherently excited to the switching threshold by a train of infrared pulses polarized along the zigzag axis of MoS2. The results of this study indicate the possibility of ultrafast ferroelectricity in stacked two-dimensional materials from the control of stacking sequence.

, its realization in the experiment is challenging because MoS 2 is stable in semiconducting 2H phase rather than metallic 1T structure 14 .
The ferroelectricity is likely to appear in stacked TMDC rather than in single layer form. It was recently shown that the horizontal mirror symmetry of individual layers is broken by the stacking in 3R structure (P3m1 for finite layers and R3m for bulk), hence the vertical electric polarization manifests itself in accordance with the global polar symmetry 15 . The direction of polarization depends on the stacking sequence, hence is reversed by the interlayer translation between the AB and AC stackings 15,16 as shown in Fig. 1a. The possibility of ferroelectric switching in the 3R structure via the interlayer translation has not been explored to date. Meanwhile, the multilayer distorted 1T WTe 2 showed switching of the polarization 17,18 and the topological phase 12 in the recent experiments, probably via the interlayer translation. The stability of the 3R structure MoS 2 is comparable to that of the 2H structure and can be selectively synthesized among competing polytypes 19 . It is hence a viable candidate for 2D ferroelectrics in which the fascinating phenomena such as high electron mobility 16 and valleytronics 15,19 can be explored altogether.
The interlayer translation universally manifests as Raman active low-frequency lattice vibrations in layered 2D crystals due to the weak van der Waals bonds (vdW) between layers [20][21][22] . Therefore, the ferroelectric switching using an optical field based on ionic Raman scattering is considerably appealing [23][24][25] . A particular mechanism, called nonlinear phononics, relies on the anharmonic phonon coupling between infrared active and targeted secondary vibrational modes, which displaces the crystal toward the reversal of polarization upon the irradiation of intense terahertz pulse 26,27 . The ionic Raman scattering is distinguished from the conventional Raman scattering which has been used to detect structural characteristics of 2D materials such as stacking structures 28 and local bonding chemistry 29 . The use of pulse with the mid-infrared frequency within short duration allows the exploration of an extremely intense light field (peak electric field reaching ~600 MVcm −1 ) without the material damage [30][31][32] . The optical ferroelectric switching is a rapidly growing topic, which will enable the ultrafast and nondestructive way to achieve coherent switching 26,27,[33][34][35] .
In this work, we theoretically show the possibility of the ferroelectric switching of the bilayer 3R MoS 2 using the intense light pulse through the anharmonic phonon coupling. Density functional theory (DFT) calculations demonstrate that a large amplitude vibration of infrared mode can effectively lower the ferroelectric switching barrier, and induce an unidirectional anharmonic force on the interlayer shear mode along the switching direction. This effect depends on the polarization angle of the incident light pulse with respect to the crystallographic axis of MoS 2 according to the selection rule. Lattice dynamics simulations indicate the possibility of dynamical ferroelectric switching through the coherent amplification of the interlayer shear mode and the lowering of the energy barrier under the repetitive pulses within a few picoseconds.

Results
Electric polarization and switching in bilayer 3R MoS 2 . The bilayer 3R MoS 2 has the polar point group symmetry of C 3v with the polar axis along the z-axis. The 3R structure can be constructed by either AB or AC stacking sequence, which develops the spontaneous electric polarization in the opposite direction of each other as shown in Fig. 1a. The magnitude of electric polarization of the bilayer structure was calculated from the Berry phase method as P = 0.24 μCcm −2 , in agreement with previous reports 15,16 . Figure 1b shows the total energy of the bilayer structure (primitive cell consisting of 6 atoms) depending on the stacking sequence calculated from the nudged elastic band (NEB) method. Both the AB and AC stackings are stable and energetically degenerate structures. The AC stacking is obtained from the AB stacking by sliding the upper B layer along the +y direction by 1.82 Å. The interlayer translation over the weak vdW interaction results in the modest energy barrier Φ NEB = 15.0 meV.
The optical switching mechanism was investigated based on the nonlinear phononics 23,24 . The bilayer 3R MoS 2 has 18 zone-center phonon modes which are decomposed into Γ = 6A 1 + 6E representations. The singly degenerate A 1 mode involves out-of-plane motion of atoms, while the doubly degenerate E mode involves in-plane motion of atoms. Figure 1c shows two kinds of E modes relevant to the nonlinear phononics mechanism. The low-frequency mode (Ω LS = 0.6 THz in phonon dispersion in Fig. 1d) referred to as the interlayer shear mode (denoted by Q LS ) involves the relative motion between adjacent layers along the in-plane polarization axis r. Therefore, the Q LS mode is related to the AB ↔ AC stacking change. The infrared activity of the v-th mode is proportional to the square of mode effective charge Z v

⁎36
. Due to the almost rigid relative ionic motion, the Q LS mode does not produce net dipole moment as the calculated effective charge ⁎ Z LS = 0.00 eμ −1/2 (where e is the electronic charge and μ is the atomic mass unit). The vanishing infrared activity indicates that it is almost impossible to directly excite the Q LS mode with large amplitude to induce the stacking change. Nonetheless, the anharmonic coupling of Q LS mode with other infrared active modes can provide an alternative route to control this mode, and the consequent ferroelectric switching. Among the other in-plane modes, only the high-frequency mode (denoted by Q IR ) at Ω IR = 11.4 THz shows finite effective charge Z IR ⁎ = 0.23 eμ −1/2 , and is the solely infrared active mode under vertical incidence of light.
The normal-mode coordinate where m i is the atomic mass of i-th atom and e i v is the normalized eigenvector of the dynamical matrix. The orthogonal basis sets were chosen to represent the degenerate Q LS and Q IR modes as {Q LSx , Q LSy } and {Q IRx , Q IRy }, respectively. They correspond to the linear polarization along the zigzag (x-axis) and armchair (y-axis) axes shown in the top view of the AB stacking of the 3R MoS 2 in Fig. 2a. The AB stacking deformed by the positive amplitude of Q LSy = +16.29 Å μ (atomic displacement of the each of the adjacent layers in the opposite direction by 0.91 Å) corresponds to the AC stacking. Meanwhile, the deformation by the negative amplitude of the Q LSy = −16.29 Å μ changes the AB stacking into the unstable AA stacking. The positive amplitude of Q LS along three crystallographically equivalent directions, r 1 (+y direction), r 2 and r 3 directions (−120° and +120° from the +y direction), equally change the AB stacking into AC stacking as shown in Fig. 2a.
The possibility of the ferroelectric switching hinges on how much the Q LS mode can be amplified along the desired direction for the AB ↔ AC stacking change by the coupling with Q IR . Here, the anharmonic coupling property was investigated from the potential energy surfaces as a function of normal-mode coordinates. The potential energy surfaces V(Q IR , Q LSx , Q LSy ) for each Q IRx and Q IRy were calculated using DFT on 21 × 21 × 23 points with steps of 0.82 Å μ for Q IR and 1.63 Å μ for Q LS modes. The energy surface was then fitted to the polynomial function as where Q IR is either Q IRx or Q IRy , and c lmn is the anharmonic coefficient, and Q IR l , Q Q and LSx m LSy n denote the l, m and n powers of the normal-mode coordinates, respectively. Using this expression, we analyze the effect of the irradiating light pulse with the linear polarization along the x-or y-axis, thus exciting Q IRx or Q IRy mode, respectively. Note that the normal-modes in the cartesian basis are classified into the odd parity modes (Q LSx and Q IRx ) and even parity modes (Q LSy and Q IRy ) under the mirror symmetry σ 1 shown in Fig. 2a coupling such that c lmn is nonzero only for l + m = even (m = even). The terms were included up to 15 th power terms (l + m + n = 15) in the polynomial function, which fits the DFT potential energy surface accurately. The representative coupling terms are displayed in Table 1. Figure 2b shows the potential energy surface V(Q IR , Q LSx , Q LSy ) represented on (Q LSx , Q LSy ) coordinates when the amplitude of Q IR is zero. The energy contour shows the directional dependence inherited from the C 3v symmetry. The energy barriers for the AB → AC change along the equivalent r 1, r 2 and r 3 directions in this potential energy surface are the same as Φ 0 = 17.3 meV. The difference between Φ 0 and Φ NEB for the AB → AC stacking change is because the deformation by the in-plane Q LS mode does not include any out-of-plane relaxation, while the NEB path includes the relaxation from the slight increase (~1.6%) of the interlayer distance, reducing the barrier. It is worth to note that the Φ 0 rather than the Φ NEB is relevant to the ultrafast switching in the picosecond time scale, while the latter is relevant to the conventional switching in a longer time scale.
The anharmonic coupling effect can be seen from the modulation of the potential energy landscape under the large Q IR amplitude. Figure 2c,d show the potential energy landscapes when the amplitude of Q IR was set to ±5.00 Å μ along the x-and y-axis, respectively. This amplitude corresponds to the displacement of Mo atoms by ~0.23 Å and that of S atoms by ~0.34 Å in opposite direction along the polarization axis. Due to the deformation, the C 3v symmetry of the potential surface on the (Q LSx , Q LSy ) coordinates was broken, and the energy barriers along the three equivalent directions became different. Under the negative amplitude of Q IRx = −5.00 Å μ , the energy barrier along the r 1 direction decreased to 8.6 meV, but that along the other directions increased to 19.4 meV (r 2 direction) and 21.3 meV (r 3 direction), respectively. The energy landscape for the positive amplitude Q IRx = +5.00 Å μ is essentially the same with that for negative amplitude, except for the fact that the energy contour is flipped with respect to the mirror σ 1 . For both signs of Q IRx , the coordinate of the potential energy minimum is slightly shifted along the r 1 direction from the origin (Q LSx = Q LSy = 0 Å μ ).
By contrast, the amplitude of Q IRy largely increases the energy barrier along the r 1 direction (31.7 meV at Q IRy = −5.00 Å μ , and 23.8 meV at Q IRy = +5.00 Å μ ). This is accompanied by a slight shift of the potential minimum along the −r 1 (−y) direction. The energy landscape is symmetric with respect to the mirror plane σ 1 , and the energy barrier along the r 2 and r 3 directions are reduced (13.4 meV at Q IRy = −5.00 Å μ , and 10.1 meV at Q IRy = +5.00 Å μ ). The change in energy barriers and the shift of the potential minimum indicate that the coupling of Q LS and Q IR modes exerts an anharmonic force on the Q LS mode.

Dynamics of coupled normal-modes under light pulse. The dynamical behavior of the normal-modes
was investigated under light pulse with a specific polarization direction. Since the motion of the Q IR mode is much faster than that of the Q LS , the Q LS modes experience the effective potential asserted by the rapidly oscillating Q IR mode; i.e., time-averaged potential energy surface depending on Q IR (t). The dynamics of the nonlinearly coupled modes were simulated by the following coupled equations of motion,̈ where γ IR and γ LS are the damping coefficients for each mode, and F(t) is the optical driving force on the Q IR mode. We used Gaussian pulse F t Z E t e ( ) sin( ) / 2 , where E 0 is the amplitude of the electric field, σ is the duration of the pulse and Ω is the frequency. Figure 3 shows the evolution of potential energy curve on the Q LS coordinate along the AB → AC switching directions when the Q IR mode is resonantly pumped by a pulse with Ω = Ω IR , E 0 = 34 MVcm −1 and σ = 100 fs. Such high intensity of the pulse is required to achieve the large amplitude of Q IR (~5 Å μ ) in MoS 2, in order to explore the strong anharmonic coupling effect. The pulse intensity used in this study is comparable to the that used in the experiment on the high harmonic generation of the single layer MoS 2 32 . The energy curve (bold black line) corresponds to the static case (Q IRx = 0 Å μ ), where the energy of the AC stacking is slightly higher than that of the AB stacking because the layer-shearing by the Q LS mode is not perfectly rigid.
The pulse polarized along the x-axis induces the oscillation of the Q IRx mode with the amplitude between ± 5.35 Å μ by which the potential curve changes (grey line). The time-averaged potential energy (orange line) results in the effective barrier 〈Φ〉 = 11.5 meV along the r 1 direction (Fig. 3a), which is a significant reduction from the 17.3 meV for the static case. The coordinate of the potential minimum was shifted by 0.72 Å μ along the r 1 direction, and the energy of the AC stacking slightly increased compared to the static case. Meanwhile, a slight increase of the barrier to 19.2 meV along the r 2 and r 3 directions was observed (Fig. 3b). In contrast, the www.nature.com/scientificreports www.nature.com/scientificreports/ Q IRy mode under the y-polarized pulse shows asymmetric vibration between −4.69 Å μ and +5.43 Å μ in the anharmonic potential due to the lack of the mirror plane perpendicular to the y-axis. This results in an increase of the effective barrier along the r 1 direction to 22.5 meV, and the shift of the potential minimum by 0.36 Å μ along the −r 1 direction (Fig. 3c). On the other hand, the effective energy barrier along the r 2 and r 3 directions diminishes to 14.4 meV (Fig. 3d).
The polarization-dependent modulation of the effective potential energy can be explained by the characteristics of the anharmonic coupling terms. The overall trend is captured by the coupling terms in the form of Q Q determines the sign of F anh . The calculated F anh by the Q IRx has a positive value of 0.92 meVÅ −1 μ −1/2 , hence unidirectionally drives the Q LSy along the +y (r 1 ) direction (as indicated by the orange arrow in Fig. 3a). The F anh decreases the effective energy barrier along the r 1 direction by 33%, but increases the energy barrier by 9% along the r 2 and r 3 directions (according to the factor π = − F F cos(2 /3) (1/2) anh a nh ). Compared to Q IRx , the F anh from Q IRy is in the opposite direction with a slightly smaller magnitude (−0.87 meVÅ −1 μ −1/2 ). This explains the increase in the energy barrier along the r 1 direction by 30% and the decreases along the r 2 and r 3 directions by 16%.
The polarization-dependent direction of F anh has a geometrical origin related to the Mo and S sublattices, which are displaced by the Q IR in the opposite direction (Fig. 1c). The motion of Q IR modulates the interlayer interaction, which is approximated by the springs connecting the Mo and the S atoms (S 1 , S 2 , S 3 ) in the adjacent layers as illustrated in Fig. 3e. The associated interaction energy is ∑ ∆ = k d i i 1,2,3 2 , where k is spring constant and ∆d i is the change in distances between the Mo and S atoms. The contour plot of the interaction energy in Fig. 3f exhibits an anisotropy arising from the triangular geometry of the atomic arrangement. Particularly, the gradient of contour (orange arrow) indicates the force component along the +y direction when the Mo sublattice oscillates along the x-axis with respect to the S sublattice. On the contrary, the x-component of the force is canceled upon the rapid motion of Q IRx . This simple picture explains the F anh along the +y direction, and agrees with the selection rule (l + m = even) for Q Q Q IRx l LSx m LSy n coupling. In contrast, the Q IRy motion induces net forces along the −y direction due to the imbalance of the force (see the length of orange arrows).
Neither Q IRx nor Q IRy imparts such an unidirectional force on Q LSx since the relevant coupling terms (the Q Q IR l LSx with even l) are absent due to the odd parity of Q LSx . It prohibits the excitation of the interlayer shear along the r 2 and r 3 directions. Although the y-polarized light pulse lowers the energy barrier along the r 2 and r 3 directions, it cannot induce the ferroelectric switching along these directions. Therefore, the most effective way to realize the ferroelectric switching is to use the x-polarized light pulse which induces both the interlayer shear motion and energy barrier lowering along the r 1 direction for the ferroelectric switching to occur.
Next, we analyze the dynamics of ferroelectric switching based on the Q IRx -Q LSy coupling under the x-polarized pulse. First, we considered a case neglecting the damping of normal modes to simply show the essential consequences of the Q IRx -Q LSy coupling on the dynamics of Q LSy mode. Figure 4a,b show the motions of Q IRx and Q LSy modes at 0 K and 300 K, respectively, without damping under the x-polarized pulse with E 0 = 34 MVcm −1 and σ = 100 fs. The initial vibration amplitudes were set as the mean-square-displacement according to the Bose-Einstein distribution at each temperature. The initial vibration of Q LS was assumed to be aligned to the y-axis by setting the initial coordinate as Q LSx = 0 Å μ . This results in the initial amplitudes of Q IRx = 0.21 Å μ and Q LSy = 0.91 Å μ at 0 K, while Q IRx = 0.28 Å μ and Q LSy = 5.70 Å μ at 300 K, respectively. In Fig. 4a, the Q IRx and Q LSy modes oscillate with the harmonic frequencies before the arrival of the pulse at 0 ps. The initial kinetic energy of Q LSy mode was  Q /2 LSy 2 = 0.6 meV. When Q IRx was pumped, Q LSy started to oscillate with larger amplitude with respect to the shifted minimum at Q LSy = 0.87 Å μ (in good agreement with aforementioned 0.72 Å μ shift in the effective potential minimum). The pumping does not affect the motion of Q LSx (the value remains as ~0 Å μ ) as there are no forcing terms on it. The kinetic energy of Q LSy mode was increased to 1.6 meV by the anharmonic energy flow from the pumped Q IRx mode. It is noted that the pulse and Q LSy should be in-phase because the anharmonic force is unidirectional. The vibration of Q LSy is restricted in a small region because the kinetic energy is still smaller than the effective barrier of 〈Φ〉 = 11.5 meV under the oscillating Q IRx . In contrast, the oscillatory curve of Q LSy mode at 300 K in Fig. 4b shows slight modulations in shape and frequency by the onset of anharmonicity of Q LSy . The kinetic energy of Q LSy mode was 13.4 meV which is yet below the static energy barrier Φ 0 = 17.3 meV, but higher than the effective barrier 〈Φ〉 = 11.5 meV under the pulse. When Q IRx mode was pumped, Q LSy mode jumped over the barrier and oscillated with colossal amplitude between −3.79 Å μ and +19.19 Å μ . The vibration corresponds to the repetitive interconversion between AB and AC stackings, due to the absence of damping.
Secondly, a more realistic model that includes the damping of the normal modes was considered. The damping coefficients of γ IR and γ LS were taken as 2% of the harmonic frequencies, which are similar to the experimental values 21 . In Fig. 4c, the Q IRx and Q IRy modes initially oscillate with small amplitudes at 0 K until the arrival of the first pulse at 0 ps (damping is turned on at 0 ps). The eight sequencial pulses are applied to substantially amplify the Q LSy mode from the zero-point vibration at 0 K. The time interval between the subsequent pulses is gradually increased by ~4% from the 1/Ω LS ~ 1.6 ps for the phase matching between the pulse and Q LSy mode, considering the increase in the period of Q LS mode due to anharmonicity. Upon each cycle of pulse irradiation, the Q LSy mode is coherently amplified by gaining kinetic energy. After the eight pulses are irradiated, the Q LSy mode has sufficient kinetic energy and jump over the effective barrier which is reduced by the Q IRx mode. Once the initial AB stacking sequence changes to the AC stacking, it maintains the AC stacking due to the dissipation of the kinetic energies of the vibrations. This corresponds to the AB → AC ferroelectric switching. The opposite switching operation, AC → AB, can be performed by the same optical input as illustrated in Fig. 4d. The direction of F anh on Q LSy mode www.nature.com/scientificreports www.nature.com/scientificreports/ in the AC stacking is reversed (−r 1 direction) with respect to that (r 1 direction) in the AB stacking. The optical parameters of pulses (e.g. E 0 and σ) used in this study might be optimized further for more efficient switching, for instance, via pulse shaping techniques 37 .

Conclusion
In summary, the polarization switching mechanism of the bilayer 3R MoS 2 whose direction of the polarization is reversed by the change of the stacking sequence was investigated. The ferroelectric switching was achieved by driving the interlayer shear mode through the anharmonic energy flow from the optically pumped infrared mode. Remarkably, due to the selection rule from the crystal symmetry of MoS 2 , the degenerate interlayer shear mode can only be driven along its armchair axis whether the infrared mode is pumped along the zigzag or armchair axis. However, the optical pulse should be polarized along the zigzag axis for successful switching since the direction of anharmonic force is aligned with the switching direction. The coherent light pulses can amplify the interlayer shear mode substantially and unidirectionally, displacing the stacking sequence into the opposite polarization. The scheme for optical modulation of the stacking structure can be applied to other 2D materials exhibiting the interlayer shear mode to explore various stacking-dependent properties in a dynamical manner.

Methods
The density functional theory calculations were performed using Vienna Ab-initio Simulation Package (VASP) 38,39 . The projector-augmented wave (PAW) method 40 and a cut-off energy of 500 eV were used with the valence electron configurations of Mo[4s 2 4p 6 5s 2 4d4] and S[3s 2 3p4], respectively. The generalized gradient approximation 41 with Grimme's D3 scheme 42 was used to describe the van der Waals interaction. The bilayer structure was simulated by the supercell containing ~40 Å of vacuum layer to avoid the artificial interaction between periodic images. The structures were fully relaxed using 24 × 24 × 1 k-mesh until the residual forces on the atoms were less than 0.001 eVÅ −1 . Spontaneous polarization was calculated using Berry phase method 43 . The phonon calculation was performed using PHONOPY code 44 using 3 × 3 × 1 supercell and 8 × 8 × 1 k-mesh. Amplification of the Q LSy mode through the anharmonic force by the pulse at 0 K. (b) At 300 K, the Q LSy has sufficient kinetic energy to overcome the effective energy barrier (〈Φ〉 = 11.5 meV at |Q IRx | = 5.35 Å μ ) after the pumping, and oscillates back and forth between AB (Q LSy = 0 Å μ ) and AC (Q LSy = 16.29 Å μ ) stackings. (c) In the presence of damping, the Q LSy mode at 0 K overcomes the energy barrier after eight sequential pulses. The AB stacking changes to the AC stacking and does not return due to the dissipation of kinetic energy. (d) Schematics of ferroelectric switching through the Q IRx -Q LSy coupling. The orange arrows indicates the directions of interlayer shear induced by the x-polarized pulse, which are opposite in the AB and AC stackings.