# Van der Waals chain: a simple model for Casimir forces in dielectrics

Helmut Hörner1, Lukas M. Rachbauer1, Stefan Rotter1, and Ulf Leonhardt2
1Institute for Theoretical Physics, Vienna University of Technology (TU Wien),
A–1040 Vienna, Austria  2Department of Physics of Complex Systems, Weizmann
Institute of Science, Rehovot 761001, Israel

###### Abstract

The Casimir force between dielectric bodies is well–understood, but not the
Casimir force inside a dielectric, in particular its renormalization. We
develop and analyse a simple model for the Casimir forces inside a medium that
is completely free of renormalization and show then how renormalization
emerges. We consider a one–dimensional chain of point particles interacting
with each other by scattering the zero–point fluctuations of the
electromagnetic field confined to one dimension. We develop a fast, efficient
algorithm for calculating the forces on each particle and apply it to study
the macroscopic limit of infinitely many, infinitely weak scatterers. The
force density converges for piece–wise homogeneous media, but diverges in
inhomogeneous media, which would cause instant collapse in theory. We argue
that short–range counter forces in the medium prevent this collapse in
reality. Their effect appears as the renormalization of the Casimir stress in
dielectrics. Our simple model also allows us to derive an elementary analogue
of the trace anomaly of quantum fields in curved space.

##  I Introduction

The Casimir effect Milonni ; Milton ; BKMM ; Rodriguez ; Buhmann ; Forces ;
State  surely is one of the most fascinating and “one of the least intuitive
consequences of quantum electrodynamics” Schwinger . Yet it all began with a
mundane, practical question LamoreauxPhysToday : under which conditions are
colloidal suspensions Israelachvili  — such as liquid paint — stable? A
colloidal suspension is a fluid containing particles that repel each other on
short distances by steric forces and attract each other on longer distances by
van der Waals forces Hermann . The balance of forces determines whether the
suspension is stable Israelachvili , for example, whether liquid paint stays
liquid.

Casimir was confronted with the following problem. Measurements indicated that
the interparticle forces fall off like $r^{-7}$ with distance $r$, but
London’s theory London  of static van der Waals forces predicted a $r^{-6}$
dependance. Overbeck suggested CasimirPolder  that retardation weakens the van
der Waals force. Following this idea, Casimir and Polder CasimirPolder
deduced the $r^{-7}$ asymptotics from a full quantum–electrodynamical theory
ShahmoonvdW . Their theory CasimirPolder  was complicated and cumbersome, but
their result was surprisingly simple and only depending on natural constants.
Upon showing it to Bohr, he “mumbled something about zero–point energy” Letter
, which put Casimir on the track to discover the Casimir effect: Casimir  two
planar perfect conductors in empty space attract each other by modifying the
zero–point energy of the electromagnetic field. That energy is infinite, but
the modifications are finite.

Lifshitz Lifshitz  generalized the theory from perfect conductors to planar
dielectrics separated by vacuum. Dzyaloshinskii and Pitaevskii joined in,
replaced the vacuum by a homogeneous fluid in thermal equilibrium, and still
managed to extend the theory to this case DLP . Schwinger et al. confirmed
their theory using a different starting point Schwinger . The fluid gives rise
to rich physical phenomena observed relatively recently: repulsive Casimir
forces Milling ; Meurk ; Lee ; Levitation  and stable Casimir equilibria
CasimirEquilibrium . Two planar dielectrics separated by a homogeneous fluid
constitute a piece-wise homogeneous planar medium. Dzyaloshinskii and
Pitaevskii thought DP  their theory applies to inhomogeneous dielectrics as
well, but they were wrong SimpsonThesis .

The problem is the renormalization required to extract the finite, physically
valid force density from the infinite zero–point energy. It took modern
computational tools to discover Simpson  that the renormalization procedure DP
of Dzyaloshinskii and Pitaevskii fails in general. Some partial solutions of
the renormalization problem are known Renormalization ; Itay2 ; Parashar ; Li
; Itai ; Shayit ; Vinas  but lack a sound physical justification.

Here we address the problem of renormalization in inhomogeneous media with a
theory inherently free of renormalization. We take as our starting point
Casimir’s initial problem of neutral particles in a fluid and simplify it
further, inspired by the scattering approach to Casimir forces Lambrecht ;
Wirzba ; Ingold . We consider a one–dimensional chain of point scatterers
interacting with each other by virtual electromagnetic waves confined to one
dimension Fiedler  (by, e.g., a waveguide Shahmoon ). Varying the scattering
strengths of the particles we implement inhomogeneous dielectrics. The
scatterers are immersed in a fluid of uniform refractive index, but by
rescaling spatial coordinates we can replace the fluid by vacuum and calculate
the van der Waals forces in the new coordinates as vacuum forces.

The force on each scatterer is always finite. We find, however, that in the
macroscopic limit of infinitely many, infinitely weak scatterers, the force
density tends to infinity if they form an inhomogeneous medium. If, on the
other hand, the scatterers are arranged to be piece–wise homogeneous, the
force density converges to finite values, except at interfaces. Perturbations
of homogeneity, however small they might be, would produce infinite forces in
the macroscopic limit. Liquid paint would instantly solidify, unless it is
absolutely homogenous, which does not happen in reality. Something is missing
in the theory. What is it?

We follow and refine an argument by Pitaevskii Pitaevskii  that local
restoring forces of the background fluid counteract the van der Waals forces
on the scatterers. The net result is finite in the macroscopic limit. These
counter forces are naturally occurring in colloidal suspensions. It is
well–known Israelachvili  that they may stabilize homogeneous suspensions,
here we argue that they also prevent a catastrophic collapse in inhomogeneous
fluids. We identify how the counter forces are expressed in terms of
macroscopic dielectric quantities. Mathematically, they appear as
renormalizers, in agreement with previous renormalization procedures in
inhomogeneous dielectrics Renormalization ; Itay2 ; Parashar ; Li ; Itai ;
Shayit ; Vinas  and in agreement with Lifshitz theory DLP , but physically
they are quite something else. Instead of seeing renormalization as a
mathematical procedure to remove infinities, we regard it as the result of a
physical process, the local equilibrium of forces.

Inhomogeneities in fluids play a role in surface forces, and a better theory
of van der Waals forces in inhomogeneous liquids may solve some of the
mysteries there Sivan . Realizing the limitations of Lifshitz theory may also
help explaining the disagreement between theory and experiment in the Casimir
force at long distances LifshitzDrude . Finally, the renormalization of
Casimir forces in inhomogeneous media is essential for understanding the
quantum vacuum in curved space BD . This is because curved space–time LL2
acts just like an inhomogeneous dielectric medium on electromagnetic fields
Gordon ; Plebanski ; Schleich ; LeoPhil . We show how in our model the trace
anomaly Itai ; Wald  arises that plays the role of the cosmological constant
Annals ; Tkatchenko . We thus have reasons to believe that our model, despite
its simplicity, elucidates the essentials of some diverse physics, ranging
from liquid paint to the cosmos.

##  II The model

Figure 1:  The model: one–dimensional chain of $N$ point particles (here
$N=5$) with varying polarizabilities $\alpha_{i}$ (varying shades of gray)
interacting with each other by scattering electromagnetic waves in the vacuum
state, freely propagating between them. The propagation is described by P
matrices [Eq. (13)] and the scattering by R matrices [Eq. (19)].

Consider a one–dimensional chain of $N$ point scatterers at positions $x_{i}$
with constant distance $\delta$ from each other (Fig. 1). The electric
polarizability $\alpha_{i}$ of each scatterer may vary such that they
constitute a medium in the macroscopic limit $\delta\rightarrow 0$ with

$\chi(x)=\frac{\alpha_{i}}{\delta}$ (1)

and $\chi$ kept finite. The polarizabilities may depend on frequency
(establishing a dispersive medium). We will show in Sec. III that $\chi$
represents the electric susceptibility Fiedler  and

$n=\sqrt{1+\chi}$ (2)

the refractive index. One could also assume a varying distance $\delta$ and
constant polarizability $\alpha$ or both $\delta$ and $\alpha$ varying, but
here we take $\delta$ as constant for simplicity. The scatterers shall
interact with each other by the vacuum forces Rodriguez ; Buhmann ; Forces ;
State  of the electromagnetic field confined to one–dimensional propagation.
These are the van der Waals forces with both retardation and multiple
scattering taken into account. In the following we work out what those forces
are (Fig. 2).

Figure 2:  van der Waals forces. a: homogeneous block of $N=51$ particles
(with $\chi$ illustrated by the gray filled curve) placed between $x_{1}=0$
and $x_{21}=1$. The dots show the numerically calculated van der Waals forces
in units of $\hbar c$ and the solid black curve the analytic formula (81) in
the macroscopic limit. They agree well for $N\geq 20$. The force peaks at the
outer particles (not shown) and is attractive and antisymmetric; the net force
is zero. b: $N=101$ particles with inhomogeneous profile [Eq. (86) with the
parameters $\chi_{0}=1$, $a=0.15$]. The force density diverges in the
macroscopic limit ($N\rightarrow\infty$).

###  II.1 Vacuum forces

According to Hamilton’s equations, the force on each scatterer is given by the
negative derivative of the Hamiltonian with respect to the position $x_{i}$.
For the Hamiltonian we take the zero–point energy of the electromagnetic field
expressed as the integral of all vacuum energies $E/2$ times the density of
states $D(E)$. Here $E$ represents the energy $E=\hbar\omega$ corresponding to
the frequency $\omega$. The density of states is given by the
Friedel–Krein–Lloyd formula Friedel ; Krein ; Lloyd ; Faulkner  in terms of
the scattering phase $\varphi$ as $D(E)=(2\pi)^{-1}\partial_{E}\varphi$. We
integrate by parts and obtain for the van der Waals force

$F_{i}=\frac{1}{4\pi}\int_{0}^{\infty}\partial_{i}\varphi\,dE\,,\quad
E=\hbar\omega\,.$ (3)

The scattering phase is the sum of all the phases of the eigenvalues of the
classical scattering matrix $\rm{S}$. Since $\det\rm{S}$ is the product of all
those eigenvalues we may write:

$\varphi=-i\ln\,\det\rm{S}\,.$ (4)

With Eqs. (3) and (4) we have expressed the quantum force on each scatterer in
terms of classical quantities. Quantum mechanics enters though $E=\hbar\omega$
and the fact that the zero–point energy is $E/2$.

To proceed, we note that the scattering matrix is unitary, and that all
unitary matrices can be expressed via the Cayley transform as

${\rm S}=\frac{\mathbb{1}-i{\rm K}}{\mathbb{1}+i{\rm K}}\quad\mbox{ with
Hermitian $\rm{K}$}.$ (5)

We extend ${\rm S}(E)$ to negative energies corresponding to negative
frequencies. As the classical electromagnetic field is real, the
negative–frequency component is the complex conjugate of the positive
frequency component and runs backwards in time. Consequently, ${\rm S}(-E)$
must be the inverse of ${\rm S}(E)$. From this follows that

${\rm K}(-E)=-{\rm K}(E)\,.$ (6)

We thus obtain

$F_{i}=\frac{i}{4\pi}\left(\int_{0}^{+\infty}-\int_{-\infty}^{0}\right)\frac{\partial_{i}\det(\mathbb{1}+i\rm{K})}{\det(\mathbb{1}+i\rm{K})}\,dE$
(7)

and by deforming the integration contours on the complex plane (Fig. 3) to the
positive imaginary axis, and expressing them as an integral over imaginary
wave numbers $\kappa$ with

$\omega=ic\kappa$ (8)

and $c$ being the speed of light in vacuum, we get a rapidly converging
formula for the force:

$F_{i}=-\frac{\hbar
c}{2\pi}\int_{0}^{\infty}\frac{\partial_{i}\det(\mathbb{1}+i\rm{K})}{\det(\mathbb{1}+i\rm{K})}\,d\kappa\,.$
(9)

We are going to show (Sec. IIC) that $\det(\mathbb{1}+i\rm{K})$ is real on the
imaginary axis. Note that we have not performed any renormalization to
calculate the force of the quantum vacuum on the point scatterers.

Figure 3:  Contour deformation. We represent the integrals
$\int_{0}^{+\infty}$ and $-\int_{-\infty}^{0}$ over $\omega$ [Eq. (7)] by
contours (shown as arrows) running to the right or the left, respectively, and
move both contours to the positive imaginary axis ($\omega=ic\kappa$) where
they coincide, which gives Eq. (9).

###  II.2 Transfer matrix

The only inconvenient feature of formula (9) is the appearance of the
scattering matrix (via $\rm{K}$). However, for one–dimensional propagation we
may express $\det(\mathbb{1}+i\rm{K})$ in terms of the transfer matrix
$\rm{T}$. The transfer matrix is much easier to calculate than the scattering
matrix, because the total transfer matrix of the chain of scatterers is simply
the product of the individual transfer matrices Transfer .

Consider real frequencies first. A right–moving wave $e^{ikx}$ on the left of
the structure is transferred to the superposition
$\alpha\,e^{ikx}+\beta\,e^{-ikx}$ with complex $\alpha$ and $\beta$ on the
right (not to be confused with the polarizabilities $\alpha_{i}$) while a
left–moving wave $e^{-ikx}$ on the left gets transferred to
$\beta^{*}\,e^{ikx}+\alpha^{*}\,e^{-ikx}$ on the right. The transfer matrix
has thus the structure Transfer

${\rm T}=\begin{pmatrix}\alpha&\beta^{*}\\\ \beta&\alpha^{*}\end{pmatrix}$
(10)

with $\det\rm{T}=|\alpha|^{2}-|\beta|^{2}>0$ as the outgoing waves must not
exceed the ingoing waves. The scattering matrix ${\rm S}$ describes how
ingoing waves are turned into outgoing waves. The ingoing waves enter as
$e^{ikx}$ on the left and $e^{-ikx}$ on the right, while the outgoing waves
emerge as $e^{-ikx}$ on the left and $e^{ikx}$ on the right. One verifies that
the coefficients of the outgoing and ingoing waves are related by the matrix

${\rm
S}=\frac{1}{\alpha^{*}}\begin{pmatrix}\sqrt{|\alpha|^{2}-|\beta|^{2}}&\beta^{*}\\\
-\beta&\sqrt{|\alpha|^{2}-|\beta|^{2}}\end{pmatrix}.$ (11)

We see from Eq. (11) that $\det S=\alpha/\alpha^{*}$ and infer from Eq. (10)
and $\det S=\det(\mathbb{1}+i{\rm K})^{*}/\det(\mathbb{1}+i{\rm K})$ the
relation

$\det(\mathbb{1}+i{\rm K})=T_{22}\,.$ (12)

This relation, once established on the real axis of $E$, is extended by
complex continuation to the imaginary axis (Fig. 3) where we need it.

Having made the connection between the transfer matrix and the force, we use
the product rule of transfer matrices to express $\rm{T}$ in terms of the
individual transfer matrices. Between the scatterers at positions $x_{i}$ the
field propagates freely (Fig. 1). The transfer matrices of free propagation
are given by

${\rm P}_{i}=\begin{pmatrix}e^{-\kappa\delta_{i}}&0\\\
0&e^{\kappa\delta_{i}}\end{pmatrix}\quad\mbox{with}\quad\delta_{i}=x_{i+1}-x_{i}$
(13)

for purely imaginary wavenumbers where waves decay exponentially in
propagation direction. For the sake of calculating the force, we keep the
distance $\delta_{i}$ variable, as we need to differentiate with respect to
the coordinates. The beginning and the end of the propagation region outside
of the scatterers we put at the fixed points $x_{0}$ and $x_{N+1}$. The
precise values of these points will not influence the force. It is wise to
break the total transfer matrix into two parts ${\rm A}_{j}$ and ${\rm B}_{j}$
describing the transfer from $x_{0}$ to $x_{j-1}$ and from $x_{j}$ to
$x_{N+1}$, respectively. We denote the transfer matrices of the scatterers by
${\rm R}_{i}$ and obtain for the partial transfer matrices the recurrence
relations

${\rm A}_{j+1}={\rm P}_{j}\,{\rm R}_{j}\,{\rm A}_{j}\quad\mbox{and}\quad{\rm
B}_{j-1}={\rm B}_{j}\,{\rm R}_{j}\,{\rm P}_{j-1}$ (14)

with the initial values

${\rm A}_{1}={\rm P}_{0}\quad\mbox{and}\quad{\rm B}_{N}={\rm P}_{N}\,.$ (15)

With the help of the A and B matrices we get for the total transfer matrix:

${\rm T}={\rm B}_{j}\,{\rm R}_{j}\,{\rm A}_{j}\,.$ (16)

As $\partial_{j}=\partial_{\delta_{j-1}}-\partial_{\delta_{j}}$ we obtain from
Eq. (13) for the derivatives of the transfer matrix

$\partial_{j}{\rm T}=\kappa{\rm B}_{j}\left[\sigma_{z},{\rm R}_{j}\right]{\rm
A}_{j}$ (17)

written in terms of the commutator with the Pauli matrix $\sigma_{z}={\rm
diag}(1,-1)$.

We have arrived at expressions for the force that are convenient for numerical
calculations (Appendix A) and, as we shall see in Sec. III, for deducing the
force density in the macroscopic limit.

###  II.3 Point scatterers

It remains to specify the transfer matrices ${\rm R}_{i}$ of the point
scatterers (Fig. 1). We describe a point scatterer at position $x_{i}$ as a
delta–function susceptibility in the Helmholtz equation of the electromagnetic
field. We only consider one polarization component of the field and denote it
by $\psi$ (for the other polarization we would get the same). For purely
imaginary wavenumbers we thus have

$\partial_{x}^{2}\psi=\kappa^{2}\left[1+\alpha_{i}\delta(x-x_{i})\right]\psi\,.$
(18)

The field to the left of $x_{i}$ we write as the superposition
$a_{-}e^{-\kappa x}+a_{+}e^{\kappa x}$ while the field on the right is given
by the (different) superposition $b_{-}e^{-\kappa x}+b_{+}e^{\kappa x}$. To
satisfy Eq. (18) the field $\psi$ must be continuous at $x_{i}$ while the
derivative $\psi^{\prime}$ is discontinuous and changes by
$\kappa^{2}\alpha_{i}\psi(x_{i})$. From this we obtain the connection between
the coefficients, i.e. the transfer matrix:

${\rm R}_{i}=\mathbb{1}-\frac{\kappa\alpha_{i}}{2}(\sigma_{z}+i\sigma_{y})$
(19)

expressed here in terms of the Pauli matrices

$\sigma_{x}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\quad\sigma_{y}=\begin{pmatrix}0&-i\\\
i&0\end{pmatrix},\quad\sigma_{z}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}$ (20)

for mathematical convenience. Note that the ${\rm R}$ matrices are real, and
so are the ${\rm P}$ matrices (13) describing the free propagation between the
scatterers. Therefore the total transfer matrix and its derivatives are real,
too, and so is $\det(\mathbb{1}+i\rm{K})$ as we have stated in Sec. IIA. From
the commutation relations of the Pauli matrices (20) we get the relation
$[\sigma_{z},{\rm R}_{i}]=-\kappa\alpha_{i}\sigma_{x}$ and hence from Eq.
(17):

$\partial_{j}{\rm T}=-\kappa^{2}\alpha_{j}{\rm B}_{j}\sigma_{x}{\rm A}_{j}\,.$
(21)

This completes the calculation of the van der Waals force on a chain of point
particles using scattering theory and the transfer matrix.

##  III Macroscopic limit

In the macroscopic limit we assume that the inter–particle distance $\delta$
goes to zero with finite ratio (1) between polarizability $\alpha_{i}$ and
$\delta$. We translate the recurrence relations (14) into differential
equations and express the solution for the force density in terms of the
macroscopic Green function.

###  III.1 Differential equations

Take the partial transfer matrices A and B as continuous functions of $x$
evaluated at the points $x_{j}$ of the scatterers. Substituting
$\alpha_{i}=\chi\delta$ in Eq. (19) and using the recurrence relations (14)
for the A and B matrices we linearize the differences ${\rm A}_{j+1}-{\rm
A}_{j}$ and ${\rm B}_{j}-{\rm B}_{j-1}$ taking also into account the
$\delta$–dependence in the propagation matrices (13). We obtain the
differential equations

$\displaystyle\partial_{x}{\rm A}=-\kappa{\rm V}{\rm
A}\quad\mbox{and}\quad\partial_{x}{\rm B}=\kappa{\rm B}{\rm
V}\quad\mbox{with}$ $\displaystyle{\rm
V}=\sigma_{z}+\frac{\chi}{2}(\sigma_{z}+i\sigma_{y})$ (22)

in terms of the Pauli matrices (20). We obtain from the initial values (15) in
the macroscopic limit (where ${\rm P}\rightarrow\mathbb{1}$) the initial
conditions

${\rm A}(x_{0})={\rm B}(x_{N+1})=\mathbb{1}\,.$ (23)

According to Eqs. (9) and (12) the force is given in terms of the matrix
element $T_{22}$ and its spatial derivative $T_{22}^{\prime}$. Since ${\rm
R}\rightarrow\mathbb{1}$ for $\delta\rightarrow 0$ we obtain from Eq. (16)
${\rm T}={\rm BA}$, and so

$T_{22}=B_{21}A_{12}+B_{22}A_{22}\,.$ (24)

From Eq. (21) we get

$T_{22}^{\prime}=-\kappa^{2}\chi\delta\left(B_{21}A_{22}+B_{22}A_{12}\right).$
(25)

One verifies that for the relevant matrix elements the solutions of the
differential equations (22) have the structure

$\displaystyle A_{12}=(1-\kappa^{-1}\partial_{x})\psi_{+}$ , $\displaystyle
A_{22}=(1+\kappa^{-1}\partial_{x})\psi_{+}$ $\displaystyle
B_{21}=(1+\kappa^{-1}\partial_{x})\psi_{-}$ , $\displaystyle
B_{22}=(-1+\kappa^{-1}\partial_{x})\psi_{-}\quad$ (26)

where the $\psi_{\pm}$ satisfy the Helmholtz equation

$\partial_{x}^{2}\psi=\kappa^{2}(1+\chi)\psi\,.$ (27)

From the initial condition (23) and the structure (26) of the solution follows
the asymptotics

$\psi_{\pm}\propto e^{\pm\kappa x}\quad\mbox{for}\quad
x\rightarrow\mp\infty\,.$ (28)

Taken all this together we obtain

$\frac{T_{22}^{\prime}}{T_{22}}=-\kappa^{2}\chi\delta\,\frac{(\psi_{+}\psi_{-})^{\prime}}{W}$
(29)

where $W$ denotes the Wronskian

$W=\psi_{+}\psi_{-}^{\prime}-\psi_{-}\psi_{+}^{\prime}\,.$ (30)

Note that the Wronskian (30) is constant in space as a consequence of the
Helmholtz equation (27).

We have seen that the transfer across the scattering structure (and hence the
force) is governed by the macroscopic Helmholtz equation (27). This shows that
the $\chi$ defined in Eq. (1) is indeed the susceptibility and Eq. (2) gives
the refractive index. In our one–dimensional case Fiedler  the susceptibility
depends linearly on the density $\delta^{-1}$. In three dimensions, this is
normally only the case for dilute fluids.

###  III.2 Green function

In the following we relate the force density in the macroscopic limit to the
Green function. The Green function $g(x,x_{0})$ describes the field
propagation from a monochromatic point source of unit strength at position
$x_{0}$ to the observation point $x$, and thus obeys the differential equation

$(\partial_{x}^{2}-\kappa^{2}n^{2})g=\delta(x-x_{0})$ (31)

with the asymptotics of a purely outgoing wave at infinity:

$g\propto e^{\pm\kappa x}\quad\mbox{for}\quad x\rightarrow\mp\infty\,.$ (32)

Away from the point of emission ($x\neq x_{0}$) the Green function must be a
solution of the Helmholtz equation (27). At the emission point $g$ must be
continuous;

$g(x_{0}+0,x_{0})=g(x_{0}-0,x_{0})$ (33)

while the derivative of $g$ needs to jump by $1$ in order to generate the
delta function on the right–hand side of the Green equation (31):

$\left.\partial_{x}g(x,x_{0})\right|_{x=x_{0}+0}-\left.\partial_{x}g(x,x_{0})\right|_{x=x_{0}-0}=1\,.$
(34)

The well–known solution of these requirements is the expression

$g(x,x_{0})=\frac{1}{W}\begin{cases}\psi_{+}(x)\,\psi_{-}(x_{0})&:\quad x\leq
x_{0}\,,\\\ \psi_{+}(x_{0})\,\psi_{-}(x)&:\quad x\geq x_{0}\,.\end{cases}$
(35)

We turn to our calculation of the force in the macroscopic limit, and define
the force density as

$f=\frac{F}{\delta}\,.$ (36)

We obtain from Eqs. (9), (12), (29) and (35) the expressions

$f=\frac{\hbar
c}{2\pi}\int_{0}^{\infty}\widetilde{f}\,d\kappa\,,\quad\widetilde{f}=\kappa^{2}\chi\,\partial_{x}g(x,x)\,.$
(37)

This is one of the central results of this paper. Similar to the Lifshitz
theory DLP  of the Casimir force we have related the force density to the
electromagnetic Green function, but in our case this is the force density
inside dielectric media, not the force between dielectric bodies as in
Lifshitz’s case DLP . Another important difference is the fact that our $g$ is
the unrenormalized Green function — at least for the time being. We shall
discuss in Sec. V how renormalization emerges.

Figure 2a shows a simple example: the van der Waals forces inside a block of
uniform material of length $a$ and refractive index $n_{0}$ surrounded by
vacuum. We see that the force density peaks at the edges of the block,
strongly pulling the material in due to van der Waals attraction, whereas near
the middle of the block the van der Waals forces from the different parts of
material nearly compensate each other. We also see the good agreement between
the numerical calculation with the prediction (37) of the macroscopic theory,
even for quite a small number of point scatterers ($\geq 50$).

###  III.3 Asymptotics

In order to understand whether and when formula (37) for the force density
converges in the macroscopic limit, we apply the Madelung representation
Madelung  for the fields $\psi_{\pm}$. The Helmholtz equation (27) plays the
role of the Schrödinger equation in optics, and Madelung has reduced the
Schrödinger dynamics to fluid mechanics Madelung . Following Madelung Madelung
we make the ansatz:

$\psi_{\pm}=\frac{1}{\sqrt{2k}}\,\exp\left(\pm\int k\,dx\right)$ (38)

where the Wronskian (30) is exactly

$W=-1\,.$ (39)

From the Helmholtz equation (27) follows an equation for $k$ that is
independent on the $\pm$ sign in the exponent,

$k^{2}=n^{2}\kappa^{2}-2\beta\,,$ (40)

and appears as a generalized dispersion relation with the extra term

$\beta=\frac{\sqrt{k}}{2}\,\partial_{x}^{2}\,\frac{1}{\sqrt{k}}=-\frac{1}{4}\left(\frac{k^{\prime}}{k}\right)^{\prime}+\frac{1}{8}\left(\frac{k^{\prime}}{k}\right)^{2}.$
(41)

For the quantum–mechanical Schrödinger equation, the $\beta$–term corresponds
to Bohm’s quantum potential Bohm , in optics $\beta$ describes the scattering
in an inhomogeneous medium. Mathematically, $\beta$ corresponds to the
Schwarzian derivative AF  of the phase.

For large frequencies ($\kappa\rightarrow\infty$) relation (40) reduces to the
ordinary dispersion relation $k^{2}\sim n^{2}\kappa^{2}$ of geometrical
optics. We take the positive sign for $k$ such that the amplitude
$(2k)^{-1/2}$ in the Madelung representation (38) is real. In subdominant
order we take $k^{2}\sim[n\kappa-\beta/(n\kappa)]^{2}$ and have

$k\sim n\kappa-\frac{\beta_{0}}{n\kappa}\,.$ (42)

Here $\beta_{0}$ denotes the dominant contribution to $\beta$, i.e. Eq. (41)
evaluated for $k=n\kappa$. We obtain

$\beta_{0}=-\frac{1}{4}\left(\frac{n^{\prime}}{n}\right)^{\prime}+\frac{1}{8}\left(\frac{n^{\prime}}{n}\right)^{2}=-\frac{n^{\prime\prime}}{4n}+\frac{3n^{\prime
2}}{8n^{2}}$ (43)

that does not explicitly depend on $\kappa$. The term $\beta_{0}$ quantifies
the deviation from geometrical optics due to inhomogeneities in the medium,
the amplitude of geometrical scattering.

In the Madelung representation (38) we get for the Green function of Eq. (35)
evaluated at $x_{0}=x$ the simple expression $g(x,x)=-(2k)^{-1}$ and hence for
the spectral force density of Eq. (37):

$\widetilde{f}=-\kappa^{2}(n^{2}-1)\,\partial_{x}\frac{1}{2k}\,.$ (44)

For large frequencies (large imaginary wavenumbers $\kappa$) we obtain from
Eq. (42) the asymptotics:

$\widetilde{f}\sim-\frac{n^{2}-1}{2}\,\partial_{x}\left(\frac{\kappa}{n}+\frac{\beta_{0}}{\kappa
n^{3}}\right).$ (45)

For piece–wise homogeneous media (Figs. 2a and 6a) the derivative in Eq. (45)
vanishes, the integral (37) converges, and the force density is finite in the
macroscopic limit. At dielectric interfaces (discontinuities of $n$) the force
density does not converge, though. Moreover, any inhomoneity, however small,
would produce an infinite force density. The suspension of scatterers would
instantly and catastrophically collapse. This does not happen in reality. So
which mechanism — not considered in our model — prevents the catastrophe?

##  IV Helmholtz force

It is wise to take a step back and relate our result (37) to the classical
optical forces in dielectrics LL8 . We begin with a brief review on the
Helmholtz force in three–dimensional electromagnetism. Then we turn to our
one–dimensional model and identify the quantum forces and stresses there.

###  IV.1 Dielectric deformation

Consider three–dimensional electromagnetism in a dielectric medium of electric
permittivity $\varepsilon=n^{2}$ (and magnetic permeability $\mu=1$). Forces
in the medium cause the dielectric to deform. Consider an infinitesimal
deformation $\delta\bm{x}$. If we can write the variation of the free energy
$F$ as the volume integral

$\delta F=-\int\bm{f}\cdot\delta\bm{x}\,dV$ (46)

then $\bm{f}$ gives the force density, according to Hamilton’s equations. The
part of the free energy modified by the dielectric is the volume integral of
the electric energy density $U=\frac{1}{2}\bm{E}\cdot\bm{D}$ with
$\bm{D}=\varepsilon_{0}\varepsilon\bm{E}$ in SI units. Gauss’ law
$\nabla\cdot\bm{D}=\rho_{\rm ext}$ gives $\bm{D}$ in terms of the external
electrical charges that are not modified by the deformation of the medium. We
thus take $\bm{D}$ as fixed in the energy density and obtain $\delta
U=-U\delta\varepsilon/\varepsilon$. Two terms contribute to
$\delta\varepsilon$, the infinitesimal shift of the medium and the change in
density $\rho$, as $\varepsilon$ may depend on density:

$\delta\varepsilon=-\nabla\varepsilon\cdot\delta\bm{x}+\frac{\partial\varepsilon}{\partial\rho}\,\delta\rho$
(47)

Differential geometry Poisson  gives
$\delta\rho=-\rho\nabla\cdot\delta\bm{x}$. Integrating by parts in the energy
produces the structure of Eq. (46) with the force density LL8

$\bm{f}=-\frac{\nabla\varepsilon}{\varepsilon}\,\frac{\bm{E}\cdot\bm{D}}{2}+\nabla
p_{\rm Ab}$ (48)

and

$p_{\rm
Ab}=\frac{\partial\varepsilon}{\partial\rho}\,\rho\,\frac{\bm{E}\cdot\bm{D}}{2\varepsilon}\,.$
(49)

The pressure $p_{\rm Ab}$ — known as the Abraham pressure LL8  — originates
from the effect of the density variation on the dielectric: The total force
(48) on the medium is called the Helmholtz force LL8 . For dilute media and in
our one–dimensional case, $\varepsilon=1+\chi$ with
$\chi=\alpha\rho/\varepsilon_{0}$ such that we get for the Helmholtz force:

$\bm{f}=\chi\nabla\,\frac{\varepsilon_{0}E^{2}}{2}=\rho\nabla\,\frac{\alpha
E^{2}}{2}\,,$ (50)

which is the dipole or gradient force familiar from optical tweezing Ashkin1 ;
Ashkin2 ; Ashkin3 ; Neuman ; Pesce .

###  IV.2 Quantum force

Let us relate our case of van der Waals forces on suspended particles to the
Helmholtz force in dielectrics. For this we need to treat the electromagnetic
field as a quantum field with the operator $\widehat{A}$ of the vector
potential as fundamental quantity. The electric field strength $\widehat{E}$
is given by $-\partial_{t}\widehat{A}$ in Coulomb gauge. Consider the
correlation functions Annals

$\displaystyle K$ $\displaystyle=$
$\displaystyle\frac{\varepsilon_{0}c}{2\hbar}\langle\widehat{A}_{1}\widehat{A}_{0}+\widehat{A}_{0}\widehat{A}_{1}\rangle=\frac{\varepsilon_{0}c}{\hbar}\,{\rm
Re}\,\langle\widehat{A}_{1}\widehat{A}_{0}\rangle\,,$ (51)
$\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle\frac{\varepsilon_{0}c}{2i\hbar}\langle\widehat{A}_{1}\widehat{A}_{0}-\widehat{A}_{0}\widehat{A}_{1}\rangle=\frac{\varepsilon_{0}c}{\hbar}\,{\rm
Im}\,\langle\widehat{A}_{1}\widehat{A}_{0}\rangle$ (52)

where the indices refer to the positions $x_{i}$ and times $t_{i}$. Note that
$\langle\widehat{A}_{1}\widehat{A}_{0}\rangle$ is an analytic function of the
time difference $t=t_{1}-t_{0}$ if the field is in the vacuum state (or
thermal) Annals . For an analytic function the real and the imaginary part are
connected AF . Furthermore, $\Gamma$ corresponds to the difference between the
retarded and the advanced Green function Annals , the difference between
outgoing and ingoing radiation: the dissipation. The correlation $K$ of the
field fluctuations is thus related to the dissipation $\Gamma$, which is a
form of the fluctuation–dissipation theorem Scheel .

In particular, in terms of our time–independent Green function $g$ and for
real frequencies $\omega$ we have for the Fourier transform: Annals

$\widetilde{\Gamma}=\frac{g(\omega)-g(-\omega)}{2c}\,.$ (53)

We see that $\widetilde{\Gamma}$ is an odd function of $\omega$. Consequently,
we get for the inverse Fourier transform:

$\Gamma(t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\widetilde{\Gamma}e^{-i\omega
t}\,d\omega=-\frac{i}{\pi}\int_{0}^{\infty}\widetilde{\Gamma}\sin\omega
t\,d\omega\,.$ (54)

For finding the real part $K$ to the imaginary part $\Gamma$ we only need to
replace $-\sin\omega t$ by $\cos\omega t$, because they represent the
imaginary and real part of the analytic $e^{-i\omega t}$ of time evolution.
Remark  Using Eq. (53) and deforming the integration contours to the positive
imaginary axis with wavenumbers $\kappa$ (Fig. 3) we obtain:

$K(t)=-\frac{1}{\pi}\int_{0}^{\infty}g\cosh(c\kappa t)\,d\kappa\,.$ (55)

Strictly speaking, formula (55) is only valid when $g$ falls off more rapidly
than $\cosh(c\kappa t)$, which is the case outside of the light cone,
$|x_{1}-x_{0}|>c|t|$, but only there we will need $K$.

From the correlation function $K$ we can calculate the electromagnetic field
correlations by differentiation. In particular, we get for the electric fields
at equal times, but different positions:

$\frac{\varepsilon_{0}}{2}\langle\widehat{E}_{1}\widehat{E}_{0}\rangle=\left.\frac{\hbar}{2c}\,\partial_{t_{1}}\partial_{t_{0}}K\right|_{t_{1}=t_{0}}=\frac{\hbar
c}{2\pi}\int_{0}^{\infty}\kappa^{2}g(x_{1},x_{0})\,d\kappa\,.$ (56)

Compare this with the result (37) for the force density in the macroscopic
limit and Eq. (50) for the Helmholtz force: the macroscopic van der Waals
force is nothing else but the gradient force Ashkin1 ; Ashkin2 ; Ashkin3 ;
Neuman ; Pesce . It is a peculiar gradient force, though. The force is not
generated by external fields, but rather by the zero–point fluctuations of the
field inside the medium.

###  IV.3 Electromagnetic stresses

Having established the force density, we are going to relate it now to the
electric and magnetic stress tensors LL8

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

From Maxwell’s equations and the linear constitutive equations
$\bm{D}=\varepsilon_{0}\varepsilon\bm{E}$ and
$\bm{H}=\varepsilon_{0}c^{2}\bm{B}$ follows the Abraham identity LL8  for
stationary fields:

$\nabla\cdot(\sigma_{E}+\sigma_{M})=-\frac{\nabla\varepsilon}{\varepsilon}\,\frac{\bm{E}\cdot\bm{D}}{2}$
(58)

and from this we get for the Helmholtz force (48):

$\bm{f}=\nabla\cdot(\sigma_{E}+\sigma_{M})+\nabla p_{\rm Ab}\,.$ (59)

The force density is thus the divergence of the total electromagnetic stress
plus the gradient of the Abraham pressure.

So far we considered the classical stress, but similarly to Sec. IVB we can
identify the quantum electromagnetic stresses as well. These stresses are the
vacuum expectation values of expressions (57) for the field operators,
symmetrized as in definition (51) of the correlation function $K$. From
$\widehat{E}=-\partial_{t}\widehat{A}$ and
$\widehat{B}=\partial_{x}\widehat{A}$ (in one dimension) we obtain:

$\widetilde{\sigma}_{E}=-n^{2}\kappa^{2}g(x,x)\,,\quad\widetilde{\sigma}_{M}=\left.\partial_{x}\partial_{x_{0}}g\right|_{x_{0}=x}$
(60)

for the spectral stresses defined as in Eq. (37). In terms of the Madelung
representation [Eqs. (35), (38) and (39)] we get:

$\widetilde{\sigma}_{E}=\frac{n^{2}\kappa^{2}}{2k}\,,\quad\widetilde{\sigma}_{M}=\frac{k}{2}-\frac{k^{\prime
2}}{8k^{3}}$ (61)

with the wavenumber $k$ satisfying Eqs. (40) and (41). We also express the
quantum Abraham pressure (49) in terms of the Green function. In our
one–dimensional case we have $\rho\,\partial\varepsilon/\partial\rho=\chi$ and
so we obtain from Eq. (56) for the spectral Abraham pressure:

$\widetilde{p}_{\rm
Ab}=\chi\kappa^{2}g(x,x)=-\frac{n^{2}-1}{n^{2}}\,\widetilde{\sigma}_{E}\,.$
(62)

One verifies that from Eqs. (40) and (41) follows the Abraham identity for the
quantum stresses (61):

$\partial_{x}(\widetilde{\sigma}_{E}+\widetilde{\sigma}_{M})=\frac{2n^{\prime}}{n}\,\widetilde{\sigma}_{E}\,.$
(63)

Formulas (61) reveal the asymptotics of the stresses for large frequencies
($\kappa\rightarrow\infty$). Equation (42) gives

$\widetilde{\sigma}_{E}\sim\frac{\kappa n}{2}+\frac{\beta_{0}}{2\kappa
n}\,,\quad\widetilde{\sigma}_{M}\sim\frac{\kappa
n}{2}-\frac{\beta_{0}}{2\kappa n}-\frac{n^{\prime 2}}{8\kappa n^{3}}$ (64)

such that we obtain for the total spectral stress,

$\widetilde{\sigma}=\widetilde{\sigma}_{E}+\widetilde{\sigma}_{M}\,,$ (65)

the asymptotics:

$\widetilde{\sigma}\sim\kappa n-\frac{n^{\prime 2}}{8\kappa n^{3}}\,.$ (66)

Like in the three–dimensional planar case Renormalization  the asymptotics of
the total stress depends only on the local dielectric functions and their
first derivatives, in our case on the refractive index $n$ and its derivative
$n^{\prime}$ (Fig. 4). We also see that the integrated stresses diverge even
in the case of piece–wise homogeneous media when the force density (37)
converges.

Figure 4:  Asymptotics of the stress. The figure illustrates the asymptotics
of the spectral stress [Eq. (65)] as a function of $\kappa$ for the
susceptibility profile [Eq. (86)] displayed in Fig. 2b (evaluated at $x=0.1$).
Parameters as in Fig. 2b. The gray curve shows $\widetilde{\sigma}-\kappa n$,
the renormalization proposed by Dzyaloshinskii and Pitaevskii DP . One sees
that the curve gently decays (as $\kappa^{-1}$). The total renormalized stress
is given by the integral over $\kappa$ and diverges. For the black curve the
asymptotics (66) is subtracted; the curve decays as $\kappa^{-3}$ and this
renormalization (Fig. 5) converges.

##  V Renormalization

Following the scattering approach to Casimir forces Lambrecht ; Wirzba ;
Ingold  we have developed a theory that directly gives the actual forces on
point scatterers, without the need of renormalization. We found that the
density of the van der Waals force and the electromagnetic stresses diverge in
the macroscopic limit, in direct contradiction to the empirical fact that some
colloidal suspensions are stable — that liquid paint may remain liquid. Some
additional feature must prevent this catastrophic collapse pure van der Waals
forces would cause.

###  V.1 Argument

Pitaevskii Pitaevskii  argued that the pressure of the fluid naturally
counteracts the van der Waals forces. These counter forces depend on the
density and are typically much stronger than the van der Waals forces, such
that a small variation in density is sufficient to reach a local mechanical
equilibrium. As this small density variation hardly affects the van der Waals
force, we ignore it in calculating the force, but we do use the fact that the
counter pressure it causes partly compensates the van der Waals force. This
counter pressure appears as the renormalization of the Casimir force.

Which contributions to the van der Waals force are compensated? Pitaevskii
Pitaevskii  points out the following. As the counter pressure reacts to
density variations, it immediately takes out the part of the quantum force
that is also generated by density variations (Sec. IVA): the Abraham pressure
[Eq. (49)]. What remains is the divergence of the electromagnetic stress in
the material Pita ; Burger ; Burger2 . But, as we argue here, this is not all.
The counter pressure can also react to local variations of the electromagnetic
stresses, establishing local mechanical equilibria everywhere.

We need to define what we mean by local. The electromagnetic stresses are
created by zero–point fluctuations of the field in inhomogeneous media. We may
characterize the length scale $\ell$ of the inhomogeneity by
$\ell^{-1}=|n^{\prime}/n|$. Only frequencies with wavelengths $\lambda\ll\ell$
will contribute to the local stresses. As $\lambda$ is inversely proportional
to $n$ we obtain the criterion

$|\nabla\lambda|\ll 1\quad\mbox{with}\quad\lambda=\frac{2\pi}{\kappa n}\,,$
(67)

which is the well–known condition for the validity of geometrical optics or of
the semiclassical approximation in quantum mechanics LL3 .

We can only define local stresses for fields consistent with geometrical
optics. The field at zero frequency will never satisfy this condition, so
whenever zero–frequency features are important, the renormalization procedure
is bound to fail. This is the case for thermal Casimir forces dominated by the
zero Matsubara frequency LL7 , which might explain why Lifshitz theory DLP
does not agree with experiments there LifshitzDrude .

###  V.2 Local stresses

Next we assume criterion (67) and work out the local electromagnetic stresses.
According to Eq. (60) the stresses are generated by the Green function. We
must define a new Green function that captures their local contributions. This
$g_{0}$ we only need in the vicinity of the point of emission $x_{0}$ and
there we assume Renormalization  that it describes a purely outgoing wave
unaffected by the refractive index profile away from $x_{0}$. Here is the
argument for $g_{0}$:

Green functions must satisfy two conditions: the propagation equation (31) and
the boundary conditions. The total Green function $g$ is required by the
boundary condition (32) to be purely outgoing at $\pm\infty$. For $g_{0}$ we
require that it is purely outgoing in an infinitesimal interval around $x_{0}$
where outgoing we define in the sense of geometrical optics. To the left or
right of $x_{0}$ it should propagate like the $\psi_{-}$ or $\psi_{+}$ waves
in the Madelung representation (38) with $k=n\kappa$. They obey the condition
$\mathbb{D}_{\pm}\psi_{\pm}=0$ with the differential operator

$\mathbb{D}_{\pm}=\partial_{x}\pm\kappa n+\frac{n^{\prime}}{2n}\,.$ (68)

The $\partial_{x}\pm\kappa n$ term captures the propagation of the phase,
while the $n^{\prime}/(2n)$ term takes care of the amplitude variation in
geometrical optics. We thus encode the boundary condition as

$\mathbb{D}_{\pm}g_{0}=0\quad\mbox{for}\quad x=x_{0}\pm 0\,.$ (69)

Being a Green function, $g_{0}$ must solve the propagation equation (31). For
$x\neq x_{0}$ it is therefore a solution of the Helmholtz equation (27), i.e.
a superposition of the fundamental system $\psi_{\pm}$ of wave functions:

$g_{0}=\begin{cases}c_{+-}\psi_{+}(x)+c_{--}\psi_{-}(x)&:\quad x<x_{0}\,,\\\
c_{++}\psi_{+}(x)+c_{-+}\psi_{-}(x)&:\quad x>x_{0}\,.\end{cases}$ (70)

Note that the coefficients depend on $x_{0}$. At the point of emission,
$x=x_{0}$, the Green function $g_{0}$ must obey the continuity condition (33)
and the discontinuity (34) of the derivative. The four conditions of Eqs.
(33-34) and (68-69) determine the four coefficients of $g_{0}$ uniquely. We
obtain

$c_{\pm a}=\left.\pm\frac{\mathbb{D}_{a}\psi_{\mp}}{2\kappa
n}\right|_{x_{0}}\quad\mbox{for}\quad a=\pm\,.$ (71)

With these coefficients for $g_{0}$ we calculate the spectral stresses (60)
using the Helmholtz equation (27) for expressing $\psi_{\pm}^{\prime\prime}$
as $\kappa^{2}n^{2}\psi_{\pm}$. We obtain the exact result

$\widetilde{\sigma}_{E0}=\frac{\kappa
n}{2}\,,\quad\widetilde{\sigma}_{M0}=\frac{\kappa
n}{2}+\partial_{x}\frac{n^{\prime}}{4\kappa n^{2}}\,.$ (72)

We see that $\widetilde{\sigma}_{E0}$ and $\widetilde{\sigma}_{M0}$ depend
only on the local refractive index $n$ and its first and second derivative.
The local Green function $g_{0}$ does indeed capture local contributions to
the electromagnetic stresses. Appendix B compares $g_{0}$ with the previous
renormalizer Itai  and finds good agreement.

###  V.3 Anomaly

We have seen that Eq. (72) describes local stresses, but as Eq. (64)
indicates, they do not completely account for the asymptotics of the full
quantum electromagnetic stresses at large frequencies
($\kappa\rightarrow\infty$). We use Eq. (43) and write the full asymptotics
(64) in terms of the local stresses (72) as

$\widetilde{\sigma}_{E}\sim\widetilde{\sigma}_{E0}+\frac{\beta_{0}}{2\kappa
n}\,,\quad\widetilde{\sigma}_{M}\sim\widetilde{\sigma}_{M0}+\frac{\beta_{0}}{2\kappa
n}\,.$ (73)

We see that in both the electric and the magnetic stress there is a
characteristic difference by $\beta_{0}/(2\kappa n)$. For dispersive media
$\beta_{0}$ will decay sufficiently fast such that the integral of the
difference converges, but there is another, deeper problem: the stresses (72)
do not satisfy the Abraham identity (63). This is because the full stresses
do, and so do the asymptotic stresses (73), but not the terms
$\beta_{0}/(2\kappa n)$. Now, the Abraham identity describes momentum
conservation and it is a consequence of Maxwell’s equations. So why is it
violated here?

In our one–dimensional case, Maxwell’s equations reduce to the Helmholtz
equation. One can derive the Abraham identity (63) for the stresses (60) from
the Helmholtz equation for the Green function [Eq. (31) at $x\neq x_{0}$] and
the reciprocity of the Green function:

$g(x,x_{0})=g(x_{0},x)\,.$ (74)

The reciprocity appears in the structure (35) of the Green function, and it
stems from the boundary conditions (32). They state that the field is purely
outgoing at infinity. Therefore, the only source of radiation is the emission
at point $x_{0}$. If there are no other sources from infinity we can swap the
positions of emitter and receiver and get the same. The local Green function,
on the other hand, does satisfy the Helmholtz equation as well, but not the
reciprocity, because in order to sustain its boundary conditions (68) and (69)
sources at $\pm\infty$ are required. Another, direct argument is the
following. For $g_{0}$ the boundary conditions (68-69) depend on the
dielectric environment at the point of emission $x_{0}$. So the emission of
radiation depends on $n(x_{0})$ and $n^{\prime}(x_{0})$, but not the
reception, which breaks the reciprocity (74).

In order to restore momentum conservation [the Abraham identity (63)] we would
need to add $\beta_{0}/(2\kappa n)$ to the spectral stresses
$\widetilde{\sigma}_{E0}$ and $\widetilde{\sigma}_{M0}$. The extra term is
called an anomaly Wald . Wald Wald  predicted it first for quantum fields in
curved space BD . He derived a trace anomaly from the lack of reciprocity of
the renormalizer if there is space–time curvature. Inhomogeneous dielectric
media correspond to space–time geometries and vice versa Gordon ; Plebanski ;
Schleich ; LeoPhil . So it is perhaps not surprising that we obtain an anomaly
as well. Note however that the van der Waals anomaly is quantitatevely
different from Wald’s anomaly Wald . It was first deduced Itai  for
three–dimensional media with spherical symmetry; here we have derived it for
the simple one–dimensional case as well.

For quantum fields in dielectric media, the anomaly describes the recoil
pressure in the material LondonProceedings . The picture is this: Each
infinitesimal cell of the material is emitting and receiving virtual
electromagnetic waves described by the Green function $g$, which creates the
electromagnetic stresses. The local contributions to the stresses depend on
the local Green function $g_{0}$, but for $g_{0}$ emission and reception are
no longer interchangeable. The lack of reciprocity creates an imbalance in the
recoil of emission and reception that appears as an additional pressure: the
anomaly.

Equations (72) derived from the local Green function plus the anomaly capture
the divergent contributions to the electromagnetic stresses. As they are local
the counter pressure of the suspension should be able to compensate for them.
What remains are the effective electromagnetic stresses:

$\widetilde{\sigma}_{F}^{\rm
eff}=\widetilde{\sigma}_{F}-\widetilde{\sigma}_{F0}-\frac{\beta_{0}}{2\kappa
n}\quad\mbox{for}\quad F=E,M$ (75)

with $\widetilde{\sigma}_{F0}$ given by Eq. (72) and $\beta_{0}$ by Eq. (43).
We have thus deduced the renormalizing counter stresses and the anomaly from
condition (69). Note that we should introduce a cutoff at low frequencies for
being consistent with the criterion (67) for locality. In three dimensions
Itai  such a cutoff would not be necessary as the density of states goes with
$\kappa^{2}$, which sufficiently suppresses the low–frequency contribution.
There the criterion (67) only matters if the zero frequency becomes singled
out and dominant, as in the thermal Casimir effect LifshitzDrude .

Figure 5:  Renormalized stress. The figure shows the renormalized total
spectral stress $\widetilde{\sigma}^{\rm eff}=\widetilde{\sigma}_{E}^{\rm
eff}+\widetilde{\sigma}_{M}^{\rm eff}$ with the effective stresses given by
Eq. (75) for the susceptibility profile $\chi$ (gray filled curve) of Eq. (86)
for several values of $\kappa$ (from top to bottom: $\kappa$ from 3 to 10 in
steps of 1). Figure 2b shows the bare van der Waals forces, whereas here we
see the effect of local counter forces (for the same parameters). The
effective stress resembles the box–like profile [Eq. (83)] for the stress in a
homogeneous layer, apart from swings at the edges. As the stress rises the
force is attractive (whereas in Fig. 6b the stress drops and the Casimir force
is repulsive).

##  VI Examples

Let us consider two examples to illustrate our ideas: first a piece–wise
homogeneous dielectric of three layers and then a genuinely inhomogeneous
medium with sech–squared profile of the susceptibility. The three layers
constitute the classic case DLP  of Casimir forces. We will revisit this case
and see in detail how the renormalization in dielectrics emerges there.

###  VI.1 Three layers

Consider a homogenous dielectric layer of thickness $a$ and refractive index
$n_{2}$ sandwiched between two infinitely thick layers of homogeneous media
with $n_{1}$ on the left and $n_{3}$ on the right. We put the coordinates such
that the left interface lies at $x=0$ and the right interface at $x=a$. For
calculating the force density according to Eq. (37) we need the Green function
at equal positions. In the central layer we have (Appendix C):

$\displaystyle g(x,x)$ $\displaystyle=$
$\displaystyle-\frac{1}{2n_{2}\kappa}-\frac{\varrho_{l}\,e^{-2n_{2}\kappa
x}+\varrho_{r}\,e^{2n_{2}\kappa(x-a)}}{2n_{2}\kappa\,(1-\varrho_{l}\varrho_{r}\,e^{-2n_{2}\kappa
a})}$ (76)
$\displaystyle-\frac{1}{n_{2}\kappa}\,\frac{1}{(\varrho_{l}\varrho_{r})^{-1}e^{2n_{2}\kappa
a}-1}$

in terms of the Fresnel reflection coefficients

$\varrho_{l}=\frac{n_{2}-n_{1}}{n_{2}+n_{1}}\quad\mbox{and}\quad\varrho_{r}=\frac{n_{2}-n_{3}}{n_{2}+n_{3}}\,,$
(77)

whereas in the left layer we obtain (Appendix C):

$g(x,x)=-\frac{1}{2n_{1}\kappa}+\frac{1}{2n_{1}\kappa}\,\frac{\varrho_{l}-\varrho_{r}\,e^{-2n_{2}\kappa
a}}{1-\varrho_{l}\varrho_{r}\,e^{-2n_{2}\kappa a}}\,e^{2n_{1}\kappa x}\,.$
(78)

For the Green function in the right layer we only need to replace in Eq. (78)
the index $n_{1}$ by $n_{3}$ and the term $e^{2n_{1}\kappa x}$ by
$e^{2n_{3}\kappa(a-x)}$ and to exchange $\varrho_{l}$ and $\varrho_{r}$. We
now have all expressions ready for calculating the force density (37).

Assume that the refractive indices do not depend on frequency (no dispersion).
In this case we can analytically perform the integration in formula (37) in
terms of the Lerch transcendent [Sec. 1.11 of Ref. Erdelyi, ]:

$\Phi(z,s,x)=\sum_{m=0}^{\infty}\frac{z^{m}}{(m+x)^{s}}\,,$ (79)

because $\Phi$ has the integral representation 1.11(3) of Ref. Erdelyi, :

$\Phi(z,s,x)=\frac{1}{\Gamma(s)}\int_{0}^{\infty}\frac{t^{s-1}e^{-xt}}{1-ze^{-t}}\,dt\,.$
(80)

We obtain in the three regions

$\displaystyle f_{1}$ $\displaystyle=$ $\displaystyle\frac{\hbar
c}{\pi}\frac{n_{1}^{2}-1}{(2n_{2}a)^{3}}\left[\varrho_{l}\,\Psi\Big{(}-\frac{n_{1}x}{n_{2}a}\Big{)}-\varrho_{r}\,\Psi\Big{(}1-\frac{n_{1}x}{n_{2}a}\Big{)}\right],$
$\displaystyle f_{2}$ $\displaystyle=$ $\displaystyle\frac{\hbar
c}{\pi}\frac{n_{2}^{2}-1}{(2n_{2}a)^{3}}\left[\varrho_{l}\,\Psi\Big{(}\frac{x}{a}\Big{)}-\varrho_{r}\,\Psi\Big{(}1-\frac{x}{a}\Big{)}\right],$
(81) $\displaystyle f_{3}$ $\displaystyle=$ $\displaystyle\frac{\hbar
c}{\pi}\frac{n_{3}^{2}-1}{(2n_{2}a)^{3}}\left[\varrho_{l}\,\Psi\Big{(}1-\frac{n_{3}x_{a}}{n_{2}a}\Big{)}-\varrho_{r}\,\Psi\Big{(}-\frac{n_{3}x_{a}}{n_{2}a}\Big{)}\right]$

with $x_{a}=a-x$ while $\Psi(\xi)$ abbreviates

$\Psi(\xi)=\Phi(\varrho_{l}\varrho_{r},3,\xi)\sim\xi^{-3}\quad\mbox{for}\quad\xi\rightarrow
0\,.$ (82)

The force density is sharply peaked at the dielectric interfaces where the
pull of the van der Waals forces is strongest (Figs. 2a and 6a). However,
since $\chi$ is spatially constant, these forces [Eq. (37)] are entirely given
by the Abraham pressure [Eq. (62)]. We have argued in Sec. VA that counter
forces in the material completely compensate this pressure.

What is left are the stress differences at the interfaces. From the Green
equation (31) follows that $g$ is some superposition of $e^{\pm n\kappa x}$ in
each of the three layers and for $x\neq x_{0}$. Reciprocity [Eq. (74)] implies
that $g(x,x_{0})$ consists of terms proportional to $e^{\pm n\kappa(x-x_{0})}$
and $e^{\pm n\kappa(x+x_{0})}$. The $e^{\pm n\kappa(x+x_{0})}$ terms have
produced the position–dependent parts in $g(x,x)$ and hence the Abraham
pressure. We see from Eq. (60) that these terms cancel each other in the sum
(65) of the electromagnetic stresses, whereas the $e^{\pm n\kappa(x-x_{0})}$
terms generate equal electric and magnetic stresses. In the limit
$x_{0}\rightarrow x$ these are the position–independent parts of Eqs. (76) and
(78). For them, twice the electric stress $-n^{2}\kappa^{2}g$ gives the total
stress:

$\displaystyle\widetilde{\sigma}_{1}=n_{1}\kappa\,,\quad\widetilde{\sigma}_{3}=n_{3}\kappa\,,$
$\displaystyle\widetilde{\sigma}_{2}=n_{2}\kappa+\frac{2n_{2}\kappa}{(\varrho_{l}\varrho_{r})^{-1}e^{2n_{2}\kappa
a}-1}\,.$ (83)

The stress is constant in each of the dielectric layers, but discontinuous at
the interfaces, generating there the localized force densities

$f=F_{l}\,\delta(x)+F_{r}\,\delta(x-a)\,.$ (84)

The dominant contributions to the stress (83) are the $n\kappa$ terms in
agreement with the asymptotics (64). But these are entirely local terms that
only depend on the local refractive index. We argued that they are also
compensated by the local counter pressure. What remains are the non–local
contributions that depend on the distance $a$ and the reflections
coefficients, the Casimir–Lifshitz forces DLP :

$F_{l}=-F_{r}=\frac{\hbar
c}{\pi}\int_{0}^{\infty}\frac{n_{2}\kappa\,d\kappa}{(\varrho_{l}\varrho_{r})^{-1}e^{2n_{2}\kappa
a}-1}\,.$ (85)

The van der Waals forces of Eq. (81) are attractive towards regions of higher
refractive index and peak at the dielectric interfaces (Figs. 2a and 6a). The
counter forces in the material compensate them completely, except at the
interfaces. The resulting Casimir–Lifshitz force of Eq. (85) may be attractive
or repulsive, depending on the sign of the product $\varrho_{l}\varrho_{r}$ of
the Fresnel reflection coefficients (77). For $n_{1}<n_{2}<n_{3}$ or
$n_{1}>n_{2}>n_{3}$ the force is repulsive DLP  (Fig. 6b). One of the layers
may even levitate on quantum fluctuations Levitation ; CasimirEquilibrium .

For the classic case of three homogeneous layers — and from our perspective on
renormalization — we have derived Eq. (85): Lifshitz’ formula.

Figure 6:  Three layers. a: The piece–wise homogeneous susceptibility profile
$\chi$ (filled steps) generates the van der Waals forces $F$ (dots) on each
point particle (in units of $\hbar c$). The solid curve shows the macroscopic
limit of the forces [Eq. (81)]. The forces are attractive towards regions of
higher susceptibility and peak at the interfaces (see also Fig. 2a). b:
Counter forces in the material remove most of the van der Waals forces and
produce the stress $\sigma$ that is zero in the outer layers and reaches the
Lifshitz value [Eq. (85)] in the central layer. The discontinuities of the
stress generate the Casimir–Lifshitz force at the interfaces [Eq. (84)] — a
repulsive force here.

###  VI.2 Inhomogeneous profile

Let us now turn to an inhomogeneous medium where standard Lifshitz
renormalization fails Simpson . We take an example where we can still
calculate the Green function analytically, and compare the resulting forces
and stresses with the asymptotics for large wavenumbers and with numerical
results for a finite number of particles. Consider the susceptibility

$\chi=\chi_{0}\sech^{2}(x/a)$ (86)

as a function of the spatial coordinate $x$ in units of the characteristic
scale $a$ of the profile. The Helmholtz equation (27) has the solution

$\psi_{+}=e^{\kappa
x}{}_{2}F_{1}(\nu_{-},\nu_{+},\kappa+1,\zeta)\,,\quad\zeta=\frac{1+\tanh(x/a)}{2}$
(87)

with the parameters

$\nu_{\pm}=\frac{1}{2}\left(1\pm\sqrt{1-4\chi_{0}\,a^{2}\kappa^{2}}\right)$
(88)

in terms of Gauss’ hypergeometric function Erdelyi  ${}_{2}F_{1}$ as $\zeta$
reduces Eqs. (27) and (86) to the hypergeometric differential equation Erdelyi
. Note that the $\nu_{\pm}$ are complex for $\kappa>1/(2a\sqrt{\chi_{0}})$
where $\nu_{\mp}$ becomes the complex conjugate of $\nu_{\pm}$. However, as
the hypergeometric series Erdelyi  goes equally with $\nu_{-}$ and $\nu_{+}$
all the terms in the series are real and so Eq. (87) always describes a real
function. For $x\rightarrow-\infty$ the variable $\zeta$ tends to zero. Since
all hypergeometric functions are unity at zero argument the wave $\psi_{+}$
approaches $e^{\kappa x}$ and so satisfies the boundary condition (28) as
required. For constructing the second solution $\psi_{-}$ of the wave equation
and of the opposite asymptotics ($\psi_{-}\sim e^{-\kappa x}$ for
$x\rightarrow+\infty$) we take advantage of the mirror symmetry of the profile
(86) and put

$\psi_{-}(x)=\psi_{+}(-x)\,.$ (89)

For the Green function (35) we need the Wronskian (30). The Wronskian is
constant throughout $x$ and so we may calculate it where the calculation is
most convenient, for $x\rightarrow+\infty$ where $\psi_{-}\sim e^{-\kappa x}$
and

$\displaystyle\psi_{+}$ $\displaystyle\sim$ $\displaystyle\frac{\Gamma(-\kappa
a)\,\Gamma(-\kappa a-1)}{\Gamma(\kappa a+\nu_{+})\,\Gamma(\kappa
a+\nu_{-})}\,e^{-\kappa x}$ (90) $\displaystyle+\frac{\Gamma(\kappa
a)\,\Gamma(\kappa a+1)}{\Gamma(\kappa a+\nu_{+})\,\Gamma(\kappa
a+\nu_{-})}\,e^{\kappa x}$

according to Eq. e7 of Ref. LL3, . Here $\Gamma$ denotes the Gamma function
Erdelyi  with $z\Gamma(z)=\Gamma(z+1)$. We thus obtain

$W=-\frac{2[\Gamma(\kappa a+1)]^{2}}{a\Gamma(\kappa a+\nu_{+})\,\Gamma(\kappa
a+\nu_{-})}\,.$ (91)

Figure 2b shows the bare van der Waals forces on a finite number $N$ of
suspended particles. The macroscopic limit does only exist for the spectral
force (Fig. 7) and the spectral stresses $\widetilde{\sigma}_{F}$, but not for
the total force density $f$ of Eq. (37) nor for the total stresses integrated
over $\kappa$. Figure 4 compares the spectral stresses in the macroscopic
limit with the asymptotics (64) for large wavenumbers, confirming our result
(66). Figure 5 illustrates the main result of this paper, the renormalized
electromagnetic stresses, evaluated according to Eq. (75) for the example
considered here.

Figure 7:  Spectral forces. Comparison of the spectral force (dots)
$-\partial_{j}\ln T_{22}$ [the integrand of the total force (97)] with the
macroscopic limit (solid curve) [Eq. (37)] and the asymptotics (dashed curve)
[Eq. (45)] for the susceptibility profile $\chi$ of Eq. (86) (gray filled
curve). Parameters as in Fig. 2b, $\kappa=8.0$. The curves eventually merge
for large $\kappa$, but are distinct yet at the moderate $\kappa$ used, for
better visibility.

##  VII Conclusions

We have analyzed a simple model for the Casimir force in dielectrics: a chain
of point particles interacting with each other by scattering virtual
electromagnetic waves. The scattering creates forces Buhmann ; Rodriguez ;
Forces  on the particles that appear as the familiar van der Waals forces
ShahmoonvdW  (with retardation taken into account CasimirPolder ). The point
scatterers may represent the molecules of a dielectric medium LL8  and they do
represent colloidal suspensions Israelachvili  such as liquid paint. In our
model, we assume both the particles and the field to be confined to one
dimension, which has allowed us to develop an efficient numerical algorithm
for calculating the forces, and to deduce the macroscopic limit in
unprecedented detail and clarity.

In the macroscopic limit of infinitely many, infinitely weak scatterers, we
found a result that seems both satisfying and paradoxical: the van der Waals
force becomes the gradient force of electric field fluctuations, which nicely
connects the scattering approach to Casimir forces Lambrecht ; Wirzba ; Ingold
with macroscopic quantum electrodynamics Buhmann ; KSW  and the rich field of
optical forces Ashkin1 ; Ashkin2 ; Ashkin3 ; Neuman ; Pesce . If the
scatterers are arranged to form a piece–wise homogeneous medium the force
density is finite, except at interfaces. Yet any inhomogeneity, however
minute, will produce infinite forces and bring about the immediate collapse of
the suspension. The scatterers appear to anticipate whether they are aligned
to form pieces of a homogeneous or of an inhomogeneous medium, which seems
paradoxical. It also disagrees with the empirical fact that some colloidal
suspensions are stable — liquid paint may remain liquid even if it is not
perfectly smooth. Homogeneous colloidal suspensions are known to be stable
when the long–range attractive van der Waals forces are compensated by other
short–range repulsive forces Israelachvili .

We have argued that the same mechanism that makes homogeneous colloidal
suspensions stable does also prevent the collapse of inhomogeneous media. The
counter forces compensate the local contributions to the van der Waals stress,
and the remaining non–local contributions constitute the Casimir force. This
local balance between the van der Waals and the counter force amounts to the
renormalization of the quantum–electromagnetic stresses. We have shown [Eqs.
(72) and (75)] that the counter stress depends on the local macroscopic
dielectric response functions, the refractive index and its first and second
derivative, independent of the microscopic physics involved.

Our model describes a suspension of particles in a fluid, but we could also
apply it to solids where the inter–atomic forces are even stronger. In either
case, the result of the counter pressure is independent of the material
details. Only the refractive–index profile counts. Suppose we replace the
dielectric material by empty space. Even in this case we need to remove the
local stress. For example (Sec.VIA) in each homogenous layer we should remove
$n\kappa$ in the spectral stress [Eq. (83)] — also in empty space where $n=1$.
Only then the discontinuity of the stress gives the correct Casimir–Lifshitz
force Lifshitz  [Eq. (85)]. One might speculate whether there is some
mechanism of counter forces in empty space as well, acting not on the
molecular scale of $10^{-9}\mathrm{m}$, but rather on the Planck scale of
$10^{-35}\mathrm{m}$. The precise microscopic nature of such a space–time
medium is unknown, but as in ordinary dielectric materials, the
renormalization should not depend on microscopic details.

In curved space–time, i.e. in gravitational fields LL2  the relation between
space and media is well–known Gordon ; Plebanski ; Schleich ; LeoPhil  in
classical electromagnetism: there, Maxwell’s equations are the equations of
electromagnetic fields in impedance–matched magneto–electric media Plebanski .
In our model of quantum electromagnetism in media, we have reproduced an
intriguing feature of renormalization in curved space–time, the van der Waals
anomaly (Sec. VC). This dielectric analogue of the trace anomaly Wald  of
quantum fields in curved space BD  has been derived before in spherical and
planar media Itai , here we have deduced it in the most elementary way. In
general relativity, the trace anomaly appears in Einstein’s field equations
LL2  as the cosmological constant LondonProceedings , it has the correct order
of magnitude and is consistent with astronomical data Berechya , so perhaps it
simply is the cosmological constant Annals ; Tkatchenko .

Our simple model has come a long way. It has shown to contain rich physics
ranging from liquid paint to the cosmological constant. The model is simple
enough for elementary calculations, but sufficiently sophisticated for deep
cross–field connections. We hope to have made it “as simple as possible, but
not simpler”.

## Acknowledgements

We thank Yael Avni, Gregory Falkovich, Jonathan Kogman, Ziv Landau, Alexander
Poddubny, and René Sedmik for discussions. L.M.R. was supported by the
Austrian Science Fund (FWF) through Project No. P 32300-N27 (WaveLand). U.L.
was supported by the Israel Science Foundation and the Murray B. Koffler
Professorial Chair. For calculating the quantum forces between large numbers
of particles we used the Vienna Scientific Cluster (VSC).

##  Appendix A Numerical procedures

Our method of using the transfer matrix for calculating Casimir forces on
individual scatterers is efficient and fast, but two problems arise for large
wavenumbers $\kappa$ (becoming relevant when the number $N$ of scatterers is
large): the products of the exponential propagation matrices (13) in the
recurrence relation (14) may cause overflow and the required $\kappa$
integrations in the force formula (9) may be slow or insufficiently accurate.
We recommend the remedies described here. We also indicate how the
calculations can be done in multiple parallel threads.

But first we make the following simplification. We do not need the full
transfer matrix and its derivative, but only the elements $T_{22}$ and
$\partial_{j}T_{22}$ for calculating the force [Eq. (9)]. Equations (16) and
(25) and the recurrence relations (14) imply that both $T_{22}$ and
$\partial_{j}T_{22}$ depend only on the second column of A and the second row
of B. We may thus replace the A and B matrices by $N$ two–dimensional vectors:

${\bm{a}}_{j}=\begin{pmatrix}A_{12}\\\
A_{22}\end{pmatrix}\,,\quad{\bm{b}}_{j}=\begin{pmatrix}B_{21}\\\
B_{22}\end{pmatrix}.$ (92)

We may give the ${\bm{a}}$ and ${\bm{b}}$ the initial values

${\bm{a}}_{1}={\bm{b}}_{N}=\begin{pmatrix}0\\\ 1\end{pmatrix}$ (93)

as the $e^{\kappa\delta}$ factors in the initial conditions (15) drop out in
the ratio $T^{\prime}_{22}/T_{22}$. Physically this means that the force is
independent of the propagation to and from the chain of scatterers.

###  A.1 Recurrence

From the macroscopic limit (Sec. IIIA) we know that the ${\bm{a}}$ and
${\bm{b}}$ vectors approach the structure (26) and go with solutions of the
Helmholtz equation (27). Furthermore, from the Madelung representation (38)
with asymptotics (42) follows that ${\bm{a}}$ and ${\bm{b}}$ go as the
exponential of $\int\\!n\kappa\,dx$, which may grow exponentially large,
causing overflow. However, as the force [Eqs. (9) and (12)] depends only on
the ratio between $T^{\prime}_{22}$ and $T_{22}$ we may multiply the
${\bm{a}}$ and ${\bm{b}}$ vectors by some factor, and get the same result. To
avoid overflow, we recommend to employ the inverse of the macroscopic
propagation factor as compensation. This means one should replace the
recurrence relations (14) by

$\displaystyle{\bm{a}}_{j+1}=e^{-n\kappa\delta}\,{\rm P}\,{\rm
R}_{j}\,{\bm{a}}_{j}\,,\quad{\bm{b}}_{j-1}=e^{-n\kappa\delta}\,{\rm P}\,{\rm
R}_{j}^{\top}{\bm{b}}_{j}$ $\displaystyle\mbox{with}\quad n=\sqrt{1+\chi_{j}}$
(94)

where the ${\rm P}$ and ${\rm R}_{j}$ denote the matrices

${\rm P}=\begin{pmatrix}e^{-\kappa\delta}&0\\\
0&e^{\kappa\delta}\end{pmatrix},\quad{\rm
R}_{j}=\mathbb{1}-\frac{\kappa\delta\chi_{j}}{2}\begin{pmatrix}1&1\\\
-1&-1\end{pmatrix}.$ (95)

We obtain from Eqs. (16) and (25) for the $j$–th particle:

$T_{22}={\bm{b}}_{j}\cdot{\rm
R}_{j}\,{\bm{a}}_{j}\,,\quad\partial_{j}T_{22}=-\delta\kappa^{2}\chi_{j}\,{\bm{b}}_{j}\cdot\sigma_{x}{\bm{a}}_{j}$
(96)

in terms of the $\sigma_{x}$ Pauli matrix defined in Eq. (20). Note that for
the calculation of the force on one particle the required transfer–matrix
elements for all other particles are computed. It is wise to cache them such
that they may be used in the force integral for other particles as well.

###  A.2 Integration

For determining the force on the $j$–th particle we need to compute the
integral

$F_{j}=-\frac{\hbar
c}{2\pi}\int_{0}^{\infty}\frac{\partial_{j}T_{22}}{T_{22}}\,d\kappa\,.$ (97)

For numerical integration, we use a modified variant of the Romberg
integration method NumRec  of order $2\times 5$. The Romberg method converges
much faster than many other numerical integration algorithms, and it also has
the advantage that more and more function evaluations in–between the already
evaluated values can be added in each iterative step, until the required
precision has been reached. For integrating from zero to infinity, we apply
the “Mixed Rule” transformation NumRec  $\kappa=\exp[\nu-\exp(-\nu)]$ turning
Eq. (97) into an integral from $-\infty$ to $+\infty$ that needs to be cut
off. Both cutoff values of $\nu$ are determined by an adaptive search
algorithm. Further, to make use of the caching–method mentioned above (Sec.
A1) we have modified the algorithm such that all possible floating–point
values of $\nu$ required for the Romberg integration get unique integer values
assigned. These values serve as lookup index for the cache. For checking the
accuracy of our calculations we used arbitrary–precision arithmetics packages.

###  A.3 Parallelization

In our software implementation, multiple of the above–described routines are
executed in parallel threads. This has to be considered in the caching
mechanism, which not only must be thread–safe, but also must be able to deal
with the possibility that a requested result is in the limbo state of being
currently in calculation by another thread. Multi–threading has allowed us to
run accurate calculations of the van der Waals force on millions of particles.

##  Appendix B Previous renormalization

How is the renormalization procedure of this paper (Sec. V) related to the
previous one Itai ? In the previous procedure Itai , an outgoing Green
function is subtracted from the full Green function as well, but this Green
function is not given as the exact solution of the Green equation (31) with
appropriate conditions at the source, Eqs. (69) and (68). Rather it comes from
the geometrical–optics approximation of the original Green function [Eq. (35)]
in the vicinity of the source where, in the renormalizer, $n(x)$ is taken as
the quadratic expansion

$n(x)=n_{0}+n_{0}^{\prime}(x-x_{0})+\frac{n_{0}^{\prime\prime}}{2}(x-x_{0})^{2}$
(98)

with $n_{0}=n(x_{0})$. We have confirmed [Eqs. (43)] and (64)] that the local
stresses do only depend on maximally the second derivative of $n$, which
justifies the quadratic expansion (98). In our case, this is not an
assumption, but a consequence of the boundary conditions (69) and (68).

The outgoing Green function $g_{0}$ is represented in terms of the optical
length $s$, the amplitude ${\cal A}$ and the first-scattering coefficient
$\beta_{1}$ as

$g_{0}=\mathrm{e}^{-\kappa s}{\cal
A}\left(1+\frac{\beta_{1}s}{\kappa}\right).$ (99)

The optical length is a solution of the eikonal equation Renormalization  (in
one dimension $s^{\prime 2}=n^{2}$) as $(x-x_{0})$ times a term up to
quadratic order. Integrating the eikonal equation we obtain

$s=(x-x_{0})\left[n_{0}+\frac{n_{0}^{\prime}}{2}(x-x_{0})+\frac{n_{0}^{\prime\prime}}{6}(x-x_{0})^{2}\right].$
(100)

The amplitude satisfies the conservation law Renormalization  $({\cal
A}^{2}s^{\prime})^{\prime}=0$ (in one dimension). We get by integration:

${\cal A}=-\frac{1}{2\kappa\sqrt{nn_{0}}}$ (101)

with $n$ given by the quadratic expansion (98). The integration constant gives
the prefactor of the amplitude. It is chosen such that $g_{0}$ satisfies the
Green equation (35) in the immediate vicinity of the delta function where $n$
can be taken as constant $n_{0}$. The first–scattering coefficient
Renormalization  $\beta_{1}$ obeys $2n\beta_{1}={\cal A}^{\prime\prime}/{\cal
A}$ evaluated at $x_{0}$. We express it in terms of the $\beta_{0}$ defined by
Eq. (43) and obtain:

$\beta_{1}=\frac{\beta_{0}}{n_{0}}\,.$ (102)

We substitute these expressions in Eq. (60) for the stresses. Up to
$\kappa^{-1}$ order we obtain exactly the same result [Eq. (64)] as with the
renormalizer advocated in this paper, including the need for the anomaly. The
differences lie in the convergent part, but it is the convergent part of
$\sigma-\sigma_{0}$ that gives the actual physical force. Another difference
is conceptual: the new procedure is much more clear and elegant. However, it
would be very useful to test the theory in experiments, including computer
experiments, bearing in mind the theorem: for every difficult problem there is
a solution that is simple, elegant, and wrong.

##  Appendix C Three layers

In this appendix we consider the case of three homogeneous dielectric layers
(as in Sec. VIA) and derive formulas (76-78) for the Green function.
Expressions for this Green functions are scattered in the literature, with
inconsistent notations, and not all what we need can be found at one place.
For the convenience of the reader we derive them here with an efficient
method, the transfer–matrix method (Sec. IIB).

Consider the transfer from the one homogenous layer with refractive index
$n_{1}$ to another layer with index $n_{2}$. In order to satisfy the Helmholtz
equation (27) the field $\psi$ and its first derivative $\psi^{\prime}$ must
be continuous. From this follows the transfer matrix ${\rm T}_{l}$ across the
interface:

${\rm T}_{l}=\frac{1}{2n_{2}}\begin{pmatrix}n_{2}+n_{1}&n_{2}-n_{1}\\\
n_{2}-n_{1}&n_{2}+n_{1}\end{pmatrix}.$ (103)

In our case, it describes the transfer across the left interface of the
central layer. We will also need the transfer matrix at the right interface:

${\rm T}_{r}=\frac{1}{2n_{3}}\begin{pmatrix}n_{3}+n_{2}&n_{3}-n_{2}\\\
n_{3}-n_{2}&n_{3}+n_{2}\end{pmatrix}.$ (104)

We express the Green function as the product (35) of wave functions
$\psi_{\pm}$ divided by their Wronskian (30). The $\psi_{\pm}$ we determine
with transfer matrices.

Consider first the case where both the point of emission $x_{0}$ and the
observation point $x$ lie in the central layer. There we make the ansatz
$\psi_{+}=e^{n_{2}\kappa x}+\varrho\,e^{-n_{2}\kappa x}$. In the left layer
$\psi_{+}\propto e^{n_{1}\kappa x}$. Consequently,

${\rm T}_{l}\begin{pmatrix}0\\\ 1\end{pmatrix}\propto\begin{pmatrix}\varrho\\\
1\end{pmatrix},$ (105)

which gives for $\varrho$ the Fresnel coefficient $\varrho_{l}$ of Eq. (77).
For $\psi_{-}$ we make the ansatz $\psi_{-}=e^{-n_{2}\kappa
x}+\varrho\,e^{n_{2}\kappa x}$ in the central layer and require
$\psi_{-}\propto e^{-n_{3}\kappa x}$ after propagating to the right layer and
crossing the interface:

${\rm T}_{r}\begin{pmatrix}e^{-n_{2}\kappa a}&0\\\ 0&e^{n_{2}\kappa
a}\end{pmatrix}\begin{pmatrix}1\\\
\varrho\end{pmatrix}\propto\begin{pmatrix}1\\\ 0\end{pmatrix},$ (106)

which gives $\varrho=\varrho_{r}\,e^{-2n_{2}\kappa a}$ in terms of the other
Fresnel coefficient (77). We thus have in the central layer:

$\displaystyle\psi_{+}$ $\displaystyle=$ $\displaystyle e^{n_{2}\kappa
x}+\varrho_{l}\,e^{-n_{2}\kappa x}\,,$ $\displaystyle\psi_{-}$
$\displaystyle=$ $\displaystyle e^{-n_{2}\kappa
x}+\varrho_{r}\,e^{n_{2}\kappa(x-2a)}\,.$ (107)

The Wronskian (30) we read off from Eq. (107) as

$W=-2n_{2}\kappa\left(1-\varrho_{l}\varrho_{r}\,e^{-2n_{2}\kappa a}\right).$
(108)

Expressions (107) and (108) give the Green function (35). Note that one can
interpret this $g(x,x_{0})$ as the summed–up geometrical series of all
multiple reflections Essential . In Sec. VIA we need
$g(x,x)=W^{-1}\psi_{+}(x)\psi_{-}(x)$ and get Eq. (76) from the expressions
derived here.

Now turn to the left layer. There $\psi_{+}$ naturally propagates without
scattering as $e^{n_{1}\kappa x}$ while $\psi_{-}$ is a superposition of the
$e^{\pm n_{1}\kappa x}$. Therefore we write

$\psi_{+}=e^{n_{1}\kappa x}\,,\quad\psi_{-}=e^{-n_{1}\kappa
x}+\varrho\,e^{n_{1}\kappa x}$ (109)

and have the Wronskian

$W=-2n_{1}\kappa\,.$ (110)

We require that $\psi_{-}\propto e^{-n_{3}\kappa x}$ after arriving at the
right layer:

${\rm T}_{r}\begin{pmatrix}e^{-n_{2}\kappa a}&0\\\ 0&e^{n_{2}\kappa
a}\end{pmatrix}{\rm T}_{l}\begin{pmatrix}1\\\
\varrho\end{pmatrix}\propto\begin{pmatrix}1\\\ 0\end{pmatrix},$ (111)

solve for $\varrho$ and express the result in terms of the Fresnel
coefficients (77):

$\varrho=-\frac{\varrho_{l}-\varrho_{r}e^{-2n_{2}\kappa
a}}{1-\varrho_{l}\varrho_{r}\,e^{-2n_{2}\kappa a}}\,.$ (112)

Inserting the result (112) in expressions (109) for the waves $\psi_{\pm}$
with Wronskian (110) we obtain the desired Green function
$W^{-1}\psi_{+}(x)\psi_{-}(x)$ and arrive at formula (78).

The coefficient (112) describes the reflection at the entire central layer.
The physical interpretation of $\varrho$ becomes apparent when we
mathematically transform formula (112) to

$\varrho=-\varrho_{l}+(1-\varrho_{l}^{2})\,\varrho_{r}\,\frac{e^{-2n_{2}\kappa
a}}{1-\varrho_{l}\varrho_{r}\,e^{-2n_{2}\kappa a}}$ (113)

and employ the geometric series

$\frac{1}{1-\varrho_{l}\varrho_{r}\,e^{-2n_{2}\kappa
a}}=\sum_{m=0}^{\infty}\left(\varrho_{l}\varrho_{r}\,e^{-2n_{2}\kappa
a}\right)^{m}\,.$ (114)

Equations (113) and (114) describes the following scenario: a wave comes in
from the left. At the left interface a part of the wave is reflected with
reflectivity $-\varrho_{l}$ (the minus sign because $\varrho_{l}$ describes
reflection from the right to the left, not from the left to the right). The
other part is transmitted with transmittance $\sqrt{1-\varrho_{l}^{2}}$. The
transmitted part propagates to the right interface with factor
$e^{-n_{2}\kappa a}$. If the wave makes it back (such that it counts in the
overall reflection coefficient $\varrho$) it must get reflected at the right
interface with coefficient $\varrho_{r}$ and propagate back to the left
interface with factor $e^{-n_{2}\kappa a}$ and be transmitted out of the
central layer with transmittance $\sqrt{1-\varrho_{l}^{2}}$. This gives the
factor $(1-\varrho_{l}^{2})\,\varrho_{r}\,e^{-2n_{2}\kappa a}$ in Eq. (113).
The wave may also get multiply reflected in the central layer, bouncing back
and forth multiple times, as the series (114) describes.

,

## References

  * (1) P. Milonni, The Quantum Vacuum: An Introduction to Quantum Electrodynamics (Academic Press, Boston, 1994).
  * (2) K. A. Milton, The Casmir Effect (World Scientific, Singapore, 2001).
  * (3) M. Bordag, G. L. Klimchitskaya, U. Mohideen, and V. M. Mostepanenko, Advances in the Casimir Effect (Oxford University Press, Oxford, 2009).
  * (4) A. W. Rodriguez, F. Capasso, and S. G. Johnson, The Casimir effect in microstructured geometries, Nat. Photon. 5, 211 (2011).
  * (5) S. Y. Buhmann, Dispersion Forces (Springer, Berlin, 2012).
  * (6) W. M. R. Simpson and U. Leonhardt (eds.) Forces of the quantum vacuum (World Scientific, Singapore, 2015).
  * (7) K. A. Milton (ed.) The State of the Quantum Vacuum (World Scientific, Singapore, 2023).
  * (8) J. Schwinger, L. L. DeRaad, and K. A. Milton, Casimir Effect in Dielectrics, Ann. Phys. (New York) 115, 1 (1978).
  * (9) S. K. Lamoreaux, Casimir forces: Still surprising after 60 years, Phys. Today 60, 40 (2007).
  * (10) J. N. Israelachvili, Intermolecular and Surface Forces, (Elsevier, Amsterdam, 2011).
  * (11) J. Hermann, R. A. DiStasio, and A. Tkatchenko, First-Principles Models for van der Waals Interactions in Molecules and Materials: Concepts, Theory, and Applications, Chem. Rev. 117, 4714 (2017).
  * (12) F. London, Zur Theorie und Systematik der Molekularkräfte, Z. Phys. 63, 245 (1930).
  * (13) H. B. G. Casimir and D. Polder, The Influence of Retardation on the London–van der Waals Forces, Phys. Rev. 73, 360 (1948).
  * (14) E. Shahmoon, Van der Waals and Casimir Polder dispersion forces, in Ref. Forces, .
  * (15) P. W. Milonni, The Casimir Effect: Physical Manifestations of Zero-Point Energy, Phys. Today 56, 49 (2003).
  * (16) H. B. G. Casimir, On the attraction between two perfectly conducting plates, Koninkl. Ned. Akad. Wetenschap. 51, 793 (1948).
  * (17) E. M. Lifshitz, The Theory of Molecular Attractive Forces between Solids, Sov. Phys. JETP 2, 73 (1956).
  * (18) I. E. Dzyaloshinskii, E. M. Lifshitz, and L. P. Pitaevskii, The general theory of van der Waals forces, Adv. Phys. 10, 165 (1961).
  * (19) A. Milling, P. Mulvaney, and I. Larson, Direct Measurement of Repulsive van der Waals Interactions Using an Atomic Force Microscope, J. Colloid Interface Sci. 180, 460 (1996).
  * (20) A. Meurk, P. F. Luckham, and L. Bergström, Direct Measurement of Repulsive and Attractive van der Waals Forces between Inorganic Materials, Langmuir 13, 3896 (1997).
  * (21) S.-W. Lee and W. M. Siegmund, Repulsive van der Waals Forces for Silica and Alumina, J. Colloid Interface Sci. 243, 365 (2001).
  * (22) J. N. Munday, F. Capasso, and V. A. Parsegian, Measured long–range repulsive Casimir–Lifshitz forces, Nature 457, 170 (2009).
  * (23) 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).
  * (24) I. E. Dzyaloshinskii and L. P. Pitaevskii, Van der Waals Forces in an Inhomogeneous Dielectric, Sov. Phys. JETP 9, 1282 (1959).
  * (25) W. M. R. Simpson, Surprises in Theoretical Casimir Physics (Springer, Berlin, 2014).
  * (26) 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).
  * (27) I. Griniasty and U. Leonhardt, Casimir stress inside planar materials, Phys. Rev. A 96, 032123 (2017) and references cited therein.
  * (28) I. Griniasty and U. Leonhardt, Casimir stress in materials: Hard divergency at soft walls, Phys. Rev. B 96, 205418 (2017).
  * (29) P. Parashar, K. A. Milton, Y. Li, H. Day, X. Guo, S. A. Fulling, and I. Cavero-Peláez, Quantum electromagnetic stress tensor in an inhomogeneous medium, Phys. Rev. D 97, 125009 (2018).
  * (30) Y. Li, K. A. Milton, X. Guo, G. Kennedy, and S. A. Fulling, Casimir forces in inhomogeneous media: Renormalization and the principle of virtual work, Phys. Rev. D 99, 125004 (2019).
  * (31) I. Y. Efrat and U. Leonhardt, Van der Waals anomaly: Analog of dark energy with ultracold atoms, Phys. Rev. B 104, 235432 (2021).
  * (32) A. Shayit, S. A. Fulling, T. E. Settlemyre, and J. Merritt, Vacuum energy density and pressure inside a soft wall, Int. J. Mod. Phys. A 37, 2241007 (2022).
  * (33) S. A. Franchino-Vinas, M. N. Mantinan, and F. D. Mazzitelli, Quantum vacuum fluctuations and the principle of virtual work in inhomogeneous backgrounds, Phys. Rev. D 105, 085023 (2022).
  * (34) This is different from planar media that are effectively one–dimensional, see e.g. J. Fiedler, F. Spallek, P. Thiyam, C. Persson, M. Boström, M. Walter, and S. Y. Buhmann, Phys. Rev. A 99, 062512 (2019).
  * (35) E. Shahmoon, I. Mazets, and G. Kurizki, Giant vacuum forces via transmission lines, Proc. Natl. Acad. Sci. USA 111, 10485 (2014).
  * (36) A. Lambrecht, P. A. Maia Neto, and S. Reynaud, The Casimir effect within scattering theory, New J. Phys. 8, 243 (2006).
  * (37) A. Wirzba, The Casimir effect as a scattering problem, J. Phys. A: Math. Theor. 41, 164003 (2008).
  * (38) G.-L. Ingold and A. Lambrecht, Casimir effect from a scattering approach, Am. J. Phys. 83, 156 (2015).
  * (39) L. P. Pitaevskii, On the Problem of van der Waals Forces in Dielectric Media in D. Dalvit, P. Milonni, D. Roberts, and F. Rosa (eds.) Casimir Physics (Springer, Berlin, 2011).
  * (40) See e.g. I. Schlesinger and U. Sivan, New Information on the Hydrophobic Interaction Revealed by Frequency Modulation AFM, Langmuir 33, 2485 (2017).
  * (41) G. L. Klimchitskaya and V. M. Mostepanenko, Current Status of the Problem of Thermal Casimir Force, Int. J. Mod. Phys. 37, 241002 (2022).
  * (42) N. D. Birrell and P. C. W. Davies. Quantum fields in curved space (Cambridge University Press, Cambridge, 1984).
  * (43) L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields (Butterworth-Heinemann, Amsterdam, 2003).
  * (44) W. Gordon, Zur Lichtfortpflanzung nach der Relativitätstheorie, Ann. Phys. (Leipzig) 72, 421 (1923).
  * (45) J. Plebanski, Electromagnetic Waves in Gravitational Fields, Phys. Rev. 118, 1396 (1960).
  * (46) 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).
  * (47) U. Leonhardt and T. G. Philbin, Geometry and Light: the Science of Invisibility, (Dover, Mineola, 2010).
  * (48) R. M. Wald, Trace anomaly of a conformally invariant quantum field in curved spacetime, Phys. Rev. D 17, 1477 (1978).
  * (49) U. Leonhardt, Lifshitz theory of the cosmological constant, Ann. Phys. (New York) 411, 167973 (2019).
  * (50) For an alternative, see A. Tkatchenko and D. V. Fedorov, Casimir Self-Interaction Energy Density of Quantum Electrodynamic Fields, Phys. Rev. Lett. 130, 041601 (2023).
  * (51) J. Friedel, Electronic structure of primary solid solutions in metals, Adv. Phys. 3, 446 (1954).
  * (52) M. Sh. Birman and M. G. Krein, On the theory of wave operators and scattering operators, Dokl. Akad. Nauk SSSR 144, 475 (1962).
  * (53) P. Lloyd, Wave propagation through an assembly of spheres: II. The density of single-particle eigenstates, Proc. Phys. Soc. London 90, 207 (1967).
  * (54) J. S. Faulkner, Scattering theory and cluster calculations, J. Phys. C: Solid State Phys. 10, 4661 (1977).
  * (55) L. L. Sanchez-Soto, J. J. Monzon, A. G. Barriuso, and J. F. Carinena, The transfer matrix: A geometrical perspective, Phys. Rep. 513, 191 (2012).
  * (56) E. Madelung, Eine anschauliche Deutung der Gleichung von Schrödinger, Naturwiss. 14, 1004 (1926).
  * (57) D. Bohm, A Suggested Interpretation of the Quantum Theory in Terms of ”Hidden” Variables. II, Phys. Rev. 85, 180 (1952).
  * (58) M. J. Ablowitz and A. S. Fokas, Complex Variables: Introduction and Applications (Cambridge University Press, Cambridge, 1997).
  * (59) L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media (Pergamon, Oxford, 1984).
  * (60) E. Poisson, A Relativist’s Toolkit (Cambridge University Press, Cambridge, 2009).
  * (61) A. Ashkin, Acceleration and Trapping of Particles by Radiation Pressure, Phys. Rev. Lett. 24, 156 (1970).
  * (62) A. Ashkin and J. Dziedzic, Optical levitation by radiation pressure, Appl. Phys. Lett. 18, 283 (1971).
  * (63) A. Ashkin, Optical trapping and manipulation of neutral particles using lasers, Proc. Natl. Acad. Sci. USA 94, 4853 (1997).
  * (64) K. C. Neuman and S. M. Block, Optical trapping, Rev. Sci. Instr. 75, 2787 (2004).
  * (65) G. Pesce, P. H. Jones, O. M. Marago, and G. Volpe, Optical tweezers: theory and practice, Eur. Phys. J. Plus 135, 949 (2020).
  * (66) S. Scheel, The Casimir stress in real materials, in Ref. Forces .
  * (67) Note that this simple substitution relies on the fact that $K+i\Gamma$ is analytic and decaying on the lower half plane of complex time difference, which is only true for the vacuum state. In the thermal state $K$ is periodic (Kubo–Martin–Schwinger relation Annals ; LL7 ).
  * (68) L. P. Pitaevskii, Comment on “Casimir force acting on magnetodielectric bodies embedded in media”, Phys. Rev. A 73, 047801 (2006).
  * (69) 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).
  * (70) F. A. Burger, R. W. Corkery, S. Y. Buhmann, and J. Fiedler, Comparison of Theory and Experiments on van der Waals Forces in Media – A Survey, J. Phys. Chem. C 124, 24179 (2020).
  * (71) L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Pergamon, Oxford, 1977).
  * (72) L. D. Landau and E. M. Lifshitz, Statistical Physics, Part 2 (Pergamon, Oxford, 1980).
  * (73) U. Leonhardt, The case for a Casimir cosmology, Phil. Trans. R. Soc. A 378, 20190229 (2020).
  * (74) U. Leonhardt, Essential Quantum Optics: From Quantum Measurements to Black Holes, Appendix A (Cambridge University Press, Cambridge, 2010).
  * (75) A. Erdélyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi, Higher Transcendental Functions (McGraw-Hill, New York, 1981).
  * (76) L. Knöll, S. Scheel, and D.-G. Welsch, QED in dispersing and absorbing media, in J. Peřina (ed.) Coherence and Statistics of Photons and Atoms (Wiley, New York, 2001); arXiv:quant-ph/0006121.
  * (77) D. Berechya and U. Leonhardt, Lifshitz cosmology: quantum vacuum and Hubble tension, Month. Not. Roy. Astron. Soc. 507, 3473 (2021).
  * (78) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge, 2007).