Terahertz Driven Reversible Topological Phase Transition of Monolayer Transition Metal Dichalcogenides

Abstract This paper shows how terahertz light can drive ultrafast topological phase transitions in monolayer transition metal dichalcogenides (TMDs). The phase transition is induced by the light interaction with both electron and phonon subsystems in the material. The mechanism of such a phase transition is formulated by thermodynamics theory: the Gibbs free energy landscape can be effectively modulated under light, and the relative stability between different (meta‐)stable phases can be switched. This mechanism is applied to TMDs and reversible phase transitions between the topologically trivial 2H and nontrivial 1T′ phases are predicted, providing appropriate light frequency, polarization, and intensity are applied. The large energy barrier on the martensitic transformation path can be significantly reduced, yielding a small energy barrier phase transition with fast kinetics. Compared with other phase transition schemes, light illumination has great advantages, such as its non‐contact nature and easy tunability. The reversible topological phase transition can be applicable in high‐resolution fast data storage and in‐memory computing devices.


Introduction
The concept of topology in band structure has revolutionized condensed-matter physics and materials science. Different phases of materials can be classified by their electronic band topological indices, which serve as indicators of various exotic properties. [1] For example, the topologically protected surface/edge states of topological insulators (TIs) are promising platforms for fault-tolerant quantum computing and high-performance electronic or spintronic devices. Notably, some materials possess phases with contrasting topologies, and the topological phase transitions between these phases have evoked tremendous attention, owing to their scientific and technological importance. [2] Until now, one of the most widely studied materials that could exhibit phases with contrasting topologies is the 2D group-VI monolayer transition metal dichalcogenides (TMDs), in the chemical formula of MX 2 (M = Mo and W, X = S, Se, and Te). [3] Atomically, the M atom layer sits in the center, sandwiched by two X atom layers. Two (meta-)stable phases, 2H and 1T′, could exist in TMDs. In the 2H phase, the three atom layers X-M-X show an A-B-A stacking pattern in a hexagonal lattice (space group of P6m2, no. 187). If one X layer shifts about 1/3 lattice constant, the system would have an A-B-C stacking pattern, and is denoted as 1T phase. Such 1T phase is dynamically unstable and would spontaneously undergo a Peierls distortion, and finalizes to a distorted structure, 1T′ phase (space group of P2 1 /m, no. 11). Both 2H and 1T′ phases of monolayer group-VI TMDs possess exotic properties. In the 2H phase, TMDs are semiconducting with large bandgaps (≈1 eV) and trivial electronic topology. Remarkably, the lack of inversion symmetry in the 2H phase induces an out-of-plane effective "magnetic" field, which lifts the spin degeneracy of both the valence and conduction bands. This leads to valley-spin locking and a perfect valley contrasting circular dichroism. [4] On the other hand, in the 1T′ phase, TMDs are theoretically predicted and experimentally demonstrated to be 2D Z 2 -TIs, [5] which exhibit topologically protected nontrivial edge states.
Both 2H and 1T′ TMDs are under intensive research. The transitions between these two phases are attracting particular attention. Technologically, the transition between 2H and 1T′ phases only requires a shuffle of one X layer. Thus, it is a crystalline-tocrystalline, diffusionless, and martensitic phase transition with much faster kinetics than diffusive phase changes (such as that Figure 1. a) Schematic plot of THz irradiation on a monolayer 2H-TMD, which could transform into 1T′-TMD. b) Potential energy landscape of 2H, 1T′, and transition SP of TMD without THz irradiation. c) and d) are schematic energy profiles under different THz frequencies, which can well-control the relative stability of the two phases. Note that the energy barrier could be significantly reduced (or even eliminated) if THz laser intensity, polarization, and frequency are appropriately selected. e) Gibbs free energy change under light illumination in the polarization (P) general coordinate system. Thick and long dashed orange and green curves represent intrinsic (no light irradiation) 2H and 1T′ structures, respectively. Under alternating THz illumination (red waves), electric field reduces Gibbs free energy (thin and short dashed curves), and the phase transition occurs (solid orange and green curves). f) Atomic geometries of 2H, SP, and 1T′ phases. Atoms in the 2H phase are slightly tilted for clarity reason. Note that even though the Mo1 and Mo2 atoms are periodic images, they possess different bonding nature with the lower Te layer (denoted as Te1). The armchair and zigzag directions are denoted as x and y, respectively.
in GeSbTe alloys [6] ). Scientifically, such a phase transition could be an excellent platform to observe and investigate the process of trivial-to-nontrivial topological phase transitions. In addition, during the transition process, 1/6 of the M-X chemical bonds need to break and reform (see below). Hence, the energy barrier is sufficiently high (over 1 eV per MX 2 ), [5a] and the two phases are fairly stable with long lifetime once a thorough phase transition occurs. Therefore, the monolayer TMDs and their phase transition could provide a good platform to observe and study quantum phase transition and martensitic phase transition in lowdimensional materials, and may serve as potential data storage and in-memory computing material.
However, the high energy barrier also makes it challenging to manipulate the different phases of monolayer TMDs. Energetically, 2H is usually more stable than 1T′ (except for monolayer WTe 2 ). Previous theoretical works have predicted a few strategies to aid and facilitate phase transition from 2H to 1T′, including Li ion intercalation, [7] mechanical tension, [8] carrier injection, [9] electron-hole pair excitation, [10] and Ar-plasma bombardment. [11] Experimentally, visible laser irradiation on MoTe 2 was proposed to achieve the 2H to 1T′ phase transition, [12] but experiments thus far indicated that the new phase is a Temetalloid-like phase, rather than 1Tʹ. [13] Later, Wang et al. used ionic liquid gating to inject carriers into monolayer MoTe 2 and observed a 2H to 1T′ transition. [14] However, the phase transition easily reverses back once the gating is turned off. These facts also suggest that the phase transitions in monolayer MoTe 2 is challenging. The metastable 1T′-MoTe 2 can be directly fabricated via chemical vapor deposition method [15] by controlling Te vacancy concentration, but 1T′ to 2H phase transition could occur. [16] Most of the above-mentioned strategies require direct mechanical or electrochemical contacts between the sample and the tips or electrodes. Such contacts could introduce unwanted chemical impurities or additional interactions onto the monolayer TMDs; hence it may affect the sample quality and transition reversibility after many cycles. [17] In the current work, we propose a new mechanism that is free of these problems. We theoretically predict that terahertz (THz) frequency laser can trigger ultrafast phase transitions between 2H and 1T′ (Figure 1a). THz irradiation is a noncontact, noninvasive, athermic, and most-often nondestructive technique, [18] thus it could minimize lattice damages during operation. [19] In addition, it can be deeperpenetrating than visible lasers and the laser frequency, polarization, and intensity can be easily controlled, providing good flexibility. Here, we first theoretically demonstrate that under THz light irradiation, the potential energy landscape profile can be effectively altered, so that light irradiation would trigger phase transitions thermodynamically. To illustrate this theory, we apply first-principles density functional theory (DFT, see Experimental Section) to evaluate the electron and phonon contributions to the optical susceptibility of monolayer TMDs in the THz range, and show that THz can effectively tune and flip relative stabilities between 2H and 1T′ phases (Figure 1b-d). In addition, the energy barrier on the phase transition path can be greatly reduced with well-controlled THz frequency and intensity. In this case, ultrafast and reversible (2H↔1T′) topological phase transitions are expected. In this work, we mainly focus on monolayer MoTe 2 , because its energy difference between the 2H and 1T′ is the lowest among all monolayer TMDs. We also show that the THz driven phase transition is similarly applicable in other monolayer TMDs.

General Thermodynamic Theory
We will first illustrate our proposal with thermodynamic theory.
The optical alternating electric field can be written as In the low frequency range of interest ( ≈THz), the wavelength of the light (a few hundred µm) is much larger than the dimension of a unit cell, therefore the optical electric field can be regarded as almost homogeneous (wavevector k≅0). For a 2D material or surface of a 3D material, when the light propagates along z (normal to the material surface), ⇀  is in the xy-plane. Then, one could apply closed boundary condition [20] and study the Gibbs free energy with the electric field ⇀  as external forcing Here A denotes the unperturbed structure with coordinates 0 A , and it can be either 2H or 1T′ phases of TMD in our case. On the other hand, the internal Helmholtz free energy F has leading positive-quadratic dependence on , which is roughly 2 , where is the interatomic force constant. Such formalism also gives rise to the phonons. Then, the outcome of this standard linear-quadratic minimization problem, to the leading order in expansion, is which is illustrated in Figure 1e, the Gibbs free energy minimum position and value shifts linearly and quadratically with ⇀ , respectively. Since the induced electric polarization is Integrating this over the electric , where 0 is the vacuum permittivity, V f.u. is the volume of the formula unit, and ) (see Supporting Information for detailed discussions). The total Gibbs free energy is thus Note that microscopically the electric field linearly interacts with the polarization, while the macroscopic Gibbs free energy depend on the E 2 (Equation (3)), which is schematically plotted in Figure 1e. Since neither 2H nor 1T′ phase contains an intrinsic in-plane polarization, before light irradiation they both locate at | ⇀ P | = 0. Here ⇀ P serves as a generalized coordinate. When THz light is turned on, these two systems shift positions, depending on their optical susceptibility responses, and phase transition occurs when the 1T′ has lower Gibbs free energy than 2H. Note that if the system is in-plane ferroelectric (with finite static polarization ⇀ P s ) and is under a slowly alternating electric field (very small ), additional linear interaction term between ⇀ P s and ⇀  needs to be included. 21 From thermodynamic laws, in equilibrium states the total Gibbs free energy  E (A, ) should be minimized with respect to different phases A. Without laser irradiation, only the intrinsic part F 0 ( 0 A ) needs to be considered. For , that is, 2H-MoTe 2 is more stable than 1T′-MoTe 2 . When an external laser is applied, the second term in (Eq. 3) comes into play. Through the dependence of , and the relative stability between different phases can be switched. Note that ↔′ ( , ) corresponds to the Raman scattering process, during which photons from the laser field excites coherent phonons in the material system, and the material is thus driven out of equilibrium. As a result, with selected laser polarization, frequency, and intensity, an inverted order  E (1T ′ , ) <  E (2H, ) could be achieved. This indicates a topological phase transition under THz irradiation.

DFT Simulation Results on MoTe 2
The transition saddle point (SP) structure of MoTe 2 is obtained through cell-variable nudged-elastic band calculation. Geometrically, it almost remains the rectangular shape unit cell, and keeps the M y mirror symmetry as in both 2H and 1T′ phases. Compared with the 2H phase, which is the starting point of the phase transition, the transition SP lattice is elongated in the x-direction by 4% and shrinks in the y-direction by 5%, toward the 1T′ phase. As for the internal chemical bonding structure, we denote the relevant atoms as Mo1, Mo2, and Te1, as indicated in Figure 1f. Here the Te1 refers to the Te atom on the lower Te layer, and the Mo1 and Mo2 are periodic images of a single Mo atom in the simulation supercell (containing 2 Mo and 4 Te atoms). One could clearly see that from the 2H phase to the transition SP, the main atomic displacement is the movement of Te1 along the +x direction, and the chemical bond Mo1-Te1 breaks. Hence, the total energy of the system increases as the system approaches the SP. The Mo1 (Mo2) correspondingly moves to +x (−x), but the displacements are much smaller. After passing the transition state, the Te1 forms chemical bond with Mo2 and the system becomes 1T′ phase. This reduces the total energy. Therefore, during the whole process, 1/6 chemical bonds (one Mo-Te bond among all six bonds) in the system break and reform, which makes the energy barrier high enough to protect the 2H and 1T′ phases under ambient environment. Now we are ready to apply the above theory to monolayer MoTe 2 under linearly polarized terahertz laser (LPTL). According to Equation (3), we calculate the susceptibility ↔ ( ) at THz frequencies, which is contributed by both electron and ion subsystems. First, we explore the electron subsystem. The calculated electron band dispersions are shown in Supporting Information. One knows that the photon-electron interaction is mainly determined by two factors, namely, the joint density of states (jDOS) and the transition dipole matrix. The jDOS reads where the nk represents the energy of band-n at momentum k. c and v denote conduction and valence bands, respectively. The integral is taken in the first Brillouin zone (BZ). The jDOS of 2H and 1T′-MoTe 2 are plotted in Figure 2a. One clearly sees that the 2H-MoTe 2 possesses a much smaller jDOS than that of 1T′-MoTe 2 for ℏ < 1.5 eV, suggesting weaker optical response. We also plot the k-resolved jDOS (the integrand in Equation (4) We then examine the transition dipole matrix between the valance and conduction bands, | cv ii (k)| 2 = |⟨u ck |∇ k i |u vk ⟩| 2 , where |u nk 〉 is the periodic part of Bloch wavefunction. For 2H-MoTe 2 (left panel of Figure 2b), the transition dipole has a maximum value of 12 Å 2 around ±K and an average value of 2.4 Å 2 over the first BZ. On the contrary, for 1T′-MoTe 2 , the transition dipoles are significantly larger (middle and right panels of Figure 2b). One can see | cv xx (k)| 2 is mainly contributed from the states around ±Λ, with a maximum value of 386 Å 2 . | cv yy (k)| 2 have major contributions from both ±Λ and the −X ↔ +X path, with a maximum value of 260 Å 2 . Here x and y directions are perpendicular to and parallel with the dimerized Mo chain in the 1T′ structure, respectively. The average value of | cv xx (k)| 2 and | cv yy (k)| 2 of 1T′-MoTe 2 in the first BZ are 5.4 and 9.1 Å 2 , respectively. These suggest that the optical responses of 1T′-MoTe 2 should be stronger than that of 2H-MoTe 2 , consistent with previous studies that TIs have enhanced optical responses than normal insulators, due to band mixing between the chalcogen p-band and transition-metal d-band. [22] In addition, 1T′-MoTe 2 should have stronger optical response to the y-polarized LPTL than to the x-polarized LPTL.
The arguments above are verified by direct calculations of the susceptibility of 2H and 1T′-MoTe 2 . According to random phase approximation, the electron contributed susceptibility can be expressed as [23] el where ii ( ). For a 2D material computed with a 3D periodic supercell, the artificial contribution from vacuum space along the z-direction needs to be eliminated. In practice, the integration over the BZ is carried out by a k -mesh sampling, ∫ BZ where S is in-plane (x-y plane) area, h is the thickness of the 3D simulation supercell, and w k is the weight of each k point. For the in-plane components of the susceptibility tensor, according to a parallel capacitor model [24] one could use a scaling relationship ii ( ) are supercell calculated values and the rescaled values of the real part of the susceptibility tensor, respectively. d is an effective thickness of the 2D material, which is taken to be the distance between the centers of two layers when the same 2D material is van der Waals stacked in its bulk phase. Note that this d is a phenomenological quantity (similar as gauge choice in field theory). Fortunately, when we apply Equation (3) to estimate the LPTL-induced Gibbs free energy, the susceptibility is multiplied with the total volume V f.u. . Hence, the final results are independent of both h and d. As for the calculated imaginary part We plot the rescaled (with d = 7.7 Å) real part of electron contributed susceptibility and absorbance in Figure 2c,d. One observes that the absorbance of 2H-MoTe 2 starts at ≈1 eV, owing to its large direct bandgap at ±K points. Whereas in 1T′-MoTe 2 , the absorbance occurs at ≈0.1 eV due to its smaller direct bandgap. The absorbance for the y-LPTL is larger than that for the x-LPTL, consistent with previous analysis based on the transition dipole matrix. According to the Kramers-Kronig relation, if the absorbance occurs at lower energy with larger value, the real part susceptibility would be higher at low frequency. This agrees with our numerical results. One sees that at lowfrequency regime (e.g., THz, gray shaded areas in Figure 2c Next, we analyze the phonon contribution to the susceptibility. According to lattice dynamics theory and group theory analysis, only IR-active vibration modes can directly interact with LPTL and contribute to the first-order susceptibility. Also, since the wavevector of photons are much smaller than the dimensions of the BZ, only phonons near Γ-point need to be considered. We calculate the phonon dispersion (see Supporting Information), and label the optical phonon modes according to the irreducible representations of the corresponding crystalline symmetry group. For the 2H phase, one has and only the two dimensional E′ mode is IR-active. Since the 2H phase is inversion asymmetric, this E′ mode would split into LO and TO modes, where the LO mode contributed to optical responses. For the 1T′ phase, one has Γ 1T′ = 6A g ⊕3B g ⊕2A u ⊕4B u and the one dimensional B u and A u modes are IR-active under the x-and y-LPTL, respectively.
The phonon contributed susceptibility can be calculated via [25] ph is the oscillator strength of the m-th phonon mode (i,j = x,y). Here Z * ,ik is the Born effective charge on ion-, and u m ( k) is mass-normalized displacement of the m-th mode on the -th ion in the k-th direction. m and m are the frequency and lifetime of the near Γ-point phonons. The phonon lifetime comes into play because it influences the photon absorption linewidth and strength by IR-active phonons. When the frequency of the photon is close to that of a phonon, the photon can be absorbed by the phonon and then later dissipated into heat. According to the uncertainty principle, the absorption peak frequency and linewidth are m and 1/ m , respectively. In other words, the photon absorption mainly occurs in the frequency regime between m −1/ m and m +1/ m . Out of this window, photons are unlikely to be directly absorbed by the system ( ′′ ,ph ≅0). In this case, the photons can be treated as a pure alternating electric field with small dissipation. The phonon lifetime is dominated by three-phonon scattering (see Supporting Information), which originates in the third-order anharmonicity of the atomic potential. In Figure 3 we plot our calculated lifetime ( m ) of IR-active phonon modes. As temperature rises, the three-phonon scattering accelerates because more phonons are excited, thus phonon lifetime reduces. In detail, the lifetime of E′-LO mode in 2H-MoTe 2 is on the order of 10 1 ps at room temperature, suggesting a linewidth of the order of 0.1 meV. For 1T′-MoTe 2, the low frequency B u and A u modes (below 5 THz) has a much longer lifetime (thus a smaller linewidth).
Besides lifetime, another important factor that influences the phonon contributed susceptibility is the pattern of the ion displacements of vibrational modes. For example, a phonon can contribute to ph xx ( ) only when it is IR-active, and it vibrates strongly along the x-direction. The vibrational displacements of IR-active modes are also plotted in Figure 3b. One sees that the E′-LO of 2H phase has two 90°-rotation related degenerate modes that could couple to x-and y-LPTL. Three B u modes (at 3.74, 4.00, and 6.22 THz) of the 1T′ phase have small x-direction but large z-direction vibration components, so they do not contribute much to ph xx ( ).
Only the B u mode at 8.65 THz vibrates strongly along x. The A u mode of the 1T′ phase at 3.17 THz vibrates along both y and z, while the mode at 5.55 THz is mainly along y.
We can calculate the phonon contributed susceptibility. Note that the phonon lifetime is temperature dependent. Hence, we plot the results under several selected temperatures, in Figure 4. It is clearly seen that the imaginary parts are composed of Lorentzian-shape peaks sitting on the frequencies of IR-active phonon modes. As temperature increases, the absorption peaks get shorter and wider because the phonon lifetime reduces (linewidth increases). For the real part of susceptibility, such temperature variation only affects their values near the resonant regime (| − m | ≲ 1/ m ). Out of this regime (| − m | ≳ 1/ m ), the real part susceptibility response remains nearly unchanged as temperature varies. At low frequency (lower than first absorbance peak), the calculated susceptibilities are  Figure 4. When the phonon linewidth is small (on the order of 0.01 THz), the absolute value of ′,ph can be as large as a few hundred around the absorption peak. This offers us opportunities to manipulate the potential energy landscape effectively when the THz frequency is carefully selected, according to Equation (3). Now we add the electron and phonon contributed susceptibilities together (Figure 5a). One can see that for < 3 THz, the maximum values of ′ yy (1T ′ , ) and ′ xx (1T ′ , ) are 60.9 and 34.5, respectively. While for 2H-MoTe 2 , one has ′ xx (2H, ) = ′ yy (2H, ) = 19.5. If one uses a y-LPTL of = 1 THz, the susceptibility difference between the 1T′-MoTe 2 and 2H-MoTe 2 is then 41.4. Hence, from Equation (3), such a laser tends to push the system toward 1T′ phase and stabilize it. When the laser intensity in free space I reaches 8.6 × 10 10 W⋅cm −2 (E = 0.80 V⋅nm −1 , from I = 1 2 0 cE 2 ), the 1T′-MoTe 2 would have lower Gibbs free energy than that of 2H-MoTe 2 , as shown in Figure 5b. This would lead to a thermodynamic phase transition from 2H to 1T′. If we move closer to the frequencies of IR-active phonon modes of 1T′-MoTe 2 , the susceptibility difference between 2H and 1T′ phases would dramatically increase. For example, at = 5.45 THz, ′ yy (2H, ) and ′ yy (1T ′ , ) are 21.5 and 519.0, respectively. Hence a = 5.45 THz y-LPTL with its laser intensity as small as 7.6 × 10 9 W⋅cm −2 (E = 0.24 V⋅nm −1 ) is sufficient to make the 1T′ phase more stable than 2H. Our result agrees well with the recent experimental observations [26] that a 0.3−1.5 THz (centered at 0.5 THz) irradiation could drive a 2H→1T′ phase transi-tion in monolayer MoTe 2 . The free space field amplitude used in the experiment is 0.03 V⋅nm −1 , which is further field enhanced by 20−50 times in most regions; hence the field amplitude irradiated onto the samples is around 0.6−1.5 V⋅nm −1 , consistent with our estimates.
The above prediction suggests that y-LPTL irradiation triggers ion movement along x, which is the 2H→1T′ structural change path, which is counterintuitive. In order to confirm and visualize this effect in real space, we perform an ab initio molecular dynamics simulation by applying an alternating electric field along y. The simulation results clearly indicate an ion equilibrium position movement along x. Around the new equilibrium state, the ions vibrate with m = 7.1 THz (see Supporting Information for details). These are consistent with macroscopic thermodynamic theory prediction.

Thermal Effect under THz
We estimate the temperature rise of the system under THz irradiation. When the laser frequency is chosen out of direct photon absorbance regime ( ′′≅0), the energy transfer from photon absorption can be omitted. We only consider optical alternating electric field effect here. First, before phase transition, the system oscillates according to the driven oscillator equation where x µ , m µ , and Z * ,xy are average values of displacement, mass, and Born effective charge of ion-. m = 7.1 THz and are the phonon frequency and its effective damping, and E 0 = 1 V⋅nm −1 and = 1 THz are the THz light magnitude and frequency, respectively. The displacement amplitude is then The amplitudes of Mo and Te are estimated to be 0.0085 and 0.0031 Å, respectively, well consistent with ab initio molecular dynamics simulation results (Supporting Information). We use the heat capacity of MoTe 2 (18.4 cal⋅K −1 ⋅mol −1 ) [27] and estimate the temperature rise when these driven oscillating motions fully convert into heat, c p T = 1 2 ∑ m 2 m (x max ) 2 , which gives 6 K. In addition, the work done under THz irradiation during phase transition could also dissipate into heat, which is W = The temperature rise from this source is similarly estimated to be 8 K. Therefore, during THz driven phase transition of MoTe 2 , the temperature rise effect is small enough and can be dissipated out quickly. The waste heat problem in conventional phase change materials can be eliminated here.

Reduced Energy Barrier Phase Transition
In order to trigger fast phase transitions, one has to evaluate the kinetics along the transition path. We use the SP structure and its optical responses (see Supporting Information) to estimate the energy barrier change under LPTL. Without laser irradiation, the energy barrier from 2H to 1T′ is found to be 1.06 eV per f.u. (corresponding to 160 µJ⋅cm −2 ). We find that in the THz range, the electron contributed susceptibility of the SP is The total susceptibility of the transition SP structure is plotted in Supporting Information. Hence, even though a y-LPTL with = 1 THz and E = 0.80 V⋅nm −1 is large enough to induce a phase transition from 2H to 1T′ thermodynamically, according to Equation (3), the energy barrier is still 0.97 eV per f.u. (corresponding to 145.5 µJ⋅cm −2 ). One can increase its intensity to I = 1.1 × 10 12 W⋅cm −2 (E = 2.84 V⋅nm −1 ), and the energy difference between the SP and the 2H phase may be wiped out. However, thorough evaluation the free energy landscape under THz irradiation is a high dimensional problem. Here, this Gibbs free energy comparison is a simple estimate. Nevertheless, we expect that under this THz intensity, the energy barrier can be significantly reduced. According to the Arrhenius law, phase transition kinetics exponentially increases as the energy barrier reduces, hence it corresponds to a fast transition under such intense light irradiation. Considering the large contrasting optical susceptibilities between these phases, it is possible to erase the energy barrier, once the optical frequency and intensity are carefully selected, but an exact prediction is very challenging. If the energy barrier is completely erased, the topological phase transition could happen on the timescale of a few picoseconds. [28] Even in this situation, the system needs additional time to dissipate its kinetic energy and fully relax to the new equilibrium state (1T′), which could occur on the order of tens of picoseconds. [29] Note that this is a rough estimate, and the dissipation time also depends on the boundary condition of the material. In addition, one has to note that our model is based on a unit cell calculation, which corresponds to a coherent phase transition in the whole sample. In reality, the material may contain different domains, and the domain boundary energy (on the order of a few tens of meV⋅Å −1 ) [8] could change the free energy profile. Then the nucleation barrier associated with the domain boundary needs to be overcome during phase transition. In order to further explore the phase transition kinetics and its relationship with light frequency and strength, more experiments with stringent condition need to be conducted. In the recent experiment, [26] the field acting onto the sample is around 0.6-1.5 V⋅nm −1 . According to our model, this can still have an energy barrier, hence, a transition timescale on the order of nanosecond, as observed in the experiment, is quite reasonable. Even though the 2H phase is energetically more stable without any external stimulus, once a large area of 1T′ phase is fabricated, it is difficult to overcome the high energy barrier and transit back to 2H. This non-volatile phase transition guarantees data safety for long-term storage. As for data writing, one needs another stimulus to boost a 1T′→2H phase transition. This can be done by selecting the frequency regime where ′(2H, ) is much larger than ′(1T′, ). From Figure 5a, one sees that a y-LPTL with = 5.63 THz can be used. This is slightly above the IR-active mode of 1T′-MoTe 2 at = 5.55 THz. At this frequency, ′ yy (2H, ), ′ yy (SP, ), and ′ yy (1T ′ , ) are 21.5, 89.9, and −305.1, respectively. Hence, as the y-LPTL intensity increases, the 2H and 1T′ energy difference is also enhanced (Figure 5c). Particularly, a LPTL with I = 3.0 × 10 11 W⋅cm −2 (E = 1.53 V⋅nm −1 ) is sufficient to greatly reduce the transition energy barrier. This completes the ultrafast topological phase transition cycle (2H↔1T′) in monolayer MoTe 2 . The optical susceptibility and Gibbs free energy for typical frequencies and electric field strengths are listed in Table 1.
2D materials with atomic-scale thickness in the out-of-plane direction may reduce the data unit volume (L × L × d in the real space), compared with conventional 3D materials (L × L × L). Here L is characteristic lateral size (usually a few tens to a few hundred nanometers) and d is effective thickness (subnanometer). In addition, the topological metallic edge state in 1T′ TMD monolayers is a few nanometers, [5a] which allows minimum horizontal feature size to be a few tens nanometers. In this way, the single data storage unit of 2D monolayer MoTe 2 may occupy smaller space than 3D bulk materials, and the information storage density (over the real space) can be high. Since the 2H-MoTe 2 is a semiconductor with sizable bandgap, while the 1T′-MoTe 2 has a small topological bandgap, the electric resistance and optical reflection could be highly distinct in the two phases. [30] Thus, one may use electrical or optical approach to distinguish the 2H and 1T′ phases. If the sample quality is good and the phase transition occurs completely, the data reading resolution should be high enough for practical applications. We thus propose that the monolayer MoTe 2 can be used as potential memory devices with large information storage density.

Other Monolayer TMD Systems
We also compute the optical responses of other monolayer TMDs with small energy difference and transition barrier, as shown in Figure 6. One could see that in both MoSe 2 and WSe 2 , the ′ yy (1T ′ , ) is higher than the other two in most frequency regimes. Thus, a y-LTPL is able to drive a 2H→1T′ phase transition in these systems. Specifically, if the y-LTPL frequency is selected to be 4.4 THz (for MoSe 2 ) and 4.5 THz (for WSe 2 ), intensity of 1.3 × 10 11 (MoSe 2 ) and 9.2 × 10 10 W⋅cm −2 (WSe 2 ) would be enough for the transition, respectively. If the frequency is low (below the IR-active mode, e.g., 2 THz), then the required intensity would be higher, which is 1.2 × 10 12 W⋅cm −2 for MoSe 2 and 2.8 × 10 11 W⋅cm −2 for WSe 2 . On the other hand, if one wants to facilitate phase transition from 1T′ to 2H, the frequency has to be chosen carefully (4.6 THz for both MoSe 2 and WSe 2 ). As for the monolayer WTe 2 , since its ground state is 1T′, we can apply a y-LPTL with = 3.4 THz and I = 6.4 × 10 10 W⋅cm −2 to drive a transition from 1T′ to 2H. Comparing Figures 5a and 6, the monolayer MoTe 2 might be the best platform to realize ultrafast topological phase transition between the 2H and 1T′ phases back and forth, because it allows wider frequency windows that have both large ′ yy (2H, ) − ′ yy (1T ′ , ) and ′ yy (1T ′ , ) − ′ yy (2H, ).

Discussion
One has to note that defects and disorders always exist in real samples. If the defects arise in the material with a considerably high concentration, they may strongly affect the optical response functions, by significantly reducing the carrier lifetime and introducing localized doping levels in the phonon and electron band structures. [31] For example, if the electronic doping level is shallow, it will enhance the THz optical absorption, and the unwanted optical loss would occur. In this case, the imaginary part of optical susceptibility needs to be considered. Actually, the optical absorption could also change the energy landscape and trigger phase transition. [10] When the defect or impurity concentration is not high enough, we can phenomenologically incorporate its effect by adopting finite lifetime in the response formulae, Equations (5) and (6). The current theoretical work is based on the perfect crystal, and a complete evaluation on defects and impurities is out of scope here. We will discuss phase transitions in defective systems elsewhere. There are several 3D bulk materials possessing different geometric phases with contrast topology. For example, bulk SnSe could have topologically trivial semiconducting phase Pnma and nontrivial phase Fm3m under ambient environment, and the phase transformation between them is also martensitic. [1b] The present thermodynamic theory is applicable to these systems, as long as their spontaneous polarization P s is zero. If the system has finite P s , additional terms need to be considered, as the THz oscillating electric field may flip P s over time. As for the dimensionality effects, the free-standing 2D materials are sandwiched by vacuum in the out-of-plane z direction, while the conventional 3D materials usually lack such free space. One has to note that during phase transformation, the associated transformation strain can be a few percent. For 2D materials, the transformation stress can be easily released by z-motion, if they are freely suspended with slight pre-buckling, so that long-range elastic energy penalty can be minimal. In contrast, the accommodation of martensitic transformation strain in 3D bulk materials is a well-known materials science problem, [32] that causes long-range elastic interactions and generally requires several orientation variants to cooperate, [33] which constrains the speed and versatility. The residual back-stress would also limit the spatial extent of the transformation and reduces the transformation reversibility. Hence, martensitic transitions in 2D materials can be greatly advantageous compared to 3D materials in speed and reversibility.

Conclusion
In summary, our theory on THz driven topological phase transition provides another route for optomechanical manipulations of the materials phase, by incorporating both photon-electron and photon-phonon coupling. In particular, this mechanism can be a novel way to realizing and understanding the topological phase transitions in monolayer TMDs under THz illumination, which can be potentially used in next-generation data storage and computing devices based on 2D materials. In addition, such an optical scheme can be widely applicable to other materials for ultrafast martensitic transformations.

Experimental Section
Density Functional Theory: First-principles calculations were based on DFT with spin-orbit coupling included self-consistently, as implemented in the Vienna ab initio simulation package. [34] The solidstate Perdew-Burke-Enzerhof [35] functional was adopted to evaluate exchange-correlation potential. Projector-augmented wave method [36] and planewave basis set were used to treat the core and valence electrons, respectively. Planewave cutoff energy was set as 350 eV. In order to eliminate the artificial image interactions under 3D periodic boundary condition, a vacuum space of 18 Å was adopted in the z-direction. Monkhorst-Pack k-mesh [37] of (15 × 15 × 1) and (9 × 15 × 1) grids were used to optimize the unit cell geometry of 2H and 1T′ phase, and meshes of (45 × 45 × 1) and (25 × 45 × 1) were used to calculate the electron contributed frequency-dependent susceptibility ↔ el ( ). Total energy and force convergence criteria were set to be 1 × 10 −7 eV and 0.001 eV Å −1 , respectively. The convergence of these parameters were well tested. Phonon dispersion and zone center phonon mode lifetime calculations were performed using density functional perturbation theory (DFPT) [38] as implemented in the Phonopy [39] and Phono3py [40] codes. LO-TO splitting effect was included, and the LO frequency near the Γ-point was used to evaluate optical responses. Note that finite displacement method yielded similar results as in DFPT.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.