Terahertz optics-driven phase transition in two-dimensional multiferroics

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₃ (refs ^1,2), BiFeO₃ (ref. ³), and SrTiO₃ (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⁶.

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 }$ (ω₀ 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⁷. 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₃, 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¹⁰ and α-SnTe¹¹ 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⁻¹) 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, ε₀ 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_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.

Fig. 1: LPTL inducing ferroic order switch.

(Left panel) Different ferroic orders (FE and FE′) are subject to structural operator $\hat {\cal{O}}$, which are degenerate in energy. (Right panel) Under LPTL, FE′ has lower GFE, corresponding to a FE→FE′ transition.

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¹⁵ 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²¹.

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 theory²², 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²³

$$\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 (τ = Γ⁻¹) of the vibration. For simplicity, we choose a universal value of Γ = 4 cm⁻¹, 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 Pmn2₁ 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²⁸, 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.

Fig. 2: Atomic and THz optical responses of monolayer β-GeSe.

a Geometric structure of monolayer β-GeSe, with dashed rectangle indicating simulation unit cell. b Calculated real part of optical susceptibility and absorbance function with respect to incident energy. The gray shaded region indicates phononic range (<8 THz), above which the responses are mainly contributed from electronic subsystem. The absorbance in the phononic range is enlarged 20 times for clarity reasons. The subscripts “1” and “2” denote the Cartesian coordinates “x” and “y”, respectively.

In addition to ferroelastic order, the lack of centrosymmetry indicates a spontaneous electric polarization P_s. Consistent with previous works²⁷, we compute its value to be 0.16 nC m⁻¹ (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²⁶, 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₁₁(ω = 5.38 THz) = 0.07%] for the x-LPTL, and two absorbance peaks [A₂₂(ω = 1.2 THz) = 0.26% and A₂₂(ω = 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⁻¹ (corresponding to an intermediate LPTL intensity of 5 × 10⁹ W cm⁻², 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⁻²) (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² (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

$$\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)

Fig. 3: Direction-dependent EELS spectra of FE0-β-GeSe.

The 65 keV electron beam is traveling at a distance of b = 25 nm parallel to the β-GeSe plane, along x- and y-directions (inset).

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⁻². Now we consider another 2D time-reversal invariant multiferroic material, monolayer α-SnTe, whose parental geometry is C₄ 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 Pmn2₁ 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⁻¹, along the puckered direction. These results agree well with previous works^39,40.

Fig. 4: Atomic and THz optical responses of monolayer α-SnTe.

a Geometric structure of FE1-α-SnTe monolayer. FE2 is 90°-rotated. b Calculated optical susceptibility and absorbance function with respect to incident energy. The gray shaded region indicates phononic range. The phononic absorbance is enlarged for clarity reasons. c EELS spectra when an electron beam with kinetic energy of 65 keV travels along the x- and y-directions, at a distance of 25 nm from the α-SnTe plane.

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⁴¹, 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₁₁ = 0.3%), and two peaks for the y-LPTL (at 1.4 and 1.9 THz, with A₂₂ = 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⁹ W cm⁻² 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⁻²). When one increases the LPTL intensity to 2.2 × 10¹⁰ W cm⁻² (0.42 V nm⁻¹), 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⁴².

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⁻¹), 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

$$\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 functions⁴⁶ 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⁴⁷. 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⁻¹ (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⁻¹ (Fig. 5e). From its momentum distribution, the dominate contribution comes from ±Λ [=(0, ±0.55, 0) Å⁻¹] 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$.

Fig. 5: DFT calculated second harmonic generation.

Real part, imaginary part, and total frequency-dependent in-plane electronic contributed SHG spectra of monolayers a, b FE0-β-GeSe and d, e FE1-α-SnTe. Horizontal axis indicates incident energy ℏω. Inset shows k-resolved total SHG susceptibility distribution at ℏω = 0.82 eV (GeSe) and 0.36 eV (SnTe), indicated by the red circles. Polarization anisotropy of absolute SHG susceptibilities of c GeSe and f SnTe, where red and green curves indicate SHG susceptibility that are parallel and perpendicular to the polarization direction, respectively.

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⁴⁸. Unfortunately, this is extremely computational challenging because of nonlocal interactions included and very dense k-mesh needed. In the previous calculations⁴⁷, 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⁻¹ (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⁴⁹ 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³⁰. 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₃ thin films deposited on transition metal dichalcogenide monolayers⁵², 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₃, which is controlled by the synergistic effects of excited-state energy profile, atomic motion velocity, and the dephasing of excited states⁵³.

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². Compared with 3D materials or thick films (such as BaTiO₃ 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⁵⁴. 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)⁵⁵. 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¹⁴ and performed in the Vienna Ab initio Simulation Package (VASP)¹⁵ with generalized gradient approximation (GGA) treatment of exchange and correlation functional in the solid-state PBE form (PBEsol)¹⁶. 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¹⁷ 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¹⁸. The electron and force component convergence criteria are set to be 10^–7 eV and 0.001 eV Å⁻¹, 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¹⁹ is used to include van der Waals interactions, which has been demonstrated to yield good results on the energy barrier³⁹. 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²⁰ to calculate vibrational modes and Born effective charge of each ion, with the help of the Phonopy code²¹.