Computation and visualization of Casimir forces in arbitrary geometries:
  non-monotonic lateral forces and failure of proximity-force approximations
Computation and visualization of Casimir forces in arbitrary geometries:
non-monotonic lateral-wall forces and failure of proximity-force approximations
Alejandro Rodriguez,1 Mihai Ibanescu,1 Davide Iannuzzi,2
Federico Capasso,3 J. D. Joannopoulos,1 and Steven G. Johnson1
1Center for Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139
2Faculty of Sciences, Department of Physics and Astronomy, Vrije Universiteit Amsterdam, The Netherlands
3Department of Physics, Harvard University, Cambridge, MA 02139
We present a method of computing Casimir forces for arbitrary geometries, with any desired ac-
curacy, that can directly exploit the efficiency of standard numerical-electromagnetism techniques.
Using the simplest possible finite-difference implementation of this approach, we obtain both agree-
ment with past results for cylinder-plate geometries, and also present results for new geometries.
In particular, we examine a piston-like problem involving two dielectric and metallic squares sliding
between two metallic walls, in two and three dimensions, respectively, and demonstrate non-additive
and non-monotonic changes in the force due to these lateral walls.
PACS numbers:
Casimir forces arise between macroscopic objects due
to changes in the zero-point energy associated with quan-
tum fluctuations of the electromagnetic field [1]. This
spectacular effect has been subject to many experimental
validations, as reviewed in Ref. 2. All of the experiments
reported so far have been based on simple geometries
(parallel plates, crossed cylinders, or spheres and plates).
For more complex geometries, calculations become ex-
tremely cumbersome and often require drastic approx-
imations, a limitation that has hampered experimental
and theoretical work beyond the standard geometries.
In this letter, we present a method to compute Casimir
forces in arbitrary geometries and materials, with no un-
controlled approximations, that can exploit the efficient
solution of well-studied problems in classical computa-
tional electromagnetism. Using this method, which we
first test for geometries with known solutions, we predict
a non-monotonic change in the force arising from lateral
side walls in a less-familiar piston-like geometry (Fig. 2).
Such a lateral-wall force cannot be predicted by “addi-
tive” methods based on proximity-force or other purely
two-body–interaction approximations, due to symmetry,
and it is difficult to find a simple correction to give a
non-monotonic force. We are able to compute forces for
both perfect metals and arbitrary dispersive dielectrics,
and we also obtain a visual map of the stress tensor that
directly depicts the interaction forces between objects.
The Casimir force was originally predicted for parallel
metal plates, and the theory was subsequently extended
to straighforward formulas for any planar-multilayer di-
electric distribution ε(x, ω) via the generalized Lifshitz
formula [3]. In order to handle more arbitrary geome-
tries, two avenues have been pursued. First, one can em-
ploy approximations derived from limits such as that of
parallel plates; these methods include the proximity-force
approximation (PFA) and its refinements [4], renormal-
ized Casimir-Polder [5] or semi-classical interactions [6],
multiple-scattering expansions [7], classical ray optics [8],
and various perturbative techniques [9, 10]. Such meth-
ods, however, involve uncontrolled approximations when
applied to arbitrary geometries outside their range of ap-
plicability, and have even been observed to give quali-
tatively incorrect results [11, 12]. Therefore, researchers
have instead sought numerical methods applicable to ar-
bitrary geometries that converge to the exact result given
sufficient computational resources. One such method
uses a path-integral representation for the effective ac-
tion [13], and has predicted the force between a cylin-
der and a plane or between corrugated surfaces. Ref. 13
uses a surface parameterization of the fields coupled via
vacuum Green’s functions, requiring O(N2) storage and
O(N3) time for N degrees of freedom, making scaling
to three dimensions (3d) problematic. Another exact
method is the “world-line approach” [12, 14, 15], based
on Monte-Carlo path-integral calculations. (The scaling
of the world-line method involves a statistical analysis,
determined by the relative feature sizes in the geome-
try, that is beyond the scope of this Letter.) Further-
more, the methods of Ref. 13 and Ref. 14 have currently
only been demonstrated for perfect-metallic z-invariant
structures—in this case, the vector unknowns can be de-
composed into TE (E · ẑ = 0) and TM (H · ẑ = 0) scalar
fields with Dirichlet or Neumann boundary conditions—
although generalizations have been proposed [16]. Here,
we propose a method based on evaluation of the mean
stress tensor via the fluctuation-dissipation theorem,
which only involves repeated evaluation of the electro-
magnetic imaginary-frequency Green’s function. For a
volume discretization with N degrees of freedom and
an efficient iterative solver, this requires O(N) stor-
age and O(N2−1/d) time in d dimensions. Further-
more, because evaluation of the Green’s function is such
a standard problem in classical computational electro-
magnetism, it will be possible to exploit many develop-
ments in fast solvers, such as finite-element, or boundary-
element methods [17]. To illustrate the method, our
initial implementation is based on the simplest-possible
finite-difference frequency-domain (FDFD) method, as
described below.
As derived by Dzyaloshinskĭı et al. [1], the net Casimir
force on a body can be expressed as an integral over any
closed surface around the body of the mean electromag-
netic stress tensor 〈Tij〉, integrated over imaginary fre-
quencies ω = iw:
surface
〈Tij(r, iw)〉 dSj . (1)
For a 3d z-invariant structure, the z integral is replaced
by an integral over the corresponding wavevector, result-
ing in a net force per unit length. The stress tensor is
defined as usual by:
〈Tij(r, iw)〉 = 〈Hi(r)Hj(r)〉 −
〈Hk(r)Hk(r)〉
+ ε(r, iw)
〈Ei(r)Ej(r)〉 −
〈Ek(r)Ek(r)〉
The connection to quantum mechanics arises from the
correlation functions of the fluctuating fields, such as
〈EiEj〉, given via the fluctuation-dissipation theorem in
terms of the imaginary-ω Green’s function Gij(iw; r−r′):
〈Ei(r)Ej(r′)〉 =
Gij(iw; r− r′) (3)
〈Hi(r)Hj(r′)〉 = −(∇×)i`(∇′×)jmG`m(iw; r− r′) , (4)
where the Green’s function Gij solves the equation:[
∇×∇×+
ε(r, iw)
Gj(iw; r−r′) = êjδ(r−r′) (5)
for a unit vector êj in the j direction, and obeys the usual
boundary conditions on the electric field from classical
electromagnetism. (The above expressions are at zero
temperature; the nonzero-temperature force is found by
changing
dw in Eq. 1 into a discrete summation [1].)
Although the Green’s function (and thus Tij) is formally
infinite at r = r′, this divergence is conventionally re-
moved by subtracting the vacuum Green’s function; in
a numerical method with discretized space, as below,
there is no divergence and no additional regularization
is required. (The vacuum Green’s function gives zero net
contribution to the dS integral, and therefore need not
be removed as long as the integrand is finite.)
Historically, this stress-tensor expression was used to
derive the standard Lifshitz formula for parallel plates,
where Gij is known analytically. However, it also forms
an ideal starting point for a computational method, be-
cause the Green’s function for arbitrary geometries is
routinely computed numerically by a variety of tech-
niques [17]. Furthermore, the problem actually becomes
easier for an imaginary ω. First, for an imaginary ω, the
linear operator in Eq. 5 is real-symmetric and positive-
definite for w 6= 0, since the dielectric function ε(ω) is
purely real and positive along the imaginary-ω axis for
physical materials without gain, due to causality. Second,
the imaginary-ω Green’s function is exponentially decay-
ing rather than oscillating, leading to a well-behaved non-
oscillatory integrand in Eq. 1.
To illustrate this method, we employed the simplest
possible computational technique: we perform a FDFD
discretization of Eq. 5 with a staggered Yee grid [18]
and periodic boundaries, inverting the linear operator by
a conjugate-gradient method. The presence of discon-
tinuous material interfaces degrades second-order finite-
difference methods to only first-order accuracy, and the
uniform spatial resolution is also suboptimal, but we
found FDFD to be nevertheless adequate for small 2d
geometries. The periodicity leads to artificial “wrap-
around” forces that decay rapidly with cell size L (at least
as 1/L3 in 2d and 1/L4 in 3d); we chose cell sizes large
enough to make these contributions negligible (< 1%).
The computational process is as follows: pick some sur-
face/contour around a given body, evaluate the Green’s
function for every grid point on this surface in order to
compute the surface integral of the stress tensor, which
is then integrated over w by adaptive quadrature.
1  10 
FIG. 1: Casimir force between a radius-R cylinder and a
plate (inset), relative to the proximity-force approximation
FPFA, vs. normalized separation a/R. The solid lines are
the Casimir force computed in Ref. 19 for TE (gray) and
TM (blue) polarizations, along with results computed by our
method with a simple finite-difference discretization (gray
squares). Error bars were estimated for some data points
by using computations at multiple spatial resolutions. Inset
shows interaction stress tensor ∆〈Txx〉 at a typical imaginary
frequency w = 2πc/a, where red indicates attractive stress.
Before we attempt to study new geometries with our
method, it is important to check it against known re-
sults. The simplest cases, of parallel metallic or dielec-
tric plates, of course match the known result from the
Liftshitz formula and are not reproduced here. A more
complicated geometry, consisting of a perfect metallic
cylinder adjacent to a perfect metallic plate in 3d, was
solved numerically by Ref. 19, to which our results are
compared in Fig. 1. Ref. 19 used a specialized Fourier-
Bessel basis specific to this cylindrical geometry, which
should have exponential (spectral) convergence. Our use
of a simple uniform grid was necessarily much less effi-
cient, especially with the first-order accuracy, but was
able to match the Ref. 19 results within ∼ 3% using rea-
sonable computational resources. A simple grid has the
advantage of being very general, as illustrated below, but
other general bases with much greater efficiency are pos-
sible using finite-element or boundary-element methods;
the latter, in particular, could use a spectral Fourier ba-
sis analogous to Ref. 19 and exploit a fast-multipole or
similar O(N logN) solver technique [17].
Also shown, in the inset of Fig. 1, is a plot of the
interaction stress-tensor component ∆〈Txx〉 at a typical
imaginary frequency w = 2πc/a. By “interaction” stress-
tensor ∆〈Tij〉, we mean the total 〈Tij〉 of the full geome-
try minus the sum of the 〈Tij〉’s computed for each body
in isolation. Here, the stress tensors of the isolated cylin-
der and plate have been subtracted, giving us a way to
visualize the force due to the interaction. As described
below, such stress plots reveal the regions in which two
objects most strongly affect one another, and therefore
reveal where a change of the geometry would have the
most impact. (In contrast, Ref. 12 plots an interaction-
energy density that does not directly reveal the force,
since the force requires the energy to be differentiated
with respect to a. For example, Ref. 12’s subtracted en-
ergy density apparently goes nearly to zero as a metallic
surface is approached, whereas the stress tensor cannot
since the stress integration surface is arbitrary.)
We now consider a more complicated geometry in
which there are interactions between multiple bodies: a
3d “piston”-like structure, shown in the inset of Fig. 2,
consisting of two z-invariant metal s × s squares sepa-
rated by a distance a from one another (here, s = a)
and separated by a distance h from infinite metal plates
on either side. We then compute the Casimir force per
unit z between the two squares as a function of the sep-
aration h. The result for perfect conductors is shown
in Fig. 2, plotted for the TE and TM polarizations and
also showing the total force. (Error bars are not shown
because the estimated error is < 1%.) In the limit of
h → 0, this structure approaches the “Casimir piston,”
which has been solved analytically [20, 21] (and also in
2d for the TM polarization [22]). Our results, extrapo-
lated to h = 0, agree with these results to within ≈ 2%
(although we have computational difficulties for small h
due to the high resolution required to resolve a small
0  0.2 0.4 0.6 0.8 1
total
FIG. 2: Casimir force per unit length between metal squares
F/FPFA, vs. distance from metal plate h (inset), normalized
by the total TE+TM force per unit length obtained using the
PFA, FPFA = ~csζ(3)/480πa4. The total force is plotted
(black squares) along with the TE (red dots) and TM (blue
circles) contributions.
feature in FDFD). For h > 0, however, the result is sur-
prising in at least two ways. First, the total force is
non-monotonic in h, due to a competition between the
TE and TM contributions to the forces. Second, the h
dependence of the force is a lateral effect of the parallel
plates on the squares, which would be zero by symmetry
in PFA or any other two-body–interaction approxima-
tion. Although lateral-wall effects can clearly arise qual-
itatively in various approximations, such as in ray optics
or in PFA restricted to “line-of-sight” interactions, non-
monotonicity is more surprising[24]. Also, in the large-h
limit, the force remains different from PFA due to finite-s
“edge” effects [12], which are captured by our method.
Our method is also capable, without modification, of
handling dielectric materials. This is demonstrated in
Fig. 3, where the Casimir force is shown for the case
where the squares are made of gold, using the experimen-
tal Drude ε(ω) from Ref. 23 for a separation a = 1µm.
Here, our calculation is for a purely 2d geometry (equiv-
alently, 3d restricted to z-invariant fields/currents), and
for comparison we also plot the corresponding 2d force
for perfect-metal squares (although the two cases are
normalized differently as described in the caption). As
might be expected, the dielectric squares have a weaker
interaction than the perfect-metal squares, but are still
non-monotonic. Note also the qualitative similarity be-
tween the perfect-metal results of Figs. 2 and 3, reflecting
the fact that the force contributions are dominated by
the zero-wavevector (z-invariant) fields. Extrapolated to
h = 0, the perfect-metal TM force agrees with the known
analytical result [22] to within ≈ 3%.
0 0.2 0.4 0.6 0.8 1 1.2
sΑu a x
FIG. 3: Solid lines: Casimir force between 2d gold squares
F/FPFA, vs. distance from metal plate h (inset), using ex-
perimental ε(ω) [23], normalized by the total force obtained
using the PFA. (Here, the PFA force is computed for x-infinite
gold slabs). The total force is plotted (black squares) along
with the TE (red dots) and TM (blue circles) contributions.
Dashed lines: force for 2d perfect-metal squares, normalized
by the perfect-metal PFA force FPFA = ~csζ(3)/8πa3.
To further explore the source of the h-dependence, we
plot the TM interaction-stress maps ∆〈Txx〉 and ∆〈Txy〉
in Fig. 4, for the 2d perfect-metal squares from Fig. 3.
The stress plots of Fig. 4 are computed at a typical fre-
quency w = 2πc/a, and for varying distances from the
metal plates (h = 0.5, 1.0, 2.0). As shown, the mag-
nitudes of both the xx (a–c) and xy (d–f) components
of the stress tensor change dramatically as the metal
plates are brought closer to the squares. For example,
one change in the force integral comes from Txy, which
for isolated squares has an asymmetric pattern at the
four corners that will contribute to the attractive force,
whereas the presence of the plates induces a more sym-
metric pattern of stresses at the four corners that will
have nearly zero integral. This results in a decreasing
TM force with decreasing h. Because stress maps indi-
cate where bodies interact and with what signs, it may
be useful in future work to explore whether they can be
used to design unusual behaviors such as non-additive,
non-monotonic, or even repulsive forces.
This work was supported in part by the Nanoscale Sci-
ence and Engineering Center (NSEC) under NSF con-
tract PHY-0117795, by the Materials Research Science
and Engineering Center program of the NSF under award
DMR-9400334, and by a DOE Comp. Science Grad. Fel-
lowship under grant DE–FG02-97ER25308. D. I. grate-
fully acknowledges support from the Netherlands Organ-
isation for Scientific Research (NWO), under the IRI
Scheme Vernieuwingsimpuls VIDI-680-47-209.
Txx Txy
FIG. 4: (a–f): TM stress map of the 2d-analogous geome-
try of Fig. 2 for various h. The intearaction stress tensors
〈Txx〉 (left) and 〈Txy〉 (right) for: (a),(d): h = 0.5a; (b),(e):
h = a; and (c),(f): h = 2a, where blue/white/red = repul-
sive/zero/attractive.
[1] E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics:
Part 2 (Pergamon, Oxford, 1980).
[2] R. Onofrio, New J. Phys. 8, 237 (2006).
[3] M. S. Tomaš, Phys. Rev. A 66, 052103 (2002).
[4] M. Bordag, Phys. Rev. D 73, 125018 (2006).
[5] M. Tajmar, Intl. J. Mod. Phys. C 15, 1387 (2004).
[6] M. Schaden and L. Spruch, Phys. Rev. A 58, 935 (1998).
[7] R. Balian and B. Duplantier, Ann. Phys. 112, 165 (1978).
[8] R. L. Jaffe and A. Scardicchio, Phys. Rev. Lett. 92,
070402 (2004).
[9] T. Emig, A. Hanke, R. Golestanian, and M. Kardar,
Phys. Rev. Lett. 87, 260402 (2001).
[10] R. B. Rodrigues, P. A. M. Neto, A. Lambrecht, and
S. Reynaud, Phys. Rev. Lett. 96, 100402 (2006).
[11] O. Kenneth, I. Klich, A. Mann, and M. Revzen, Phys.
Rev. Lett. 89, 033001 (2002).
[12] H. Gies and K. Klingmuller, Phys. Rev. Lett. 97, 220405
(2006).
[13] R. Büscher and T. Emig, Phys. Rev. Lett. 94, 133901
(2005).
[14] H. Gies, K. Langfeld, and L. Moyaerts, J. High Energy
Phys. p. 018 (2003).
[15] H. Gies and K. Klingmuller, Phys. Rev. Lett. 96, 220401
(2006).
[16] T. Emig and R. Büscher, Nucl. Phys. B 696, 468 (2004).
[17] W. C. Chew, J. Jian-Ming, E. Michielssen, and S. Jiming,
Fast and Efficient Algorithms in Computational Electro-
magnetics (Artech, Norwood, MA, 2001).
[18] A. Christ and H. L. Hartnagel, IEEE Trans. Microwave
Theory Tech. 35, 688 (1987).
[19] T. Emig, R. L. Jaffe, M. Kardar, and A. Scardicchio,
Phys. Rev. Lett. 96, 080403 (2006).
[20] M. P. Hertzberg, R. L. Jaffe, M. Kardar, and A. Scardic-
chio, Phys. Rev. Lett. 95, 250402 (2005).
[21] V. M. Marachevsky, PRD 75, 085019 (2007).
[22] R. M. Cavalcanti, Phys. Rev. D 69, 065015 (2004).
[23] I. Brevik, J. B. Aarseth, J. S. Hoye, and K. A. Milton,
Phys. Rev. E 71 (2005).
[24] A publication with R. L. Jaffe is in preparation that
shows non-monotonicity arising in ray optics.
	References
