2021

These authors contributed equally to this work.

[2]Dominik M. Juraschek These authors contributed equally to this work.

[1,5]Wanzheng Hu

1]Department of Physics, Boston University, 590 Commonwealth Avenue, Boston,
02215, Massachusetts, USA

2]School of Physics and Astronomy, Tel Aviv University, 6997801, Tel Aviv,
Israel

3]Max Planck Institute for the Structure and Dynamics of Matter, Center for
Free-Electron Laser Science (CFEL), Luruper Chaussee 149, Hamburg, 22761,
Germany

4]Rutgers Center for Emergent Materials and Department of Physics and
Astronomy, Rutgers University, 136 Frelinghuysen Road, Piscataway, 08854, New
Jersey, USA

5]Division of Materials Science and Engineering, Boston University, 590
Commonwealth Avenue, Boston, 02215, Massachusetts, USA

# Ultrafast simultaneous manipulation of multiple ferroic orders through
nonlinear phonon excitation

Daniel A. Bustamante Lopez<EMAIL_ADDRESS><EMAIL_ADDRESS>Michael
Fechner<EMAIL_ADDRESS>Xianghan Xu<EMAIL_ADDRESS>Sang-
Wook Cheong<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [  [

###### Abstract

Recent experimental studies have demonstrated the possibility of utilizing
strong terahertz pulses to manipulate individual ferroic orders on pico- and
femtosecond timescales. Here, we extend these findings and showcase the
simultaneous manipulation of multiple ferroic orders in BiFeO3, a material
that is both ferroelectric and antiferromagnetic at room temperature. We find
a concurrent enhancement of ferroelectric and antiferromagnetic second-
harmonic generation (SHG) following the resonant excitation of a high-
frequency fully-symmetric phonon mode. Based on first-principles calculations
and phenomenological modeling, we ascribe this observation to the inherent
coupling of the two ferroic orders to the nonequilibrium distortions induced
in the crystal lattice by nonlinearly driven phonons. Our finding highlights
the potential of nonlinear phononics as a technique for manipulating multiple
ferroic order parameters at once. In addition, this approach provides a
promising avenue to studying the dynamical magnetization and polarization
behavior, as well as their intrinsic coupling, on ultrashort timescales.

###### keywords:

multiferroics, nonlinear phononics, ultrafast dynamics

##  1 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,
and significant research efforts have been made to develop this concept within
the last decade Eerenstein ; SpaldinRamesh . BiFeO3, 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 ZhaoRamesh ;
LebeugleGukasov . While this magnetoelectric coupling is appealing for low-
energy consumption spintronics HeronRamesh2011 ; HeronRamesh2014 , the
timescale of the switching process is limited by the ramp-up time of the
applied electromagnetic fields, generally to nanosecond timescales.

In contrast, ultrashort laser pulses provide the possibility to generate large
electric field strengths on sub-picosecond timescales, and recent experimental
realizations of intense terahertz (THz) and mid-infrared (MIR) pulses with
peak electric fields exceeding several megavolts per centimetre (MV/cm) have
opened up new pathways for dynamical materials engineering Liu2017 ;
Vicario2020 . 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 Foerst2011 ; subedi:2014 ; fechner:2016 ;
juraschek:2017 , creating transient crystal geometries that are not accessible
in equilibrium. Examples of nonlinear phononic rectification have achieved to
manipulate and induce ferroelectric subedi:2015 ; Mankowsky_2:2017 ; Nova2019
; Li2019 ; Chen2022  and magnetic order Fechner2018 ; Khalsa2018 ; Rodriguez-
Vega2020 ; Disa2020 ; Afanasiev2021 ; Stupakiewicz2021 ; Disa2021
individually. It has so far remained unexplored however, whether light can
control multiple ferroic order parameters simultaneously in a single domain.

Here, we coherently excite a high-frequency fully-symmetric ($A_{1}$) phonon
mode in BiFeO3 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_{1}$ 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, which we interpret as a stabilization of the antiferromagnetic
state (Fig. 1b).

Figure 1: Control of multiferroic order through nonlinear phonon dynamics.
a, Rhombohedral crystal structure of BiFeO3 (space group $R3c$), in which the
ferroelectric polarization ($\mathbf{P}_{\mathrm{FE}}$) is along the stretched
body-diagonal of the pseudocubic unit cell, identified as $[1\bar{1}1]$ here.
The spins between nearest-neighbor iron ions are antiferromagnetically aligned
(red arrows) and form a spin cycloid with a wavelength of 62 nm along the
normal vector $\boldsymbol{\tau}$ of the mirror plane (yellow) Lee2008 ;
Ratcliff2016 . In a single-magnetic domain crystal, the spin cycloid breaks
the three-fold rotation symmetry around $\mathbf{P}_{\mathrm{FE}}$. b, (Left)
Initial state with the equilibrium ferroelectric polarization and a thermal
broadening of the iron-spin alignment. (Right) The lattice distortion
following the resonant excitation of a $A_{1}$ phonon leads to a change in
polarization and to an improved alignment of the iron spins.

##  2 Results

###  2.1 Phonon-driven enhancement of second-harmonic generation

We study a single-crystal sample of BiFeO3, cut and polished along the
pseudocubic (110) plane. The ferroelectric polarization
$\mathbf{P}_{\mathrm{FE}}$ is oriented along the body diagonal, [1$\bar{1}$1],
which is parallel to the polished surface, as illustrated in Fig. 1a. Previous
works show that the lattice, electronic, and magnetic degrees of freedom in
BiFeO3 are highly sensitive to external fields Lee2008 ; Choi2009 ; Kundys2010
; Chen2015 ; Rubio-Marcos2018 ; Liou2019 . The large optical nonlinearities in
BiFeO3 Kumar2008  enables optical rectification that leads to coherent THz
emission, which is a sensitive probe of the ferroelectric polarization
Talbayev2008 . We observed a pronounced angular dependence of the emitted THz
field amplitude that verifies that our sample is predominantly in a single
domain (Supplementary Information Section 1).

To control structural parameters critical to the ferroic properties, we used
intense MIR pulses to excite the high-frequency $A_{1}$ mode in BiFeO3 at 15.7
THz Palai2010 . We used frequency-tunable MIR pulses with an excitation
fluence of 12 mJ/cm2 to couple to the electric dipole moment of the $A_{1}$
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\omega)\propto\lvert
P(2\omega)\rvert^{2}$. The transient state is characterized by the light-
induced change in SHG intensity, $\Delta I(2\omega)$, with the probe
polarization oriented parallel to the light scattering plane (angle
$\phi=0^{\circ}$) 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_{1}$ mode that are
primarily composed of oxygen-ion motion. In Fig. 2b, we show the
characterization of the MIR pump pulse and the polarization dependence of
$\Delta I(2\omega)$. When the polarization of the MIR pulse is oriented along
$\mathbf{P}_{\mathrm{FE}}$, the SHG intensity is enhanced by 1.5%, whereas a
polarization of the MIR pulse perpendicular to $\mathbf{P}_{\mathrm{FE}}$
yields no pump-induced change. When tuning the center frequency of the MIR
pump pulse away from the $A_{1}$ phonon frequency, the pump-induced $\Delta
I(2\omega)$ 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
Neugebauer2021 .

Figure 2: Phonon-driven enhancement of second-harmonic generation. a, We used
MIR pulses with a center frequency of 16.9 THz (69 meV) to drive the high-
frequency $A_{1}$ phonon in BiFeO3 at normal incidence. A probe beam at 800 nm
was used for SHG with a 45∘ incidence in reflection geometry. b, (Top panel)
Electric field component of the MIR pump and Fourier spectrum (inset). (Middle
panel) The MIR pump induces a 1.5% transient enhancement of the SHG intensity
when polarized along $\mathbf{P}_{\mathrm{FE}}$. (Bottom panel) Rotating the
MIR pump polarization by $90^{\circ}$ results in no pump-probe signal. c,
Dependence of the SHG signal on the center frequency of the MIR pump pulse.
The signal becomes nonzero as the frequency of the pump pulse approaches that
of the high-frequency $A_{1}$ mode, indicating a phonon-driven mechanism. d,
Linear dependence of the light-induced change in second-harmonic polarization
on the pump fluence, consistent with nonlinear phononic rectification and
ionic Raman scattering Neugebauer2021 .

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

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

where $\mathbf{E}(\omega)$ is the electric field component of the laser pulse
and $\chi^{(i)}$ and $\chi^{(c)}$ respectively are the ferroelectric and
magnetic contributions to the SHG susceptibility tensor Denev2011 ; Fiebig2005
. The nonzero elements of the SHG tensors are determined by the point-group
symmetry of the crystal. At room temperature, the crystal lattice of BiFeO3 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
$(\boldsymbol{\tau}_{1},\boldsymbol{\tau}_{2},\boldsymbol{\tau}_{3})$. The
wave vectors $\boldsymbol{\tau}_{i}$ are perpendicular to
$\mathbf{P}_{\mathrm{FE}}$ Ratcliff2016 ; Lee2008 . Our sample is
predominantly in a single ferroelectric domain, as well as a single spin-
cycloid domain Lee2008 . The uniqueness of the spin-cycloid wave vector lifts
the three-fold rotation symmetry around $\mathbf{P}_{\mathrm{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.

Figure 3: Static and time-resolved SHG polarimetry. Polar plots of the
equilibrium SHG intensity, $I$, are shown for a, $p$-polarization detection
and b, $s$-polarization detection (open circles). Polar plots of the pump-
induced change of the SHG intensity, $\Delta I$, are shown for c,
$p$-polarization detection and d, $s$-polarization detection at the maximum of
the pump-probe response (open squares). In all plots, solid lines are fits of
the SHG intensities, composed of the ferroelectric contribution (grey areas)
and the antiferromagnetic contribution (red areas). Both the ferroelectric and
antiferromagnetic contributions show an enhancement following the resonant
phonon excitation.

In Figs. 3a and 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 fitted using nearly
exclusively the $\chi^{(i)}$ 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 $\chi^{(c)}$ tensor components (antiferromagnetic
contribution, red area) arising from the symmetry lowering in the single spin-
cycloid domain. Similar features were observed in BiFeO3 thin films, in which
the orientation of the asymmetric SHG component was used to identify different
antiferromagnetic domains Chauleau2017 . In Figs. 3c and 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^{max}_{p}/I_{p}$ and $\Delta
I^{max}_{s}/I_{s}$ is approximately 1.5%. Intriguingly, the dominant
contributions to the pump-induced $p$\- and $s$-signal changes are from
different $\chi$ 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_{1}$ mode by the MIR pulse.

###  2.2 Simulations of nonlinear phonon dynamics

We evaluate the nonlinear vibrational response and the corresponding lattice-
polarization dynamics in a semi-classical oscillator model, for which we
compute the input parameters from first principles using density functional
theory. The point-group symmetry of BiFeO3 allows for quadratic-linear three-
phonon couplings of the resonantly driven high-frequency $A_{1}$ mode
(calculated eigenfrequency 15.3 THz) to the three low-frequency $A_{1}$ modes
(calculated eigenfrequencies 8.8, 7.4, and 4.8 THz) present in the system,
which we denote as $A^{h}_{1}$ and $A^{l}_{1}$ in the following. The coupling
is described by an interaction potential of
$V=cQ_{A^{h}_{1}}^{2}Q_{A^{l}_{1}}$, where $Q$ is the phonon amplitude and $c$
the coupling coefficient, and can be utilized for nonlinear phononic
rectification: when the $A^{h}_{1}$ mode is coherently driven, it enacts a
unidirectional force on the $A^{l}_{1}$ modes. This force induces a
quasistatic displacement of the atoms along the eigenvectors of the
$A^{l}_{1}$ modes that is proportional to the mean-squared amplitude of the
$A^{h}_{1}$ mode, $\langle Q_{A^{l}_{1}}\rangle\propto\langle
Q^{2}_{A^{h}_{1}}\rangle$. Furthermore, as the high-frequeny $A_{1}$ 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^{h}_{1}}^{3}$. Because all $A_{1}$ 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 Methods section and in the Supplementary
Information Section 4.

Figure 4:  Simulations of nonlinear phonon dynamics and phonon-induced
ferroelectric SHG. a, Transient displacements of the ions along the rectified
components of all four $A_{1}$ modes. Displacements are exaggerated for
clarity. Besides the change in the oxygen octahedron, the two cations, bismuth
and iron, move in opposite directions along the [1$\bar{1}$1] direction,
marked by black and white arrows. b, Polarization induced by the rectified
components of the four $A_{1}$ modes, and percentile change of the electric
polarization with respect to the equilibrium value ($90$ $\mu$C/cm2).
Calculated pump-induced changes in the SHG signal for the ferroelectric
contribution are shown for c, $p$-polarization detection and d,
$s$-polarization detection at the maximum response of the phonon amplitudes.
The shape and magnitude of the polar plot in $p$-polarization fits excellently
to the experimental data. In $s$-polarization, the shape of the experimental
data for the ferroelectric contribution could be reproduced, however with an
overestimation of the magnitude.

In Fig. 4a, we show the displacements of the ions for the sum of the rectified
components of all four $A_{1}$ 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_{1}$
modes, as well as the relative change of lattice polarization compared to the
equilibrium ferroelectric polarization of $90$ $\mu$C/cm2. The rectifications
of the individual $A_{1}$ 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.

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 (see Supplementary Information Section 4 for
computational details). In Figure 4c and d, we show the polar plots of the
change in SHG intensity, $\Delta 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 and d), which confirms the ferroelectric
contribution to the SHG enhancement of $\Delta I_{p}$, and supports the
experimental finding that the enhancement of $\Delta I_{s}$ along 60∘ and 240∘
must stem from the antiferromagnetic $\chi^{(c)}$ tensor (red area in Fig.
3d).

##  3 Discussion

We have demonstrated a simultaneous enhancement of ferroelectric and
antiferromagnetic second-harmonic generation in BiFeO3 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 Lee2008 . The SHG signals from minor domains are
negligible and do not match the experimental data, see Supplementary
Information Section 3. We can further exclude the contribution from MIR-
induced SHG, which involves the third-order susceptibility tensor,
$\chi^{(3)}$.The $\chi^{(3)}$ signal is spectrally filtered (see Methods) and
the $\chi^{(3)}$ process is spatially separated from the $\chi^{(2)}$ process
due to the 45∘ angle between the SHG probe and the MIR pump beams. In
addition, the $\chi^{(3)}$ process would result in a quadratic pump fluence
dependence for $\Delta P(2\omega)$, not the linear relation observed here
(Fig. 2d).

The enhancement of the ferroelectric contribution can be well reproduced by
our calculations of the $\chi^{(i)}$ tensor-component modulations due to the
transient lattice distortion induced by nonlinear phononic rectification.
Intriguingly, the lattice polarizations of the $A_{1}$ phonon have opposite
signs due to the counter-phase motion of the cations along the direction of
ferroelectric polarization, $\mathbf{P}_{\mathrm{FE}}$, and nearly cancel each
other out. In turn, they enhance the SHG tensor components directly, therefore
leading to a larger $\mathbf{P}(2\omega)$ in comparison to the equilibrium
state that our measurements pick up. While our theory captures the
ferroelectric SHG features, the microscopic mechanism responsible for the
enhancement of the antiferromagnetic contribution remains an open question,
and we provide a number of possible explanations below that will have to be
tested in the future.

A first possibility concerns phonon-induced changes of the exchange
interaction, which have been described for a variety of magnetic systems in
the past Fechner2018 ; Khalsa2018 ; Rodriguez-Vega2020 ; Padmanabhan2022 . The
quasistatic distortion induced through nonlinear phononic rectification
modifies the effective bond lengths and therefore potentially the exchange
coupling that is inherently dependent on the distance between the atoms. As a
result, this modification could lead to changes of the antiferromagnetic
$\chi^{(c)}$ tensor components, similar to the ferroelectric case. The fully-
symmetric $A_{1}$ modes that we drive by definition preserve the symmetry of
the system and furthermore leave the relative distances between the iron ions
mostly unchanged, however the relative distances between the iron and oxygen
ions change, which affects the Dzyaloshinskii-Moriya interaction that is
responsible for the noncollinear arrangement of spins in BiFeO3. A second
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 Disa2023 . Due to
the fully-symmetric nature of the $A_{1}$ modes, we can likely rule out this
mechanism in our case. A third possibility concerns coherent magnon excitation
through coupling to the phonons, which has in the past been described as ionic
Raman scattering by magnons nova:2017 ; juraschek2:2017 ; Juraschek2020_3 .
While ionic Raman scattering of the $A_{1}$ modes by the magnons might be
possible in BiFeO3, only circularly polarized phonons would lead to a
unidirectional force on the spins Dominik2022 . All $A_{1}$ 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 BiFeO3 DeSousa2008 ; Cazayous2008 , 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 a modification of the exchange interactions leads to a
stabilization of the antiferromagnetic order and to a change of the
$\chi^{(c)}$ components.

##  4 Conclusion

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 BiFeO3, but we expect
that this approach will be applicable generically to materials exhibiting
magnetoelectric coupling. The ferroic orders in BiFeO3 form independently
below the respective Curie and Néel temperatures at 1100 and 640 K. In so-
called type-II multiferroics, in contrast, the ferroic orders are
interdependent and emerge together, but these materials often suffer from very
low multiferroic ordering temperatures Cheong2007 . An intriguing question
therefore arises, whether both ferroic orders can be stabilized simultaneously
above the multiferroic ordering temperature through coherent phonon
excitation. This problem can be regarded in the same broad context as the
light-induced stabilization of superconductivity, ferroelectricity, and
magnetism demonstrated in recent experiments mankowsky:2014 ; Mitrano2016 ;
Nova2019 ; Li2019 ; Disa2023 .

##  5 Methods

###  5.1 Crystal growth

The BiFeO3 single crystal was grown using a traveling solvent technique in a
Laser Floating Zone (LFZ) furnace Ito2011 . Bi2O3 (99.975%, Alfa Aesar) and
Fe2O3 (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 hours with one intermediated grinding,
and the solvent was sintered at 700∘C for 5 hours. 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.

###  5.2 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. 5. A half-wave plate (HWP) was
used to rotate the polarization of the fundamental beam. The beam was focused
SHG signal was spectrally filtered by a band-pass filter, spatially filtered
by an iris, and then detected by a photo-multiplier tube. 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.

MIR pump pulses with tunable wavelengths between 13 and 18.5 micron, 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
micron and a maximum excitation fluence of 12 mJ/cm2.

Figure 5: Schematic of the MIR pump-SHG probe experiment, with the
orientations of the $(x,y,z)$ laboratory frame and the $(X_{1},X_{2},X_{3})$
crystal frame indicated. A half-wave plate (HWP) rotates the polarization of
the 800 nm pulses. A dichroic mirror (not shown) reflects the second-harmonic
generation signal and removes the residual 800 nm signal. A band-pass filter
(10 nm bandwidth) removes signals from higher-order optical processes. An
analyzer selects either the vertically polarized ($s$) or horizontally
polarized ($p$) SHG signal.

###  5.3 Fits of the SHG signal

At room temperature, BiFeO3 crystallizes in rhombohedral symmetry with point
group ${3m}$ and space group $R3c$. We set $X_{3}$ as the ferroelectric axis
and $X_{2}$ as the normal to the mirror plane of the trigonal structure (Fig.
5). In a single antiferromagnetic domain, the long-range cycloidal modulation
of the spin breaks the three-fold rotation symmetry around
$\mathbf{P}_{\mathrm{FE}}$, reducing the symmetry of the material to
monoclinic with point group $m$. Details of the SHG tensors can be found in
Supplementary Information Section 2. The fit reveals that the cycloidal spin
modulation direction is along $X_{2}$.

###  5.4 Theory of nonlinear phonon dynamics

We use a phenomenological oscillator model to evaluate the nonlinear
vibrational response of the $A_{1}$ 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 $\kappa_{i}$
the phenomenological linewidths of the $A_{1}$ modes in BiFeO3,
$i\in\\{A_{1}\\}$. $V=V_{0}+V_{\mathrm{int}}$ is the potential energy,
consisting of an equilibrium potential, $V_{0}$, and an interaction potential,
$V_{\mathrm{int}}$. We expand $V_{0}$ up to fourth order in anharmonic and
nonlinear coupling terms,

$V_{0}=\sum\limits_{i=\\{A_{1}\\}}\frac{\Omega_{i}^{2}}{2}Q_{i}^{2}+\sum\limits_{i,j,k=\\{A_{1}\\}}c_{ijk}Q_{i}Q_{j}Q_{k}+\sum\limits_{i,j,k,l=\\{A_{1}\\}}d_{ijkl}Q_{i}Q_{j}Q_{k}Q_{l},$
(3)

where $\Omega_{i}$ is the phonon eigenfrequency, and $c_{ijk}$ and $d_{ijkl}$
are the anharmonic and nonlinear coupling coefficients. Because the $A_{1}$
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_{1}$(15.3 THz) mode, which is resonantly driven, to the other $A_{1}$ modes
of the system. See Supplementary Information Section 4 for the density
functional theory calculations of the coupling coefficients.

The interaction potential, $V_{\mathrm{int}}$, contains the coupling of the
electric dipole moment of the infrared-active phonon $A_{1}$ modes, $p_{i}$,
to the electric field component of the laser pulse, $E(t)$, and can be written
as

$V_{\mathrm{int}}=\sum\limits_{i=\\{A_{1}\\}}p_{i}E(t)=\sum\limits_{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{8\ln 2})^{2})]\cos(\omega_{0}t)$ and
set its parameters to the experimental values: a full width at half maximum
(FWHM) duration of $\tau_{0}=250$ fs, a central frequency of $\omega_{0}=16.9$
THz, and a fluence of 12 mJ/cm2 (peak electric field $E_{0}=5.3$ MV/cm).
Because $\omega_{0}$ is far off resonance from the three low-frequency $A_{1}$
modes, only the contribution of the $A_{1}$(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_{1}$
modes, is given by

$P(t)=\sum\limits_{i=\\{A_{1}\\}}p_{i}(t)/V_{c},$ (5)

where $V_{c}$ is the volume of the unit cell.

We finally 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 $\chi^{(2)}$ as

$\chi_{\text{pumped}}^{(2)}(t)=\chi_{\text{equilibrium}}^{(2)}+\sum_{j=\\{A_{1}\\}}\frac{\partial\chi_{j}^{(2)}}{\partial
Q_{j}}Q_{j}(t).$ (6)

Here, $\chi_{\text{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, $\chi^{(i)}$, from first principles, see
Supplementary Information Section 4.

Supplementary Information

Details of the SHG analysis and density functional theory calculations can be
found in the Supplementary Information.

Acknowledgments We acknowledge the helpful discussions with Liuyan Zhao, Ankit
Disa, and Rui Zu. D.A.B.L. and W.H acknowledge support from the U.S.
Department of Energy, Office of Science, Office of Basic Energy Sciences Early
Career Research Program under Award Number DE-SC-0021305. The work at Rutgers
was supported by the W.M. Keck Foundation. D.M.J. was supported by Tel Aviv
University.

Author contributions

D.A.B.L. and D.M.J. contributed equally to this work. W.H. conceived the
project, designed and performed the pump-probe measurements, and analyzed the
data. D.A.B.L. performed data fitting and modeling. D.M.J. performed the
density functional theory calculations and simulations of the nonlinear phonon
dynamics. M.F. performed the density functional theory calculations of the SHG
tensor components. X.X. and S-W.C. grew the sample. W.H. and D.M.J. wrote the
paper with input from all other authors.

## Declarations

All authors declare no competing interests.

## References

  *   * (1) Eerenstein, W., Mathur, N.D., Scott, J.F.: Multiferroic and magnetoelectric materials. Nature 442, 759–765 (2006). https://doi.org/10.1038/nature05023
  * (2) Spaldin, N.A., Ramesh, R.: Advances in magnetoelectric multiferroics. Nat. Mater. 18, 203–212 (2019). https://doi.org/10.1038/s41563-018-0275-2
  * (3) Zhao, T., Scholl, A., Zawaliche, F., Lee, K., Barry, M., Doran, A., Cruz, M.P., Chu, Y.H., Ederer, C., Spaldin, N.A., Das, R.R., Kim, D.M., Baek, S.H., Eom, C.B., Ramesh, R.: Electrical control of antiferromagnetic domains in multiferroic BiFeO3 films at room temperature. Nat. Mater. 5, 823–829 (2006). https://doi.org/10.1038/nmat1731
  * (4) Lebeugle, D., Colson, D., Forget, A., Viret, M., Bataille, A.M., Gukasov, A.: Electric-field-induced spin flop in BiFeO3 single crystals at room temperature. Phys. Rev. Lett. 100, 227602 (2008). https://doi.org/10.1103/PhysRevLett.100.227602
  * (5) Heron, J.T., Trassin, M., Ashraf, K., Gajek, M., He, Q., Yang, S.Y., Nikonov, D.E., Chu, Y.-H., Salahuddin, S., Ramesh, R.: Electric-field-induced magnetization reversal in a ferromagnet-multiferroic heterostructure. Phys. Rev. Lett. 107, 217202 (2011). https://doi.org/10.1103/PhysRevLett.107.217202
  * (6) Heron, J.T., Bosse, J.L., He, Q., Gao, Y., Trassin, M., Ye, L., Clarkson, J.D., Wang, C., Liu, J., Salahuddin, S., Ralph, D.C., Schlom, D.G., Íñiguez, J., Huey, B.D., Ramesh, R.: Deterministic switching of ferromagnetism at room temperature using an electric field. Nature 516, 370–373 (2014). https://doi.org/10.1038/nature14004
  * (7) Liu, B., Bromberger, H., Cartella, A., Gebert, T., Först, M., Cavalleri, A.: Generation of narrowband, high-intensity, carrier-envelope phase-stable pulses tunable between 4 and 18 THz. Opt. Lett. 42, 129–131 (2017). https://doi.org/10.1364/OL.42.000129
  * (8) Vicario, C., Trisorio, A., Allenspach, S., Rüegg, C., Giorgianni, F.: Narrow-band and tunable intense terahertz pulses for mode-selective coherent phonon excitation. Appl. Phys. Lett. 117, 101101 (2020). https://doi.org/10.1063/5.0015612
  * (9) Först, M., Manzoni, C., Kaiser, S., Tomioka, Y., Tokura, Y., R., M., Cavalleri, A.: Nonlinear phononics as an ultrafast route to lattice control. Nat. Phys. 7, 854–856 (2011). https://doi.org/10.1038/nphys2055
  * (10) Subedi, A., Cavalleri, A., Georges, A.: Theory of nonlinear phononics for coherent light control of solids. Phys. Rev. B 89, 220301 (2014). https://doi.org/10.1103/PhysRevB.89.220301
  * (11) Fechner, M., Spaldin, N.A.: Effects of intense optical phonon pumping on the structure and electronic properties of yttrium barium copper oxide. Phys. Rev. B 94, 134307 (2016). https://doi.org/10.1103/PhysRevB.94.134307
  * (12) Juraschek, D.M., Fechner, M., Spaldin, N.A.: Ultrafast Structure Switching through Nonlinear Phononics. Phys. Rev. Lett. 118, 054101 (2017). https://doi.org/10.1103/PhysRevLett.118.054101
  * (13) Subedi, A.: Proposal for ultrafast switching of ferroelectrics using midinfrared pulses. Phys. Rev. B 92, 214303 (2015). https://doi.org/10.1103/PhysRevB.92.214303
  * (14) Mankowsky, R., von Hoegen, A., Först, M., Cavalleri, A.: Ultrafast Reversal of the Ferroelectric Polarization. Phys. Rev. Lett. 118, 197601 (2017). https://doi.org/10.1103/PhysRevLett.118.197601
  * (15) Nova, T., Disa, A., Fechner, M., Cavalleri, A.: Metastable ferroelectricity in optically strained SrTiO3. Science 364, 1075–1079 (2019). https://doi.org/10.1126/science.aaw4911
  * (16) Li, X., Qiu, T., Zhang, J., Baldini, E., Lu, J., Rappe, A.M., Nelson, K.A.: Terahertz field–induced ferroelectricity in quantum paraelectric SrTiO3. Science 364, 1079–1082 (2019). https://doi.org/10.1126/science.aaw4913
  * (17) Chen, P., Paillard, C., Zhao, H.J., Íñiguez, J., Bellaiche, L.: Deterministic control of ferroelectric polarization by ultrafast laser pulses. Nat. Commun. 13, 2566 (2022). https://doi.org/%****␣BFO.tex␣Line␣675␣****10.1038/s41467-022-30324-5
  * (18) Fechner, M., Sukhov, A., Chotorlishvili, L., Kenel, C., Berakdar, J., Spaldin, N.A.: Magnetophononics: Ultrafast spin control through the lattice. Phys. Rev. Materials 2, 064401 (2018). https://doi.org/10.1103/PhysRevMaterials.2.064401
  * (19) Khalsa, G., Benedek, N.A.: Ultrafast optically induced ferromagnetic/anti-ferromagnetic phase transition in GdTiO3 from first principles. npj Quantum Mater. 3, 15 (2018). https://doi.org/10.1038/s41535-018-0086-3
  * (20) Rodriguez-Vega, M., Lin, Z.-X., Leonardo, A., Ernst, A., Chaudhary, G., Vergniory, M.G., Fiete, G.A.: Phonon-mediated dimensional crossover in bilayer CrI3. Phys. Rev. B 102, 081117 (2020). https://doi.org/10.1103/physrevb.102.081117
  * (21) Disa, A.S., Fechner, M., Nova, T.F., Liu, B., Först, M., Prabhakaran, D., Radaelli, P.G., Cavalleri, A.: Polarizing an antiferromagnet by optical engineering of the crystal field. Nat. Phys. 16, 937–941 (2020). https://doi.org/10.1038/s41567-020-0936-3
  * (22) Afanasiev, D., Hortensius, J.R., Ivanov, B.A., Sasani, A., Bousquet, E., Blanter, Y.M., Mikhaylovskiy, R.V., Kimel, A.V., Caviglia, A..D.: Ultrafast control of magnetic interactions via light-driven phonons. Nat. Mater. 20, 607 (2021). https://doi.org/10.1038/s41563-021-00922-7
  * (23) Stupakiewicz, A., Davies, C.S., Szerenos, K., Afanasiev, D., Rabinovich, K.S., Boris, A.V., Caviglia, A., Kimel, A.V., Kirilyuk, A.: Ultrafast phononic switching of magnetization. Nat. Phys. 17, 489 (2021). https://doi.org/10.1038/s41567-020-01124-9
  * (24) Disa, A.S., Nova, T.F., Cavalleri, A.: Engineering crystal structures with light. Nat. Phys. 17, 1087–1092 (2021). https://doi.org/10.1038/s41567-021-01366-1
  * (25) Lee, S., Choi, T., Ratcliff II, W., Erwin, R., Cheong, S.-W., Kiryukhin, V.: Single ferroelectric and chiral magnetic domain of single-crystalline BiFeO3 in an electric field. Phys. Rev. B 78, 100101 (2008). https://doi.org/%****␣BFO.tex␣Line␣825␣****10.1103/PhysRevB.78.100101
  * (26) Ratcliff, W., Lynn, J.W., Kiryukhin, V., Jain, P., Fitzsimmons, M.R.: Magnetic structures and dynamics of multiferroic systems obtained with neutron scattering. npj Quantum Materials 1, 16003 (2016). https://doi.org/10.1038/npjquantmats.2016.3
  * (27) Choi, T., Lee, S., Choi, Y.J., Kiryukhin, V., Cheong, S.-W.: Switchable ferroelectric diode and photovoltaic effect in BiFeO3. Science 324, 63 (2009). https://doi.org/10.1126/science.116863
  * (28) Kundys, B., Viret, M., Colson, D., Kundys, D.O.: Light-induced size changes in BiFeO3 crystals. Nat. Mater. 9, 803 (2010). https://doi.org/10.1038/nmat2807
  * (29) Chen, F., Goodfellow, J., Liu, S., Grinberg, I., Hoffmann, M.C., Damodaran, A.R., Zhu, Y., Zalden, P., Zhang, X., Takeuchi, I., Rappe, A.M., Martin, L.W., Wen, H., Lindenberg, A.M.: Ultrafast terahertz gating of the polarization and giant nonlinear optical response in BiFeO3 thin films. Adv. Mater. 27, 6371 (2015). https://doi.org/10.1002/adma.201502975
  * (30) Rubio-Marcos, F., Ochoa, D.A., Del Campo, A., García, M.A., Castro, G.R., Fernández, J.F., García, J.E.: Reversible optical control of macroscopic polarization in ferroelectrics. Nature Photon. 12, 29–32 (2018). https://doi.org/10.1038/s41566-017-0068-1
  * (31) Liou, Y.-D., Chiu, Y.-Y., Hart, R.T., Kuo, C.-Y., Huang, Y.-L., Wu, Y.-C., Chopdekar, H.-J. R. V.and Liu, Tanaka, A., Chen, C.-T., Chang, C.-F., Tjeng, L.H., Cao, Y., Nagarajan, V., Chu, Y.-H., Chen, Y.-C., Yang, J.-C.: Deterministic optical control of room temperature multiferroicity in BiFeO3 thin films. Nat. Mater. 18, 580–587 (2019). https://doi.org/10.1038/s41563-019-0348-x
  * (32) Kumar, A., Rai, R.C., Podraza, N.J., Denev, S., Ramirez, M., Chu, Y.-H., Martin, L.W., Ihlefeld, J., Heeg, T., Schubert, J., Schlom, D.G., Orenstein, J., Ramesh, R., Collins, R.W., Musfeldt, J.L., Gopalan, V.: Linear and nonlinear optical properties of BiFeO3. Appl. Phys. Lett. 92, 121915 (2008). https://doi.org/10.1063/1.2901168
  * (33) Talbayev, D., Lee, S., Cheong, S.-W., Taylor, A.J.: Terahertz wave generation via optical rectification from multiferroic BiFeO3. Appl. Phys. Lett. 93, 212906 (2008). https://doi.org/10.1063/1.3036526
  * (34) Palai, R., Schmid, H., Scott, J.F., Katiyar, R.S.: Raman spectroscopy of single-domain multiferroic BiFeO3. Phys. Rev. B 81, 064110 (2010). https://doi.org/10.1103/PhysRevB.81.064110
  * (35) Neugebauer, M.J., Juraschek, D.M., Savoini, M., Engeler, P., Boie, L., Abreu, E., Spaldin, N.A., Johnson, S.L.: Comparison of coherent phonon generation by electronic and ionic Raman scattering in LaAlO3. Phys. Rev. Research 3, 13126 (2021). https://doi.org/10.1103/physrevresearch.3.013126
  * (36) Denev, S.A., Lummen, T.T.A., Barnes, E., Kumar, A., Gopalan, V.: Probing ferroelectrics using optical second harmonic generation. J. Am. Ceram. Soc. 94, 2699–2727 (2011). https://doi.org/10.1111/j.1551-2916.2011.04740.x
  * (37) Fiebig, M., Pavlov, V.V., Pisarev, R.V.: Second-harmonic generation as a tool for studying electronic and magnetic structures of crystals: review. J. Opt. Soc. Am. B 22, 96–118 (2005). https://doi.org/10.1364/JOSAB.22.000096
  * (38) Chauleau, J.-Y., Haltz, E., Carrétéro, C., Fusil, S., Viret, M.: Multi-stimuli manipulation of antiferromagnetic domains assessed by second-harmonic imaging. Nat. Matter. 16, 803–807 (2017). https://doi.org/10.1038/nmat4899
  * (39) Padmanabhan, H., Poore, M., Kim, P.K., Koocher, N.Z., Stoica, V.A., Puggioni, D., Wang, H., Shen, X., Reid, A.H., Gu, M., Wetherington, M., Lee, S.H., Schaller, R.D., Mao, Z., Lindenberg, A.M., Wang, X., Rondinelli, J.M., Averitt, R.D., Gopalan, V.: Interlayer magnetophononic coupling in MnBi2Te4. Nat. Commun. 13, 1929 (2022). https://doi.org/10.1038/s41467-022-29545-5
  * (40) Disa, A.S., Curtis, J., Fechner, M., Liu, A., Hoegen, A., Först, M., Nova, T.F., Narang, P., Maljuk, A., Boris, A.V., Keimer, B., Cavalleri, A.: Photo-induced high-temperature ferromagnetism in YTiO3. Nature 617, 73 (2023). https://doi.org/10.1038/s41586-023-05853-8
  * (41) Nova, T.F., Cartella, A., Cantaluppi, A., Först, M., Bossini, D., Mikhaylovskiy, R.V., Kimel, A.V., Merlin, R., Cavalleri, A.: An effective magnetic field from optically driven phonons. Nat. Phys. 13, 132–137 (2017). https://doi.org/10.1038/nphys3925
  * (42) Juraschek, D.M., Fechner, M., Balatsky, A.V., Spaldin, N.A.: Dynamical multiferroicity. Phys. Rev. Mater. 1, 014401 (2017). https://doi.org/10.1103/PhysRevMaterials.1.014401
  * (43) Juraschek, D.M., Narang, P., Spaldin, N.A.: Phono-magnetic analogs to opto-magnetic effects. Phys. Rev. Research 2, 043035 (2020). https://doi.org/10.1103/PhysRevResearch.2.043035
  * (44) Juraschek, D.M., Neuman, T., Narang, P.: Giant effective magnetic fields from optically driven chiral phonons in 4f paramagnets. Phys. Rev. Research 4, 013129 (2022). https://doi.org/10.1103/PhysRevResearch.4.013129
  * (45) De Sousa, R., Moore, J.E.: Optical coupling to spin waves in the cycloidal multiferroic BiFe O3. Phys. Rev. B 77, 012406 (2008). https://doi.org/10.1103/PhysRevB.77.012406
  * (46) Cazayous, M., Gallais, Y., Sacuto, A., De Sousa, R., Lebeugle, D., Colson, D.: Possible observation of cycloidal electromagnons in BiFeO3. Phys. Rev. Lett. 101, 037601 (2008). https://doi.org/10.1103/PhysRevLett.101.037601
  * (47) Cheong, S.-W., Mostovoy, M.: Multiferroics: a magnetic twist for ferroelectricity. Nat. Mater. 6, 13–20 (2007). https://doi.org/10.1038/nmat1804
  * (48) Mankowsky, R., Subedi, A., Först, M., Mariager, S.O., Chollet, M., Lemke, H.T., Robinson, J.S., Glownia, J.M., Minitti, M.P., Frano, A., Fechner, M., Spaldin, N.A., Loew, T., Keimer, B., Georges, A., Cavalleri, A.: Nonlinear lattice dynamics as a basis for enhanced superconductivity in YBa2Cu3O6.5. Nature 516, 71–73 (2014). https://doi.org/10.1038/nature13875
  * (49) Mitrano, M., Cantaluppi, A., Nicoletti, D., Kaiser, S., Perucchi, A., Lupi, S., Di Pietro, P., Pontiroli, D., Riccò, M., Clark, S.R., Jaksch, D., Cavalleri, A.: Possible light-induced superconductivity in K3C60 at high temperature. Nature 530, 461–464 (2016). https://doi.org/10.1038/nature16522
  * (50) Ito, T., Ushiyama, T., Yanagisawa, Y., Kumai, R., Tomioka, Y.: Growth of highly insulating bulk single crystals of multiferroic BiFeO3 and their inherent internal strains in the domain-switching process. Cryst. Growth Des. 11, 5139–5143 (2011). https://doi.org/10.1021/cg201068m