Local Fluctuations in Cavity Control of Ferroelectricity

Control of quantum matter through resonant electromagnetic cavities is a promising route towards establishing control over material phases and functionalities. Quantum paraelectric insulators -- materials which are nearly ferroelectric -- are particularly promising candidate systems for this purpose since they have strongly fluctuating collective modes which directly couple to the electric field. In this work we explore this possibility in a system comprised of a quantum paraelectric sandwiched between two high-quality metal mirrors, realizing a Fabry-Perot type cavity. By developing a full multimode, continuum description we are able to study the effect of the cavity in a spatially resolved way for a variety of system sizes and temperatures. Surprisingly, we find that once a continuum of transverse modes are included the cavity ends up suppressing ferroelectric correlations. This effect arises from the screening out of transverse photons at the cavity boundaries and as a result is confined to the surface of the paraelectric sample. We also explore the temperature dependence of this effect and find it vanishes at high temperatures, indicating it is a purely quantum mechanical effect. We connect our result to calculations of Casimir and Van der Waals forces, which we argue are closely related to the dipolar fluctuations in the quantum paraelectric. Our results are based on a general formalism and are expected to be widely applicable, paving the way towards studies of the quantum electrodynamics of heterostructures featuring multiple materials and phases.


I. INTRODUCTION
Optically engineering properties of quantum materials may potentially allow for the design and development of novel technologies and the creation of phases of matter which are otherwise difficult to obtain and study.Ideally, one would like to think of this as adding a new "control knob" to the toolbox of solid-state physics just like temperature, pressure, external field, and twist-angle, allowing for new explorations of physical systems and device structures [1].For instance, intense electromagnetic radiation can induce nonequilibrium phases of matter and generate new phase diagrams, with sometimes no counterpart in equilibrium [2][3][4][5][6][7][8][9][10][11][12][13].However, the nonequilibrium route towards optical control has a number of drawbacks which impede its practical application, chief among them are problems related to heating, optical access, and complicated theoretical modeling.Therefore, it would be desirable to obtain a similar degree of optical control without leaving thermal equilibrium.Recently, it has been argued that this may be done by instead shaping the environment of electromagnetic fluctuations through the use of optical cavities, resonators, and metamaterials [14,15].Many systems have recently been proposed to be amenable to control in this way, including superconductors [6,[16][17][18], excitonic insulators [19], antiferromagnets [20,21], spin-liquids [22], semiconductors [23], quantum Hall fluids [24], and ferroelectrics [25][26][27][28][29], with a great deal more proposed to exhibit strong coupling between material and optical excitations [30,31].Recent * jon.curtis.94@gmail.comexperiments on the metal-insulator transition in 1T -TaS 2 even seem to have seen promising signatures of cavity control on the transition temperature [32].Experiments have also seen fascinating phenomena occur when unconventional superconductors are strongly coupled to the quantum electromagnetic bath [33,34].
Ferroelectrics are particularly promising candidates since the relevant fluctuations-phonon polaritons-directly couple strongly to the electromagnetic field via the electric dipole transition even down to atomic scale [31,[35][36][37].Furthermore, there are a number of appealing candidate systems, such as SrTiO 3 [38][39][40][41][42] and various moiré and Van der Waals materials [43][44][45][46][47] which may be suitable for proof-ofprinciple experiments.Intrinsic SrTiO 3 is believed a quantum paraelectric (QPE) [38][39][40][41][42], lying right at the border of the ferroelectric phase, with long-range order suppressed by quantum fluctuations.Strain, chemical, and isotope substitution have all been shown to tip the system over the edge in to the ordered phase [39], and recently resonant optical excitation of the lattice [40,48] have also been shown to seemingly induce a transition into the ordered phase [49,50], making this a prime candidate for demonstrating cavity control over the phase diagram [25,26,28].Previous theoretical investigations have largely been limited to single-mode, dipole coupling, and translationally invariant approximations.In this Article we show that going beyond these simplifying assumptions can lead to a qualitatively different behavior making our approach necessary in order to describe realistic experiments [51].
In this work, we extend our analysis of fluctuating quantum paraelectric to a fully multimode [23], spatially resolved system beyond the standard dipole approximation -the standard optical approximation which treats the electromagnetic field as homogeneous across the sample.In fact, this is a crucial technical development, giving a number of predictions which starkly differ from previous simplified models.By making use of connections to the study of Casimir-Polder and Van der Waals forces, we are able to efficiently solve the full multimode problem including a continuum of electromagnetic modes.In doing so, we find that in fact the presence of the cavity suppresses ferroelectric fluctuations in the system-the opposite of what is expected based on a single-or fewmode calculations.This surprising result then has important implications for future experiments on cavity control of ferroelectricity and potentially other phases of matter [32][33][34].
The key insight is using the fluctuation dissipation relation to reformulate the problem in terms of the dielectric response and its variational dependence on material parameters, thereby encapsulating the effect of electromagnetic fluctuations in terms of well-known frequencydependent electric-field correlation function.The behavior of this correlations function is very well studied, dating back to seminal work on the Casimir force [52][53][54][55] and is by now well documented and experimentally verified.In fact, cavity control over the QPE fluctuations in a material is closely related to the problem of using the cavity to modify the Van-der Waals forces between virtual dipolar excitations in the cavity [56].Furthermore, we argue our technique can be easily extended to include phonon loss, anisotropy, and mode couplings provided the dielectric constant dispersion ϵ(ω) is known well enough, and may even be extended to include more complicated heterostructure geometries such as interfaces between quantum paraelectrics and metals, air, or more exotic twodimensional systems via characterization in terms of the reflection coefficients at interfaces [55].In a rough sense, the problem is similar to calculating the Casimir force but instead of computing the energy as a function of the boundary separation, we keep the boundary conditions fixed and directly compute the photon and phonon fluctuations inside the cavity.
After obtaining these general relations, we use our method to compute the local behavior of the QPE fluctuations by considering a Fabry-Perot type system with perfect metallic boundaries sandwiching a QPE, illustrated in Fig. 1 in a cross-sectional view.We find as our key result that actually towards the boundaries of the sample the phonon fluctuations ⟨Q 2 (r, t)⟩ increase, resulting in a blue-shift of the soft-mode transverse frequency Ω T due to the anharmonic coupling of the phonon mode, characterized by the local displacement field Q(r, t).This leads to an overall thickness dependence which may be pronounced for thinner cavities and results in a diminished low-frequency dielectric constant ϵ(0), signaling a hardening of the soft polar mode.
We also find that this effect-the difference between the surface and bulk fluctuations-vanishes at higher The system has translational symmetry in the xy plane, which leads to a conserved in-plane momentum q, which can be taken along q ∥ êx, and obeys ideal metallic boundary conditions at z = ±L/2.Due to metallic boundary conditions, the photonic part of the wavefunction is suppressed near the boundary, leading to enhanced phonon fluctuations.
temperatures, indicating the origin of this effect is of a truly quantum origin, and should drop off once T ≳ ℏΩ T /K B .For materials such as SrTiO 3 , with Ω T ∼ 2THz, this provides an important ceiling on the temperature of experiments which is roughly of order 100K.
The remainder of our Paper is structured as follows.In Sec.II we use the fluctuation-dissipation theorem to connect the phonon fluctuations to the dielectric response of the material.We first do this in real-time formalism in Sec.II A, followed by a reformulation on the Matsubara axis for finite-temperature calculations in Sec.II B. In Sec.III we demonstrate how our results connect to the more familiar method based on phonon-polaritons in the case of a bulk translationally invariant system.In particular, in Sec.III A we show how at high-temperatures the photons decouple highlighting that the cavity control is a manifestation of quantum effects.In Sec.IV we apply this technique to derive the local QPE fluctuations in a Fabry-Perot geometry which does not preserve translational symmetry.We then conclude by discussing general aspects of our results beyond our cavity-QPE model and highlighting potential directions for future study in Sec.V.

II. FLUCTUATION-DISSIPATION THEOREM
In the following we will focus on the case of a local [57], isotropic, polar phonon mode.In particular, since the phonon group velocity is much slower than the photon, the approximation of a non-dispersion phonon mode should be suitable for studying electrodynamic effects.Going beyond this approximation to include the phonon dispersion would be an interesting direction for future studies and may be important very close to the ferroelectric critical point or in the ordered phase, which we will not study in this work.
We thus consider a model of a local polar phonon mode Q(r) with transverse optical (TO) mode frequency Ω T and effective charge η coupled to the electromagnetic field, described by E and B. This system is most simply described in terms of a Lagrangian which generates the equations of motion.We have The Maxwell Lagrangian is In terms of the gauge potentials, the electromagnetic fields are expressed as Here and throughout we use units where The phonon Lagrangian, in the absence of dispersion, is purely local and is simply given by: where Ω 0 is the bare TO mode frequency, and λ is the phonon-phonon interaction strength.Due to symmetry, these are the only local terms allowed at up to quartic order and second order in time-derivatives.Finally, we have the dipole-coupling between the phonons and the electric field.The phonon displacement field Q generates a polarization which then couples to E, via Here, η is the light-matter interaction constant and it sets, among other things, the size of the splitting between the longitudinal optical (LO) and transverse optical frequency splitting due to the Coulomb part of the electromagnetic interaction (the so-called LO-TO splitting).Before proceeding, in order to perform calculations we must fix a gauge.In this work, we will employ the "Weyl gauge", which is obtained by demanding A 0 = 0. We proceed by treating the phonon nonlinearity in the Hartree approximation, such that the system is essentially linear, albeit with a renormalized phonon frequency of In general, we allow for spatially varying fluctuations of Q, and therefore the phonon frequency may be renormalized in an inhomogeneous way, which is the subject of this investigation.Therefore, our primary objective is to compute the spatially resolved phonon fluctuations, ⟨Q 2 (r, t)⟩.We do this by the familiar linear-response formalism, obtaining the equilibrium fluctuations of Q by solving for the causal response to an external perturbation.We thus introduce a source field F(x) which couples to the phonon mode via such that ) for causal response function.Here and throughout, when confusion is not likely, we will use x = (r, t) to represent a spacetime four-coordinate, while r as a spatial three-coordinate.From this we can apply the fluctuationdissipation relation [58] to obtain This then allows to characterize the local density of phonon fluctuations.

A. Real-Time Equations of Motion
The relevant equations of motion can be written down, including the source term which acts on the phonon field.This gives us the equations in the frequency domain We compute the response function of Q as follows.Let us first introduce the bare response function for Q(r, ω) of such that we can obtain the response of Q as Our job is not done because we need the response of the phonon not to the total force, which is ηE + F, but only to the external force F, which is partly screened by the electromagnetic field.This screening is found by solving the equations of motion for the electromagnetic field.We have The second equation contains the dependence on the forcing field; we can eliminate the magnetic field to derive a closed equation for the response of E, which we use to find the depolarizing field, from which we find the effective force acting on the phonon mode due to the external force.We get Here we have introduced the dielectric constant We now can use this to eliminate the electric field formally as where ĜR (r, We can evaluate this Greens's function in a manner of our choosing; in a bulk system it makes sense to use momentum space functions.
The full response of the phonons, dressed by the photons, is given as a function of the photon Green's function: In the absence of phonon dispersion the first term is completely local, and exhibits a resonance at the bare TO mode frequency, while the second term involves dispersion of the phonon polaritons in the system, as it involves ϵ(ω), and therefore is sensitive to the cavity geometry.
We note we can write this in an elegant way by using to obtain

B. Matsubara Formalism
The result presented is derived using a real-time linearresponse formalism, which is the most transparent presentation.In Appendix A, we equivalently derive this result using the equilibrium Matsubara frequency representation, as well as yet a third way of deriving this result based on a variational procedure for the Matsubara freeenergy functional in Appendix B.
The Matsubara frequency representation is particularly useful since it allows for a more efficient implementation of the result in terms of a well-behaved, convergent sum over Matsubara frequencies.For more details, we refer the reader to the relevant Appendices.However, we present the end formula here as it is relevant for the results to follow.
The dielectric function is analytically continued to bosonic Matsubara frequencies ω → iω m = 2πimT as allowing for a locally varying TO mode frequency.Likewise, the unscreened phonon propagator becomes and the photon propagator becomes Ĝ (r, The local phonon fluctuations as a function of temperature are calculated by combining equations ( 9),( 21), which is conveniently expressed in terms of a Matsubara sum as Note that ∂ϵ(iω m , r) is negative semidefinite, while the other terms in the equation are positive semidefinite.This leads to a natural conclusion that the electric field actually has the effect of reducing phonon fluctuations, since it subtracts spectral weight away from the otherwise unscreened phonons.We also see that the result can be implemented by a sum over terms which are non-singular, greatly aiding the numerical evaluation of this expression.This is somewhat similar to the recent approach of Ref. [28], though extended to the multimode formalism.Nevertheless, the concept of a 1/N expansion would be useful to apply in this context as well.
We can then ultimately express quantities in terms of the effective frequency renormalization as compared to the bulk value For each configuration, this is done by extracting the bulk value as the value computed as system size tends to infinity.We can then visualize the predicted spatial deviations of the frequency away from the bulk value.We now roughly estimate the size of the coupling λ, the momentum space cutoff Λ, and other relevant parameters.

C. Parameters
Following Ref. [48], we first extract the value of the anharmonic potential U (2) as a function of the ionic displacement coordinate u, for the case of SrTiO 3 , estimated from spectroscopy to give ν T ∼ 35THz 2 /pm 2 .This is obtained from nonlinear terahertz spectroscopy by measuring ).However, this is not directly related to the long-wavelength order parameter [59], which in general involves a complicated coarsegraining of the microscopic degrees of freedom.By dimensional analysis, we see that the long-wavelength coupling constant can be related to a microscopic parameter through λ ∼ ν T V eff /M eff where M eff is an effective mass scale and V eff is an effective volume scale.For the purposes of our crude model, we use a single effective mass to characterize the mode, which we take to be the titanium ionic mass Finally, in the spirit of a real-space renormalization group procedure, we anticipate that the effective volume V eff which appears in the relation between the orderparameter and the microscopic degrees of freedom should be related to the momentum space cutoff we place on our model, Λ ∼ 1/a with a the size of the coarse-grained "blocks."As a result, we estimate that V eff ∼ 1/Λ 3 ∼ a 3 .We thus identify the long-wavelength coupling as λ = a 3 ν T /M Ti .This leads to the cutoff-dependent estimate of which we use in this work.Finally, by ignoring the gradient terms ∼ ∇Q, our model treats the phonon fluctuations as purely local.In reality, we expect this approximation to breakdown below a length scale a, roughly related to the correlation length of the ferroelectric order parameter Q(r).On the one hand, near the critical point this length scale may become large as correlations develop across longer length scales.One the other hand, treating physics at a length scale l < a requires a more complicated theory than the one we develop here.In this work we are interested in macroscopic physics, with samples of order L ∼ 50µm or larger; as a result, we will consider a cutoff of a = 5µm.
In the future it would be particularly interesting to consider going even closer to the critical point, so that a ∼ L, in which case electrodynamic control may be even more important.This is because the photon dispersion is much faster than the phonon, and thus the region in reciprocal space amenable to cavity control is roughly q 3 QED ∼ (1THz/c) 3 ∼ (.3mm) −3 , which is to be compared against the volume of the Brillouin zone which exhibits fluctuations, which goes roughly as 1/a 3 .As such, we expect the physics we consider to be most important materials which are close to the ferroelectric instability, such that a −1 ∼ q QED .However, in this regime a more elaborate model which considers the role of phonon dispersion and spatial gradients is required.Evidently, a more elaborate scaling theory of the transition is needed in the future.We elaborate on this slightly later, in Sec.V.
Finally, we must fix the phonon TO and LO mode frequencies.For the bulk TO frequency Ω T and LO-TO splitting ∼ η, we take Ω T = .5THzto emulate the ferroelectric soft mode, and η = 10THz, such that ϵ(0) ∼ η 2 /Ω 2 T ∼ 400 is large, as is the case for many incipient ferroelectrics.We also ignore the temperature dependence of these parameters for simplicity, though this could be accounted for more careful if we were modeling a specific material.
Finally, for numerical purposes we take a Matsubara frequency cutoff of ω c = 40THz, well above all other frequency scales; this is not expected to be an important parameter, provided it is large enough.In order to make the calculations simpler and more transparent we also relax the self-consistency demand on the TO frequency, and instead simply evaluate the perturbative correction to the bulk value.In principle, this should be addressed though if the shift is small it is likely to be qualitatively correct.

III. HOMOGENEOUS SYSTEM
Before proceeding on to our main calculation, we briefly outline how the calculation works in the case of a system with full translational symmetry.In this case, momentum is a good quantum number in all directions.We can then go to the plane-wave basis. Using we find the electromagnetic correlation function splits into two decoupled subspaces: the transverse modes and longitudinal modes, with We easily recognize the first term as the dynamical Coulomb portion of the electric field, while the second term comes from the transverse photon modes.The unscreened phonon propagator is trivial in this case, with We find for the overall result, with UV cutoff on momentum Λ, This simplifies into one longitudinal mode and d−1 transverse modes (with d the spatial dimensionality), with the longitudinal modes contributing (per momentum space mode) The momentum space integral will then simply yield a factor of Λ 3 /6π 2 , with the cubic divergence due to the dispersionless nature of the model.The transverse modes are more difficult to evaluate.Evaluating the Matsubara sums and performing the analytical continuation to real frequencies reveals that the overall result including the longitudinal and two transverse modes is The dispersion of the upper polariton branch (Ω +,q ) and lower polariton branch (Ω −,q ) are found to be This is in accordance with the result one would expect from Bogoliubov transformation and diagonalizing the quasiparticle Hamiltonian.

A. High-Temperature Limit
We also can look at the temperature dependence of this effect, and in particular at high-temperatures we obtain, after simplifying terms This is exact and independent of momentum q.
In this case we see the electromagnetic field only serves to split the longitudinal modes away from the transverse modes, which then essentially decouple form the electromagnetic field fluctuations.All the non-trivial electrodynamic effects end up vanishing as 1/T , which can be interpreted as them being a true manifestation of quantum effects.This is evident from examining the relevant Matsubara sums, which end up going as ω 2 m /q 2 at low frequencies, and therefore vanish in the static ω m = 0 2. Schematic depiction of temperature dependence of phonon fluctuations ⟨Q 2 ⟩ as a function of temperature in the bulk and near the surface of the sample-cavity boundary.At high temperatures the fluctuations are independent of electrodynamic details and therefore are ignorant of the proximity to the surface.At lower temperatures, the electrodynamic contribution becomes important and leads to a more pronounced fluctuations near the surface as compared to the bulk.
limit.The finite Matsubara frequencies are in turn arising from dynamical quantum fluctuations of the field, which are ultimately the source of the polaritonic splittings.We emphasize that this does not mean that the fluctuations of ⟨Q 2 ⟩ diminish with increasing temperature, but rather that the contribution from the electromagnetic field falls off.As a result, if one performs a hightemperature expansion of ⟨Q 2 ⟩ in the bulk one would get ⟨Q 2 ⟩ bulk ∼ a bulk T + b bulk /T + ..., with the leading term the classical equipartition result and the subleading terms coming from the high-temperature expansion of coth(ωβ/2) (which is odd, so the series has only odd powers of temperature).Performing the same expansion for the fluctuations near the surface (which as we will show, differs in that it doesn't couple to the transverse photons due to boundary conditions) we get ⟨Q 2 ⟩ surf ∼ a surf T + b surf /T + ... and find that a bulk = a surf , so that the effects of the boundary conditions are subleading in T and due to truly quantum effects.This is illustrated in Fig. 2 schematically, which shows that at high temperatures the classical result dominates and only at low temperatures is the interaction with photons relevant.This is reminiscent of the Bohr-Van Leeuwen theorem in classical mechanics, and can be heuristically understood as a consequence of the current and displacement being canonically conjugate for the phonon.This is explored more technically in detail Appendix C; here we provide a simple heuristic.Essentially, photons don't actually couple directly to the phonon displacement Q, but rather couple to the displacement current J ∼ ∂ t Q associated to transverse oscillations in the charge distribution.This current is canonically conjugate to the ob-ject of interest, Q, via [Q(r), J(r ′ )] ∝ iℏδ 3 (r − r ′ ) and therefore in the quantum system (at low temperatures) they are not independently fluctuating.Thus, modifying the fluctuations of the current J can in turn induce changes in the fluctuations of Q.However, this coupling is purely due to the canonical commutation relations between the two fields, and therefore at high-temperatures (i.e. in the classical limit), the two variables become independently fluctuating, just as q and p are independent in a classical system.Therefore, the photon decouples from the fluctuations of Q since this is now independent from the current.We also see that this is not the case for the longitudinal part, which instead directly couples the electrostatic potential to the induced phonon charge ρ ∼ ∇ • Q, and indeed the LO-TO splitting due to this interaction remains unaffected in the high-temperature classical limit.
To summarize, we see directly from this calculation that the coupling to the electric field reduces the fluctuations of the phonons, both by dressing the eigenfrequencies, and also by reducing the projection of the quasiparticle wavefunction onto the phononic subsystem.Therefore, we uncover a counter-intuitive heuristic that coupling to the quantum electromagnetic field actually facilitates ferroelectric order, by virtue of reducing the fluctuation-induced shift of phonon frequency.In the following section we will apply this formalism to the case of an inhomogeneous slab structure, and demonstrate how this is potentially visible in terms of a local shift in the phonon mode frequency.

IV. PLANAR GEOMETRY
We now evaluate the correlation functions needed to compute the frequency shift.In particular, we focus on the electromagnetic Green's function since the phonon correlation function is local and trivial to solve.More over, this term will essentially be a bulk background contribution, and in this work we are focused on the contribution from the degrees of freedom we can change by the boundary conditions (the photons).
In the planar geometry we can utilize in-plane translation symmetry to reduce the problem to a one-dimensional differential equation in terms of the spatial coordinate z.It turns out this can still be solved analytically utilizing the transfer matrix method, at least in the case of a homogeneous dielectric constant, as was originally done a long time ago for the purposes of calculating Casimir-Polder forces (which turns out to be a closely related problem) [53][54][55].We can take the momentum to lie in the x-direction, such that the Green's function for the gauge potential G reads (37) This system is depicted schematically in Fig. 1, illustrating the metal-paraelectric-metal geometry.In this work, we take the physically reasonable limit of infinite plasma frequency in the metal, such that the metal mirrors can be modeled by Dirichlet boundary conditions on the tangential components of E and normal component of B.
We take the cavity plates to be located at z = ±L/2 such that L is the total size of the cavity, and also the full extent of the paraelectric is all the way up to the boundaries.This is solved in detail in Appendix D, we will only present the final result here.We find that the trace of the Green's function evaluated at coincident spatial points trG (z, z) is given in closed form as Here κ = ω 2 m ϵ(iω m ) + q 2 governs the length-scale for the recovery to the bulk value for a given frequency and in-plane momentum.We see that this involves the divergent quantity δ(0), which is understood as lim z ′ →z δ(z − z ′ ).This quantity is in fact independent of system-size and geometry and thus ends up getting renormalized away by the counter-term Ω 2 0 .In particular, we only compare this result to result obtained for L → ∞ (with z finite), which ultimately gives the closed-form formula for the renormalization of the phonon-frequency shift due to the cavity, as a function of position, as We note this expression does still depend on the cutoff Λ = π/a.While this is much easier to evaluate than the full Green's function, it is still not completely trivial due to the subtle behavior of the integrand at ω m → 0, which governs to the high-temperature limit of the effect.For all non-zero Matsubara frequencies, this is easily evaluated by numerically summing, for ω m = 0 there is an issue about how the limit of zero frequency is taken; on the one hand, the numerator involves a power of ω 2 m which then vanishes quadratically at small frequency.On the other T = 10K FIG. 3. First order correction to phonon frequency due to cavity boundary conditions at z = ±L/2.We here fix the temperature to T = 10K, which is lower than the coherence temperature for the phonon frequency of ΩT = .5THz(we note the coherence frequency involves a factor 2π so that 2πT ≲ ΩT ).We explicitly note that this is the correction to the phonon frequency due to the cavity and in particular, this blue shifts at lower temperature, unlike the typical phonon anharmonicity which blue shifts at higher temperatures.The blue shift due to the anharmonicity is renormalized away in this treatment; here we only plot the additional shift observed between the bulk and finite systems at the same temperature.
Careful examination of the limit reveals that the numerator ends up winning and thus, we are to discard altogether the zeroth frequency contribution.This is again a manifestation of the quantum mechanical origin of the effect, as explained in Sec.III A and shown in greater detail in the Appendix C. We are now able to efficiently evaluate this effect in order to obtain not just the phonon frequency shift but its entire spatial dependence.We also comment that this expression is clearly manifestly positive, owing to the inequality sinh κ(L/2 − z) sinh κ(L/2 + z)/ sinh κL ≤ 1 2 .As a result, we find the frequency shift is always positive, with ∆Ω 2 T (z) ≥ 0, such that the cavity in fact suppresses the ferroelectric order.We will comment on this more later.
In particular, it is interesting to analyze the dependence of the renormalized frequnecy on system size L and temperature T (we consider the value in the midpoint of the cavity).This is shown in Fig. 3.We indeed see that for "bulk systems" with large system size L, there is no effect due to the boundary conditions (which is expected on grounds of locality).For smaller systems however, the typical phonon frequency strongly blue shifts as the surface effects set in.For very small systems, the entire system is essentially "surface" and thus we see the first characteristic prediction which is a significant blue-shifting of the soft-mode frequency for thin samples.
We next confirm the temperature dependence of this effect; namely, that at high-temperatures any signature of the cavity should disappear.This is seen explicitly in Fig. 4. Indeed we see that in all system sizes, the frequency shift vanishes at high temperature irrespective of FIG. 4. First order correction to phonon frequency due to the cavity as a function of temperature T .Here we consider two cases: a small system with L = 50µm (a), and a larger system with L = 70µm (b).In the inset we depict schematic profiles of the renormalized frequency, illustrating how the effect is larger for smaller systems since the surface effects are still dominant, whereas in a larger system the frequency converges to the bulk value.
the system size L. At lower temperatures, the phonon frequency shift sets in, but for larger systems the effect is small since it ultimately must recover to the bulk as L → ∞.For smaller systems however, the shift may become sizeable at the cavity midpoint.
These results are best understood by simply looking at the spatial profile for the renormalized frequency Ω 2 T (z), depicted in Fig. 5.For a cavity size of 100µm we see that the phonon frequency significantly blue-shifts near the boundary at low temperatures, while it remains essentially equal to the bulk value at high temperatures and in the center of the cavity z ∼ 0).This is inline with the previous arguments we have given, namely that the cavity renormalization is (i) a purely quantum effect setting in once 2πT ∼ Ω T , and (ii) that it is a surface effect in response to the boundary conditions at z = ±L/2.We also always see that the frequency blue-shifts-this is in some tension with a number of previous theoretical investigations [25,26,28,29], which all seem to find that cavities tend to enhance ferroelectric order.We now reconcile our calculation with these previous studies.
We believe the origin of this tension is in the way that the UV cutoff on the electromagnetic fluctuations is imposed, and ultimately scaled to the continuum limit (or not).In particular, in a real physical system (with an appropriate UV cutoff) the number of electromagnetic modes is proportional to system size, with a finite number of modes (e.g.lattice points) per unit volume, albeit of potentially high frequency and thus far from resonance.We again emphasize that it is qualitatively important to take into account the multimode nature of the electromagnetic fluctuations [23,60].
In our new formalism, which up to, and including Eqn. ( 39) is an analytically exact solution to the problem, the UV cutoff is essentially scaled in this fashion, so that the density of modes per unit volume is constant.This is also in congruence with known and experimentally verified results on Casimir forces [52][53][54][55].In contrast, previous studies have essentially imposed a cutoff on number of cavity eigenfunctions taken, selecting the lowest N c modes to include in an eigenmode expansion for the Green's function G regardless of system size.This assumption leads to the unphysical result that number of electromagnetic modes per unit volume ∼ N c /L 3 is scaling to zero for a bulk system!It turns out that, unfortunately this has exactly the opposite behavior as what we find using the more complicated analytical solution [61].
While our approach draws upon parallels with Casimir forces, it is however investigating a distinct phenomenon and thus it still warrants experimental confirmation.We now briefly discuss the prospects for this in the next section.

V. DISCUSSION
As we saw in the previous sections, the effect of the cavity ends up being largely confined to surface of the material, with the net result being a blue shift of the soft phonon mode near the boundaries.This would result in an apparent diminishing of the dielectric constant ε(ω = 0) ∼ η 2 /Ω 2 T for a small sample in comparison to the bulk value.In fact, this effect has already been observed and known for quite some time as the "deadlayer" effect [62].Canonically, this effect was found in exactly this sort of system, comprised of a thin layer of SrTiO 3 sandwiched between two electrodes to form a capacitor [63].Historically, the dead-layer effect has been attributed largely to strain effects [64] between the two interfaces which also acts to blue-shift the soft mode, however it seems the matter has not been entirely settled [65].It is quite interesting then in the light of this new work and interpretation to re-open the investigation and determine if any effects due to quantum electrodynamics may be relevant.To this end, experiments may want to try experiments with different metallic interfaces that induce varying levels of strain, to see how sensitive the dead-layer is to the level of strain, as opposed to the dielectric environment.Indeed, experiments investigating the interplay between electromagnetic response and the surrounding dielectric environment [66] have already shown promising results [67][68][69].Theoretically, this would also warrant further calculations which have open (air) boundary conditions for the SrTiO 3 layer, rather than the metallic interface conditions.
Our results also point to another interesting connection which ties our results to studies of the Casimir force [52,54,55].This is most apparent if one returns to Eqn. 39, which closely mirrors the results for Casimir forces (see, e.g.Ref. [55]).This makes sense since our setup is very similar to the one considered originally by Casimir except in this case the fluctuation force does work on the dielectric constant of the QPE (which is determined self-consistently) rather than the plates of the cavity [41].It may be interesting to turn this around and try to utilize Casimir force spectroscopy to probe incipient or critical ferroelectric fluctuations via the effect of fluctuations on radiative forces.
A closely related phenomenon is that of the Van der Waals force, which is also an entropic force due to virtual electromagnetic fluctuations.Van der Waals forces typically act between polarizeable molecules and is notably attractive, leading to the total energy being lowered as a result of the coupling to electromagnetic fields [52][53][54].Similarly, the infrared-active phonon field can be thought of as a lattice of polarizeable molecules undergoing virtual polar fluctuations.The primary difference between the two cases is that in the case of the phonon polaritons the molecules are arranged in a regular lattice, whereas in a fluid the molecules are spatially disordered.When placed in a cavity, the metallic boundary serves to screen the electromagnetic field, leading to a net increase in the free-energy as compared to the bulk system, since the Van der Waals forces which get screened out are attractive.
In fact, the extended many-mode nature of this problem is paramount.This is seen by contrasting our results with the simpler case of a single fluctuating dipole localized near a metallic surface.In the localized case, the interaction between the single dipole and its image charge (which is a manifestation of the cavity screening) is attractive and reduces the cost of a fluctuating dipolar moment.However, when we consider an extended distribution of dipoles (as realized by the bulk paraelectric), the result is very different.In particular, one can check that in the long-wavelength, longitudinally polarized q → 0 modes remain unaffected by the screening.To see this, consider an infinite line of longitudinallypolarized dipoles interacting with its image, which is also an infinite line with opposite in-plane component of the dipole moment.In contrast to the localized case, now the overall interaction is zero since the attractive tip-to-tail interaction for small transverse separations is cancelled by the interaction between distant dipoles, which are essentially tip-to-tip oriented.This naturally manifests in our calculations as the LO-TO splitting of the longwavelength modes remaining unchanged by the boundary conditions.Transverse modes on the other hand, do interact with the screening effect and this ends up most easily diagnosed by recasting the problem as determining ϵ(ω, r) rather than directly determining the dipole-dipole interactions, which become very complicated once quantum dynamics are included.Our result is also consistent with various recent no-go theorems pertaining to superradiance in cavities [70,71].
One may alternatively view this renormalization of the phonon stiffness as a realization of the concept of dynamical localization [25,26,51].This renormalization comes because when the phonon mode oscillates it also linearly couples to the electromagnetic vacuum, and thus the energy and inertia of that mode receive contributions also from the cloud of photons which are attached to the phonon.By using the cavity to suppress the electromagnetic field, we are essentially decoupling the phonon from its photonic cloud, and this is reflected in the evaluation of the correlation functions for the given geometry.What is somewhat unanticipated is that this cloud of photons actually makes the phonon mode "lighter," delocalizing the phonon coordinate.
We also can understand the purely quantum origin of the effect in this way.Cavity photons directly couple to the transverse component of the current, J ∼ ∂ t Q and by using the cavity one can change the radiative renormalization of the current fluctuations.This can only influence the actual phonon displacement ⟨Q 2 ⟩ via the canonical commutation relations, which couple the charge and current together.Therefore, observation of this experimentally would truly demonstrate the principle of "quantum control of quantum materials." As discussed earlier, in Sec.II C, we expect our theory to be valid not too close to the critical point, due to the local nature of the approximation we use for the phonon fluctuations.We now elaborate slightly on how this may breakdown, and what a more complete theory may look like.In particular, the key approximation we made is that the phonon correlation function is purely local, so that ⟨Q(x)Q(x ′ )⟩ ∼ δ(r−r ′ ).This then allowed us to integrate the phonons out in favor of a theory purely in terms of the local dielectric response, ϵ(iω m ).This approximation is motivated by the observation that since the photon group velocity is much larger than the phonon velocity, all phonon dispersion in the vicinity of the light-cone can be ignored.Naively, this is a very good approximation, but we do have reason to expect this to breakdown close to the critical point, since at the critical point lattice fluctuations become correlated over much larger length scales than the naive lattice estimate indicates, as one expects from the general theory of critical phenomena.
In order to properly accommodate these correlations, one must introduce phonon dispersion in to the model, with the simplest modification to our current theory given by the imaginary-time Lagrangian (using the same conventions we use in the rest of the paper) We have emphasize the new terms to be added to accommodate the spatial correlations of the phonon field, which by symmetry in an isotropic medium are characterized by two sound velocities; a transverse mode velocity v ⊥ and a longitudinal mode velocity v ∥ .Unfortunately, including these terms makes the solution technique we employ in this paper more difficult to apply, since it relied on the local nature of the phonon fluctuations.In principle, since this modification is only needed near the critical point anyways, it then also makes the mean-field treatment we use for the phonon frequency shift somewhat inapplicable as well.Instead, in order to proceed we propose employing a renormalization-group (RG) type procedure to handle this theory.This would then be able to quantitatively assess the degree to which the phonon correlations become important and potentially search for modifications to the critical phenomena due to the long-wavelength modifications due to the cavity geometry, though this is certainly beyond the scope of this paper.Physically though, it may be interesting to look for modifications to the lattice structure factor and diffuse x-ray scattering due to the presence of the cavity near the critical point, since this may reflect changes to the range of spatial correlations in the incipient lattice distortion.
We also comment on the relation between what we propose here and recent experiments which seem to indicate optically induced ferroelectricity in SrTiO 3 using strong terahertz driving [49,50].It seems likely that in both the present paper and in the terahertz driven experiments [49,50], it is important that there be strong precursor fluctuations of a polar soft-mode present in equilibrium.However, since a conclusive theoretical explanation for the induced ferroelectricity in these systems is still pending, it is difficult to draw a concrete connection.
In particular, owing to the larger magnitudes of electric field accessible in an optically driven material, nonlinear phononic effects may become important [40] which can complicate the picture by effectively inducing higher order couplings between the electric field and the softmode fluctuations.On the other hand, since the cavity fluctuations are much weaker in amplitude, it is expected that these will more selectively couple to the soft-mode directly, and do so through a predominantly linear polaritonic interaction.Therefore, it is hard to say with certainty how the mechanism of driven terahertz control is related to the mechanism of cavity control we outline here, other than the fact that both require preexisting fluctuations.
To summarize, we have developed an approach to studying local phonon fluctuations in a QPE interacting with a quantized cavity electromagnetic field.Rather than making a single-mode approximation or invoking the dipole approximation, we have included a complete continuum of modes which all couple locally to the QPE material.This allowed us to study the variation in the blue shift of the phonon modes in a spatially resolved way, and we found that near to the cavity walls the fluctuations increase due to the screening of the electric field by the cavity.This approach was then connected to the study of Casimir and Van-der Waals forces, both of which are manifestations of forces induced by electromagnetic fluctuations.
In the future it would be extremely interesting to extend this approach to include more complex heterostructures which feature multiple types material [72,73] including metals [74], insulators [75], superconductors [76,77], semiconductors [23], ferroelectrics [43][44][45][46][47], magnets [78][79][80][81][82][83][84], and multiferroics [85]-all of which are phases of matter which can be characterized by their couplings to electromagnetic fields and which are now realizeable down to the two-dimensional limit [86].In the future it would be very interesting to consider the mutual coupling of different stacked phases through their shared quantum electrodynamic environment.It is also interesting to extend these calculations into the ordered phase, where the order parameter as well as the fluctuations become important.Finally, probably the most important direction for future research are experimental realizations and confirmations of this theory.To that end, it is likely necessary to use more accurately obtained parameters and electromagnetic solvers.We therefore would envision interfacing this framework with ab initio calculations of microscopic parameters [87][88][89], as well as finite-element Maxwell-equation solvers.Conceptually, this is relatively straightforward but likely requires a large degree of technical work before it can be widely applied.
The first term is simply the free energy of the noninteracting ansatz, while the second term becomes ) This can be evaluated using Wick's theorem.We now vary the parameter Ω 2 T (r) to find the best approximation to the free energy.The evaluation is aided by performing a canonical transformation, as done in Appendix A. We shift Q = δQ − ηωm ω 2 m +Ω 2 T (r) A. This then allows us to express the free energy derivative as and the dielectric constant is found as We note as well that by direct examination, we have B9) which implies that the local TO frequency is a variational parameter which is conjugate to the local phonon fluctuations.
We can now prove that this is equivalent to the previous two approaches.We first use Wick's theorem to derive the quartic term in terms of the quadratic propagator.In the isotropic, paraelectric phase, we have ) where d is the spatial dimension.In the large d limit this is simply the Hartree term, with the Fock term being subleading in that limit.
We can now write our variational functional, using the expression for the interaction derived above and the relation between ⟨Q(x) 2 ⟩ and the derivative of W eff to get the exact result Now, we take a derivative with respect to the parameter Ω 2 T (r).After applying the chain rule, one finds T (r) is the local phonon fluctuation density.The first two terms cancel and the derivative of the phonon density is in general dependent on the TO frequency, so this leaves the variational equation simplified as We note that other than the factor of 2/d due to the Fock correction at finite d, this exactly matches our condition from the Hubbard-Stranovich method.
Appendix C: High-Temperature Behavior Here we provide another argument for why the transverse electromagnetic modes decouple at hightemperature.We begin with the quantum finitetemperature Lagrangian describing the phonons and their coupling to the electromagnetic field in a gauge independent form We can make the argument more clear if we do partial integration on the coupling term, so express it directly in terms of the gauge fields We choose the Coulomb gauge, where ∇•A = 0, such that A clearly couples to the transverse currents induced by the phonon modes, as these are the true radiative degrees of freedom.We can then separate the longitudinal and transverse parts of the radiation coupling as This clearly endows the LO phonon mode with the longrange Coulomb interaction.It also clearly couples to the phonon coordinate Q and makes perfect sense in the classical limit, where Q and ϕ become "time-independent" and the integral over Matsubara time becomes simply β 0 dτ → 1/T .On the other hand, the transverse modes couple via We note a number of important points.The first is that the radiative electric field has both retardation (due to the first term), and intrinsic dynamics due to the inductive response of the second term.Second, we see that the field couples to the transverse phonon current.In the classical limit high-temperature limit the retardation due to (∂ τ A) 2 goes away as the cost of a thermal tunneling event ∼ T becomes suppressed, and we are left with In this case the current J ⊥ = η∂ τ Q is essentially the canonical momentum associated to the dielectric polarization, but at high-temperatures this becomes uncorrelated with Q.Therefore, the photon field indeed dresses the expectation value for the phonon current fluctuations ⟨(J ⊥ ) 2 ⟩, but this now is independent of the phonon fluctuations ⟨Q 2 ⟩.In this way, we see that this is a quantum effect since the photons only couple to the phonon displacement through the commutation relations between the phonon displacement and current [Q, J] ̸ = 0, and at high-temperatures this vanishes.
trG(z, z) = sinh κ(L/2 − z) sinh κ(L/2 + z) sinh κL The quantity in brackets can be simplified as We therefore obtain the result for the electric-field fluctuations, which are weighted by an additional ω 2 m factor.This gives (D21) Therefore, the result ends up being relatively simple (we still do need to integrate over q and sum over ω m , to be clear).
(D22) This is therefore manifestly negative; now if we recall that the overall contribution goes as −⟨E 2 ⟩ due to the change in dielectric constant being inverse to the phonon frequency, we find that this will lead to a mode hardening at the boundary.
In total, we find that the final result including the Matsubara sum and momentum integral is This is to be numerically evaluated.We do need to be careful since the integrand is singular in small q, ω.We therefore numerically compute The first terms are the quantum corrections, which are evaluated by summing over the m = ±1, ±2, ... Matsubara frequencies, and this can be evaluated relatively easily numerically.The last term is the classical contribution which is tricky to evaluate due to the singular limit as ω → 0. the last term does indeed vanish in this limit.This is seen by first manipulating the integrand via x = κL and utilizing trig identities to get

FIG. 1 .
FIG.1.Schematic depiction of cavity.Two perfect metal plates are located in the xy plane at z = ±L/2, and the electromagnetic field and phonon modes coexist within the interior.The system has translational symmetry in the xy plane, which leads to a conserved in-plane momentum q, which can be taken along q ∥ êx, and obeys ideal metallic boundary conditions at z = ±L/2.Due to metallic boundary conditions, the photonic part of the wavefunction is suppressed near the boundary, leading to enhanced phonon fluctuations.

FIG. 5 .
FIG. 5. Perturbatively renormalized phonon frequency ΩT (z)including the electrodynamic correction due to cavity boundary conditions as a function of spatial coordinate z.System size is fized at L = 100µm, with coordinate z referenced from the midpoint, as shown in the inset.Temperature varies from high temperature of T = 100K (with little effect) down to low tmeperature of T = 10K, with a pronounced blue-shift in the local phonon frequency occurring at the boundary.