{'fraction_non_alphanumeric': 0.029302736098852603, 'fraction_numerical': 0.014916151809355693, 'mean_word_length': 4.763224821973551, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlocal nonlinear phononics

Nonlinear phononics relies on the resonant optical excitation of infrared-active lattice vibrations to coherently induce targeted structural deformations in solids. This form of dynamical crystal-structure design has been applied to control the functional properties of many interesting systems, including magneto-resistive manganites, magnetic materials, superconductors, and ferroelectrics. However, phononics has so far been restricted to protocols in which structural deformations occur locally within the optically excited volume, sometimes resulting in unwanted heating. Here, we extend nonlinear phononics to propagating polaritons, effectively separating in space the optical drive from the functional response. Mid-infrared optical pulses are used to resonantly drive an 18 THz phonon at the surface of ferroelectric LiNbO3. A time-resolved stimulated Raman scattering probe reveals that the ferroelectric polarization is reduced over the entire 50 micron depth of the sample, far beyond the ~ micron depth of the evanescent phonon field. We attribute the bulk response of the ferroelectric polarization to the excitation of a propagating 2.5 THz soft-mode phonon-polariton. For the highest excitation amplitudes, we reach a regime in which the polarization is reversed. In this this non-perturbative regime, we expect that the polariton model evolves into that of a solitonic domain wall that propagates from the surface into the materials at near the speed of light.

Generally, optical radiation tuned to frequencies immediately above the TO frequency of an infraredactive vibration does not propagate into a material and is evanescently screened by the induced polarization. The interaction of the light field with the material, and the ability to manipulate its microscopic properties remains therefore inherently local. Here, we show that if the material breaks inversion symmetry, the optically excited mode can couple anharmonically to other symmetry-odd, infrared active excitations of the solid, which hybridize with the electromagnetic field and form a propagating polariton. This type of nonlinear phonon-polariton coupling between an evanescent drive and a propagating mode opens up a broad class of phenomena that we term nonlocal nonlinear phononics.
Nonlocal polariton excitation builds on the extensive work performed in the past in which polaritons were excited using impulsive stimulated Raman scattering from propagating near-infrared pulses 11,12 In these studies, steering, focusing, and amplification of THz radiation was shown to occur in regions of the crystal beyond those traversed by the excitation pulse 13,14 . However, even in cases in which phase matching was optimized to drive phonon polaritons of the largest amplitudes, these experiments never reached the deep nonlinear regime that is afforded by nonlinear phononics. In the present paper we show that nonlinear phononics can bring this class of phenomena to the more consequential regime of ferroelectric switching, predicted 15 but never observed experimentally through polaritons.
The principle of nonlocal nonlinear phononics is sketched in Fig. 1a for the model ferroelectric LiNbO3, which is studied in this work. Both infrared active mode QIR and polariton QP are connected to lattice vibrations. The c-axis polarized 18.3 THz A1 mode (QIR) is nonlinearly coupled to the 7.3 THz A1 "soft mode" (QP), which in turn modulates the ferroelectric polarization. Previous work has shown that driving the high-frequency vibration QIR with a resonant sub-picosecond mid-infrared pulse results in a transient reversal of the ferroelectric polarization. However, the reversal was only detected below the surface of the material by time resolved second-harmonic generation 16 . The spatial extent of this effect was not accessible by that method, as it was limited by the micrometer coherence length of the detection process. Here, we demonstrate that the ferroelectric polarization in LiNbO3 can be controlled non-locally throughout the bulk. Resonant excitation of QIR at 18.3 THz also launches an anharmonically coupled low-frequency polariton of the soft mode QP, which, contrary to the evanescent driven mode, propagates throughout the crystal (Figure 1b). Femtosecond stimulated Raman scattering measurements, which employ short pulses of visible light below the optical bandgap, interrogate the lattice of LiNbO3 and reveal a reduction of the Raman cross section for a number of polar modes, which disappear transiently as the ferroelectric polarization vanishes. At the highest fields, we even observe that the Raman susceptibility of these modes changes sign, indicating a reversal of the ferroelectric polarization throughout the bulk of the solid, and in fact far beyond the evanescent field of the mid-infrared optical drive.
To describe these physics quantitatively, we consider the following potential for the dynamics of the crystal lattice: are high-order nonlinearities along the phonon coordinates and , relevant during high amplitude excitation 17 . The last term describes the lowest-order anharmonic coupling between and which enables energy transfer between the two modes. The polarization in the crystal is given by: where is the vacuum permittivity and ∞ is the high frequency dielectric constant. = * + * is the ionic contribution to the polarization resulting from the effective charges * and * of the modes and , respectively. = * + * is the polarization which results from the coupling of and to the total electric field E via the respective Raman tensors * and * . Assuming a light-matter interaction given by • , the respective equations of motion are: An incident electric field thus imparts two forces upon and : one that can couple directly to the modes through their effective charges and a "Raman" force, which we found to be much weaker by  THz) mode of the polariton branch shown in Figure 1(b), resulting in a distinct polarization wave propagating throughout the material. As discussed below, we assert that the control of the ferroelectric polarization in the bulk of LiNbO3 is intertwined with the propagation of this ferroelectric soft mode polariton.
In our experiments, we studied the response of the Raman-active lattice modes throughout a 50 micron thick LiNbO3 single crystal using femtosecond stimulated Raman spectroscopy (FSRS) [18][19][20] . The theory behind FSRS has been described elsewhere 21,22 . Here, we briefly summarize the general concept.  THz). To probe the response of the A1 modes, both and were polarized along the c-axis, and the time delay between these two pulses was kept fixed. An FSRS spectrum was recorded for each pumpprobe delay between the FSRS pulse pair and the mid-infrared pulse as depicted in Figure 3(a).
Further details on the experimental setup are provided in the supplementary information 23 . THz. This is effectively a form of second-order Raman scattering which results from the excitation of atomic motions along the coordinates just below the TO frequency . The group velocity at the 400 nm wavelength of the FSRS probe matches the phase velocity of the polariton at 16 THz 23 , thus allowing for the effective detection of effects on the Raman tensor at this modulation frequency.
More important is the reduction in the peak amplitudes of and . This effect was well-modelled under the assumption that the equilibrium components of both * and * become transiently reduced upon excitation with the mid-IR pulse. As shown in Figure 4(b), the group velocity of the probe around 400 nm wavelengths approximately matches the soft mode polariton phase velocity over the course of the 50 µm sample volume. Thus, reductions in the Raman tensor correlated to the propagating soft mode polariton can be detected.
To produce the calculated pumped spectra shown in Figure 4 The final curve, shown in black in Figure 4(c), models the effective * ( ) seen by the probe and was used to produce the spectra shown in Figure 4(a). Complete details regarding the calculations are provided in the supplementary information 23 .
The fluence dependence of the peak reductions are shown in Figure 5(a) for a fixed time delay  = 0.08 ps. As the mid-infrared excitation influence is increased, the FSRS peaks decrease monotonically in amplitude and then reverse sign at the maximum fluence of 185 mJ/cm 2 used in the experiment. First principles calculations presented in Figure 5(b) show that if the ferroelectric polarization is reversed, the signs of both Raman tensors * and * would also reverse. The amplitudes of the FSRS peaks corresponding to and as a function of mid-IR pump fluence are shown in Figure 5(c). For fluences above ~ 60 mJ/cm 2 , the peak integrals scale nonlinearly with fluence and display a saturation-like behaviour before they reverse sign. This is likely due to a softening of the frequency of the driven mode QIR, an effect which occurs when a vibration is driven to high amplitudes and effectively reduces its response to the driving field 17  It should be stated explicitly here that the manipulation of the ferroelectric polarization in the bulk is not described comprehensively by the theory of propagating polaritons which is illustrated in Figure   1. When a ferroelectric domain is either reduced or reversed at the surface, a solitonic phase front of altered domains will also presumably propagate into the bulk at near-luminal speeds. Such a phenomenon breaches the highly interesting regime of driven phase transitions that should be explored further, both theoretically and experimentally.
In summary, by measuring the dynamical Raman tensor elements in LiNbO3 using a stimulated Raman scattering probe, we show control of the ferroelectric polarization well-beyond the excitation region.
We extend the ideas of nonlinear phononics to propagating polaritons, which provides a foundation for explaining this phenomenon. Our experiments also provide evidence for a transient reversal of the ferroelectric polarization throughout the bulk of the material, a striking observation which requires further theoretical pursuit.
We foresee that further steps forward could involve integrating nonlinear phononics with nonlocal optical platforms such as metamaterials. This would introduce the ability to refract 24,25 , tightly focus 26 , and even accelerate 27,28 nonlinearly driven modes. Furthermore, nonlocal nonlinear phononics broadens the scope of optical control in quantum materials and is naturally applicable to other types of propagating polaritons 29 . Finally, we note that the physics discussed here is likely connected to propagating phase changes induced by phononics across hetero-interfaces, [30][31][32] which to date have remained unexplained.   As discussed in the text, the features in the pumped spectrum can be explained under the assumption that the Raman tensors for the A1 modes change dynamically.

(S2)
Here, n is the index of refraction, c is the vacuum speed of light, and is the vacuum permittivity.
Eqns. (S1) -(S2) describe a four-wave mixing process in which is ultimately emitted in the same direction as [34]. When sent to a spectrometer, the two fields interfere in the frequency domain, with the resulting intensity given by: As shown in Figure 2 of the manuscript, the frequency components of at − and + interfere with the incident field constructively and destructively, respectively. All FSRS spectra shown in the manuscript are normalized to the power spectrum of according to:

S.2) Sample preparation
All measurements were performed at room temperature and in ambient air. The sample investigated was a 50µm x-cut LiNbO3 single crystal which was coated on the back side with an anti-reflection layer designed for 400 nm wavelengths. This reduced the prominence of fringes in the spectrum of the white light continuum due to Fabry-Perot interferences.

S.3) Experimental Setup a) mid-IR pump.
A regeneratively amplified Ti:sapphire oscillator operating at a 1 kHz repetition rate was used for the experiments, delivering pulses of ~100 fs duration at 800 nm carrier wavelength. 4 mJ of the pulse energy was used to pump an optical parametric amplifier which produced 1.2 mJ of combined signal and idler when tuned to a signal wavelength of 1520 nm. The two near-infrared beams were then used to generate mid-infrared pulses with 12 µJ of energy and 20 THz central frequency (7 THz FWHM) by difference frequency mixing in a GaSe crystal.

b) FSRS probe.
To construct the FSRS probe, a portion of the remaining 800 nm output from the regenerative amplifier was split into two paths. The first path employed roughly 1µJ in order to generate a white light continuum pulse in a CaF2 substrate. Using a motorized stage, the substrate was kept under constant translational motion in order to avoid damage. The spectral weight of the continuum pulse extended to 350 nm wavelengths; however, the edge pass filter used to remove the residual 800 nm light reduced the shortest usable wavelength to 388 nm. The continuum pulse was collimated and then focused onto the sample using reflective optics.
To generate the narrowband Raman pump field, 10 µJ pulses of the 800 nm was incident upon a 1.5 cm thick type-I BBO, generating 1.4 µJ of second harmonic with central wavelength of 400 nm (750 THz). The pulse envelope was characterized by cross correlation with the white light continuum pulse in a time-resolved two-photon absorption measurement in the 50 µm-thick LiNbO3. The temporal width of the pulse envelope, which is shown in Figure S1, is determined by the group velocity mismatch between the fundamental 800 nm pulse and the generated 400 nm pulse in the BBO. The asymmetry in the intensity envelope results from the depletion of the 800 nm pulse as it passes through the crystal. The temporal width between the rising and falling edge of the pulse was roughly 3 ps.

c) FSRS data acquisition
In order to facilitate single-shot detection of the white light continuum pulses, the CCD of the visible spectrometer was triggered at the 1 kHz repetition rate of the regenerative amplifier. The narrowband Raman pump pulses were mechanically chopped at half the frequency (500 Hz), facilitating a consecutive narrowband pulse on/off sequence impinging upon the sample. To produce a given FSRS spectrum, we collected a series of shots (minimum 50,000) and binned the data such that the sum of the "on" shots was divided by the sum of the "off" shots. The delay between the white light continuum and narrowband pulses was fixed to 1.85 ps with respect to the time trace shown in Figure S1.
When performing time-resolved FSRS measurements, the mid-IR pump illuminated the sample during both "on" and "off" shots. The sum of the "on" shots was divided by the sum of the "off" shots, in the same manner as the equilibrium FSRS measurements. Thus, effects which are dependent on the carrierenvelope phase of the mid-IR pulse, which can vary over modulo 2π from shot to shot, are averaged out in both the "on" and "off" shots. This was essential for averaging out parasitic carrier-envelope phase dependent effects of the impinging mid-IR pulses on the white light continuum spectrum.

S.4) Consideration of the polariton and probe velocities
The velocity of the probe pulses relative to the excitations arising from the mid-IR pump was determined by considering the group refractive index of LiNbO3 at 400 nm ( 0 = 2.86) [35] as well as the = 30° angle between the mid-IR pump and the white light continuum. In the rest frame of the polariton, the probe moves slower than it would in the case that it was collinear with the polariton; the effective group velocity is given by = / = c/( 0 cos (30°)) where is the speed of light in vacuum. Figure S2 shows the mapping of 0 and onto the dispersion relation for the c-axis modes, which was generated by using parameters from infrared reflectivity measurements in Ref. 36. Notably, the group velocity line intersects the polariton branch of the 18.3 THz TO phonon at a frequency of 16 THz. This indicates that the probe is phase-matched to effects which result from the propagating 16 THz polariton.
The phase index at 2.5 THz is = 5.15, yielding a ratio of / = 1.56. The corresponding velocity line was then used for the plot of Figure 4(b) and in the calculations which account for the velocity mismatch between the probe and the propagating reduction of the Raman tensor (see section S.5).

S.5) Time-domain model of the FSRS process
Here, we discuss the model that was used to produce the calculations shown in Figures 4(a)  (S10) To model the equilibrium FSRS spectrum, the parameters * and were chosen so that eq. S10 fits the experimental data.
Our model for the dynamic FSRS spectra resulting from mid-IR excitation assumes that * becomes time-dependent and consists of two components: * ( ) = − ( ) + 16 ( ) (S11) Here, − ( ) corresponds to a transient reduction of the equilibrium Raman tensor for the n th c-axis mode, and 16 ( ) models the oscillatory component responsible for the second-order scattering. The shape of − ( ) was determined by first assuming that the ferroelectric polarization P is reduced over the same time scale as the data reported in Ref. [37]. First principles calculations of the Raman tensors * ( ) and * ( ) (see Figure 5(b) of the main text) were then mapped onto the time-dependent polarization in order to produce * ( ) and * ( ). Our model then assumes that these Raman tensors propagate throughout a 50 µm sample volume with a group velocity equal to the soft mode polariton phase velocity. To account for the velocity mismatch between the propagating * ( ) and the FSRS probe, the velocity line of the probe was projected onto the space-time maps for the propagating Raman tensors, and an integral along the velocity line was carried out for each time step. This produces the effective − ( ) and − ( ) which are seen in the rest frame of the probe. For each data set, the degree of polarization reduction was adjusted in order to produce a given − ( ) and − ( ) which produced calculations that best fit the data. The Raman tensors of the two additional c-axis modes which appear in the Raman spectrum at 8.3 THz and 10 THz were assumed follow the same dynamics as − ( ). 16 ( ) was assumed to have the form: where [ ] is an error function with rise parameter a that was chosen by integrating the Gaussian envelope of the mid-IR pulse used in the experiments. Additionally, is the lifetime of the 18.3 THz mode extracted from infrared measurements reported in Ref. 36, f was taken to be 16 THz, A was chosen as a fit for each data set, and is the phase of the oscillation. Generally, the phase of a linearly driven polariton is dependent on the carrier envelope phase (CEP) of the excitation pulse. In the experiments, the CEP of the mid-IR source fluctuates randomly from shot to shot, hence the timeresolved FSRS measurements were averaged over ~ 100,000 shots. To account for the phase averaging, equations S1 -S5 were evaluated from = 0 2 in steps of π/2 for each delay between the FSRS probe and * ( + ). The expression which describes the simulated CEP-averaged FSRS spectrum for a given is then:

S.6) DFT calculations of the polarization dependent Raman tensor
The dependence of the Raman tensor elements on the polarization was mapped out by utilizing the first-principle framework of density functional theory (DFT) calculations. We performed our computations with the Vienna ab-initio simulation package VASP.6.2 [38]. Moreover, for the phonon calculation, we used the Phonopy software package [39]. Our computations further utilized pseudopotentials generated within the Projected Augmented Wave (PAW) [40] method. Specifically, we take the following configurated default potentials: Li 2s 1 , Nb 4p 6 5s 1 4d 4 , and O 2s 2 2p 4 . We applied the PBEsol [41] approximation for the exchange-correlation potential. As a numerical setting, we used a 9x9x9 Monkhorst [42] generated k-point-mesh sampling of the Brillouin zone and a plane-wave energy cutoff of 750 eV. The self-consistent calculations were reiterated until the change in total energy becomes less than 10 -9 eV. Last, to compute the Raman tensor elements, we followed the approach of Ref. [43] and determined the R* elements from the induced change in the permittivity due to phonon displacement.
The starting point of our investigation was the ferroelectric ground state of LiNbO3. To create a forcefree reference configuration, we first structurally relaxed the LNO unit-cell in DFT. We initialized the computation with the experimental structure with symmetry R3c of Ref. [44]. We then minimized the system until the Hellman Feynman forces on the atoms and the pressure on the unit cell became less than 0.1 eV/Ang. Tab. S1 contains the final structure. For this atomic arrangement, we then computed the phonon spectrum at the zone center and identified the four polar A1 modes at the frequencies noted in the lower part of Tab S1. All results agree with experimental findings.
Next, we constructed the reversed ferroelectric (FE) polarization configuration by inverting the positions of all atoms. Further, we built the centrosymmetric state, which lies in between the two FE states. Utilizing these three states, we finally computed the atomic displacement field, which connects all three states. In Figure S3, we show the evolution of the FE polarization, calculated by the Berry-Phase Approach [45,46] along this displacement field. Note that λ=0 corresponds to the centrosymmetric state.
Finally, we applied the approach of Ref. [43] to compute the Raman tensor for the structures along this trajectory. Note that we recalculated the phonon spectrum for each structural configuration to obtain the corresponding correct eigendisplacements of the phonons. Fig. S3 (b) shows the resulting R* dependence. An apparent reduction and sign change become evident upon passing through the centrosymmetric state, with exception to the 10.1 THz mode, which remains finite through the transition. When restoring inversion symmetry, the phonon modes recapture their unique character to be either infrared or Raman active. In the centrosymmetric state, LNO exhibits three polar modes and one Raman active mode, which become the c-axis polar modes in the ferroelectric phase. The 10.1 THz mode corresponds to the original Raman mode, which maintains its Raman activity, and thus a nonzero R*, in the centrosymmetric state.     Table S1. Structural parameters in the hexagonal setting and phonon frequencies from DFT (this work) and experimental results.