Unifying microscopic and continuum treatments of van der Waals and Casimir interactions
Prashanth S. Venkataram Jan Hermann Alexandre Tkatchenko Alejandro W. Rodriguez




Department of Electrical Engineering, Princeton University, Princeton, New Jersey 08544, USA

Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4–6, 14195, Berlin, Germany

Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4–6, 14195, Berlin, Germany
Physics and Materials Science Research Unit, University of Luxembourg, L-1511 Luxembourg

Department of Electrical Engineering, Princeton University, Princeton, New Jersey 08544, USA





  We present an approach for computing long-range van der Waals (vdW)
  interactions between complex molecular systems and arbitrarily
  shaped macroscopic bodies, melding atomistic treatments of
  electronic fluctuations based on density functional theory in the
  former, with continuum descriptions of strongly shape-dependent
  electromagnetic fields in the latter, thus capturing many-body and
  multiple scattering effects to all orders. Such a theory is
  especially important when considering vdW interactions at mesoscopic
  scales, i.e. between molecules and structured surfaces with features
  on the scale of molecular sizes, in which case the finite sizes,
  complex shapes, and resulting nonlocal electronic excitations of
  molecules are strongly influenced by electromagnetic retardation and
  wave effects that depend crucially on the shapes of surrounding
  macroscopic bodies. We show that these effects together can modify
  vdW interactions by orders of magnitude compared to previous
  treatments based on Casimir–Polder or non-retarded approximations,
  which are valid only at macroscopically large or atomic-scale
  separations, respectively.


Unifying microscopic and continuum treatments of van der Waals and Casimir interactions
    Alejandro W. Rodriguez
    December 30, 2023
=======================================================================================




Van der Waals (vdW) interactions play an essential role in
non-covalent phenomena throughout biology, chemistry, and
condensed-matter physics <cit.>. It has long been known that vdW
interactions among a system of polarizable atoms are not
pairwise-additive but instead strongly depend on geometric and
material properties <cit.>. However, only recently developed theoretical methods
have made it possible to account for short-range quantum interactions
in addition to long-range many-body screening in molecular
ensembles <cit.>, demonstrating that nonlocal many-body effects
cannot be captured by simple, pairwise-additive descriptions; these
calculations typically neglect electromagnetic retardation effects in
molecular systems. Simultaneously, recent theoretical and experimental
work has characterized dipolar Casimir–Polder interactions between
macroscopic metallic or dielectric objects and atoms, molecules, or
Bose–Einstein condensates, further extending to nonzero temperatures,
dynamical situations, and fluctuations in excited states (as in
so-called Rydberg atoms) <cit.>. Yet,
while theoretical treatments have thus far accounted for the full
electrodynamic response of macroscopic bodies (including retardation),
they often treat molecules as point dipoles of some effective bulk
permittivities or as collections of noninteracting atomic dipoles,
ignoring finite size and other many-body electromagnetic effects.

In this paper, motivated by the aforementioned theoretical
developments <cit.>, we describe an approach that
seamlessly connects atomistic descriptions of large molecules to
continuum descriptions of arbitrary macroscopic bodies,
characterizing their mutual vdW interactions. In particular, while
molecules that are very close to macroscopic objects require atomistic
descriptions of the latter, and very large molecules that are far from
macroscopic objects require consideration of the contributions of
vibrational (in addition to electronic) resonances to the vdW
interaction energy, we focus on a mesoscopic regime involving
molecular sizes and separations on the order of 1–100 nm,
where macroscopic objects can be treated continuously for the purposes
of computing electromagnetic field responses (and molecular
vibrational resonances can be neglected), yet electromagnetic
retardation in conjunction with the finite sizes, nontrivial shapes,
and nonlocal electronic correlations of large molecules need to be
self-consistently considered to accurately characterize vdW
interactions. We specifically investigate interactions among various
large molecules and gold surfaces, and show that the effect of
nonlocal polarization correlations, encapsulated in the ratio of
retarded, many-body (RMB) to pairwise vdW energies (or forces), causes
relative deviations from pairwise treatments ranging from 20% to over
3 orders of magnitude; further differences of over an order of
magnitude are observed when retardation or finite size effects are
neglected. 


The basis of our work is an equation for the long-range dispersive vdW
energy of a system of polarizable bodies, consisting of N
microscopic bodies (molecules), labeled by k and described by
electric susceptibilities 𝕍_k, and a collection of
continuum bodies (an environment) described by a collective,
macroscopic susceptibility 𝕍_env, displayed
schematically in schematic. The energy of such a collection
of bodies can be obtained from the scattering
framework <cit.> and written as an integral over imaginary
frequency ω = ξ,

    ℰ = ħ/2π∫_0^∞dξ ln[(𝕋_∞𝕋^-1)],

in terms of T-operators that depend on the bodies' susceptibilities as
well as on the homogeneous electric Green's function
𝔾_0(ξ, x⃗, x⃗') = (∇⊗∇ -
ξ^2/c^2𝕀)e^-ξ |x⃗ -
      x⃗'|/c/4π |x⃗ - x⃗'| (including
retardation) mediating electromagnetic interactions; they encode the
scattering properties of the various bodies, and are given by,

    𝕋 = (𝕀 - (𝕍 + 𝕍_env) 𝔾_0)^-1 (𝕍 + 𝕍_env),

where 𝕍 = ∑_k 𝕍_k; 𝕋_∞ =
𝕋_env∏_k𝕋_k, written in terms
of 𝕋_k(env) = (𝕀 -
𝕍_k(env)𝔾_0)^-1𝕍_k(env), encodes the scattering response of the
bodies in isolation from one another <cit.>.



The energy in (<ref>) treats microscopic and macroscopic bodies on
an equal footing, yet the key to its accurate evaluation lies in
appropriately representing the degrees of freedom (DOFs) of each
entity. Typically, macroscopic environments are well described by
continuum susceptibilities 𝕍_env, whose response
can be expanded in a basis of incoming and outgoing propagating
planewaves, as is typical of the scattering
framework <cit.>, or via localized functions,
e.g. tetrahedral mesh elements, in brute-force
formulations <cit.>. Microscopic bodies,
on the other hand, generally require quantum descriptions, but recent
work has shown that one can accurately represent their response
𝕍_k = ∑_p α_p |f_p⟩⟨f_p| through bases
{|f_p⟩} of either exponentially localized (for insulators)
or polynomially delocalized (for metals) functions <cit.>,
that accurately capture multipolar interactions among electronic
wavefunctions <cit.>. For molecules with finite
electronic
gaps, 
the bare response is well described by sums over dipolar ground-state
oscillator densities <cit.>,

    f_p (iξ, x⃗) = (√(2π)σ_p(iξ))^-3exp(- (x⃗ - x⃗_p)^2 /2σ_p^2(iξ)),

centered at the locations x⃗_p of each atom p, normalized
such that ∫d^3x⃗ f_p = 1, and featuring a
Gaussian width that, rather than being
phenomenological <cit.>, depends
on the atomic polarizability via σ_p (iξ) =
(α_p
    (iξ)/√(72π^3))^1/3 <cit.>. The isotropic atomic polarizabilities
α_p are computed via density functional theory, as in recent
works <cit.>, which include
short-range electrostatic, hybridization, and quantum exchange
effects.

Since microscopic and macroscopic bodies are assumed to be disjoint,
it is more efficient to partition the T-operators into blocks
belonging to either molecules or macroscopic objects, allowing a trace
over the macroscopic DOFs. The definitions of
𝕋_k(env) imply
𝕋_k(env)^-1 = 𝕍_k(env)^-1
- 𝔾_0, which means that the relevant T-operators can be
written as:

    𝕋^-1   = [ 𝕋_mol^-1     -𝔾_0;     -𝔾_0 𝕋_env^-1 ],    𝕋_∞ = [ 𝕋_mol,∞       0;       0   𝕋_env ]

thus partitioning the molecular and macroscopic (environmental)
DOFs. These depend on the molecular T-operators

    𝕋_mol^-1   = [  𝕋_1^-1    -𝔾_0       …    -𝔾_0;    -𝔾_0  𝕋_2^-1       …    -𝔾_0;       ⋮       ⋮       ⋱       ⋮;    -𝔾_0    -𝔾_0       … 𝕋_N^-1, ]

with 𝕋_mol,∞ = ∏_k 𝕋_k, which
are in turn partitioned into blocks for each of the N molecular
bodies. Given this, the product in the determinant can be evaluated
as:

    (𝕋_∞𝕋^-1)    =
      (𝕋_mol,∞𝕋_mol^-1) (𝕀 - 𝔾_0𝕋_env𝔾_0𝕋_mol) 
       = (𝕋_mol,∞𝕋_mol^-1) (𝕀 -
      𝔾_env𝕍) 
       × (𝕀 - 𝔾_0𝕍)^-1

where we used the property 𝔾_0𝕋_k(env) = (𝕀 - 𝔾_0𝕍_k(env))^-1 - 𝕀, and consolidated the
scattering properties of the macroscopic bodies into the operator
𝔾_env = 𝔾_0 (𝕀 -
𝕍_env𝔾_0)^-1, which solves

    [∇×∇× +
        ξ^2/c^2(𝕀 + 𝕍_env)]
      𝔾_env = -ξ^2/c^2𝕀

for an imaginary frequency ω = ξ, thereby encoding the
macroscopic DOFs purely in the electric field response; this can be
solved via any number of state-of-the-art analytical or numerical
classical electrodynamic techniques <cit.>, including but not limited to
scattering <cit.> and
finite-difference <cit.> methods. Moreover, as the molecules are all
disjoint, then (𝕋_mol,∞𝕋_mol^-1) = (𝕀 - 𝔾_0𝕍) ∏_k (𝕀 - 𝔾_0𝕍_k)^-1. Putting all of these identities together yields
the following expression for the energy:

    ℰ = ħ/2π∫_0^∞dξ ln[(𝕄𝕄_∞^-1)]

where 𝕄 = 𝕀 - 𝔾_env𝕍
and 𝕄_∞ = ∏_k (𝕀 - 𝔾_0𝕍_k).

The above log-determinant formula for the energy includes retardation
by construction and accounts for many-body screening and multiple
scattering to all orders, thereby ensuring full consideration of
finite size, complex shape effects, and collective polarization
excitations (see supplement for an alternate equivalent derivation
including all of these effects). Moreover, existing sophisticated
techniques for modeling molecular and electromagnetic-field responses
come together in the operator products 𝔾𝕍_k;
when represented in the p-dimensional molecular basis
{|f_p⟩}, their block matrix elements are of the form:

    f_p | 𝔾𝕍_k  f_q =
      α_q∫d^3x⃗ d^3x⃗' f_p
      (x⃗) 𝔾 (x⃗, x⃗') f_q (x⃗')

(see supplement for more details). The equivalence of (<ref>)
and (<ref>) captures the seamless unification of ideas and methods
previously confined to either atomistic vdW or continuum Casimir
physics <cit.>: (<ref>) is similar to prior
log-determinant expressions used to describe molecular interactions in
vacuum <cit.>, except that 𝔾_0 and
𝔾_env are replaced by nonretarded (quasistatic)
vacuum fields 𝔾_0(ξ = 0).

We demonstrate the importance of all of these effects by comparing the
vdW energies (or forces) obtained from (<ref>) to those from
pairwise or other approximate treatments in a number of
configurations, consisting of one or two molecules above either a gold
half-space or a conical gold tip. While the Green's function of the
half-plate can be computed analytically <cit.>, the latter
is computed using brute-force numerical techniques <cit.>, with the
dielectric function of gold taken
from <cit.>. 
We specifically study a C_500-fullerene of radius
1 nm, a 250 atom 30 nm-long linear carbyne wire,
and a 1944 atom-large
2.6 nm × 2.9 nm × 5.5 nm
protein associated with human Huntington's
disease <cit.>.

We further compare the RMB energy from (<ref>) to typical
approximations used in the literature: the non-retarded vdW energy
ℰ_0, obtained by evaluating (<ref>) with 𝔾_0
and 𝔾_env replaced by their respective quasistatic
(ξ = 0) responses, and the Casimir–Polder (CP) energy,

    ℰ_CP = -ħ/2π∫_0^∞ dξ [α·𝔾_env·(𝕀 + 1/2α·𝔾_env) ]

which ignores finite size effects by instead contracting the dressed
susceptibility of the molecular ensemble into effective dipolar
polarizabilities,

    α = ⊕_k ∑_p,q⟨ f_p |
      (𝕀 - 𝕍_k𝔾_0)^-1𝕍_k f_q⟩,

thus neglecting higher-order many-body interactions among the
different molecules and surfaces. Finally, we define a pairwise
interaction energy,

    ℰ_PWS = -ħ/2π∫_0^∞dξ [∑_k𝕍_k𝔾_env(𝕀 + 1/2∑_l ≠ k𝕍_l𝔾_env) ]

which, as in (<ref>), is obtained as a lowest-order expansion of
(<ref>) in the scattering; this captures both finite size and
retardation but ignores all high-order many-body interactions, with
the sums over k, l running over either individual or pairs of
molecules. When comparing non-retarded and CP energies to their
corresponding pairwise approximations, it suffices to take the
quasistatic limit in (<ref>) and to let (𝕀 -
𝕍_k 𝔾_0)^-1→𝕀 for the effective
polarizability α in (<ref>), respectively.



onemol shows the RMB to pairwise energy ratio
ℰ/ℰ_PWS of various
configurations (insets), with the fullerene interaction (blue line)
found to vary only slightly, attaining a maximum of 1.16 at z ≈
10 nm; such a small discrepancy stems from the small size
and isotropic shape of the fullerene, which limits possible nonlocal
correlations in its polarization response. Even weaker relative
correlations are observed in the case of the protein (green line),
which despite its greater size, number of atoms, and chemical
complexity, has a reduced response compared to semi-metallic carbon
allotropes <cit.>. To separate
the various many-body effects, the inset of onemol compares
the RMB power law ∂ln(ℰ)/∂ln(z)
of the fullerene interaction to its counterparts when neglecting
either finite size or retardation. As expected, both approximations
become accurate in their corresponding regimes of validity, with the
power law asymptoting to -4 and -1.9 at large and small z,
respectively, but fail in the intermediate, mesoscopic regime z
≈ 10 nm. Even larger discrepancies arise in the case
of the wire, whose large size and highly anisotropic shape support
long-wavelength collective fluctuations. We find that the absolute
values of both ℰ_0 (dashed red) and
ℰ_CP (dotted red) for the parallel wire
overestimate ℰ by factors of 3–7 [onemol(a)] due
to the slower decay of the Green's function in the former and lack of
screening over the length (or modes) of the wire in the latter. The
corresponding energy ratios, however, behave differently in that the
effect of screening is strongest in the quasistatic limit, which ends
up greatly dampening the many-body excitations relative to pairwise
approximations and hence leads to smaller non-retarded energy ratios;
in contrast, by construction CP ignores many-body interactions with
the surface and thus screening has a much weaker impact relative to
the pairwise approximation, leading to larger CP energy ratios. At
intermediate z ≈ 10 nm of the order of the wire
length, ℰ/ℰ_PWS≈ 30, with the
approximate energy ratios deviating by 20%. Similar results are
observed in the case of a wire in the perpendicular orientation (black
lines), with the pairwise energy leading to slightly larger
discrepancies at short separations due to the screening and decreasing
impact of atoms farther away from the plate.








































We now investigate the mutual vdW interactions among two fullerenes or
parallel wires oriented either parallel or perpendicular to the gold
plate [twomol], focusing primarily on horizontal separations
d on the order of molecular sizes, where many-body and finite size
effects are strongest. Especially in the case of two wires, the
pairwise approximation is shown to fail by many orders of magnitude,
with the largest energy ratios occurring at asymptotically large z,
i.e. for two molecules in vacuum, while at small z a decreasing
ratio reflects the dominant interactions (and screening) of the
individual molecules with the plate. The transition and competition
between the two limiting behaviors occurs at mesoscopic z ∼ d,
and is more clearly visible from the plots in twomol(lower
inset), which show ℰ/ℰ_PWS
versus d at several values of z. In particular, in the case of
parallel wires at mesoscopic z = 10 nm, the competition
leads to a nonmonotonic energy ratio, with the maximum of 70 occurring
at intermediate d ≈ 3 nm. Comparisons against
non-retarded and CP approximations illustrate behaviors similar to the
previous case of a single wire, with each under- and over-estimating
the ratios by approximately 20% and 30%, respectively. Also shown in
twomol(upper inset) is the ratio of the horizontal force F_y
= -∂ℰ/∂ y on the wires to its
pairwise counterpart, plotted against z for parallel wires at d =
10 nm. Note that by construction, F_y,PWS is
independent of z and thus, the system experiences an absolute
decrease in the force due to the screening induced by the
plate. Comparing F_y,0 and F_y,CP, one finds the
surprising result that in contrast to the energy ratio of a single
molecule, the screening by the plate makes retardation more rather
than less relevant to the force at small z, leading to an ≈
10% decrease in the force magnitude.



Finally, we consider the energy of a molecule above a gold conical
tip [cone] by comparing it to that of a gold plate at the same
vertical separation z, with 𝔾_env in the
former computed through the use of a free, surface-integral Maxwell
solver, SCUFF-EM <cit.>. The finite cone has a
base diameter of 54 nm and a height of 50 nm
from the base to the bottom of a hemispherical tip of diameter
20 nm.  The ratio decreases with increasing z, with the
energy scaling as z^-6 at asymptotically large separations (not
shown) as the finite sizes of the cone and molecule become irrelevant
and their interactions dipolar. (Note that a decreasing ratio is
expected also for a semi-infinite cone due to its smaller effective
area and hence stronger decay compared to a plate.) The ratios at
small z for the fullerene and perpendicular wire approach 1 since in
this limit, their small horizontal sizes allow the hemispherical tip,
which effectively acts like a plate at such short separations, to
dominate the interaction. By contrast, the ratio in the case of a
parallel wire is non-monotonic, decreasing with at short separations
since in this configuration, the wire excitations in the limit z →
0 still sample the finite curvature of the tip and conical slope,
leading to a different asymptotic power law.

In conclusion, we have demonstrated a unifying approach to computing
vdW interactions among molecules and macroscopic bodies that accounts
for many-body and multiple-scattering effects to all orders. By
comparing against commonplace pairwise, CP, and non-retarded
approximations, we quantified the impact of nonlocality, finite size,
and retardation on the vdW energy between molecules and either a
planar or conical macroscopic body. We have consistently found larger
deviations in approximate interactions for long, semi-metallic
molecules such as carbyne wires, whereas compact, insulating molecules
such as many proteins are reasonably well-described as effectively
dilute dielectric particles, allowing these low-order approximations
to be more valid. In the future, one might consider more complex
macroscopic bodies, such as periodic gratings <cit.> that may elicit larger differences between RMB and
approximate interactions even for compact biomolecules, as well as
extend these results to incorporate the effects of infrared molecular
resonances <cit.>.

This material is based upon work supported by the National Science
Foundation under Grant No. DMR-1454836 and by the National Science
Foundation Graduate Research Fellowship Program under Grant No. DGE
1148900.
















































§ APPROXIMATIONS


Our calculations above make two related approximations related to
Gaussian damping. First, we approximate (<ref>) as

    ⟨ f_p|𝔾𝕍_k f_q⟩≈α_q∫ d^3x⃗' 𝔾 (ξ,
      x⃗_p, x⃗') f_p+q (ξ, x⃗')
                                                               
where f_p + q is the same as f_q, but with √(2)σ_p replaced with √(σ_p^2 + σ_q^2), in
line with <cit.>; this effectively approximates the
Galerkin discretization by a collocation method, with the basis
Gaussian functions f_p + q acquiring modified widths. Secondly,
for computational convenience, we consider only scattered fields
(Green's functions) from dipolar rather than Gaussian sources, which
is justified so long as the atoms are several widths (angstroms) away
from the macroscopic surfaces.

















































































§ VDW ENERGY VIA FLUCTUATION–DISSIPATION THEOREM


We provide a heuristic derivation of the retarded, many-body (RMB) vdW
energy of a general collection of molecular or macroscopic bodies,
requiring only that they be disjoint and have no correlations in the
polarization response between bodies. Each body k is
described by an electric susceptibility
𝕍_k, relating its polarization to the total electric
field via |P⃗_k⟩ = 𝕍_k|E⃗⟩; these
susceptibilities account for short-range quantum and electrostatic
correlations, allowing us to focus solely on long-range electrodynamic
correlations when considering the vdW energy. [When
  neglecting retardation, the charge density and electric
  potential are more frequently used, so the density response is
  written in terms of the susceptibility as χ(ω, x⃗,
  x⃗') = ∑_i,j∂_i∂_j
  (𝕍_k)_ij (ω, x⃗, x⃗').] Our derivation
follows analysis <cit.>
based on the fluctuation–dissipation theorem; we note previous
demonstrations <cit.> of its
equivalence to the summation of ground-state energies of the coupled
molecular system <cit.>.

Following RosaPRA2011, the assembly of the constituents
of all bodies from infinite separation into the final configuration
defining 𝕍 = ∑_k𝕍_k can be considered the
result of an adiabatic change in the particle–field coupling strength
λ∈ [0, 1], in which case the energy of the system can be
written as,

    ℰ = -∫_0^∞ dω ∫_0^1 dλ/λ P⃗|E⃗_ZP,

per the Feynman–Hellmann theorem <cit.>. Here, |E⃗⟩ denotes zero-point
fluctuating electric fields, |P⃗⟩ = 𝕍|E⃗⟩ is the induced polarization, and
_ZP denotes the quantum statistical average
over zero-point fluctuations. The connection to scattering problems
comes from the well-known fluctuation–dissipation
theorem <cit.>,

    |E⃗⟩⟨E⃗|_ZP = 
      ħ/π𝔾,

which expresses field fluctuations in terms of the Green's function
𝔾 of the system. The latter solves Maxwell's equations and
can be written in terms of the susceptibility as 𝔾 =
(𝕀 - 𝔾_0𝕍)^-1𝔾_0 <cit.>. Exploiting the analyticity of
𝕍 and 𝔾_0 in the complex-ω
plane <cit.> and performing a
Wick rotation of the energy integral from real to imaginary frequency
ω = ξ, leads to a simplified expression for the
energy [The entirety of this derivation is identical to that
  of past work employing the so-called adiabatic connection
  fluctuation–dissipation (ACFD) framework, but using the vector
  polarization, tensor electric susceptibility, and tensorial vacuum
  Green's function instead of the scalar charge density, density
  response, or Coulomb potential, in order to account for
  retardation. It is therefore not a coincidence that the
  log-determinant frequency integrand is so similar in form to past
  expressions for the vdW energy of a single body.],

    ℰ = ħ/2π∫_0^∞ dξ ln((𝕀 - 𝔾_0𝕍))

where we rescaled the response functions 𝔾_0 and
𝕍 by the coupling constant λ and integrated over
λ. The net interaction energy among the bodies is found by
subtracting self-energies of the form in (<ref>), replacing
𝕍 by 𝕍_k separately for each k. This allows
recasting the net vdW interaction energy as (<ref>) in terms of
scattering operators:

    𝕄   = 𝕀 - 𝔾_0𝕍
    𝕄_∞   = ∏_k (𝕀 - 𝔾_0𝕍_k).


If the system considered consists of N molecular bodies labeled k,
and an arbitrary number of macroscopic bodies collectively described
by 𝕍_env, then one can write

    𝕄   = 𝕀 - 𝔾_0𝕍_env
      - 𝔾_0𝕍
    𝕄_∞   = (𝕀 - 𝔾_0𝕍_env) ∏_k (𝕀 - 𝔾_0𝕍_k),

where 𝕍 = ∑_k𝕍_k only runs over the
molecular bodies. Multiplying 𝕄𝕄_∞^-1
produces terms of the form

    (𝕀 - 𝔾_0𝕍_env)^-1𝔾_0≡𝔾_env,

which is just the electric field response due to
𝕍_env alone and can be computed via analytical
or numerical formulations of continuum electrodynamics. Redefining

    𝕄'    = 𝕀 - 𝔾_env𝕍
    𝕄'_∞   = ∏_k (𝕀 - 𝔾_0𝕍_k),

and dropping primes, these new operators can then be substituted
into (<ref>) to obtain the net vdW interaction energy among N
molecules and a general macroscopic environment.

















































































































































































































