## 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 induced electronic, phononic, and mechanical responses of two experimentally fabricated multiferroic (with both ferroelastic and ferroelectric) materials, *β*-GeSe and *α*-SnTe monolayer. We show that the relative stability of different orientation variants can be effectively manipulated via the polarization direction of the terahertz laser, which is selectively and strongly 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.

## Introduction

Optical control of material geometries is receiving rapidly growing attention in very recent years. For example, some ferroelectric/multiferroic perovskites, such as BaTiO_{3} (refs ^{1,2}), BiFeO_{3} (ref. ^{3}), and SrTiO_{3} (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 scales^{6}.

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 practically^{7}. 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 BaTiO_{3}, 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 interactions^{8,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 *β*-GeSe^{10} and *α*-SnTe^{11} 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 **E**_{max} 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 condition^{12,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 **P**_{s} 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 **P**_{s} and −**P**_{s}. Here *E*_{max} indicates the magnitude of electric field vector, |**E**_{max}|. When the *E*_{ma}_{x} 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)

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

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 details^{14,15,16,17,18,19,20,21}). We use the Vienna ab initio simulation package^{15} 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 package^{21}.

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

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

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 theory^{22}, 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 to^{23}

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 couplings^{24,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*2_{1} space group (no. 31) without centrosymmetry, consistent with previous works^{26,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 monolayers^{28}, 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 **P**_{s}. Consistent with previous works^{27}, 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 model^{29,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 works^{26}, 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

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 × 10^{9} 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 signatures^{32,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

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

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} symmetric^{39,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*2_{1} 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 works^{39,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 monolayers^{41}, 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 × 10^{9} 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 × 10^{10} 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 picoseconds^{42}.

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 **P**_{s}) and *y*- (parallel to **P**_{s}) 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 as^{43,44,45}

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

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 functions^{46} 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 wavefunctions^{47}. 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*_{2v}, 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 calculations^{48}. Unfortunately, this is extremely computational challenging because of nonlocal interactions included and very dense **k**-mesh needed. In the previous calculations^{47}, 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 strains^{49} to the systems, and phase transition may occur^{50,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 small^{30}. 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 (**P**_{s} and −**P**_{s}). 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 BaTiO_{3} thin films deposited on transition metal dichalcogenide monolayers^{52}, 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 PbTiO_{3}, which is controlled by the synergistic effects of excited-state energy profile, atomic motion velocity, and the dephasing of excited states^{53}.

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 movement^{2}. Compared with 3D materials or thick films (such as BaTiO_{3} 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 bands^{54}. 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 DFT^{14} 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 method^{17} 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 scheme^{18}. 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 scheme^{19} is used to include van der Waals interactions, which has been demonstrated to yield good results on the energy barrier^{39}. We use modern theory of polarization based on Berry phase approach^{56,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 DFPT^{20} to calculate vibrational modes and Born effective charge of each ion, with the help of the Phonopy code^{21}.

## 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

- 1.
Akamatsu, H. et al. Light-activated gigahertz ferroelectric domain dynamics.

*Phys. Rev. Lett.***120**, 096101 (2018). - 2.
Rubio-Marcos, F. et al. Reversible optical control of macroscopic polarization in ferroelectrics.

*Nat. Photonics***12**, 29–32 (2018). - 3.
Liou, Y.-D. et al. Deterministic optical control of room temperature multiferroicity in BiFeO

_{3}thin films.*Nat. Mater.***18**, 580–587 (2019). - 4.
Nova, T. F., Disa, A. S., Fechner, M. & Cavalleri, A. Metastable ferroelectricity in optically strained SrTiO

_{3}.*Science***364**, 1075–1079 (2019). - 5.
Li, X. et al. Terahertz field–induced ferroelectricity in quantum paraelectric SrTiO

_{3}.*Science***364**, 1079–1082 (2019). - 6.
Yeats, A. L. et al. Persistent optical gating of a topological insulator.

*Sci. Adv.***1**, e1500640 (2015). - 7.
Salén, P. et al. Matter manipulation with extreme terahertz light: progress in the enabling THz technology.

*Phys. Rep.***836–837**, 1–74 (2019). - 8.
Mankowsky, R., von Hoegen, A., Först, M. & Cavalleri, A. Ultrafast reversal of the ferroelectric polarization.

*Phys. Rev. Lett.***118**, 197601 (2017). - 9.
Qi, T., Shin, Y.-H., Yeh, K.-L., Nelson, K. A. & Rappe, A. M. Collective coherent control: synchronization of polarization in ferroelectric PbTiO

_{3}by shaped THz fields.*Phys. Rev. Lett.***102**, 247603 (2009). - 10.
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). - 11.
Chang, K. et al. Discovery of robust in-plane ferroelectricity in atomic-thick SnTe.

*Science***353**, 274–278 (2016). - 12.
Stengel, M., Spaldin, N. A. & Vanderbilt, D. Electric displacement as the fundamental variable in electronic-structure calculations.

*Nat. Phys.***5**, 304–308 (2009). - 13.
Song, W., Fei, R. & Yang, L. Off-plane polarization ordering in metal chalcogen diphosphates from bulk to monolayer.

*Phys. Rev. B***96**, 235420 (2017). - 14.
Hohenberg, P. & Kohn, W. Inhomogeneous electron gas.

*Phys. Rev.***136**, B864 (1964). - 15.
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). - 16.
Perdew, J. P. et al. Restoring the density-gradient expansion for exchange in solids and surfaces.

*Phys. Rev. Lett.***100**, 136406 (2008). - 17.
Blöchl, P. E. Projector augmented-wave method.

*Phys. Rev. B***50**, 17953 (1994). - 18.
Monkhorst, H. J. & Pack, J. D. Special points for Brillouin-zone integrations.

*Phys. Rev. B***13**, 5188 (1976). - 19.
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). - 20.
Baroni, S. & Resta, R. Ab initio calculation of the macroscopic dielectric constant in silicon.

*Phys. Rev. B***33**, 7017 (1986). - 21.
Togo, A., Oba, F. & Tanaka, I. First-principles calculations of the ferroelastic transition between rutile-type and CaCl

_{2}-type SiO_{2}at high pressures.*Phys. Rev. B***78**, 134106 (2008). - 22.
Born, M. & Huang, K.

*Dynamical theory of crystal lattices*(Oxford University Press, 1954). - 23.
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). - 24.
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). - 25.
Winta, C. J., Wolf, M. & Paarmann, A. Low-temperature infrared dielectric function of hyperbolic α-quartz.

*Phys. Rev. B***99**, 144308 (2019). - 26.
Luo, N. et al. Saddle‐point excitons and their extraordinary light absorption in 2D β‐phase group‐IV monochalcogenides.

*Adv. Funct. Mater.***28**, 1804581 (2018). - 27.
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). - 28.
Li, W. & Li, J. Ferroelasticity and domain physics in two-dimensional transition metal dichalcogenide monolayers.

*Nat. Commun.***7**, 10843 (2016). - 29.
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). - 30.
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). - 31.
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). - 32.
Krivanek, O. L. et al. Vibrational spectroscopy in the electron microscope.

*Nature***514**, 209–212 (2014). - 33.
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). - 34.
Ibach, H. & Mills, D. L.

*Electron energy loss spectroscopy and surface vibrations*(Academic Press, London, 1982). - 35.
Radtke, G., Taverna, D., Lazzeri, M. & Balan, E. First-principles vibrational electron energy loss spectroscopy of

*β*-guanine.*Phys. Rev. Lett.***119**, 027402 (2017). - 36.
Radtke, G. et al. Polarization selectivity in vibrational electron-energy-loss spectroscopy.

*Phys. Rev. Lett.***123**, 256001 (2019). - 37.
Kociak, M. et al. Experimental evidence of surface-plasmon coupling in anisotropic hollow nanoparticles.

*Phys. Rev. Lett.***87**, 075501 (2001). - 38.
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). - 39.
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). - 40.
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). - 41.
Wang, H. & Qian, X. Two-dimensional multiferroics in monolayer group IV monochalcogenides.

*2D Mater.***4**, 015042 (2017). - 42.
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). - 43.
Aversa, C. & Sipe, J. Nonlinear optical susceptibilities of semiconductors: results with a length-gauge analysis.

*Phys. Rev. B***52**, 14636–14645 (1995). - 44.
Sipe, J. E. & Ghahramani, E. Nonlinear optical response of semiconductors in the independent-particle approximation.

*Phys. Rev. B***48**, 11705–11722 (1993). - 45.
Hughes, J. L. P. & Sipe, J. E. Calculation of second-order optical response in semiconductors.

*Phys. Rev. B***53**, 10751–10763 (1996). - 46.
Mostofi, A. A. et al. wannier90: A tool for obtaining maximally-localised Wannier functions.

*Comput. Phys. Commun.***178**, 685–699 (2008). - 47.
Wang, H. & Qian, X. Giant Optical second harmonic generation in two-dimensional multiferroics.

*Nano Lett.***17**, 5027–5034 (2017). - 48.
Fei, R., Tan, L. Z. & Rappe, A. M. Shift-current bulk photovoltaic effect influenced by quasiparticle and exciton.

*Phys. Rev. B***101**, 045104 (2020). - 49.
Hu, H. et al. Quantum electronic stress: density-functional-theory formulation and physical manifestation.

*Phys. Rev. Lett.***109**, 055501 (2012). - 50.
Peng, B. et al. Sub-picosecond photo-induced displacive phase transition in two-dimensional MoTe

_{2}.*npj 2D Mater. Appl.***4**, 14 (2020). - 51.
Si, C. et al. Photoinduced vacancy ordering and phase transition in MoTe

_{2}.*Nano Lett.***19**, 3612–3617 (2019). - 52.
Li, T. et al. Optical control of polarization in ferroelectric heterostructures.

*Nat. Commun.***9**, 3344 (2018). - 53.
Chen, N.-K. et al. Optical subpicosecond nonvolatile switching and electron-phonon coupling in ferroelectric materials.

*Phys. Rev. B***102**, 184115 (2020). - 54.
Mieghem, P. V. Theory of band tails in heavily doped semiconductors.

*Rev. Mod. Phys.***64**, 755–793 (1992). - 55.
Zhou, J., Mao, S. & Zhang, S. Noncontacting optostriction driven anisotropic and inhomogeneous strain in two-dimensional materials.

*Phys. Rev. Res.***2**, 022059(R) (2020). - 56.
King-Smith, R. D. & Vanderbilt, D. Theory of polarization of crystalline solids.

*Phys. Rev. B***47**, 1651–1654 (1993). - 57.
Resta, R. Macroscopic polarization in crystalline dielectrics: the geometric phase approach.

*Rev. Mod. Phys.***66**, 899–915 (1994).

## 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

### Affiliations

### 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

## 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

## 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/.

## 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

Received:

Accepted:

Published: