Cavity control of nonlinear phononics

Nonlinear interactions between phonon modes govern the behavior of vibrationally highly excited solids and molecules. Here, we demonstrate theoretically that optical cavities can be used to control the redistribution of energy from a highly excited coherent infrared-active phonon state into the other vibrational degrees of freedom of the system. The hybridization of the infrared-active phonon mode with the fundamental mode of the cavity induces a polaritonic splitting that we use to tune the nonlinear interactions with other vibrational modes in and out of resonance. We show that not only can the efficiency of the redistribution of energy be enhanced or decreased, but also the underlying scattering mechanisms may be changed. This work introduces the concept of cavity control to the field of nonlinear phononics, enabling nonequilibrium quantum optical engineering of new states of matter.


I. INTRODUCTION
Phonon modes describe the vibrational behavior of atoms that experience a harmonic increase of the potential energy when they move away from their equilibrium positions in a crystal lattice. When the dipole moment of a phonon mode is resonantly driven with the electric field component of an ultrashort laser pulse, the amplitudes of the vibrations become so large that anharmonic corrections to the potential energy are essential. The vibrational dynamics of the solid are then governed by nonlinear interactions between the phonon modes which allow for the electronic correlations in the material to be modified, leading to emergent states of matter that do not exist in equilibrium [1,2]. Two of the most prominent examples of nonlinear phonon induced states include superconductivity far above the equilibrium critical temperature [3][4][5][6][7][8] and light-induced ferroelectricity and ferroelectric switching [9][10][11][12][13][14]. These mechanisms are based on the energy redistribution from an initially coherently excited infrared (IR)-active phonon mode to a Ramanactive one that couples to the electronic system through transient modulations of the crystal structure, exchange interaction, or the electron-phonon coupling. The efficiency of the energy transfer is however often hindered due to a frequency mismatch of the phonon modes participating in the interaction, as illustrated in Fig. 1(b), limiting the exploitation of nonlinear phonon induced states to a few materials in which these resonances occur naturally.
At the same time, developments in the field of strong light-matter coupling have made it possible to engineer the optical properties of materials in optical cavities. An electromagnetic mode sustained in the cavity hybridizes * djuraschek@seas.harvard.edu † tomasneuman@g.harvard.edu ‡ prineha@seas.harvard.edu with optically active states of the material to form polariton states as shown schematically in Fig. 1(a). This polaritonic engineering has been used to control the chemical reactivity of molecules in their ground and excited states through a modulation of the chemical landscape [15][16][17][18][19][20][21][22]. Recently, various studies for solids have suggested an influence of strong light-matter coupling on superconductivity by exciton-photon or phonon-photon coupling [23][24][25][26][27][28]. Despite the tremendous progress in theoretical methods for strong coupling in molecular and solid-state systems, the influence of strong light-matter coupling on nonlinear phonon-phonon interactions has Cavity-mediated resonance in nonlinear phononics. (a) Schematic setup of a sample material in a terahertz cavity. The energy level diagram shows the formation of two phonon-polariton modes with frequencies Ω+ and Ω− from the hybridization between an IR-active phonon mode, ΩIR, and the fundamental cavity mode, ω1. (b) Frequency mismatch between an IR-active and two Raman-active phonon modes, ΩIR and ΩR, hinders the transfer of energy between them. (c) Polaritonic splitting matches the energy of one of the branches, Ω−, with the Raman-active phonon modes. remained unexplored, motivating this work. Strong lightmatter coupling, like nonlinear phonon interactions, relies on the matching of energies to phonon resonances in solids.
We therefore synthesize the physics of nonlinear phononics and strong light-matter coupling in this work to overcome the frequency mismatch in nonlinear phononic processes. Specifically, we demonstrate that we can tune the polaritonic splitting of an IR-active phonon mode in the cavity such that one of the polariton branches is pushed into resonance with nonlinearly coupled Raman-active phonon modes. As a result, we are able to greatly enhance the efficiency of energy redistribution through the nonlinear interaction and even change the underlying process entirely. We sketch the principle of our technique for the example of a parametric downconversion process in Figs. 1(b) and (c). In Fig. 1(b), an IR-active phonon mode with frequency Ω IR is out of resonance with two Raman-active phonon modes with frequencies Ω R , preventing efficient coupling between them. By splitting the IR-active phonon mode in a terahertz-frequency cavity in Fig. 1(c), one of the polariton branches is shifted into resonance with the Ramanactive phonon modes, Ω − = 2Ω R , enabling the redistribution of energy between them.

II. NONLINEAR PHONONIC PROCESSES
We begin by analyzing the dominating nonlinear phononic mechanisms of energy redistribution from a coherently excited IR-active phonon mode to a Ramanactive one, of which we provide an overview in Table I. The potential energy V including the most prominent phonon couplings beyond the harmonic approximation can be written as [29] V (Q) = Ω 2 Here, Ω IR and Ω R are the eigenfrequencies of the IR and Raman-active phonon modes in 2πTHz, and Q IR and Q R are their respective amplitudes inÅ √ amu, where amu is the atomic mass unit. c, d, and e are the phonon coupling coefficients in meV/(Å √ amu) n , where n is the order of expansion in the phonon amplitudes [30].
Nonlinear phononic mechanisms can be classified by those linear and squared in the Raman-active phonon amplitude. In the linear case, the square of the IRactive phonon amplitude Q 2 IR acts as an effective driving force on the Raman-active phonon mode. If Q IR ∼ A IR sin(Ω IR t) follows a sinusoidal shape with a decaying amplitude A IR , then Q 2 IR ∼ A 2 IR [1 − 2 cos(2Ω IR t)] possesses one static component and another oscillating at twice the frequency. Depending on whether the IRactive phonon mode has a higher or smaller frequency than the Raman-active one, this coupling leads to either difference-frequency (DF) or sum-frequency ionic Raman scattering (SF-IRS), see Table I. Difference-frequency ionic Raman scattering is a nonresonant process that creates quasi-static transient distortions in the crystal structure due to its static force component on the Ramanactive phonon mode, and it is held responsible for the enhancement of superconductivity in yttrium barium copper oxide [3,8,31,32] and for ferroelectric switching in lithium niobate [11,33]. Sum-frequency ionic Raman scattering, in contrast, is a two-phonon absorption process resonant at Ω IR = Ω R /2, which has been predicted to induce Raman-active phonon mode oscillations with large amplitudes [30,34,35]. In the quadratic case, the IR-active phonon amplitude acts as an effective modulation of the Raman-active phonon frequency, This coupling leads to a parametric amplification of the Raman-active phonon mode in different orders, and it is resonant when either Table I. Parametric amplification could be responsible for the enhancement of superconductivity through an enhancement of electron-phonon coupling [36][37][38][39], and is a well-established method to generate squeezed phonon states [40][41][42][43]. The symmetry of the parametric coupling allows an IR-active phonon mode, in addition to Raman-active phonon modes, to also couple to silent phonon modes that are otherwise inaccessible for coherent excitation [44]. In the remainder of the manuscript, we will focus on the leading-order contributions from the three-phonon processes in regimes I and II and set e = 0.

III. THEORETICAL FORMALISM
In order to evaluate the interaction of the phonon modes with the cavity, we consider a quasi-1D model, in which a slab sample is placed at position z 0 into a cavity formed by two parabolic mirrors at positions z = 0 and L, see Fig. 1(b). As only Brillouin-zone center phonon modes are relevant to the dynamics, we consider translational invariance in directions parallel to the cavity mirrors. These assumptions guarantee the transversality of the fields that enter the dynamics of the system, which simplifies their theoretical treatment.
We begin by considering Ampere's law in the Lorenz gauge expressed via the transverse vector potential A: where µ 0 is the vacuum permeability, c the speed of light, and J is a volume current density generated by the IRactive phonon mode. J is related to the macroscopic polarization P generated by the IR-active phonon mode by where we have assumed the mode effective charge Z IR of the IR-active phonon mode, polarized along the unit vector e p , to be time independent, and where V c is the volume of the unit cell associated with Z IR . (Note that for very strong driving, the mode effective charge can be dependent on the phonon amplitude [45].) Assuming perfectly reflecting mirrors, the vector potential of the cavity mode m of polarization p can be expressed as: where k m = mπ/L, m is a natural number, and A m,p are time-dependent, but spatially invariant amplitudes. The two transverse polarizations p = {x, y} of the cavity mode are degenerate. In the following, we pick the relevant polarization direction aligned with the dipole moment of the IR-active phonon mode and take into account only on the fundamental cavity mode, m = 1, whose frequency ω 1 = πc/L is assumed to be close to Ω IR . We justify this assumption a posteriori, as our results turn out to remain unchanged when we take into account higherorder cavity modes. Integrating over the spatial coordinate z, the equation of motion coupling the IR-active and cavity mode can be written in a scalar form, with A 1,x → A 1 , as where ε 0 is the vacuum permittivity, and ∆z is the thickness of the slab. More details of the derivation and the general form including multiple cavity modes are provided in the Supplementary Material [46]. Defining G 1 ≡ sin(k 1 z 0 )Z IR , the coupled equations of motion for the IR-active, Raman-active, and cavity modes driven by an ultrashort terahertz pulse yield Here is the electric field component of the terahertz pulse, where E 0 is the peak electric field, ω 0 is the center frequency, and τ is the full width at half maximum duration of the pulse [30]. κ i denote phenomenological linewidths of the respective phonon and cavity modes.

IV. NONLINEAR PHONON DYNAMICS
We now investigate how the coherent nonlinear phonon dynamics change when we tune the cavity mode in and out of resonance with the IR-active phonon mode in the geometry shown in Fig. 1(a). As input parameters to Eqs. . In order to explore the efficiency of energy redistribution from the IR-active to the Raman-active phonon mode, we vary the frequencies of the terahertz pulse and the cavity mode, as it would be possible in an experimental setup. We look at the two different regimes, in which the resonant three-phonon processes work: Regime I, Ω IR > Ω R , for difference-frequency ionic Raman scattering and cubic-order parametric amplification, and regime II, Ω IR < Ω R , for sum-frequency ionic Raman scattering.

A. Regime I: ΩIR > ΩR
We demonstrate our technique at the example of a 15 THz IR and a 6.75 THz Raman-active phonon mode that are detuned by ∆Ω IR /(2π) = 1.5 THz such that Ω IR −∆Ω = 2Ω R . We numerically solve Eqs. (7)  phonon modes in response to the excitation by a terahertz pulse with a peak electric field of E 0 = 5 MV/cm and a duration of τ = 250 fs in Fig. 2. In Figs. 2(a) and (b), we show maps of the integrated spectral densities, S i = ∞ 0 |F i | 2 dω that we obtain from Fast Fourier Transformation (FFT) spectra F i of the time-dependent phonon amplitudes of the IR (i=IR) and Raman-active (i=R) phonon modes, respectively. We plot the maps as a function of the pulse frequency ω 0 (vertical axis) and cavity frequency ω 1 (horizontal axis). The values are normalized to the maximum value that we obtain for the IR-active phonon mode, S 0 = max[S IR (ω)]. A clear signature of phonon-photon strong coupling emerges in the spectral map of the IR-active phonon mode in Fig. 2(a). The maxima form two continuous branches that feature an avoided crossing as the cavity frequency ω 1 is varied. When the frequency of the cavity is far detuned from the IR-active phonon frequency, ω 1 Ω IR and ω 1 Ω IR , the integrated spectral density peaks when the terahertz pulse is resonant with the frequency of the IR-active phonon mode, ω 0 = Ω IR . The integrated spectral density reduces as ω 1 is approaching the value of Ω IR towards the center of the map, as the cavity mode hybridizes with the IR-active phonon mode and the energy is shared between the two phonon-polariton branches. In contrast, the integrated spectral density of the Ramanactive phonon mode Fig. 2(b) exhibits a sharp peak at ω 1 /(2π) ≈ 16.8 THz and ω 0 /(2π) ≈ 13.9 THz. Here, the frequency of the lower phonon-polariton branch is shifted into resonance with the Raman-active phonon modes, Ω − = 2Ω R , strongly enhancing the efficiency of energy redistribution. The marker "OFF" denotes a point at which the IR-active phonon mode is driven resonantly by the terahertz pulse, but the energy redistribution to the Raman-active phonon mode is hindered, and the marker "ON" denotes the resonance point of the nonlinear phononic process.
In order to separate the contributions from the mechanisms of difference-frequency ionic Raman scattering and parametric downconversion, we show the FFT spectra F IR and F R of the respective phonon modes at the points marked "OFF" and "ON" in Figs. 2(c) and (d). The values are normalized to the maximum value that we obtain for the IR-active phonon mode, F 0 = max[F IR (ω)]. The FFT spectrum of the IR-active phonon mode in Fig. 2(c) shows a peak at its eigenfrequency of Ω IR /(2π) = 15 THz, when the cavity is out of resonance, and a peak at the frequency of the lower phonon polariton Ω − /(2π) = 13.5 THz, when the cavity is tuned into resonance, thus fulfilling the resonance condition for the cubic-order parametric downconversion process, Ω − = 2Ω R . The FFT spectrum of the Raman-active phonon mode in Fig. 2(d) shows only a static component arising from difference-frequency ionic Amplitude, 3.

Laser frequency
[THz] Cavity frequency [THz] Raman scattering when the cavity is out of resonance. When the cavity is tuned into resonance, this static component decreases, but a large peak at its eigenfrequency of Ω R /(2π) = 6.75 THz emerges.
The cavity-induced resonance becomes even more distinct when looking at the time-dependent evolutions of the phonon amplitudes Q IR and Q R in Figs. 2(e) and (f). The values are normalized to the maximum value that we obtain for the IR-active phonon mode, . The IR-active phonon mode shows the common evolution of a resonantly driven and slowly decaying phonon mode in Fig. 2(e) when the cavity is out of resonance. In resonance, the amplitude exhibits an additional decrease and recovery around 4 ps from the redistribution of energy to the Raman-active phonon mode. Q R in turn exhibits a transient shift of its equilibrium position and only tiny oscillatory amplitude in Fig. 2(f) when the cavity is out of resonance, a common feature of difference-frequency ionic Raman scattering. When the cavity is tuned into resonance, the static shift decreases, but a large-amplitude oscillation builds up over time through the cubic-order parametric downconversion process. The time delay of the build-up depends on various factors, such as the relative frequencies of the phonon modes, their dampings, and the strength of their nonlinear coupling.

B. Regime II: ΩIR < ΩR
We now investigate the case of a 7.5 THz IR and a 16.5 THz Raman-active phonon mode that are detuned by ∆Ω IR /(2π) = 0.75 THz such that 2(Ω IR − ∆Ω IR ) = Ω R . We show the numerical results for an excitation of the phonon modes by a terahertz pulse with a peak electric field of E 0 = 2 MV/cm and a duration of τ = 500 fs in Fig. 3, which is structured in the same way as Fig. 2. Analog to the behavior in regime I, the spectral map of the IR-active phonon mode in Fig. 3(a) shows the formation of phonon-polariton branches as a result of phonon-photon strong coupling. The integrated spectral density of the Raman-active phonon mode Fig. 3(b) exhibits a peak around ω 0 /(2π) ≈ 7.8 THz and ω 1 /(2π) ≈ 5.2 THz, where the frequency of the upper phonon-polariton branch is shifted into resonance with the Raman-active phonon mode, Ω + = Ω R /2, strongly enhancing the efficiency of energy redistribution. The peaks is not as sharp as in the case of parametric downconversion, but the redistribution of energy is so large, that it visibly imprints on the spectral map of the IR-active phonon mode.
We find that the term dQ IR Q 2 R is negligible in this regime, as the dynamics remain unchanged when we set d = 0. The features of the FFT spectra F IR and F R of the respective phonon modes in Figs. 3(c) and (d) therefore arise from sum-frequency ionic Raman scattering entirely. In this case, the frequency range of the cavity is close enough to the IR-active phonon mode in order to induce a shift even at the point marked as "OFF" resonant, as can be seen by a shift of the frequency in Fig. 3(c) away from the eigenfrequency Ω IR . In the resonant case, "ON", the peak of the FFT spectrum lies around 16.2 THz, which does not yet perfectly match the eigenfrequency of the Raman-active phonon mode, Ω R /(2π) = 16.5 THz. The consequences for the Raman-active phonon mode can be seen in Fig. 3(d). "OFF" resonance, the Raman-active phonon mode shows a small peak at 2Ω IR in addition to a tiny static component. In the resonant case, however, a large double peak emerges. This is a result of the Raman-active phonon mode simultaneously following the two-frequency component of the upper phonon polariton, 2Ω + , and gaining significant spectral weight at its eigenfrequency Ω R /(2π) = 16.5 THz with the same magnitude.
The IR-active phonon amplitude Q IR in Fig. 2(e) shows the common evolution of a resonantly driven and slowly decaying phonon mode when the cavity is out of resonance, as in regime I. With the cavity tuned into resonance, Q IR exhibits a small beating signal in addition, arising from the mutual energy exchange with the Raman-active phonon mode. The Raman-active phonon amplitude Q R is oscillating with small amplitude in the "OFF"-resonant case in Fig. 2(f) and gets strongly enhanced when the cavity is tuned into resonance with the IR-active phonon mode ("ON"). While the beating is small compared to the amplitude of the IR-active phonon mode, it is significant for the Raman-active one. This example shows that even when the IR-active phonon mode is not perfectly shifted into resonance with the Ramanactive one, the phonon amplitude can still be strongly enhanced. The increase of Q R by a factor of four in this case increases the efficiency of energy redistribution by ×16.

V. CONCLUSION
The ultimate goal of the technique presented here is to broaden the range of materials in which resonant nonlinear phononic processes can be exploited to yield nonequilibrium states of matter. Concepts in cavity control of nonlinear processes enable a new pathway for quantum optical engineering of new states of matter. The analysis presented here is applicable to resonant coupling mechanisms between IR-active phonon modes and other fundamental excitations in solids and molecules. In particular, we anticipate that it will be used to enhance processes involving the amplification of Goldstone excitations [44,52], the coupling to magnetic degrees of freedom [53,54], and the control of superconducting order parameters [55].
This work was supported by the DARPA DRINQS Program under award number D18AC00014 as well as the Swiss National Science Foundation (SNSF) under project ID184259. T.N. acknowledges the 'Photonics at Thermodynamic Limits' Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number de-sc0019140 that supported the computational approaches used here. J. F. acknowledges partial financial support from the Deutsche Forschungsgemeinschaft (DFG) under contract No. FL 997/1-1. P. N. is a Moore Inventor Fellow supported by the Gordon and Betty Moore Foundation and is a CIFAR Azrieli Global Scholar.