I

# Van der Waals Anomaly

Itai Y. Efrat and Ulf Leonhardt  Department of Physics of Complex Systems,
Weizmann Institute of Science, Rehovot 761001, Israel

###### Abstract

In inhomogeneous dielectric media the divergence of the electromagnetic stress
is related to the gradients of $\varepsilon$ and $\mu$, which is a consequence
of Maxwell’s equations. Investigating spherically symmetric media we show that
this seemingly universal relationship is violated for electromagnetic vacuum
forces such as the generalized van der Waals and Casimir forces. The stress
needs to acquire an additional anomalous pressure. The anomaly is a result of
renormalization, the need to subtract infinities in the stress for getting a
finite, physical force. The anomalous pressure appears in the stress in media
like dark energy appears in the energy-momentum tensor in general relativity.
We propose and analyse an experiment to probe the van der Waals anomaly with
ultracold atoms. The experiment may not only test an unusual phenomenon of
quantum forces, but also an analogue of dark energy, shedding light where
nothing is known empirically.

##  I Introduction

Van der Waals forces DLP ; Buhmann ; Shahmoon  dominate the microcosm of the
nanoworld: they cause electrically neutral atoms and molecules to interact
with each other on submicron scales. The force is generated by vacuum
fluctuations that electrically polarize these particles. Other vacuum forces
Forces  such as the Casimir force Casimir  we can understand Rodriguez  as the
net effect of van der Waals interactions between the constituents of
dielectric media (taking retardation into account CasimirPolder ). There they
cause stresses in the material described by the electric and the magnetic
components $\sigma_{E}$ and $\sigma_{M}$ of Abraham’s stress tensor Pita ;
Burger . Differences in these stresses create forces.

Van der Waals and Casimir forces in dielectric media need to be renomalized,
for otherwise their energy and stress were infinite. The renormalization is
local: it depends on the electric permittivity $\varepsilon$ and the magnetic
permeability $\mu$ as functions of space $\bm{r}$. We know this from
piece–wise homogeneos materials DLP  where the renomalizer depends on the
$\varepsilon$ and $\mu$ of each piece. For sandwiches of several materials,
the renormalization procedure DLP  does not only give finite results, but its
quantitative predictions DLP  have agreed with high–precision measurements
Levitation ; CasimirEquilibrium ; Decca . There is thus empirical evidence for
the locality of renormalization.

The question is: how local is renormalization? This question arises in
inhomogeneous media SimpsonSurprise  where $\varepsilon(\bm{r})$ and
$\mu(\bm{r})$ vary gradually. Suppose we would approximate a given
inhomogeneous medium with a sequence of homogeneous pieces, making each piece
finer and finer. The renormalized stress on each piece is always finite, but
the limit is not: the sequence of stresses diverges Simpson . From this
follows that, if macroscopic electromagnetism can account for vacuum forces in
media at all, the renormalizer must also depend on derivatives of
$\varepsilon(\bm{r})$ and $\mu(\bm{r})$. It should be still sufficiently
local, it cannot depend on all derivatives of $\varepsilon$ and $\mu$, for
otherwise the difference between the bare stress and the renormalizer would
vanish, the stress would get lost in renormalization. So how local is
renormaization?

We know from planar inhomogeneous media Grin1  (where $\varepsilon$ and $\mu$
vary in one direction) that the renormalizing Green function must depend on
the derivatives of $\varepsilon$ and $\mu$ up to second order, it is not
enough to take the gradients of $\varepsilon$ and $\mu$ into account.
Furthermore, second–order locality is not only necessary, but sufficient for
getting finite stresses, provided $\varepsilon$ and $\mu$ depend on frequency
and tend to unity sufficiently fast for large frequencies Grin1  — as is the
case for real materials LL8 . Experimental tests of the results of this
second–order renormalization procedure have been proposed Grin2 , but not
carried out yet.

In this paper, we take the next step and study spherically symmetric media.
Like in the planar case Grin1 , the symmetry preserves the polarizations of
the electromagnetic field that would normally get mixed in inhomogeneous
media, which considerably simplifies the problem. Spherically piece–wise
homogeneous media introduce problems of their own Milton80 ; Brevik82 ;
Milton97 ; Brevik98 ; Brevik99 ; Barton99 ; Bordag99 ; Avni ; Comment  but
here we avoid them by taking $\varepsilon(r)$ and $\mu(r)$ as gradually
varying with radius $r$.

In our paper we use Lifshitz theory DLP ; LL9 ; Lifshitz ; Scheel  to
calculate vacuum forces. This is the theory that agrees best with experiments
on the Casimir force Rodriguez ; Decca . Lifshitz theory uses the
fluctuation–dissipation theorem Scheel  to relate the quantum stress of the
vacuum to classical electromagnetic Green functions. The renormalization is
carried out by subtracting from the total Green function the outgoing part
such that only the scattered part remains. The physical picture behind this
renormalization procedure is the idea that van der Waals or Casimir forces DLP
are caused by the scattering of virtual electromagnetic waves at the
boundaries or inhomogeneities of media. The outgoing Green function depends on
the local dielectric environment, on $\varepsilon(\bm{r})$ and $\mu(\bm{r})$,
and so renormalization is local.

Like in the planar case Grin1 , we postulate second–order locality for the
renormalizing Green function. However, we find that we need to subtract not
only the zeroth–order Green function $D_{0}$, but also the first-order $D_{1}$
to get a finite vacuum stress. The Green function $D_{0}$ describes a purely
outgoing wave, whereas $D_{1}$ captures the first scattering at
inhomogeneities around the point of emission. As renormalization is local, the
choice of the renormalizing Green function must not depend on the global
symmetry of the material — whether it is planar or spherically symmetric must
not matter there. We thus conclude that one should always subtract both the
outgoing wave and first–order scattering in renormalization.

It is not entirely clear what the meaning of this modified renormalization is.
Subtracting the outgoing wave has a definite physical meaning: what is left
are the scattered virtual electromagnetic waves, and those are the ones
causing vacuum forces. But what is the meaning of the first–order scattering
we need to subtract as well? Maybe the best description for renormalization is
this: in the context of van der Waals and Casimir forces, renormalization is
the removal of the near field. Each infinitesimal cell of the medium acts like
an emitter and a receiver of virtual electromagnetic waves, but the receiver
should not interact with the near field of the emitter. In our procedure we
identify and subtract this field.

In adopting this renormalization procedure, we deduce a remarkable feature of
generalized van der Waals forces in inhomogeneous media: we find that the
electric and magnetic stresses $\sigma_{E}$ and $\sigma_{M}$ should be
replaced as

$\sigma_{F}\rightarrow\sigma_{F}+\frac{p}{2}\mathbb{1}$ (1)

with $p$ given by Eqs. (30) and (35) for spherically symmetric media. This
phenomenon we call van der Waals anomaly, because it corresponds to the trace
anomaly of conformally invariant fields in curved space–time Wald  where the
energy–momentum tensor $T^{\alpha}_{\beta}$ gets replaced by

$T^{\alpha}_{\beta}\rightarrow
T^{\alpha}_{\beta}+\varepsilon_{\Lambda}\delta^{\alpha}_{\beta}\,.$ (2)

The energy–momentum contribution
$\varepsilon_{\Lambda}\delta^{\alpha}_{\beta}$ has exactly the form of the
cosmological constant Weinberg ; HawkingInflation  written on the right–hand
side of Einstein’s field equations LL2 , but it is not necessarily constant.
Such a contribution to $T^{\alpha}_{\beta}$ is called dark energy Dark .
Following other analogues of cosmology in condensed–matter physics Volovik ;
Visser ; FF ; Kolo1 ; Kolo2 ; Kolo3 ; Stein21  the van der Waals anomaly thus
establishes an analogue of dark energy.

A problem of our modified renormalization procedure is the fact that first
scattering, described by $D_{1}$ in Eq. (29), is only defined in the
asymptotic expansion of the Green function for large frequencies. This
expansion captures the infinities of the stress, because those appear in the
high–frequency limit, but it does not necessarily describe correctly the
finite part of the stress. However, the central quantity of this paper, the
anomalous pressure $p$, we calculate to leading order in frequency.
Corrections to $p$ will be practically negligible.

What exactly is the quantity we calculate? As renormalization is local, it may
affect identities for the electromagnetic field one would normally take for
granted, as they follow from Maxwell’s equations, such as

$\nabla\cdot(\sigma_{E}+\sigma_{M})=\frac{\nabla\varepsilon}{\varepsilon}\,\mathrm{tr}\,\sigma_{E}+\frac{\nabla\mu}{\mu}\,\mathrm{tr}\,\sigma_{M}$
(3)

for the electric and magnetic components of Abrahams’s stress tensor with

$\sigma_{E}=\bm{E}\otimes\bm{D}-\frac{\bm{E}\cdot\bm{D}}{2}\mathbb{1}\,,\quad\sigma_{M}=\bm{B}\otimes\bm{H}-\frac{\bm{B}\cdot\bm{H}}{2}\mathbb{1}$
(4)

where in isotropic dielectric media and in SI units:

$\bm{D}=\varepsilon_{0}\varepsilon\bm{E}\,,\quad\bm{B}=\mu_{0}\mu\bm{H}\,,\quad\varepsilon_{0}\mu_{0}=c^{-2}\,.$
(5)

The Abraham identity of Eq. (3) relates the divergence of the stress to
gradients in $\varepsilon$ and $\mu$. In our paper we show that the
electromagnetic stresses need to be replaced according to Eq. (1) for the
identity to hold. This is the phenomenon we call van der Waals anomaly. The
additional pressure $p$ in the stress we call the anomalous pressure.

The anomaly is a direct consequence of the locality of our near–field
renormalization that takes first scattering into account. In particular, it is
the result of the violation of reciprocity by renormalization Wald , as
follows. Since the renormalizing Green function depends on the local
dielectric environment, the emitted and the received Green functions will be
different, violating reciprocity. The exchange of virtual electromagnetic
waves creates the van der Waals and Casimir forces DLP ; Buhmann , and so Case
the recoil imbalance causes a momentum imbalance: the pressure $p$. This
anomalous pressure depends only on local dielectric properties, and we have
found explicit expressions, Eqs. (30) and (35), for the planar and the
spherical case.

Our analysis indicates that the dilute limit of vacuum forces (the limit of
vanishing density) is dominated by local dielectric quantities. This is
because the van der Waals anomaly dominates the dilute limit. We see this as
follows: Eq. (3) with replacement (1) implies that
$\nabla\cdot(\sigma_{E}+\sigma_{M})+\nabla p$ vanishes quadratically in
density, because $\nabla\varepsilon$ and $\nabla\mu$ go at least linearly with
density and the renormalized stresses are also at least linear in density.
Now, $p$ depends on the local values of $\varepsilon$ and $\mu$ and their
derivatives up to second order, which implies that $p$ is linear in density.
We thus conclude that $\nabla\cdot(\sigma_{E}+\sigma_{M})\sim-\nabla p$ in the
dilute limit and that the vacuum forces are local. From this also follows that
$\nabla\cdot(\sigma_{E}+\sigma_{M})$ is linear in density in the dilute limit
Avni , in contrast to what has been widely believed Comment  without
questioning.

The dilute limit opens up the possibility of experimental tests of the van der
Waals anomaly with ultracold atoms Pethick . In our paper we outline the basic
ideas for such tests and estimate the required measurement precision. Note
that such experiments would not only probe the anomaly, but also additional
assumptions we need to make: we assume that the anomaly gives the electric and
the magnetic stress directly and equally, and we assume that the dielectric
force is the Helmholtz force LL8  even for vacuum fluctuations. Our estimation
indicates that these measurements are within the range of current experiments.

Such experiments may test more than just a curious phenomenon in Casimir
physics, as the van der Waals anomaly is an analogue of dark energy in
condensed–matter physics. While dark energy amounts to the most abundant form
of mass ($70\%$) in the universe CMBPlanck , its theoretical understanding has
been riddled with problems Weinberg . Testing an analogue of dark energy in
the laboratory may give some much–needed empirical guidance to elucidate the
matter.

How exactly is the van der Waals anomaly in dielectric media related to the
trace anomaly Wald  in space–time? Dielectric media appear to electromagnetic
fields like curved spaces Gordon ; Plebanski ; Schleich ; LeoPhil . In the
case of impedance–matched media,

$\varepsilon=\mu=n\,,$ (6)

the effective geometry is exact, with line element

$ds^{2}=n^{2}d\bm{r}^{2}$ (7)

and, consequently LL2 , spatial metric tensor $g_{ij}=n^{2}\mathbb{1}$, its
inverse $g^{ij}=n^{-2}\mathbb{1}$ and determinant $g=n^{6}$. For impedance
matching, the electric and magnetic vacuum stresses are the same,
$\sigma_{E}=\sigma_{M}$, and the total stress is a tensor density LL2 :
$\sigma=\sigma_{E}+\sigma_{M}=\sqrt{g}\,\tau$ with $\tau$ being a symmetric
tensor $\tau^{i}_{j}$. We obtain for the covariant divergence of a symmetric
tensor LL2 :

$\displaystyle\tau^{i}_{j;i}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{g}}\,\partial_{i}\sqrt{g}\,\tau^{i}_{j}-\frac{1}{2}\left(\partial_{j}g_{ik}\right)\tau^{ik}$
(8) $\displaystyle=$
$\displaystyle\frac{1}{n^{3}}\left(\nabla\cdot\sigma-\frac{\nabla
n}{n}\,\mathrm{tr}\,\sigma\right)\,,$ (9)

which gives Eq. (3) in the impedance–matched case of Eq. (6) if we require
$\tau^{i}_{j;i}=0$. Therefore, in impedance–matched media where the medium
corresponds to an exact geometry, the Abraham identity is equivalent to the
covariant conservation of momentum. In curved space–time, the trace anomaly
(2) plays the same role: it serves to restore the covariant conservation of
energy–momentum Wald . The van der Waals anomaly in dielectric media is thus
an analogue of the trace anomaly in space–time, and the trace anomaly appears
as dark energy in Einstein’s field equations LL2 .

The van der Waals anomaly might be more than just an analogue, it might
actually be the archetype of a mechanism Annals  that explains dark energy
without the need of new physics, as an effect of quantum electromagnetism.
This explanation appears to be in agreement with astronomical data Dror
resolving the tension between the measured Hubble constant and its prediction
from the Cosmic Microwave Background Verde . It is therefore possible that the
same force that rules the microcosm of the nanoworld does also rule the
macrocosmos of the expanding universe. Time and tests will tell.

##  II Spherical symmetry

Our theory is based on the Lifshitz theory LL9 ; Lifshitz ; Scheel  of vacuum
forces in dielectric media Forces  where renormalization is performed by
subtracting the outgoing wave from the bare Green function. The
electromagnetic Green function $\mathrm{G}$ describes the vector potential at
position $\bm{r}$ created by a point dipole at $\bm{r}_{0}$ oscillating with
frequency $\omega$. It is a bi–tensor with one index referring to the
direction of the vector potential and the other to the dipole direction. The
vacuum stress is calculated from $\mathrm{G}$ in the limit
$\bm{r}_{0}\sim\bm{r}$ by integration along positive imaginary frequencies and
differentiations with respect to $\bm{r}$ and $\bm{r}_{0}$. If there is a
symmetry that preserves the two polarizations of light, the Green bi–tensor
gets reduced to two scalar Green functions. In the planar case, these are the
Green functions $g_{\mathrm{E}}$ and $g_{\mathrm{M}}$ of the electric and the
magnetic polarizations. There either the electric or the magnetic field
remains orthogonal to the propagation direction and the direction in which the
medium varies. In a general inhomogeneous medium, the polarizations mix; it
takes a high degree of symmetry to preserve them.

###  II.1 Setting the scene

Consider a rotationally symmetric medium. We employ spherical coordinates
$\\{r,\theta,\phi\\}$ and require

$\varepsilon=\varepsilon(r)\,,\quad\mu=\mu(r)\,.$ (10)

Without writing this explicitly, we also assume that $\varepsilon$ and $\mu$
depend on frequency (which will become important for renormalization Grin1 ).
In spherically–symmetric media, the stress depends only on $r$ and has no
shear,

$\sigma=\mathrm{diag}(\sigma_{r}^{r},\sigma_{\theta}^{\theta},\sigma_{\phi}^{\phi})$
(11)

with

$\sigma_{\theta}^{\theta}=\sigma_{\phi}^{\phi}\,.$ (12)

The simplest of such media is a homogeneous dielectric sphere embedded in a
uniform background Shell  (with $\varepsilon=\varepsilon_{1}$, $\mu=\mu_{1}$
inside and $\varepsilon=\varepsilon_{2}$, $\mu=\mu_{2}$ outside the radius of
the sphere). Despite extensive effort Boyer ; MDS ; Balian ; Nesterenko ;
Milton80 ; Brevik82 ; Milton97 ; Brevik98 ; Brevik99 ; Barton99 ; Bordag99 ;
Avni ; Comment  there is still no general solution for the Casimir force in
this case. One of the major difficulties of this seemingly simple problem is
posed by the sharp interface between the sphere and the background. The issue
is that even renormalized vacuum fields tend to infinity at interfaces LeoBook
. For piece–wise homogeneous planar media these infinities do not matter, but
for spherically symmetric media they do. To see this, consider the divergence
of the stress (that gives the force density in mechanical equilibrium Pita ).
We get in spherical coordinates Divergence

$(\nabla\cdot\sigma)_{r}=\partial_{r}\sigma_{r}^{r}+\frac{2}{r}\,\sigma_{r}^{r}-\frac{2}{r}\,\sigma_{\theta}^{\theta}\,.$
(13)

At a spherical interface, $(\nabla\cdot\sigma)_{r}$ must be a delta function,
and so the radial stress $\sigma_{r}^{r}$ must be discontinuous. But both
$\sigma_{r}^{r}$ and $\sigma_{\theta}^{\theta}$ vary and tend to infinity at
the interface Avni . These diverging contribution cancel each other in the
force density Avni , but their cancellation is subtle Avni  and not fully
understood when $\varepsilon$ and $\mu$ depend on frequency. For planar media,
on the other hand, the stress $\sigma_{z}^{z}$ orthogonal to the interface is
finite and piece–wise constant, and the diverging stresses in the other
directions, $\sigma_{x}^{x}=\sigma_{y}^{y}$, do not contribute to the force,
as $(\nabla\cdot\sigma)_{z}=\partial_{z}\sigma_{z}^{z}$. There this problem
does not occur. Here we avoid it by assuming $\varepsilon(r)$ and $\mu(r)$ to
be continuous.

###  II.2 Spectral stresses

Lifshitz theory relates the stress to the electromagnetic Green function in
the form of an integral along positive imaginary frequencies Scheel :

$\sigma=\left.-\frac{\hbar
c}{2\pi}\int_{0}^{\infty}W\,d\kappa\right|_{\bm{r}_{0}\rightarrow\bm{r}}$ (14)

written here in terms of the imaginary wavenumber $\kappa$ for the frequency
$\omega=ic\kappa$. The spectral stresses $W$ depend on the Green function, and
on the source and observation points $\bm{r}_{0}$ and $\bm{r}$. The limit
$\bm{r}_{0}\rightarrow\bm{r}$ is taken in the integral of Eq. (14) but not
before. Spherical symmetry implies that the $W$ are functions of $r$ and
$r_{0}$, and of the angle $\gamma$ between the origin, $\bm{r}$ and
$\bm{r}_{0}$. The limit $\bm{r}_{0}\rightarrow\bm{r}$ corresponds then to
$r_{0}\rightarrow r$ and $\gamma\rightarrow 0$.

We use the symbol $W_{jF}^{jP}$ to denote the contribution of the
$P$-polarization to the $F$–component of the stress in $j$-direction, and
obtain (Appendix A):

$\displaystyle W_{rE}^{rE}$ $\displaystyle=$
$\displaystyle\varepsilon\kappa^{2}g_{E}\,,$ $\displaystyle W_{rM}^{rE}$
$\displaystyle=$
$\displaystyle-\frac{\partial_{r}r\,\partial_{r_{0}}r_{0}\,g_{E}}{\mu
rr_{0}}-W_{\theta M}^{\theta E}\,,$ $\displaystyle W_{\theta M}^{\theta E}$
$\displaystyle=$
$\displaystyle\frac{\partial_{\gamma}\sin\gamma\,\partial_{\gamma}g_{E}}{\mu
r^{2}\sin\gamma}$ (15)

for the electric polarization, and the equivalent expressions for the magnetic
polarization with $\varepsilon\leftrightarrow\mu$ and $E\leftrightarrow M$.
All other spectral stresses vanish. In the limit of large radii the
expressions (15) reproduce the planar case Grin1  where
$W^{zE}_{z}=\varepsilon\kappa^{2}g_{E}-\mu^{-1}\nabla\cdot\nabla_{0}\,g_{E}$
for the electric polarization and
$W^{zM}_{z}=\mu\kappa^{2}g_{M}-\varepsilon^{-1}\nabla\cdot\nabla_{0}\,g_{M}$
for the magnetic polarization.

###  II.3 Scalar Green functions

The two polarizations are described by the two scalar Green functions $g_{P}$.
Like in the planar case Grin1  we write them as

$g=\nu\nu_{0}D$ (16)

(dropping the polarization–index) with $\nu_{0}=\nu(r_{0})$ and

$\nu_{E}=\mu\,,\quad\nu_{M}=\varepsilon\,.$ (17)

We show in Appendix A that $D$ obeys the wave equation

$\frac{\nabla\cdot\nu\nabla
D}{n^{2}\nu}-\frac{R_{r}^{r}}{2}D-\kappa^{2}D=\frac{\delta(\bm{r}-\bm{r}_{0})}{n^{2}\nu}$
(18)

in terms of the refractive index $n$ with

$n^{2}=\varepsilon\mu$ (19)

and the abbreviation

$R_{r}^{r}=\frac{2}{n^{2}\nu}\left(\frac{\nu^{\prime
2}}{\nu}-\nu^{\prime\prime}-\frac{\nu^{\prime}}{\nu}\right)$ (20)

where the primes denote differentiations with respect to $r$. In the
impedance–matched case (6) when the medium establishes an exact geometry with
line element (7) the term $n^{-3}\nabla\cdot n\nabla$ in the wave equation
(18) is the Laplacian of a scalar and $R_{r}^{r}$ becomes the radial component
of the Ricci tensor LL2 . In the planar case Grin1 , one obtains the same wave
equation with $R_{r}^{r}$ replaced by $R_{z}^{z}$. The wave equations for $D$
thus assume an entirely geometrical form. Finally, for spherical symmetry we
have:

$\frac{\nabla\cdot\nu\nabla}{\nu}=\frac{\partial_{r}r^{2}\nu\,\partial_{r}}{r^{2}\nu}+\frac{\partial_{\gamma}\sin\gamma\,\partial_{\gamma}}{r^{2}\sin\gamma}\,.$
(21)

Expressions (16-21) determine the scalar Green functions and hence the
spectral stresses.

###  II.4 Abraham identity

The Abraham identity, Eq. (3), follows from Maxwell’s equations in general.
Let us see how it follows from the wave equation (18) in our case. We obtain
from the expression for $W_{rE}^{rE}$ in Eq. (15) by straightforward
differentiation

$\left(\partial_{r}+\partial_{r_{0}}+\frac{2}{r}-\frac{\varepsilon^{\prime}}{\varepsilon}\right)W_{rE}^{rE}=\mathbb{D}_{0}D_{E}$
(22)

in the limit $\bm{r}_{0}\sim\bm{r}$ where $\mathbb{D}_{0}$ denotes the
differential operator

$\mathbb{D}_{0}=\frac{\kappa^{2}}{r_{0}^{2}}\,\partial_{r_{0}}r_{0}\mu_{0}\,n^{2}+\frac{\kappa^{2}}{r^{2}}\,\partial_{r}r\mu\,n_{0}^{2}\,.$
(23)

We also obtain for $W_{rM}^{rE}$ and $W_{\theta M}^{\theta E}$ in the limit
$\bm{r}_{0}\sim\bm{r}$:

$\left(\partial_{r}+\partial_{r_{0}}+\frac{2}{r}-\frac{\mu^{\prime}}{\mu}\right)W_{rM}^{rE}-\left(\frac{2}{r}+\frac{2\mu^{\prime}}{\mu}\right)W_{\theta
M}^{\theta E}=\mathbb{D}_{1}D_{E}$ (24)

with the differential operator

$\displaystyle\mathbb{D}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{r_{0}^{2}}\,\partial_{r_{0}}r_{0}\mu_{0}\,n^{2}\left(-\frac{\nabla\cdot\mu\nabla}{n^{2}\mu}+\frac{R_{r}^{r}}{2}\right)$
(25)
$\displaystyle+\frac{1}{r^{2}}\,\partial_{r}r\mu\,n_{0}^{2}\left(-\frac{\nabla_{0}\cdot\mu_{0}\nabla_{0}}{n_{0}^{2}\mu_{0}}+\frac{R_{0r}^{r}}{2}\right)$

using Eq. (21) for the Laplacians. For the magnetic polarization we obtain the
corresponding expressions by interchanging $E\leftrightarrow M$ and
$\varepsilon\leftrightarrow\mu$.

Inspecting the differential operators, we see that the sum
$\mathbb{D}_{0}+\mathbb{D}_{1}$ applied to $D_{E}$ contains the left–hand side
of the wave equation (18) for $\bm{r}$ and also for $\bm{r}_{0}$. For
$\bm{r}\neq\bm{r}_{0}$ the right–hand side of Eq. (18) vanishes. So the sum of
the left–hand sides of Eqs. (22) and (24) vanishes, too, provided $g_{E}$
satisfies the same wave equation in $\bm{r}_{0}$ as it does in $\bm{r}$. This
is the case if the Green functions are reciprocal,

$D(\bm{r},\bm{r}_{0})=D(\bm{r}_{0},\bm{r})\,.$ (26)

The stresses are calculated after the limit $\bm{r}_{0}\rightarrow\bm{r}$ is
taken in Eq. (14), which implies that $\partial_{r}\sigma$ corresponds to the
sum of $\partial_{r}W$ and $\partial_{r_{0}}W$ in the spectral stresses. These
and the $2/r$ terms in Eqs. (22) and (24) form the divergence of the stress
according to Eq. (13), and
$\sigma_{r}^{r}+2\sigma_{\theta}^{\theta}=\sigma_{r}^{r}+\sigma_{\theta}^{\theta}+\sigma_{\phi}^{\phi}$
gives $\mathrm{tr}\,\sigma$. We thus obtain the Abraham identity (3), provided
the Green functions are reciprocal. Yet renormalization breaks reciprocity, as
we show next.

##  III Renormalization

The field of Casimir physics began with Casimir’s three–page paper Casimir  in
1948. There he did two things: calculate the force between two perfect mirrors
and handwave the justification for renormalization. The latter has haunted the
field ever since. Renormalization is necessary, because the bare stress
calculated, for example, using the procedure of Sec. II is infinite. One may
view the infinity of the bare vacuum energy and stress as an inevitable
feature of quantum field theory or as an artefact of the theory. In either
case, the infinity needs to be subtracted to get meaningful results. This
process of subtraction is called renormalization. In Lifshitz renormalization
DLP  the subtraction is performed on the Green functions. There the diverging
part stems from high frequencies and is therefore well–captured in geometrical
optics.

###  III.1 Geometrical optics

Let us briefly summarize geometrical optics in spherically–symmetric media
(further details in Appendix B). In the asymptotic limit of large $\kappa$ the
Green function $D$ approaches

$D\sim D_{0}+D_{1}$ (27)

with

$D_{0}={\cal A}_{0}\,e^{-\kappa s}$ (28)

and (for $\bm{r}\sim\bm{r}_{0}$)

$D_{1}=\beta_{1}d_{1}\,,\quad d_{1}=\frac{{\cal A}_{0}s}{\kappa}\,e^{-\kappa
s}$ (29)

where $s$ denotes the optical path length with increment $ds$ given by Eq. (7)
and ${\cal A}_{0}$ the amplitude. The contribution $D_{0}$ describes a purely
outgoing wave, while $D_{1}$ captures the first reflection in inhomogeneities
near the point of emission Grin1 .

Renormalization depends on the local values of $\varepsilon$ and $\mu$.
Consider the vicinity of $\bm{r}_{0}$ up to quadratic order. One verifies that
$D_{0}+D_{1}$ satisfies the wave equation (18) up to $O(\kappa^{-1})$. One
also finds that both $s$ and ${\cal A}_{0}$ are reciprocal up to quadratic
order (Appendix B). From this follows, according to Eqs. (28) and (29), that
$D_{0}$ and $d_{1}$ are reciprocal. However, the prefactor $\beta_{1}$ of
$d_{1}$ turns out to be nonreciprocal: it is given by (Appendix B):

$\beta_{1}=\frac{1}{24n^{2}}\left(\frac{n^{\prime
2}-2n\nabla^{2}n}{n^{2}}+\frac{6\nu\nu^{\prime\prime}-9\nu^{\prime
2}}{\nu^{2}}\right)$ (30)

evaluated at $\bm{r}=\bm{r}_{0}$. Equation (30) closely resembles the result
Grin1  for planar media where $\nabla^{2}n=n^{\prime\prime}$ and the primes
denote $\partial_{z}$ whereas here
$\nabla^{2}n=n^{\prime\prime}+(2/r)n^{\prime}$ and the primes stand for
$\partial_{r}$. Unlike the planar case Grin1 , the spherical $D_{1}$ does
contribute to the stress, and this contribution diverges (Appendix C). For
getting a finite result one therefore needs to subtract from the bare Green
function $D$ both $D_{0}$ and $D_{1}$. Appendix D proves the convergence of
the stress calculated according to this procedure. As renormalization is
local, it must not depend on the global symmetry of the medium — whether it is
planar or spherically symmetric should not play any role. We thus conclude
that we should renormalize with $D_{0}+D_{1}$ in general, even in the planar
case.

###  III.2 Reciprocity violation

The reciprocity–violating first–scattering amplitude $\beta_{1}$ is going to
be responsible for the violation of the Abraham identity (3). Consider the
spectral stresses $W$ in Eq. (14). Here the violation of the Abraham identity
is described by the sum of the left–hand sides of Eqs. (22) and (24) with
definitions (23) and (25). We call this sum $Q$. The wave equation (18) would
remove all of the terms in $Q$ if not only $D_{0}$ were reciprocal, but also
$D_{1}=\beta_{1}d_{1}$. However, in
$\nabla\cdot\nu\nabla\,\beta_{1}(r_{0})d_{1}$ the amplitude $\beta_{1}$
appears as a constant, whereas in
$\nabla_{0}\cdot\nu_{0}\nabla_{0}\,\beta_{1}(r_{0})d_{1}$ we get the extra
terms
$d_{1}\,\nabla_{0}\cdot\nu_{0}\nabla_{0}\,\beta_{1}+2\nu_{0}(\nabla_{0}\beta_{1})\cdot(\nabla_{0}d_{1})$.
We thus obtain:

$Q=-\frac{1}{r}\,\partial_{r}\,r\nu\left(\frac{(\partial_{r_{0}}\,r_{0}^{2}\nu_{0}\,\partial_{r_{0}}\beta_{1})}{r_{0}^{2}\nu_{0}}+2(\partial_{r_{0}}\beta_{1})\partial_{r_{0}}\right)d_{1}$
(31)

evaluated for $\bm{r}_{0}\rightarrow\bm{r}$. Note that Eq. (31) is only valid
up to $O(\kappa^{-1})$ since $D_{0}+D_{1}$ solves the wave equation only up to
this order. For large $\kappa$ the dominant term in $Q$ is given by the
highest power in $\kappa$. According to Eq. (29) this term is produced by
$\partial_{r}\partial_{r_{0}}e^{-\kappa s}$, which is quadratic in $\kappa$
for $\bm{r}_{0}\rightarrow\bm{r}$. In this limit the optical length $s$
approaches $n\rho$ where $\rho$ denotes the distance
$\rho=|\bm{r}-\bm{r}_{0}|$. Consequently, we obtain:

$\partial_{r}\partial_{r_{0}}e^{-\kappa
s}\sim\kappa^{2}n^{2}e^{-n\rho\,\kappa}\,.$ (32)

Furthermore, the amplitude ${\cal A}_{0}$ approaches $-(4\pi\nu\rho)^{-1}$,
because near the point of emission the wave equation (18) reduces to the
Laplace equation $\nabla^{2}D=\delta(\bm{r}-\bm{r}_{0})/\nu$ with this
solution. Taking all this into account, we get

$Q\sim\kappa\,\frac{n^{3}}{2\pi}\,\beta_{1}^{\prime}e^{-n\rho_{0}\kappa}$ (33)

where we kept a finite distance $\rho_{0}$ between emitter and receiver that
should be of atomic size.

###  III.3 Anomalous pressure

Equation (33) describes the violation of the Abraham identity for the spectral
stress of the renormalizer. We give it a minus sign to account for the
subtraction in renormalization, $W=-Q$ in Eq. (14), and obtain the result

$(\nabla\cdot\sigma)_{r}=\frac{\varepsilon^{\prime}}{\varepsilon}\,\mathrm{tr}\,\sigma_{E}+\frac{\mu^{\prime}}{\mu}\,\mathrm{tr}\,\sigma_{M}-n^{3}\partial_{r}\frac{p}{n^{3}}$
(34)

for the total stress $\sigma=\sigma_{E}+\sigma_{M}$ with the pressure

$p=\frac{\hbar
c}{(2\pi)^{2}}\int_{0}^{\infty}n^{3}\left(\beta_{1E}+\beta_{1M}\right)e^{-n\rho_{0}\kappa}\kappa\,d\kappa\,.$
(35)

In Eq. (34) we pulled $\varepsilon^{\prime}/\varepsilon$ and
$\mu^{\prime}/\mu$ out of the spectral integral (14) for simplicity. In
reality, $\varepsilon^{\prime}/\varepsilon$ and $\mu^{\prime}/\mu$ depend on
$\kappa$ and therefore should be left inside the integral or regarded as
operators outside. Equation (35) shows that the anomalous pressure depends on
the sum of the scattering amplitudes of the two polarizations. Finally, as
$n^{3}(p/n^{3})^{\prime}=p^{\prime}-(3/2)(\varepsilon^{\prime}/\varepsilon+\mu^{\prime}/\mu)$,
$\mathrm{tr}\,\mathbb{1}=3$, and
$(\nabla\cdot\sigma)_{\theta}=(\nabla\cdot\sigma)_{\phi}=0$, we can cast Eq.
(34) in a remarkably simple form: as the van der Waals anomaly (1) in the
Abraham identity (3).

###  III.4 Cosmological considerations

Without dispersion (frequency dependance of $\varepsilon$ and $\mu$) the
anomalous pressure of Eq. (35) would diverge like $(n\rho_{0})^{-2}$. This
does not occur in ordinary dielectric media where $\varepsilon$ and $\mu$ fall
off like $\kappa^{-2}$ beyond the last resonance of the material. Therefore
the pressure integral (35) does only logarithmically depend on the atomic size
$\rho_{0}$ (and the precise value of the cutoff $\rho_{0}$ does not matter in
comparison with the other contributions to $p$). But now think of the curved
space of general relativity LL2  as a medium Gordon ; Plebanski ; Schleich ;
LeoPhil . The equivalence principle LL2  requires that space acts the same on
all objects, in particular on all frequencies of the electromagnetic field —
up to the scale where classical general relativity does no longer hold. It is
generally believed Case  that this scale is given by the Planck length

$\ell_{\mathrm{P}}=\sqrt{\frac{\hbar G}{c^{3}}}$ (36)

with $G$ being Newton’s gravitational constant. Our analogy would thus suggest
that $n\rho_{0}$ for the medium of space lies in the order of the Planck
length $\ell_{\mathrm{P}}$ of Eq. (36). In this case, the pressure (35) may
have a cosmologically relevant gravitational effect $Gp$, as both $G$ and
$\hbar$ drop out. It depends on the details whether it does.

According to the cosmological principle LL2  space is homogeneous and
isotropic, as empirically verified in astronomical surveys Survey . For
homogeneous and isotropic space, the refractive index is given by LL2

$n=\frac{2n_{1}}{1+k(r/a)^{2}}$ (37)

with constant $n_{1}$ and scale factor $a$ where $a$ gives the radius of
curvature and $k\in\\{-1,0,1\\}$ its sign. Space acts as an impedance–matched
medium according to Eq. (6). We obtain from Eqs. (30) and (37):

$\beta_{1E}=\beta_{1M}=\frac{1}{6n^{3}}\left(n^{\prime\prime}-\frac{n^{\prime}}{r}-\frac{2n^{\prime
2}}{n}\right)=0\,.$ (38)

There is no scattering in maximally symmetric spaces Fish3D , and as our
anomalous pressure is given by the scattering amplitude, there is no anomaly,
provided $a$ is indeed constant.

However, the universe is expanding, the scale factor $a$ is growing with time.
The time–dependence of $a$ generates a dynamical Casimir effect GH  with a
trace anomaly Annals  that also goes as $(n\rho_{0})^{-2}$. If we take for
$n\rho_{0}$ the Planck length $\ell_{\mathrm{P}}$ we obtain the correct order
of magnitude of the cosmological constant Annals . Moreover, a recent
comparison Dror  with astronomical data has found that the anomaly Annals
appears to agree not only on the order of magnitude, but with the details of
the cosmic expansion, as measured by supernova explosions, Baryon Acoustic
Oscillations and the Cosmic Microwave Background. The theory may resolve the
Hubble tension Verde  without invoking new, untested physics DiValentino  (but
rather by revisiting aspects of quantum electromagnetism that have been
neglected for decades). The van der Waals anomaly, although not acting in the
homogeneous and isotropic space of cosmology, is the archetype of this theory.

##  IV Possible experiment

Trace anomalies have been theoretically predicted a considerable time ago Wald
, but they have never been experimentally tested. Here we describe the idea
for an experimental test with Bose–Einstein condensates of alkali atoms
Pethick . At first glance, alkali condensates might not appear as the most
likely candidates for demonstrating delicate dielectric effects, as they are
dilute gases with low density and hence low dielectric response, but the
possibility to optically excite and probe elementary excitations with
well–defined spatial structures may very well outweigh this disadvantage.

###  IV.1 Dilute gases

Let us first discuss some consequences of the diluteness of Bose Einstein
condensates within the context of the van der Waals anomaly that apply to
other dilute gases as well. The renormalized vacuum stresses $\sigma_{E}$ and
$\sigma_{M}$ are caused by the scattering of virtual electromagnetic waves on
inhomogeneities. As these are proportional to the number density $\rho$, they
too must be at least linear in $\rho$, but so are also $\varepsilon^{\prime}$
and $\mu^{\prime}$ in Eq. (3). From this and the van der Waals anomaly (1)
follows in the limit of vanishing density:

$\nabla\cdot(\sigma_{E}+\sigma_{M})\sim-\nabla p\,.$ (39)

Assuming that the stress is isotropic in this limit, we integrate:

$\sigma_{E}+\sigma_{M}\sim-p\mathbb{1}\,.$ (40)

As $p$ enters both $\sigma_{E}$ and $\sigma_{M}$ equally in the van der Waals
anomaly (1) we further assume that the stresses are equal in the dilute limit,
and get

$\sigma_{F}\sim-\frac{p}{2}\mathbb{1}\,.$ (41)

These ideas suggest that in the dilute limit the electric and magnetic
stresses are entirely given by the anomalous pressure and hence by local
dielectric properties (by $\varepsilon$ and $\mu$ and their derivatives up to
second order).

Consider a dilute gas with $\mu=1$ (such as the alkali Bose–Einstein
condensates Pethick ). In the dilute limit we have:

$\varepsilon=1+\chi\,,\quad\chi=\frac{\alpha}{\varepsilon_{0}}\,\rho$ (42)

where $\alpha$ denotes the polarisability Pethick . Neglecting in Eq. (30) all
terms of higher order in $\rho$ we obtain in the planar Grin1  and the
spherical case:

$\beta_{1E}+\beta_{1M}=\frac{\alpha}{6\varepsilon_{0}}\begin{cases}\partial_{z}^{2}\rho&\quad\mbox{(planar)}\\\
\left(\partial_{r}^{2}-\frac{\partial_{r}}{r}\right)\rho&\quad\mbox{(spherical).}\end{cases}$
(43)

In either case, the resulting van der Waals anomaly is linear in density. This
and Eq. (40) confirms the previous, contraversal result Avni  that the vacuum
stress on a homogeneous sphere in a uniform background is linear in the dilute
limit, despite the pairwise nature of van der Waals interactions — and the
very existence of the van der Waals anomaly invalidates the criticism Comment
as a whole. However, the divergence of the stress does not necessarily give
the force density, at least not in mechanical equilibrium Pita .

The force in dielectrics is the Helmholtz force obtained by varying the free
energy with respect to the density, keeping all external charges and currents
constant LL8 . In this way one gets the force density LL8 :

$\bm{f}=\nabla\,\frac{\partial\varepsilon}{\partial\rho}\,\rho\,\frac{\bm{E}\cdot\bm{D}}{2}-\frac{\bm{E}\cdot\bm{D}}{2}\,\frac{\nabla\varepsilon}{\varepsilon}\,.$
(44)

In the dilute limit, the Helmholtz force becomes the dipole force Bjorkholm
with

$\bm{f}=-\rho\,\nabla V$ (45)

in terms of the dipole potential

$V=-\frac{\alpha}{2}\,E^{2}=\frac{\alpha}{\varepsilon_{0}}\,\mathrm{tr}\,\sigma_{E}\,.$
(46)

We assume this expression also for vacuum forces (and integrate over the
entire spectrum of the polarisability as a function of imaginary wavenumber
$\kappa$).

###  IV.2 Estimation

Let us first estimate the effect of the van der Waals anomaly on Bose–Einstein
condensates. Consider a simple toy model: the particle in a box. Suppose the
quantum gas is confined by optical forces to a box of height $a$ (in
$z$–direction) and width and length $b$. Assume further that the gas is made
non–interacting by Feshbach resonance Pethick . One can create excited states
in $z$–direction, for example by applying optical pulses Spectroscopy .
Suppose, for simplicity, that the gas is homogeneous in the other directions.
We thus have for the number density

$\rho=2\rho_{0}\sin^{2}k_{l}z\,,\quad k_{l}=\frac{\pi}{a}\,l$ (47)

with positive integer $l$, and average density

$\rho_{0}=\frac{N}{ab^{2}}$ (48)

where $N$ denotes the total number of atoms. Without atom–atom collisions —
and without van der Waals anomaly — the energies of these states are just the
ones of a quantum particle in a box:

$E_{l}=\frac{\hbar^{2}k_{l}^{2}}{2m}$ (49)

where $m$ denotes the atomic mass. The van der Waals anomaly will create a
perturbation of the energies we calculate by first–order perturbation theory,
averaging the dipole potential of Eq. (46) over the spatial probability
distribution of an individual atom $\rho/N$. We obtain from Eqs. (35), (41),
(43) and (47-49):

$\frac{\langle V\rangle}{E_{l}}=\delta_{0}$ (50)

with

$\delta_{0}=\frac{\rho_{0}}{4\pi\lambda_{C}}\int_{0}^{\kappa_{0}}\left(\frac{\alpha}{\varepsilon_{0}}\right)^{2}\kappa\,d\kappa$
(51)

where $\lambda_{C}$ denotes the Compton wavelength

$\lambda_{c}=\frac{2\pi\hbar}{mc}\,.$ (52)

In the integral (35) over the imaginary wavenumbers $\kappa$ we have replaced
the exponential by a hard cutoff $\kappa_{0}$. We see that the van der Waals
anomaly appears as a density–dependent contribution to the effective mass,
independent of quantum number $l$.

Let us estimate the magnitude of this contribution. The most frequently used
atom for Bose–Einstein condensates is Rubidium. It has a Compton wavelength of
$\lambda_{c}\approx 2.5\times 10^{-18}\mathrm{m}$. As a rough approximation we
assume the polarizibility as being constant and of Pethick
$\alpha/\varepsilon_{0}\approx 300\times 4\pi a_{B}^{3}$ in atomic units (with
Bohr radius $a_{B}$). We assume the cutoff in the optical range:
$\kappa_{0}=2\pi/\lambda_{0}$ with $\lambda_{0}=0.5\mu\mathrm{m}$. Finally we
suppose that $10^{6}$ atoms are confined in a box of $a=100\mu\mathrm{m}$ and
$b=10\mu\mathrm{m}$. This gives $\delta_{0}=10^{-4}$ as a rough estimate of
the relative contribution of the van der Waals anomaly to the energy. The
details depend on the actual polarizibilities as a function of imaginary
frequency. The effect is small, but not out of the range of precision
experiments.

The modified spectrum can be probed by exciting motional states of the trapped
condensate. The excitation lines should vary with density as indicated in Eqs.
(50) and (51). The atom–atom collision energy does also vary with density,
though. We have assumed that the collisions are reduced to zero by Feshbach
resonance Pethick , but a residual interaction will always remain in practice.
However, its contribution to the energy does depend differently on $l$. As the
collision term in the Gross–Pitaevskii equation Pethick  is proportional to
$\rho$, the relative contribution to the energy goes with
$\langle\rho\rangle/E_{l}\propto l^{-2}$ for the density profile of Eq. (47).
The different scaling with quantum number $l$ may serve to discriminate
between the residual atom–atom collisions and the van der Waals anomaly.

###  IV.3 Harmonic trap

Consider now a spherically symmetric harmonic trap of frequency $\Omega$. As
before, atomic collisions shall be switched off by Feshbach resonance Pethick
. Without atom–atom interaction, the spherically–symmetric eigenstates
$\psi_{l}$ obey the Schrödinger equation:

$E_{l}\,\psi_{l}=-\frac{\hbar^{2}}{2m}\left(\partial_{r}^{2}+\frac{2}{r}\,\partial_{r}\right)\psi_{l}+\frac{m\Omega^{2}}{2}\,r^{2}\,\psi_{l}\,.$
(53)

We scale the position $\bm{r}$ in terms of the characteristic length $a$ as

$\bm{r}=a\bm{\xi}\quad\mbox{with}\quad a^{2}=\frac{\hbar}{m\Omega}\,,$ (54)

and obtain the energies

$E_{l}=\hbar\Omega\left(2l+\frac{3}{2}\right)$ (55)

with non–negative integer $l$, and the wave functions

$\psi_{l}=\frac{{\cal
N}_{l}}{\xi}\,\exp\left(-\frac{\xi^{2}}{2}\right)\,H_{2l+1}(\xi)$ (56)

in terms of the Hermite polynomials Erdelyi . The wave functions are
normalized to unity with respect to $\bm{\xi}$ by

${\cal N}_{l}^{2}=\frac{1}{\pi^{3/2}\,2^{2l+2}\,(2l+1)!}\,.$ (57)

Writing the number density $\rho$ of the condensate as

$\rho=\rho_{0}\,\psi_{l}^{2}$ (58)

with average density

$\rho_{0}=\frac{N}{a^{3}}$ (59)

we obtain from Eq. (43) in the spherical case the result:

$\frac{\langle
V\rangle}{\hbar\Omega}=\delta_{0}\int_{0}^{\infty}\left(\partial_{\xi}^{2}\psi_{l}^{2}-\frac{\partial_{\xi}\psi_{l}^{2}}{\xi}\right)4\pi\,\xi^{2}\,d\xi$
(60)

with $\delta_{0}$ given by Eqs. (51) and (52). The integral is dominated by
small $\xi$. We apply the asymptotics of the Hermite polynomials Erdelyi
$H_{2l+1}(\xi)$ for large index and small variable,

$\psi_{l}\sim(-1)^{l}\,\frac{\sin(\sqrt{4l+3}\,\xi)}{(4l+2)^{1/4}\pi\xi}\,,$
(61)

and arrive at an approximation of Eq. (60) that is quite accurate for $l\geq
1$:

$\frac{\langle
V\rangle}{\hbar\Omega}\sim\frac{\delta_{0}}{3\pi^{2}}\,\frac{(4l+3)^{3/2}}{2l+1}\,.$
(62)

For $l=0$ we have

$\frac{\langle
V\rangle}{\hbar\Omega}=\frac{\delta_{0}}{\sqrt{2}\,\pi^{3/2}}\,.$ (63)

In the harmonic trap, the relative contribution $\langle V\rangle/E_{l}$ of
the van der Waals anomaly to the energy varies with $l$ as
$(4l+3)^{1/2}/(2l+1)\sim l^{-1/2}$ in contrast to the condensate in a box
where $\langle V\rangle/E_{l}=\mathrm{const}$. This is because
harmonic–oscillator states are not rigidly confined, but reach to the
classical turning points that depend on the energy, which reduces the density.
We can also work out the effect of collisions: the collision term in the
Gross–Pitaevskii equation Pethick  scales with

$\displaystyle\langle\rho\ \rangle$ $\displaystyle\sim$
$\displaystyle\frac{\rho_{0}}{2\pi^{2}}\,\frac{(4l+3)^{1/2}}{2l+1}\quad\mbox{for}\quad
l\geq 1\,,$ (64) $\displaystyle\langle\rho\ \rangle$ $\displaystyle=$
$\displaystyle\frac{\rho_{0}}{(2\pi)^{3/2}}\quad\mbox{for}\quad l=0\,.$ (65)

Like for the condensate in a box, the ratio $\langle
V\rangle/\langle\rho\rangle$ goes with the energy $E_{l}$, which allows to
discriminate the van der Waals anomaly from residual collisions.

Overall we found the effect of the anomaly to be in the range of $10^{-4}$ of
the energy. The van der Waals anomaly can be distinguished from the influence
of residual collisions, and it can be measured by probing the transition
frequencies of a trapped Bose–Einstein condensate, following Schawlow’s motto
for precision measurements Bres : “never measure anything but frequency”.

## Acknowledgements

We thank Y. Avni, D. Berechya, N. Davidson, N. Ebel, E. Efrati, I. Griniasty,
and E. Shahmoon for discussions. The paper has been supported by the Israel
Science Foundation and the Murray B. Koffler Professorial Chair.

##  Appendix A Green function and stress

In this appendix we derive the expression for the Abraham stress as a function
of the scalar Green functions $g_{E}$ and $g_{M}$. Throughout, we use the
vector spherical harmonic notation introduced and discussed in Appendix E.

###  A.1 Scalar Green functions

We start with $\mathrm{G}(\bm{r},\bm{r}_{0})$, the matrix Green function (the
Green bi–tensor) defined as the solution of the inhomogeneous wave equation
LeoBook :

$\nabla\times\frac{1}{\mu}\nabla\times\mathrm{G}+\varepsilon\kappa^{2}\mathrm{G}=\mathbb{1}\,\delta(\bm{r}-\bm{r}_{0})$
(66)

for imaginary wavelengths $\kappa$. Equation (66) is solved by decomposing the
matrix Green function into two contributions $\mathrm{G}_{E}$ and
$\mathrm{G}_{M}$ that originate from the two electromagnetic polarisations
preserved in spherically symmetric media, and a contact term $\mathrm{G}_{C}$:

$\mathrm{G}=\mathrm{G}_{E}+\mathrm{G}_{M}+\mathrm{G}_{C}$ (67)

where each polarisation corresponds to one of the scalar Green functions
$g_{E}$ and $g_{M}$ as

$\displaystyle\mathrm{G}_{E}$ $\displaystyle=$
$\displaystyle-\sum_{lm}\frac{g^{lm}_{E}(r,r_{0})}{l(l+1)}\,\bm{\Phi}_{lm}(\theta,\phi)\otimes\bm{\Phi}^{*}_{lm}(\theta_{0},\phi_{0})\,,$
$\displaystyle\mathrm{G}_{M}$ $\displaystyle=$
$\displaystyle\frac{1}{\kappa^{2}\varepsilon(r)\varepsilon(r_{0})}\sum_{lm}\nabla\times\frac{g^{lm}_{M}(r,r_{0})}{l(l+1)}\,\bm{\Phi}_{lm}(\theta,\phi)$
(68)
$\displaystyle\quad\quad\quad\quad\quad\quad\otimes\,\bm{\Phi}^{*}_{lm}(\theta_{0},\phi_{0})\times\overleftarrow{\nabla}_{0}$

in terms of the vector spherical harmonics $\bm{\Phi}_{lm}(\theta,\phi)$
described in Appendix E. We are going to specify the contact term later on.
Let us first focus on the component $\mathrm{G}_{E}$ of the $E$–polarization.
Consider Eq. (66) applied to $\mathrm{G}_{E}$ only. From the curl identities
of the vector spherical harmonics [Eqs. (3.12) of Barrera et al. Barrera ]
follows

$\nabla\times\frac{1}{\mu}\nabla\times
f(r)\bm{\Phi}_{lm}=\frac{1}{r}\partial_{r}\frac{1}{\mu}\partial_{r}rf(r)\,\bm{\Phi}_{lm}$
(69)

for an arbitrary function $f(r)$. We see that for $\mathrm{G}_{E}$ the
left–hand side of the inhomogeneous wave equation (66) consists entirely of
$\bm{\Phi}_{lm}\otimes\bm{\Phi}^{*}_{lm}$ terms. We decompose the right–hand
side into vector spherical harmonics:

$\displaystyle\delta(\bm{r}-\bm{r}_{0})\mathbb{1}$ $\displaystyle=$
$\displaystyle\frac{\delta(r-r_{0})}{r^{2}}\sum_{lm}\bigg{(}\bm{Y}_{lm}(\theta,\phi)\otimes\bm{Y}^{*}_{lm}(\theta_{0},\phi_{0})$
(70)
$\displaystyle\quad\quad+\frac{\bm{\Phi}_{lm}(\theta,\phi)\otimes\bm{\Phi}^{*}_{lm}(\theta_{0},\phi_{0})}{l(l+1)}$
$\displaystyle\quad\quad+\frac{\bm{\Psi}_{lm}(\theta,\phi)\otimes\bm{\Psi}^{*}_{lm}(\theta_{0},\phi_{0})}{l(l+1)}\bigg{)}$

and see that Eq. (66) is satisfied for the
$\bm{\Phi}_{lm}\otimes\bm{\Phi}^{*}_{lm}$ terms if we require:

$\frac{1}{r}\partial_{r}\frac{1}{\nu}\,\partial_{r}rg^{lm}-\frac{l(l+1)}{r^{2}\nu}g^{lm}-\frac{n^{2}\kappa^{2}}{\nu}\,g^{lm}=\frac{\delta(r-r_{0})}{r^{2}}$
(71)

dropping the polarization label and using $\nu=\mu$ according to Eq. (17) and
the refractive index $n$ of Eq. (19). This settles the electric polarization.

Now turn to the magnetic polarization. It should take care of the remaining
$\bm{Y}_{lm}\otimes\bm{Y}^{*}_{lm}$ and
$\bm{\Psi}_{lm}\otimes\bm{\Psi}^{*}_{lm}$ terms in Eq. (70) for
$\delta(\bm{r}-\bm{r}_{0})\mathbb{1}$. Let us define the contact term

$\mathrm{G}_{C}=\frac{\delta(r-r_{0})}{\kappa^{2}rr_{0}\,\varepsilon(r_{0})}\sum_{lm}\left(\bm{Y}_{lm}\otimes\bm{Y}^{*}_{lm}+\frac{\bm{\Psi}_{lm}\otimes\bm{\Psi}^{*}_{lm}}{l(l+1)}\right)$
(72)

where we omit the dependancies on the various angles for simplicity of
writing. The contact term does not contribute to the stress for
$\bm{r}_{0}\neq\bm{r}$, but it is required for mathematical consistency as we
shall see below. When substituted in Eq. (66)
$\varepsilon\kappa^{2}\mathrm{G}_{C}$ immediately gives the
$\bm{Y}_{lm}\otimes\bm{Y}^{*}_{lm}$ and
$\bm{\Psi}_{lm}\otimes\bm{\Psi}^{*}_{lm}$ terms for
$\delta(\bm{r}-\bm{r}_{0})\mathbb{1}$. Furthermore, we obtain from the curl
identities [Eqs. (3.12) of Barrera et al. Barrera ]:

$\displaystyle\nabla\times\mathrm{G}_{C}$ $\displaystyle=$
$\displaystyle\frac{1}{\kappa^{2}\varepsilon(r_{0})}\sum_{lm}\bm{\Phi}_{lm}\otimes\bigg{(}-\frac{\delta(r-r_{0})}{r^{2}r_{0}}\bm{Y}_{lm}^{*}$
(73)
$\displaystyle\quad\quad\quad\quad\,+\frac{1}{r}\partial_{r}\frac{\delta(r-r_{0})}{r_{0}}\frac{\bm{\Psi}_{lm}^{*}}{l(l+1)}\bigg{)}$

and writing $r^{-1}\partial_{r}r_{0}^{-1}\delta(r-r_{0})$ as
$-r_{0}^{-1}\partial_{r_{0}}\delta(r-r_{0})$ and using the curl identities
once more, but this time for $\bm{r}_{0}$:

$\nabla\times\mathrm{G}_{C}=\frac{\delta(r-r_{0})}{\kappa^{2}rr_{0}\,\varepsilon(r_{0})}\sum_{lm}\bm{\Phi}_{lm}\otimes\bm{\Phi}_{lm}^{*}\times\overleftarrow{\nabla}_{0}\,.$
(74)

Requiring Eq. (71) for $g_{M}^{lm}$ (with $\nu=\varepsilon$) we obtain along
similar lines as in the case of the electric polarization:

$\nabla\times\frac{1}{\mu}\nabla\times(\mathrm{G}_{M}+\mathrm{G}_{C})+\varepsilon\kappa^{2}\mathrm{G}_{M}=0\,.$
(75)

This proves that, in Eq. (66) for the Green tensor, $\mathrm{G}_{M}$ plus
$\mathrm{G}_{C}$ generates the remaining terms in Eq. (70) for the right–hand
side.

Equation (71) thus defines the scalar Green functions for both polarizations.
As Eq. (71) does not depend on $m$ the Green function components $g^{lm}$ do
not depend on $m$ either. The scalar Green functions are given in terms of the
components $g^{lm}$ with respect to the standard spherical harmonics. In
position space they read as

$g(\bm{r},\bm{r_{0}})=\sum_{lm}g^{lm}(r,r_{0})\,Y_{lm}(\theta,\phi)\,Y^{*}_{lm}(\theta_{0},\phi_{0})$
(76)

and the scalar wave equation (71) takes the form of Eq. (18) with definition
(16).

###  A.2 Electric and magnetic stresses

The fluctuation–dissipation theorem Scheel  relates the vacuum stress to the
matrix Green function. In particular, defining

$\displaystyle\tau_{E}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\langle
0|\bm{D}(\bm{r})\otimes\bm{E}(\bm{r}_{0})+\bm{D}(\bm{r}_{0})\otimes\bm{E}(\bm{r})|0\rangle\,,$
$\displaystyle\tau_{M}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\langle
0|\bm{B}(\bm{r})\otimes\bm{H}(\bm{r}_{0})+\bm{B}(\bm{r}_{0})\otimes\bm{H}(\bm{r})|0\rangle$
(77)

we arrive at the expressions Scheel ; LeoBook :

$\displaystyle\tau_{E}$ $\displaystyle=$ $\displaystyle-\frac{\hbar
c}{\pi}\int_{0}^{\infty}\kappa^{2}\varepsilon\,(\mathrm{G}_{E}+\mathrm{G}_{M})\,d\kappa\,,$
$\displaystyle\tau_{M}$ $\displaystyle=$ $\displaystyle\frac{\hbar
c}{\pi}\int_{0}^{\infty}\frac{1}{\mu}\nabla\times(\mathrm{G}_{E}+\mathrm{G}_{M})\times\overleftarrow{\nabla}_{0}\,d\kappa\,.$
(78)

The stress tensors are then given by

$\sigma_{F}=\tau_{F}-\frac{1}{2}\,\mathrm{tr}\,\tau_{F}\,\mathbb{1}$ (79)

in the limit $\bm{r}_{0}\rightarrow\bm{r}$. In order to derive manageable
expressions for the stress–tensor contributions, we write them in the spectral
representation (14) for the components $i$, fields $F$ and polarizations $P$
as

$W^{iP}_{iF}=\sum_{lm}{\cal
W}_{iF}^{iP}\,Y_{lm}(\theta,\phi)Y^{*}_{lm}(\theta_{0},\phi_{0})\,.$ (80)

The sum of the angular coefficients ${\cal W}_{iF}^{iP}$ over $P$ corresponds
to the Fourier component ${\cal W}$ in the planar case Grin1 . To illustrate
our method of determining the ${\cal W}_{iF}^{iP}$ consider the simplest case:
${\cal W}_{iE}^{iE}$. We get for the corresponding part of the correlation
function $\tau_{E}$:

$\displaystyle\tau_{E}^{(E)}$ $\displaystyle=$ $\displaystyle-\frac{\hbar
c}{\pi}\int_{0}^{\infty}\kappa^{2}\varepsilon\,\mathrm{G}_{E}\,d\kappa$ (81)
$\displaystyle=$ $\displaystyle\frac{\hbar
c}{\pi}\sum_{lm}\int_{0}^{\infty}\frac{\kappa^{2}\varepsilon
g^{lm}_{E}}{l(l+1)}\,\bm{\Phi}_{lm}\otimes\bm{\Phi}^{*}_{lm}$

from Eq. (68) and assuming the limit $\theta_{0}\rightarrow\theta$ and
$\phi_{0}\rightarrow\phi$. Then we use the fact that $g^{lm}$ is independent
of $m$ and apply the last one of the sum identities (140). From this and Eq.
(79) follows

${\cal W}_{rE}^{rE}=\kappa^{2}\varepsilon\,g_{E}^{lm}\,.$ (82)

For calculating ${\cal W}_{iM}^{iE}$ we define the corresponding
$\tau_{M}^{(E)}$ according to Eq. (78), use Eq. (68), apply the curl
identities of the vector spherical harmonics [Eqs. (3.12) of Barrera et al.
Barrera ] and finally use the sum identities (140). We obtain from Eq. (79):

$\displaystyle{\cal W}_{rM}^{rE}$ $\displaystyle=$
$\displaystyle-\frac{\partial_{r_{0}}r_{0}\,\partial_{r}r\,g_{E}^{lm}}{\mu
rr_{0}}-{\cal W}_{\theta M}^{\theta E}\,,$ $\displaystyle{\cal W}_{\theta
M}^{\theta E}$ $\displaystyle=$ $\displaystyle{\cal W}_{\phi M}^{\phi
E}=-\frac{l(l+1)}{\mu rr_{0}}\,g^{lm}_{E}\,.$ (83)

Equations (82) and (83) give the expressions in Eq. (15) when represented in
position space according to Eqs. (76) and (80). In order to determine the
components for the magnetic polarization we may proceed similarly, but
applying also the wave equation (71) to reduce the number of derivatives
MagRemark . Alternatively, we may just take advantage of the electromagnetic
duality and interchange the fields and the polarizations $E\leftrightarrow M$
and $\mu\leftrightarrow\varepsilon$.

##  Appendix B Geometrical optics

In this appendix we provide some more details on the geometrical–optics theory
for the renormalizing scalar Green functions $D_{0}$ and $D_{1}$. We closely
follow the planar case Grin1  where the general expressions have been derived,
and state the results for the spherical case.

###  B.1 Quadratic expansion

As renormalization is local, the renormalizing Green function $D$ should
depend on the local $\varepsilon$ and $\mu$ and their derivatives. We
postulate that these are maximally second derivatives. This means that we
expand the dielectric functions up to second order around the point of
emission $\bm{r}_{0}$:

$n=n_{0}+n_{0}^{\prime}(r-r_{0})+\frac{n_{0}^{\prime\prime}}{2}(r-r_{0})^{2}$
(84)

with $n_{0}=n(r_{0})$, and similarly for $\varepsilon$ and $\mu$. The
renormalizing Green function thus perceives the medium as if the dielectric
functions depend on the point of emission at $r_{0}$, which potentially breaks
reciprocity.

###  B.2 Optical length

The central quantity of geometrical optics is the optical length $s$ that
satisfies the eikonal equation:

$(\nabla s)^{2}=n^{2}\,.$ (85)

Spherical symmetry implies that $s$ depends on $r$, $r_{0}$ and the angle
$\gamma$ between the origin, $\bm{r}$ and $\bm{r}_{0}$. The distance between
$\bm{r}$ and $\bm{r}_{0}$ is given by:

$\rho=\sqrt{r^{2}+r_{0}^{2}-2rr_{0}\cos\gamma}\,.$ (86)

One verifies that the solution of Eq. (85) is:

$\displaystyle s$ $\displaystyle=$
$\displaystyle\rho\left(n_{0}+\frac{n_{0}^{\prime}}{2}(r-r_{0})+\frac{n_{0}^{\prime\prime}}{6}(r-r_{0})^{2}\right.$
(87)
$\displaystyle\left.\quad\quad-\frac{n_{0}^{\prime}r_{0}\,(n_{0}^{\prime}r_{0}+2n_{0})}{24n_{0}}\,\sin^{2}\gamma\right)$

up to quadratic order in $r-r_{0}$ and $\gamma$. One also verifies the
reciprocity of $s$: if we replace $r\leftrightarrow r_{0}$ and
$n\leftrightarrow n_{0}$ in Eq. (87) and use expansion (84) we obtain the same
$s$, apart from third–order corrections. So despite the fact that $n$ depends
on $r_{0}$, the optical length $s$ is reciprocal up to quadratic order.

###  B.3 Amplitude

The amplitude ${\cal A}_{0}$ of the wave in geometrical optics satisfies the
continuity equation Grin1 :

$\nabla\cdot\left(\nu{\cal A}_{0}^{2}\nabla s\right)=0$ (88)

with $\nu$ defined in Eq. (17). One verifies that up to quadratic order the
solution of Eq. (88) is the same as in the planar case Grin1 :

${\cal
A}_{0}=-\frac{1}{4\pi\sqrt{\nu_{0}\nu}}\left(\frac{1}{\rho}+\frac{n_{0}^{2}R_{0}}{48}\rho\right)$
(89)

where $R_{0}$ denotes the 3D curvature scalar LeoPhil :

$R_{0}=-\frac{4\nabla^{2}n_{0}}{n_{0}^{3}}+\frac{2n_{0}^{\prime
2}}{n_{0}^{4}}\,.$ (90)

The only difference is that in the planar case Grin1
$\nabla^{2}=\partial_{z}^{2}$ whereas in the spherical case we have
$\nabla^{2}=\partial_{r}^{2}+(2/r)\partial_{r}$. One verifies again that
${\cal A}_{0}$ is reciprocal up to quadratic order.

###  B.4 First scattering

As $s$ and ${\cal A}_{0}$ are reciprocal, the outgoing wave $D_{0}$ defined in
Eq. (28) is reciprocal, too. The contribution $D_{1}$ of Eq. (29) depends on
$s$ and ${\cal A}_{0}$ as well, but also on the amplitude $\beta_{1}$. This
quantity describes the first scattering in inhomogeneities around the point of
emission Grin1 . The scattering amplitude is given by Grin1 :

$\beta_{1}=\frac{\nabla\cdot\nu\nabla{\cal A}_{0}}{2n^{2}\nu{\cal
A}_{0}}-\frac{R_{r}^{r}}{4}$ (91)

with $R_{r}^{r}$ defined in Eq. (20). For $\nu=n$ when the medium establishes
an exact geometry, $R_{r}^{r}$ is the radial component of the 3D Ricci tensor
LeoPhil . Evaluating Eq. (91) we obtain Eq. (30) for $\beta_{1}$. One verifies
that $D=D_{0}+D_{1}$ with the expressions given satisfies the wave equation
(18) around $\bm{r}_{0}$ within $\mathrm{O}(\kappa^{-1})$. The scattering
amplitude $\beta_{1}$ depends entirely on $r_{0}$ and therefore clearly breaks
reciprocity, which causes the van der Waals anomaly.

##  Appendix C Asymptotics

In order to see whether the renormalizing Green functions remove the infinity
of the vacuum stress, we need to identify and characterize this divergency in
the first place. The stress is determined in Eq. (14) by an integration over
spectral stresses, and those are given in Eq. (15) by maximally second
derivatives of the scalar Green functions $g_{P}$. We thus need to analyse the
asymptotics of the $g_{P}(\bm{r},\bm{r}_{0})$ for large $\kappa$ and
$\bm{r}_{0}\sim\bm{r}$. The Green functions $g_{P}$ we represent in Eq. (76)
by the Green coefficients $g^{lm}$ dropping from now on the polarization label
$P$. The required asymptotics of $g$ is then given by the asymptotics of the
$g^{lm}$ for large $l$ and large $\kappa$. This is the asymptotics we
determine in this appendix.

###  C.1 Langer transformation

Similar to the planar case Grin1  we apply a WKB technique. However, for the
spherical case, it is wise to adopt the Langer transformation Langer

$r=e^{x}\,,\quad g^{lm}=e^{-x/2-x_{0}/2}\,u(x)$ (92)

dropping in $u$ the $lm$–indices for simplicity. We have put the constant
$x_{0}/2$ in the exponent for keeping the reciprocity of $g^{lm}$ in $u$.
Equation (71) appears then in the form:

$\partial_{x}\frac{1}{\nu}\,\partial_{x}u-\frac{\nu^{\prime}u}{2\nu^{2}}-\frac{p^{2}+e^{2x}\kappa^{2}n^{2}}{\nu}u=\delta(x-x_{0})$
(93)

where the primes denote differentiations with respect to $x$ (such that
$\nu^{\prime}=r\partial_{r}\nu$) and the parameter $p$ abbreviates

$p=l+\frac{1}{2}$ (94)

that should not be confused with the pressure. In the following we determine
the solution of Eq. (93) for large $\kappa$ and $p$.

###  C.2 Wronskian representation

First, we write down an expression for the exact solution in terms of the
homogeneous solutions $h_{\pm}$ of Eq. (93):

$u(x,x_{0})=\frac{1}{W}\begin{cases}h_{+}(x)h_{-}(x_{0}):&x>x_{0}\\\
h_{+}(x_{0})h_{-}(x):&x<x_{0}\end{cases}$ (95)

where $W$ denotes the Wronskian:

$W=\frac{h^{\prime}_{+}(x)\,h_{-}(x)-h^{\prime}_{-}(x)\,h_{+}(x)}{\nu}\,.$
(96)

The Wronskian is constant, as a consequence of Eq. (93). Since the $h_{\pm}$
solve Eq. (93) for $r\neq r_{0}$, $u$ does this as well. One verifies that the
jump at $x=x_{0}$ specified in Eq. (95) generates the delta function on the
right hand side of Eq. (93), which proves that $u$ solves Eq. (93). The
$h_{\pm}$ need to be chosen such that they decay at $\pm\infty$, respectively,
which determines $u$ and hence the Green coefficient $g^{lm}$ uniquely.

###  C.3 WKB asymptotics

Now we use the standard WKB technique for the homogeneous solutions $h_{\pm}$
in the limit $\kappa,p\rightarrow\infty$. For convenience, we replace $\kappa$
and $p$ by $\kappa/q$ and $p/q$ in terms of the formal parameter $q$ we regard
as small. We write

$h_{\pm}=\exp\left(-\frac{1}{q}\sum_{m=0}^{\infty}q^{m}s_{m}\right)$ (97)

and solve, order by order, the homogeneous equation

$\partial_{x}\frac{1}{\nu}\,\partial_{x}h-\frac{\nu^{\prime}h}{2\nu^{2}}-\frac{p^{2}+e^{2x}\kappa^{2}n^{2}}{q^{2}\nu}u=0\,.$
(98)

Inserting our ansatz (97) into Eq. (98), we get the relation:

$\displaystyle\frac{1}{q^{2}}\sum_{m=0}^{\infty}q^{m}\sum_{k=0}^{m}\frac{s^{\prime}_{k}s^{\prime}_{m-k}}{\nu}+\frac{1}{q}\sum_{m=0}^{\infty}q^{m}\left(-\frac{1}{\nu}s^{\prime\prime}_{m}+\frac{\nu^{\prime}}{\nu^{2}}s^{\prime}_{m}\right)$
$\displaystyle=\frac{p^{2}+e^{2x}\kappa^{2}n^{2}}{q^{2}\nu}+\frac{\nu^{\prime}}{2\nu^{2}}\,.$
(99)

In lowest order ($q^{-2}$) we have

$s_{0}^{\prime}=\pm\sqrt{\kappa^{2}n^{2}e^{2x}+p^{2}}$ (100)

where the $\pm$ corresponds to the two choices $h_{\pm}$ that should and do
vanish for $x\rightarrow\pm\infty$ according to Eq. (97). In the second–lowest
order ($q^{-1}$) we have

$s_{1}^{\prime}=\frac{1}{2s_{0}^{\prime}}\left(s^{\prime\prime}_{0}-\frac{\nu^{\prime}}{\nu}s^{\prime}_{0}\right),$
(101)

whereas in the next order ($q^{0}$) we get

$s^{\prime}_{2}=\frac{1}{2s^{\prime}_{0}}\left(s^{\prime\prime}_{1}-\frac{\nu^{\prime}}{\nu}s^{\prime}_{1}-{s_{1}^{\prime}}^{2}-\frac{\nu^{\prime}}{2\nu}\right).$
(102)

For all other orders, Eq. (99) reduces to

$s^{\prime}_{m}=\frac{1}{2s^{\prime}_{0}}\left(s^{\prime\prime}_{m-1}-\frac{\nu^{\prime}}{\nu}s^{\prime}_{m-1}-\sum_{k=1}^{m-1}s^{\prime}_{k}s^{\prime}_{m-k}\right)$
(103)

that includes Eq. (101) as a special case. Note that Eq. (103) is exactly the
same recurrence relation as in the planar case Grin1 . Note also that the
$\pm$ in Eq. (100) does not change the sign of $s_{1}^{\prime}$ in Eq. (101).
The $\pm$ sign carries then over from $s_{0}^{\prime}$ to $s_{2}^{\prime}$ in
Eq. (102). In fact, all even orders carry the $\pm$ sign, while all odd orders
are unchanged, as one sees from Eq. (103). In terms of our ansatz (97) we thus
have

$h_{\pm}=e^{\mp s_{\cal E}}e^{-s_{\cal O}}$ (104)

where the $s_{\cal E}$ and $s_{\cal O}$ collect the even and orders in $s$,
respectively:

$s_{\cal E}\equiv\frac{1}{q}\sum_{m=0}^{\infty}q^{2m}s_{2m}\,,\quad s_{\cal
O}\equiv\frac{1}{q}\sum_{m=0}^{\infty}q^{2m+1}s_{2m+1}.$ (105)

Given our result (104) for $h_{\pm}$, we write down the Wronskian of Eq. (96):

$W=-\frac{2}{\nu}\,s_{\cal E}^{\prime}\,e^{-2s_{\cal O}}.$ (106)

As $W=\mathrm{const}$, we may put $W=\sqrt{W(x)W(x_{0})}$ and have

$W=-2\sqrt{\frac{s_{\cal E}^{\prime}(x)s_{\cal
E}^{\prime}(x_{0})}{\nu(x)\nu(x_{0})}}\,e^{-s_{\cal O}(x)-s_{\cal
O}(x_{0})}\,,$ (107)

which gives the compact expression

$u=-\frac{1}{2}\sqrt{\frac{\nu(e^{x})\nu(e^{x_{0}})}{s_{\cal
E}^{\prime}(x)s_{\cal
E}^{\prime}(x_{0})}}\exp\left(-\absolutevalue{\int_{x_{0}}^{x}s_{\cal
E}^{\prime}\,dx}\right)$ (108)

that only depends on the even orders in the asymptotic expansion (97). Formula
(108) is explicitly reciprocal, and so is the Green function given by the
Langer transformation, Eq. (92). Note that Eq. (108) is identical with the
expression for the Fourier–transformed Green functions in the planar case
Grin1 . However, the starting point of the recurrence and the expression for
$s_{2}$, Eqs. (100) and (102), are different.

###  C.4 Divergency of the stress

Now we are going to characterize the divergency of the vacuum stress using as
cutoff $\Lambda=1/q$ and the expressions derived for the scalar Green
function. But first we note that due to dispersion LL8  the refractive index
and impedance have the following asymptotic form in the high–frequency limit:

$n(r;\kappa)\sim 1+\frac{n_{\infty}(r)}{\kappa^{2}}\,,\quad Z(r;\kappa)\sim
1+\frac{Z_{\infty}(r)}{\kappa^{2}}$ (109)

where the impedance is defined as $Z=\sqrt{\mu/\varepsilon}$.

In the following we express our results in terms of the spectral stresses
defined in Eq. (80). We take dispersion into account according to Eq. (109)
and calculate the bare stresses given by Eqs. (82) and (83) using our WKB
technique. We obtain for ${\cal W}^{r}_{r}$ for each order of divergency:

$\displaystyle\Lambda^{4}$ $\displaystyle:$
$\displaystyle-\frac{4pw}{r^{3}}\,,$ $\displaystyle\Lambda^{2}$
$\displaystyle:$
$\displaystyle\frac{p\left(-8r^{2}w^{4}n_{\infty}+p^{4}+w^{4}\right)}{2r^{3}w^{5}}\,,$
$\displaystyle\ln\Lambda$ $\displaystyle:$
$\displaystyle\frac{p}{32r^{3}w^{11}\left(w^{2}-p^{2}\right)}\Bigl{(}32p^{4}r^{3}w^{6}n_{\infty}^{\prime}$
(110)
$\displaystyle-16r^{2}w^{4}n_{\infty}\left(w^{2}-p^{2}\right)\left(5p^{4}+w^{4}\right)$
$\displaystyle-\left(w^{2}-p^{2}\right)^{2}\bigl{(}105p^{6}-63p^{4}w^{2}+7p^{2}w^{4}-w^{6}\bigr{)}$
$\displaystyle-\,64p^{2}r^{4}w^{8}n_{\infty}^{2}\Bigr{)}$

with the abbreviation

$w=\sqrt{p^{2}+\kappa^{2}r^{2}}\,.$ (111)

We obtain for ${\cal W}^{\theta}_{\theta}$:

$\displaystyle\Lambda^{4}$ $\displaystyle:$
$\displaystyle\frac{2p^{3}}{r^{3}w}\,,$ $\displaystyle\Lambda^{2}$
$\displaystyle:$
$\displaystyle\frac{-8p^{3}r^{2}w^{4}n_{\infty}+5p^{7}-6p^{5}w^{2}+p^{3}w^{4}-2pw^{6}}{4r^{3}w^{7}}\,,$
$\displaystyle\ln\Lambda$ $\displaystyle:$
$\displaystyle\frac{p}{64r^{3}w^{13}\left(w^{2}-p^{2}\right)}\Bigl{(}-(w^{2}-p^{2})^{2}(1155p^{8}$
(112) $\displaystyle\quad-1617p^{6}w^{2}+553p^{4}w^{4}-47p^{2}w^{6}+4w^{8})$
$\displaystyle+16r^{2}w^{4}(2p^{4}rw^{2}(rw^{2}n_{\infty}^{\prime\prime}-5(w^{2}-p^{2})n^{\prime})$
$\displaystyle\quad+n(35p^{8}-65p^{6}w^{2}+33p^{4}w^{4}-5p^{2}w^{6}+2w^{8})$
$\displaystyle\quad+4p^{2}r^{2}(2w^{2}-3p^{2})w^{4}n_{\infty}^{2})\Bigr{)}\,.$

These are the divergencies of the vacuum stress. In Appendix D we prove that
they are removed by renormalization.

##  Appendix D Convergence

In this appendix we prove that the renormalization converges. By this we mean
that the infinite vacuum stress calculated from the renormalizer, $D_{0}$
+$D_{1}$, compensates for the infinite vacuum stress calculated using the
asymptotics of the Green function (Appendix C). What is left is finite.

###  D.1 Spherical harmonics decomposition

In order to compare the vacuum stress of the renormalizer with the bare stress
given by the asymptotics of the Green function, we need to represent $D_{0}$
and $D_{1}$ in terms of spherical harmonics. Since we are only interested in
the behavior up to quadratic order, we approximate and decompose $D_{0}$ and
$D_{1}$ according to Eqs. (28), (29), (84), (87) and (89):

$\displaystyle D_{0}$ $\displaystyle=$ $\displaystyle{\cal A}_{0}\,e^{-\kappa
s}$ (113) $\displaystyle\sim$
$\displaystyle-\frac{1}{4\pi\rho\sqrt{\nu\nu_{0}}}\left[1+\frac{n_{0}^{2}R_{0}}{48}\rho^{2}\right](1+\kappa\rho\alpha_{0}\gamma^{2})\,e^{-\kappa\rho\chi}$
$\displaystyle=$ $\displaystyle\sum_{m=1}^{4}D_{0}^{m}\quad\mbox{where}$
$\displaystyle\alpha_{0}$ $\displaystyle=$
$\displaystyle\frac{n_{0}^{\prime}r_{0}\,(n_{0}^{\prime}r_{0}+2n_{0})}{24n_{0}}\,,$
(114) $\displaystyle\chi$ $\displaystyle=$ $\displaystyle
n_{0}+\frac{n_{0}^{\prime}}{2}(r-r_{0})+\frac{n_{0}^{\prime\prime}}{6}(r-r_{0})^{2}\,,$
(115) $\displaystyle D^{1}_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi\rho\sqrt{\nu\nu_{0}}}\,e^{-\kappa\rho\chi}\,,$
(116) $\displaystyle D^{2}_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi\sqrt{\nu\nu_{0}}}\frac{n_{0}^{2}R_{0}}{48}\rho\,e^{-\kappa\rho\chi}\,,$
(117) $\displaystyle D^{3}_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi\sqrt{\nu\nu_{0}}}\kappa\alpha_{0}\gamma^{2}\,e^{-\kappa\rho\chi}\,,$
(118) $\displaystyle D^{4}_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi\sqrt{\nu\nu_{0}}}\frac{n_{0}^{2}R_{0}}{48}\rho^{2}\kappa\alpha_{0}\gamma^{2}\,e^{-\kappa\rho\chi}\,,$
(119) $\displaystyle D_{1}$ $\displaystyle=$ $\displaystyle{\cal
A}_{0}\frac{\beta_{1}s}{\kappa}e^{-\kappa
s}\sim-\frac{1}{4\pi\sqrt{\nu\nu_{0}}}\frac{\beta_{1}}{\kappa}\,e^{-\kappa\rho\chi}\,.$
(120)

These terms generate all divergencies in the vacuum stress, with further terms
being of order $\mathrm{O}(\Lambda^{-2})$ or below in their contribution to
the stress.

Equations (113-120) reduce our problem to figuring out how to decompose the
functions $f$ and $\gamma^{2}f_{\alpha}$ into spherical harmonics with

$f_{\alpha}(\bm{r};\bm{r}_{0})=\frac{\rho^{\alpha}}{4\pi}e^{-k\rho}$ (121)

where $k$ abbreviates $\kappa\chi$ here. The $f_{\alpha}$ depend via $\rho$ on
the angles $\theta$ and $\phi$ of the spherical coordinates, whereas $k$ does
not depend on them. The spherical harmonics $f_{-1,lm}$ of $f_{-1}$ are given
by the well–known expression Erdelyi :

$f_{-1,lm}(r,r_{0};k)=\frac{1}{\sqrt{rr_{0}}}\begin{cases}K_{p}(kr_{0})I_{p}(kr)&r<r_{0}\\\
I_{p}(kr_{0})K_{p}(kr)&r>r_{0}\end{cases}$ (122)

where $I_{p}$ and $K_{p}$ are the modified Bessel functions Erdelyi  with
indices $p=l+\frac{1}{2}$. From definition (121) follows then

$f_{\alpha,lm}=(-k)^{-\alpha-1}\partial_{k}^{\alpha+1}f_{-1,lm}\,.$ (123)

Given the decomposition of $f_{\alpha}$ one can also determine the spherical
harmonics for $\gamma^{2}f_{\alpha}$ and $\gamma\sim 0$. For this we take
$\gamma^{2}\sim 2(1-\cos\gamma)$ and rotate our coordinate system such that
$\gamma=\theta$. From the recurrence relations of the Legendre
polynomialsErdelyi  in the spherical harmonics Erdelyi  follows then:

$(\cos\gamma f)_{lm}=\frac{(l+1)f_{l+1,m}-lf_{l-1,m}}{2l+1}\,.$ (124)

Now we need to determine the asymptotics of these expressions for large
$\kappa$ and $p$ with $p=l+\frac{1}{2}$.

###  D.2 Asymptotics

As the vacuum stress is given by a linear differential operator on the Green
function, its asymptotic behavior for $\bm{r}_{0}\rightarrow\bm{r}$ is
determined by this linear differential operator applied on functions of type
$IK$ and their derivatives. We take the choice $r>r_{0}$ (or $r<r_{0}$) in Eq.
(122). Using the Bessel identities Erdelyi

$\displaystyle\frac{2\alpha}{x}I_{\alpha}(x)$ $\displaystyle=$ $\displaystyle
I_{\alpha-1}(x)-I_{\alpha+1}(x)\,,$ $\displaystyle 2I_{\alpha}^{\prime}(x)$
$\displaystyle=$ $\displaystyle I_{\alpha-1}(x)+I_{\alpha+1}(x)\,,$
$\displaystyle-\frac{2\alpha}{x}K_{\alpha}(x)$ $\displaystyle=$ $\displaystyle
K_{\alpha-1}(x)-K_{\alpha+1}(x)\,,$ $\displaystyle-2K_{\alpha}^{\prime}(x)$
$\displaystyle=$ $\displaystyle K_{\alpha-1}(x)+K_{\alpha+1}(x)$ (125)

we can always reduce the stress to the form:

$\displaystyle\sigma$ $\displaystyle\sim$ $\displaystyle R_{1}I_{p}(\kappa\chi
r_{0})K_{p}(\kappa\chi r)$ (126) $\displaystyle+R_{2}I_{p+1}(\kappa\chi
r_{0})K_{p}(\kappa\chi r)$ $\displaystyle+R_{3}I_{p}(\kappa\chi
r_{0})K_{p+1}(\kappa\chi r)$ $\displaystyle+R_{4}I_{p+1}(\kappa\chi
r_{0})K_{p+1}(\kappa\chi r)$

where the $R_{i}$ are some rational functions of $r,r_{0},\kappa$ and $p$. We
only need the asymptotics of the products $IK$ for large index and large
argument.

The required uniform asymptotic behaviour of the Bessel functions is well
understood. For large $\alpha$ and finite $z$ we have Olver :

$\displaystyle I_{\alpha}(\alpha z)$ $\displaystyle\sim$
$\displaystyle\frac{e^{\alpha\xi}}{(2\pi\alpha)^{1/2}(1+z^{2})^{1/4}}\sum_{s=0}^{\infty}\frac{U_{s}(w)}{\alpha^{s}}\,,$
$\displaystyle K_{\alpha}(\alpha z)$ $\displaystyle\sim$
$\displaystyle\left(\frac{\pi}{2\alpha}\right)^{1/2}\frac{e^{-\alpha\xi}}{(1+z^{2})^{1/4}}\sum_{s=0}^{\infty}(-1)^{s}\frac{U_{s}(w)}{\alpha^{s}}\,,$
$\displaystyle\xi$ $\displaystyle=$
$\displaystyle(1+z^{2})^{1/2}+\ln\frac{z}{1+(1+z^{2})^{1/2}}\,,$
$\displaystyle w$ $\displaystyle=$ $\displaystyle(1+z^{2})^{-1/2}$ (127)

in terms of the polynomials $U_{s}$ defined by the relations

$\displaystyle U_{s+1}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}x^{2}(1-x^{2})U_{s}^{\prime}(x)+\frac{1}{8}\int_{0}^{x}(1-5t^{2})U_{s}(t)\,dt\,,$
$\displaystyle U_{0}(x)$ $\displaystyle=$ $\displaystyle 1\,.$ (128)

This gives the asymptotics of the $IK$ products:

$\displaystyle
I_{\alpha_{1}}(\alpha_{1}z_{1})K_{\alpha_{2}}(\alpha_{2}z_{2})\sim\frac{e^{\alpha_{1}\xi_{1}-\alpha_{2}\xi_{2}}}{2(\alpha_{1}\alpha_{2})^{1/2}\,(1+z_{1}^{2})^{1/4}\,(1+z_{2}^{2})^{1/4}}$
$\displaystyle\qquad\times\sum_{\begin{subarray}{c}0\leq s_{1}<\infty\\\ 0\leq
s_{2}\leq
s_{1}\end{subarray}}(-1)^{s_{2}}\frac{U_{s_{1}-s_{2}}(w_{1})}{\alpha_{1}^{s_{1}-s_{2}}}\frac{U_{s_{2}}(w_{2})}{\alpha_{2}^{s_{2}}}\,.$
(129)

We use this equation to get the asymptotic forms of all the needed $IK$
combinations, by taking $\alpha_{1},\alpha_{2}=p,p+1$ as necessary, and
$z_{1}=(\kappa\chi r)/p$ for the $p$ case and $z_{1}=(\kappa\chi r_{0})/(p+1)$
for the $(p+1)$ case. From here we expand around $\Lambda\to\infty$, take
$r\to r_{0}$, and obtain the divergency. The calculations are best done with
computer algebra; in the following we state the main result.

###  D.3 Result

We find that the $D_{0}$ contribution to the stress almost takes care of all
the infinities described in Eqs. (110) and (112) with only one term left:

$\left.\mathcal{W}^{\theta}_{\theta}-\mathcal{W}^{\theta}_{\theta}\right|_{D_{0}}=\frac{p^{3}}{w^{3}}\,\frac{rn^{\prime\prime}_{\infty}(r)-n^{\prime}_{\infty}(r)}{3\kappa^{2}r^{2}}\,.$
(130)

If we think of the $p$–summation as an integration and write
$p=w\cos\vartheta$ and $\kappa r=w\sin\vartheta$ such that we get a
two–dimensional integral in polar coordinates with (large) radius $w$, the
result (130) leads to a logarithmic divergence in $\Lambda$. Using Eq. (109)
we may also express the result as:

$\left.\mathcal{W}^{\theta}_{\theta}-\mathcal{W}^{\theta}_{\theta}\right|_{D_{0}}\sim\frac{p^{3}}{w^{3}}\,\frac{2(\beta_{E}+\beta_{M})}{r}$
(131)

for large $\kappa$ in terms of the scattering amplitudes $\beta_{E}$ and
$\beta_{M}$ defined in Eq. (30). These scattering amplitudes characterize the
$D_{1}$ waves of Eq. (29). In fact, we find that the left–over term
$\left.\mathcal{W}^{\theta}_{\theta}-\mathcal{W}^{\theta}_{\theta}\right|_{D_{0}}$
is entirely accounted for by the contribution of $D_{1}$ to the stress. So, if
we include $D_{1}$ in the renormalization procedure, the stress becomes
finite.

##  Appendix E Vector spherical harmonics

In this appendix we state the definitions, main relations and sum formulas for
the vector spherical harmonics we use in this paper.

###  E.1 Definition

According to Barrera et al. Barrera  the vector spherical harmonics are
defined as:

$\displaystyle\bm{Y}_{lm}(\theta,\phi)$ $\displaystyle=$
$\displaystyle\hat{\mathbf{r}}Y_{lm}(\theta,\phi)\,,$
$\displaystyle\bm{\Phi}_{lm}(\theta,\phi)$ $\displaystyle=$
$\displaystyle\bm{r}\times\nabla Y_{lm}(\theta,\phi)\,,$
$\displaystyle\bm{\Psi}_{lm}(\theta,\phi)$ $\displaystyle=$ $\displaystyle
r\nabla Y_{lm}(\theta,\phi)$ (132)

where the $Y_{lm}$ are the usual spherical harmonics Erdelyi  and
$\hat{\mathbf{r}}=\bm{r}/r$. The vector spherical harmonics establish three
orthogonal vectors in three–dimensional space for equal $lm$, and they are
orthogonal functions on the unit sphere for different $lm$.

###  E.2 Decomposition

Any vector function $\bm{V}(\theta,\phi)$ on the unit sphere can be decomposed
into vector spherical harmonics:

$\displaystyle\bm{V}(\theta,\phi)$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\Bigl{(}A_{lm}\bm{Y}_{lm}(\theta,\phi)+B_{lm}\bm{\Phi}_{lm}(\theta,\phi)$
(133) $\displaystyle\qquad\qquad+\,C_{lm}\bm{\Psi}_{lm}(\theta,\phi)\Bigr{)}$

with the coefficients given by Barrera

$\displaystyle A_{lm}$ $\displaystyle=$
$\displaystyle\int\bm{V}(\theta,\phi)\cdot\bm{Y^{*}}_{lm}(\theta,\phi)\,d\Omega\,,$
$\displaystyle B_{lm}$ $\displaystyle=$
$\displaystyle\frac{1}{l(l+1)}\int\bm{V}(\theta,\phi)\cdot\bm{\Phi^{*}}_{lm}(\theta,\phi)\,d\Omega\,,$
$\displaystyle C_{lm}$ $\displaystyle=$
$\displaystyle\frac{1}{l(l+1)}\int\bm{V}(\theta,\phi)\cdot\bm{\Psi^{*}}_{lm}(\theta,\phi)\,d\Omega\,.$
(134)

###  E.3 Sum identities

In order to derive our results for the vacuum stress we make use of several
sum identities of the vector spherical harmonics that are nigh impossible to
find in the literature. Therefore we list them here. They originate from the
familiar sum identity Erdelyi

$\frac{4\pi}{2l+1}\sum_{m=-l}^{l}Y_{lm}(\theta_{1},\phi_{1})\,Y^{*}_{lm}(\theta_{2},\phi_{2})=P_{l}(\cos\gamma)$
(135)

where $\gamma$ denotes the angle between the two points
$(\theta_{1},\phi_{1})$ and $(\theta_{1},\phi_{2})$ on the unit sphere, with

$\cos\gamma=\cos\theta_{1}\cos\theta_{2}+\sin\theta_{1}\sin\theta_{2}\cos(\phi_{1}-\phi_{2})\,.$
(136)

Let us explain our procedure by giving an example: the sum of
$\bm{\Psi}_{lm}(\theta,\phi)\otimes\bm{\Psi}_{lm}(\theta,\phi)^{*}$ over $m$.
From definition (132) follows

$\displaystyle\bm{\Psi}_{lm}(\theta_{1},\phi_{1})\otimes\bm{\Psi}_{lm}(\theta_{2},\phi_{2})^{*}$
$\displaystyle=\begin{bmatrix}0&0&0\\\
0&\partial_{\theta_{1}}\partial_{\theta_{2}}&\frac{1}{\sin\theta_{2}}\partial_{\theta_{1}}\partial_{\phi_{2}}\\\
0&\frac{1}{\sin\theta_{1}}\partial_{\phi_{1}}\partial_{\theta_{2}}&\frac{1}{\sin\theta_{1}\sin\theta_{2}}\partial_{\phi_{1}}\partial_{\phi_{2}}\\\
\end{bmatrix}Y_{lm1}Y^{*}_{lm2}$ (137)

with $Y_{lmi}\equiv Y_{lm}(\theta_{i},\phi_{i})$. Hence we obtain

$\displaystyle\sum_{m=-l}^{l}\bm{\Psi}_{lm}(\theta_{1},\phi_{1})\otimes\bm{\Psi}_{lm}(\theta_{2},\phi_{2})^{*}$
$\displaystyle=\begin{bmatrix}0&0&0\\\
0&\partial_{\theta_{1}}\partial_{\theta_{2}}&\frac{1}{\sin\theta_{2}}\partial_{\theta_{1}}\partial_{\phi_{2}}\\\
0&\frac{1}{\sin\theta_{1}}\partial_{\phi_{1}}\partial_{\theta_{2}}&\frac{1}{\sin\theta_{1}\sin\theta_{2}}\partial_{\phi_{1}}\partial_{\phi_{2}}\\\
\end{bmatrix}\frac{2l+1}{4\pi}\,P_{l}(\cos\gamma)$ (138)

from Eq. (135). We take the limit $\gamma\rightarrow 0$ and use for the $k$–th
derivatives $P^{(k)}_{l}(1)$ the formula

$P^{(k)}_{l}(1)=\frac{\prod_{q=-k+1}^{k}(l+q)}{\prod_{q=1}^{k}2q}$ (139)

that is easy to derive by induction. The result is given below, and also the
results for similar sum identities we have derived with the same method:

$\displaystyle\sum_{m=-l}^{l}\bm{Y}_{lm}\otimes\bm{Y}_{lm}^{*}$
$\displaystyle=\begin{bmatrix}1&0&0\\\ 0&0&0\\\ 0&0&0\\\
\end{bmatrix}\frac{2l+1}{4\pi}$
$\displaystyle\sum_{m=-l}^{l}\bm{Y}_{lm}\otimes\bm{\Psi}_{lm}^{*}$
$\displaystyle=0\,,$
$\displaystyle\sum_{m=-l}^{l}\bm{Y}_{lm}\otimes\bm{\Phi}_{lm}^{*}$
$\displaystyle=0\,,$
$\displaystyle\sum_{m=-l}^{l}\bm{\Psi}_{lm}\otimes\bm{Y}_{lm}^{*}$
$\displaystyle=0\,,$
$\displaystyle\sum_{m=-l}^{l}\bm{\Psi}_{lm}\otimes\bm{\Psi}_{lm}^{*}$
$\displaystyle=\begin{bmatrix}0&0&0\\\ 0&1&0\\\ 0&0&1\\\
\end{bmatrix}\frac{l(l+1)(2l+1)}{8\pi}\,,$
$\displaystyle\sum_{m=-l}^{l}\bm{\Psi}_{lm}\otimes\bm{\Phi}_{lm}^{*}$
$\displaystyle=\begin{bmatrix}0&0&0\\\ 0&0&1\\\ 0&-1&0\\\
\end{bmatrix}\frac{l(l+1)(2l+1)}{8\pi}\,,$
$\displaystyle\sum_{m=-l}^{l}\bm{\Phi}_{lm}\otimes\bm{Y}_{lm}^{*}$
$\displaystyle=0\,,$
$\displaystyle\sum_{m=-l}^{l}\bm{\Phi}_{lm}\otimes\bm{\Psi}_{lm}^{*}$
$\displaystyle=\begin{bmatrix}0&0&0\\\ 0&0&-1\\\ 0&1&0\\\
\end{bmatrix}\frac{l(l+1)(2l+1)}{8\pi}\,,$
$\displaystyle\sum_{m=-l}^{l}\bm{\Phi}_{lm}\otimes\bm{\Phi}_{lm}^{*}$
$\displaystyle=\begin{bmatrix}0&0&0\\\ 0&1&0\\\ 0&0&1\\\
\end{bmatrix}\frac{l(l+1)(2l+1)}{8\pi}$ (140)

for $\theta_{0}=\theta$ and $\phi_{0}=\phi$.

## References

  * (1) I. E. Dzyaloshinskii, E. M. Lifshitz, and L. P. Pitaevskii, The general theory of van der Waals forces, Adv. Phys. 10, 165 (1961).
  * (2) S. Y. Buhmann, Dispersion Forces (Springer, Berlin, 2012),
  * (3) E. Shahmoon, Van der Waals and Casimir–Polder dispersion forces, in Ref. Forces  below.
  * (4) W. M. R. Simpson and U. Leonhardt (eds.) Forces of the quantum vacuum (World Scientific, Singapore, 2015).
  * (5) H. B. G. Casimir, On the attraction between two perfectly conducting plates, Koninkl. Ned. Akad. Wetenschap. 51, 793 (1948).
  * (6) A. W. Rodriguez, F. Capasso, and S. G. Johnson, The Casimir effect in microstructured geometries, Nat. Photon. 5, 211 (2011).
  * (7) H. B. G. Casimir and D. Polder, The Influence of Retardation on the London–van der Waals Forces, Phys. Rev. 73, 360 (1948).
  * (8) L. P. Pitaevskii, Comment on “Casimir force acting on magnetodielectric bodies embedded in media”, Phys. Rev. A 73, 047801 (2006).
  * (9) F. A. Burger, J. Fiedler, and S. Y. Buhmann, Zero-point electromagnetic stress tensor for studying Casimir forces on colloidal particles in media, Europhys. Lett. 121, 24004 (2018).
  * (10) J. N. Munday, F. Capasso, and V. A. Parsegian, Measured long–range repulsive Casimir–Lifshitz forces, Nature 457, 170 (2009).
  * (11) R. Zhao, L. Li, S. Yang, W. Bao, Y. Xia, P. Ashby, Y. Wang, and X. Zhang, Stable Casimir equilibria and quantum trapping, Science 364, 984 (2019).
  * (12) R. S. Decca, Measuring Casimir Phenomena, in Ref. Forces .
  * (13) W. M. R. Simpson, Surprises in Theoretical Casimir Physics (Springer, Berlin, 2014).
  * (14) W. M. R. Simpson, S. A. R. Horsley, and U. Leonhardt, Divergence of Casimir stress in inhomogeneous media, Phys. Rev. A 87, 043806 (2013); ibid. 88, 059901 (2013).
  * (15) I. Griniasty and U. Leonhardt, Casimir stress inside planar materials, Phys. Rev. A 96, 032123 (2017).
  * (16) I. Griniasty and U. Leonhardt, Casimir stress in materials: Hard divergency at soft walls, Phys. Rev. B 96, 205418 (2017).
  * (17) K. A. Milton, Casimir Energy for a Spherical Cavity in a Dielectric: Toward a Model for Sonoluminescence? Ann. Phys. (New York) 127, 49 (1980).
  * (18) I. Brevik and H. Kolbenstvedt, The Casimir effect in a solid ball when $\varepsilon\mu=1$, Ann. Phys. (New York) 143, 179 (1982).
  * (19) K. A. Milton and Y. J. Ng, Casimir energy for a spherical cavity in a dielectric: Applications to sonoluminescence, Phys. Rev. E 55, 4207 (1997).
  * (20) I. Brevik, V. V. Nesterenko and I. G. Pirozhenko, Direct mode summation for the Casimir energy of a solid ball, J. Phys. A 31, 8661 (1998).
  * (21) I. Brevik, V. N. Marachevsky, and K. A. Milton, Identity of the van der Waals Force and the Casimir Effect and the Irrelevance of These Phenomena to Sonoluminescence, Phys. Rev. Lett. 82, 3948 (1999).
  * (22) G. Barton, Perturbative check on the Casimir energies of nondispersive dielectric spheres, J. Phys. A 32, 525 (1999).
  * (23) M. Bordag, K. Kirsten and D. Vassilevich, Ground state energy for a penetrable sphere and for a dielectric ball, Phys. Rev. D 59, 85011 (1999).
  * (24) Y. Avni and U. Leonhardt, Casimir self–stress in a dielectric sphere, Ann. Phys. (New York) 395, 326 (2018).
  * (25) K. A. Milton, P. Parashar, I. Brevik, G. Kennedy, Self–stress on a dielectric ball and Casimir–Polder forces, Ann. Phys. (New York) 412, 168008 (2020).
  * (26) L. D. Landau and E. M. Lifshitz, Statistical Physics, Part 2 (Pergamon, Oxford, 1980).
  * (27) E. M. Lifshitz, The Theory of Molecular Attractive Forces between Solids, J. Exper. Theoret. Phys. USSR 29, 94 (1954) [English translation: Sov. Phys. JETP 2, 73 (1956)].
  * (28) S. Scheel, The Casimir stress in real materials, in Ref. Forces .
  * (29) R. M. Wald, Trace anomaly of a conformally invariant quantum field in curved spacetime, Phys. Rev. D 17, 1477 (1978).
  * (30) S. Weinberg, The cosmological constant problem, Rev. Mod. Phys. 61, 1 (1989).
  * (31) S. W. Hawking, T. Hertog, and H. S. Reall, Trace anomaly driven inflation, Phys. Rev. D 63, 083504 (2001).
  * (32) L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields (Butterworth-Heinemann, Amsterdam, 2003).
  * (33) L. Amendola and S. Tsujikawa, Dark Energy: Theory and Observations (Cambridge University Press, Cambridge, 2010).
  * (34) G. Volovik, The Universe in a Helium Droplet (Oxford University Press, Oxford, 2003).
  * (35) C. Barcelo, S. Liberati and M. Visser, Analogue Gravity, Living Rev. Relativity 8, 12 (2005).
  * (36) P. O. Fedichev and U. R. Fischer, Gibbons–Hawking Effect in the Sonic de Sitter Space–Time of an Expanding Bose–Einstein–Condensed Gas, Phys. Rev. Lett. 91, 240407 (2003).
  * (37) E. B. Kolomeisky, Natural analog to cosmology in basic condensed matter physics, Phys. Rev. B 100, 140301(R) (2019).
  * (38) E. B. Kolomeisky, Normal modes of vibrations around Hubble flow in jellium, Phys. Rev. B 101, 174304 (2020).
  * (39) E. B. Kolomeisky, Analog de Sitter space in a controlled Coulomb explosion, Phys. Rev. A 103, 043101 (2021).
  * (40) J. Steinhauer, M. Abuzarli, T. Aladjidi, T. Bienaimé, C. Piekarski, W. Liu, E. Giacobino, A. Bramati, and Q. Glorieux, Analogue cosmological particle creation in an ultracold quantum fluid of light, arXiv:2102.08279.
  * (41) U. Leonhardt, The case for a Casimir cosmology, Phil. Trans. R. Soc. A 378, 20190229 (2020).
  * (42) C. Pethick and H. Smith, Bose–Einstein condensation in dilute gases (Cambridge University Press, Cambridge, 2008).
  * (43) L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media (Pergamon, Oxford, 1984).
  * (44) Planck Collaboration, Planck 2018 results. VI. Cosmological parameters, Astron. & Astrophys. 641, A6 (2020).
  * (45) W. Gordon, Zur Lichtfortpflanzung nach der Relativitätstheorie, Ann. Phys. (Leipzig) 72, 421 (1923).
  * (46) J. Plebanski, Electromagnetic Waves in Gravitational Fields, Phys. Rev. 118, 1396 (1960).
  * (47) W. Schleich and M. O. Scully, General relativity and modern optics, in New trends in atomic physics: Les Houches, session XXXVIII, 1982 by G. Grynberg and R. Stora (eds.) (Elsevier, Amsterdam, 1984).
  * (48) U. Leonhardt and T. G. Philbin, Geometry and Light: the Science of Invisibility, (Dover, Mineola, 2010).
  * (49) U. Leonhardt, Lifshitz theory of the cosmological constant, Ann. Phys. (New York) 411, 167973 (2019).
  * (50) D. Berechya and U. Leonhardt, Lifshitz cosmology: Hubble tension and quantum vacuum, MNRAS 507, 3473 (2021).
  * (51) L. Verde, T. Treu, and A. G. Riess, Tensions between the early and late Universe, Nat. Astron. 3, 891 (2019).
  * (52) The famous, solved problem of the Casimir force on a perfectly conducting shell Boyer ; MDS ; Balian ; Nesterenko  can be reduced to the problem of the Casimir force on the homogeneous sphere in a uniform background Avni : the force on the shell is the force on a solid, perfectly conducting sphere plus the force on a hollow sphere in a uniform, perfectly conducting background.
  * (53) T. H. Boyer, Quantum electromagnetic zero-point energy of a conducting spherical shell and the Casimir model for a charged particle, Phys. Rev. 174, 1764 (1968).
  * (54) K. A. Milton, L. L. De Raad Jr, J. Schwinger, Casimir self-stress on a perfectly conducting spherical shell, Ann. Phys. (New York) 115, 388 (1978).
  * (55) R. Balian and B. Duplantier, Electromagnetic waves near perfect conductors. I. Multiple scattering expansions. Distribution of modes, Ann. Phys. (New York) 104, 300 (1977).
  * (56) V. V. Nesterenko, Simple method for calculating the Casimir energy for a sphere, Phys. Rev. D 57, 1284 (1998).
  * (57) U. Leonhardt, Essential Quantum Optics: From Quantum Measurements to Black Holes, (Cambridge University Press, Cambridge, 2010).
  * (58) Relation (13) follows from Eq. (8) in spherical coordinates with $ds^{2}=dr^{2}+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta\,d\phi^{2}$.
  * (59) J. R. Gott III, M. Jurić, D. Schlegel, F. Hoyle, M. Vogeley, M. Tegmark, N. Bahcall, and J. Brinkmann, A Map of the Universe, Astrophys. J. 624, 463 (2005).
  * (60) U. Leonhardt and T. G. Philbin, Perfect imaging with positive refraction in three dimensions, Phys. Rev. A 81, 011804 (2010).
  * (61) U. Leonhardt, Cosmological horizons radiate, Europhys. Lett. 135, 10002 (2021).
  * (62) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, In the Realm of the Hubble tension — a Review of Solutions, Class. Quantum Grav. 38, 153001 (2021).
  * (63) J. E. Bjorkholm, R. H. Freeman, A. Ashkin, and D. B. Pearson, Observation of Focusing of Neutral Atoms by the Dipole Forces of Resonance-Radiation Pressure, Phys. Rev. Lett. 41, 1361 (1978).
  * (64) R. Ozeri, N. Katz, J. Steinhauer, and N. Davidson, Colloquium: Bulk Bogoliubov excitations in a Bose-Einstein condensate, Rev. Mod. Phys. 77, 187 (2005).
  * (65) A. Erdélyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi, Higher Transcendental Functions (McGraw-Hill, New York, 1981).
  * (66) C.-S. Brés, With a fine-tooth comb, Nat. Phys. 16, 600 (2020).
  * (67) R. G. Barrera, G. A. Estevez, and J. Giraldo, Vector spherical harmonics and their application to magnetostatics, Eur. J. Phys. 6, 287 (1985).
  * (68) One should use the homogeneous wave equation here, i.e. Eq. (66) with right–hand side equal to zero, because $\bm{r}_{0}\neq\bm{r}$ (before the limit $\bm{r}_{0}\rightarrow\bm{r}$ is taken).
  * (69) R. E. Langer, On the connection formulas and the solutions of the wave equation, Phys. Rev. 51, 669 (1937).
  * (70) F. Olver, Asymptotics and special functions, (CRC Press, London, 1997).