# Local Fluctuations in Cavity Control of Ferroelectricity

Jonathan B. Curtis<EMAIL_ADDRESS>College of Letters and Science,
University of California, Los Angeles, CA 90095, USA  Department of Physics,
Harvard University, Cambridge, MA 02138, USA  Marios H. Michael  Max Planck
Institute for the Structure and Dynamics of Matter, 22761 Hamburg, Germany
Eugene Demler  Institute for Theoretical Physics, ETH Zurich, 8093 Zurich,
Switzerland

###### Abstract

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.

###### pacs:

##  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 [16, 17, 18, 6],
excitonic insulators [19], antiferromagnets [20, 21], spin-liquids [22],
quantum Hall fluids [23], and ferroelectrics [24, 25, 26, 27, 28], with a
great deal more proposed to exhibit strong coupling between material and
optical excitations [29]. Recent experiments on the metal-insulator transition
in 1$T$-TaS2 even seem to have seen promising signatures of cavity control on
the transition temperature [30]. Experiments have also seen fascinating
phenomena occur when unconventional superconductors are strongly coupled to
the quantum electromagnetic bath [31, 32].

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 [33, 34,
35]. Furthermore, there are a number of appealing candidate systems, such as
SrTiO3 [36, 37, 38, 39, 40] and various moiré and Van der Waals materials [41,
42, 43, 44, 45] which may be suitable for proof-of-principle experiments.
Intrinsic SrTiO3 is believed a quantum paraelectric (QPE) [36, 37, 38, 39,
40], 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 [37], and recently resonant optical excitation of the lattice
[46, 38] have also been shown to seemingly induce a transition into the
ordered phase [47, 48], making this a prime candidate for demonstrating cavity
control over the phase diagram [24, 25, 27]. 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 [49].

In this work, we extend our analysis of fluctuating quantum paraelectric to a
fully multimode, 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 few-mode calculations. This surprising result then has
important implications for future experiments on cavity control of
ferroelectricity and potentially other phases of matter [30, 31, 32].

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 frequency-dependent electric-field
correlation function. The behavior of this correlations function is very well
studied, dating back to seminal work on the Casimir force [50, 51, 52, 53] 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 [54]. Furthermore, we argue our technique can be
easily extended to include phonon loss, anisotropy, and mode couplings
provided the dielectric constant dispersion $\epsilon(\omega)$ 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 two-dimensional systems via characterization in terms of the
reflection coefficients at interfaces [53]. 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
$\langle\mathbf{Q}^{2}(x)\rangle$ increase, resulting in a blue-shift of the
soft-mode frequency due to the anharmonic coupling of the phonon mode. This
leads to an overall thickness dependence which may be pronounced for thinner
cavities and results in a diminished low-frequency dielectric constant
$\epsilon(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 temperatures, indicating the origin of this
effect is of a truly quantum origin, and should drop off once
$T\gtrsim\hbar\Omega_{T}/K_{B}$. For materials such as SrTiO3, with
$\Omega_{T}\sim 2$THz, this provides an important ceiling on the temperature
of experiments which is roughly of order 100K.

Figure 1: Schematic depiction of cavity. Two perfect metal plates are located
in the $xy$ plane at $z=\pm L/2$, and the electromagnetic field and phonon
modes coexist within the interior. The system has translational symmetry in
the $xy$ plane leading to momentum which can be taken along
$\mathbf{q}\parallel\hat{\mathbf{e}}_{x}$, and obeys ideal metallic boundary
conditions at $z=\pm L/2$. Due to metallic boundary conditions, the photonic
part of the wavefunction is suppressed near the boundary, leading to enhanced
phonon fluctuations.

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.1, followed by a reformulation on the Matsubara axis for finite-
temperature calculations in Sec. II.2. 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.1
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 localized, isotropic, polar
phonon mode. In particular, since the phonon group velocity is much slower
than the photon, the local phonon approximation 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 $\mathbf{Q}(\mathbf{r})$
with TO mode frequency $\Omega_{T}$ and effective charge $\eta$ coupled to the
electromagnetic field, described by $\mathbf{E}$ and $\mathbf{B}$. This system
is most simply described in terms of a Lagrangian which generates the
equations of motion; quantum mechanically these enter into the appropriate
partition function. We have

$\mathcal{L}=\mathcal{L}_{\rm EM}+\mathcal{L}_{\rm ph}+\mathcal{L}_{\rm int}.$
(1)

The Maxwell Lagrangian is

$\mathcal{L}_{\rm EM}=\frac{1}{2}\mathbf{E}^{2}-\frac{1}{2}\mathbf{B}^{2}.$
(2)

In terms of the gauge potentials, the electromagnetic fields are expressed as

$\displaystyle\mathbf{E}=-\nabla A_{0}-\partial_{t}\mathbf{A}$ (3a)
$\displaystyle\mathbf{B}=\nabla\times\mathbf{A}.$ (3b)

Here and throughout we use units where $\hbar=c=k_{B}=\epsilon_{0}=1$.

The phonon Lagrangian, in the absence of dispersion, is purely local and is
simply given by:

$\mathcal{L}_{\rm
ph}=\frac{1}{2}\left[\left(\partial_{t}\mathbf{Q}\right)^{2}-\Omega_{0}^{2}\mathbf{Q}^{2}\right]-\frac{\lambda}{4}\left(\mathbf{Q}^{2}\right)^{2},$
(4)

where $\Omega_{0}$ is the bare tranverse optical (TO) mode frequency, and
$\lambda$ is the phonon-phonon interaction strength. Due to symmetry, these
are the only terms allowed at ${\bf q}=0$ 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 $\bf Q$ generates a polarization which
then couples to $\bf E$, via

$\mathcal{L}_{\rm int}=+\eta\mathbf{Q}\cdot\left[-\nabla
A_{0}-\partial_{t}\mathbf{A}\right].$ (5)

Here, $\eta$ is the light-matter interaction constant and it sets, among other
things, the size of the splitting between the longitudinal optical and
transverse optical (LOTO) frequency splitting due to the Coulomb part of the
electromagnetic interaction. 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

$\Omega_{T}^{2}=\Omega_{0}^{2}+\lambda\langle\mathbf{Q}^{2}(\mathbf{r})\rangle.$
(6)

In general, we allow for spatially varying fluctuations of $\mathbf{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, $\langle\mathbf{Q}^{2}(\mathbf{r})\rangle$. We do this by the
familiar linear-response formalism, obtaining the equilibrium fluctuations of
$\mathbf{Q}$ by solving for the causal response to an external perturbation.
We thus introduce a source field $\mathbf{F}(\mathbf{r},t)$ which couples to
the phonon mode via

$\mathcal{L}_{\rm source}=\mathbf{F}(x)\cdot\mathbf{Q}(x),$ (7)

such that

$\mathbb{\hat{D}}^{R}(x,x^{\prime})=-\frac{\delta\langle\mathbf{Q}(x)\rangle}{\delta\mathbf{F}(x^{\prime})}\bigg{|}_{\mathbf{F}=0}=-i\theta(t-t^{\prime})\langle[\mathbf{Q}(x),\mathbf{Q}(x^{\prime})]\rangle$
(8)

for causal response function. From this we can apply the fluctuation-
dissipation relation to obtain

$\langle\mathbf{Q}(x)\cdot\mathbf{Q}(x)\rangle=-\int\frac{d\omega}{2\pi}\coth\frac{\beta\omega}{2}\Im\left[\mathbf{tr}\mathbb{\hat{D}}^{R}(\mathbf{r},\mathbf{r};\omega)\right].$
(9)

This then allows to characterize the local density of phonon fluctuations.

###  II.1 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

$\displaystyle\left[-\omega^{2}+\Omega_{T}^{2}\right]\mathbf{Q}(\mathbf{r},\omega)=\eta\mathbf{E}(\mathbf{r},\omega)+\mathbf{F}(\mathbf{r},\omega)$
(10a)
$\displaystyle+i\omega\mathbf{B}(\mathbf{r},\omega)=\nabla\times\mathbf{E}(\mathbf{r},\omega)$
(10b)
$\displaystyle-i\omega\left[\mathbf{E}(\mathbf{r},\omega)+\eta\mathbf{Q}(\mathbf{r},\omega)\right]=\nabla\times\mathbf{B}(\mathbf{r},\omega)$
(10c)

We compute the response function of $\mathbf{Q}$ as follows. Let us first
introduce the bare response function for $\mathbf{Q}(\mathbf{r},\omega)$ of

$\chi_{0}(\omega)=\frac{1}{-\omega^{2}+\Omega_{T}^{2}},$ (11)

such that we can obtain the response of $\mathbf{Q}$ as

$\mathbf{Q}(\mathbf{r},\omega)=\chi_{0}(\omega)\left[\eta\mathbf{E}(\mathbf{r},\omega)+\mathbf{F}(\mathbf{r},\omega)\right].$
(12)

Our job is not done because we need the response of the phonon not to the
total force, which is $\eta\mathbf{E}+\mathbf{F}$, but only to the external
force $\mathbf{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

$\displaystyle+i\omega\mathbf{B}(\mathbf{r},\omega)=\nabla\times\mathbf{E}(\mathbf{r},\omega)$
(13a)
$\displaystyle-i\omega\left(\mathbf{E}(\mathbf{r},\omega)+\eta\chi_{0}(\omega)\left[\eta\mathbf{E}(\mathbf{r},\omega)+\mathbf{F}(\mathbf{r},\omega)\right]\right)=\nabla\times\mathbf{B}(\mathbf{r},\omega)$
(13b)

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
$\mathbf{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

$\omega^{2}\epsilon(\omega)\mathbf{E}(\mathbf{r},\omega)-\nabla\times\nabla\times\mathbf{E}(\mathbf{r},\omega)=-\omega^{2}\eta\chi_{0}(\omega)\mathbf{F}(\mathbf{r},\omega).$
(14)

Here we have introduced the dielectric constant

$\epsilon(\omega)=1+\eta^{2}\chi_{0}(\omega).$ (15)

We now can use this to eliminate the electric field formally as

$\mathbf{E}(\mathbf{r},\omega)=\int
d^{3}r^{\prime}\mathbb{\hat{G}}(\mathbf{r},\mathbf{r}^{\prime};\omega)\left[-\eta\omega^{2}\chi_{0}(\omega)\mathbf{F}(\mathbf{r}^{\prime},\omega)\right]$
(16)

where

$\mathbb{\hat{G}}(\mathbf{r},\mathbf{r}^{\prime};\omega)=\left[\omega^{2}\epsilon(\omega)\mathds{1}-\nabla\times\nabla\times\right]^{-1}.$
(17)

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,

$\mathbf{Q}(\mathbf{r},\omega)=\int
d^{3}r^{\prime}\mathbb{\hat{D}}(\mathbf{r},\mathbf{r}^{\prime};\omega)\mathbf{F}(\mathbf{r}^{\prime},\omega),$
(18)

is given as a function of the photon Green’s function:

$\mathbb{\hat{D}}^{R}(\mathbf{r},\mathbf{r}^{\prime};\omega)=\chi_{0}(\omega)\left[\delta^{3}(\mathbf{r}-\mathbf{r}^{\prime})\mathds{1}-\eta^{2}\omega^{2}\chi_{0}(\omega)\mathbb{\hat{G}}^{R}(\mathbf{r},\mathbf{r}^{\prime};\omega)\right].$
(19)

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
$\epsilon(\omega)$, and therefore is sensitive to the cavity geometry.

We note we can write this in an elegant way by using

$-\eta^{2}\left(\chi_{0}(\omega)\right)^{2}=\frac{\delta\epsilon(\omega)}{\delta\Omega_{T}^{2}}$
(20)

to obtain

$\begin{split}\mathbb{\hat{D}}^{R}(\mathbf{r},\mathbf{r}^{\prime};\omega)=&\chi_{0}(\omega)\delta^{3}(\mathbf{r}-\mathbf{r}^{\prime})\mathds{1}+\\\
&\omega^{2}\frac{\delta\epsilon(\omega)}{\delta\Omega_{T}^{2}}\mathbb{\hat{G}}^{R}(\mathbf{r},\mathbf{r}^{\prime};\omega).\end{split}$
(21)

###  II.2 Matsubara Formalism

The result presented is derived using a real-time linear-response formalism,
which is the most transparent presentation. In Appendix A, we equivalently
derive this result using the equilibrium Matsubara frequency representation,
and in Appendix B we present a third way of deriving this result based on a
variational procedure for the Matsubara free-energy functional.

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 details, we refer the reader to
the relevant Appendix A. 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 $\omega_{m}=2\pi mT$ as

$\epsilon(\omega_{m};\mathbf{r})=1+\frac{\eta^{2}}{\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r})},$
(22)

allowing for a locally varying TO mode frequency. Likewise, the unscreened
phonon propagator becomes

$\mathscr{\hat{D}}_{0}(\mathbf{r},\mathbf{r}^{\prime};\omega_{m})=\frac{1}{\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r})}\delta^{3}(\mathbf{r}-\mathbf{r}^{\prime}),$
(23)

and the photon propagator becomes

$\mathscr{\hat{G}}(\mathbf{r},\mathbf{r}^{\prime};\omega_{m})=\left[\omega_{m}^{2}\epsilon(\omega_{m},\mathbf{r})+\nabla\times\nabla\times\right]^{-1}.$
(24)

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

$\langle\mathbf{Q}^{2}(\mathbf{r},t)\rangle=T\sum_{\omega_{m}}\\\
\mathbf{tr}\left[\mathscr{\hat{D}}_{0}(\mathbf{r},\mathbf{r};\omega_{m})+\omega_{m}^{2}\frac{\partial\epsilon(\omega_{m},\mathbf{r})}{\partial\Omega_{T}^{2}}\mathscr{\hat{G}}(\mathbf{r},\mathbf{r};\omega_{m})\right]$
(25)

Note that

$\frac{\partial\epsilon(\omega_{m},\mathbf{r})}{\partial\Omega_{T}^{2}}=-\frac{\eta^{2}}{(\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r}))^{2}}$
(26)

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. [27], 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

$\delta\Omega_{T}^{2}(\mathbf{r})=\lambda\left[\langle\mathbf{Q}^{2}(\mathbf{r})\rangle-\langle\mathbf{Q}^{2}\rangle_{\rm
bulk}\right].$ (27)

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 $\lambda$, the momentum space cutoff
$\Lambda$, and other relevant parameters.

###  II.3 Parameters

Following Ref. [46], we first extract the value of the anharmonic potential
$U^{(2)}$ as a function of the ionic displacement coordinate $\mathbf{u}$, for
the case of SrTiO3, estimated from spectroscopy to give $\nu_{T}\sim
35$\mathrm{T}\mathrm{Hz}^{2}\mathrm{/}\mathrm{p}\mathrm{m}^{2}$$. This is
obtained from nonlinear terahertz spectroscopy by measuring
$\Omega_{T}^{2}(\mathbf{u})\sim\Omega_{T}^{2}+\nu_{T}\mathbf{u}^{2}+O(\mathbf{u}^{4})$.
However, this is not directly related to the long-wavelength order parameter
111This is similar to the case of Ginzburg-Landau theory where, in general the
microscopic parameters such as density-of-states and quasiparticle gap are not
directly related to the long-wavelength order parameter. Rather, only certain
ratios of coupling constants are fixed in terms of microscopic parameters.,
which in general involves a complicated coarse-graining 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
$\lambda\sim\nu_{T}V_{\rm eff}/M_{\rm eff}$ where $M_{\rm eff}$ is an
effective mass scale and $V_{\rm 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 $M_{\rm eff}\sim M_{\rm
Ti}\sim 50$\mathrm{u}$$.

Finally, in the spirit of a real-space renormalization group procedure, we
anticipate that the effective volume $V_{\rm eff}$ which appears in the
relation between the order-parameter and the microscopic degrees of freedom
should be related to the momentum space cutoff we place on our model,
$\Lambda\sim 1/a$ with $a$ the size of the coarse-grained “blocks.” As a
result, we estimate that $V_{\rm eff}\sim 1/\Lambda^{3}\sim a^{3}$. We thus
identify the long-wavelength coupling as $\lambda=a^{3}\nu_{T}/M_{\rm Ti}$.
This leads to the cutoff-dependent estimate of

$\lambda=32\times 10^{4}$\mathrm{T}\mathrm{Hz}^{3}$a^{3},$ (28)

which we use in this work.

Finally, by ignoring the gradient terms $\sim\nabla\mathbf{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
$\mathbf{Q}(\mathbf{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\sim
50$\mathrm{\SIUnitSymbolMicro m}$$ or larger; as a result, we will consider a
cutoff of $a=5$\mathrm{\SIUnitSymbolMicro m}$$.

In the future it would be particularly interesting to consider going even
closer to the critical point, so that $a\sim 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_{\rm
QED}^{3}\sim(1$\mathrm{THz}$/c)^{3}\sim(.3$\mathrm{mm}$)^{-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}\sim q_{\rm 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.

H

Finally, we must fix the phonon TO and LO mode frequencies. For the bulk TO
frequency $\Omega_{T}$ and LO-TO splitting $\sim\eta$, we take
$\Omega_{T}=.5$THz to emulate the ferroelectric soft mode, and $\eta=10$THz,
such that $\epsilon(0)\sim\eta^{2}/\Omega_{T}^{2}\sim 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
$\omega_{c}=40$THz, 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

$\nabla\times\nabla\to-\mathbf{q}\times\mathbf{q}\times=\mathbf{q}^{2}\mathds{1}-\mathbf{q}\otimes\mathbf{q}$
(29)

we find the electromagnetic correlation function splits into two decoupled
subspaces: the transverse modes and longitudinal modes, with

$\omega_{m}^{2}\mathscr{\hat{G}}(\omega_{m},\mathbf{q})=\\\
\frac{\mathbf{q}\otimes\mathbf{q}}{\mathbf{q}^{2}}\frac{1}{\epsilon(\omega_{m})}+\left(\mathds{1}-\frac{\mathbf{q}\otimes\mathbf{q}}{\mathbf{q}^{2}}\right)\frac{\omega_{m}^{2}}{\omega^{2}_{m}\epsilon(\omega_{m})+\mathbf{q}^{2}}.$
(30)

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

$\mathscr{\hat{D}}_{0}(\omega_{m},\mathbf{q})=\frac{1}{\omega_{m}^{2}+\Omega_{T}^{2}}.$
(31)

We find for the overall result, with UV cutoff on momentum $\Lambda$,

$\langle\mathbf{Q}^{2}\rangle=T\sum_{\omega_{m}}\int_{|{\bf q}|<\Lambda}\\\
\left[\frac{3}{\omega_{m}^{2}+\Omega_{T}^{2}}+\frac{\partial\epsilon(\omega_{m})}{\partial\Omega_{T}^{2}}\left(\frac{1}{\epsilon(\omega_{m})}+\frac{2\omega_{m}^{2}}{\omega_{m}^{2}\epsilon(\omega_{m})+\mathbf{q}^{2}}\right)\right].$
(32)

This simplifies into one longitudinal mode and $d-1$ transverse modes, with
the longitudinal modes contributing (per momentum space mode)

$T\sum_{\omega_{m}}\left[\frac{1}{\omega_{m}^{2}+\Omega_{T}^{2}}+\frac{\partial\epsilon(\omega_{m})}{\partial\Omega_{T}^{2}}\frac{1}{\epsilon(\omega_{m})}\right]\\\
=\frac{\coth(\beta\Omega_{L}/2)}{2\Omega_{L}}.$ (33)

The momentum spacec integral will then simply yield a factor of
$\Lambda^{3}/6\pi^{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

$\langle\mathbf{Q}^{2}\rangle=\int_{|{\bf
q}|<\Lambda}\bigg{[}\frac{\coth(\beta\Omega_{L}/2)}{2\Omega_{L}}\\\
+2\times\frac{\coth(\beta\Omega_{+,q}/2)}{2\Omega_{+,q}}\frac{\Omega_{+,q}^{2}-q^{2}}{\Omega_{+,q}^{2}-\Omega_{-,q}^{2}}\\\
+2\times\frac{\coth(\beta\Omega_{-,q}/2)}{2\Omega_{-,q}}\frac{\Omega_{-,q}^{2}-q^{2}}{\Omega_{-,q}^{2}-\Omega_{+,q}^{2}}\bigg{]}.$
(34)

The dispersion of the upper polariton branch ($\Omega_{+,q}$) and lower
polariton branch ($\Omega_{-,q}$) are found to be

$\Omega_{\pm,q}=\sqrt{\frac{\Omega_{L}^{2}+q^{2}}{2}\pm\sqrt{\left(\frac{\Omega_{L}^{2}+q^{2}}{2}\right)^{2}-\Omega_{T}^{2}q^{2}}}.$
(35)

This is in accordance with the result one would expect from Bogoliubov
transformation and diagonalizing the quasiparticle Hamiltonian.

###  III.1 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

$\langle\mathbf{Q}^{2}\rangle=T\int_{|{\bf
q}|<\Lambda}\bigg{[}\frac{1}{\Omega_{L}^{2}}+\frac{2}{\Omega_{T}^{2}}\bigg{]}.$
(36)

This is exact and independent of momentum $\mathbf{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
$\omega_{m}^{2}/q^{2}$ at low frequencies, and therefore vanish in the static
$\omega_{m}=0$ 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
$\langle\mathbf{Q}^{2}\rangle$ diminish with increasing temperature, but
rather that the contribution from the electromagnetic field falls off. As a
result, if one performs a high-temperature expansion of
$\langle\mathbf{Q}^{2}\rangle$ in the bulk one would get
$\langle\mathbf{Q}^{2}\rangle_{\rm bulk}\sim a_{\rm bulk}T+b_{\rm
bulk}/T+...$, with the leading term the classical equipartition result and the
subleading terms coming from the high-temperature expansion of
$\coth(\omega\beta/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
$\langle\mathbf{Q}^{2}\rangle_{\rm surf}\sim a_{\rm surf}T+b_{\rm surf}/T+...$
and find that $a_{\rm bulk}=a_{\rm 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.

Figure 2: Schematic depiction of temperature dependence of phonon fluctuations
$\langle\mathbf{Q}^{2}\rangle$ 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.

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 can be
understood more technically in Appendix C; here we provide a simple heuristic.
Essentially, photons don’t actually couple directly to the phonon displacement
$\mathbf{Q}$, but rather couple to the displacement current
$\mathbf{J}\sim\partial_{t}\mathbf{Q}$ associated to transverse oscillations
in the charge distribution. This current is canonically conjugate to the
object of interest, $\mathbf{Q}$, via
$[\mathbf{Q}(\mathbf{r}),\mathbf{J}(\mathbf{r}^{\prime})]\propto
i\hbar\delta^{3}(\mathbf{r}-\mathbf{r}^{\prime})$ and therefore in the quantum
system (at low temperatures) they are not independently fluctuating. Thus,
modifying the fluctuations of the current $\mathbf{J}$ can in turn induce
changes in the fluctuations of $\mathbf{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 $\mathbf{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
$\rho\sim\nabla\cdot\mathbf{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) [51, 52, 53]. We can take the momentum to lie in the $x$-direction,
such that the Green’s function for the gauge potential $\mathscr{G}$ reads

$\left[\omega_{m}^{2}\epsilon(\omega_{m})\mathds{1}+\begin{pmatrix}-\partial_{z}^{2}&0&iq\partial_{z}\\\
0&q^{2}-\partial_{z}^{2}&0\\\ iq\partial_{z}&0&q^{2}\\\
\end{pmatrix}\right]\mathscr{G}(z,z^{\prime})=\mathds{1}\delta(z-z^{\prime}).$
(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 $\mathbf{E}$ and normal component of $\bf B$. We take the cavity plates to
be located at $z=\pm 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 $\mathbf{tr}\mathscr{G}(z,z)$ is given in closed form as

$\mathbf{tr}\mathscr{G}(z,z)=\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{2\kappa\sinh\kappa
L}\\\ +\frac{1}{\omega_{m}^{2}\epsilon(\omega_{m})}\delta(0).$ (38)

Here $\kappa=\sqrt{\omega_{m}^{2}\epsilon(\omega_{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 $\delta(0)$,
which is understood as $\lim_{z^{\prime}\to z}\delta(z-z^{\prime})$. This
quantity is in fact independent of system-size and geometry and thus ends up
getting renormalized away by the counter-term $\Omega_{0}^{2}$. In particular,
we only compare this result to result obtained for $L\to\infty$ (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

$\left(\Delta\Omega_{T}(z)\right)_{\rm cav}^{2}=\lambda
T\sum_{\omega_{m}}\int^{\Lambda}\frac{d^{2}q_{\parallel}}{(2\pi)^{2}}\frac{\partial\epsilon(\omega_{m})}{\partial\Omega_{T}^{2}}\frac{\omega_{m}^{2}}{2\kappa}\left[\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}-\frac{1}{2}\right].$ (39)

We note this expression does still depend on the cutoff $\Lambda=\pi/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
$\omega_{m}\to 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 $\omega_{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
$\omega_{m}^{2}$ which then vanishes quadratically at small frequency. On the
other hand, the denominator involves $\kappa\sinh\kappa L$, which for
$q_{\parallel}\to 0$ also vanishes quadratically as
$|\omega_{m}|\sinh|\omega_{m}|\sqrt{\epsilon(0)}L$.

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.1 and shown in greater detail in App. 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\kappa(L/2-z)\sinh\kappa(L/2+z)/\sinh\kappa L\leq\frac{1}{2}$. As a
result, we find the frequency shift is always positive, with
$\Delta\Omega_{T}^{2}(z)\geq 0$, such that the cavity in fact suppresses the
ferroelectric order. We will comment on this more later.

Figure 3: First order correction to phonon frequency due to cavity boundary
conditions at $z=\pm L/2$. We here fix the temperature to $T=10K$, which is
lower the coherence temperature for the phonon frequency of $\Omega_{T}=.5$THz
(we note the coherence frequency involves a factor $2\pi$ so that $2\pi
T\lesssim\Omega_{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.

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.

Figure 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$\mathrm{\SIUnitSymbolMicro m}$$ (a), and a larger system with
$L=70$\mathrm{\SIUnitSymbolMicro 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.

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 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\to\infty$. For smaller systems however, the shift may become sizeable at
the cavity midpoint.

Figure 5: Perturbatively renormalized phonon frequency $\Omega_{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$\mathrm{\SIUnitSymbolMicro m}$$, with coordinate $z$ referenced from
the midpoint, as shown in the inset. Temperature varies from high temperature
of $T=100$K (with little effect) down to low tmeperature of $T=10$K, with a
pronounced blue-shift in the local phonon frequency occurring at the boundary.

These results are best understood by simply looking at the spatial profile for
the renormalized frequency $\Omega_{T}^{2}(z)$, depicted in Fig. 5. For a
cavity size of $100$\mathrm{\SIUnitSymbolMicro 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\sim 0$). This is inline with the previous
arguments we can given, namely that the cavity renormalization is (i) a purely
quantum effect setting in once $2\pi T\sim\Omega_{T}$, and (ii) that it is a
surface effect in response to the boundary conditions at $z=\pm L/2$. We also
always see that the frequency blue-shifts—this is in some tension with a
number of previous theoretical investigations [24, 25, 27, 28], 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.

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 [50, 51, 52, 53]. 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 $\mathscr{G}$ regardless of system size. This assumption leads to the
unphysical result that number of electromagnetic modes per unit volume $\sim
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 222The authors would like to especially
acknowledge private correspondence with Yuto Ashida on this point, who
patiently and carefully worked through this subtlety with us, and who
ultimately helped uncover the origin of the discrepancies. Discussions with
Ataç İmamoğlu were also valuable in this regard..

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
$\varepsilon(\omega=0)\sim\eta^{2}/\Omega_{T}^{2}$ 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 “dead-layer” effect [57]. Canonically,
this effect was found in exactly this sort of system, comprised of a thin
layer of SrTiO3 sandwiched between two electrodes to form a capacitor [58].
Historically, the dead-layer effect has been attributed largely to strain
effects [59] between the two interfaces which also acts to blue-shift the soft
mode, however it seems the matter has not been entirely settled [60]. 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 [61] have already shown promising results [62, 63, 64].
Theoretically, this would also warrant further calculations which have open
(air) boundary conditions for the SrTiO3 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 [50, 52, 53]. This is most apparent if
one returns to Eqn. 39, which closely mirrors the results for Casimir forces
(see, e.g. Ref. [53]). 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 [39]. It may be investigate
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 [50, 51, 52]. 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 $\mathbf{q}\to 0$ modes remain unaffected by the
screening. To see this, consider an infinite line of longitudinally-polarized
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 long-
wavelength 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
$\epsilon(\omega,\mathbf{r})$ rather than directly determining the dipole-
dipole interactions, which become very complicated once quantum dynamics are
included.

One may alternatively view this renormalization of the phonon stiffness as a
realization of the concept of dynamical localization [25, 24, 49]. 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,
$\mathbf{J}\sim\partial_{t}\mathbf{Q}$ and by using the cavity one can change
the radiative renormalization of the current fluctuations. This can only
influence the actual phonon displacement $\langle\mathbf{Q}^{2}\rangle$ 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.”

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
[65, 66] including metals [67], insulators [68], superconductors [69, 70],
ferroelectrics [41, 42, 43, 44, 45], magnets [71, 72, 73, 74, 75, 76, 77], and
multiferroics [78]—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 [79]. 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 [80, 81, 82], 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.

###### Acknowledgements.

We would like to acknowledge invaluable discussions with Yuto Ashida, Ankit
Disa, Urs Staub, Prineha Narang, Ataç İmamoğlu, David Hsieh, Pavel Dolgirev,
Aaron Lindenberg, N. Peter Armitage, John Sous, Maximilian Daschner, Andrea
Cavalleri, Jérôme Faist, Dieter Jaksch, Misha Fogler, Ilya Esterlis, and John
Philbin. J.B.C. is an HQI Prize Postdoctoral Fellow and gratefully
acknowledges support from the Harvard Quantum Initiative. E.D. and J.B.C. were
supported by Harvard-MIT CUA, AFOSR-MURI: Photonic Quantum Matter Award No.
FA95501610323, the ARO grant Control of Many-Body States Using Strong Coherent
Light-Matter Coupling in Terahertz Cavities.

## References

  * Basov  _et al._ [2017] D. Basov, R. Averitt, and D. Hsieh, Towards properties on demand in quantum materials, Nat. Mat. 16, 1077 (2017).
  * Oka and Aoki [2009] T. Oka and H. Aoki, Photovoltaic Hall effect in graphene, Phys. Rev. B 79, 081406(R) (2009).
  * Kitagawa  _et al._ [2010] T. Kitagawa, E. Berg, M. Rudner, and E. Demler, Topological characterization of periodically driven quantum systems, Phys. Rev. B 235114, 82 (2010).
  * Lindner  _et al._ [2011] N. Lindner, G. Refael, and V. Galitski, Floquet topological insulator in semiconductor quantum wells, Nat. Phys. 7, 490 (2011).
  * Cavalleri [2018] A. Cavalleri, Photo-induced superconductivity, Contemp. Phys. 59, 31 (2018).
  * Gao  _et al._ [2020] H. Gao, F. Schlawin, M. Buzzi, A. Cavalleri, and D. Jaksch, Photoinduced Electron Pairing in a Driven Cavity, Phys. Rev. Lett. 125, 053602 (2020).
  * Michael  _et al._ [2022] M. H. Michael, M. Först, D. Nicoletti, S. R. U. Haque, Y. Zhang, A. Cavalleri, R. D. Averitt, D. Podolsky, and E. Demler, Generalized fresnel-floquet equations for driven quantum materials, Phys. Rev. B 105, 174301 (2022).
  * Michael  _et al._ [2020] M. H. Michael, A. von Hoegen, M. Fechner, M. Först, A. Cavalleri, and E. Demler, Parametric resonance of Josephson plasma waves: A theory for optically amplified interlayer superconductivity in YBa2Cu3O6+x, Phys. Rev. B 102, 174505 (2020).
  * Dolgirev  _et al._ [2020] P. E. Dolgirev, M. H. Michael, A. Zong, N. Gedik, and E. Demler, Self-similar dynamics of order parameter fluctuations in pump-probe experiments, Phys. Rev. B 101, 174306 (2020).
  * Dolgirev  _et al._ [2021] P. E. Dolgirev, A. Zong, M. H. Michael, J. B. Curtis, D. Podolsky, A. Cavalleri, and E. Demler, Periodic dynamics in superconductors induced by an impulsive optical quench (2021), arXiv:2104.07181 .
  * Baumann  _et al._ [2010] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Dicke quantum phase transition with a superfluid gas in an optical cavity, Nature 464, 1301 (2010).
  * Babadi  _et al._ [2017] M. Babadi, K. M., I. Martin, G. Refael, and E. Demler, Theory of parametrically amplified electron-phonon superconductivity, Phys. Rev. B 96, 014512 (2017).
  * Chaudhary  _et al._ [2020] S. Chaudhary, A. Haim, Y. Peng, and G. Refael, Phonon-induced Floquet topological phases protected by space-time symmetries, Phys. Rev. Res. 2, 043431 (2020).
  * Bloch  _et al._ [2022] J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Strongly correlated electron–photon systems, P.N.A.S. 118, e2105618118 (2022).
  * Schlawin  _et al._ [2022] F. Schlawin, D. Kennes, and M. Sentef, Cavity quantum materials, Appl. Phys. Rev. 9, 011312 (2022).
  * Sentef  _et al._ [2018] M. Sentef, M. Ruggenthaler, and A. Rubio, Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity, Sci. Adv. 4, eaau6969 (2018).
  * Curtis  _et al._ [2019] J. Curtis, Z. Raines, A. Allocca, M. Hafezi, and V. Galitski, Cavity Quantum Eliashberg Enhancement of Superconductivity, Phys. Rev. Lett. 122, 167002 (2019).
  * Schlawin  _et al._ [2019] F. Schlawin, A. Cavalleri, and D. Jaksch, Cavity-Mediated Electron-Photon Superconductivity, Phys. Rev. Lett. 122, 133602 (2019).
  * Mazza and Georges [2019] G. Mazza and A. Georges, Superradiant Quantum Materials, Phys. Rev. Lett. 122, 017401 (2019).
  * Curtis  _et al._ [2022] J. Curtis, A. Grankin, N. Poniatowski, V. Galitski, P. Narang, and E. Demler, Cavity magnon-polaritons in cuprate parent compounds, Phys. Rev. Res. 4, 013101 (2022).
  * Parvini  _et al._ [2021] T. Parvini, V. Bittencourt, and S. Kusminskiy, Antiferromagnetic cavity optomagnonics, Phys. Rev. Res. 2, 022027(R) (2021).
  * Chiochetta  _et al._ [2021] A. Chiochetta, D. Kiese, C. Zelle, F. Piazza, and S. Diehl, Cavity-induced quantum spin liquids, Nat. Comm. 12, 5901 (2021).
  * Paravicini-Bagliani  _et al._ [2019] G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist, Magneto-transport controlled by Landau polariton states, Nature Physics 15, 186 (2019), 1805.00846 .
  * Ashida  _et al._ [2020] Y. Ashida, A. İmamoğlu, J. Faist, D. Jaksch, A. Cavalleri, and E. Demler, Quantum Electrodynamic Control of Matter: Cavity-Enhanced Ferroelectric Phase Transition, Physical Review X 10, 041027 (2020), 2003.13695 .
  * Latini  _et al._ [2021] S. Latini, D. Shin, S. Sato, C. Schäfer, U. De Giovanni, H. Hübener, and A. Rubio, The ferroelectric photo ground state of SrTiO3: Cavity materials engineering, P.N.A.S. 118, e2105618118 (2021).
  * Zhang  _et al._ [2019] Q. Zhang, Z. Zhen, Y. Yang, G. Gan, D. Jariwala, and X. Cui, Hybrid phonon-polaritons at atomically-thin van der Waals heterointerfaces for infrared optical modulation, Optics Express 27, 18585 (2019).
  * Lenk  _et al._ [2022a] K. Lenk, J. Li, P. Werner, and M. Eckstein, Collective theory for an interacting solid in a single-mode cavity (2022a), arXiv:2205.05559 .
  * Lenk  _et al._ [2022b] K. Lenk, J. Li, and M. Eckstein, Dynamical mean-field study of a photon-mediated ferroelectric phase transition (2022b), arXiv:2209.04246 .
  * Basov  _et al._ [2020] D. Basov, A. Asenjo-Garcia, P. Schuck, X. Zhu, and A. Rubio, Polariton panorama, Nanophotonics 10, 549 (2020).
  * Jarc  _et al._ [2022] G. Jarc, S. Mathengattil, A. Montanaro, F. Giusti, E. Rigoni, F. Fassioli, S. Winnerl, S. Zilio, D. Mihailovic, P. Prelovšek, M. Eckstein, and D. Fausti, Cavity control of the metal-to-insulator transition in 1T-TaS2 (2022), arXiv:2210.02346 .
  * Thomas  _et al._ [2019] A. Thomas, E. Devaux, K. Nagarajan, T. Chervy, M. Seidel, D. Hagenmüller, S. Schütz, J. Schachenmayer, C. Genet, G. Pupillo, and T. Ebbesen, Exploring Superconductivity under Strong Coupling with the Vacuum Electromagnetic Field (2019), arXiv:191101459 .
  * Thomas  _et al._ [2021] A. Thomas, E. Devaux, K. Nagarajan, G. Rogez, M. Seidel, F. Richard, C. Genet, M. Drillon, and T. Ebbesen, Large Enhancement of Ferromagnetism under a Collective Strong Coupling of YBCO Nanoparticles, Nano Lett. 21, 4365 (2021).
  * Rivera  _et al._ [2019] N. Rivera, T. Christensen, and P. Narang, Phonon Polaritonics in Two-Dimensional Materials, Nano Letters 19, 2653 (2019), 1901.05308 .
  * Berte  _et al._ [2018] R. Berte, C. R. Gubbin, V. D. Wheeler, A. J. Giles, V. Giannini, S. A. Maier, S. D. Liberato, and J. D. Caldwell, Sub-nanometer Thin Oxide Film Sensing with Localized Surface Phonon Polaritons, ACS Photonics 5, 2807 (2018).
  * Dai  _et al._ [2019] S. Dai, W. Fang, N. Rivera, Y. Stehle, B.-Y. Jiang, J. Shen, R. Tay, C. Ciccarino, Q. Ma, D. Rodan-Legrain, P. Jarillo-Herrero, E. Teo, M. Fogler, P. Narang, J. Kong, and D. Basov, Phonon Polaritons in Monolayers of Hexagonal Boron Nitride, Adv. Mater. 31, 1806603 (2019).
  * Yamanaka  _et al._ [2000] A. Yamanaka, M. Kataoka, Y. Inaba, K. Inoue, B. Hehlen, and E. Courtens, Evidence for competing orderings in strontium titanate from hyper-Raman scattering spectroscopy, Europhysics Letters 50, 688 (2000).
  * Rowley  _et al._ [2014] S. E. Rowley, L. J. Spalek, R. P. Smith, M. P. M. Dean, M. Itoh, J. F. Scott, G. G. Lonzarich, and S. S. Saxena, Ferroelectric quantum criticality, Nature Physics 10, 367 (2014), 0903.1445 .
  * Kozina  _et al._ [2019] M. Kozina, M. Fechner, P. Marsik, T. v. Driel, J. M. Glownia, C. Bernhard, M. Radovic, D. Zhu, S. Bonetti, U. Staub, and M. C. Hoffmann, Terahertz-driven phonon upconversion in SrTiO3, Nature Physics 15, 387 (2019), 1807.10788 .
  * Pálová  _et al._ [2009] L. Pálová, P. Chandra, and P. Coleman, Quantum critical paraelectrics and the Casimir effect in time, Physical Review B 79, 075101 (2009), 0803.1517 .
  * Esswein and Spaldin [2022] T. Esswein and N. Spaldin, Ferroelectric, quantum paraelectric, or paraelectric? Calculating the evolution from BaTiO3 to SrTiO3 to KTaO3 using a single-particle quantum mechanical description of the ions, Phys. Rev. Res. 4, 033020 (2022).
  * Zheng  _et al._ [2020] Z. Zheng, Q. Ma, Z. Bi, S. d. l. Barrera, M.-H. Liu, N. Mao, Y. Zhang, N. Kiper, K. Watanabe, T. Taniguchi, J. Kong, W. A. Tisdale, R. Ashoori, N. Gedik, L. Fu, S.-Y. Xu, and P. Jarillo-Herrero, Unconventional ferroelectricity in moiré heterostructures, Nature 588, 71 (2020).
  * Stern  _et al._ [2021] M. V. Stern, Y. Waschitz, W. Cao, I. Nevo, K. Watanabe, T. Taniguchi, E. Sela, M. Urbakh, O. Hod, and M. B. Shalom, Interfacial ferroelectricity by van der Waals sliding, Science 372, 1462 (2021), 2010.05182 .
  * Yasuda  _et al._ [2021] K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Stacking-engineered ferroelectricity in bilayer boron nitride, Science 372, 1458 (2021), 2010.06600 .
  * Woods  _et al._ [2021] C. R. Woods, P. Ares, H. Nevison-Andrews, M. J. Holwill, R. Fabregas, F. Guinea, A. K. Geim, K. S. Novoselov, N. R. Walet, and L. Fumagalli, Charge-polarized interfacial superlattices in marginally twisted hexagonal boron nitride, Nature Communications 12, 347 (2021).
  * Moore  _et al._ [2021] S. Moore, C. Ciccarino, D. Halbertal, L. McGilly, N. Finney, K. Yao, Y. Shao, G. Ni, A. Sternbach, E. Telford, B. Kim, S. Rossi, K. Watanabe, T. Taniguchi, A. Pasupathy, C. Dean, J. Hone, P. Schuck, P. Narang, and D. Basov, Nanoscale lattice dynamics in hexagonal boron nitride moiré superlattices, Nat. Comm. 12, 5741 (2021).
  * Katayama  _et al._ [2012] I. Katayama, H. Aoki, J. Takeda, H. Shimosato, M. Ashida, R. Kinjo, I. Kawayama, M. Tonouchi, M. Nagai, and K. Tanaka, Ferroelectric Soft Mode in a SrTiO3 Thin Film Impulsively Driven to the Anharmonic Regime Using Intense Picosecond Terahertz Pulses, Physical Review Letters 108, 097401 (2012).
  * Nova  _et al._ [2019] T. F. Nova, A. S. Disa, M. Fechner, and A. Cavalleri, Metastable ferroelectricity in optically strained SrTiO3, Science 364, 1075 (2019), 1812.10560 .
  * Li  _et al._ [2019] X. Li, T. Qiu, J. Zhang, E. Baldini, J. Lu, A. M. Rappe, and K. A. Nelson, Terahertz field–induced ferroelectricity in quantum paraelectric SrTiO3, Science 364, 1079 (2019), 1812.10785 .
  * Ashida  _et al._ [2021] Y. Ashida, A. İmamoğlu, and E. Demler, Cavity Quantum Electrodynamics at Arbitrary Light-Matter Coupling Strengths, Phys. Rev. Lett. 126, 153603 (2021).
  * Casimir and Polder [1948] H. Casimir and D. Polder, The Influence of Retardation on the London-van der Waals Forces, Phys. Rev. 73, 360 (1948).
  * Lifschitz [1956] E. Lifschitz, The Theory of Molecular Attractive Forces between Solids, Sov. Phys. J.E.T.P. 2, 73 (1956).
  * Dzyaloshinskii  _et al._ [1961] I. Dzyaloshinskii, E. Lifshitz, and L. Pitaevskii, GENERAL THEORY OF VAN DER WAALS’ FORCES, Sov. Phys. Usp. 4, 153 (1961).
  * Kenneth and Klich [2006] O. Kenneth and I. Klich, Opposites Attract: A Theorem about the Casimir Force, Phys. Rev. Lett. 160401, 97 (2006).
  * Philbin  _et al._ [2022] J. Philbin, T. Haugland, T. Ghosh, E. Ronca, M. Chen, P. Narang, and H. Koch, Molecular van der Waals fluids in cavity quantum electrodynamics (2022), arXiv:2202.07956 .
  * Note [1] This is similar to the case of Ginzburg-Landau theory where, in general the microscopic parameters such as density-of-states and quasiparticle gap are not directly related to the long-wavelength order parameter. Rather, only certain ratios of coupling constants are fixed in terms of microscopic parameters.
  * Note [2] The authors would like to especially acknowledge private correspondence with Yuto Ashida on this point, who patiently and carefully worked through this subtlety with us, and who ultimately helped uncover the origin of the discrepancies. Discussions with Ataç İmamoğlu were also valuable in this regard.
  * Zhou and Newns [1997] C. Zhou and D. Newns, Intrinsic dead layer effect and the performance of ferroelectric thin film capacitors, J. of App. Phys. 82, 3081 (1997).
  * Sirenko  _et al._ [2000] A. A. Sirenko, C. Bernhard, A. Golnik, A. M. Clark, J. Hao, W. Si, and X. X. Xi, Soft-mode hardening in SrTiO3 thin films, Nature 404, 373 (2000).
  * Katayama  _et al._ [2008] I. Katayama, H. Shimosato, D. S. Rana, I. Kawayama, M. Tonouchi, and M. Ashida, Hardening of the ferroelectric soft mode in SrTiO3 thin films, Applied Physics Letters 93, 132903 (2008).
  * Yang  _et al._ [2016] Q. Yang, J. Cao, Y. Zhou, L. Sun, and X. Lou, Dead layer effect and its elimination in ferroelectric thin film with oxide electrodes, Acta Mater. 112, 216 (2016).
  * Welakuh and Narang [2022] D. Welakuh and P. Narang, Nonlinear optical processes in centrosymmetric systems by strong-coupling-induced symmetry breaking (2022), arXiv:2202.11117 .
  * Khatoniar  _et al._ [2022] M. Khatoniar, R. Bushati, A. Mekawy, F. Dirnberger, A. Alù, and V. M. Menon, Relaxing Symmetry Rules for Nonlinear Optical Interactions in Van der Waals Materials via Strong Light–Matter Coupling, ACS Photonics 9, 503 (2022).
  * Kumah  _et al._ [2014] D. P. Kumah, A. Malashevich, A. S. Disa, D. A. Arena, F. J. Walker, S. Ismail-Beigi, and C. H. Ahn, Effect of Surface Termination on the Electronic Properties of LaNiO3 Films, Physical Review Applied 2, 054004 (2014).
  * Mannebach  _et al._ [2017] E. Mannebach, C. Nyby, F. Ernst, Y. Zhao, J. Tolsma, Y. Li, M.-J. Sher, I.-C. Tung, Z. H., Q. Zhang, K. Seyler, G. Clark, Y. Lin, D. Zhu, J. Glownia, M. Kozina, S. Song, S. Nelson, A. Mehta, Y. Yu, A. Pant, O. Aslan, A. Raja, Y. Guo, A. DiChiara, W. Mao, L. Cao, S. Tongay, J. Sun, D. Singh, T. Heinz, X. Xu, A. MacDonald, E. Reed, H. Wen, and A. Lindenberg, Dynamic Optical Tuning of Interlayer Interactions in the Transition Metal Dichalcogenides, Nano Lett. 17, 7761 (2017).
  * Basov  _et al._ [2016] D. N. Basov, M. M. Fogler, and F. J. G. d. Abajo, Polaritons in van der Waals materials, Science 354, aag1992 (2016).
  * Liu  _et al._ [2016] M. Liu, A. J. Sternbach, and D. N. Basov, Nanoscale electrodynamics of strongly correlated quantum materials, Reports on Progress in Physics 80, 014501 (2016).
  * Jing  _et al._ [2021] R. Jing, Y. Shao, Z. Fei, C. F. B. Lo, R. A. Vitalone, F. L. Ruta, J. Staunton, W. J.-C. Zheng, A. S. Mcleod, Z. Sun, B.-y. Jiang, X. Chen, M. M. Fogler, A. J. Millis, M. Liu, D. H. Cobden, X. Xu, and D. N. Basov, Terahertz response of monolayer and few-layer WTe2 at the nanoscale, Nature Communications 12, 5594 (2021).
  * Caldwell  _et al._ [2019] J. D. Caldwell, I. Aharonovich, G. Cassabois, J. H. Edgar, B. Gil, and D. N. Basov, Photonics with hexagonal boron nitride, Nature Reviews Materials 4, 552 (2019).
  * Xi  _et al._ [2015] X. Xi, Z. Wang, W. Zhao, J.-H. Park, K. Law, H. Berger, L. Forró, J. Shan, and K. Mak, Ising pairing in superconducting NbSe2 atomic layers, Nat. Phys. 12, 139 (2015).
  * Yu  _et al._ [2019] Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. Gu, X. Chen, and Y. Zhang, High-temperature superconductivity in monolayer Bi2Sr2CaCuO8+δ, Nature 156, 575 (2019).
  * Rizzo  _et al._ [2022] D. Rizzo, A. McLeod, C. Carnahan, E. Telford, A. Dismukes, R. Wiscons, Y. Dong, C. Nuckolls, C. Dean, A. Pasupathy, X. Roy, D. Xiao, and D. Basov, Visualizing Atomically Layered Magnetism in CrSBr, Adv. Mater. 34, 2201000 (2022).
  * Huang  _et al._ [2018] B. Huang, G. Clark, D. Klein, D. MacNeill, E. Navarro-Moratalla, K. Seyler, N. Wilson, M. McGuire, D. Cobden, D. Xiao, W. Yao, P. Jarillo-Herrero, and X. Xu, Electrical control of 2D magnetism in bilayer CrI3, Adv. Mater. 34, 2201000 (2018).
  * Lee  _et al._ [2021] K. Lee, A. Dismukes, E. Telford, R. Wiscons, J. Wang, X. Xu, C. Nuckolls, C. Dean, X. Roy, and X. Zhu, Magnetic Order and Symmetry in the 2D semiconductor CrSBr, Nano Lett. 21, 3511 (2021).
  * Klein  _et al._ [2022a] J. Klein, Z. Song, B. Pingault, F. Dirnberger, H. Chi, J. B. Curtis, R. Dana, R. Bushati, J. Quan, L. Dekanovsky, Z. Sofer, A. Alù, V. M. Menon, J. S. Moodera, M. Lončar, P. Narang, and F. M. Ross, Sensing the local magnetic environment through optically active defects in a layered magnetic semiconductor (2022a), arXiv:2207.02884 .
  * Klein  _et al._ [2022b] J. Klein, B. Pingault, M. Florian, M.-C. Heißenbüttel, A. Steinhoff, Z. Song, K. Torres, F. Dirnberger, J. B. Curtis, T. Deilmann, R. Dana, R. Bushati, J. Quan, J. Luxa, Z. Sofer, A. Alù, V. M. Menon, U. Wurstbauer, M. Rohlfing, P. Narang, M. Lončar, and F. M. Ross, The bulk van der Waals layered magnet CrSBr is a quasi-1D quantum material (2022b), arXiv:2205.13456 .
  * Bedoya-Pinto  _et al._ [2021] A. Bedoya-Pinto, J.-R. Ji, A. K. Pandeya, P. Gargiani, M. Valvidares, P. Sessi, J. M. Taylor, F. Radu, K. Chang, and S. S. P. Parkin, Intrinsic 2D-XY ferromagnetism in a van der Waals monolayer, Science 374, 616 (2021), 2006.07605 .
  * Klein  _et al._ [2021] J. Klein, T. Pham, J. D. Thomsen, J. B. Curtis, M. Lorke, M. Florian, A. Steinhoff, R. A. Wiscons, J. Luxa, Z. Sofer, F. Jahnke, P. Narang, and F. M. Ross, Atomistic spin textures on-demand in the van der Waals layered magnet CrSBr (2021), arXiv:2107.00037 .
  * Song  _et al._ [2022] Q. Song, C. A. Occhialini, E. Ergeçen, B. Ilyas, D. Amoroso, P. Barone, J. Kapeghian, K. Watanabe, T. Taniguchi, A. S. Botana, S. Picozzi, N. Gedik, and R. Comin, Evidence for a single-layer van der Waals multiferroic, Nature 602, 601 (2022).
  * Novoselov  _et al._ [2016] K. Novoselov, A. Mishchenko, A. Carvalho, and A. Castro Neto, 2D materials and van der Waals heterostructures, Science 353, aac9439 (2016).
  * Flick  _et al._ [2015] J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Kohn-Sham approach to quantum electrodynamical density-functional theory: Exact time-dependent effective potentials in real space, P.N.A.S. 112, 15285 (2015).
  * Flick  _et al._ [2019] J. Flick, D. Welakuh, M. Ruggenthaler, H. Appel, and A. Rubio, Light-Matter Response in Nonrelativistic Quantum Electrodynamics, A.C.S. Photonics 6, 2757 (2019).
  * Ruggenthaler  _et al._ [2014] M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. Tokatly, and A. Rubio, Quantum-electrodynamical density-functional theory: Bridging quantum optics and electronic-structure theory, Phys. Rev. A 90, 012508 (2014).
  * Abrikosov  _et al._ [1963] A. Abrikosov, L. Gor’kov, and I. Dzyaloshinskii,  _Methods of Quantum Field Theory in Statistical Physics_ (Dover Publications, 1963).

##  Appendix A Matsubara Formalism

In the Matsubara approach, we introduce the partition function via functional
integral

$Z=\int\mathcal{D}[Q,A,\chi]e^{-S[Q,A,\chi]}.$ (40)

The Matsubara action is written in the Weyl gauge with $\phi=0$. There is a
subtlety about writing the electromagnetic field in the Matsubara formalism,
which is that since the electric field $\mathbf{E}$ is in fact a canonical
momentum (it is the time derivative of the coordinate $\mathbf{A}$), it
carries an additional factor of $i$ when coupling to the polarization. Thus,
the action is written as the integral of the Lagrangian

$\mathcal{L}=\frac{1}{2}\left[\left(\frac{\partial\mathbf{A}}{\partial\tau}\right)^{2}+\left(\nabla\times\mathbf{A}\right)^{2}+\left(\frac{\partial\mathbf{Q}}{\partial\tau}\right)^{2}+\Omega_{0}^{2}\mathbf{Q}^{2}\right]\\\
+i\eta\mathbf{Q}\cdot\frac{\partial\mathbf{A}}{\partial\tau}+i\frac{1}{2}\chi(x)\mathbf{Q}^{2}+\frac{1}{4\lambda}\chi^{2}$
(41)

over space and imaginary time $\tau\in[0,\beta]$. Here we see the
hybridization between the phonon and photon through the electric dipole
coupling $\mathbf{Q}\cdot
i\partial_{\tau}\mathbf{A}\to\mathbf{Q}\cdot\mathbf{E}$ upon returning to real
time. Since the electromagnetic field experiences dispersion (due to the
magnetic induction), it must also be subjected to boundary conditions which in
this case are the same perfect-metal conditions we used previously.

The last term is a Hubbard-Stratonovich term used to decouple the quartic
phonon-phonon interaction in the Hartree channel. Here the factor of $i$
reflects the interaction is repulsive. Integration over $\chi$ can be
performed, and this returns the theory to its standard form in terms of
$\mathbf{A}$ and $\mathbf{Q}$ with a nonlinearity
$\sim\lambda(\mathbf{Q}^{2})^{2}$. Qualitatively, $\chi$ represents the
renormalization of the phonon frequency from bare value $\Omega_{0}$ to the
physical value $\Omega_{T}$, which is renormalized by the phonon-phonon
interactions.

We now show that up to one-loop order this also yields the same result as the
FDT calculation above. In particular, we expect the saddle-point in $\chi$ to
correspond to the self-consistent Hartree approximation. We can obtain this
formally exactly by integrating out the phonon and photon, which now appear
quadratically, to get

$\frac{\beta}{2\lambda}\chi(x)+\frac{\delta W[\chi]}{\delta\chi(x)}=0$ (42)

with functional determinant

$W[\chi]=\frac{1}{2}\mathbf{Tr}\log\mathbb{K}[A]$ (43)

where the kernel can be identified as the Gaussian part of the Matsubara
action.

This expression is simplified if we make an ansatz that the self-energy $\chi$
is time-independent in the saddle-point. Then, we can perform a shift to
“complete the square” for $\mathbf{Q}$, writing in the frequency domain

$\mathbf{Q}(x)=\delta\mathbf{Q}(x)-\frac{\eta\omega_{m}}{\omega_{m}^{2}+\Omega_{0}^{2}+i\chi(\mathbf{r})}\mathbf{A}(x).$
(44)

The first term characterizes the “unscreened” fluctuations, as we argued in
the first section, while the second term describes the screening due to the
electromagnetic field. By writing $\chi(\mathbf{r})$ we have allowed for the
possibility of an inhomogeneous shift in the phonon self-energy, as we expect
near the boundary of the system.

With this transformation the kernel decouples into the part due to
$\delta\mathbf{Q}(x)$ and the part due to the dielectric energy. We also see
that the variation with respect to $\chi$ can be framed as a variation with
respect to the dressed phonon frequency, since they are related by

$\Omega_{T}^{2}(\mathbf{r})=\Omega_{0}^{2}+i\chi(\mathbf{r})\Rightarrow\frac{\delta}{\delta\chi(\mathbf{r})}=i\frac{\delta}{\delta(\Omega_{T}^{2}(\mathbf{r}))}.$
(45)

We find, after performing the transformation that the functional has two
contributions;

$W_{0}[\chi]=\frac{1}{2}\sum_{\omega_{m}}\mathbf{Tr}\log\left[\mathds{1}\delta^{3}(\mathbf{r}^{\prime}-\mathbf{r})\left(\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r})\right)\right]$
(46)

from the unscreened response (in the absence of phonon dispersion this is
purely local and diverges with UV cutoff as $\Lambda^{d}$), and

$W_{scr}[\chi]=\frac{1}{2}\sum_{\omega_{m}}\mathbf{Tr}\log\left[\delta^{3}(\mathbf{r}^{\prime}-\mathbf{r})\left(\epsilon(\omega_{m},\mathbf{r})\omega_{m}^{2}\mathds{1}-\nabla^{2}+\nabla\nabla\cdot\right)\right]^{-1}$
(47)

from the dielectric screening. Now, all of the dependence on $\chi$ is
captured through

$\epsilon(\omega_{m},\mathbf{r})=1+\frac{\eta^{2}}{\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r})}.$
(48)

We therefore can easily evaluate the derivative in terms of the Matsubara
Green’s functions (note the minus sign is different than the usual definition
here)

$\displaystyle\mathscr{D}_{0}(\mathbf{r},\mathbf{r}^{\prime};\omega_{m})=\left[\mathds{1}\delta^{3}(\mathbf{r}^{\prime}-\mathbf{r})\left(\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r})\right)\right]^{-1}$
(49a)
$\displaystyle\mathscr{G}(\mathbf{r},\mathbf{r}^{\prime};\omega_{m})=\left[\delta^{3}(\mathbf{r}^{\prime}-\mathbf{r})\left(\epsilon(\omega_{m},\mathbf{r})\omega_{m}^{2}\mathds{1}-\nabla^{2}+\nabla\nabla\cdot\right)\right]^{-1}$
(49b)

as

$\frac{\delta
W[\chi]}{\delta\chi(\mathbf{r})}=i\frac{1}{2}\sum_{\omega_{m}}\mathbf{tr}\left[\mathscr{D}_{0}(\mathbf{r},\mathbf{r};\omega_{m})+\omega_{m}^{2}\frac{\delta\epsilon(\omega_{m})}{\delta\Omega_{T}^{2}}\mathscr{G}(\mathbf{r},\mathbf{r};\omega_{m})\right].$
(50)

We therefore find an equation for
$i\chi(\mathbf{r})=\Omega_{T}^{2}(\mathbf{r})-\Omega_{0}^{2}$ of

$\frac{1}{\lambda}\left(\Omega_{T}^{2}(\mathbf{r})-\Omega_{0}^{2}\right)\\\
=T\sum_{\omega_{m}}\mathbf{tr}\left[\mathscr{D}_{0}(\mathbf{r},\mathbf{r};\omega_{m})+\omega_{m}^{2}\frac{\delta\epsilon(\omega_{m},\mathbf{r})}{\delta\Omega_{T}^{2}}\mathscr{G}(\mathbf{r},\mathbf{r};\omega_{m})\right].$
(51)

$\Omega_{0}^{2}$ is a counter term which is set by the renormalization
condition that $\Omega_{T}^{2}(\mathbf{r})$ match the bulk value at a given
temperature.

We now recover the result from the previous section if we evaluate the right-
hand side to lowest order (i.e. not self-consistently) in $\lambda$, taking
$\Omega_{T}^{2}(\mathbf{r})=\Omega_{TO}^{2}$ to be the bulk TO mdoe frequency.
Then the right-hand side is nothing but
$\langle\mathbf{Q}(\mathbf{r},t)^{2}\rangle$ evaluated as a Matsubara sum.
This can now be evaluated efficiently, and in an unbiased manner, provided one
can evaluate the Green’s functions.

##  Appendix B Variational Approach

As a final sanity check, we also provide a derivation based on a variational
approach for the thermodynamic free energy. We again work in Matsubara
formalism, but instead of using a Hubbard-Stratonovich transformation we
employ the Feynman-Gibbs-Bogoliubov inequality. To this end, we write the
partition function as

$Z=\int\mathcal{D}[\mathbf{Q},\mathbf{A}]e^{-S}=\int\mathcal{D}[\mathbf{Q},\mathbf{A}]e^{-S_{\rm
eff}}e^{-(S-S_{\rm eff})}.$ (52)

The effective action is then chosen to constitute the variational ansatz. We
use

$S_{\rm eff}=\int
d^{3}r\sum_{\omega_{m}}\left[\frac{1}{2}\mathbf{Q}_{-\omega}(\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r}))\mathbf{Q}_{\omega}+\eta\omega_{m}\mathbf{Q}_{-\omega}\cdot\mathbf{A}_{\omega}\right]\\\
+S_{\rm Maxwell},$ (53)

which essentially replaces the bare TO frequency and interactions with a
local, effective TO frequency. Note this differs from the analysis of Ref.
[24] so far only in the choice of a local variational TO frequency as opposed
to a single global parameter. With respect to this ansatz, correlation
functions may be calculated easily using the above frameworks, since the
effective action used in the ansatz is quadratic and time independent.

We then obtain an upper bound on the free energy as a functional of our ansatz
via the Feynman-Bogoliubov-Gibbs inequality as

$W[\Omega_{T}^{2}(\mathbf{r})]\leq W_{\rm eff}+\langle S-S_{\rm
eff}\rangle_{\rm eff}.$ (54)

The first term is simply the free energy of the noninteracting ansatz, while
the second term becomes

$\langle S-S_{\rm eff}\rangle_{\rm eff}=\int
d^{4}x\left[\frac{\lambda}{4}\langle(\mathbf{Q}^{2}(x))^{2}\rangle+\frac{1}{2}(\Omega_{0}^{2}-\Omega_{T}^{2}(\mathbf{r}))\langle\mathbf{Q}^{2}(x)\rangle\right].$
(55)

This can be evaluated using Wick’s theorem. We now vary the parameter
$\Omega_{T}^{2}(\mathbf{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
$\mathbf{Q}=\delta\mathbf{Q}-\frac{\eta\omega_{m}}{\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r})}\mathbf{A}$.
This then allows us to express the free energy derivative as

$\frac{\delta W_{\rm
eff}}{\delta\Omega_{T}^{2}(\mathbf{r})}=\frac{1}{2}\sum_{\omega_{m}}\mathbf{tr}\mathscr{\hat{D}}_{0}(\mathbf{r},\mathbf{r};\omega_{m})+\frac{\delta\epsilon(\omega_{m},\mathbf{r})}{\delta\Omega_{T}^{2}(\mathbf{r})}\omega_{m}^{2}\mathbf{tr}\mathscr{\hat{G}}(\mathbf{r},\mathbf{r};\omega_{m})$
(56)

where the two Green’s functions are

$\mathscr{\hat{D}}_{0}(\mathbf{r},\mathbf{r}^{\prime};\omega_{m})=\langle\delta\mathbf{Q}_{-\omega}(\mathbf{r})\delta\mathbf{Q}_{\omega}(\mathbf{r}^{\prime})\rangle_{\rm
eff}$ (57)

and

$\mathscr{\hat{G}}(\mathbf{r},\mathbf{r}^{\prime};\omega_{m})=\langle\mathbf{A}_{-\omega}(\mathbf{r})\mathbf{A}_{\omega}(\mathbf{r}^{\prime})\rangle_{\rm
eff}$ (58)

and the dielectric constant is found as

$\epsilon(\omega_{m},\mathbf{r})=1+\frac{\eta^{2}}{\omega_{m}^{2}+\Omega_{T}^{2}(\mathbf{r})}.$
(59)

We note as well that by direct examination, we have

$\frac{\delta W_{\rm eff}}{\delta\Omega_{T}^{2}(\mathbf{r})}=\int
d\tau\frac{1}{2}\langle\mathbf{Q}^{2}(\mathbf{r})\rangle_{\rm
eff}=\frac{1}{2}\sum_{\omega_{m}}\langle\mathbf{Q}_{-\omega}(\mathbf{r})\cdot\mathbf{Q}_{\omega}(\mathbf{r})\rangle_{\rm
eff},$ (60)

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

$\langle
Q^{a}(x)Q^{a}Q^{b}(x)Q^{b}(x)\rangle=\left(\langle\mathbf{Q}(x)^{2}\rangle\right)^{2}+\frac{2}{d}\left(\langle\mathbf{Q}(x)^{2}\rangle\right)^{2},$
(61)

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
$\langle\mathbf{Q}(x)^{2}\rangle$ and the derivative of $W_{\rm eff}$ to get
the exact result

$W=W_{\rm eff}+\int
d^{3}r^{\prime}(\Omega_{0}^{2}-\Omega_{T}^{2}(\mathbf{r}))\frac{\delta W_{\rm
eff}}{\delta\Omega_{T}^{2}(\mathbf{r}^{\prime})}+\frac{u}{\beta}(1+\frac{2}{d})\left(\frac{\delta
W_{\rm eff}}{\delta\Omega_{T}^{2}(\mathbf{r}^{\prime})}\right)^{2}.$ (62)

Now, we take a derivative with respect to the parameter
$\Omega_{T}^{2}(\mathbf{r})$. After applying the chain rule, one finds

$\frac{\delta
W}{\delta\Omega_{T}^{2}(\mathbf{r})}=F-F+\left[\Omega_{0}^{2}-\Omega_{T}^{2}(\mathbf{r})\right]\frac{\delta
F}{\delta\Omega_{T}^{2}}+2\frac{u}{\beta}(1+\frac{2}{d})F\frac{\delta
F}{\delta\Omega_{T}^{2}},$

where $F=\frac{\delta W_{\rm eff}}{\delta\Omega_{T}^{2}(\mathbf{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

$\left[\Omega_{0}^{2}-\Omega_{T}^{2}(\mathbf{r})\right]+2uT(1+\frac{2}{d})\frac{\delta
W_{\rm eff}}{\delta\Omega_{T}^{2}(\mathbf{r})}=0.$ (63)

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 high-temperature. We begin with the quantum finite-temperature
Lagrangian describing the phonons and their coupling to the electromagnetic
field in a gauge independent form

$\mathscr{L}=\frac{1}{2}\left[\left(\frac{\partial\mathbf{Q}}{\partial\tau}\right)^{2}+\Omega_{0}^{2}\mathbf{Q}^{2}\right]+\frac{1}{2}\left[\left(\nabla\phi+\frac{\partial\mathbf{A}}{\partial\tau}\right)^{2}+(\nabla\times\mathbf{A})^{2}\right]+i\eta\mathbf{Q}\cdot\left(\frac{\partial\mathbf{A}}{\partial\tau}+\nabla\phi\right).$
(64)

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

$\mathscr{L}=\frac{1}{2}\left[\left(\frac{\partial\mathbf{Q}}{\partial\tau}\right)^{2}+\Omega_{0}^{2}\mathbf{Q}^{2}\right]+\frac{1}{2}\left[\left(\nabla\phi+\frac{\partial\mathbf{A}}{\partial\tau}\right)^{2}+(\nabla\times\mathbf{A})^{2}\right]-i\eta\left[\phi\nabla\cdot\mathbf{Q}+\mathbf{A}\cdot\frac{\partial\mathbf{Q}}{\partial\tau}\right].$
(65)

We choose the Coulomb gauge, where $\nabla\cdot\mathbf{A}=0$, such that
$\mathbf{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

$\mathscr{L}^{\parallel}_{\rm
int}=\frac{1}{2}\left(\nabla\phi\right)^{2}-i\eta\phi\nabla\cdot\mathbf{Q}^{\parallel}.$
(66)

This clearly endows the LO phonon mode with the long-range Coulomb
interaction. It also clearly couples to the phonon coordinate $\mathbf{Q}$ and
makes perfect sense in the classical limit, where $\mathbf{Q}$ and $\phi$
become “time-independent” and the integral over Matsubara time becomes simply
$\int_{0}^{\beta}d\tau\to 1/T$.

On the other hand, the transverse modes couple via

$\mathscr{L}^{\perp}_{\rm
int}=\frac{1}{2}\left(\frac{\partial\mathbf{A}}{\partial\tau}\right)^{2}+\frac{1}{2}(\nabla\times\mathbf{A})^{2}-i\eta\mathbf{A}\cdot\frac{\partial\mathbf{Q}^{\perp}}{\partial\tau}.$
(67)

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 $(\partial_{\tau}\mathbf{A})^{2}$
goes away as the cost of a thermal tunneling event $\sim T$ becomes
suppressed, and we are left with

$\mathscr{L}^{\perp}_{\rm
int}=\frac{1}{2}(\nabla\times\mathbf{A})^{2}-i\mathbf{A}\cdot\mathbf{J}^{\perp}.$
(68)

In this case the current $\mathbf{J}^{\perp}=\eta\partial_{\tau}\mathbf{Q}$ is
essentially the canonical momentum associated to the dielectric polarization,
but at high-temperatures this becomes uncorrelated with $\mathbf{Q}$.
Therefore, the photon field indeed dresses the expectation value for the
phonon current fluctuations $\langle(\mathbf{J}^{\perp})^{2}\rangle$, but this
now is independent of the phonon fluctuations $\langle\mathbf{Q}^{2}\rangle$.
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 $[\mathbf{Q},\mathbf{J}]\neq 0$, and at
high-temperatures this vanishes.

##  Appendix D Evaluation of Green’s Function in Cavity

Here we elaborate on the details for the calculations of the Green’s function
in the slab geometry. This largely follows Ref. [83], which calculates nearly
the same quantity we need but only evaluates at the boundary of the Fabry-
Perot geometry; we want to evaluate it locally in the bulk. We must obtain the
Matsubara Green’s function for the vector potential in the Weyl gauge. This
satisfies the equation

$\left[\omega_{m}^{2}\epsilon(\omega_{m})+\begin{pmatrix}-\partial_{z}^{2}&0&iq\partial_{z}\\\
0&q^{2}-\partial_{z}^{2}&0\\\
iq\partial_{z}&0&q^{2}\end{pmatrix}\right]\mathbb{\hat{G}}(z,z^{\prime})=\delta(z-z^{\prime})$
(69)

where $q$ is the in-plane momentum, which we have taken to lie along
$\mathbf{\hat{e}}_{x}$ (this can be done provided there is in-plane isotropy.

Clearly, the problem decoules in to the transverse electric (TE) modes, which
solve

$\left[\omega_{m}^{2}\epsilon(\omega_{m})+q^{2}-\partial_{z}^{2}\right]G_{yy}(z,z^{\prime})=\delta(z-z^{\prime})$
(70)

(they have the electric field polarized along $\mathbf{\hat{e}}_{y}$ which is
transverse to the momentum), and the transverse magnetic (TM) and longitudinal
(L) modes which are coupled and solve

$\left[\omega_{m}^{2}\epsilon(\omega_{m})+\begin{pmatrix}-\partial_{z}^{2}&iq\partial_{z}\\\
iq\partial_{z}&q^{2}\end{pmatrix}\right]\mathbb{\hat{G}}(z,z^{\prime})=\delta(z-z^{\prime}).$
(71)

We also have boundary conditions on the $x,y$ components at $\pm L/2$.

We begin with the TE modes. These can be solved for analytically by the method
of matching. We can write down the solution almost immediately (after some
deep introspection)

$G_{yy}(z,z^{\prime})=\begin{cases}A\sinh\kappa(L/2-z)\sinh\kappa(z^{\prime}+L/2)&z>z^{\prime}\\\
A\sinh\kappa(L/2-z^{\prime})\sinh\kappa(z+L/2)&z<z^{\prime}.\\\ \end{cases}$
(72)

We have introduced $\kappa=\sqrt{\omega_{m}^{2}\epsilon(\omega_{m})+q^{2}}$.
What now remains is to impose continuity of the derivative at $z=z^{\prime}$.
Integrating across the singularity and utilizing hyperbolic trig identities we
obtain the relation

$A\kappa\sinh\kappa L=1.$ (73)

We therefore obtain the TE mode Green’s function in analytical form of

$G_{yy}(z,z^{\prime})=\frac{1}{\kappa\sinh\kappa L}\\\
\times\begin{cases}\sinh\kappa(L/2-z)\sinh\kappa(z^{\prime}+L/2)&z>z^{\prime}\\\
\sinh\kappa(L/2-z^{\prime})\sinh\kappa(z+L/2)&z<z^{\prime}.\\\ \end{cases}$
(74)

The case of the remaining two modes is harder, in particular since the
equations are second order and coupled, allowing for the potential for a
fourth order equation upon decoupling. Let us write out all the equations
explicitly

$\displaystyle(\omega_{m}^{2}\epsilon(\omega_{m})-\partial_{z}^{2})G_{xx}(z,z^{\prime})+iq\partial_{z}G_{zx}(z,z^{\prime})=\delta(z-z^{\prime})$
(75a)
$\displaystyle\kappa^{2}G_{zz}(z,z^{\prime})+iq\partial_{z}G_{xz}(z,z^{\prime})=\delta(z-z^{\prime})$
(75b)
$\displaystyle(\omega_{m}^{2}\epsilon(\omega_{m})-\partial_{z}^{2})G_{xz}(z,z^{\prime})+iq\partial_{z}G_{zz}(z,z^{\prime})=0$
(75c) $\displaystyle
iq\partial_{z}G_{xx}(z,z^{\prime})+\kappa^{2}G_{zx}(z,z^{\prime})=0.$ (75d)

We can locally eliminate $G_{zx}$ to get

$G_{zx}(z,z^{\prime})=-\frac{1}{\kappa^{2}}iq\partial_{z}G_{xx}(z,z^{\prime})$
(76)

giving

$(\omega_{m}^{2}\epsilon(\omega_{m})-\partial_{z}^{2}+\frac{q^{2}}{\kappa^{2}}\partial_{z}^{2})G_{xx}(z,z^{\prime})=\delta(z-z^{\prime})$
(77)

for the in-plane component. This simplifies slightly to produce

$(\kappa^{2}-\partial_{z}^{2})G_{xx}(z,z^{\prime})=\frac{\kappa^{2}}{\omega_{m}^{2}\epsilon(\omega_{m})}\delta(z-z^{\prime})$
(78)

In fact, this can be solved just as trivially as the TE modes since we need
only replace the constant prefactor $A$. We obtain

$G_{xx}(z,z^{\prime})=\frac{\kappa}{\omega_{m}^{2}\epsilon(\omega_{m})\sinh\kappa
L}\\\
\times\begin{cases}\sinh\kappa(L/2-z)\sinh\kappa(z^{\prime}+L/2)&z>z^{\prime}\\\
\sinh\kappa(L/2-z^{\prime})\sinh\kappa(z+L/2)&z<z^{\prime}.\\\ \end{cases}$
(79)

All in all we then obtain for the contributions $G_{xx}(z,z)+G_{yy}(z,z)$
together

$\mathbf{tr}\mathbb{G}_{\parallel}(z,z)=\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}\left[\frac{1}{\kappa}+\frac{\kappa}{\omega_{m}^{2}\epsilon(\omega_{m})}\right].$
(80)

While this is not itself singular, we are reminded that we still must perform
the integral over $q$ and this will in general incur a cutoff dependence.

The last step is the most difficult; we must determine the normal component of
$G$. The normal component is found following Ref. [83], whereby we first
obtain the off-diagonal component

$G_{zx}(z,z^{\prime})=-\frac{iq}{\kappa^{2}}\partial_{z}G_{xx}(z,z^{\prime}).$
(81)

It then seems that an assumption has been made in Ref. [83], which seems we
must make as well, which is that the Green’s function is reciprocal in the
sense that we can interchange $G_{zx}$ and $G_{xz}$, so as to close the
equations. This seems valid provided time-reversal, or more importantly,
reciprocity of the system is preserved.

If we make this simplification, then we can obtain $G_{zz}$ from $G_{zx}$,
ultimately giving

$\kappa^{2}G_{zz}(z,z^{\prime})=(\omega_{m}^{2}\epsilon(\omega_{m})-\partial_{z}^{2})G_{xx}(z,z^{\prime}).$
(82)

We can in turn utilize the equation for $G_{xx}$ to remove the derivative,
which is very singular acting on the Green’s function. We then instead recover
an inhomogeneous equation of

$\kappa^{2}G_{zz}(z,z^{\prime})=-q^{2}G_{xx}(z,z^{\prime})+\frac{\kappa^{2}}{\omega_{m}^{2}\epsilon(\omega_{m})}\delta(z-z^{\prime}),$
(83)

such that we have

$G_{zz}(z,z^{\prime})=-\frac{q^{2}}{\kappa^{2}}G_{xx}(z,z^{\prime})+\frac{1}{\omega_{m}^{2}\epsilon(\omega_{m})}\delta(z-z^{\prime}).$
(84)

In fact, the last term is the origin of the most severe singularity, since the
Green’s function itself is singular at $z=z^{\prime}$, not merely non-
differentiable. However, we are saved by the fact that this term is
essentially constant and independent of the geometry. We therefore obtain the
complete expression for the local vector-potential fluctuations as

$\mathbf{tr}\mathbb{G}(z,z)=\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}\left[\frac{1}{\kappa}+\frac{\kappa}{\omega_{m}^{2}\epsilon(\omega_{m})}-\frac{q^{2}}{\kappa\omega_{m}^{2}\epsilon(\omega_{m})}\right]+\frac{1}{\omega_{m}^{2}\epsilon(\omega_{m})}\delta(0).$
(85)

The quantity in brackets can be simplified as

$\left[\frac{1}{\kappa}+\frac{\kappa}{\omega_{m}^{2}\epsilon(\omega_{m})}-\frac{q^{2}}{\kappa\omega_{m}^{2}\epsilon(\omega_{m})}\right]=\frac{2}{\kappa}.$
(86)

We therefore obtain the result for the electric-field fluctuations, which are
weighted by an additional $\omega_{m}^{2}$ factor. This gives

$\omega_{m}^{2}\mathbf{tr}\mathbb{G}(z,z)=\frac{\omega_{m}^{2}}{2\kappa}\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}+\frac{1}{\epsilon(\omega_{m})}\delta(0).$ (87)

We now must subtract the part which is independent of geometry or system size.
In particular, what matters for the Casimir formula is

$\Delta\langle
EE\rangle=\omega_{m}^{2}\left[\mathbf{tr}\mathbb{G}(z,z)-\mathbf{tr}\mathbb{G}(0,0)\bigg{|}_{L\to\infty}\right].$
(88)

We find a simple result of

$\Delta\langle
EE\rangle=\frac{\omega_{m}^{2}}{2\kappa}\left(\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}-\frac{1}{2}\right).$ (89)

Therefore, the result ends up being relatively simple (we still do need to
integrate over $q$ and sum over $\omega_{m}$, to be clear).

We can already see however that this is going to be positive comparing the
bulk and surface. In particular, we have

$\Delta\langle EE\rangle(L/2)-\Delta\langle
EE\rangle(0)=-\frac{\omega_{m}^{2}}{2\kappa}\sinh^{2}(\kappa L/2)/\sinh\kappa
L.$ (90)

This is therefore manifestly negative; now if we recall that the overall
contribution goes as $-\langle E^{2}\rangle$ 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

$\Delta\Omega_{T}^{2}(z)=\lambda
T\sum_{\omega_{m}}\int^{\Lambda}\frac{d^{2}q}{(2\pi)^{2}}\frac{\partial\epsilon(\omega_{m})}{\partial\Omega_{T}^{2}}\frac{\omega_{m}^{2}}{2\kappa}\left(\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}-\frac{1}{2}\right).$ (91)

This is to be numerically evaluated. We do need to be careful since the
integrand is singular in small $q$, $\omega$.

We therefore numerically compute

$\Delta\Omega_{T}^{2}(z)=\frac{\lambda}{2\pi}T\sum_{m=1}^{\frac{\omega_{c}}{2\pi
T}}\left[-\frac{\omega_{m}^{2}\eta^{2}}{(\omega_{m}^{2}+\Omega_{T}^{2})^{2}}\right]\int_{\sqrt{\omega_{m}^{2}\epsilon(\omega_{m})}}^{\sqrt{\omega_{m}^{2}\epsilon(\omega_{m})+\Lambda^{2}}}d\kappa\frac{1}{2\kappa}\left(\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}-\frac{1}{2}\right)\\\ +\frac{\lambda}{4\pi
T}\left[-\frac{\eta^{2}}{\Omega_{T}^{4}}\right]\lim_{\omega\to
0}\omega^{2}\int_{\omega\sqrt{\epsilon(0)}}^{\Lambda}d\kappa\frac{1}{2\kappa}\left(\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}-\frac{1}{2}\right).$ (92)

The first terms are the quantum corrections, which are evaluated by summing
over the $m=\pm 1,\pm 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 $\omega\to 0$. the
last term does indeed vanish in this limit. This is seen by first manipulating
the integrand via $x=\kappa L$ and utilizing trig identities to get

$\lim_{\omega\to
0}\omega^{2}\int_{\omega\sqrt{\epsilon(0)}}^{\Lambda}d\kappa\frac{1}{2\kappa}\left(\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}-\frac{1}{2}\right)=\lim_{\omega\to
0}\frac{\omega^{2}}{4}\int_{L\omega\sqrt{\epsilon(0)}}^{L\Lambda}\frac{dx}{x}\left[-1+\frac{\sinh
x(\frac{1}{2}-s)\sinh x(\frac{1}{2}+s)}{\sinh x}\right].$ (93)

We can now apply L’Hôspital’s rule to this to determine the result. Applying
this once we find

$\lim_{\omega\to
0}\frac{\omega^{2}}{4}\int_{L\omega\sqrt{\epsilon(0)}}^{L\Lambda}\frac{dx}{x}\left[-1+\frac{\sinh
x(\frac{1}{2}-s)\sinh x(\frac{1}{2}+s)}{\sinh x}\right]=\lim_{\omega\to
0}\frac{1}{2/\omega^{3}}\frac{1}{4}\frac{1}{x}\left[-1+\frac{\sinh
x(\frac{1}{2}-s)\sinh x(\frac{1}{2}+s)}{\sinh
x}\right]\bigg{|}_{x=L\sqrt{\epsilon(0)}\omega}.$ (94)

The first term in the brackets is non-singular so we find

$\lim_{\omega\to
0}\frac{\omega^{2}}{4}\int_{L\omega\sqrt{\epsilon(0)}}^{L\Lambda}\frac{dx}{x}\left[-1+\frac{\sinh
x(\frac{1}{2}-s)\sinh x(\frac{1}{2}+s)}{\sinh x}\right]=\lim_{\omega\to
0}\frac{\omega^{2}}{8}\frac{1}{L\sqrt{\epsilon(0)}}\left[\frac{\sinh
x(\frac{1}{2}-s)\sinh x(\frac{1}{2}+s)}{\sinh
x}\right]\bigg{|}_{x=L\sqrt{\epsilon(0)}\omega}.$ (95)

To proceed further we write

$\sinh x(\frac{1}{2}-s)\sinh x(\frac{1}{2}+s)=\frac{1}{2}(\cosh x-\cosh 2xs).$

Once again, the first term is not singular enough to possibly cancel the
$\omega^{2}$ numerator. so we find the only possible way out is from
evaluating

$\lim_{\omega\to
0}\frac{\omega^{2}}{4}\int_{L\omega\sqrt{\epsilon(0)}}^{L\Lambda}\frac{dx}{x}\left[-1+\frac{\sinh
x(\frac{1}{2}-s)\sinh x(\frac{1}{2}+s)}{\sinh x}\right]=\lim_{\omega\to
0}\frac{\omega^{2}}{16}\frac{1}{L\sqrt{\epsilon(0)}}\frac{\cosh
2\omega\sqrt{\epsilon(0)}z}{\sinh L\sqrt{\epsilon(0)}\omega}.$ (96)

Since $z$ is ultimately bounded by $L/2$ this also cannot scale with $z$ fast
enough to possibly yield a singular limit, so we see that applying L’Hôspitals
rule once more we will find this limit vanishes. Thus, we can disregard the
classical component altogether, since it does not contribute to this quantity.

This establishes that the electrostatic fluctuations do not end up
contributing, and we find that the result is simply due to the quantum
fluctuations vis a vis

$\Delta\Omega_{T}^{2}(z)=\frac{\lambda}{2\pi}T\sum_{m=1}^{\frac{\omega_{c}}{2\pi
T}}\left[-\frac{\omega_{m}^{2}\eta^{2}}{(\omega_{m}^{2}+\Omega_{T}^{2})^{2}}\right]\int_{\sqrt{\omega_{m}^{2}\epsilon(\omega_{m})}}^{\sqrt{\omega_{m}^{2}\epsilon(\omega_{m})+\Lambda^{2}}}d\kappa\frac{1}{2\kappa}\left(\frac{\sinh\kappa(L/2-z)\sinh\kappa(L/2+z)}{\sinh\kappa
L}-\frac{1}{2}\right).$ (97)