ABSTRACT Displacive martensitic phase transition is potentially promising in semiconductor-based data storage applications with fast switching speed. In addition to traditional phase

transition materials, the recently discovered two-dimensional ferroic materials are receiving a lot of attention owing to their fast ferroic switching dynamics, which could tremendously

boost data storage density and enhance read/write speed. In this study, we propose that a terahertz laser with an intermediate intensity and selected frequency can trigger ferroic order

switching in two-dimensional multiferroics, which is a damage-free noncontacting approach. Through first-principles calculations, we theoretically and computationally investigate optically

coupled with the transverse optical phonon modes. The transition from one orientation variant to another can be barrierless, indicating ultrafast transition kinetics and the conventional

nucleation-growth phase transition process can be avoidable. SIMILAR CONTENT BEING VIEWED BY OTHERS CHARACTERIZING G-TYPE ANTIFERROMAGNETISM QUANTITATIVELY WITH OPTICAL SECOND HARMONIC

GENERATION Article Open access 22 April 2025 LIGHT-DRIVEN ELECTRODYNAMICS AND DEMAGNETIZATION IN FE_N_GETE2 (_N_ = 3, 5) THIN FILMS Article Open access 12 November 2024 SHEAR-STRAIN-MEDIATED

LARGE NONVOLATILE TUNING OF FERROMAGNETIC RESONANCE BY AN ELECTRIC FIELD IN MULTIFERROIC HETEROSTRUCTURES Article Open access 22 January 2021 INTRODUCTION Optical control of material

geometries is receiving rapidly growing attention in very recent years. For example, some ferroelectric/multiferroic perovskites, such as BaTiO3 (refs 1,2), BiFeO3 (ref. 3), and SrTiO3 (refs

4,5) would experience their ferroic order switch under laser pulse illumination (with the optical frequency well below their electronic bandgaps). It offers a remarkable and promising

scheme to manipulate the structure of materials, avoiding mechanical or electrochemical contacts with these samples, which might slow down the effect and introduce unwanted impurities or

disorders. Thus, this fast-paced noncontacting opto-mechanical strategy is less susceptible to lattice damage and provides an ultrafast manipulation of materials on picosecond time scales

and sub-micrometer length scales6. According to optical response theory, there are four regimes of optical (angular) frequency, namely, low frequency regime where optical frequency \(\omega

\,<\, \omega _0 - \frac{1}{\tau }\) (_ω_0 is energy difference between two optical-transition allowed states and _τ_ is lifetime), absorption frequency regime where \(\omega _0 -

\frac{1}{\tau } \,<\, \omega \,<\, \omega _0 + \frac{1}{\tau }\), reflection regime where \(\omega _0 + \frac{1}{\tau } \,<\, \omega \,<\, \omega _{\mathrm{p}}\) (where _ω_p is

plasmon frequency), and transparent regime where _ω_ > _ω_p. The second frequency regime may introduce significant heat, the light is highly reflected in the third frequency regime, and

the optical response in transparent regime is usually small. Hence, low optical energy light would be intriguing to control the material behaviors practically7. Low optical energy light can

be further classified into two categories: near or mid-infrared optics with its photon energy above phonon but below electronic bandgap of semiconducting materials (such as experiments on

BaTiO3, refs 1,2), and far infrared optics with its (terahertz) frequency strongly and directly coupled with phonons. Atomically thin two-dimensional (2D) materials, with their extremely

large surface-area-to-volume ratio, are more optically addressable and accessible. Therefore, noncontacting optically driven ferroic order transition in 2D ferroic materials will be

promising with their easy and damage-free manipulation, large information storage density, and ultrafast kinetics. Previous theoretical studies mainly use parameterized model including

phonon–phonon nonlinear interactions8,9. In this work, we theoretically and computationally evaluate the interactions between light and phonons, as well as light and electrons in 2D

time-reversal invariant multiferroic (with both ferroelectric and ferroelastic orders) materials, from a first-principles approach. We choose two experimentally fabricated systems, i.e.,

monolayers _β_-GeSe10 and _α_-SnTe11 which represent the mostly observed group-IV monochalcogenide monolayer family, to illustrate our theory. We predict that linearly polarized terahertz

laser (LPTL) pulses with intermediate intensity (around 1–2 MV cm−1) can trigger ferroic order switch in these systems. A contact-free direction-dependent vibrational electron energy loss

spectroscopy (EELS) is also theoretically calculated, which can be used to detect and measure the structural signature in a high resolution. In addition, we calculate the second harmonic

generation (SHG) signals, which also couple to ferroic orders. These strategies could provide powerful and non-invasive tools to characterize ferroicity that is indistinguishable by the

traditional diffraction methods. We analyze the light–matter interaction thermodynamically. When LPTL is illuminated onto a semiconductor (with optical bandgap larger than ~40 meV), the

electron-hole pair formation is eliminated owing to low photon energy. Therefore, only optical electromagnetic field effects need to be considered. Here we are interested in the

time-reversal invariant systems, in which the magnetic field interaction is very weak and can be omitted. Hence, we focus on the alternating electric field component (\({\mathbf{{\cal{E}}}}

= {\mathbf{E}}_{{\mathrm{max}}}e^{i\omega t}\), _ω_ ~ THz) effect, here Emax is the maximum value (amplitude) of the electric field. The LPTL would accelerate both electronic and ionic

subsystems in the material, and its work done per volume can be evaluated by \(du = {\mathrm{Re}}\left( {{\mathbf{{\cal{E}}}} \cdot d{\mathbf{{\cal{P}}}}^ \ast } \right)\), where

\({\mathbf{{\cal{P}}}}\) is the time-dependent polarization density. Here, closed boundary condition12,13 is applied, since electric polarization is in the _xy_-plane. Using Legendre

transformation, the Gibbs free energy (GFE) density change is then \(dg = - {\mathrm{Re}}\left\langle {{\mathbf{{\cal{P}}}}^ \ast \cdot d{\mathbf{{\cal{E}}}}} \right\rangle\), and

\(\left\langle \cdot \right\rangle\) indicates time average. Note that \({\mathbf{{\cal{P}}}} = {\mathbf{P}}_{\mathrm{s}} + \varepsilon _0{\mathrm{Re}}\mathop{\chi}\limits^{\leftrightarrow}

\left( \omega \right) \cdot {\mathbf{{\cal{E}}}}\), where Ps is spontaneous electric polarization, _ε_0 is vacuum permittivity, and \(\mathop{\chi}\limits^{\leftrightarrow} \left( \omega

\right)\) is optical susceptibility tensor at the frequency _ω_, containing electronic and phononic contributions (\(\mathop{\chi}\limits^{\leftrightarrow} =

{\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}}} + {\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{ph}}}\)). Under LPTL illumination (along the _i_-direction), note that the light

frequency is on the THz order, which could be lower than the phonon Debye frequency of the system. Hence, if _E_max is large enough (much stronger than the coercive field _E_c to reverse

polarization), the spontaneous polarization may follow the electric field and switch back and forth between Ps and −Ps. Here _E_max indicates the magnitude of electric field vector, |Emax|.

When the _E_max is much stronger than _E_c (corresponds to the situation in our current discussion), we can show that it contributes to the time-averaged GFE density in the form of (see

Supplementary Note 1 for details) $$g_0 = - \frac{{{E}_{{\mathrm{max}}}{P}_{\mathrm{s}}\cos \theta }}{2}$$ (1) where \(\theta \in [0,\frac{\pi }{2}]\) is the angle between the LPTL

polarization direction and the spontaneous polarization direction. This is a first-order interaction that measures between light and polarization. If the light frequency is highly above the

phonon frequency range (e.g. several tens THz and above), then the ion position change cannot follow the light field (off-resonant). Then this term would become zero as the

\({\mathbf{P}}_{\mathrm{s}}\) is time-independent. We also include the second-order interactions by incorporating optical susceptibility. This can be considered by the averaged GFE density

change through integrating time and electric field in \(dg = - \varepsilon _0{\mathbf{{\cal{E}}}}^ \ast (t) \cdot {\mathrm{Re}}[

{{\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}}}\left( \omega \right) + {\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{ph}}}\left( \omega \right)} ] \cdot

d{\mathbf{{\cal{E}}}}(t)\), and the total GFE density contributed from optical linear responses can be written as $$\begin{array}{ccccc}\\ g_1\left( \omega \right) = & -

\frac{1}{4}\varepsilon _0{\mathbf{E}}_{{\mathrm{max}}} \cdot {\mathrm{Re}}\left[ {{\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}}}\left( \omega \right) +

{\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{ph}}}\left( \omega \right)} \right] \cdot {\mathbf{E}}_{{\mathrm{max}}}\cr \\ & = - \frac{1}{4}\varepsilon _0{\mathrm{Re}}\left[ {\chi

_{ii}^{{\mathrm{el}}}\left( \omega \right) + \chi _{ii}^{{\mathrm{ph}}}\left( \omega \right)} \right]{E}_{{\mathrm{max}},i}^2\\ \end{array}.$$ (2) Note that the integration of \(g_1\) over

electric field space gives one 1/2 factor, and the time average of sinusoidal electric field gives another 1/2. Hence, there is a 1/4 coefficient in Eq. (2). The total GFE density change

under light can then be estimated by \(g = g_0 + g_1\). In a ferroic material with multiple ferroic orders, different orders are subject to a simple structural operator (such as a proper or

improper rotation \(\hat {\cal{O}}\)), and hence are energetically degenerate. Certain directional LPTL could break such degeneracy (Eqs. 1 and 2), and one can expect optically induced

ferroic order switch from a metastable order with higher GFE to a stable order with lower GFE. This mechanism is illustrated in Fig. 1. The optical susceptibility can be evaluated via ab

initio density functional theory (DFT) calculations (see “Methods” for details14,15,16,17,18,19,20,21). We use the Vienna ab initio simulation package15 to calculate the electron behaviors

self-consistently, and compute forces on ions to evaluate phonon dispersions using density functional perturbation theory (DFPT), as implemented in the Phonopy package21. According to the

linear response theory with random phase approximation (RPA), the electronic part of susceptibility \({\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}}}\left( {{\mathbf{q}} = 0,\omega

} \right)\) is the Fourier transformation of real space density response function \({\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}}}\left( {{\mathbf{r}},{\mathbf{r}}^{\prime},\omega

} \right)\), which is solved from a Dyson equation $${\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}}}\left( {{\mathbf{r}},{\mathbf{r}}^{\prime},\omega } \right) =

{\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}},\left( 0 \right)}\left( {{\mathbf{r}},{\mathbf{r}}^{\prime},\omega } \right) + {\int}

{d{\mathbf{r}}_1d{\mathbf{r}}_2{\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}},\left( 0 \right)}\left( {{\mathbf{r}},{\mathbf{r}}_1,\omega } \right)\frac{1}{{\left| {{\mathbf{r}}_1 -

{\mathbf{r}}_2} \right|}}{\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}},\left( 0 \right)}\left( {{\mathbf{r}}_2,{\mathbf{r}}^{\prime},\omega } \right)},$$ (3) where

\({\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}},\left( 0 \right)}\) is bare density response function (independent particle) contributed by electron transitions

$${\mathop{\chi}\limits^{\leftrightarrow}}^{{\mathrm{el}},\left( 0 \right)}\left( {{\mathbf{r}},{\mathbf{r}}^{\prime},\omega } \right) = \mathop {\sum}\limits_{m \, \ne \, n} {\left( {f_n -

f_m} \right)\frac{{\varphi _m^ \ast \left( {\mathbf{r}} \right)\varphi _n^ \ast \left( {{\mathbf{r}}}^{\prime} \right)\varphi _n\left( {\mathbf{r}} \right)\varphi _m\left(

{{\mathbf{r}}}^{\prime} \right)}}{{\hbar \omega - {\it{\epsilon }}_m + {\it{\epsilon }}_n + i\xi }}}.$$ (4) Here, _f__i_, _ϵ__i_, _φ__i_ are the Fermi-Dirac distribution, energy, and

wavefunction of the _i_-th level, and _ξ_ is the Lorentzian phenomenological damping parameter (taken to be 0.025 eV in our calculations), representing disorder, finite temperature, and

impurity effects. When the frequency of LPTL lies in the range of phonon frequency (a few THz), vibrational phononic contribution to optical susceptibility needs to be included. According to

lattice dynamics theory22, the LPTL couples with infrared-active transverse optical (TO) mode of phonons at the Γ-point of the Brillouin zone. The phononic contributions to susceptibility

is calculated according to23 $$\chi _{\alpha \beta }^{{\mathrm{ph}}}\left( \omega \right) = \frac{1}{{\Omega }}\mathop {\sum}\nolimits_m {\frac{{{\cal{Z}}_{m,\alpha }^ \ast

{\cal{Z}}_{m,\beta }^ \ast }}{{\omega _m^2 - \left( {\omega + i{{\Gamma }}} \right)^2}}},$$ (5) where Ω is the unit cell volume, _α_ and _β_ (= 1, 2, 3) are the Cartesian-coordinates

indices, and _m_ is the TO phonon mode index. The Born effective charge of each phonon mode is evaluated as \({\cal{Z}}_{m,\alpha }^ \ast = \mathop {\sum}\nolimits_{\kappa ,\alpha ^{\prime}}

{Z_{\kappa ,\alpha \alpha ^{\prime}}^ \ast u_{m,\kappa \alpha ^{\prime}}}\), where _κ_ is the ion index, _u_ is the mass normalized displacement, and \(Z_{\kappa ,\alpha \alpha^{\prime}}^

\ast\) is the Born effective charge component on each ion. The numerator accounts for the vibrational mode oscillator strengths. Finally, Γ is linewidth parameter accounting for the finite

lifetime (_τ_ = Γ−1) of the vibration. For simplicity, we choose a universal value of Γ = 4 cm−1, which is comparable with the linewidth determined by phonon–phonon couplings24,25. RESULTS

MONOLAYER _Β_-GESE We now apply these analyses to two experimentally fabricated 2D multiferroic materials. Figure 2a shows the atomic structure of _β_-GeSe monolayer, with relaxed lattice

constants to be _a_ = 3.59 Å and _b_ = 5.73 Å. This structure belongs to _Pmn_21 space group (no. 31) without centrosymmetry, consistent with previous works26,27. One clearly observes that

it looks like a distorted honeycomb lattice (such as h-BN, compressed along the _y_-direction). Actually, the honeycomb structure is serving as a parental phase (Supplementary Note 2), and

the _β_-GeSe shown in Fig. 2a is one of its ferroelastic orientation variants (denoted as FE0). The (_a_ × _b_) rectangle unit cell in _β_-GeSe corresponds to the (1 × √3) supercell of the

high symmetric parental structure, from which a spontaneous structural transformation occurs, forming different orientation variants. Similar as the 1T and 1Tʹ phases in transition-metal

dichalcogenide monolayers28, there are three symmetrically equivalent ferroelastic _β_-GeSe (FE0, FE1, and FE2, subjecting to 120°-rotations). The 2D transformation strain tensor of these

orientation variants are \({\mathop{\eta}\limits^{\leftrightarrow}}_0 = \left( {\begin{array}{*{20}{c}} {0.042} & 0 \\ 0 & { - 0.040} \end{array}} \right)\),

\({\mathop{\eta}\limits^{\leftrightarrow}}_1 = \left( {\begin{array}{*{20}{c}} { - 0.016} & { - 0.041} \\ { - 0.041} & {0.021} \end{array}} \right)\), and

\({\mathop{\eta}\limits^{\leftrightarrow}}_2 = \left( {\begin{array}{*{20}{c}} { - 0.016} & {0.041} \\ {0.041} & {0.021} \end{array}} \right)\). Therefore, the switch strain from one

orientation variant to another is small, which could occur under intermediate energy input in experiments (Supplementary Note 2). Since they are symmetrically equivalent, we will only

calculate physical properties of the FE0. One could perform rotational operation for the other two orientation variants. In addition to ferroelastic order, the lack of centrosymmetry

indicates a spontaneous electric polarization Ps. Consistent with previous works27, we compute its value to be 0.16 nC m−1 (see Supplementary Note 3 for details), along the armchair (_y_−)

direction (Fig. 2a). This polarization is switchable under intermediate in-plane static electric field. Thus, it also possesses a ferroelectric order, making it a 2D time-reversal invariant

multiferroic material. Next, we calculate the electronic and phonon band dispersions (Supplementary Note 4), based on which we can compute the optical susceptibility of FE0-_β_-GeSe,

according to Eqs. (3) and (4). Note that all calculations are performed in a 3D periodic supercell with a large vacuum space separating different images. In order to eliminate the vacuum

effects and obtain 2D values, we adopt a widely used scaling method for the real part susceptibility, \({\mathop{\chi}\limits^{\leftrightarrow}}_{||}^{3{\mathrm{D}}}d =

{\mathop{\chi}\limits^{\leftrightarrow}}_{||}^{2{\mathrm{D}}}h\), according to a parallel capacitor model29,30,31. Here, superscripts “3D” and “2D” refer to susceptibility in the supercell

and in the 2D form, || indicates that only _xy_-plane component can be scaled, and _d_ and _h_ are simulation supercell and 2D material thickness, respectively. We use the separation

distance in the van der Waals stacked bulk structure to estimate _h_, which gives 5.5 Å. The calculated \({\mathop{\chi}\limits^{\leftrightarrow}}_{||}^{2{\mathrm{D}}}\) are shown in Fig.

2b. Detailed electronic and phononic properties can be seen in Supplementary Note 4. We find that above the phonon dispersion region (~8 THz, corresponds to 0.033 eV) and below the direct

bandgap (1.6 eV), the \(\chi _{11}^{{\mathrm{el}}} = 25.2 > \chi _{22}^{{\mathrm{el}}} = 7.8\). This can be understood from anisotropic electronic transition strength, reflected by the

imaginary part of susceptibility at the direct bandgap, which determines the absorbance. The absorbance is calculated by \(A_{ii}\left( \omega \right) = 1 - \exp \left( {{{ - \omega \chi

}}_{ii}^{\prime\prime} d/c} \right)\), where _c_ is the speed of light and _χ_″ is imaginary part of susceptibility. The absorbance for the _x_-polarized and _y_-polarized light at 1.6 eV

are 0.63% and 6.7%, respectively, owing to the saddle-like exciton feature and large anisotropic joint density of states. This is consistent with previous works26, and quantitative

differences come from different formulae. In the phonon frequency region, the TO phonon modes interact with LPTL strongly. The optical branches of vibrational modes at the Γ-point can be

decomposed according to the irreducible representations as $${{\Gamma }}_{{\mathrm{op}}.} = 3A_1 \oplus 2A_2 \oplus 3B_1 \oplus 1B_2.$$ (6) We find one absorbance peak [_A_11(_ω_ = 5.38 THz)

 = 0.07%] for the _x_-LPTL, and two absorbance peaks [_A_22(_ω_ = 1.2 THz) = 0.26% and _A_22(_ω_ = 1.87 THz) = 0.11%] for the _y_-LPTL. According to Kramers-Kronig relation, these modes

could produce a macroscopic polarization and contribute to nonzero (real part) susceptibility components. Since the _y_-absorbance peaks are stronger and lower in frequency than that of the

_x_-absorbance peak, the corresponding real part of \(\chi _{22}^{{\mathrm{ph}}}\) is also much larger than \(\chi _{11}^{{\mathrm{ph}}}\) at the low frequency region. For example, at the

frequency of 1 THz, the real part of \(\chi _{22}^{{\mathrm{ph}}} = 713.4\) and \(\chi _{11}^{{\mathrm{ph}}} = 3.3\). According to Eqs. (1) and (2), when an _x_-LPTL (_y_-LPTL) is applied,

the FE0 phase would have largest (lowest) GFE density among the three symmetry equivalent phases. For example, on an FE0 phase, one shines an LPTL (_ω_ = 1 THz) along 60° (or 120°) with

respect to the _x_-axis, then the FE0 could transit to FE1 (or FE2) phase, aligning the armchair direction parallel to the light polarization. We plot energy difference versus polarization

angle in Supplementary Note 5. For example, if the alternative electric field strength _E_max is chosen to be 0.2 V nm−1 (corresponding to an intermediate LPTL intensity of 5 × 109 W cm−2,

which is easily accessible experimentally) and frequency is 1 THz, according to Eqs. (1) and (2), the GFE difference between the FE1/FE2 and FE0 under _x_- or _y_-LPTL illumination will be

10.66 meV per f.u. (=1.6 μJ cm−2) (Supplementary Note 2). Thus, using LPTL, one could drive a ferroic order transition schematically plotted in Fig. 1. This energy difference is proportional

to the _E_2 (or the intensity of LPTL). Increasing light intensity could boost transition kinetics. ANISOTROPIC EELS The different ferroic order geometries can be detected via

ultrahigh-resolution EELS in the (scanning) transmission electron microscope, especially for the vibrational signatures32,33,34. Here, we will show that this approach can be employed to

distinguish different phases, owing to their anisotropic optical response. A simple approximation to describe the experimental setup is sketched inset of Fig. 3: when an electron travels

along a rectilinear trajectory (velocity _v_) parallel to the material surface (with a finite distance _b_). The local dielectric tensor (_ε_ = 1 + _χ_) is the key component entering the

classical description of electron scattering by a surface. According to a classical and quasistatic approach (refs 35,36,37,38 and Supplementary Note 6), energy loss probability per unit

angular frequency and per unit electron path length is $$\frac{{d^2P\left( {\omega ,b} \right)}}{{d\omega dr_{||}}} = \frac{{e^2}}{{\left( {2\pi } \right)^2v^2\varepsilon _0\hbar }} \times

{\int}_{ - \infty }^{ + \infty } {{\mathrm{Im}}\left[ {\frac{{\alpha \left( {\omega ,q_{||},q_ \bot } \right)\left( {e^{2Qh} - 1} \right)}}{{e^{2Qh} - \alpha ^2\left( {\omega ,q_{||},q_ \bot

} \right)}}} \right]\frac{{e^{ - 2Qb}}}{Q}dq_ \bot }.$$ (7) Here, _q_|| and _q_⊥ are wavevectors parallel and perpendicular to the electron movement direction (_q_|| = _ω_/_v_),

respectively, \(Q^2 = q_{||}^2 + q_ \bot ^2\), and $$\alpha \left( {\omega ,q_{||},q_ \bot } \right) = \frac{{\sqrt {\varepsilon _{33}\left( {\varepsilon _{||}q_{||}^2 + \varepsilon _ \bot

q_ \bot ^2} \right)/Q} - 1}}{{\sqrt {\varepsilon _{33}\left( {\varepsilon _{||}q_{||}^2 + \varepsilon _ \bot q_ \bot ^2} \right)/Q} + 1}}$$ (8) is the angle-dependent polarizability. When

the system is isotropic, the polarizability reduces to its well-known form \(\alpha = \frac{{\varepsilon - 1}}{{\varepsilon + 1}}\). We plot the EELS spectra when the electron is moving

along the _x_- and _y_-directions (Fig. 3). One could clearly observe large direction-dependent EELS feature, especially in the phonon region. This anisotropic vibrational EELS originates

from different infrared-active phonon characters of the _x_- and _y_-LPTL in the _β_-GeSe monolayer. This result provides a high-resolution damage-free approach to distinguish the geometric

structure and anisotropy when the ferroic order switches. MONOLAYER _Α_-SNTE The monolayer _β_-GeSe has three orientation variants with 120° rotation to each other, owing to the \(\hat C_3\)

character of the parental geometry. Thus, even though the LPTL response is largely anisotropic (~66 times difference in the _x_- and _y_-directions at incident frequency of 1 THz), the

energies separating different orientation variants are on the order of 1 μJ cm−2. Now we consider another 2D time-reversal invariant multiferroic material, monolayer _α_-SnTe, whose parental

geometry is _C_4 symmetric39,40. As shown in Fig. 4a, the _α_-SnTe also shows a puckered structure, with slight expansion (compression) along the _y_- (_x_-)direction. Note that even though

bulk and multilayered _α_-phase other group IV–VI compounds (such as _α_-GeS, _α_-GeSe, _α_-SnS, and _α_-SnSe) have been experimentally seen, their monolayer remains to be fabricated.

Hence, we only focus on the _α_-SnTe monolayer, and similar results for those analogs can be obtained. Our relaxation reveals that the structure also belongs to be the _Pmn_21 space group,

and the deformation strain tensor of this ferroelastic structure (FE1) is \(\eta _1 = \left( {\begin{array}{*{20}{c}} { - 0.011} & 0 \\ 0 & {0.011} \end{array}} \right)\), and the

other ferroelastic structure (FE2) is subject to 90°-rotation with \(\eta _2 = \left( {\begin{array}{*{20}{c}} {0.011} & 0 \\ 0 & { - 0.011} \end{array}} \right)\). According to the

modern theory of polarization, we find that its spontaneous polarization is 0.13 nC m−1, along the puckered direction. These results agree well with previous works39,40. We employ DFT and

DFPT methods to calculate the electron and phonon dispersions (Supplementary Note 4), and compute the optical susceptibility (Fig. 4b). We find the electronic contributed optical response

anisotropy in _α_-SnTe monolayer is not as large as that in the _β_-GeSe. At the direct optical bandgap (0.9 eV), the absorbances of the _x_-polarized and _y_-polarized light are 0.9% and

1.6%, respectively. The electronic contributed real part of susceptibility is \(\chi _{11}^{{\mathrm{el}}} = 44.7 \, > \, \chi _{22}^{{\mathrm{el}}} = 38.3\) below the direct optical

bandgap. The anisotropic excitonic absorption at different wavelength have been reported for other multiferroic 2D group-IV monochalcogenide monolayers41, which should yield polarization

direction-dependent imaginary part of optical susceptibility. Then the anisotropic electron contributed real part of susceptibility at low frequency can be obtained, according to

Kramers-Kronig relation. As for the phononic contribution, there is one clear absorbance peak for the _x_-LPTL (at 1.2 THz, _A_11 = 0.3%), and two peaks for the _y_-LPTL (at 1.4 and 1.9 THz,

with _A_22 = 0.03% and 0.2%, respectively). Hence, the phononic contributed real part of susceptibility also possesses a large anisotropy. In total, when the incident energy is 0.004 eV

(1.1 THz), the \(\chi _{11} = 731.5 \, > \, \chi _{22} = 137.4\). For example, if a _y_-LPTL (with 1.1 THz frequency and 5 × 109 W cm−2 intensity), the FE1 orientation variant (as shown

in Fig. 4a) would have higher GFE density than the FE2 orientation variant by 12.3 meV per f.u. (3.7 μJ cm−2). When one increases the LPTL intensity to 2.2 × 1010 W cm−2 (0.42 V nm−1), this

ferroic order switch (FE1 to FE2) can be _barrier-free_ (Supplementary Note 7). A similar _x_-LPTL can drive the FE2 to FE1 switch, suggesting its good reversibility. Such barrierless phase

transition does not require latent heat and can occur anywhere the LPTL is shined. This indicates a spinodal-decomposition in the reaction coordinate space, avoiding the conventional

nucleation-and-growth kinetics. Such process only requires one or a few vibrational oscillatory processes, which is ultrafast and could occur on the order of several picoseconds42. Note that

this system is interesting as its spontaneous polarization favors to align perpendicular to the optical polarization direction. This is even correct when the LPTL frequency reduces to

zero—a static electric field. Since \(\chi _{11}\left( {\omega = 0} \right) \, > \, \chi _{22}\left( {\omega = 0} \right)\), under _x_-directional electric field (magnitude larger than

0.3 V nm−1), the spontaneous polarization prefers to align along _y_, counterintuitive with the conventional E//P picture. The direction-dependent EELS spectra are also evaluated (Fig. 4c).

We again observe that the vibrational EELS shows a large anisotropy. The main EELS spectra peak positions differ by about 1 THz when the electron is moving along the _x_- (perpendicular to

spontaneous polarization Ps) and _y_- (parallel to Ps) directions. This signal shows larger directional contrast than that in the _β_-GeSe monolayer, boosting an ultrahigh-resolution

characterization of ferroic orders of _α_-SnTe monolayer. DIRECTION-DEPENDENT SHG RESPONSE We calculate the electronic contribution of SHG susceptibility \(\chi _{abc}^{\left( 2

\right)}\left( { - 2\omega ;\omega ,\omega } \right)\) of FE0-_β_-GeSe and FE1-_α_-SnTe, where _a_, _b_, _c_ are Cartesian coordinates, to measure the ferroicity via nonlinear optical

response. Under the electric field with angular frequency of _ω_, the second-order nonlinear polarization takes the form \(P_a\left( {2\omega } \right) = \chi _{abc}\left( { - 2\omega

;\omega ,\omega } \right){\cal{E}}_b\left( \omega \right){\cal{E}}_c\left( \omega \right)\), which is correlated with SHG emitted field. The SHG susceptibility can be expressed as43,44,45

$$\chi _{abc}^{\left( 2 \right)}\left( { - 2\omega ;\omega ,\omega } \right) = \chi _{abc}^{{\mathrm{II}}}\left( { - 2\omega ;\omega ,\omega } \right) + \eta _{abc}^{{\mathrm{II}}}\left( { -

2\omega ;\omega ,\omega } \right) + \sigma _{abc}^{{\mathrm{II}}}\left( { - 2\omega ;\omega ,\omega } \right),$$ (9) where \(\chi _{abc}^{{\mathrm{II}}}\left( { - 2\omega ;\omega ,\omega }

\right)\), \(\eta _{abc}^{{\mathrm{II}}}\left( { - 2\omega ;\omega ,\omega } \right)\), and \(\sigma _{abc}^{{\mathrm{II}}}\left( { - 2\omega ;\omega ,\omega } \right)\) are interband

transition, modulation of the linear susceptibility due to intraband motions of electrons, and modification by the polarization energy associated with interband motions, respectively.

Specially, they can be evaluated by $$\chi _{abc}^{{\mathrm{II}}}\left( { - 2\omega ;\omega ,\omega } \right) = \frac{{e^3}}{{{\hbar}^2}}\mathop {\sum}\limits_{nml} {{\int}

{\frac{{d^3\mathbf{k}}}{{\left( {2\pi } \right)^3}}\frac{{r_{nm}^ar_{ml}^br_{ln}^c}}{{\omega _{ln} - \omega _{ml}}} \times \left( {\frac{{f_{ml}}}{{\omega _{ml} - \omega }} +

\frac{{f_{ln}}}{{\omega _{ln} - \omega }} + \frac{{2f_{nm}}}{{\omega _{mn} - 2\omega }}} \right)} },$$ (10) $$\eta _{abc}^{{\mathrm{II}}}\left( { - 2\omega ;\omega ,\omega } \right) =

\frac{{e^3}}{{{\hbar}^2}}{\int} {\frac{{d^3\mathbf{k}}}{{\left( {2\pi } \right)^3}}} \left\{ {\mathop {\sum}\limits_{nml} {\omega _{mn}r_{nm}^a\left\{ {r_{ml}^br_{ln}^c} \right\}} \left[

{\frac{{f_{nl}}}{{\omega _{ln}^2\left( {\omega _{ln}^2 - \omega } \right)}} + \frac{{f_{lm}}}{{\omega _{ml}^2\left( {\omega _{ml}^2 - \omega } \right)}}} \right] - 8i\mathop

{\sum}\limits_{nm} {\frac{{f_{nm}r_{nm}^a}}{{\omega _{mn}^2\left( {\omega _{mn} - 2\omega } \right)}}} \left\{ {{{\Delta }}_{mn}^br_{mn}^c} \right\} - 2\frac{{\mathop {\sum }\limits_{nml}

f_{nm}r_{nm}^a\left\{ {r_{ml}^br_{ln}^c} \right\}\left( {\omega _{ln} - \omega _{ml}} \right)}}{{\omega _{mn}^2\left( {\omega _{mn} - 2\omega } \right)}}} \right\},$$ (11) $$\sigma

_{abc}^{{\mathrm{II}}}\left( { - 2\omega ;\omega ,\omega } \right) = \frac{{e^3}}{{2{\hbar}^2}}{\int} {\frac{{d^3\mathbf{k}}}{{\left( {2\pi } \right)^3}}} \left\{ {\mathop

{\sum}\limits_{nml} {\frac{{f_{nm}}}{{\omega _{mn}^2\left( {\omega _{mn} - \omega } \right)}}} \left[ {\omega _{nl}r_{lm}^a\left\{ {r_{mn}^br_{nl}^c} \right\} - \omega _{lm}r_{nl}^a\left\{

{r_{lm}^br_{mn}^c} \right\}} \right] + i\mathop {\sum}\limits_{nm} {\frac{{f_{nm}r_{nm}^a\left\{ {{{\Delta }}_{mn}^br_{mn}^c} \right\}}}{{\omega _{mn}^2\left( {\omega _{mn} - \omega }

\right)}}} } \right\},$$ (12) where the position matrix element is defined as \({r_{nm}^a}\left(k \right) = \frac{\langle{n{\mathbf{k}}}{|} {\hat v^a}{|}{m}{\mathbf{k}}\rangle}{{i\omega

_{nm}}}\left({n \ne m} \right)\), the energy and Fermi-Dirac distribution difference between state-_n_ and state-_m_ at K being \(\omega _{nm}\left( \mathbf{k} \right) = \omega _n\left(

\mathbf{k} \right) - \omega _m\left( \mathbf{k} \right)\), \(f_{nm}\left( \mathbf{k} \right) = f_n\left( \mathbf{k} \right) - f_m\left( \mathbf{k} \right)\). In addition, \(\left\{

{r_{ml}^b}{r_{ln}^c} \right\} = \frac{1}{2}\left( {{r_{ml}^b}r_{ln}^c + r_{ml}^cr_{ln}^b} \right)\) and \({{\Delta }}_{mn}^b = v_{mm}^b - v_{nn}^b\) are velocity differences. Here the clear

dependence on wavevector K is omitted. We fit the DFT calculated electronic structure with maximally localized Wannier functions46 and use a dense K-mesh sampling of (81 × 81 × 1) to perform

the integration. Once the bands around the Fermi levels are sufficiently fitted, this localized basis set should yield very similar results as directly using DFT wavefunctions47. Note that

when more bands are included, the better converged results could be obtained. In our calculations, we fit 32 bands by including all s and p orbitals of Ge, Se, Sn, and Te atoms, which could

reproduce the bands between –5 and 5 eV around the Fermi level. We then assume that those band sets can serve as a good approximate to a full complete wavefunction basis set and perform

calculations. We calculate the \(\chi _{yxx}\left( { - 2\omega ;\omega ,\omega } \right)\) which reflects the _y_-polarized response under _x_-polarized light, as shown in Fig. 5. The

previously mentioned re-scaling scheme is also applied. The real and imaginary parts of \(\chi _{yxx}\left( { - 2\omega ;\omega ,\omega } \right)\) are shown as red and green curves (Fig.

5a, d), which satisfies Kramers-Kronig relationship. The first peak of GeSe lies at ℏ_ω_ = 0.82 eV, with \(|\chi _{yxx}\left( { - 2\omega ;\omega ,\omega } \right)|\) is 0.21 (=0.11 − 0.17 ×

 _i_) nm V−1 (Fig. 5b). The momentum distribution of SHG is shown in the inset. One clearly observes that it is mainly contributed around the R (=1/2, 1/2, 0) point, consistent with its

direct electronic bandgap at R (Supplementary Note 4). As for the SnTe, its smaller bandgap yields larger SHG susceptibility. At ℏ_ω_ = 0.36 eV, \(| {\chi _{yxx}\left( { - 2\omega ;\omega

,\omega } \right)} |\) becomes as large as 24.4 (=14.54 + 19.53 × _i_) nm V−1 (Fig. 5e). From its momentum distribution, the dominate contribution comes from ±Λ [=(0, ±0.55, 0) Å−1] point,

which are the direct bandgap position of its electronic band structure (Supplementary Note 4). One also notes that \(\chi _{xxx}\left( { - 2\omega ;\omega ,\omega } \right) = \chi

_{yxy}\left( { - 2\omega ;\omega ,\omega } \right) = \chi _{yyx}\left( { - 2\omega ;\omega ,\omega } \right) = 0\) due to mirror symmetry \({\cal{M}}_x\). We calculate other components

(Supplementary Note 8) and plot polar distribution of anisotropic SHG susceptibility \(| {\chi _{||}\left( \theta \right)}|\) (red curve) and \(\left| {\chi _ \bot \left( \theta \right)}

\right|\) (green curve) in Fig. 5c, f, where _θ_ is angle that between the light polarization and _x_-axis. This also corresponds to ferroicity-dependent SHG susceptibility. For example, for

90°-rotated FE2-_α_-SnTe, the \(\chi _{xxx}\left( { - 2\omega ;\omega ,\omega } \right) \ne \, 0\) and \(\chi _{yxx}\left( { - 2\omega ;\omega ,\omega } \right) = 0\). We briefly compare

our SHG calculations with previous works on other centrosymmetry broken 2D time-reversal multiferroic materials. The point groups of _β_-GeSe and _α_-SnTe are both _C_2_v_, which suggests

that their SHG anisotropy is the same as those of other _α_-phase group IV–VI monolayers, such as monolayers GeS and SnSe. It is well-known that the standard DFT calculations underestimate

the electronic bandgap, hence in principles one has to adopt more accurate approaches, such as hybrid functional or many-body calculations48. Unfortunately, this is extremely computational

challenging because of nonlocal interactions included and very dense K-mesh needed. In the previous calculations47, Wang and Qian adopted a scissor operator scheme to shift the DFT

calculated bandgap to be consistent with the optical bandgap calculated through many-body theory with exciton interaction correction. In our current study, we only focus on the relative

strength of SHG responses and its anisotropy, and the scissor correction is not applied here. Therefore, in experiments, the peaks of the SHG responses would shift to higher energy compared

with our calculation results, usually by ~0.5–1 eV. According to Eqs. (10)–(12), the SHG peak magnitude roughly scales as _E_g−1 (_E_g is the electronic bandgap). Hence, the experimentally

measured SHG peak values would also be smaller than our calculation results. Nevertheless, the overall SHG shape and anisotropy trend should be similar with our calculations. DISCUSSION In

our current LPTL-driven phase transition mechanism, we only elaborate the real part of optical susceptibility and assume its imaginary part to be small. Actually when the optical frequency

is chosen around finite Im(_χ__ii_), direct optical absorption occurs. This corresponds to a conversion from photons to infra-red phonon or electron-hole pairs (electron interband

excitation). Actually this could also trigger mechanical strains49 to the systems, and phase transition may occur50,51. However, they would generate unwanted heat into the system via

phonon–phonon interaction or non-radiative electron-hole recombination so that the device reversibility will reduce. In our mechanism, only electric field work done is applied to the

systems. Even though all of the work done converts to be internal energy, the temperature rise is still very small30. Hence, in the current discussion, we only focus on the real part of

optical susceptibility contribution and avoid such direct photon absorption process. In addition, we note that in the current setup only ferroelastic order degeneracy can be broken under

LPTL irradiation. One could not break the degeneracy between the ferroelectric phases with opposite polarization states (Ps and −Ps). Thus, one usually needs to apply additional fields (such

as zero frequency static electric field or introducing surface/interface effects) to break such degeneracy. Actually, optical control of the polarization direction (reversal by 180°) has

been observed experimentally in BaTiO3 thin films deposited on transition metal dichalcogenide monolayers52, where surface effects are adopted. This suggests that when additional spatial

broken interactions exist, one could further break the polarization degeneracy and control/manipulate the ferroelectric order. Very recently, Chen et al. have shown that photo-induced

polarization flip could be realized in 3D ferroelectric GeTe and PbTiO3, which is controlled by the synergistic effects of excited-state energy profile, atomic motion velocity, and the

dephasing of excited states53. In the current model, we use a single unit cell to perform calculations. This indicates a coherent phase transition, with all unit cells rotate their ferroic

orders simultaneously. In reality, the materials could contain domain walls, which spatially separate different ferroic orientation variants. The phase transition usually occurs along with

the domain wall movement2. Compared with 3D materials or thick films (such as BaTiO3 perovskite) in which the domain wall is a 2D interface, in 2D materials, the domain wall is actually in

quasi-1D. This would make the phase transition in 2D materials much easier than the 3D materials or thick films, with smaller residue stress during transition. Hence, the phase transition

would cost low energy and can be nonvolatile. Previous works have shown that the domain wall migration energy is around a few to a few tens meV per Å28,41, much smaller than the formation

energy of 1D defects in 3D materials, indicating a fast domain-wall-assisted order switching. Therefore, the existence of domain wall may facilitate the phase transition. In addition, light

usually has a large spot size, on the order of μm or larger. Once the light polarization, intensity, and frequency are carefully selected, it can trigger ultrafast and barrier-free phase

transition wherever it irradiates. The conventional nucleation and growth kinetics can be avoided, and a coherent phase transition can be expected. However, if the material contains a large

number of point, line or area defects, they may strongly affect the optical response functions, by significantly reducing the carrier lifetime (_τ_) and introducing doping levels into the

phonon and electron band structures. For example, defects usually bring both shallow and deep energy levels into semiconductor electronic bands54. Then the low frequency light may be

absorbed and the imaginary part of optical susceptibility needs to be considered. When the defect or impurity concentration is not high enough, their effects can be phenomenologically

incorporated by the finite lifetime in the response formulae, in Eqs. (4) and (5). In order to avoid such uncertainty, we choose the light frequency away from resonant absorption

frequencies, then the exact value of lifetime does not affect the estimate too much. However, an exact and complete evaluation is very complicated, and is out of the scope of our current

study. The phase transition related strain in the current study is within 8%, which is usually sustainable for 2D materials. However, one may notice that if the sample size is a few to a few

tens of nanometers, then such strain would cause a big lattice mismatch in the system, and may even affect the chemical bonds at the boundary between the transformed and untransformed

domains. A direct numerical simulation of such process in a large area is very computationally challenging and memory demanding. Actually one may allow a freestanding 2D material to be

slightly slack in the _z_-direction, or carefully select a surface to support them with weak interactions (such as van der Waals)55. This allows the atoms to move in the _z_-direction with

small energy cost, which could effectively release the in-plane strains during martensitic phase transitions. This is different from 3D bulk materials or thick films, where such

strain-induced damage can be significantly large and is fatal to their reversible usage. We implement a theoretically and computationally combined approach to study the electronic, phononic,

and mechanical responses of 2D time-reversal invariant multiferroic monolayers under LPTL illumination. In a microscopic picture, light irradiation could trigger a strong and coherent Raman

vibration in the system (Supplementary Note 9). Taking two experimentally fabricated multiferroic _β_-GeSe and _α_-SnTe monolayers as examples, we find that they both show a large

anisotropic optical response, especially at the terahertz region, owing to their selective vibrational infrared-active vibrational modes. According to the thermodynamic theory, we predict

that LPTL can efficiently drive phase transition with an ultrafast kinetics (or even a barrierless GFE profile). This noncontacting optomechanical approach to switch the ferroelastic order

of 2D materials can be easily controlled experimentally. In order to detect different orders, we propose to measure vibrational EELS spectra, which is direction dependent and has an

ultrahigh resolution experimentally. Anisotropic SHG response calculations are also provided. Different from the parameterized phonon–phonon coupling models in the optically induced phase

transition, we provide a first-principles quantitative estimation on the terahertz optics effects. Such mechanism can also apply to other frequency range, as long as the direct optical

absorption is eliminated. Owing to the rapid developments of using terahertz laser triggering (topological or structural) phase transition, our theory provides a route to explain these

experiments and predict unexplored phenomena from a precise first-principles approach. METHODS DENSITY FUNCTIONAL THEORY The first-principles calculations are based on DFT14 and performed in

the Vienna Ab initio Simulation Package (VASP)15 with generalized gradient approximation (GGA) treatment of exchange and correlation functional in the solid-state PBE form (PBEsol)16. In

order to simulate 2D materials, a vacuum distance of 15 Å in the out-of-plane _z_-direction is applied to eliminate layer interactions. The projector augmented wave method17 and planewave

basis set are applied to treat the core and valence electrons, respectively. The kinetic cutoff energy of planewave is set to be 350 eV. The reciprocal space is represented by Monkhorst-Pack

K-mesh scheme18. The electron and force component convergence criteria are set to be 10–7 eV and 0.001 eV Å−1, respectively. Spin-orbit coupling (SOC) interactions are included

self-consistently. For the SnTe monolayer, _d_-electrons are incorporated in the Sn valence electrons, and Grimme’s D3 scheme19 is used to include van der Waals interactions, which has been

demonstrated to yield good results on the energy barrier39. We use modern theory of polarization based on Berry phase approach56,57 to evaluate spontaneous polarization (Supplementary Note

3). The optical dielectric function contributed from the electron subsystems are calculated according to RPA. In order to evaluate dielectric function contributed from the ionic subsystem

(phonons), we use DFPT20 to calculate vibrational modes and Born effective charge of each ion, with the help of the Phonopy code21. DATA AVAILABILITY All data generated or analyzed during

this study are included in this published article (and its Supplementary Information files), and are available from the authors upon reasonable request. CODE AVAILABILITY The related codes

are available from the authors upon reasonable request. REFERENCES * Akamatsu, H. et al. Light-activated gigahertz ferroelectric domain dynamics. _Phys. Rev. Lett._ 120, 096101 (2018).

Article  Google Scholar  * Rubio-Marcos, F. et al. Reversible optical control of macroscopic polarization in ferroelectrics. _Nat. Photonics_ 12, 29–32 (2018). Article  CAS  Google Scholar 

* Liou, Y.-D. et al. Deterministic optical control of room temperature multiferroicity in BiFeO3 thin films. _Nat. Mater._ 18, 580–587 (2019). Article  CAS  Google Scholar  * Nova, T. F.,

Disa, A. S., Fechner, M. & Cavalleri, A. Metastable ferroelectricity in optically strained SrTiO3. _Science_ 364, 1075–1079 (2019). Article  CAS  Google Scholar  * Li, X. et al.

Terahertz field–induced ferroelectricity in quantum paraelectric SrTiO3. _Science_ 364, 1079–1082 (2019). Article  CAS  Google Scholar  * Yeats, A. L. et al. Persistent optical gating of a

topological insulator. _Sci. Adv._ 1, e1500640 (2015). Article  CAS  Google Scholar  * Salén, P. et al. Matter manipulation with extreme terahertz light: progress in the enabling THz

technology. _Phys. Rep._ 836–837, 1–74 (2019). Article  Google Scholar  * Mankowsky, R., von Hoegen, A., Först, M. & Cavalleri, A. Ultrafast reversal of the ferroelectric polarization.

_Phys. Rev. Lett._ 118, 197601 (2017). Article  CAS  Google Scholar  * Qi, T., Shin, Y.-H., Yeh, K.-L., Nelson, K. A. & Rappe, A. M. Collective coherent control: synchronization of

polarization in ferroelectric PbTiO3 by shaped THz fields. _Phys. Rev. Lett._ 102, 247603 (2009). Article  CAS  Google Scholar  * von Rohr, F. O. et al. High-pressure synthesis and

characterization of β-GeSe—a six-membered-ring semiconductor in an uncommon boat conformation. _J. Am. Chem. Soc._ 139, 2771–2777 (2017). Article  CAS  Google Scholar  * Chang, K. et al.

Discovery of robust in-plane ferroelectricity in atomic-thick SnTe. _Science_ 353, 274–278 (2016). Article  CAS  Google Scholar  * Stengel, M., Spaldin, N. A. & Vanderbilt, D. Electric

displacement as the fundamental variable in electronic-structure calculations. _Nat. Phys._ 5, 304–308 (2009). Article  CAS  Google Scholar  * Song, W., Fei, R. & Yang, L. Off-plane

polarization ordering in metal chalcogen diphosphates from bulk to monolayer. _Phys. Rev. B_ 96, 235420 (2017). Article  Google Scholar  * Hohenberg, P. & Kohn, W. Inhomogeneous electron

gas. _Phys. Rev._ 136, B864 (1964). Article  Google Scholar  * Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis

set. _Phys. Rev. B_ 54, 11169 (1996). Article  CAS  Google Scholar  * Perdew, J. P. et al. Restoring the density-gradient expansion for exchange in solids and surfaces. _Phys. Rev. Lett._

100, 136406 (2008). Article  CAS  Google Scholar  * Blöchl, P. E. Projector augmented-wave method. _Phys. Rev. B_ 50, 17953 (1994). Article  Google Scholar  * Monkhorst, H. J. & Pack, J.

D. Special points for Brillouin-zone integrations. _Phys. Rev. B_ 13, 5188 (1976). Article  Google Scholar  * Grimme, S., Antony, J., Ehrlich, S. & Krieg, S. A consistent and accurate

ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _J. Chem. Phys._ 132, 154104 (2010). Article  CAS  Google Scholar  * Baroni, S. &

Resta, R. Ab initio calculation of the macroscopic dielectric constant in silicon. _Phys. Rev. B_ 33, 7017 (1986). Article  CAS  Google Scholar  * Togo, A., Oba, F. & Tanaka, I.

First-principles calculations of the ferroelastic transition between rutile-type and CaCl2-type SiO2 at high pressures. _Phys. Rev. B_ 78, 134106 (2008). Article  CAS  Google Scholar  *

Born, M. & Huang, K. _Dynamical theory of crystal lattices_ (Oxford University Press, 1954). * Gonze, X. & Lee, C. Dynamical matrices, Born effective charges, dielectric permittivity

tensors, and interatomic force constants from density-functional perturbation theory. _Phys. Rev. B_ 55, 10355 (1997). Article  CAS  Google Scholar  * Leite Alves, H. W., Neto, A. R. R.,

Scolfaro, L. M. R., Myers, T. H. & Borges, P. D. Lattice contribution to the high dielectric constant of PbTe. _Phys. Rev. B_ 87, 115204 (2013). Article  CAS  Google Scholar  * Winta, C.

J., Wolf, M. & Paarmann, A. Low-temperature infrared dielectric function of hyperbolic α-quartz. _Phys. Rev. B_ 99, 144308 (2019). Article  CAS  Google Scholar  * Luo, N. et al.

Saddle‐point excitons and their extraordinary light absorption in 2D β‐phase group‐IV monochalcogenides. _Adv. Funct. Mater._ 28, 1804581 (2018). Article  CAS  Google Scholar  * Guan, S.,

Liu, C., Lu, Y., Yao, Y. & Yang, S. A. Tunable ferroelectricity and anisotropic electric transport in monolayer β-GeSe. _Phys. Rev. B_ 97, 144104 (2018). Article  Google Scholar  * Li,

W. & Li, J. Ferroelasticity and domain physics in two-dimensional transition metal dichalcogenide monolayers. _Nat. Commun._ 7, 10843 (2016). Article  CAS  Google Scholar  * Jiang, Z.,

Liu, Z., Li, Y. & Duan, W. Scaling universality between band gap and exciton binding energy of two-dimensional semiconductors. _Phys. Rev. Lett._ 118, 266401 (2017). Article  Google

Scholar  * Zhou, J., Xu, H., Li, Y., Jaramillo, R. & Li, J. Opto-mechanics driven fast martensitic transition in two-dimensional materials. _Nano Lett._ 18, 7794–7800 (2018). Article 

CAS  Google Scholar  * Laturia, A., Van de Put, M. L. & Vandenberghe, W. G. Dielectric properties of hexagonal boron nitride and transition metal dichalcogenides: from monolayer to bulk.

_npj 2D Mater. Appl._ 2, 6 (2018). Article  CAS  Google Scholar  * Krivanek, O. L. et al. Vibrational spectroscopy in the electron microscope. _Nature_ 514, 209–212 (2014). Article  CAS 

Google Scholar  * Lagos, M. J., Trügler, A., Hohenester, U. & Batson, P. E. Mapping vibrational surface and bulk modes in a single nanocube. _Nature_ 543, 529–532 (2017). Article  CAS 

Google Scholar  * Ibach, H. & Mills, D. L. _Electron energy loss spectroscopy and surface vibrations_ (Academic Press, London, 1982). * Radtke, G., Taverna, D., Lazzeri, M. & Balan,

E. First-principles vibrational electron energy loss spectroscopy of _β_-guanine. _Phys. Rev. Lett._ 119, 027402 (2017). Article  CAS  Google Scholar  * Radtke, G. et al. Polarization

selectivity in vibrational electron-energy-loss spectroscopy. _Phys. Rev. Lett._ 123, 256001 (2019). Article  CAS  Google Scholar  * Kociak, M. et al. Experimental evidence of

surface-plasmon coupling in anisotropic hollow nanoparticles. _Phys. Rev. Lett._ 87, 075501 (2001). Article  CAS  Google Scholar  * Chen, C. H. & Silcox, J. Calculations of the

electron-energy-loss probability in thin uniaxial crystals at oblique incidence. _Phys. Rev. B_ 20, 3605 (1979). Article  CAS  Google Scholar  * Ye, Q.-J., Liu, Z.-Y., Feng, Y., Gao, P.

& Li, X.-Z. Ferroelectric problem beyond the conventional scaling law. _Phys. Rev. Lett._ 121, 135702 (2018). Article  CAS  Google Scholar  * Liu, K., Lu, J., Picozzi, S., Bellaiche, L.

& Xiang, H. Intrinsic origin of enhancement of ferroelectricity in SnTe ultrathin films. _Phys. Rev. Lett._ 121, 027601 (2018). Article  Google Scholar  * Wang, H. & Qian, X.

Two-dimensional multiferroics in monolayer group IV monochalcogenides. _2D Mater._ 4, 015042 (2017). Article  Google Scholar  * Xu, H., Zhou, J., Li, Y., Jaramillo, R. & Li, J.

Optomechanical control of stacking patterns of h-BN bilayer. _Nano Res._ 12, 2634–2639 (2019). Article  CAS  Google Scholar  * Aversa, C. & Sipe, J. Nonlinear optical susceptibilities of

semiconductors: results with a length-gauge analysis. _Phys. Rev. B_ 52, 14636–14645 (1995). Article  CAS  Google Scholar  * Sipe, J. E. & Ghahramani, E. Nonlinear optical response of

semiconductors in the independent-particle approximation. _Phys. Rev. B_ 48, 11705–11722 (1993). Article  CAS  Google Scholar  * Hughes, J. L. P. & Sipe, J. E. Calculation of

second-order optical response in semiconductors. _Phys. Rev. B_ 53, 10751–10763 (1996). Article  CAS  Google Scholar  * Mostofi, A. A. et al. wannier90: A tool for obtaining

maximally-localised Wannier functions. _Comput. Phys. Commun._ 178, 685–699 (2008). Article  CAS  Google Scholar  * Wang, H. & Qian, X. Giant Optical second harmonic generation in

two-dimensional multiferroics. _Nano Lett._ 17, 5027–5034 (2017). Article  CAS  Google Scholar  * Fei, R., Tan, L. Z. & Rappe, A. M. Shift-current bulk photovoltaic effect influenced by

quasiparticle and exciton. _Phys. Rev. B_ 101, 045104 (2020). Article  CAS  Google Scholar  * Hu, H. et al. Quantum electronic stress: density-functional-theory formulation and physical

manifestation. _Phys. Rev. Lett._ 109, 055501 (2012). Article  CAS  Google Scholar  * Peng, B. et al. Sub-picosecond photo-induced displacive phase transition in two-dimensional MoTe2. _npj

2D Mater. Appl._ 4, 14 (2020). Article  CAS  Google Scholar  * Si, C. et al. Photoinduced vacancy ordering and phase transition in MoTe2. _Nano Lett._ 19, 3612–3617 (2019). Article  CAS 

Google Scholar  * Li, T. et al. Optical control of polarization in ferroelectric heterostructures. _Nat. Commun._ 9, 3344 (2018). Article  CAS  Google Scholar  * Chen, N.-K. et al. Optical

subpicosecond nonvolatile switching and electron-phonon coupling in ferroelectric materials. _Phys. Rev. B_ 102, 184115 (2020). Article  CAS  Google Scholar  * Mieghem, P. V. Theory of band

tails in heavily doped semiconductors. _Rev. Mod. Phys._ 64, 755–793 (1992). Article  Google Scholar  * Zhou, J., Mao, S. & Zhang, S. Noncontacting optostriction driven anisotropic and

inhomogeneous strain in two-dimensional materials. _Phys. Rev. Res._ 2, 022059(R) (2020). Article  Google Scholar  * King-Smith, R. D. & Vanderbilt, D. Theory of polarization of

crystalline solids. _Phys. Rev. B_ 47, 1651–1654 (1993). Article  CAS  Google Scholar  * Resta, R. Macroscopic polarization in crystalline dielectrics: the geometric phase approach. _Rev.

Mod. Phys._ 66, 899–915 (1994). Article  CAS  Google Scholar  Download references ACKNOWLEDGEMENTS This work is supported under the National Key Research and Development Program (Grant No.

2019YFA0210600), the National Natural Science Foundation of China (NSFC) under Grant Nos. 21903063, 11974270, and 11904350, and the Young Talent Startup Program of Xi’an Jiaotong University.

S.Z. also acknowledges the support of Anhui Provincial Natural Science Foundation (Grant. No. 2008085QA30). AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Center for Alloy Innovation and

Design, State Key Laboratory for Mechanical Behavior of Materials, Xi’an Jiaotong University, Xi’an, 710049, China Jian Zhou * International Center for Quantum Design of Functional Materials

(ICQD), Hefei National Laboratory for Physical Sciences at Microscale, and CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of

China, Hefei, 230026, Anhui, China Shunhong Zhang Authors * Jian Zhou View author publications You can also search for this author inPubMed Google Scholar * Shunhong Zhang View author

publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS J.Z. conceived the concept. J.Z. and S.Z. performed the calculations. J.Z. analyzed the data. J.Z. and

S.Z. wrote the manuscript. CORRESPONDING AUTHORS Correspondence to Jian Zhou or Shunhong Zhang. ETHICS DECLARATIONS COMPETING INTERESTS The authors declare no competing interests. ADDITIONAL

INFORMATION PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. SUPPLEMENTARY INFORMATION

SUPPLEMENTARYMATERIAL RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation,

distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and

indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to

the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will

need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE

CITE THIS ARTICLE Zhou, J., Zhang, S. Terahertz optics-driven phase transition in two-dimensional multiferroics. _npj 2D Mater Appl_ 5, 16 (2021). https://doi.org/10.1038/s41699-020-00189-7

Download citation * Received: 02 September 2020 * Accepted: 02 December 2020 * Published: 22 January 2021 * DOI: https://doi.org/10.1038/s41699-020-00189-7 SHARE THIS ARTICLE Anyone you

share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to clipboard Provided by the

Springer Nature SharedIt content-sharing initiative