Origin of Terahertz Soft-Mode Nonlinearities in Ferroelectric Perovskites

Soft modes are intimately linked to structural instabilities and are key for the understanding of phase transitions. The soft modes in ferroelectrics, for example, map directly the polar order parameter of a crystal lattice. Driving these modes into the nonlinear, frequency-changing regime with intense terahertz (THz) light fields is an efficient way to alter the lattice and, with it, the physical properties. However, recent studies show that the THz electric-field amplitudes triggering a nonlinear soft-mode response are surprisingly low, which raises the question on the microscopic picture behind the origin of this nonlinear response. Here, we use linear and two-dimensional terahertz (2D THz) spectroscopy to unravel the origin of the soft-mode nonlinearities in a strain-engineered epitaxial ferroelectric SrTiO 3 thin film. We find that the linear dielectric function of this mode is quantitatively incompatible with pure ionic or pure electronic motions. Instead, 2D THz spectroscopy reveals a pronounced coupling of electronic and ionic-displacement dipoles. Hence, the soft mode is a hybrid mode of lattice (ionic) motions and electronic interband transitions. We confirm this conclusion with model calculations based on a simplified pseudopotential concept of the electronic band structure. It reveals that the entire THz nonlinearity is caused by the off-resonant nonlinear response of the electronic interband transitions of the lattice-electronic hybrid mode. With this work, we provide fundamental insights into the microscopic processes that govern the softness that any material assumes near a ferroic phase transition. This knowledge will allow us to gain an efficient all-optical control over the associated large nonlinear effects.


I. INTRODUCTION
Soft modes are excitations that are intrinsically linked to a phase instability and mediate the phase transition, where they approach zero frequency [1,2]. Polar soft modes in ferroelectrics, for instance, carry an electric dipole moment and are therefore susceptible to external light fields. Thus, there has been extensive effort towards driving soft modes into the strongly nonlinear regime (a regime where the higher-order terms in the expansion of polarization or current density in the material dominate its response) with intense THz light fields, either to photoinduce ferroelec-tricity [3][4][5] or to alter lattice structures all-optically [6,7]. While most of these studies exploit the soft-mode nonlinearities, direct experimental insight into the very origin of soft-mode nonlinear response is rather limited. Some of the rare attempts towards this goal reveal that the nonlinear response of THz soft modes in ferroelectrics [8] or molecular crystals [9] occurs at surprisingly low THz electric-field amplitudes of approximately 50 kV=cm. Upon resonant excitation with a THz light field, this low-field nonlinear response manifests itself as a pronounced blueshift of the vibrational resonance and photon-echo signals at negative coherence times [9]. In the literature, two different microscopic pictures are used to explain the observations associated with the soft-mode nonlinearities. These pictures are diametrically opposed to each other-while in the first picture, the ionic motions govern the nonlinear response, in the second picture, it is the electronic motions. Both pictures, however, treat the respective motions in a classical manner and hence suffer from discrepancies, such that the fundamental origin of the THz-induced soft-mode nonlinearities in ferroelectric perovskites has remained a topic of intense debate for decades.
Before we tackle this issue, we briefly recall the two pictures and their discrepancies.
The ionic picture is based on the purely phononic origin of the soft-mode nonlinearity and relies unconditionally on the Born-Oppenheimer approximation. Therein, the lowfrequency motion of charged particles in a crystal or a molecule is entirely performed by ions on the potential energy surface of the electronic ground state. In such an approach, any THz nonlinearity is connected to the anharmonicity of the potential energy surface [10]. Such anharmonicity is the basis of classical molecular-dynamics simulations, which can describe the dynamics of nonlinearly coupled phonon modes and has been measured by using femtosecond x-ray diffraction methods [11]. When the coupling to the electromagnetic field is included, however, the purely phononic approach to the soft-mode dynamics gives unphysical values to certain parameters. These include non-natural Born effective charges [12] on the ions [8,13] and large zero-point motions [8]. Further, the purely phononic approach fails to reproduce the observed nonlinear signals in 2D THz experiments, such as the presence of photon-echo signals at negative coherence times [9].
The electronic picture, in contrast, involves both ions and electrons and is based on Cochran's core-shell model [1]. In this model, the electronic equations of motion are introduced in a classical manner, similar to the pure lattice motion of the ions in the crystal. Close to the phase transition, the soft-mode frequency gradually approaches zero, leading to an instability of the entire crystal. Concomitantly with this softening, the electronic contribution to the oscillator strength of the soft mode grows to values much larger than the contribution from its pure ionic motion. Such a redistribution of electronic oscillator strength into the vibrational system of the crystal is mediated by the coupling between the electronic-and ionic-displacement dipoles via the Lorentz field (cf. the Clausius-Mossotti relation [14,15]). The electronic redistribution has been demonstrated via femtosecond x-ray diffraction experiments performed both in paraelectric and/ or ferroelectric phases of various ferroelectric [16][17][18] and molecular crystals [19]. With classical equations for the electronic motion, this picture can reproduce soft-mode nonlinearities associated with slow collective motions, i.e., with frequencies far below the THz range, such as during thermal phase transitions [20,21]. However, as this picture neglects the intrinsic nonlinearity of the electronic oscillators and further ignores the kinetic energy of the electrons [22], it fails to account for the correct electronic nonlinearities associated with THz soft modes.
In this paper, our observation of nonlinear THz signalsin particular, the photon echo-shows that both the abovementioned pictures are quantitatively unphysical. This is mainly because the nonlinearity associated with the electronic motions has been neglected, which, however, is the dominant one at THz frequencies. Instead, we show that the soft mode is a hybrid mode of ionic motions and electronic interband transitions. We present an expanded hybrid picture where pure lattice motion of ions is described by perfect harmonic oscillators and the entire THz nonlinearity of the soft mode is dominated by the electronic oscillators. In our hybrid picture, the correct electronic nonlinearity is determined by the nonlinear modification of the electron's velocity, which is connected to the interband transition dipoles. The coupling between the electronic-and the ionic-displacement dipoles depends sensitively on the effective magnitude of these interband transition dipoles. Hence, our picture not only highlights the hybrid character of the soft mode but also leads us to unambiguously identify the origin of THz soft-mode nonlinearity, namely, the off-resonant nonlinear response of the electronic interband transitions. Our experimental finding that the electronic degrees of freedom in the soft-mode quasiparticle of ferroelectrics are responsible for its nonlinear response to resonant driving fields is valid for all condensed matter systems with elementary excitations around 1 THz and lower, e.g., liquid water. The role of electronic degrees of freedom for low-frequency modes in condensed matter has been underestimated so far, and our work can stimulate future research in this direction.

II. PRECISE FORMULATION OF THE OPEN QUESTIONS AND OUR EXPERIMENTAL STRATEGY
A. Mechanisms contributing to the nonlinear response of phonons in crystals We begin by commemorating the basic mechanisms determining the light-matter interaction of phonons in crystals. The electronic ground state of an (insulating) crystal determines the anharmonic potential surface on which the nuclear motions are performed. In the quantummechanical description, the Hamiltonian can be formulated by creation (a † ) and annihilation (a) operators of all normal modes of the crystal [23]. For the time being, we consider just one anharmonic mode, i.e., the soft mode whose anharmonic potential [shown in Fig. 1(a)], and its coupling to light, is described by where ℏ is the reduced Planck constant and ω soft is the softmode frequency. Here, η quantum mechanically describes the lowest order of anharmonic forces beyond the harmonic oscillator, which is that of the Morse potential [24]. Note that E loc ðtÞ is the local field that originates from the coupling between the electronic-and ionic-displacement dipoles via the Lorentz field, and P mac is the macroscopic electric polarization of phonons.
In the paraelectric phase [red potential in Fig. 1(a)], the ladder of vibrational energy eigenstates shows an increasing energy spacing for excited states that corresponds to a positive η > 0 in Eq. (1). In contrast, for the ferroelectric phase [blue potential in Fig. 1(a)], the spacing decreases, i.e., η < 0. According to the seminal work of Lax and Burstein [25], P mac is, in general, a nonlinear function of the phonon coordinate(s) that can be expanded in the latter: If several modes are considered, A, B, C, and D are tensors of the first, second, third, and fourth rank, respectively [26].
Lax and Burstein were the first to point to the fact that there is "… an appreciable deformation of the charge distribution about the atoms during lattice vibration … the deformation of the charge distribution is suggested to be a possible alternative explanation to that of anharmonic forces …." Thus, there are two different sources for the nonlinear phonon response, i.e., the anharmonic potential [ Fig. 1(a)] and the nonlinear deformation of electronic charge. The concept of nonlinear redistribution of electronic charge density is rather universal [27,28].
In the ionic picture, all nonlinear terms beyond the linear term with tensor B in the expansion [i.e., Eq. (2)] are neglected. With such an approximation, still compatible with the Born-Oppenheimer approximation, the soft mode of a paraelectric phase will show a nonlinear blueshift upon excitation, whereas that of a ferroelectric phase must show a redshift as dictated by its ladder of vibrational energy eigenstates [the blue potential curve in Fig. 1(a)]. Note that the terms in Eq. (2), which are nonlinear in the phonon coordinates (e.g., C and D), are incompatible with the Born-Oppenheimer approximation. This is because the nonlinear deformation of the electronic charge is an additional degree of freedom that is not contained in the wave function describing the nuclear motions only.

B. Coupling of ionic (nuclear) motions to electronic interband currents via local-field effects
The alternative derivation of the Clausius-Mossotti relation [14] provides a very simple physical picture of  dipole-dipole coupling in condensed matter via local-field effects. This approximative approach is essentially based on two assumptions: (i) All contributing dipoles can be approximated by point dipoles (i.e., by using δ functions), and (ii) the lattice vectors of the crystal are in a good approximation to those of a cubic crystal. The local field in our Eq. (1) is determined by the macroscopic electric field in the crystal and by all (time-dependent) localized dipoles d k ðtÞ with their respective volume density N k : There is a strong impact of local-field effects on the soft mode in perovskites. As sketched in Fig. 1(b), the electronic interband dipoles are located on the Ti-O bonds (here, of SrTiO 3 ), whereas the ionic-displacement dipoles are located on the nuclei. Both dipoles are strongly localized in real space, residing at spatially disjunct positions. Thus, Hannay's approximation [14] is fully applicable. The linear response of the material already contains valuable information on the relevance of dipole-dipole coupling between the electronic and the ionic degrees of freedom. For a measured frequency-dependent dielectric function ϵðωÞ, one can compare the scenarios, including [i.e., Eq. (4)] and neglecting [i.e., Eq. (5)] local-field effects: It turns out that the combination of a huge dielectric function at the soft-mode frequency and a considerable off-resonant electronic high-frequency ϵ ∞ leads to a predominant electronic contribution to the oscillator strength of the soft-mode quasiparticle. This leads to dramatic consequences for the nonlinear soft-mode response because the local driving field acting on dipoles [Eq. (3)] has a temporal structure that is distinctly different from the THz pulses, as illustrated by the green and orange transients in Fig. 1 Neglecting any local-field effects [i.e., the second term in Eq. (3)], the applied THz pulses themselves are the driving fields acting on all dipoles in the system, thereby selecting only the most resonant contributions to the nonlinear response-a typical situation considered in 2D spectroscopy on molecular systems in the liquid phase [29]. In such a situation, rephasing photon-echo signals can exclusively occur for those interaction sequences of THz pulses A and B with a first interaction with the electric field of the respective leading pulse followed by two interactions with the electric field of the respective trailing pulse [e.g., the ABB photon echo in the uppermost pulse sequence in Fig. 1(b)]. Note that in collinear 2D THz spectroscopy, the rephasing ABB photon echo can be uniquely demarcated by experimental means from other nonlinear signals, such as the nonrephasing signals with an ABB interaction sequence or the rephasing BAA photon echo, i.e., the second trace in Fig. 1(b). Here, the ABB photon echo vanishes because pulse A (green) cannot interact with the sample as a leading pulse anymore.
For the opposite situation, the local field experienced by individual dipoles is dominated by the electric field of all other dipoles; i.e., the second term dominates over the first term in Eq. (3). In this case, we have protracted driving fields on all localized dipoles as shown in the two lowermost traces of Fig. 1(b). Such driving fields will allow for the occurrence of both rephasing photon-echo signals, ABB and BAA, even in the respective reversed pulse sequence.

C. Large electronic nonlinearity at THz frequencies
If the oscillator strength of soft modes in perovskites is dominated by its electronic interband contribution, it is expected that its nonlinear response dominates, too. Here, the excursion of valence-band electrons in k space plays a significant role, as sketched in Fig. 1(c). As shown recently [18], the microscopic polarization density is not a uniquely defined physical quantity, whereas the microscopic current density is (in principle) a measurable physical quantity. Thus, it is convenient to consider the time integrals of electronic (j el ), Eq. (6), and ionic (j ions ) current densities, Eq. (7), as the dipole densities in Eq. (3): where N el and N ions are the densities of electronic and ionic oscillators, respectively. When driving a valence-band electron with an electric field at THz frequencies, its k vector makes an excursion in the Brillouin zone as sketched in Fig. 1(c). The k-dependent interband-transition dipole (red curve) now plays a crucial role for the off-resonant electronic interband nonlinearity. At a frequency of 1 THz, an electric-field amplitude of a few tens of kV/cm is sufficient to explore regions in k space with a distinctly smaller interband-transition dipole, thereby modulating the local-field effects in Eq. (3) dramatically. Such a reduction of dipole-dipole coupling leads to a nonlinear blueshift of the soft-mode resonance in both paraelectric and ferroelectric phases. we enlist below and briefly explain our experimental strategy to find the answers: (1) What is the portion of electronic degrees of freedom vs nuclear (ionic) motions in the soft-mode quasiparticle? An answer to this important question can be obtained from a detailed quantitative analysis of the oscillator strength of the soft mode in the linear optical regime.
(2) When entering the nonlinear regime, what is the dominating soft-mode nonlinearity in the Hamiltonian, Eq. (1): the anharmonicity of the potential or the deformation of the charge density distribution?
The nonlinear response of a ferroelectric perovskite will give a definite answer to this question: observation of redshift or blueshift upon soft-mode excitation. While the redshift implies potential anharmonicities, the blueshift is a signature of charge density deformation.
(3) Can we provide experimental evidence for the relevant local-field contributions to the soft-mode nonlinearity?
The occurrence of rephasing ABB photon-echo signals for both THz-pulse sequences (A-B vs B-A) in our 2D THz experiments will directly prove the dipole-dipole coupling between ionic soft-mode motions and electronic interband currents in the localized Ti-O bonds of SrTiO 3 , used for our experiments.
(4) Is the experimentally observed nonlinear response of such a hybrid mode in quantitative agreement with a theoretical model with reasonable parameters?
The existing models are insufficient and hence we will need a new hybrid model. We show that, in our model, even when neglecting the potential anharmonicity [i.e., setting η ¼ 0 in Eq. (1)], we can reproduce all the experimental results (i.e., the linear and nonlinear responses) in a quantitative way.

III. EXPERIMENT
A. Sample design As a model system, we need a ferroelectric perovskite that hosts soft modes at THz frequencies and is close to the ferroelectric phase transition at room temperature. We have chosen a ferroelectric perovskite, in particular, because its quasicubic unit cell facilitates the theoretical description of local-field effects, and the well-known anomalously large values of its Born effective charges [13] are nicely reproduced by state-of-the-art first-principles calculations concerning ferroelectric polarizations. This approach further ensures that our experimental findings can be reproduced for all perovskites listed in Table I of Ref. [13].
Thin films of SrTiO 3 (STO) grown on DyScO 3 (DSO) induce an epitaxial strain that makes the STO undergo a ferroelectric phase transition at room temperature and hosts THz soft modes [30][31][32][33]. Here, the thin-film setting provides the degree of freedom to control the transition temperature via epitaxial strain. The choice of STO is further motivated by its high susceptibility [8], which provides an ideal platform to demonstrate the admixture of electronic degrees of freedom to the soft-mode quasiparticle via the dipole-dipole coupling. We grow our (001)-oriented STO film of d ¼ 50 nm thickness on a (110)-oriented DSO substrate using pulsed laser deposition [34]. Figure 2(a) shows the x-ray peaks corresponding to the DSO substrate and the STO film. The substrateimposed tensile strain as verified by x-ray reciprocal space mapping, shown in Fig. 2(b), amounts to 0.99%. As a result, an in-plane ferroelectricity is induced in our film at room temperature [35] with a soft-mode frequency in the THz range. Further details on the sample preparation and strain characterization are provided in Secs. I and II of the Supplemental Material [36].

B. Linear and 2D THz spectroscopy
Linear THz spectroscopy is performed by using a single THz pulse generated by optical rectification of 120-fs pulses from a Ti:Sapphire amplified laser system at 800 nm in a 0.5-mm-thick (110)-oriented ZnTe crystal. The THz electric field transmitted through the sample is detected as a function of the so-called real (or sampling) time t by free-space electro-optic sampling [37]. A small fraction of the 800-nm pulse is used as the sampling beam. The electro-optic sampling is carried out in a 0.5-mm-thick (110)-oriented ZnTe crystal that is optically bonded on a 2-mm-thick (100)-oriented ZnTe crystal.
The 2D THz spectroscopy [38][39][40][41][42][43][44][45][46] is performed by generating two phase-locked THz pulses, A and B, separated by the delay time τ. The second THz pulse is generated by optical rectification of the 120-fs pulses at 800 nm in a 0.5-mm-thick (110) [9]. Contour plots of the individual transmitted THz pulses (A and B) are shown in Sec. III of the Supplemental Material [36]. All experiments are performed in an inert N 2 atmosphere at room temperature.

A. Electronic and ionic contributions to the soft mode in STO
To estimate the contribution of the electronic and the ionic motions to the soft-mode oscillator strength, we first evaluate the soft-mode dielectric function ϵ soft ðωÞ. The THz transients transmitted through a bare DSO (reference) and the STOjDSO (sample) are shown in Fig. 2(c). The transmittance TðωÞ is obtained by taking the ratio of the Fouriertransformed THz spectra of the sample [red curve in Fig. 2(d)] to that of the reference [black curve in Fig. 2(d)]. The transmittance spectrum, in Fig. 2(e), shows the softmode resonance [8,[31][32][33] [47,48]), the dielectric function ϵ soft ðωÞ ¼ ϵ 0 ðωÞ þ iϵ 00 ðωÞ is obtained from the complex transmittance TðωÞ using the relation [8] where c is the vacuum velocity of light and n is the refractive index of the substrate. Figure 2(f) shows the zero crossing of the real part and the peak of the imaginary part of ϵ soft ðωÞ at the soft-mode frequency. We fit the measured dielectric function to that of a damped harmonic oscillator having a Lorentzian profile [red and black dashed lines in Fig. 2(f)], with a frequency ω soft and a damping rate γ soft . This leads to where ω 2 ions ¼ Z 2 eff N=ðϵ 0 M red Þ is the ionic plasma frequency, Z eff is the Born effective charge, N is the density of soft-mode oscillators, and ϵ 0 is the vacuum permittivity. The reduced mass M red is given by where M Sr , M Ti , and M O are the atomic masses of strontium, titanium, and oxygen, respectively. The background dielectric constant ϵ ∞ ¼ 5 stands for the linear electronic polarizability of the STO crystal [49].
When neglecting for the nonce any electronic local-field effects, ω ions is the only parameter determining the softmode oscillator strength. Now, if we assume that the soft mode is essentially a relative motion between the positive (Sr 2þ and Ti 4þ ) and negative (O 2− ) ions [see Fig. 3(a)] with the reduced mass M red , as in Ref. [20], the large value of ϵ soft ðωÞ ¼ 4000 leads to a high value of Z eff ≈ 13e 0 (with e 0 as the elementary charge). However, this value is unphysical to be considered as a real charge on the ions. The high value of Z eff can be explained immediately by introducing the local-field corrections according to the Clausius-Mossotti relation [14,15]: where the total polarizability is defined as the sum of electronic, N el α elec ðωÞ, and ionic, N ions α ions ðωÞ, contributions. Within such a concept, the polarizability due to purely ionic motions shown in Fig. 3(c) is small and corresponds to realistic charges on the ions. Thus, the large value of Z eff from the linear response [8,13] results from the hybrid electronic-ionic character of the soft mode. The dipole-dipole coupling via local fields controls the admixture of electronic degrees of freedom to the soft modes and strongly renormalizes the soft-mode frequency from 20 THz down to 1 THz. Concomitantly, electronic interband currents gradually resume the dominance in the contribution to the macroscopic dielectric function. Based on Cochran's equations [1], we have calculated the portions of electronic vs ionic currents to the soft mode via the integral of the real part of the conductivity, as detailed in Sec. IV of the Supplemental Material [36]. A comparison in Fig. 3(d) clearly shows that, in the present situation, the electronic part is at least 5 times larger than the ionic part. This clearly points to a predominant contribution of electronic degrees of freedom in the soft mode, thereby answering the first open question.
To scrutinize the hybrid character of the soft mode and the predominant electronic contribution further, we discuss the observation of nonlinear signals from our 2D THz experiments in the following.

B. Overview of 2D nonlinear signals
The 2D THz spectroscopy with two phase-locked THz pulses, separated by the delay time τ, is used to study the nonlinear response of the soft mode. The total electric field transmitted through the STO film when both THz pulses are present is shown in Fig. 4(a). Figure 4(b) shows the emitted nonlinear field E NL ðt; τÞ. By performing a 2D Fourier transform of E NL ðt; τÞ, we obtain the nonlinear signals in the frequency domain as a function of the detection frequency ν t and the excitation frequency ν τ . The 2D spectrum in Fig. 4(c) is comprised of four types of nonlinear signals. These signals are also known as the χ ð3Þnonlinear signals [29,39] since they result from three-field interactions, with at least one field from each THz pulse.
All of these occur at the detection frequency ν t ¼ ν 0 ¼ 1.17 THz, corresponding to the soft-mode resonance. Within the 2D frequency map, all nonlinear signals can be expressed as a linear combination of frequency vectors [see green (ν A ) and red (ν B ) arrows in Fig. 4(d)] of the incident THz pulses [29,39], which ultimately defines the signal. The green and red arrows have one-to-one correspondence to the wave vectors k A and k B in the wave-vector space used in noncollinear 2D spectroscopy [29].
The most intense nonlinear signals in the 2D map are the so-called pump-probe signals, which are further classified as two types: (i) A pu -B pr (E AB pp ) located at ðν t ν τ Þ ¼ ðν  that they carry the information from the phase evolution of the probe fields only, while the phase evolution of the respective pump fields cancels out. For example, the signal at (ν 0 ; 0) is expressed as The so-called photon-echo signals appear less pronounced in the map and are also classified as two types: (i) the ABB photon-echo signal (E ABB pe ) located at ðν t ; ν τ Þ ¼ ðν 0 ; ν 0 Þ and (ii) the BAA photon-echo signal (E BAA pe ) located at ðν t ; ν τ Þ ¼ ðν 0 ; −2ν 0 Þ. These signals are called photon-echo signals because they contain frequency vector combinations that preserve the phase evolution from both fields, as in conventional photon-echo experiments [50,51]. The signal at (ν 0 , ν 0 ) can be expressed as ν ABB ¼ 2ν B − ν A , while the signal at (ν 0 ; −2ν 0 ) can be expressed as Our theoretical approach to model the soft-mode nonlinear response is based on the simplified pseudopotential concept of electronic band structure, elaborated in Secs. V and VI of the Supplemental Material [36]. The light-matter interactions are treated within the dipole approximation by the use of velocity gauge. Because of the coupling between the electronic-and ionic-displacement dipoles, the interaction with incident THz fields results in a nonlinear electronic motion. Such nonlinear motion of electrons further modifies the polarization of the material, eventually leading to the emission of a nonlinear electromagnetic field [see Eq. (16) and Fig. S5(a) of the Supplemental Material [36]]. The 2D Fourier transform of this nonlinear field shows a striking agreement [ Fig. 4(d)] with our experiments. We now individually scrutinize all the 2D nonlinear signals in Fig. 4(c) to illustrate their role within our hybrid picture.

C. Pump-probe signals
The separation of the different nonlinear signals in Fig. 4(c) allows us to examine their contributions to soft-mode nonlinearities in a background-free manner. We do this by applying a 2D-Gaussian spectral filter to the nonlinear signals in Fig. 4(c). The filtered signals corresponding to A pu -B pr and B pu -A pr are shown in Figs. 5(a) and 5(b), respectively. We perform a Fourierback-transform of these signals from the frequency (ν t , ν τ ) to the time (t, τ) domain. The experimental spectra of E AB pp ðt; τÞ and E BA pp ðt; τÞ signals in Figs. 5(c) and 5(d) and the corresponding theoretical spectra in Figs. 5(e) and 5(f), respectively, provide not only a qualitative but a rather quantitative agreement.
The blueshift of the soft-mode frequency is one of the most striking manifestations of soft-mode nonlinearities. To examine this further, we perform 1D Fourier transforms  Fig. 6(b) shows, the shift reaches a value as large as 0.3 THz, or 26% of the soft-mode frequency. Here, the coupling between the electronic-and ionic-displacement dipoles via the Lorentz field gives rise to local-field effects that strongly modify the character of the nonlinear response. It is quite remarkable that even a THz electric-field strength as low as 49 kV=cm, as used in our work, is sufficient to drive the system far out of the perturbative regime of light-matter interaction [9,29,39].
To characterize the spectral evolution of the soft-mode frequency over the entire delay range, we evaluate the spectrally resolved A pu -B pr signal, S AB pp ðν t ; τÞ. This signal, shown in Fig. 6(c), is obtained by performing a 1D Fourier transform of E AB pp ðt; τÞ along t and using the relation S AB pp ðν t ; Slices of S AB pp ðν t ; τÞ at several values of τ are shown in Fig. 6(d). These spectra clearly display a transmission decrease [ΔT < 0, blue shaded region in Fig. 6(d)] in the range of the original 1.17-THz soft-mode frequency and a concomitant transmission increase [ΔT > 0, red shaded region in Fig. 6(d)] at higher frequencies. Since the material is so close to the ferroelectric phase transition, the soft mode is extremely susceptible to small perturbations that can easily induce nonlinear motions. As a consequence, once the pump pulse leaves the sample, the soft-mode excitation changes its frequency during its decay [9,39], resulting in the observed behavior. Note that this is in striking contrast to the purely phononic nonlinearities [29,39], where one strictly expects a bleaching (ΔT > 0) at the original soft-mode frequency and an induced absorption (ΔT < 0) at the frequency of the excited soft mode. The spectrally resolved A pu -B pr signal along the entire delay range in Fig. 6(c) further exhibits the long-lived (Δτ > 3 ps) character of the soft-mode response, with respect to the duration of the THz pump pulse. The complementary pump-probe signal S BA pp ðν t ; τÞ with the inverted pulse sequence shows the same behavior (see Fig. S7 in the Supplemental Material [36]), in agreement with the interrelation of the pump-probe experiments for χ ð3Þ nonlinearities [39].

D. Photon-echo signals
The 2D-filtered ABB and BAA photon-echo signals are shown in Figs. 7(a) and 7(b), respectively. As before, we perform a Fourier-back-transformation of these signals to get the corresponding nonlinear signals in the time (t, τ) domain. The interaction sequence of the pulses dictates that the photon-echo signals have defined phase fronts. The A-B-B pulse interaction sequence gives phase fronts perpendicular to the green dashed line in Fig. 7, while the B-A-A pulse sequence gives phase fronts that are almost parallel to the green dashed line in Fig. 7. The contour plots of E ABB pe ðt; τÞ and E BAA pe ðt; τÞ in Figs. 7(c) and 7(d) and the corresponding theoretical 2D plots in Figs. 7(e) and 7(f), respectively, display a beautiful agreement. Our observation of photon-echo signals with a significant strength at negative coherence times τ x is striking, and yet it coheres with the expected soft-mode nonlinearities driven by electronic motions; i.e., the softness of the system near the phase transition allows for a nontrivial field overlap.
In our hybrid picture, the local-field effects lead to a situation where the effective electric field resulting from the interaction of all the dipoles (i.e., electronic and ionic) is much larger than the average macroscopic THz electric field in the crystal. This leads to a significant overlap of the THz pulses A and B, resulting in the photon-echo signals at negative coherence times [9,39,[51][52][53][54]. In addition, it characterizes the long-lived nature (>3 ps) of the driving fields, as expected for noninstantaneous contributions that are dominated by local-field effects. In contrast, for the purely phononic picture, theoretical 2D simulations show that both the nonlinear ABB and BAA photon-echo signals are only present for positive coherence times; see Figs. S9 (b) and S9(c) in Sec. IX of the Supplemental Material [36]. Evidently, the presence of photon-echo signals at negative coherence times is the very hallmark of the low-field softmode nonlinearity, which is direct experimental evidence for the relevant local-field contributions to the soft-mode nonlinearity-the answer to the third open question. Furthermore, our hybrid model has reproduced all experimental results (i.e., both the linear and nonlinear responses) in a quantitative way, thereby answering the fourth open question.

V. CONCLUSION
In conclusion, we have presented a comprehensive study on the origin of THz nonlinearities of a polar soft mode in a strain-engineered epitaxial ferroelectric thin film. We find that the soft mode is a hybrid mode of pure lattice (ionic) motions and electronic interband transitions. The dielectric function of the soft mode in our ferroelectric film shows a predominant electronic contribution, which is at least 5 times larger than the purely ionic contribution. The apparent low threshold for the soft-mode nonlinear response results from a combination of two complementary facts: first, the soft mode has a hybrid nature, and second, being close to the ferroelectric phase transition, we can easily drive electronic motions into the nonlinear regime with low THz fields unlike pure ionic motions. Our expanded hybrid picture quantitatively reproduces all the features associated with the soft-mode nonlinearities, namely, the blueshift of the soft-mode frequency and the photon-echo signals at negative coherence times, which are, however, in striking contrast to purely phononic nonlinearities. In our model, the local-field effects that originate from the coupling of electronic and ionic degrees of freedom not only dominate but also strongly modify the soft-mode nonlinearities, while switching them off leads to a vanishing nonlinear signal. We identified several important open questions and provided quantitative answers that shed light on the microscopic picture behind the THz softmode nonlinearities in ferroelectric perovskites.
At the fundamental level, these findings can be extrapolated to any polar mode that is intrinsically linked to the softness of a system near a phase transition. Going beyond ferroelectrics, analogies can be drawn for unstable spin excitations (for example, certain antiferromagnetic magnons [55,56]) and thereby associate their role during light-driven phase transitions in magnetic materials. From a more utilitarian perspective, 2D THz spectroscopy turns out to be a very powerful tool to unravel the interplay between the transient nonlinear electronic motions and the macroscopic polarization, an understanding of which is beneficial for engineering oxidebased heterostructures [57] with functional properties. On the other hand, we can make use of our understanding to drive polar phonons nonlinearly by optical straining and potentially steer many-body quantum phenomena, such as superconductivity. Strain, in particular, has been proved to be an efficient way to control the soft-modeinduced superconductivity [58]. Observations such as an increase of superconducting T c under tensile strain, broadening of the superconducting dome, and a shift of T c towards lower carrier densities are general for any superconductor where pairing is mediated by softening of ferroelectric modes [58][59][60]. While these observations are mostly theoretical, our work provides a strong experimental groundwork. Furthermore, our results have deeper implications for classical molecular-dynamics simulations that tacitly rely on the purely ionic picture. The molecular-dynamics simulations have been widely used for understanding the nonlinear dynamics of coupled phonon modes across thermal phase transitions. However, they fail to reproduce the correct electronic nonlinearities associated with the softening of phonon modes in a light-induced phase transition at ultrafast timescales-a regime where our hybrid model shows great potential.