Ultrafast simultaneous manipulation of multiple ferroic orders through nonlinear phonon excitation

Introduction

Identifying efficient pathways to control the quantum phases of matter is a key challenge for the development of new materials with functional properties. In magnetoelectric multiferroics, the coupled magnetic and ferroelectric orders enable the control of magnetization using electric fields and vice versa, promising energy-efficient switching operations for data processing. Significant research efforts have been made to develop this concept within the last decade^1,2. BiFeO₃, shown in Fig. 1a, is a well-known multiferroic material at room temperature, whose ferroelectric polarization can be switched by an applied electric field, which is accompanied by a reorientation of its antiferromagnetic domain^3,4. While this magnetoelectric coupling is appealing for low-energy consumption spintronics^5,6, the timescale of the switching process is limited by the ramp-up time of the applied electromagnetic fields, generally to nanosecond timescales.

**Fig. 1: Control of multiferroic order through nonlinear phonon dynamics.**

In contrast, ultrashort laser pulses provide the possibility to generate large electric field strengths on sub-picosecond timescales. Recent experimental realizations of intense terahertz (THz) and mid-infrared (MIR) pulses with peak electric fields exceeding several megavolts per centimeter (MV cm⁻¹) have opened up new pathways for dynamical materials engineering^7,8. When the electric field component of a laser pulse is coupled resonantly to the dipole moment of an optical phonon mode, coherent lattice vibrations can be driven with amplitudes exceeding the harmonic regime, making anharmonic contributions of the interatomic potential-energy surface integral parts of the dynamics. A particular manifestation of this anharmonicitiy is nonlinear phononic rectification, a lattice analog to optical rectification, in which nonlinear coupling between phonon modes leads to a quasistatic distortion of the crystal structure^9,10,11,12, creating transient crystal geometries that are not accessible in equilibrium. Examples of nonlinear phononic rectification have achieved to manipulate and induce ferroelectric^{13,14,15,16,17} and magnetic order^{18,19,20,21,22,23,24} individually. It has so far remained unexplored however, whether light can manipulate multiple ferroic order parameters simultaneously in a single domain.

Here, we coherently excite a high-frequency fully-symmetric (A₁) phonon mode in BiFeO₃ with intense MIR pulses and we probe the evolution of ferroelectricity and antiferromagnetism using time-resolved second-harmonic generation (SHG) at room temperature. We find that both the ferroelectric and antiferromagnetic contributions to the SHG signal are enhanced on a sub-picosecond timescale. We perform phonon dynamics simulations supported by density functional theory calculations that predict a rectification of the A₁ modes arising from both three-phonon coupling and cubic anharmonicities, whose polar atomic displacements increase the magnitude of the ferroelectric contribution to the second-order nonlinear susceptibility tensor. Comparison with experimental data confirms that the total change in SHG signal contains a significant proportion of antiferromagnetic contribution. Our calculations indicate that the rectified displacements along the A₁ modes lead to an enhancement of the exchange interaction between the iron spins, which we interpret as a stabilization of the antiferromagnetic state (Fig. 1(b)).

Results

Phonon-driven enhancement of second-harmonic generation

We study a single-crystal sample of BiFeO₃, cut and polished along the pseudocubic (110) plane. The ferroelectric polarization P_FE is oriented along the body diagonal, [1(bar{1})1], which is parallel to the polished surface, as illustrated in Fig. 1(a). Previous works show that the lattice, electronic, and magnetic degrees of freedom in BiFeO₃ are highly sensitive to external fields^{25,26,27,28,29,30}. The large optical nonlinearities in BiFeO₃³¹ enables optical rectification that leads to coherent THz emission, which is a sensitive probe of the ferroelectric polarization³². We observed a pronounced angular dependence of the emitted THz field amplitude that verifies that our sample is predominantly in a single domain (Supplementary Note 1).

To control structural parameters critical to the ferroic properties, we used intense MIR pulses to excite the high-frequency A₁ mode in BiFeO₃ at 15.7 THz³³. We used frequency-tunable MIR pulses with an excitation fluence of 12 mJ cm⁻² to couple to the electric dipole moment of the A₁ mode that is oriented along the ferroelectric polarization direction. Time-resolved SHG was used to probe the MIR-induced changes in the ferroic order parameters. The SHG intensity is proportional to the square of the light-induced nonlinear polarization, I(2ω) ∝ ∣P(2ω)∣². The transient state is characterized by the light-induced change in SHG intensity, ΔI(2ω), with the probe polarization oriented parallel to the light scattering plane (angle ϕ = 0°) and the p-polarization detection configuration. A schematic of the pump-probe geometry is shown in Fig. 2a, indicating the atomic displacements induced by the high-frequency A₁ mode that is primarily composed of oxygen-ion motion. In Fig. 2b, we show the characterization of the MIR pump pulse and the polarization dependence of ΔI(2ω). When the polarization of the MIR pulse is oriented along P_FE, the SHG intensity is enhanced by 1.5%, whereas a polarization of the MIR pulse perpendicular to P_FE yields no pump-induced change. When tuning the center frequency of the MIR pump pulse away from the A₁ phonon frequency, the pump-induced ΔI(2ω) decays rapidly, as shown in Fig. 2c, indicating a resonant phonon-driven effect. Furthermore, the maximum transient change in polarization scales linearly with the pump fluence, shown in Fig. 2d, indicating a mechanism based on nonlinear phononic rectification and ionic Raman scattering³⁴.

**Fig. 2: Phonon-driven enhancement of second-harmonic generation.**

SHG is a sensitive probe for ferroelectricity and magnetic ordering, as the symmetry of the material is reflected by the coefficients of the SHG susceptibility tensor that is given by

$${bf{P}}(2omega )propto left({chi }^{(i)}+{chi }^{(c)}right):{bf{E}}(omega )otimes {bf{E}}(omega ),$$

(1)

where E(ω) is the electric field component of the laser pulse and χ⁽ⁱ⁾ and χ^(c) respectively are the ferroelectric and magnetic contributions to the SHG susceptibility tensor^35,36,37. The nonzero elements of the SHG tensors are determined by the point-group symmetry of the crystal. At room temperature, the crystal lattice of BiFeO₃ is rhombohedral with a 3m point-group symmetry (space group R3c). In the bulk, the antiferromagnetically aligned spins exhibit an additional long-range cycloidal modulation (Fig. 1a) that can be oriented along three equivalent wave vectors, commonly denoted by (τ₁, τ₂, τ₃). The wave vectors τ_i are perpendicular to P_FE^25,38. Our sample is predominantly in a single ferroelectric domain, as well as a single spin-cycloid domain, as discussed in Supplementary Notes 1 and 3. The uniqueness of the spin-cycloid wave vector lifts the three-fold rotation symmetry around P_FE, and effectively lowers the symmetry of the system to monoclinic, described by point group m. Because the ferroelectric and magnetic properties determine the point-group symmetry, information about them is intrinsically encoded in the symmetry and magnitude of the SHG tensor components that can be probed by varying the polarization of the probe beam.

In Fig. 3a, b, we show the probe-polarization dependence of the SHG signal in equilibrium. For the p-polarization detection, the SHG polar plot has two lobes at 0° and 180°, respectively. This can be fit using nearly exclusively the χ⁽ⁱ⁾ tensor components (ferroelectric contribution, grey area). For the s-polarization detection, the SHG signal has four lobes, with a clear asymmetry along the 60° and 240° directions. This asymmetry is due to the nonzero χ^(c) tensor components (antiferromagnetic contribution, red area) arising from the symmetry lowering in the single spin-cycloid domain. Similar features were observed in BiFeO₃ thin films, in which the orientation of the asymmetric SHG component was used to identify different antiferromagnetic domains³⁹. In Fig. 3c, d, we show the SHG polar plots in the phonon-driven state at the maximum of the pump-probe response. For both the p– and s-polarization detections, (Delta {I}^{max }) is anisotropic and primarily positive with varying probe polarization. The maximum enhancement for both (Delta {I}_{{rm{p}}}^{max }/{I}_{{rm{p}}}) and (Delta {I}_{{rm{s}}}^{max }/{I}_{{rm{s}}}) is approximately (1.5 ± 0.4)%. Intriguingly, the dominant contributions to the pump-induced p– and s-signal changes are from different χ tensors. In the p-detection channel, the maximum enhancement in SHG intensity is along 0° and 180°, given by the ferroelectric SHG tensor (grey area in Fig. 3c). In contrast, in the s-detection channel, the maximum enhancement in SHG intensity is along 60° and 240°, dominated by the antiferromagnetic SHG tensor (red area in Fig. 3d). In order to shed light on the physical mechanisms that lead to this simultaneous enhancement of ferroelectric and antiferromagnetic SHG enhancement, we perform simulations of the nonlinear phonon dynamics that follow the coherent excitation of the high-frequency A₁ mode by the MIR pulse.

**Fig. 3: Static and time-resolved SHG polarimetry.**

Simulations of nonlinear phonon dynamics

We evaluate the nonlinear vibrational response and the corresponding dynamics of the lattice polarization and the magnetic interactions in a semi-classical oscillator model, for which we compute the input parameters from first principles using density functional theory. The noncentrosymmetric point-group symmetry of BiFeO₃ allows for quadratic-linear three-phonon couplings of the resonantly driven high-frequency A₁ mode (calculated eigenfrequency 15.3 THz) to the three low-frequency A₁ modes (calculated eigenfrequencies 8.8, 7.4, and 4.8 THz) present in the system, which we denote as ({A}_{1}^{h}) and ({A}_{1}^{l}) in the following. The coupling is described by an interaction potential of (V=c{Q}_{{A}_{1}^{h}}^{2}{Q}_{{A}_{1}^{l}}), where Q is the phonon amplitude and c the coupling coefficient, and can be utilized for nonlinear phononic rectification: when the ({A}_{1}^{h}) mode is coherently driven, it enacts a unidirectional force on the ({A}_{1}^{l}) modes. This force induces a quasistatic displacement of the atoms along the eigenvectors of the ({A}_{1}^{l}) modes that is proportional to the mean-squared amplitude of the ({A}_{1}^{h}) mode, (langle {Q}_{{A}_{1}^{l}}rangle propto langle {Q}_{{A}_{1}^{h}}^{2}rangle). Because the ({A}_{1}^{h}) mode is driven linearly by the electric field component of the pump, ({Q}_{{A}_{1}^{h}} sim E), phononic rectification is quadratic in the pump field and therefore linear in the pump fluence. Furthermore, as the high-frequeny A₁ mode is driven with a large amplitude, it can experience self-rectification through cubic-order anharmonicity, described by the potential term (V=tilde{c}{Q}_{{A}_{1}^{h}}^{3}). Because all A₁ modes are infrared active, the rectification induces a quasistatic electric dipole moment and therefore an electric polarization in addition to the ferroelectric polarization. Details of the evolution of the phonon amplitudes according to the oscillator model and the density functional theory calculations of the coupling coefficients are provided in the Materials and Methods and in the Supplementary Note 4.

In Fig. 4a, we show the displacements of the ions for the sum of the rectified components of all four A₁ modes. The cations, bismuth and iron, move in opposite directions along the [1(bar{1})1] direction. In Fig. 4b, we show the transient polarization induced by the rectified components of all four A₁ modes, as well as the relative change of lattice polarization compared to the equilibrium ferroelectric polarization of 90 μC cm⁻². The rectifications of the individual A₁ modes induce opposing polarizations due to the opposite signs of their mode effective charges that lead to only a small net change in polarization, which can be explained by the opposite net motion of cations shown in Fig. 4a.

**Fig. 4: Simulations of nonlinear phonon dynamics and phonon-induced ferroelectric SHG.**

In order to predict the pump-induced change in SHG intensity, we calculate and compare the ferroelectric contribution to the SHG tensor components of both the equilibrium and transiently distorted lattice structures at the maximum of the vibrational response. In Fig. 4c, d, we show the polar plots of the change in SHG intensity, ΔI, for the p– and s-detection channels that we obtain from the modified SHG tensor components due to nonlinear phononic rectification. The calculations well reproduce the experimental features (grey areas in Fig. 3c, d), which confirms the ferroelectric contribution to the SHG enhancement of ΔI_p, and supports the experimental finding that the enhancement of ΔI_s along 60° and 240° must stem from the antiferromagnetic χ^(c) tensor (red area in Fig. 3d).

Finally, we investigate the impact of the driven structural state on the magnetic order by computing the changes in the exchange interactions and magnetocrystalline anisotropy of the iron spins that arise from the quasistatic displacements of the A₁ modes. We show the time evolution of the mean values of the magnetic interactions in Fig. 5. Intriguingly, we obtain a net enhancement of the magnetic interactions in the driven state, similar to the changes in electric polarization. At the maximum rectified amplitude of the lattice (Supplementary Note 4 Fig. S4), we find that the nearest neighbor exchange (Fig. 5a) increases in magnitude by ∣ΔJ₁∣ = 7 μeV or 0.03% and the next nearest neighbor exchange (Fig. 5b) increases by ∣ΔJ₂∣ = 0.7 μeV or 0.4%, whereas the magnetocrystalline anisotropy (Fig. 5c) decreases by ∣ΔK∣ = 0.6 μeV or 0.1%. The net change of J₁ is an order of magnitude larger than the changes of J₂ and K. While the overall net changes are small, they indicate a trend towards a more robust magnetic ordering, consistent with the enhanced antiferromagnetic SHG signal that we measure. Note however that we employed a simplified model, neglecting the spin spiral in the magnetic order of BiFeO₃ that requires further terms in the spin Hamiltonian, such as the Dzyaloshinskii-Moriya interaction. Consequently, these numbers should be used with caution and more comprehensive first-principles studies will be necessary in the future to discuss the quantitative relation between the modifications of the magnetic interactions and the antiferromagnetic SHG signal, which are beyond our current methodology.

**Fig. 5: Simulations of phonon-dependent magnetic interactions.**

Discussion

We have demonstrated a simultaneous enhancement of ferroelectric and antiferromagnetic second-harmonic generation in BiFeO₃ by resonant excitation of a high-frequency fully-symmetric phonon mode. We stress that the pump-induced SHG intensity changes are not a multi-domain effect, as the single-crystal sample studied here is predominantly in a single ferroelectric and spin-cycloid domain. The SHG signals from minor domains are negligible and do not match the experimental data, see Supplementary Note 3. We can further exclude the contribution from MIR-induced SHG, which involves the third-order susceptibility tensor, χ⁽³⁾, in the form given by^40,41

$$P(2omega )propto {chi }^{(3)}{E}^{MIR}E(omega )E(omega )$$

The χ⁽³⁾ signal is spectrally filtered (see Materials and Methods) and the χ⁽³⁾ process is spatially separated from the χ⁽²⁾ process due to the 45° angle between the SHG probe and the MIR pump beams. In addition, the χ⁽³⁾ process would result in a linear pump field dependence for ΔP(2ω), corresponding to a square-root dependence on the pump fluence. However, we observe a linear pump fluence dependence (Fig. 2d).

The enhancement of the ferroelectric contribution can be well reproduced by our calculations of the χ⁽ⁱ⁾ tensor-component modulations due to the transient lattice distortion induced by nonlinear phononic rectification. Intriguingly, the lattice polarizations of the A₁ phonon have opposite signs due to the counter-phase motion of the cations along the direction of ferroelectric polarization, P_FE, and nearly cancel each other out. In turn, they enhance the SHG tensor components directly, therefore leading to a larger P(2ω) in comparison to the equilibrium state that our measurements pick up. While our theory captures the ferroelectric SHG features well, a direct calculation of the antiferromagnetic SHG tensor is beyond the scope of this work. Our calculations indicate however an overall trend towards enhanced exchange interactions by the rectified displacements of the A₁ modes. This mechanism of magnetic-interaction control through phonons has been shown for nonpolar magnets before^18,19,20,42. We expect this modification to lead to a more robust magnetic ordering and to changes of the antiferromagnetic χ^(c) tensor components, consistent with our experimental findings. The fully-symmetric A₁ modes by definition preserve the symmetry of the system and while they leave the relative distances between the iron ions mostly unchanged, they modify the relative distances between the iron and oxygen ions. This could further affect the Dzyaloshinskii-Moriya interaction that is responsible for the noncollinear arrangement of spins in BiFeO₃. Therefore, more comprehensive first-principles studies clarifying the quantitative relation between the modifications of the magnetic interactions and the antiferromagnetic SHG signal are required in the future.

Finally, we briefly touch on other possible mechanisms of phonon-induced magnetization control, which could be at play in addition to the enhancement of the exchange interaction. One possibility concerns phonon-induced symmetry breaking of the orbital configuration, which has recently been shown to stabilize magnetization above the ordering temperature in a transition-metal ferromagnet⁴³. Due to the fully-symmetric nature of the A₁ modes, we can likely rule out this mechanism in our case. Another possibility concerns coherent magnon excitation through coupling to the phonons, which has in the past been described as ionic Raman scattering by magnons^44,45,46. While ionic Raman scattering of the A₁ modes by the magnons might be possible in BiFeO₃, only circularly polarized phonons would lead to a unidirectional force on the spins⁴⁷. All A₁ modes involved in the dynamics are linearly polarized, and we therefore also rule out this mechanism for our case. Furthermore, the spectrum of our MIR pulse is far away from the electromagnon resonances in BiFeO₃^48,49, which allows us to exclude a coupling to the spins through the dipole moment of the spin cycloid. With the data and calculations given thus far, we, therefore, find it most likely that the pronounced enhancement of the nearest neighbor exchange interaction leads to a stabilization of the antiferromagnetic order and to a change of the χ^(c) components.

To conclude, our results show that ultrashort pulses in the THz and MIR spectral range can be used to address multiple ferroic order parameters at once, utilizing the coupling to nonlinearly driven phonon modes. We have demonstrated simultaneous ferroelectric and antiferromagnetic SHG enhancement for the case of the prototypical multiferroic material BiFeO₃. Our first-principles calculations confirm that the increased ferroelectric SHG signal arises from phonon-mediated increases of the ferroelectric χ⁽²⁾ tensor. They further indicate an enhancement of the nearest neighbor exchange interaction, consistent with our finding of increased antiferromagnetic SHG. We expect that this approach will be applicable generically to materials exhibiting magnetoelectric coupling. The ferroic orders in BiFeO₃ form independently below the respective Curie and Néel temperatures at 1100 and 640 K and are only weakly coupled. In contrast, the ferroic orders are interdependent and emerge together in so-called type-II multiferroics, but these materials often suffer from very low multiferroic ordering temperatures⁵⁰. An intriguing question to address in the future therefore concerns whether both ferroic orders can be stabilized simultaneously above the multiferroic ordering temperature through coherent phonon excitation, for example through an enhancement of the magnetoelectric coupling between the two ferroic orders. This problem can be regarded in the same broad context as the light-induced stabilization of superconductivity, ferroelectricity, and magnetism demonstrated in recent experiments^{15,16,43,51,52}.

Materials and methods

Crystal growth

The BiFeO₃ single crystal was grown using a traveling solvent technique in a Laser Floating Zone (LFZ) furnace⁵³. Bi₂O₃ (99.975%, Alfa Aesar) and Fe₂O₃ (99.995%, Alfa Aesar) powder are mixed in molar ratio 1:1 and 4:1 for the feed rod and the traveling solvent, respectively. The feed rod material was sintered at 800 °C for 20 h with one intermediated grinding, and the solvent was sintered at 700 °C for 5 h. A 0.15 g solvent pellet was placed on top of the seed rod before starting the floating zone growth. The growth was performed in air, and the growth rate was 1 mm h⁻¹. The as-grown crystal is oriented by a Laue X-ray camera. The ferroelectric domain structure has been characterized by optical microscope, circular interference differential contrast, and piezoresponse force microscopy imaging (Supplementary Note 3).

MIR pump-SHG probe setup

The fundamental beam for SHG is from a Ti:sapphire laser with a repetition rate of 1 kHz, energy of 30 nJ, and a duration of 40 fs. A schematic representation of the setup is shown in Fig. 6. A half-wave plate (HWP) was used to rotate the polarization of the fundamental beam. The beam was focused to 120 micron in diameter on the sample. The SHG signal was spectrally filtered by a band-pass filter (with a 400 nm center frequency and a 10 nm FWHM bandwidth), spatially filtered by an iris, and then detected by a photo-multiplier tube. The spectral and spatial filtering processes eliminate the detection of the MIR-induced SHG^40,41. A Glan-Taylor prism was used to select the desired polarization of the second-harmonic signal. All optical measurements were performed at room temperature, with a 45° incidence in reflection geometry.

**Fig. 6: Schematic of the MIR pump-SHG probe experiment, with the orientations of the (x, y, z) laboratory frame and the (X₁, X₂, X₃) crystal frame indicated.**

MIR pump pulses with tunable wavelengths between 13 and 18.5 microns, and a duration of 250 fs were generated by optical parametric amplification and subsequent difference frequency generation in a GaSe crystal. The MIR pulses were focused on to the sample at normal incidence, with a diameter of 220 microns and a maximum excitation fluence of 12 mJ cm⁻².

Fits of the SHG signal

At room temperature, BiFeO₃ crystallizes in rhombohedral symmetry with point group 3m and space group R3c. We set X₃ as the ferroelectric axis and X₂ as the normal to the mirror plane of the trigonal structure (Fig. 6). In a single antiferromagnetic domain, the long-range cycloidal modulation of the spin breaks the three-fold rotation symmetry around P_FE, reducing the symmetry of the material to monoclinic with point group m. Details of the SHG tensors can be found in Supplementary Note 2. The fit reveals that the cycloidal spin modulation direction is along X₂.

Theory of nonlinear phonon dynamics

We use a phenomenological oscillator model to evaluate the nonlinear vibrational response of the A₁ modes and the corresponding lattice-polarization dynamics following the excitation by an ultrashort MIR pulse, for which we compute all parameters from first principles using density functional theory. Specifically, we solve the equations of motion

$${ddot{Q}}_{i}+{kappa }_{i}{dot{Q}}_{i}+{partial }_{{Q}_{i}}V=0,$$

(2)

where Q_i are the normal mode coordinates (amplitudes) and κ_i the phenomenological linewidths of the A₁ modes in BiFeO₃, i ∈ {A₁}. V = V₀ + V_int is the potential energy, consisting of an equilibrium potential, V₀, and an interaction potential, V_int. We expand V₀ up to fourth order in anharmonic and nonlinear coupling terms,

$${V}_{0}=sum _{i={{A}_{1}}}frac{{Omega }_{i}^{2}}{2}{Q}_{i}^{2}+sum _{i,j,k={{A}_{1}}}{c}_{ijk}{Q}_{i}{Q}_{j}{Q}_{k}+sum _{i,j,k,l={{A}_{1}}}{d}_{ijkl}{Q}_{i}{Q}_{j}{Q}_{k}{Q}_{l},$$

(3)

where Ω_i is the phonon eigenfrequency, and c_ijk and d_ijkl are the anharmonic and nonlinear coupling coefficients. Because the A₁ modes are fully symmetric, all permutations of modes are allowed by symmetry in the anharmonic terms. We take into account all single-mode anharmonicities (c_iii, d_iiii), but only nonlinear couplings of the high-frequency A₁(15.3 THz) mode, which is resonantly driven, to the other A₁ modes of the system. See Supplementary Note 4 for the density functional theory calculations of the coupling coefficients.

The interaction potential, V_int, contains the coupling of the electric dipole moment of the infrared-active phonon A₁ modes, p_i, to the electric field component of the laser pulse, E(t), and can be written as

$${V}_{{rm{int}}}=sum _{i={{A}_{1}}}{p}_{i}E(t)=sum _{i={{A}_{1}}}{Z}_{i}{Q}_{i}E(t),$$

(4)

where the linear polarization of the laser pulse is assumed to be aligned along the ferroelectric polarization direction of the crystal. Q_i is the phonon amplitude as evaluated in Eq. (2) and Z_i is the mode effective charge given by ({Z}_{i}={sum }_{n}{Z}_{n}^{* }{q}_{n,i}/sqrt{{M}_{n}}). ({Z}_{n}^{* }) is the Born effective charge tensor, q_n,i the phonon eigenvector, and M_n the ionic mass of ion n, and the index n runs over all ions in the unit cell. We model the electric field component of the laser pulse as (E(t)={E}_{0}exp [-{t}^{2}/(2{({tau }_{0}/sqrt{8ln 2})}^{2})]cos ({omega }_{0}t)) and set its parameters to the experimental values: a full width at half maximum (FWHM) duration of 250 fs, a central frequency of 16.9 THz, and a fluence of 12 mJ cm⁻² (peak electric field E₀ = 5.3 MV cm⁻¹). Because ω₀ is far off resonance from the three low-frequency A₁ modes, only the contribution of the A₁(15.3 THz) mode plays a role in Eq. (4). Finally, the time-dependent electric polarization, P(t), induced by the displacements of the ions along the coordinates of the infrared-active A₁ modes, is given by

$$P(t)=sum _{i={{A}_{1}}}{p}_{i}(t)/{V}_{c},$$

(5)

where V_c is the volume of the unit cell.

We next compute the phonon-induced modulation of the SHG tensor components. Due to the lack of inversion symmetry in the crystal, any infrared-active mode can linearly modify the components of χ⁽²⁾ as

$${chi }_{{rm{pumped}}}^{(2)}(t)={chi }_{{rm{equilibrium}}}^{(2)}+mathop{sum}limits_{j = {{A}_{1}}}frac{partial {chi }_{j}^{(2)}}{partial {Q}_{j}}{Q}_{j}(t).$$

(6)

Here, ({chi }_{{rm{equilibrium}}}^{(2)}) is the equilibrium second-order susceptibility and ({partial }_{{Q}_{j}}{chi }_{j}^{(2)}) is the phonon-induced change. We compute both the equilibrium and nonequilibrium tensor components for the ferroelectric contribution, χ⁽ⁱ⁾, from first principles, see Supplementary Note 4.

We finally compute the time-dependent modifications of the exchange interactions, J₁ and J₂, as well as the magnetocrystalline anisotropy, K, as a result of the rectified A₁ mode displacements

$${J}_{1/2}(t)={J}_{1/2}^{{rm{equilibrium}}}+sum _{j={{A}_{1}}}frac{partial {J}_{1/2}}{partial {Q}_{j}}{Q}_{j}(t),$$

(7)

$$({E}_{{n}_{110}}-{E}_{{n}_{1bar{1}1}})(t)equiv K(t)={K}^{{rm{equilibrium}}}+sum _{j={{A}_{1}}}frac{partial K}{partial {Q}_{j}}{Q}_{j}(t).$$

(8)

We calculate ({J}_{1/2}^{{rm{equilibrium}}}), K^equilibrium, ({partial }_{{Q}_{j}}{J}_{1/2}), and ({partial }_{{Q}_{j}}K) from first principles, see Supplementary Note 4. The mean change in magnetic interactions over time is given by their dependence on the rectified phonon coordinate, e.g., 〈K(t)〉 ≡ K[〈Q(t)〉].