Kapitza stabilization of quantum critical order

Dynamical perturbations modify the states of classical systems in surprising ways and give rise to important applications in science and technology. For example, Floquet engineering exploits the possibility of band formation in the frequency domain when a strong, periodic variation is imposed on parameters such as spring constants. We describe here Kapitza engineering, where a drive field oscillating at a frequency much higher than the characteristic frequencies for the linear response of a system changes the potential energy surface so much that maxima found at equilibrium become local minima, in precise analogy to the celebrated Kapitza pendulum where the unstable inverted configuration, with the mass above rather than below the fulcrum, actually becomes stable. Our starting point is a quantum field theory of the Ginzburg-Devonshire type, suitable for many condensed matter systems, including particularly ferroelectrics and quantum paralectrics. We show that an off-resonance oscillatory electric field generated by a laser-driven THz source can induce ferroelectric order in the quantum-critical limit. Heating effects are estimated to be manageable using pulsed radiation; ``hidden"radiation-induced order can persist to low temperatures without further pumping due to stabilization by strain. We estimate the Ginzburg-Devonshire free-energy coefficients in SrTiO${}_{3}$ using density-functional theory (DFT) and the stochastic self-consistent harmonic approximation accelerated by a machine-learned force field. Although we find that SrTiO${}_{3}$ is not an optimal choice for Kapitza stabilization, we show that scanning for further candidate materials can be performed at the computationally convenient density-functional theory level. We suggest second-harmonic-generation, soft-mode-spectroscopy, and X-ray-diffraction experiments to characterize the induced order.


I. INTRODUCTION
The Kapitza effect [1-3] takes place when an anharmonic oscillator is dynamically stabilized at unstable equilibria by driving it at a frequency significantly higher than its natural frequency.The motion is divided into slow and fast components, with time averaging over the fast drive's short period enabling the analysis of the slowmode motion's stability in terms of an effective potential that depends on the drive's amplitude and frequency.The concept that fast degrees of freedom can modify slow behavior is not new, as illustrated by the thermal expansion of solids induced by rapid lattice fluctuations.This extreme rectification of fast fluctuations into a slowly varying (dc) behavior embodies the Kapitza engineering idea.
Kapitza engineering has been employed to stabilize repulsive Bose-Einstein condensates in an optical lattice [4], periodically driven spin systems [5], a many-body generalization of the Kapitza pendulum as a periodically driven sine-Gordon model [6], and a Josephson junction coupled to a nanomagnet [7].In each instance, the system is characterized by a nonlinear equation of motion for an effective degree of freedom coupled to a generalized field that induces externally controlled fast oscillations.
In this work, we suggest using Kapitza stabilization to dynamically control Quantum Matter near a quantum critical point (QCP).To develop this proposal, we use a model for the QCP state.As a prototypical example, we consider incipient ferroelectric (FE) materials like SrTiO 3 (STO) [8][9][10], where the QCP can be tuned by control parameters such as strain [11,12], Ca [13,14], and 18 O isotope doping [15].Incipient displacive FEs offer several advantages, including the availability of laserdriven electric fields E(r, t) directly coupled to the polarization.Additionally, the typical time scales associated with optical phonons responsible for displacive polarization lie within the femtosecond to picosecond range.We can therefore implement the Kapitza approach by applying fields in the accessible few-terahertz regime, as discussed below.
Observing the relatively long-lived polarization in a nominally paraelectric state due to applied time-varying electric fields would be a significant step towards dynamic control of the quantum state.Pristine STO is a quantum paraelectric.As an ionic insulator with a strong coupling between strain and polarization, the ferroelectric state in STO can have a long lifetime after the THz radiation is turned off.The relaxation of the induced FE state would occur on longer time scales associated with strain and lattice coupling.We envisage that this coupling helps in stabilizing the induced state without energy input producing Joule heating.
Our approach to QCP under drive follows the general approach to the Kapitza pendulum.We integrate the fast harmonics of P and find that new minima of the free energy for the FE order parameter can emerge, eventually leading to instability.To make a contact with realistic parameters we use ab initio analysis for free energy parameters for paraelectric state.In this work, we identify three essential features of the model required to induce the effect: i) a nonlinear potential, specifically the Devonshire-Landau free energy for the FE order parameter [8,37], ii) retardation in the response, and iii) application of the drive "off-resonance", where, as for the pendulum, the drive frequency of the applied field is chosen to be above resonance with the natural frequencies of the FE order parameter P, but not so high as to induce excitation above the insulating band gap.Due to these reasons, we refer to this approach as "Kapitza engineering", as discussed in Ref. 27.
We predict the required orientation of the optical electric field to observe the induction of QCP, contingent on the light's polarization type and the direction of the induced FE moment.Potential signatures of the induced ferroelectric order include low-frequency electrical measurements of displacement current, optical measurements of second harmonic generation, and X-ray diffraction measurements of the resulting strain.
The rest of the paper is structured as follows.In Sec.II we give a rather general exposition of the effective free energy of the FE and (linear) susceptibility renormalization in the presence of Kapitza drive from a strong, coherent, off-resonant light field.Then, focusing on particular crystal orientation, in Sec.III we analyze the critical electric field to make the PE state unstable (Sec.III A), the free-energy-gradient landscape under Kapitza drive (Sec.III B), and control of the QCP with light (Sec.III D).Our ab-initio calculations to extract the GD coefficients for STO are presented in Sec.III C. We discuss possible detection schemes in Sec.III E, with the details of the discussion moved to the corresponding Appendix.

II. MODEL AND GENERAL FRAMEWORK
We suppose that both the thermodynamics and the slow dynamics of the FE polarization P are governed by an effective action Γ[P] expressed as in integral over imaginary-time τ (sum over conjugate Matsubara frequency ν n = 2n π/T , where we choose units where k B = ℏ = 1) of the free energy [38].The effective action is composed of three terms -a harmonic term Γ 0 , an anharmonic term Γ ah , and linear coupling to the applied electric field E [39] Γ The harmonic part Γ 0 is given via the inverse susceptibility tensor of the unpolarized FE.Assuming translation invariance, it is best expressed in momentum-frequency space where the Fourier transform with a wave-vector q and a bosonic Matsubara frequency ν n , and Pj = P j /P 0 is the dimensionless polarization scaled by some convenient unit P 0 pertaining to the crystal.Starting from cubic symmetry, we scale down to tetragonal (orthorhombic) symmetry by scaling each component Pi by a suitable scaling factor ζ i .A prototypical Ansatz for the ionic inverse susceptibility tensor is that of a Lorentzian diagonal in the orthorhombic crystal axes and with transverse propagation where α(T, g) is the homogeneous, zero-frequency part, and coincides with the quadratic term of the Ginzburg-Devonshire (GD) free energy [9,[40][41][42].It is a function of temperature T , and, possibly, any parameters g that control the proximity to a QCP.Since it determines the inverse square of the correlation length, the zerotemperature characteristic exponent is 2ν.Mean-field theory predicts ν = 1/2, and other choices corresponding to a proximity to a QCP are considered in Sec.III D. We model the finite temperature dependence by another non-universal exponent n (1d) For a classical paraelectric n = 1 in accordance with the Curie-Weiss law.At low temperatures, a more convenient choice is n = 2 corresponding to a quantum paraelectric.There are more complicated cases [43], but this ansatz is sufficient for illustration purposes.The parameters c s and r are the sound propagation speed and damping, respectively, and ω 0 is a reference frequency (for example, the room-temperature value of the TO phonon frequency at the Γ point ω q=0 ).The results can be extended to include any dynamical exponent z [44,45].We choose the simplest case of a Gaussian theory with z = 1 for the rest of the paper.
The anharmonic action Γ ah is just a space-time integral over the (local) GD quartic free-energy density terms.
(1f) If α < 0, the FE would undergo a PE-FE phase transition along one of the main crystal axis.According to Eq. (1d), for g < g c , the Curie-Weiss temperature is given by T CW (g) = T 0 cs Λ ω0 (1 − g/g c ) 2ν/n .We assume the dynamics are governed by a particular optical phonon mode with a dispersion relation whose q = 0 frequency is √ α ω 0 .Because it is a polynomial, the susceptibility Eq. (1c) has a trivial analytic continuation from the experimentally measurable retarded susceptibility χret q=(q,ω) .Other Ansätze, more suitable for proximity to a quantum critical point with an anomalous dynamic exponent z may be used, but a prescription for the homogeneous q = 0, zero-Matsubara-frequency term ought to be given.The anharmonic term Eq. (1f) is the most general local quartic term consistent with cubic symmetry, and the scale factors ζ i account for lowered symmetry.Here, we focus on the dynamics in a plane perpendicular to the c-axis below the antiferrodistrotive transition [46][47][48][49][50].
The equation of motion for the q-th Fourier mode is obtained by extremizing Eq. (1a) Pq j Expanding the anharmonic term in Eq. ( 2) up to linear order in Pq around a reference static and uniform polarization ⟨P⟩, continuing analytically from Matsubara frequencies to the upper complex half-plane in ω, we obtain the effective linear susceptibility The retarded effective susceptibility determines the real-time response under an applied electric field with components along the principal crystal axes E(t) = E 0 ⟨cos (ϕ) cos (ω t), sin (ϕ) cos (ω t − δ), 0.⟩ (4) This gives rise to an oscillating component of the dimensionless polarization pi (π = cos (ϕ), sin (ϕ) e i δ , 0 ) The second moments over the rapidly oscillating polarization fields (remembering that χret (−ω) = ( χret (ω)) * ) evaluate to Using Eq. ( 6), we can express the equation of motion for the uniform and stationary solution ⟨P⟩ = T /L 3 P q=0 from Eq. ( 2) for q = 0 by averaging over the period of fast motion III. RESULTS Eqs. ( 3), ( 6), ( 7) are valid for any general linear susceptibility and anharmonic free energy density, and for ar-bitrary optical polarization π.For a quartic anharmonic free energy, the third derivative in Eq. ( 7) is linear in ⟨P⟩, and, thus, the Kapitza effect amounts to renormalization of the zero-Matsubara-frequency inverse susceptibility where We note that, for a non-zero ⟨P⟩, the second moments Eq. ( 6) are functions of ⟨P⟩, so the term ⟨P ⟩ j in the gradient of the free energy is still a highly non-linear function of ⟨P⟩.The logic behind going to Eqs. ( 8) from Eq. ( 7), and the detailed form of the third derivative are expounded in Appendix A.

A. Instability of the paraelectric state
We start by investigating the stability of the PE (⟨P⟩ = 0) point, by considering the sign of the eigenvalues of Eqs. ( 8), using Eq. ( 6).Because of Eq. ( 3), the retarded susceptibility is not renormalized in the PE state.
For definiteness, the plane of optical polarization is chosen to be the a − b plane.In that case, all the second moments ⟨ pi p3 ⟩ are zero for i = 1, 2, 3.This decouples the P 3 mode from the a − b plane modes.Furthermore, in the tetragonal phase in the c-direction is not renormalized.For the correction of the zero-Matsubara-frequency susceptibility in the a − b plane, we have a 2 × 2 sub-matrix where Ē = (ε 0 |χ ret (ω)| E 0 )/P 0 is a dimensionless measure of the amplitude of the electric field of laser light.We also note here that incoherent light with randomized ϕ (and δ) will lead to averaged matrix Eqs. ( 8) having vanishing off-diagonal components, corresponding to a zero ⃗ d ϕ in Eq. (9b), and inhibiting the possibility of Kapitza stabilization.
The eigenvalue of the tensor Eq. ( 9) that can become negative has a polarization direction pFE = ⟨cos (ψ), sin (ψ), 0⟩, where and the condition for instability is Ē2 2 The smallest magnitude Ē that fulfills the instability condition Eq. (10c) is in the direction ϕ for which the term in the square brackets acquires the largest positive value.From Eq. (9b), it is easy to see that the angulardependent part of that expression is Depending on whether |β2| α |cos (δ)| or |3β1−β2| 2α is greater, the optimal direction for the electric field is  7) under periodic drive for values of the parameters β1,2/α, together with the corresponding optimal electric field direction and critical electric field, corresponding to the white plus (a), and white cross (b), respectively, on the parameter diagram Fig. 2(a).The magnitude of the drive is E0/E0,cr = 3.0, and the frequency of the drive is ω/ω0 = 3.0 with a small damping r/ω0 = 0.05.We have checked that the choice of exact harmonic ratio of frequencies does not drastically alter the effect.The vector field is scaled according to the non-linear scaling vi = sgn(vi) 1 − e , with a suitable scaling factor v0.The white plus (cross) indicates the saddle point of the paraelectric phase, while the red points are the new stable fixed points.
ϕ opt = π/4 (direction [110]), or ϕ opt = 0 (direction [100]).Then, Eq. (10c) has two cases The condition corresponding to Eq. (10d) is fulfilled when the term in the square brackets is positive.This delimits three regions for the parameter The case for negative β 2 is unfeasible on grounds that the quartic GD free energy Eq.(1f) for the tetragonal phase is not bounded from below in the polarization direction [110] for β 2 < −β 1 (which explains the lower diagonal limit in Fig. 2(a)).We still have to verify that depicted by the green dashed line in Fig. 2(a).The first region is delimited by a dashed green line in Fig. 2.
The magnitude of the critical electric field is obtained by the value that saturates the inequality Eq. (10d) From Eqs. (10a), (10b), we see that the instability polarization direction is ψ = −π/4 (direction [ 110]), which is perpendicular to the direction of the applied electric field vector.
A similar analysis can be performed for the condition corresponding to Eq. (10e).The expression in the square brackets is positive when β 2 < 0. We still have to verify that which is outside the region of thermodynamic stability.The magnitude of the critical electric field is then and the instability polarization direction is ψ = π/2 (direction [010]), again, perpendicular to the applied electric field vector.
The intensity plot of these optimal fields is plotted in Fig. 2(a).
For directions of the applied field other than the optimal direction, the critical electric field magnitude is larger than the optimal, and Figs.2(b, c) depict the angular dependence of the magnitude of the critical electric field to the optimal one (red lines) for a representative point of the upper (lower) triangular regions in Fig. 2(a) (white cross, and white plus, respectively).We see that there are asymptotic bounds around the optimal directions, and the grey regions of the polar plot indicate combinations of direction and magnitude of the applied electric field that ought to stabilize an FE state.
We note that for sufficiently small magnitudes of the cross-coupling term β 2 in Eq. (1f) (the excluded white region in Fig. 2(a) between the two allowed triangular regions), we cannot have Kapitza stabilization.

B. Ferroelectric fixed points
Fixing the electric field beyond its critical value, we focus on the gradient field of the effective free energy Eq.(7).
For the same choice of parameters labeled in Fig. 2 by the white plus for the [100] direction and by the white cross in the [110] direction, the corresponding gradient field is depicted in Fig. 3.As predicted by the analysis in Sec.III A, Eqs.(10a)-(10b), the instability direction is in the a − b plane in a perpendicular direction to the applied electric field.We note that in both plots the PE point (P = 0) becomes a saddle-point, because in the direction of the applied electric field, there is only mode stiffening.This behavior of the gradient of the effective free energy is quite different than for the free energy in the absence of driving fields, which further reinforces the finding that cross-coupling of different polarization directions is crucial for Kapitza stabilization within this model.Additionally, the stable (saddle) points of the GD free energy, −α/β 1 [100], or −α/(β 1 + β 2 ) [110], are no longer extremal points of the effective free energy density.

C. Estimation of GD parameters from ab-initio calculations
We estimate the coefficients in the GD free energy density expression for the case of STO.The harmonic part of the GD free energy is given by the q = 0 contribution in Eq. (1b), remembering that Γ 0 = L 3 /T F 0 , and P q=0 = L 3 /T P. The anharmonic part takes the form of Eq. (1f).The fitting model then has the form We begin by using DFT to relax the internal coordinates of cubic and tetragonal (including antiferrodistortive rotations of the oxygen octahedra) unit cells, both in their centrosymmetric reference structures as well as with a displacement of the Ti atoms added in the [100] or [110] directions.The displacements break the symmetry and lead to a polarization in the respective directions.The size of the polarization is calculated based on the atomic displacements and fixed Born effective charges.The free-energy density as a function of polarization is then obtained by interpolating between the corresponding centrosymmetric and polar structures, and is shown in Fig. 4. The details of the DFT calculation are given in Appendix D.
One challenge is that standard DFT does not include temperature or quantum ion effects, and so predicts STO to be ferroelectric, yielding a negative quadratic term Eq. (11b).To capture the experimental quantum paraelectric behavior, and a positive quadratic term Eq. (11b), we use the Stochastic Self-Consistent Harmonic Approximation (SSCHA) in combination with machine learning force fields to include these effects in an abinitio description of STO at a reasonable computational cost [51].The details of these calculations are explained in Appendices E and F.
To obtain a corresponding energy versus polarization curve within the SSCHA method, for which the structure relaxes as expected to the non-polar phase, we calculate the SSCHA energies for the DFT-calculated polar/centrosymmetric structures scaled to the experimental lattice parameters [52].The results are shown in Fig. 4. We present the fit parameters in Table I.As expected, the quadratic term in the DFT case is negative indicating a characteristic double well potential with a ferroelectric instability from the non-polar phase.SSCHA on the other hand successfully yields a positive quadratic term.It also provides values for β ′ 2 and β ′ 1 which are both positive, but in a ratio of considerably less than the three required (Fig. 2(a)) for the Kapitza effect.Searches for other materials must therefore be undertaken.Interestingly Interestingly, the quartic coefficients obtained with DFT and SSCHA are similar, with β 1 almost identical, and β 2 reduced by ∼ 30% in SSCHA compared with DFT.We therefore expect that scans for further suitable materials with appropriate β values can be performed at the DFT level, with SSCHA used to check the values for the most promising candidates.

D. Control of the quantum critical point with light
The homogeneous, zero-Matsubara frequency inverse susceptibility α has characteristic scaling with temperature T for finite temperature or a control parameter g at zero temperature.Whenever it changes sign signals a critical point, either thermal or quantum-critical.Conversely, one can utilize α as a scaling variable to model the temperature or control parameter.Therefore, a helpful way to visualize the effect of the applied electric field is to plot the stable fixed point ⟨P ⟩ /P 0 as a function of α/β 1 and electric field magnitude ε 0 E 0 /P 0 for a particular choice of the ratio of quartic coefficients, as is done in Fig. 5(a).For negative α, we see that applying an electric field enhances the value of the FE order.For positive α (PE state), there is a critical region in the α-E plane where the stable order is zero.The line is simply the condition that the smaller eigenvalue of Eq. (9a) becomes zero, and we depict the power-law scaling of this curve in Fig. 5(b).The deviation of the exponent from 2 is due to the dependence of the susceptibility on α.For (b) small values of α, it is well approximated by 2.
We attribute the shift of the critical value for α with E as a shift of the critical transition point induced by light irradiation.This suggests a convenient parametrization of the dependence of the critical field E = E cr (α) on a control parameter T or g as a function of α in Eqs.(10f), (10g), and inverting by inverting Eq. (1d) with respect to T or g as a function of the parameter 1. Thermal transition at finite temperature and for the control parameter g < g c .From Eq. (1d), we can read off (note that α(0, g) < 0 for g < g c ) ∆T CW (g, α) 2. Quantum phase transition at zero temperature and for a control parameter g > g c .Again, from Eq. (1d), we can read off Both of these scalings are depicted in Fig. 5(c).For small values of α, the power-law dependence of the critical temperature is 2, regardless of the values of the non-universal exponent n, whereas the power-law dependence of the critical value of the control parameter 1/ν.

E. Experimental signatures of the polarized state
Concerning experimental observations of the proposed effect, we suggest the following experimental tests: (a) conventional dielectric measurements, (b) optical second harmonic generation, and (c) X-ray diffraction.a) In the case of an electrical measurement, one performs a poling measurement of the accumulated charge across a capacitor in an unbiased setup.The charge saturates during the irradiation and should vanish when irradiation is interrupted.The displacement current, which is the derivative of the charge, should show oppositely directed peak signals during the start and end of the irradiation.A sketch of the proposed setup and expected time traces of the signal are presented in Appendix B. b) In this case, a second harmonic to a weak signal at a frequency different than the drive frequency is generated when the sample is FE, with signal strength proportional to the induced polarization and specific angular dependence.Using the quantum effective action Eq.(1), we have obtained an expression for the non-linear susceptibility tensor, and the details are discussed in Appendix C.
(c) X-ray diffraction [53] measures directly the lattice strain associated with FE polarization and is nowadays routinely performed at free electron lasers with femtosecond time resolution.The experiment here will involve pumping the sample with THz radiation and then measuring the position of the principal Bragg peaks as a function of the time between the pump pulse and the Xray pulse.Because of the tensor character of the strain, which is mirrored in the resulting displacements of different Bragg peaks, it will be possible to reconstruct the full strain tensor as a function of experimental parameters such as THz pulse length and amplitude.
The strong coupling between strain and electric polarization is an important feature of the titanates and other (near) ferroelectrics.Generically, the latter builds up first, followed by the appearance of strain [53,54], which in turn stabilizes the polarization [22,23,55].This has two benefits: the first is that there is no need for a continuous power supply to maintain the FE state, and the second is that there is ample time to probe whether a FE state has actually been realized experimentally.
As far as an order-of-magnitude estimate for the necessary electric field magnitude, we use the following estimates for the GD free energy parameters [9, see Table I in Supplementary Information], a = α = 5.1 × 10 −5 , b = β 1 /P 2 0 .Taking β 1 /α = 0.5, implies P 0 = (β 1 /α) a/b = 1.9 µC/cm 2 .Then, using a magnitude of a complex permittivity εω = −60 + i 40 at a drive frequency of around ω/2π = 3.8 THz, which is twice the frequency of the softmode phonon [56], and a dimensionless critical field of the order of Ēcr = 2.8 (see Fig. 2), we get an absolute electric field magnitude E 0,cr = Ēcr P 0 /(ε 0 |χ ret (ω)|) = 830 kV/cm.For the estimate of the heat deposited in the sample by the laser field we use the Joule heating power density q = σ ω E 2 0 /2, where σ ω = ω ε 0 ε ′′ ω is the electric conductivity of the sample, which is estimated to be around 8.5 × 10 3 Ω −1 m −1 , giving an estimate for q = 29 W/µm 3 .To get an estimate of the heat fluence ϕ q , one has to remember that the electric field attenuates in the sample with a characteristic length scale √ εω , λ ω = 0.18 µm.The integrated power density over the sample thickness t and pulse duration τ gives ϕ q = 2 λ ω τ 1 − e − t 2λω q 0 .Assuming sample thickness much larger than 2λ ω , and pulse duration of at least 10 periods of the light field, the estimate for the fluence is ϕ q = 2.7 mJ/cm 2 .
It is important to acknowledge that the necessary power levels per unit volume may not be practically feasible in continuous-wave set-up.As previously mentioned, we must rely on coupling to strain to maintain polarization after the Kapitza "pump" has been removed.Notably, experience with dynamically induced nucleation and growth (at resonance) suggests that the relevant time scales for growth and decay lie within the femtosecond to picosecond range.Consequently, terahertz radiation pulses lasting around 1 picosecond and resulting in a net energy dissipation of 3 picowatts per cubic micrometer should be sufficient to observe the predicted effects.

IV. CONCLUSION
In this study, we extend the ideas of Kapitza engineering to the effective field theory near a quantum critical point (QCP) to show how the incipient order can be induced by the off resonance drive.As a specific example we considered the case of the strong off-resonance photon field applied to stabilize the ferroelectric (FE) phase in an incipient displacive FE.We find that it could be feasible to induce the FE phase if the material is close to a QCP and has a sufficiently large cross-coupling term in the anharmonic portion of the free energy.
We predict that the critical light field will exhibit a sensitive dependence on the direction and phase of the polarization axis.When field intensities surpass the critical field, the polarization will develop along an axis in the plane of the light's polarization and perpendicular to the applied electric field.
To make contact with materials we have developed detailed ab-initio microscopic models that enable us evaluate the realistic parameters of the Ginzburg-Devonshire action, such as β 1 , β 2 coefficients, and force fields in SrTiO 3 .The advantage of the proposed approach is in combining realistic materials parameters with the effective GD modeling of Kapitza engineering of FE state.While Kapitza engineering in SrTiO 3 would be a challenge based on the estimated parameters, we stress that we present the general framework.Our results: a) Establish a clear-cut numerical procedure that takes into account quantum-mechanical zero-point motion effects that render STO a quantum paraelectric; b) Demonstrate a deviation from isotropy (which corresponds to β 2 = β 1 ) and are not far off from the required stability criterion.The parameters are of the same order on the "wrong side" of the inequality Fig. 2(a), yet they are close.We draw the conclusion that the working regime in a similar class of quantum paraelectrics is possible.It is also possible that the choice of the functional would affect the relation between β 2 , β 1 ) and would bring them in the desired range, where Kapitza stabilization would be possible; c) These results together with the ab initio work serve as a guide for a well-defined process to be followed in different perovskite materials and under varying constraints in the search for suitable materials and optimal experimental constraints to make Kapitza stabilization possible.
We propose experimental signatures, such as poling measurements of accumulated polarization charge in an unbiased capacitor, second-harmonic generation in response to a weak probe pulse with an appropriate frequency and polarization axis, and X-ray diffraction measurements of the resulting strain.We also discuss the role of the Joule heating that would make the proposed set up impractical for the continuous-wave conditions.On the other hand, we estimate that finite-duration pulses with slow-strain dynamics could facilitate the stabilization of these phases.
Our proposed approach can be extended to other quantum-critical orders, including magnetism and superconductivity.We anticipate that factors like linear vs. quadratic coupling to drive fields will depend on the order's symmetry.The necessary resonant field strengths can be achieved with current experimental capabilities, making this mechanism a complementary method for manipulating quantum matter using terahertz light.Note: While this manuscript was in review, we became aware of a preprint [57] that is related to our work.Going from Eq. ( 7) to Eqs. (8) relies on the fact that the anharmonic free energy is a quartic function of polarization ⟨P⟩, and, so, its third derivative that enters in the third term in Eq. ( 7) is a linear function of ⟨P⟩ that we can ascribe to a renoramlization of the coefficient in terms of the linear term, which is the first term in Eq. (7).Having in mind the specific form of Eq. (1f), we can write the third derivatives in the following form This third derivative Eq. (A1a) gets contracted by ⟨ pj pk ⟩ in the second term on the right-hand side of Eq. ( 7).We get where we used a shorthand notation Π i,j for every coefficient in front of a Pj /ζ j term.These are different for i = j, and for i ̸ = j, and their form is given in Eqs.(8b), (8c) Susceptibility is a rank-three tensor that is defined as the third (functional) derivative of the cumulantgenerating function W Keeping in mind the Legendre transform relation where (P, E) is a shorthand for a space-time integral d 4 x P(x) • E(x), it is easily verified that the variational derivative of the effective action Γ[P] satisfies the equation of motion In its Fourier form, this is the same equation as Eq. ( 2) for a FE sample.The linear susceptibility is related to the second derivative.We can easily verify the following identities for the second derivative Using the chain rule, the third derivative Eq. (C1) can be expressed as where in the third line we used the matrix derivative identity Having in mind that the anharmonic free energy is local, the third variational derivative actually contains two Dirac delta functions Plugging Eq. (C6), and expanding the susceptibilities in Fourier modes, assuming the developed polarization is uniform and does not break translation symmetry, we get i,j,k (q 1 , q 2 ; ⟨P⟩), (C7) The Feynman diagram corresponding to Eq. (C8) is presented in Fig. 7.For light, the wavevector is close to q 1 = q 2 = 0, and Eq.(C8) is only a function of frequency.Because there is no frequency-dependent term in the third derivative (triangle vertes), it has a simple analytic continuation in the upper complex half-plane to get the retarded non-linear susceptibility χ i,j,k (ω 1 , ω 2 ; ⟨P⟩) = − χ ret eff (ω 1 + ω 2 ; ⟨P⟩) i,m P 0 χ ret eff (ω 1 ; ⟨P⟩) n,j χ ret eff (ω 2 , ⟨P⟩) p,k .(C9) tion set of 63 systems.The network was trained until no further improvements in the loss could be found during the last 50 epochs.Using the same test set as discussed in the supplementary material of Ref. [51] we obtain energy and force RMSEs of 0.21 meV/atom and 2 meV/ Å vs 0.18 meV/atom and 37 meV/ Å in Ref. [51].This is an 18-fold improvement in the force errors while retaining a comparable energy error.We used a cutoff radius of 6 Å, 4 message passing steps, and a maximum order of l = 2 for the equivariant message passing.The input config and output files of NequIP including the model are available at Ref. [58].
As the force field was trained using data created with slightly different DFT input parameters [51] we confirmed that the fit parameters obtained using solely the force field are sufficiently similar to the ones obtained with the DFT parameters used in this work.Using the force field for the tetragonal cells we obtain α ′ = −2.8× 10 8 J m C −2 , β ′ 1 = 2.0 × 10 10 J m 5 C −4 and β 2 = 0.57 × 10 10 J m 5 C −4 .While the double well becomes deeper with the force field the quartic coefficients and their ratio stay nearly unchanged.
Appendix G: Uncertainty in the fit parameters While the standard errors for the fit parameters are all at least one order of magnitude smaller than the parameters themselves, associating an exact uncertainty with the fit parameters is challenging, as many choices quantitatively influence the parameters.Examples include the direction of polarization chosen for the curves, the number of directions included in the fit, and whether the curves are calculated at constant or relaxed volume.Other factors are whether constant Born charges are assumed, and the choice of DFT parameters.We demonstrate the impact of some of these choices and the robustness of the results with respect to them.
Allowing for different volumes in the non-polar/polar structure and interpolating between them changes the ratio β ′ 1 /β ′ 2 to 2.5 (3.2 at fixed volume) for the DFT calculations on the tetragonal structure and to 3.5 (4.0 at fixed volume) for the SSCHA calculations.
We fitted the data in both polarization directions with a single linear fit.An alternative approach is to perform a separate fit of the [100] direction to obtain parameters for the quartic and square terms, followed by fitting the cross term for the [110] direction.For e.g., the cubic DFT calculations, this change yields α ′ = −0.94×10 8J m C −2 , β ′ 1 = 1.6 × 10 10 J m 5 C −4 , and β ′ 2 = 0.44 × 10 10 J m 5 C −4 .This difference demonstrates that more data in different polarization directions would be optimal to obtain more accurate fitting parameters, however the difference is also sufficiently small to not influence the conclusions.
In a final test, we use different Born charges calculated with the atomate2 dielectric workflow for the non-polar tetragonal structure.This changes the fitting parameters to α ′ = −1.

Figure 1 .
Figure 1.Diagram of the setup.A coherent laser field (blue line corresponding to electric field vector) is incident in the [001] crystal direction.The electric field is in the a − b crystal plane at an angle ϕ.The incident light can have a possible polarization phase δ, depicted by a spiral.

Figure 3 .
Figure3.Plots of the vector field, colored by its magnitude, of the gradient of the effective free energy Eq.(7) under periodic drive for values of the parameters β1,2/α, together with the corresponding optimal electric field direction and critical electric field, corresponding to the white plus (a), and white cross (b), respectively, on the parameter diagram Fig.2(a).The magnitude of the drive is E0/E0,cr = 3.0, and the frequency of the drive is ω/ω0 = 3.0 with a small damping r/ω0 = 0.05.We have checked that the choice of exact harmonic ratio of frequencies does not drastically alter the effect.The vector field is scaled

Figure 4 .
Figure 4. Energy versus polarization curves for the [110] (continuous line) and [100] (dashed line) directions calculated with DFT for the tetragonal (blue) and cubic (orange) phase and with SSCHA (green) for the tetragonal phase.The calculated results are presented as dots, while the fits to the data, which yielded the parameters listed in Table I, are shown as lines.The energies are plotted relative to the energy of the nonpolar structure.

Figure 5 .
Figure 5.Control of PE-FE transition with applied light.We choose a ratio of the quartic coefficients β2/β1 = 5.0 with the optimal polarization direction of the electric field vector in the [ 110] and linearly polarized light (δ = 0).(a) Stable FE order ⟨P ⟩ as a function of the magnitude of light's electric field E0 and quadratic coefficient α.The dashed line in the ⟨P ⟩ = 0 plane indicates a critical line.(b) The critical line from panel (a) presented on a log-log plot, with a linear fit indicated by a dashed line with the corresponding power.(c) Scaling of the critical temperature (left vertical axis (red)), or critical control parameter gc (right vertical axis (blue)) for several values of n and ν in Eq. (1d).

Figure 6 .
Figure 6.(a) Sketch of the unbiased poling measurement with a galvanometer (G) of an irradiated FE (blue sample) in a capacitative setup (gold electrodes).Illustrations of time traces for the accumulated charge and displacement current during irradiation

Figure 7 .
Figure 7.A Feynman diagram for the expression Eq. (C8), assuming q1 = q2 = 0 on the external labels.The triangle is the third derivative of the anharmonic free energy and is frequency-independent.The circles are linear susceptibilities at the corresponding frequency.