<figure><img src="image_1.png" /></figure>

HAL
open science

# The Casimir effect derived from scattering theory

Romain Guérout

► To cite this version:

Romain Guérout. The Casimir effect derived from scattering theory. Physics [physics]. Sorbonne Université, 2022. tel-03808005

HAL Id: tel-03808005

https://hal.science/tel-03808005

Submitted on 10 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Habilitation (HDR) Thesis

Romain GUÉROUT

Laboratoire Kastler Brossel

CNRS et Sorbonne Université campus Pierre et Marie Curie

Corridor 13-23 2nd floor, 75005 Paris, France

Email : guerout@lkb.upmc.fr

# The Casimir effect derived from scattering theory

## Composition of the Jury

Éric Charron: Reviewer

Professor at ISMO, Paris-Saclay University

Gabriel Dutier: Reviewer

Associate Professor at LPL, Sorbonne University Paris Nord

Nicolas Pavloff: Reviewer

Professor at LPTMS, Paris-Saclay University

Marie-Christine Angonin: Examinator

Professor at SYRTE, Sorbonne University

David Guéry-Odelin: Examinator

Professor at LCAR, Toulouse University

Catherine Schwob: Examinator

Professor at INSP, Sorbonne University

Defended publicly on September 13, 2022
# Curriculum Vitæ

* Romain Guérout, né le 31 Juillet 1976 (45 ans)

* Étudiant à l'Université Paris-Sud (Orsay)

* 2004 : Diplôme de doctorat du département de chimie de l'université de Bâle (Suisse) intitulé « De l'utilisation des calculs *ab initio* appliqués à la théorie moléculaire du défaut quantique »

* 2004-2011 : Différents séjours post-doctoraux

* 2004-2006 : à Steacie Institute for Molecular Sciences / CNRC (Ottawa, Canada) sous la supervision de P. Bunker. Étude théorique du spectre rovibronique de molécules triatomiques sous l'influence de l'effet Renner-Teller.

* 2006-2009 : au Laboratoire Aimé Cotton sous la supervision de Ch. Jungen et O. Dulieu. Théorie moléculaire du défaut quantique et calculs de structures électronique dans les molécules diatomiques alcalines.

* 2009-2011 : au Laboratoire Kastler Brossel sous la supervision de A. Lambrecht. Effet Casimir pour des géométries complexes.

* Depuis 2011 : recruté Chargé de Recherche au CNRS. Affecté au Laboratoire Kastler Brossel dans le groupe « Fluctuations Quantiques et Relativité » sous la supervision de A. Lambrecht et S. Reynaud.

* 50+ publications dans des revues scientifiques à comité de lecture.

## Co-encadrement d'étudiants en thèse

* **Johan Lussange**, soutenance le 10 septembre 2012 (https://tel.archives-ouvertes.fr/tel-00879989v1) intitulée « The Casimir energy and radiative heat transfer between nanostructured surfaces ». Encadrement à 50%.

Directrice de thèse : Astrid Lambrecht (50%)

* **Gabriel Dufour**, soutenance le 20 novembre 2015 (https://tel.archives-ouvertes.fr/tel-01242290v2) intitulée « Quantum reflection from the Casimir-Polder potential ». Encadrement à 25%.

Directrice de thèse : Astrid Lambrecht (50%)

Autre encadrant : Serge Reynaud (25%)

* **Axel Maury**, soutenance le 27 septembre 2016 (https://tel.archives-ouvertes.fr/tel-01478886v1) intitulée « Effet Casimir-Polder sur des atomes piégés ». Encadrement à 50%.

Directrice de thèse : Astrid Lambrecht (50%)

* **Pierre-Philippe Crépin**, soutenance le 3 Juin 2019 (https://tel.archives-ouvertes.fr/tel-02155141) intitulée « Quantum reflection of a cold antihydrogen wave packet ». Encadrement à 25%.

Directeur de thèse : Serge Reynaud (75%)

## Co-encadrement d'étudiants post doctoraux

* **Manuel Donaire**, « Effet Casimir-Polder sur des atomes piégés »
  Encadrement à 10%.

Responsable : Astrid Lambrecht (90%)

* **Antoine Canaguier-Durand**, « Nanophotonique, couplage lumière-matière cohérent »
  Encadrement à 10%.

Responsable : Serge Reynaud (90%)

## Co-encadrement de stagiaires

* Encadrement à 50% de différents stagiaires.
Habilitation à Diriger des Recherches

The Casimir effect derived from scattering theory

Romain Guérout

April 7, 2022
# Contents

1 Introduction 3

2 The Casimir effect 3

2.1 Casimir interaction energy between scatterers 3

2.1.1 The scattering operator of an ensemble of scatterers 5

2.1.2 The determinant of the scattering operator 8

2.1.3 Expression of the Casimir interaction free energy 9

2.1.4 Explicit expressions for few scatterers 10

2.1.5 Expression in terms of imaginary frequencies 11

2.2 Applications 13

2.2.1 Plane-parallel mirrors: The Lifshitz formula 13

2.2.2 Drude-plasma models for metallic mirrors 16

2.2.3 Grating geometry 22

2.2.4 Quantum reflection and gravitational quantum states 25

2.3 Miscellaneous 31

2.3.1 Van der Waals interaction between two different atoms 31

2.3.2 Multiscattering formalism 37

2.3.3 ForcaG experiment 45

3 Conclusion 53
# 1 Introduction

This manuscript describes some of my work at the Laboratoire Kastler Brossel on Casimir physics. I joined the group “Fluctuations Quantiques et Relativité” at the LKB in 2009 and got recruited in 2011 in the CNRS as Chargé de Recherche.

Casimir physics study the many phenomena related to fluctuations-induced forces. Those include the van der Waals interaction between atoms or molecules as well as the Casimir interaction between macroscopic bodies.

There are two main interpretations for Casimir interactions. In the so-called “field” interpretation, the vacuum field fluctuates and the bodies act merely as boundary conditions for those fluctuations. On the other hand, in the so-called “source” interpretation, fluctuations originates inside the matter forming the bodies and propagates through the vacuum field. The field interpretation was used by Casimir in his original paper [1] where he studied the interaction between perfectly reflecting plane-parallel plates. The fact that the plates were perfectly reflecting allowed to enforce hard boundary conditions on the vacuum field. The source interpretation is used to explain the van der Waals interaction where a fluctuating dipole in one atom creates a field which polarize the other atom which then acts back on the first atom.

There is no “correct” interpretation and source and field interpretations have been shown to be completely equivalent [2]. In fact, for the situation of explaining spontaneous emission of an excited atom, it was shown that one can go from one interpretation to the other simply by changing the ordering of atomic and field operators used to calculate expectation values [3]. The same line of reasoning has been done for the Casimir interaction [4].

There is a surprisingly large number of different methods to calculate the Casimir interaction. Most of those methods though, being based on quantum field theory, often require some regularization methods in order to get rid of infinite quantities. In this manuscript, we present a method of calculating the Casimir interaction between an arbitrary number of bodies. The method is expressed solely in terms of the scattering properties of the bodies in interaction and do not require any regularization.

We begin by carefully deriving our method of calculation based on the scattering theory and proceed to show some applications. Those applications include comparison with Casimir experiment using grating surfaces, quantum reflection which arise from the Casimir-Polder interaction between an atom and a material surface and other related phenomena in Casimir physics.

# 2 The Casimir effect

## 2.1 Casimir interaction energy between scatterers

In this section, we derive an expression for the Casimir interaction energy between an arbitrary number of scatterers in thermal equilibrium with their surrounding which we take to be vacuum. We base our demonstration on the fundamental concept which is the electromagnetic local density of states (LDOS) and its relation with
scattering operators.

As an introduction, let us first recover the energy density of the blackbody radia-
tion at temperature T. Consider a blackbody cavity of volume V. Thermodynamics
and quantification of the field inside of the cavity leads to a mean (internal) energy
$U_\omega$ for a given frequency

$$
Z^{-1} = 2 \sinh \frac{\beta \hbar \omega}{2} \tag{1}
$$

$$
U_{\omega} = - \frac{\partial \log Z}{\partial \beta} = \hbar \omega \left( n_{Pl} + \frac{1}{2} \right) \quad (2)
$$

$$
n_{Pl} = (e^{\beta\hbar\omega} - 1)^{-1} \qquad (3)
$$

with $\beta^{-1} = k_B T$. In vacuum, the density of states per unit volume per unit angular frequency is $\rho_0(\omega) = \frac{\omega^2}{\pi^2 c^3}$. The total energy inside the cavity is given by

$$
U = \int_V d^3\mathbf{r} \int_0^\infty d\omega \, \rho_0(\omega) U_\omega = V \int_0^\infty d\omega \, \frac{\hbar\omega^3}{\pi^2 c^3} \left(n_{Pl} + \frac{1}{2}\right) \quad (4)
$$

The observable blackbody energy density per unit volume per unit angular frequency $u_\omega$ is the part of the above expression without the zero-point energy. We recover the usual expression

$$
u_{\omega} = \frac{\hbar \omega^3}{\pi^2 c^3} n_{Pl} \qquad (5)
$$

From this example and eq. (4), we see that any observable can be calculated as an
integral of the spectral density of said observable weighted by the density of states.
We will now use this formalism in order to derive the Casimir interaction energy. We
place in the cavity an ensemble of *N* scatterers. The scatterers are unspecified, we
merely require them to be of finite volume. In that case, the ensemble of scatterers
can always be embedded in a (fictitious) sphere of radius *R* with respect to which a
scattering setup can be properly defined. The mere presence of the scatterers makes
the density of states ρ(r, ω) a local quantity which now depends of the position r
inside the cavity.

We define Δρ(r, ω) = ρ(r, ω) − ρ₀(ω) the difference of the LDOS with respect to
the situation without any scatterers. Importantly, Δρ(r, ω) vanishes for r far from
any scatterers surfaces. We can then increase the size of the blackbody cavity to
infinity and define

$$
\Delta D(\omega) = \int d^3\mathbf{r} \Delta\rho(\mathbf{r}, \omega) \qquad (6)
$$

and refer to *D*(ω) as the global density of states [5]. We cannot say that the quantity
Δ*D*(ω) is finite because of the behaviour of the LDOS ρ(*r*, ω) on the surface of the
scatterers but at least it does not depend explicitly on the volume of the blackbody
cavity, as it should. Next, instead of the internal energy (which can be calculated) we
are now rather interested in the Helmholtz free energy *A* or rather its modification
because of the presence of the scatterers

$$
A_{\omega} = -\frac{1}{\beta} \log Z \qquad (7)
$$

$$
\Delta A_{vac} = \int d\omega \Delta D(\omega) A_{\omega} \qquad (8)
$$
Again, the quantity $\Delta A_{vac}$ is the modification of the vacuum free energy induced by the presence of the ensemble of scatterers. We are now in position of making the connection with the scattering theory. Indeed, we now invoke the fundamental relation between global density of states and Wigner time delay [6]. Defining $\mathbf{Q}$, the Wigner-Smith time-delay operator [7], we have

$$ \mathbf{Q} = -i\mathbf{S}^{\dagger}\frac{\partial\mathbf{S}}{\partial\omega} = -i\frac{\partial}{\partial\omega}\log\mathbf{S} \quad (9) $$

$$ \Delta D(\omega) = \text{tr} \mathbf{Q}/2\pi = \frac{1}{2i\pi} \frac{\partial}{\partial\omega} \log \det \mathbf{S} \quad (10) $$

where the last equation is nothing more than the Krein-Friedel-Lloyd formula. In the above equations, we have used the unitarity of the scattering operator $\mathbf{S}^{\dagger} = \mathbf{S}^{-1}$ and the fact that it is trace-class so that $\text{tr}\log = \log\det$. As mentioned earlier, the scattering operator $\mathbf{S}$ can be defined with respect to a sufficiently large sphere enclosing all scatterers. Putting everything together, we arrive at

$$ \Delta A_{vac} = \frac{1}{2i\pi} \int d\omega \left[ \frac{\partial}{\partial\omega} \log \det \mathbf{S} \right] A_{\omega} \quad (11) $$

Finally, we now define the Casimir interaction free energy $A_{int}(\mathbf{L})$ as

$$ A_{int}(\mathbf{L}) = \Delta A_{vac}(\mathbf{L}) - \Delta A_{vac}^{\infty} \quad (12) $$

In the above equation, the quantity $\mathbf{L}$ symbolically represents the position of all scatterers with respect to each others. The quantity $\Delta A_{vac}^{\infty}$ corresponds to the situation where all scatterers are infinitely separated from each others. We see that in order to calculate the Casimir interaction energy, one has to extract from the scattering operator $\mathbf{S}$ (or rather its determinant) only the part which explicitly depends on the position of the scatterers with respect to each others. This will constitute the next section.

### 2.1.1 The scattering operator of an ensemble of scatterers

Let us briefly recall the definition of the scattering operator $\mathbf{S}$ first in the case of one scatterer. Given a bounding surface enclosing the scatterer and $\{V_{-}, V_{+}\}$ basis states for the incoming and outgoing fields. We have

$$ \psi_{\text{in}}(\mathbf{r}) = V_{-}(\mathbf{r})c_{\text{in}} \quad (13) $$

$$ \psi_{\text{out}}(\mathbf{r}) = V_{+}(\mathbf{r})c_{\text{out}} \quad (14) $$

Then, the outgoing amplitudes are related to the incoming ones by

$$ c_{\text{out}} = S c_{\text{in}} \quad (15) $$

A few words about the basis states $V_{\pm}$: those are line vectors collecting a number of basis functions $\varphi_k(\mathbf{r})$. Concretely, those could be for instance a number of plane waves or vector spherical harmonics. The index $k$ collects any number of quantum numbers defining the basis (frequency, wavector, angular momentum, etc...). In a
more abstract way, the basis states $V_+$ collects any number of scattering channels or scattering ports. In a lossless embedding medium, the outgoing channels $V_+$ can be chosen to be the time-reversal of the incoming ones $V_-$.

In the case of an ensemble of $N$ scatterers, the arguments above basically still holds. We need a bounding surface enclosing all scatterers. It is convenient to express the total scattering operator $S$ in a multicenter basis: given the position vectors of each scatterers $\mathbf{r}_i$, $1 < i < N$, the basis states $V_+$ will now consist of functions $\varphi_{\mp k}(\mathbf{r} - \mathbf{r}_i)$ centered on each scatterer. We now introduce the concept of internal and external scattering channels. Intuitively, a set of internal scattering channels $\varphi_i^+$ is one linking two different scatterers. In other words, internal channels between scatterer $n$ and $m$ are such that

$$ \varphi_{-}^{i}(\mathbf{r} - \mathbf{r}_{m}) = \varphi_{+}^{i}(\mathbf{r} - \mathbf{r}_{n})\mathbf{T}_{nm} \quad (16) $$

meaning that incoming channels on scatterer $n$ are obtained from applying the translation operator $\mathbf{T}_{nm}$ to the outgoing channels from scatterer $m$. We then collect in the external channels all other scattering channels. External scattering channels include for instance channels connecting each scatterer to the enclosing bounding surface.

<figure><img src="image_3.png" /><figcaption>Figure 1: Scattering diagrams for the case of $N = 2$ and $N = 3$ scatterers. Each diagram consists of $N$ external ports (blue), $N$ external channels (solid lines), $N(N-1)$ internal ports (red) and $N(N-1)/2$ internal channels (dashed lines).</figcaption></figure>

Obviously, the total scattering operator $S$ will be expressed in terms of the external channels only. Indeed, the derivation of $S$ is a known problem which has been tackled both in a very concrete way [8] or in a more abstract way [9, 10]. It is based on the so-called *internal channels elimination* method which is formally similar to the *closed channels elimination* method of multichannel quantum defect
theory. In the concrete way, each scatterer is a sphere and the basis states are vector spherical harmonics leading to the so-called generalized multiparticle Mie-solution formalism. In the abstract way, the total scattering operator is obtained for unspecified scatterers and basis states. We will briefly rederive the scattering operator S in a way suited to our purpose. In the following, we will work with operators in block-matrix form. We first define the operators S_n, the scattering operator of each individual scatterer. These operators act upon the basis functions $\varphi_{\mp}(\mathbf{r} - \mathbf{r}_n)$ some of which are internal channels if they link to other scatterer. This leads to the following block structure

$$S_n = \begin{pmatrix} S_n^{ee} & S_n^{ei} \\ S_n^{ie} & S_n^{ii} \end{pmatrix} \qquad (17)$$

We consider that each of the scatterer can be absorbing i.e. it is described by complex permittivities and permeabilities. This leads, a priori, to non-unitary scattering operators S_n and S. However, in the derivation of eq. (9), we have used the unitarity of S. Therefore, we have to include in the operators S_n any channels responsible for the absorption. We include these absorption processes in additional external channels. By doing so, we now work with unitary individual operators S_n which will necessarily lead to a unitary total operator S. As a consequence of unitarity, we have the following property (which we will use later) for the individual scattering operators:

$$\det S_n = \frac{\det(S_n^{ee})}{\det(S_n^{ii})^*} \quad \text{if} \quad S_n^\dagger = S_n^{-1} \qquad (18)$$

Next, we introduce two auxiliary multi-centered operators S and W. Those act on the full basis states $V_{\mp}$ which include all channels i.e. internal and external centered on all scatterers. The operator S collects all individual scattering operators as if the different scatterers were not interacting with each others. We choose to partition it according to (multi-centered) external and internal channels. Each blocks are themselves block-diagonal in terms of the different scatterers:

$$S = \begin{pmatrix} S^{ee} & S^{ei} \\ S^{ie} & S^{ii} \end{pmatrix} = \begin{pmatrix} S_1^{ee} & S_1^{ei} \\ \vdots & \vdots \\ S_N^{ee} & S_N^{ei} \\ S_1^{ie} & S_1^{ii} \\ \vdots & \vdots \\ S_N^{ie} & S_N^{ii} \end{pmatrix} \qquad (19)$$

Note that the off-diagonal blocks $S_{ie}$ and $S_{ei}$ are in general non-square since there are N external ports and $N(N-1)$ internal ones.

Then, the operator W is the weighted adjacency matrix of the graph associated to the situation. The weights over each internal channel being simply the relevant translation operator T. Obviously, by definition the operator W only acts over the
internal channels. It then has then following structure:

$$W = \begin{pmatrix} 0 & 0 \\ 0 & W^{ii} \end{pmatrix} \qquad (20)$$

$$[W^{ii}]_{nm} = \begin{cases} \mathbf{T}_{kl} & \text{if } (nm) \text{ is an internal channel } \mathbf{k} \to \mathbf{l} \\ \mathbf{0} & \text{otherwise} \end{cases} \qquad (21)$$

Above, the indexes n, m make reference to internal ports. The notation $\mathbf{k}$, $\mathbf{l}$ make reference to particular scatterer. At the end of the section, we will give explicit expressions for $W$ in simple cases of two and three scatterers.

We then have the operator $\mathcal{S}$ which only depends on the scattering properties of each scatterer, independently of their relative positions. On the other hand, the operator $W$ only depends on the relative position of each scatterer, independently of their scattering properties. Remarkably, the total scattering operator $\mathbf{S}$ of the ensemble of scatterer has a rather simple form. It is the part acting on the external channels of the operator $(1 - \mathcal{S}W)^{-1}\mathcal{S}$ i.e.

$$\mathbf{S} = \mathbf{X} (1 - \mathcal{S}W)^{-1} \mathcal{S} \mathbf{X}^{\mathrm{T}} \qquad (22)$$

if $\mathbf{X}$ denotes a projector on the external channels. Using the above partitioning of the operators $\mathcal{S}$ and $W$ and the expression of the inverse of a 2 × 2 block matrix, we obtain the explicit expression for $\mathbf{S}$

$$\mathbf{S} = \mathcal{S}^{ee} + \mathcal{S}^{ei}W^{ii}(1 - \mathcal{S}^{ii}W^{ii})^{-1}\mathcal{S}^{ie} \qquad (23)$$

### 2.1.2 The determinant of the scattering operator

Having obtained an explicit expression for the total scattering operator $\mathbf{S}$ of an ensemble of $N$ scatterers eq. (23), we are left evaluating its determinant. We will use the so-called *matrix determinant lemma* which enables us to manipulate the determinant of the sum of two matrices.

We first obtain

$$\det \mathbf{S} = \frac{\det(\mathcal{S}^{ee})}{\det(1 - \mathcal{S}^{ii}W^{ii})} \det [1 - (\mathcal{S}^{ii} - \mathcal{S}^{ie}\mathcal{S}^{ee-1}\mathcal{S}^{ei})W^{ii}] \qquad (24)$$

As a consequence of the unitarity of the individual scattering operators $\mathbf{S}_n$, we have

$$\mathbf{S}_n^{ie} \mathbf{S}_n^{ee-1} \mathbf{S}_n^{ei} = \mathbf{S}_n^{ii} - \mathbf{S}_n^{ii\dagger-1} \qquad (25)$$

So, we get

$$\det \mathbf{S} = \frac{\det(\mathcal{S}^{ee})}{\det(1 - \mathcal{S}^{ii}W^{ii})} \det(1 - \mathcal{S}^{ii\dagger-1}W^{ii}) \qquad (26)$$

After additional use of the matrix determinant lemma on the last term of the above equation, we get

$$\det \mathbf{S} = \frac{\det(\mathcal{S}^{ee})}{\det(\mathcal{S}^{ii})^*} \det(-W^{ii}) \frac{\det(1 - W^{ii-1}\mathcal{S}^{ii\dagger})}{\det(1 - \mathcal{S}^{ii}W^{ii})} \qquad (27)$$
Let us discuss each terms of this expression. The first one, thanks to the property in eq. (18) and the fact that both $S^{ee}$ and $S^{ii}$ are block-diagonal gives

$$ \frac{\det(S^{ee})}{\det(S^{ii})^*} = \prod_n \det S_n \qquad (28) $$

For the second term, we note that the dimension of $W^{ii}$, in terms of blocks, is $N(N-1)$ which is the number of internal ports and which is an even number. So

$$ \det(-W^{ii}) = \det(W^{ii}) \qquad (29) $$

For the last term, we will use the fact that in a lossless embedding medium $W^{ii}$ is unitary. Indeed, it is easy to see that the matrix $W^{ii}W^{ii\dagger}$ is block-diagonal with elements $T_{kl}T_{kl}^{\dagger} = 1$ if the medium is lossless. From the unitarity of $W^{ii}$, we get

$$ \frac{\det(1 - W^{ii-1}S^{ii\dagger})}{\det(1 - S^{ii}W^{ii})} = \frac{\det(1 - S^{ii}W^{ii})^*}{\det(1 - S^{ii}W^{ii})} \qquad (30) $$

Putting everything together, we obtain our main result

$$ \det S = \prod_n \det S_n \det(W^{ii}) \frac{\det(1 - S^{ii}W^{ii})^*}{\det(1 - S^{ii}W^{ii})} \qquad (31) $$

Under this form, it becomes evident that $\det S$ is a complex number of unit modulus providing that $S_n$ and $W^{ii}$ are unitary.

### 2.1.3 Expression of the Casimir interaction free energy

We have derive the expression of the determinant of the total scattering operator $S$ of an ensemble of $N$ scatterers. In the subtraction in eq. (12), only the distance-dependent part of this determinant will remain. It means that the first term $\prod_n \det S_n$ which represent a kind of self-energy for each scatterer drops out. We also discard the second term $\det(W^{ii})$. Indeed, this term originates because we work in a multicentered basis. The outgoing field, analyzed on the enclosing surface, is a sum of fields originating from each scatterer. Its phase-shift with respect to the incoming field must take into account the multicentered character of our basis via all the translations from one scatterer to the other encoded in $W^{ii}$. In the expression of the Casimir interaction free energy, we only keep the last term in eq. (31) to get

$$ A_{int}(\mathbf{L}) = -\frac{1}{\pi} \int d\omega \left[ \frac{\partial}{\partial\omega} \mathrm{Im} \log \det(1 - S^{ii}W^{ii}) \right] A_{\omega} \qquad (32) $$

Finally, we wish to transform a bit the above expression. It can be shown that the quantity $\mathrm{Im} \log \det(1 - S^{ii}W^{ii}) A_{\omega}$ vanishes both for $\omega = 0$ and $\omega \to \infty$. Under this condition, we can perform an integration by parts on eq. (32). Noting that $\frac{\partial A_{\omega}}{\partial \omega} = \frac{\hbar}{2} \coth \beta \hbar \omega / 2$, we obtain our final expression

$$ A_{int}(\mathbf{L}) = \frac{\hbar}{2\pi} \int d\omega \coth \beta\hbar\omega/2 \left[ \mathrm{Im} \log \det(1 - S^{ii}W^{ii}) \right] \qquad (33) $$
Other quantities of interest are the Casimir interaction force $F_{int} = -\partial A_{int}/\partial \mathbf{L}$ as the change in interaction energy as a function of a change in the arrangements of the scatterers and the interaction entropy $S_{int} = -\partial A_{int}/\partial T$. The derivative with respect to $\mathbf{L}$, the relative positions of the scatterers, acts only on $\mathcal{W}^{ii}$ and the derivative with respect of the temperature acts only on the hyperbolic cotangent (we neglect the variation with temperature of the scattering properties in $\mathcal{S}^{ii}$). Under the additional assumption that we can interchange derivatives and integrals, we get

$$F_{int}(\mathbf{L}) = \frac{\hbar}{2\pi} \int d\omega \coth \beta\hbar\omega/2 [\mathrm{Im} \, \mathrm{tr} (1 - \mathcal{S}^{ii}\mathcal{W}^{ii})^{-1} \partial_{\mathbf{L}} \mathcal{S}^{ii}\mathcal{W}^{ii}] \quad (34)$$

$$S_{int}(\mathbf{L}) = -\frac{\beta\hbar^2\omega}{4\pi T} \int d\omega \, \mathrm{csch}^2\beta\hbar\omega/2 [\mathrm{Im} \, \mathrm{log} \, \mathrm{det}(1 - \mathcal{S}^{ii}\mathcal{W}^{ii})] \quad (35)$$

The internal energy $U_{int}$ can be obtained plugging $U_{\omega}$ in eq. (32) or equivalently from the thermodynamic relation $A_{int} = U_{int} - TS_{int}$. Finally, we note that eq. (33) has a non-vanishing expression for $T \to 0$ giving the Casimir interaction energy $E_{int}$

$$E_{int}(\mathbf{L}) = \lim_{T \to 0} A_{int}(\mathbf{L}) = \frac{\hbar}{2\pi} \int d\omega \, \mathrm{Im} \, \mathrm{log} \, \mathrm{det}(1 - \mathcal{S}^{ii}\mathcal{W}^{ii}) \quad (36)$$

### 2.1.4 Explicit expressions for few scatterers

The vast majority of problem studied in Casimir physics involve the interaction between two objects. We will therefore give an explicit expression of the Casimir interaction energy in that case. We then give an explicit expression for the case of three objects, a case which is yet tractable explicitly and which illustrate concretely how to extrapolate to $N$ objects.

Consider then two objects in interaction. There is $N = 2$ external ports from which originate 2 external channels. There is also $N(N - 1) = 2$ internal ports which lead to 1 internal channel. In this simple case, we have

$$\mathcal{S}^{ii} = \begin{pmatrix} \mathbf{S}_1^{ii} & \mathbf{0} \\ \mathbf{0} & \mathcal{S}_2^{ii} \end{pmatrix} \quad (37)$$

$$\mathcal{W}^{ii} = \begin{pmatrix} \mathbf{0} & \mathbf{T}_{12} \\ \mathbf{T}_{21} & \mathbf{0} \end{pmatrix} \quad (38)$$

$$\det(1 - \mathcal{S}^{ii}\mathcal{W}^{ii}) = \begin{vmatrix} \mathbf{1} & -\mathcal{S}_1^{ii}\mathbf{T}_{12} \\ -\mathcal{S}_2^{ii}\mathbf{T}_{21} & \mathbf{1} \end{vmatrix} = \det(1 - \mathcal{S}_1^{ii}\mathbf{T}_{12}\mathcal{S}_2^{ii}\mathbf{T}_{21}) \quad (39)$$

It turns out that the case with two scatterers is somewhat of a special case because $N = N(N - 1)$ and do not illustrate how the situation evolves with $N > 2$. The reason is that in order to write down the operators $\mathcal{S}$ and $\mathcal{W}$ one has to choose how to index the internal ports i.e. which internal port is attached to which scatterer. For two scatterers, this indexing is trivial: port 1 attached to scatterer ① and port 2 attached to scatterer ②. One has to go to three scatterers at least to have a more general picture.

We now consider three objects in interaction. There are $N = 3$ external ports and $N(N - 1) = 6$ internal ones. We want to keep the block-diagonal structure of the operator $\mathcal{S}$. A natural indexing method is to sequentially associate internal ports
with scatterers i.e. ports 1 and 2 attached to scatterer ①, ports 3 and 4 attached to scatterer ② and so on. Of course, the end result will not depend on how we choose that indexing method. With this indexing method, we have

$$S^{ii} = \begin{pmatrix} S_1^{ii} & 0 & 0 \\ 0 & S_2^{ii} & 0 \\ 0 & 0 & S_3^{ii} \end{pmatrix} \qquad (40)$$

$$W^{ii} = \begin{pmatrix} 0 & 0 & 0 & 0 & T_{13} \\ 0 & 0 & T_{12} & 0 & 0 \\ 0 & T_{21} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & T_{23} \\ 0 & 0 & 0 & T_{32} & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & \tilde{T}_{12} & \tilde{T}_{13} \\ \tilde{T}_{21} & 0 & \tilde{T}_{23} \\ \tilde{T}_{31} & \tilde{T}_{32} & 0 \end{pmatrix} \qquad (41)$$

$$\det(1 - S^{ii}W^{ii}) = \begin{vmatrix} \mathbf{1} & -S_1^{ii}\tilde{T}_{12} & -S_1^{ii}\tilde{T}_{13} \\ -S_2^{ii}\tilde{T}_{21} & \mathbf{1} & -S_2^{ii}\tilde{T}_{23} \\ -S_3^{ii}\tilde{T}_{31} & -S_3^{ii}\tilde{T}_{32} & \mathbf{1} \end{vmatrix} \qquad (42)$$

$$= \det(\mathbf{1} - S_1^{ii}\tilde{T}_{12}S_2^{ii}\tilde{T}_{21}) \det(\mathbf{1} - S_1^{ii}\tilde{T}_{13}S_3^{ii}\tilde{T}_{31}) \\ \times \det[\mathbf{1} - (S_2^{ii}\tilde{T}_{23} + S_2^{ii}\tilde{T}_{21}S_1^{ii}\tilde{T}_{13})(\mathbf{1} - S_3^{ii}\tilde{T}_{31}S_1^{ii}\tilde{T}_{13})^{-1} \\ (S_3^{ii}\tilde{T}_{32} + S_3^{ii}\tilde{T}_{31}S_1^{ii}\tilde{T}_{12})(\mathbf{1} - S_2^{ii}\tilde{T}_{21}S_1^{ii}\tilde{T}_{12})^{-1}] \qquad (43)$$

As we see, an explicit formula becomes quickly intractable. Nevertheless, the above expression is interesting : the first two terms describe the Casimir interaction between scatterers ① ↔ ② and ① ↔ ③. The last term though cannot be so simple and reflect the fact that Casimir (or van der Waals) interactions are non-additive. It is the interaction between scatterers ② ↔ ③ taking into account the presence of scatterer ①.

### 2.1.5 Expression in terms of imaginary frequencies

It turns out that the final expression in eq. (33) is of little practical use. The integral over real frequencies $\omega$ is tedious and it is hard to obtain an accurate numerical result. It is because the end result of the integral involves near-cancellations between large positive and negative parts. In addition to that, the integral converges slowly at large frequencies because the integrand has an oscillatory behavior. Hopefully, we are able to obtain a particularly simple and easy-to-evaluate expression using contour integration in the complex plane as we will see. We introduce a complex frequency $z = \omega + i\xi$ and the function $\mathcal{E}(\omega)$:

$$\mathcal{E}: \mathbb{R} \rightarrow \mathbb{C} \\ \omega \rightarrow \log \det(1 - S^{ii}(\omega)W^{ii}(\omega)) \qquad (44)$$

where we have explicitly written down that both $S^{ii}$ and $W^{ii}$ are in general functions of the frequency $\omega$. It can be shown that the function $\mathcal{E}$ can be extended to complex arguments. We will come back to this point later. We then know how to define
the function $\mathfrak{e}(z)$. On the real axis, we have the property $\mathfrak{e}^*(\omega) = \mathfrak{e}(-\omega)$ meaning that the real part of $\mathfrak{e}$ is an even function of $\omega$ whereas its imaginary part is an odd function. In the whole complex plane, this property becomes a Schwartz reflection principle along the imaginary axis i.e. $\mathfrak{e}^*(z) = \mathfrak{e}(-z^*)$. Then, we define the function $f$:

$$f: \mathbb{C} \to \mathbb{C} \\ z \to \coth \beta \hbar z/2 [\mathfrak{e}(z)] \qquad (45)$$

Our goal is to apply Cauchy's theorem on the integral of $f(z)$ over a well chosen contour. Let us first discuss the analytic properties of $f(z)$. The total $\mathbb{S}$ matrix given by eq. (23), when extended in the complex frequency plane, has poles corresponding to the modes of the system of scatterers. This can be the modes corresponding to each individual scatterer (the poles of $\mathbb{S}$) or the modes of the complete system resulting from a particular spatial arrangements of the scatterers (the zeroes of $\det(1 - \mathbb{S}^{ii}\mathbb{W}^{ii})$). In the expression of the interaction energy, only the latter are present. Application of the logarithm function then leads to logarithmic singularities for the function $\mathfrak{e}(z)$ from which branch cuts emerge. We are not interested any further in those singularities because it can be shown that they necessarily lie in the complex lower half-plane and that the branch cuts can always be chosen to also totally lie in the complex lower half-plane. We wish to choose a contour which lie in the complex upper half-plane which will then avoid all the modes of the system of scatterers.

In addition to singularities from $\mathfrak{e}(z)$, the function $f(z)$ inherits those from the hyperbolic cotangent. Those lie on the imaginary axis at $z = in2\pi k_B T/\hbar \equiv i\xi_n$, $n$ being an integer. We call those (purely imaginary) frequencies Matsubara frequencies. The residue is

$$\mathrm{Res}(f(z), z = i\xi_n) = \frac{2k_B T}{\hbar} \mathfrak{e}(i\xi_n) \qquad (46)$$

As a contour, we choose an outer semi-circle $C_R$ of radius $R$, an inner semi-circle $C_\eta$ of radius $\eta$. Both semi-circles lie in the upper half-plane, we close our contour by parts on the real axis covering $\omega \in [-R, -\eta]$ and $\omega \in [\eta, R]$ and the whole contour is oriented counterclockwise. Lastly, we deform the contour by taking the limits $R \to \infty$ and $\eta \to 0$ and apply Cauchy's theorem. Upon deformation

* the integral over $C_R$ vanishes

* the integral over $C_\eta$ picks up the residue of $f(z)$ at $z = 0$ viewed from a $-\pi$ angle

* the result of Cauchy's theorem picks up the residues of $f(z)$ at $z = i\xi_n$, $n \ge 1$, viewed from a $+2\pi$ angle

Putting everything together, we obtain

$$\int_{-\infty}^{\infty} d\omega \coth \beta \hbar \omega/2 [\mathfrak{e}(\omega)] = 2i\pi \frac{2k_B T}{\hbar} \sum_{n=0}^{\infty} \mathfrak{e}(i\xi_n) \qquad (47)$$
We have introduced above two notations. The slashed integral reads $\int_{-\infty}^{\infty} = \lim_{\eta \to 0} \int_{-\eta}^{\eta} + \int_{\eta}^{\infty}$ and the primed sum $\sum_{n=0}^{\infty}$ means that the term with $n = 0$ is weighted with 1/2. Finally, taking advantage of the parity properties of $f(\omega)$ and if the slashed integral is equal to the regular one, we obtain a particularly simple expression for the Casimir interaction free energy:

$$A_{int}(\mathbf{L}) = k_B T \sum_{n=0}^{\infty}{}^{\prime} \mathcal{E}(i\xi_n) \iff \int_{-\infty}^{\infty} d\omega f(\omega) = \int_{-\infty}^{\infty} d\omega f(\omega) \quad (48)$$

Providing that the function $\mathcal{E}(i\xi)$ can be evaluated, this expression is now trivial to use since the series is exponentially convergent. An expression for the limit $T \to 0$ can be easily obtained by noting that $\frac{2\pi}{h} A_{int} = \frac{2\pi k_B T}{h} \sum_{n=0}^{\infty}{}^{\prime} \mathcal{E}(in\frac{2\pi k_B T}{h})$ is a Riemann sum. We then have

$$E_{int}(\mathbf{L}) = \lim_{T \to 0} A_{int}(\mathbf{L}) = \frac{\hbar}{2\pi} \int d\xi \mathcal{E}(i\xi) \quad (49)$$

which bears a striking formal resemblance with eq. (36). This is why, in numerous work and as far as $T = 0$ is concerned, this whole section on the application of Cauchy's theorem is often summed up by saying that the replacement $\omega \to i\xi$ is to be done.

## 2.2 Applications

In this section, we present a few concrete applications. The formulas presented in the previous section are basis-independent. All is left to specify a basis in which the operators are represented. This amounts to specify the quantum numbers defining a scattering channel. Obviously, this is chosen depending on which type of scatterer are considered.

In the following, we begin with the simplest application of the scattering formula which is the study of the Casimir interaction energy between two plane-parallel mirrors. This will lead to the so-called Lifshitz formula which can then be used in approximations for non-planar geometries. We continue with situations involving diffraction gratings for which a few experiments have been performed.

We start from the expression of the Casimir interaction free energy between two objects given by eqs. (33) and (39). We introduce the so-called "round-trip" operator $\rho$ as

$$\rho = S_{1}^{ii} T_{12} S_{2}^{ii} T_{21} \quad (50)$$

in terms of the objects' scattering operators $S_i$ acting among the internal channels (which we might refer to as reflection operators) and translation operator $T$ (in our case, we will always have $T_{12} = T_{21} = T$).

### 2.2.1 Plane-parallel mirrors: The Lifshitz formula

Consider now the situation of two plane-parallel mirrors. The mirrors lie in the $(xy)$ plane and are separated by a distance $L$ along the $z$ direction. This distance is the only degree of freedom and coincide with the symbolic quantity $\mathbf{L}$ entering
eq. (33). Obviously a basis adapted to this geometry consists of plane waves. Those are characterized by a wave vector **k** and a polarization **σ**. Therefore, a scattering channel is defined by the three components of the wave vector k_x, k_y and k_z and the polarization **σ**: |k_x, k_y, k_z, **σ**〉. We will not use this basis though but one very closely related to it: we wish to emphasize on the good quantum numbers which are conserved through a scattering process. Those are the norm of the wavevector k = ω/c (which correspond to the fact that the mirrors are static and therefore the reflection conserve the frequency ω) and the two components k_x and k_y (which correspond to the fact that we consider the reflection by the plane mirror to be specular). We collect the two components k_x and k_y in a vector **q** as **k** = **q** + k_z ê_z. The z-component is now deduced as k_z = ±√ω²/c² − **q**² where the sign of k_z indicates the direction of propagation. Our final basis is then |ω, **q**, ±, **σ**〉. Finally, we consider the mirrors to be made of a linear, homogeneous material whose optical properties can be described by a local relative permittivity ϵ(ω) and permeability μ(ω). Under those assumptions and the fact that the positive z axis goes from mirror ① to mirror ②, we have

$$
\mathbf{S}_1|\omega, \mathbf{q}, -, \sigma\rangle = r_1^\sigma(\omega, \mathbf{q})|\omega, \mathbf{q}, +, \sigma\rangle \quad (51)
$$

$$
\mathbf{S}_2|\omega, \mathbf{q}, +, \sigma\rangle = r_2^\sigma(\omega, \mathbf{q})|\omega, \mathbf{q}, -, \sigma\rangle \quad (52)
$$

$$
\mathbf{T}_{12}|\omega, \mathbf{q}, +, \sigma\rangle = e^{ik_z L}|\omega, \mathbf{q}, +, \sigma\rangle \quad (53)
$$

$$
\mathbf{T}_{21}|\omega, \mathbf{q}, -, \sigma\rangle = e^{ik_z L}|\omega, \mathbf{q}, -, \sigma\rangle \quad (54)
$$

$$
\rho|\omega, \mathbf{q}, +, \sigma\rangle = r_1^\sigma r_2^\sigma e^{2ik_z L} |\omega, \mathbf{q}, +, \sigma\rangle \quad (55)
$$

Above, the matrix elements of the scattering operators are Fresnel reflection amplitudes [11] for the two polarizations σ = s and σ = p which can be simply written in terms of the z-component of the wavevector in vacuum k_z and in the material K_z:

$$
r^s = \frac{\mu k_z - K_z}{\mu k_z + K_z} \qquad (56)
$$

$$
r^p = \frac{\epsilon k_z - K_z}{\epsilon k_z + K_z} \qquad (57)
$$

$$
K_z = \sqrt{\epsilon \mu \omega^2 / c^2 - q^2} \qquad (58)
$$

The plane-plane geometry is therefore a very special case where the operator **ρ** and hence also log(**1** − **ρ**) is diagonal in a basis of plane waves. The continuous variable **q** can be discretized in the usual way by imposing periodic boundary conditions on a cell of dimension 2π/L_x × 2π/L_y in the (k_x, k_y) space to obtain

$$
\mathrm{Im} \log \det(\mathbf{1} - \boldsymbol{\rho}) \xrightarrow[L_{x},L_{y} \to \infty]{\qquad} \frac{L_{x}L_{y}}{4\pi^2} \sum_{\sigma=s,p} \int d^2\mathbf{q} \operatorname{Im} \log [1 - r_1^\sigma(\omega, \mathbf{q})r_2^\sigma(\omega, \mathbf{q})e^{2ik_z L}] \qquad (59)
$$
and $L_x, L_y$ can be identified with the transverse dimensions of the plane mirrors <fnref n="1" />. We finally get

$$ \mathcal{A}(L) = \frac{\hbar}{2\pi} \int_{0}^{\infty} d\omega \coth \beta\hbar\omega/2 \left( \sum_{\sigma=s,p} \int \frac{d^2\mathbf{q}}{4\pi^2} \mathrm{Im} \log [1 - r_1^{\sigma}(\omega, \mathbf{q})r_2^{\sigma}(\omega, \mathbf{q})e^{2ik_zL}] \right) \quad (60) $$

$$ = k_B T \sum_{n=0}^{\infty} \sum_{\sigma=s,p}^{'} \int \frac{d^2 \mathbf{q}}{4\pi^2} \log [1 - r_1^{\sigma}(i\xi_n, \mathbf{q})r_2^{\sigma}(i\xi_n, \mathbf{q})e^{-2\kappa_z L}] \quad (61) $$

where now $\mathcal{A} = A_{int}/L_x L_y$ is a free energy per unit area and $\kappa_z = \sqrt{\xi_n^2/c^2 + \mathbf{q}^2}$. This is the so-called Lifshitz formula. It has been derived (for the Casimir force) by EM Lifshitz in his 1956 seminal paper [12] using a completely different approach based on fluctuating electromagnetism.

For the sake of simplicity, let us define dimensionless quantities and an energy spectral density. As a unit of frequency, we choose the free spectral range $\omega_0$ of the Fabry-Perot formed by the two plane mirrors. We also define a characteristic thermal frequency $\omega_T$ and associated thermal wavelength $\lambda_T$. Finally, we define a characteristic unit of energy $E_0$ in terms of $\hbar, c$ and $L$. We set

$$ \omega_0 = 2\pi \times \frac{c}{2L} \quad (62) $$

$$ \omega_T = \frac{2k_B T}{\hbar} = 2\pi \times \frac{c}{2\lambda_T}; \quad \lambda_T \approx 12 \, \mu\mathrm{m} \, @ \, T = 300 \, \mathrm{K} \quad (63) $$

$$ E_0 = \frac{\hbar c}{8\pi^2 L^3} \quad (64) $$

$$ \mathfrak{e}^{\sigma}(\omega) = L^2 \int d^2\mathbf{q} \log [1 - r_1^{\sigma}(\omega, \mathbf{q})r_2^{\sigma}(\omega, \mathbf{q})e^{2ik_zL}] \quad (65) $$

and we get

$$ \frac{\mathcal{A}(L)}{E_0} = \int_0^\infty d\left(\frac{\omega}{\omega_0}\right) \coth\left(\frac{\lambda_T}{L}\frac{\omega}{\omega_0}\right) \mathrm{Im}\left[\mathfrak{e}^s\left(\frac{\omega}{\omega_0}\right) + \mathfrak{e}^p\left(\frac{\omega}{\omega_0}\right)\right] \quad (66) $$

$$ = \frac{\xi_1}{\omega_0} \sum_{n=0}^{\infty} \left[ \mathfrak{e}^s \left( i \frac{\xi_n}{\omega_0} \right) + \mathfrak{e}^p \left( i \frac{\xi_n}{\omega_0} \right) \right] \quad (67) $$

In the last equation, $\xi_n = n\pi \frac{L}{\lambda_T}$ where $\xi_n$ are the Matsubara frequencies. Note that, in first approximation <fnref n="2" />, the spectral densities $\mathfrak{e}^{\sigma}(\omega)$ do not depend on the temperature but only on the geometry. The only effect of the temperature is through the hyperbolic cotangent weighting function which will enhance the effect of the low

<footnote n="1">¹This point actually deserves a discussion: if the plane mirrors have infinite transverse dimension, tr log(1 - $\rho$) is actually infinite because it is an extensive quantity. In this case, the operator log(1 - $\rho$) is not trace class but is anyway ill-defined because there cannot be an "outside from the two mirrors" from which incoming fluctuations are scattered to. The solution is to divide by the size of the mirrors to get an intensive quantity.</footnote>

<footnote n="2">²Rigorously, the optical properties of the scatterers depend of the temperature. Here, we neglect this effect.</footnote>
frequencies. We conclude with an interesting formula: by using the series representation of the hyperbolic cotangent, we can show that

$$ \coth \left( \frac{\lambda_T}{L} \frac{\omega}{\omega_0} \right) = \omega_0 \frac{L}{\lambda_T} \sum_{n=0}^{\infty} \frac{2\omega}{\omega^2 + \xi_n^2} \quad (68) $$

$$ \mathrm{Re} [\mathfrak{e}^{\sigma}(i\xi)] = \mathfrak{e}^{\sigma}(i\xi) = \frac{1}{\pi} \int_{0}^{\infty} d\omega \frac{2\omega}{\omega^2 + \xi^2} \mathrm{Im} [\mathfrak{e}^{\sigma}(\omega)] \quad (69) $$

The above formula is interesting because it is a Kramers-Kronig relation: the spectral density $\mathfrak{e}^{\sigma}(\omega)$ is actually a causal linear response function which obeys Kramers-Kronig relation. Its real and imaginary parts are therefore related through a Hilbert transform. It also emphasizes the fact that the value $\mathfrak{e}^{\sigma}(i\xi_n)$ for a single Matsubara frequency is always an integral over all frequencies of the spectral density. This integral is weighted by the function $\frac{2\omega}{\omega^2+\xi_n^2}$ which has a maximum at $\omega = \xi_n$.

### 2.2.2 Drude-plasma models for metallic mirrors

In order to illustrate an important point concerning the comparison between theory and experiments, let us now express the preceding formula when applied to a metallic plane mirror. To that end, we have to specify a form for the relative permittivity $\epsilon(\omega)$ (we consider here non-magnetic metals i.e. we take $\mu(\omega) = 1$). The complex refractive index $N = n + ik$ of usual materials is tabulated for a certain frequency range and the associated relative permittivity is given by $\epsilon_{tab}(\omega) = N^2(\omega)$. Those tabulated data have to be extrapolated both at low and high frequencies. While the extrapolation to high frequencies do not pose a problem, there is in the Casimir physics community two different approaches for the extrapolation to low frequencies.

A reasonable way of extrapolating the relative permittivity to low frequencies is to use a Drude model

$$ \epsilon_{\text{Drude}}(\omega) = 1 - \frac{\omega_p^2}{\omega(\omega + i\gamma)}. \quad (70) $$

While rather crude, this model reproduces the static conductivity $\frac{\sigma(0)}{\epsilon_0} = \frac{\omega_p^2}{\gamma}$ for suitable values of the plasma frequency $\omega_p$ and the relaxation frequency $\gamma$. Another way is to neglect the conduction electrons relaxation properties and to use a so-called plasma model

$$ \epsilon_{\text{plasma}}(\omega) = 1 - \frac{\omega_p^2}{\omega^2}. \quad (71) $$

This model obviously do not reproduce the properties of metals at low frequencies and one could wonder why use it in the first place. It turns out that the most precise Casimir experiments are in agreement with theoretical calculations using the plasma model, excluding calculations using the more physically accurate Drude model. This state of affair has been dubbed the "Drude-plasma puzzle" and has not yet been solved. The issue appeared back in 2007 [13] and has again been confirmed very recently [14].

To gain insight on the differences in using the Drude or the plasma models, let us illustrate the two important preceding formulas: the Lifshitz formula eq. (66) and the Lifshitz-Matsubara formula eq. (67) for a Drude relative permittivity. We take
a Drude permittivity at all frequencies for the sake of clarity because we are only interested into what happens at low frequencies. We take numerical values $\hbar\omega_p = 9$ eV and $\hbar\gamma = 35$ meV which are typical of good metals (those particular values are often used to model gold). As for the other parameters, we take $T = 300$ K and a distance $L = 500$ nm (the most precise experiments are done at ambient temperature for distances ranging from 300 to 800 nm). The Lifshitz formula associated with the particular plane-plane geometry lends itself into an analysis in terms of contributions from $s$ or $p$ polarizations as well as from propagative ($\omega^2/c^2 > |\mathbf{q}|$) or evanescent ($\omega^2/c^2 < |\mathbf{q}|$) waves.

We show in figures (2) and (3) the illustration of the Lifshitz formula (66) for the parameters cited above. Every features in those figures are associated with the

<figure><img src="image_4.png" /><figcaption>Figure 2: Illustration of the Lifshitz formula (66) for the $s$ polarization. The Casimir free energy is the integral of the functions shown.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Illustration of the Lifshitz formula (66) for the $p$ polarization. The Casimir free energy is the integral of the functions shown.</figcaption></figure>

different modes pertaining to the geometry of two plane-parallel metallic surfaces :

* the modes associated with the Fabry-Perot cavity formed by the two mirrors,

* the modes associated with surface plasmons and

* the modes associated with magnetic eddy currents.
Fabry-Perot modes are propagative in nature and occur for both *s* and *p* polarizations in the frequency range $0 < \omega \lesssim \omega_p$. They show up in both figures (2) and (3) as excess in free energy density near the frequencies $\omega \approx n\omega_0$ together with a default in free energy density between those frequencies. With the Drude model for the permittivity which we use, the plasma frequency $\omega_p$ acts as a rather sharp frequency cut-off for those modes. With a more realistic model based on tabulated optical data, the vanishing of those modes would be smoother. The net effect (excess or default) of those modes when integrated over all frequencies is not obvious. It turns out though that the net effect is always a default in free energy density. Surface plasmons modes occur for evanescent *p* polarization in a frequency range $0 < \omega \lesssim \omega_p/\sqrt{2}$ and therefore show up in figure (3). As already studied in ref. [15], the two plasmonic branches lead to both an excess and default in the free energy density, depending on the frequency. Finally, the magnetic eddy currents modes occur for evanescent *s* polarization in a frequency range $0 < \omega \lesssim \gamma$. They show up as the very sharp peak near $\omega = 0$ in figure (2) and they contribute to an excess of free energy density.

We now illustrate the Lifshitz-Matsubara formula (67) in figure (4). For the sake

<figure><img src="image_6.png" /><figcaption>Figure 4: Illustration of the Lifshitz-Matsubara formula (67). Note that the sign of the *s* polarization contribution has been changed for clarity. The Casimir free energy is the area of the rectangles.</figcaption></figure>

of clarity, we show the functions $-\epsilon^s(i\frac{\xi_n}{\omega_0})$ and $\epsilon^p(i\frac{\xi_n}{\omega_0})$. The Casimir free energy is the sum of the area of the rectangles in figure (4) and is obviously a negative quantity.

Now that we have some intuition about every feature in the free energy spectral density $\epsilon^\sigma(\omega)$ we can study what happens in the limiting case where the conduction electron relaxation parameter $\gamma$ vanishes; which will be the situation when using the plasma model to describe the low frequency part of the relative permittivity $\epsilon(\omega)$. In the limit $\gamma \to 0$, nothing much happens concerning the contribution to the free energy originating from the Fabry-Perot modes and the surface plasmons modes. This is not surprising as those are loosely dependent on $\gamma$ only and the Drude parameters we have chosen already correspond to a good metal having $\gamma/\omega_p \ll 1$. The situation is less clear for the contribution from the eddy currents which are restricted to the frequency range $0 < \omega \lesssim \gamma$. It turns out that the sharp peak
corresponding to the contribution of the eddy currents in fig. (2) becomes a dirac delta function $\delta(\omega)$ in the limit of $\gamma \to 0$. Therefore, the contribution of the eddy currents becomes a constant independent of $\gamma$ in the limit $\gamma \to 0$. We stress that this is a direct consequence of the finite temperature $T$ (for $T = 0$, the contribution of the eddy currents would vanish as $\gamma \to 0$).

Therefore, as far as the Lifshitz formula is concerned, the calculated free energy is continuous and smooth in the limit $\gamma \to 0$. The situation is different when using the Lifshitz-Matsubara expression eq. (67). Remember that this expression corresponds to the slashed integral as defined around eq. (47) which amounts to remove from the frequency integral the single point $\omega = 0$. As such, the Lifshitz-Matsubara expression eq. (67) fails at taking into account the contribution from the magnetic eddy currents as it becomes a Dirac delta function $\delta(\omega)$ as $\gamma \to 0$. The Casimir free energy for metallic mirrors when calculated with the Lifshitz-Matsubara formula is discontinuous at $\gamma = 0$. Going back to the "Drude-plasma puzzle", every theoretical calculation is always performed using the Lifshitz-Matsubara formula. Since it is a calculation with $\gamma = 0$ which is in agreement with the experiments and we have just shown that this calculation do not take into account the interaction between magnetic eddy currents, we have to conclude that, in the experiments, the interaction between eddy currents is somehow suppressed.

Recently, it has been recognized that by using a more sophisticated model for the permittivity of metals which include nonlocal effects can lead to a good agreement between theory and experiments [16]. However, this is still merely a phenomenological explanation rather than an *ab initio* one.

*Interpretation* - The function $\mathfrak{e}(\omega, L) = \log \det(\mathbf{1} - \rho)$ does not have an obvious physical interpretation. On the one hand, as we have seen in eq. (32) its derivative with respect to the frequency $\Delta D(\omega) \propto -\partial_\omega \mathrm{Im} \mathfrak{e}(\omega, L)$ represents the modification in electromagnetic local density of states, with respect to the vacuum DOS, induced by the mere presence of the scatterers. Furthermore, remember that the fact that we have subtracted the situation pertaining with each individual scatterer (see eq. (12)) means that we are really dealing with the modification of LDOS concentrated *between the scatterers* where the LDOS is modified by the mutual interaction between the scatterers <fnref n="3" />.

On the other hand, the derivative with respect to the distance $p(\omega) \propto -\partial_L \mathrm{Im} \mathfrak{e}(\omega, L)$ is the quantity which enters the formula for the Casimir force, or in this case pressure (see eq. (34)). This has a trivial interpretation in that when this quantity is positive, resp. negative, its contribution to the pressure is repulsive, resp. attractive.

We show in Figure (5) the two quantities $\Delta D(\omega)$ and $p(\omega)$ corresponding to the parameters of the preceding example. Those two quantities are seen to be nearly identical and lends credence to the following interpretation of the Casimir interaction: for frequencies where there is less DOS between the scatterers than on the outside (where $\Delta D(\omega) < 0$), the contribution from those frequencies is attractive. And vice versa. The vacuum outside the scatterers exerts an inward pressure while the DOS located between the scatterers exerts an outward pressure. The net pres-

<footnote n="3">³The concept of "between the scatterers" can be a bit fuzzy when dealing, for instance, with the interaction between two spheres. Here, in the plane-plane geometry it is quite clear</footnote>
sure (attractive or repulsive) depends whether there is more or less DOS between the scatterers than outside the scatterers. This interpretation of the Casimir interaction

<figure><img src="image_7.png" /><figcaption>Figure 5: Derivative with respect to frequency and distance of the quantity $e(\omega, L)$ defined in (62). Positive parts correspond to excess in density of states between the plates and a repulsive force and vice versa.</figcaption></figure>

is often used and is therefore seen to be quite accurate.

*Thermodynamics of the Casimir interaction* - Recall that the Casimir Helmholtz free energy $A$, the Casimir internal energy $U$ and the entropy $S$ are linked by the thermodynamical equation $A = U - TS$. Furthermore, since the entropy is defined as $S = -\partial_T A$ we obtain the additional equation $\partial_T U = T\partial_T S$ which would define the heat capacity. We define a reduced temperature $\mathcal{T} = \pi L/\lambda_T \approx 837 (L/1\text{m}) (T/1\text{K})$. We show in figs (6) and (7) the free energy $A$, the internal energy $U$ and the entropy $S$ as a function of the reduced temperature $\mathcal{T}$ for a system of two gold plates and for a system of gold-silica plates. The plates are separated by a distance $L = 4~\mu\text{m}$ leading to realistic absolute temperature.

In the case of the system of two gold plates, a peculiar feature appears. The Casimir entropy is negative for a range of temperature. This is obviously related to the fact that the free energy begins to *decrease in magnitude as the temperature increase*. Since $S(0) = 0$ according to Nernst heat theorem, the entropy must necessarily be negative in a certain range of temperature. The fact that the free energy begins to decrease in magnitude as the temperature increase is due to the interaction between eddy currents as we have seen already before. More precisely, it is due to the fact that this interaction is at very low frequency and is repulsive. On the other hand, in the case of the interaction between a gold plate and a silica one there is no emergence of negative entropy because the dielectric plate does not support eddy currents. This can be seen in figure (7).

Other more universal features which can be seen on both figures (6) and (7) are the fact that the entropy reaches an asymptotic value $S_\infty$ and the fact that the internal energy $U$ vanishes at large temperature. Therefore, the free energy reaches an asymptotic regime $A(T) \to -TS_\infty$ where it increases linearly with temperature as $T \to \infty$.

We conclude this part on the thermodynamics of the Casimir interaction by emphasizing that negative entropies in this context are not necessarily surprising.
<figure><img src="image_8.png" /><figcaption>Figure 6: Thermodynamics of the Casimir interaction between two gold plates separated by a distance $L = 4$ µm.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Thermodynamics of the Casimir interaction between a gold and a silica plates separated by a distance $L = 4$ µm.</figcaption></figure>
Indeed, we stress once again that the quantity we calculate here is an interaction entropy of the electromagnetic field. As such, it is the difference in entropy in the field with or without the presence of the objects. Most notably, it does not take into account any entropy in the objects themselves. If one were to take into account both the entropies in the field and in the objects, for a given geometric arrangement one would of course find a positive quantity. Another note is that the figures (6) and (7) are done by neglecting any dependence of the permittivity with the temperature. One can model the temperature dependence of the permittivity to a certain extent by, for instance for metals, relating the Drude dissipation parameter $\gamma$ to the resistivity $\rho$ and use tabulated values for $\rho(T)$. Doing so only change marginally the results shown in figures (6) and (7).

### 2.2.3 Grating geometry

In this section, we consider the Casimir interaction between nanostructured surfaces. In particular, we consider surfaces with nanoscale rectangular corrugations as a number of experiments have been performed using that kind of surfaces.

As with the previous section, we start from the basis-independent expression for the Casimir free interaction energy between two objects (eqs. (33) and (39))

$$A_{int}(\mathbf{L}) = \frac{\hbar}{2\pi} \int d\omega \coth \beta\hbar\omega/2 [\mathrm{Im} \log \det(\mathbf{1} - \mathbf{S}_1 \mathbf{T}_{12} \mathbf{S}_2 \mathbf{T}_{21})] \quad (72)$$

where the scattering operators $\mathbf{S}_1$ and $\mathbf{S}_2$ are understood to be the part acting on the internal channels. We consider the two surfaces to lie in the $(xy)$ plane so that the distance between them is along the $z$ direction i.e. $\mathbf{L} = L\hat{e}_z$.

If we consider the nanostructures to consist of corrugation lines along the $y$ direction, of period $d$, the surface will diffract an incoming plane wave with wavevector $\mathbf{k}_{in} = (k_x, k_y, k_z)$ into an infinite number of outgoing plane waves with wavevectors $\mathbf{k}_{out} = (k_x^n, k_y^n, k_z^n)$ with $k_x^n = k_x + n\frac{2\pi}{d}$, $k_z^n = \pm\sqrt{\omega^2/c^2 - (k_x^n)^2 - k_y^n}$ and $n$ is the diffraction order.

By analogy with the previous section, our final basis consists of all the different plane waves coupled by the diffraction process i.e. $|\omega, k_x^n, k_y, \pm, \sigma\rangle \equiv |\omega, \mathbf{q}, n, \pm, \sigma\rangle$ with $\mathbf{q} = (k_x, k_y)$, $\pm$ represent the direction of propagation and $\sigma$ is the polarization.

We are left to evaluate the matrix element of the scattering operators $\mathbf{S}_i$ as well as of the translation operators $\mathbf{T}$ entering the eq. (72). At this point, we stress an important advantage of using the scattering method: by expressing the Casimir interaction energy solely in terms of scattering operators, we can take advantage of the fact that those operators are usually known from other fields of research than Casimir physics. Here, obviously the expression of the scattering operator of a diffraction grating is known from classical physics and we can readily use the expressions previously derived. We won't go into the details of the derivation which is presented in [17] and heavily inspired from an earlier classical derivation [18].

The scattering operators are, in practice, truncated to a maximum number of diffraction orders $N$. In that case, the operators are of dimension $2 \times (2N + 1)$ as the diffraction order satisfy $-N < n < N$ and the factor two accounts for the polarization degree of freedom. Reflection and transmission operators are found
by solving a system of 2 × (2N + 1) coupled differential equations satisfied by the reflected and transmitted fields. Translation from one surface to the other does not diffract nor mixes polarizations ; the translation operators are therefore diagonal :

$$
\mathbf{S}_1|\omega, \mathbf{q}, n, -, \sigma\rangle = \sum_{n'\sigma'} [\mathbf{S}_1]^{\sigma\sigma'}_{nn'} |\omega, \mathbf{q}, n', +, \sigma'\rangle \quad (73)
$$

$$
\mathbf{S}_2|\omega, \mathbf{q}, n, +, \sigma\rangle = \sum_{n'\sigma'} [\mathbf{S}_2]^{\sigma\sigma'}_{nn'} |\omega, \mathbf{q}, n', -, \sigma'\rangle \quad (74)
$$

$$
\mathbf{T}_{12}|\omega, \mathbf{q}, n, +, \sigma\rangle = e^{ik_z^n L}|\omega, \mathbf{q}, n, +, \sigma\rangle \quad (75)
$$

$$
\mathbf{T}_{21}|\omega, \mathbf{q}, n, -, \sigma\rangle = e^{ik_z^n L}|\omega, \mathbf{q}, n, -, \sigma\rangle \quad (76)
$$

As in the previous section concerning the parallel plane geometry, we need to write out the trace appearing in the expression of the Casimir interaction energy (72).<fnref n="4" /> The situation is very similar to the parallel planes case except that the translational invariance has been lost because of the corrugations. The periodicity along the *x* direction means that the integration over **q** has to be restricted to the first Brillouin zone (FBZ). The first Brillouin zone restricts *k*<sub>x</sub> to an interval of length *2π/d* while the range for *k*<sub>y</sub> remains unrestricted

$$
\int_{\text{FBZ}} d^2\mathbf{q} \equiv \int_{-\pi/d}^{\pi/d} dk_x \int_{-\infty}^{\infty} dk_y \quad (77)
$$

We arrive at a final expression for the Casimir interaction energy between cor-
rugated surfaces

$$
\mathcal{A}(L) = \frac{\hbar}{2\pi} \int_{0}^{\infty} d\omega \coth \beta\hbar\omega/2 \int_{\text{FBZ}} \frac{d^2\mathbf{q}}{4\pi^2} \text{Im} \log \det(1 - \mathbf{S}_1 \mathbf{T}_{12} \mathbf{S}_2 \mathbf{T}_{21}) \quad (78)
$$

where once again $\mathcal{A}$ is an interaction energy per unit surface.

Thanks to that expression of the Casimir interaction energy, we are able to com-
pare with a series of experiments involving corrugated surfaces. Those experiments
were performed with silicon surfaces with either deep or shallow corrugations [19, 20].
We show in figure (8) scanning electron microscope images of the corrugated sur-
faces used in the experiments. They consist of corrugations with a 400 nm period
and the depth of the corrugations is either 90 nm or 1 µm. Those surfaces are made
out of doped silicon and the quality of the corrugations is very good as seen on the
SEM images.

The bottom-right part of figure (8) shows a schematics of the experiment. It
consists of two gold covered spheres of radius R = 50 µm glued on a micromechani-
cal resonator consisting of a 500 µm square silicon plate suspended by two torsional
rods. The experiments measures the interaction between the gold-covered sphere
and the corrugated surface which is brought above the micromechanical resonator.
The presence of the corrugated surface induces a shift in the frequency of the mi-
cromechanical resonator which is proportional to the force gradient F'_{sc}(L) between
the sphere and the corrugated surface. The quantity measured in the experiments is
therefore the gradient of the Casimir force whose expression is F'_{sc}(L) = -∂F_{sc}/∂L

<footnote n="4">⁴Remember that for a trace-class operator, $\log \det \equiv \text{tr} \log$.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 8: Scanning electron microscope images of the corrugated silicon surfaces used in the experiments. Top row: side (a) and top (b) views of the deep corrugations surface. Bottom row: side (a) view of the shallow corrugations surface and (b) schematics of the experiment.</figcaption></figure>

with $F_{sc}(L)$ being the Casimir force between the gold-coated sphere and the corrugated surface. Finally, the Casimir force between the gold-coated sphere and the corrugated surface can be related to our previous expression (78) using the so-called proximity force approximation (PFA) which is valid when the distance between the sphere and the surface is small compared to the radius of the sphere i.e. $L/R \ll 1$. Under the proximity force approximation, we have $F_{sc}(L) \approx 2\pi R\mathcal{A}(L)$.

Finally, the force gradient $F'_{sc}(L)$ between the sphere and the corrugated surface can be normalized by its expression under the pairwise additive approximation (PAA) to absorb most of the variation with respect to the distance $L$. Under this approximation, the interaction Casimir energy $\mathcal{A}(L)$ is given by the sum of the contributions from the top of the corrugations and from the bottom. It therefore reads $\mathcal{A}_{PAA}(L) = 0.5 [\mathcal{A}_{pp}(L) + \mathcal{A}_{pp}(L+a)]$ where $\mathcal{A}_{pp}$ is the Casimir interaction energy between plane parallel surfaces given by (60), the height of the corrugation is given by $a$ and the factor 0.5 represents the duty cycle of the corrugations. The ratio $\rho = \frac{F'_{sc}(L)}{F'_{PAA}(L)}$ is therefore a smooth quantity of order one for all measured distances $L$.

We show in figure (9) the ratio $\rho$ measured experimentally together with theoretical calculations for the experiments involving either deep or shallow corrugations. As predicted, the ratio $\rho$ stays of the order unity throughout the range of distances explored. This ratio measures in a sense the degree of “non additivity” of the Casimir interactions which is seen to reach around 30 % at most. The theoretical calculations are done both for $T = 0$ K and $T = 300$ K, which is the temperature at which the experiments are done. The agreement between theory and experiments is reasonably good while the experimental uncertainties cannot highlight any clear effect due to the temperature.
<figure><img src="image_11.png" /><figcaption>Figure 9: Comparison with theory for the experiments measuring the Casimir force gradient between a gold-coated sphere and a silicon surface with (a) deep and (b) shallow corrugations. The quantity shown is the ratio $\rho = F'_{sc}(L)/F'^{PAA}_{sc}(L)$ as a function of the distance $L$ (see text). Experimental data are the points with the error bars. Theoretical calculations are shown as solid curves and are done for $T = 0$ K (blue) and $T = 300$ K (red).</figcaption></figure>

### 2.2.4 Quantum reflection and gravitational quantum states

In this part, we are interested in the phenomenon of quantum reflection. When a slow quantum particle of energy $E$ is incident onto a material surface, it has a probability of being reflected by the Casimir-Polder interaction potential even though this interaction potential is attractive. The probability of reflection tends to unity as the particle's momentum $\hbar k = \sqrt{2mE}$ tends to zero. Quantum reflection has been observed experimentally using slow metastable helium or neon atoms [21, 22] incident upon different surfaces.

With the possibility of reaching nearly perfect quantum reflection, we then study the problem of gravitational quantum states which are quantum states of a particle above a material surface in a gravitational field. While the energies of those states are analytically known in the case of perfect reflection, we study the shifts of those energies when taking into account the Casimir-Polder potential.

*Quantum reflection* — We start by expressing the Casimir-Polder interaction
energy between an atom and a material surface. It is obtained starting from the general expression for the Casimir free interaction energy between two objects given by equation (72). From this expression, we will make a number of approximation :

* First of all, we will consider the interaction energy at zero temperature $E_{int}$ rather than the free energy $A_{int}$. Temperature has mostly the effect of populating excited states of the atom. For an atom in its ground state, which we consider here, a temperature of 300 K leads to a negligible population in the excited states which lies at optical frequencies from the ground state.

* Second of all, we will consider that the reflection operator for the atom is small. This assumes that an atom is a very poor scatterer and will be made clear when we give the expression of the atomic reflection operator. Under this assumption, we will then develop the logarithm appearing in eq.(72) to first order.

Under those two approximations, our starting expression for the Casimir-Polder interaction energy is

$$E_{int}^{CP}(\mathbf{L}) = -\frac{\hbar}{2\pi} \int d\omega \operatorname{Im} \operatorname{tr} (\mathbf{S}_1 \mathbf{T}_{12} \mathbf{S}_2 \mathbf{T}_{21}) \quad (79)$$

In a basis of plane waves $|\omega, \mathbf{q}, \pm, \sigma\rangle$, the reflection operator for the atom is [23]

$$\mathbf{S}_1|\omega, \mathbf{q}', -\sigma'\rangle = \int \frac{d^2\mathbf{q}}{4\pi^2} \sum_{\sigma} \left[ \frac{i\omega^2\alpha(\omega)}{2c^2k_z} \hat{e}_{\sigma}^{-}(\mathbf{q}) \cdot \hat{e}_{\sigma'}^{+}(\mathbf{q}') \right] |\omega, \mathbf{q}, +, \sigma\rangle \quad (80)$$

where $\hat{e}_{\sigma}^{\pm}(\mathbf{q})$ are the unitary vectors for a plane wave $|\omega, \mathbf{q}, \pm, \sigma\rangle$ and $\alpha(\omega)$ is the atomic dynamic dipole polarizability<fnref n="5" />. The expressions for the translation operators $\mathbf{T}$ and the reflection operator for the planar surface $\mathbf{S}_2$ in terms of Fresnel amplitudes $r^s$ and $r^p$ have been given in eqs. (51).

Putting everything together, we finally obtain

$$E_{int}^{CP}(z) = -\frac{\hbar}{2\pi} \int d\omega \operatorname{Im} \frac{i\omega^2\alpha(\omega)}{c^2} \int \frac{d^2\mathbf{q}}{4\pi^2} \frac{e^{2ik_z z}}{2k_z} \left[ r^s - \left(1 - \frac{2c^2q^2}{\omega^2}\right) r^p \right] \quad (81)$$

with $\mathbf{L} = z\hat{e}_z$. The Casimir-Polder interaction energy given by the above expression has two limiting behavior for small and large distances $z$ i.e.

$$E_{int}^{CP}(z) \underset{z/\lambda \ll 1}{\longrightarrow} -\frac{C_3}{z^3} \quad (82)$$

$$E_{int}^{CP}(z) \underset{z/\lambda \gg 1}{\longrightarrow} -\frac{C_4}{z^4} \quad (83)$$

where the change of behavior occurs around $\lambda$ which is a typical atomic transition wavelength appearing in the dynamical polarizability $\alpha(\omega)$.

As an example, we show in Figure (10) the Casimir-Polder potential between an hydrogen atom in its ground state and a silicon surface. In order to have an atom-

<footnote n="5">⁵The polarizability is $\alpha = \alpha_{SI}/4\pi\epsilon_0$ and has the dimension of a volume.</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 10: Casimir-Polder interaction potential $E_{int}^{CP}(z)$ between an hydrogen atom in its ground state and a silicon surface as a function of the atom-surface distance $z$. The quantity shown is $z^3|E_{int}^{CP}|$ and highlights the non-retarded behavior in $-C_3/z^3$ and the retarded behavior in $-C_4/z^4$.</figcaption></figure>

surface interaction potential at all distances, the Casimir-Polder potential must be supplemented with a short range part responsible for hydrogen atoms adsorption on the silicon surface. In the case of hydrogen on a silicon surface, some data exist in the literature. We use the data found in ref. [24] and fit a Morse potential. We obtain an atom-surface potential $V(z)$ valid at all distances that we show in Figure (11).

Now that we have the atom-surface at all distance $V(z)$, we can calculate the probability of quantum reflection. This is done by solving the 1D time-independent Schrödinger equation with appropriate initial conditions. The Schrödinger equation is

$$ -\frac{\hbar^2}{2m} \Psi''(z) + V(z) \Psi(z) = E \Psi(z) \quad (84) $$

with $E = \frac{\hbar^2 k^2}{2m} = \frac{1}{2}mv^2$ is the atom incident kinetic energy on the surface ($v$ is the component of the atom's velocity normal to the surface). This equation is solved by imposing that the wave function is an incoming wave at a $z = z_s$ close to the surface i.e. $\Psi(z_s) = \varphi^{in}(z_s)$. The wave function is then analyzed at a distance $z = z_\infty$ far from the surface and decomposed over incoming and outgoing solutions i.e. $\Psi(z_\infty) = A\varphi^{in}(z_\infty) + B\varphi^{out}(z_\infty)$. The probability of quantum reflection is then $R = |r|^2 = |B/A|^2$. The value of the distance $z_s$ close to the surface is chosen to be around 1 nm; the final value for the probability of quantum reflection should not depend on the exact value of $z_s$. As for the value of $z_\infty$, it must be chosen at a distance where the effect of the potential are negligible i.e. $E \gg V(z_\infty)$. Quantum reflection occurs around a distance $z_r$ when the local wave vector $k(z) = \hbar^{-1}\sqrt{2m(E - V(z))}$ changes behavior i.e. $|V(z_r)| = E$. We show in Figure (12) the probability of quantum reflection for an hydrogen atom in its ground state impinging on a silicon surface as a function of the incident velocity. We observe that the reflectivity goes to unity as the incident velocity vanishes and that it drops rapidly for higher velocity.

**Gravitational quantum states** — In this section, we introduce the gravitational quantum states (GQS) and study the effects of the Casimir-Polder potential on their energies. GQS are the bound states of a sufficiently cold particle above a material
<figure><img src="image_13.png" /><figcaption>Figure 11: Complete atom-surface potential V(z) for an hydrogen atom and a silicon surface. The potential consists of a short range well modeled by a Morse potential responsible for adsorption as well as the Casimir-Polder part which itself can be separated into a non-retarded ($-C_3/z^3$) and retarded part ($-C_4/z^4$).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12: Reflectivity $|r|^2$ for an hydrogen atom impinging on a silicon surface as a function of the incident velocity.</figcaption></figure>

surface and in a gravitational field. Bound states occur because the particle is bounded from below by quantum reflection and from above by the gravitational field.

Let us first introduce the energies and wave functions in the case of perfect reflection on the material surface. In this case, the problem is to solve the 1D time independent Schrödinger equation given by equation (84) with the potential $V(z) = mgz$ and subjected to the boundary condition $\Psi(0) = 0$ corresponding to a condition of perfect reflection on the material surface. The equation (84) can be re-written by introducing adimensional lengths $\tilde{z} = z/z_g$ and energies $\tilde{E} = E/\varepsilon_g$ as
follows

$$ -\Psi''(\tilde{z}) + \tilde{z}\Psi(\tilde{z}) = \tilde{E}\Psi(\tilde{z}) \qquad (85) $$

$$ z_g^3 = \frac{\hbar^2}{2m^2g} \qquad (86) $$

$$ \varepsilon_g = mgz_g \qquad (87) $$

Under this form, the solutions to the equation (85) subjected to the boundary conditions $\Psi(0) = 0$ (perfect reflection) and $\Psi(\infty) = 0$ (bound state) are proportional to Airy functions $\Psi_n(\tilde{z}) \propto \text{Ai}(\tilde{z} - \tilde{E}_n)$ with the energies $\tilde{E}_n = -\lambda_n$ where the $\lambda_n$ are the zeros of the regular Airy function $\text{Ai}(z)^6$.

To take into account the effect of the Casimir-Polder potential, we include it in the potential entering the Schrödinger equation i.e. $V(z) = E_{\text{int}}^{\text{CP}}(z) + mgz$ and solve for the energies which we write $\mathcal{E}_n = -\lambda_n\varepsilon_g + \Delta E$. The Schrödinger equation including the Casimir-Polder potential is solved in a manner similar to the problem of quantum reflection. The initial condition consists of imposing the wavefunction to be an incoming wave close to the surface at $z = z_s$ and solving for the energies which lead to a bound state i.e. $\Psi(\infty) = 0$. The energies compatible with those boundary conditions are complex valued and express the fact that because of the imperfect reflection, the GQS are merely quasistationary states with a finite lifetime above the surface.

<figure><img src="image_15.png" /><figcaption>Figure 13: Complex energy shifts $\Delta E$ for the first 100 GQS of an hydrogen atom above a silicon surface.</figcaption></figure>

We show in Figure (13) the complex energy shifts $\Delta E$ which arise when we take into account the effect of the Casimir-Polder potential between the atom and the

<footnote n="6">⁶Note that the zeros of the regular Airy function $\lambda_n$ are negative quantities leading to positive GQS energies.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 14: Probability of quantum reflection for an incident energy equal to that of the first GQS.</figcaption></figure>

surface. We see that those energy shifts are on the level of $|\Delta E/\lambda_n| \approx 10^{-4}$. This can be understood because the effect of the Casimir-Polder potential is only at short distances to the surface. Indeed, for an hydrogen atom above an silicon surface, we find that the altitude $z_\times$ for which $mgz_\times = |E_{int}^{CP}(z_\times)|$ is $z_\times \approx 0.1z_g$. This is to be compared to the typical size of a perfect GQS $\Psi_n(\tilde{z})$ which is $-\lambda_n$ (we recall that $-\lambda_1 \approx 2.3$ and $-\lambda_{100} \approx 60.4$).
## 2.3 Miscellaneous

In this section, we quickly present a few other results related to Casimir physics based on a particular publication.

### 2.3.1 Van der Waals interaction between two different atoms

In ref. [25], included hereafter, we have studied the van der Waals interaction between two different atoms, one of which is in an excited state.

This situation, at first glance a textbook example of QED calculations, has led to a long standing controversy regarding the behavior of the interaction as a function of the distance between the two atoms. Indeed, some authors predict a monotonic behavior of the interaction energy as a function of the inter-atomic distance whereas other authors predict an oscillating behavior. It has been argued that the different behaviors found from different authors stem from an ambiguity in the treatment of poles arising from energy denominators in the standard QED treatment of the problem.

We have studied this problem using time-dependent perturbation theory as there is no such ambiguities there: the poles arising from energy denominators have to be handled in a way as to make the relevant integrals converge. Using this method, we have found an interaction energy which present an oscillatory behavior as a function of the distance between the two atoms.

We mention also that the aforementioned ambiguities in treating energy denominators can be related to the use of either advanced, retarded or Feynman propagators. A similar, oscillating, behavior of the van der Waals interaction energy has been found in ref. [26] to be a direct consequence of the use of Feynman propagators.

A conclusion to that controversy seems to have been found recently. It involves reflecting on what exactly is the van der Waals interaction energy in that situation. It can be identified as the shift, as a function of the inter-atomic distance, of the ground state atom, of the excited atom or a phase shift rate of the two atoms wavefunction [27, 28].

The conclusion is quite interesting. The shift on the ground state atom is found to be monotonic as a function of the distance between the two atoms. This was, to some extent, the quantity calculated by the authors who found a monotonic van der Waals interaction. On the other hand, the shift of the excited atom is found to have an oscillatory behavior as a function of the distance and this was the quantity calculated by the other group of authors. The controversy lasted so long because both parties were right.

This surprising results raise some questions because it apparently violates the principle of action-reaction or the conservation of momentum. The explanation is as follows: during the time the atom stays in its excited state, there is a transfer of momentum between the two atoms system and the electromagnetic field. This transferred momentum is then released when the excited atom decays in its ground state via spontaneous emission.
Quasiresonant van der Waals Interaction between Nonidentical Atoms

M. Donaire,\* R. Guérout, and A. Lambrecht

Laboratoire Kastler Brossel, UPMC-Sorbonne Universités, CNRS, ENS-PSL Research University,
Collège de France, 4, place Jussieu, F-75252 Paris, France

(Received 24 March 2015; published 15 July 2015)

We present a time-dependent quantum calculation of the van der Waals interaction between a pair of dissimilar atoms, one of which is initially excited while the other one is in its ground state. For small detuning, the interaction is predominantly mediated at all distances by the exchange of doubly resonant photons between the two atoms. We find that it presents both temporal and spatial oscillations. Spatially oscillating terms depend on the resonant frequencies of both atoms, while the frequency of the time oscillations is given by their detuning. We analyze the physical content of our findings and discuss to what extent previous conflicting stationary approaches provide compatible results. A proper account of causality is found essential in order to obtain the correct result.

DOI: 10.1103/PhysRevLett.115.033201

PACS numbers: 34.35.+a, 34.20.Cf

Dispersion forces between neutral atoms are often inter-
preted as a result of the quantum fluctuation of both the
electromagnetic (EM) field and the atomic charges [1,2]. A
prominent example of those are van der Waals (vdW) forces
acting between neutral atoms and molecules, which are
important in atomic and molecular interferometry where they
influence the measured interference pattern [3]. In quantum
information, vdW forces between Rydberg atoms produce a
Rydberg blockade, which may be exploited to realize
quantum gates [4]. In biophysical and chemical processes
vdW forces are known to play a crucial role for the stability
and assembling of molecules [5].

At zero temperature, two atoms in their ground states undergo a series of virtual transitions to upper levels. It is the coupling of the charges of each atom to the quantum EM field that induces the correlation between their transient dipole moments, giving rise to a nonvanishing vdW interaction. The lifetime of the virtual atomic transitions is very short in comparison to ordinary observation times and thus, the use of stationary quantum perturbation theory is well justified for the calculation of this interaction [6]. For short interatomic distances *r* in comparison to the relevant transition wavelengths (nonretarded regime) the interaction scales as *r*<sup>-6</sup>, while for large distances (retarded regime) it goes like *r*<sup>-7</sup> [1,6–9]. The situation is different for excited atoms. First, excited states are unstable and present finite lifetimes. This implies that, generically, the interaction between excited atoms must depend on time. Second, if any of the transitions from the excited to lower atomic levels is relevant to the interaction, the exchange of resonant photons between the atoms must be considered. The energy of the interaction mediated by resonant photons is usually referred to as resonant van der Waals potential in the literature [10–13]. In the retarded regime the resonant potential overtakes by far the nonresonant one. It is in this regime that different approaches yield conflicting results concerning the spatial

oscillations of the interaction [14–21]. This long-standing
problem is the main motivation of the present Letter.

In the following, we address the time-dependent quan-
tum computation of the interaction between two dissimilar
atomic dipoles, one of which is excited. The excited atom is
taken of type *A* while the atom in its ground state is
considered of a different type *B*. Without loss of generality
we approximate the atoms by two-level systems of resonant
frequencies ω<sub>A</sub> and ω<sub>B</sub>, respectively, with respective line-
widths Γ<sub>A</sub> and Γ<sub>B</sub>. Further, in order to ensure the pertur-
bative nature of the calculation and to avoid resonant
energy transfer we set the detuning Δ<sub>AB</sub> ≡ ω<sub>A</sub> − ω<sub>B</sub> such
that |Δ<sub>AB</sub>| > (Γ<sub>A</sub> + Γ<sub>B</sub>)/2 and |Δ<sub>AB</sub>| ≫ 〈W(T)⟩/ħ, with
W(T) being the interaction Hamiltonian at the time of
observation *T*. Since the observation is made for atom *A*
excited, we must have *T* ≲ 2πΓ<sub>A</sub><sup>-1</sup>. Lastly, we assume
without much loss of generality Γ<sub>A,B</sub> < |Δ<sub>AB</sub>| ≪ ω<sub>A,B</sub>,
which is easily met by pairs of alkali atoms. We will refer
to this condition as *quasiresonant*. We will see that it allows
for a great reduction in the number of calculations and
makes the resonant potential dominant at all distances.
We will show that the interaction energy oscillates both
in time and in space. It contains time-independent terms
that oscillate in space with frequency 2ω<sub>A</sub>/c, and time-
dependent terms which oscillate in time with frequency
Δ<sub>AB</sub> and in space with frequency 2ω<sub>B</sub>/c. We compare our
results to previous conflicting approaches and discuss in
detail to which extent they provide compatible results.

We aim at computing the EM energy of atom A due to
the presence of atom B. To this end we apply standard time-
dependent quantum perturbative techniques in the electric
dipole approximation [22]. At any given time T the state
of the two-atom-vacuum system can be written as |Ψ(T)⟩ =
∪(T)|Ψ(0)⟩, where the state of the system at time 0 is
|Ψ(0)⟩ = |A<sub>+</sub>⟩ ⊗ |B<sub>−</sub>⟩ ⊗ |0<sub>γ</sub>⟩. In this expression ⟨A, B⟩<sub>+,−</sub>
label the upper or lower internal states of the atoms A and
$B$, respectively, and $|0_\gamma\rangle$ is the EM vacuum state. $\cup(T)$ denotes the time evolution operator in the Schrödinger representation,

$$ \cup(T) = T \exp\left\{-i\hbar^{-1} \int_{0}^{T} dt \left[H_A + H_B + H_{EM} + W\right]\right\}. $$

In this equation $H_A + H_B$ is the free Hamiltonian of the internal atomic states, $\hbar\omega_A|A_+\rangle\langle A_+| + \hbar\omega_B|B_+\rangle\langle B_+$, while the Hamiltonian of the free EM field is $H_{EM} = \sum_{k,\epsilon} \hbar\omega(a_{k,\epsilon}a_{k,\epsilon}^{\dagger} + 1/2)$, where $\omega = ck$ is the photon frequency, and the operators $a_{k,\epsilon}^{\dagger}$ and $a_{k,\epsilon}$ are the creation and annihilation operators of photons with momentum $\hbar k$ and polarization $\epsilon$, respectively. Finally, the interaction Hamiltonian reads $W = W_A + W_B$, with $W_{A,B} = -d_{A,B} \cdot E(R_{A,B})$. In this expression $d_{A,B}$ are the electric dipole operators of each atom and $E(R_{A,B})$ is the electric field operator evaluated at the position of each atom, which can be written in the usual manner as a sum over normal modes as [6]

$$ \begin{aligned} E(\mathbf{R}_{A,B}) &= \sum_{\mathbf{k}} \mathbf{E}_{\mathbf{k}}^{(-)}(\mathbf{R}_{A,B}) + \mathbf{E}_{\mathbf{k}}^{(+)}(\mathbf{R}_{A,B}) \\ &= i \sum_{\mathbf{k},\epsilon} \sqrt{\frac{\hbar c k}{2\mathcal{V}\epsilon_0}} [\epsilon a_{\mathbf{k}} e^{i\mathbf{k} \cdot \mathbf{R}_{A,B}} - \epsilon^{\dagger} a_{\mathbf{k}}^{\dagger} e^{-i\mathbf{k} \cdot \mathbf{R}_{A,B}}], \end{aligned} $$

where $\mathcal{V}$ is a generic volume and $E_{\mathbf{k}}^{(\tau)}$ denote the annihilation or creation electric field operators of photons of momentum $\hbar\mathbf{k}$, respectively. While the internal atomic and EM degrees of freedom are quantum variables, the position vectors $\mathbf{R}_{A,B}$ are classical variables. We emphasize here that we do not make further simplifications to these potentials. In particular, we do not replace the EM response of any of the atoms by its ordinary polarizability, as is the case in Ref. [20].

Next, considering $W$ as a perturbation to the free Hamiltonians, the unperturbed time-evolution operator for atom and free photon states is $\cup_0(t) = \exp[-i\hbar^{-1}(H_A + H_B + H_{EM})t]$. In order to make contact with a realistic setup, we imagine that atom A starts being excited at time $-\tau$ by a laser pulse of duration $\tau < |\Delta_{AB}^+|$. This fixes our temporal resolution and implies that at time $\approx 0$ the initial state $|\Psi(0)\rangle$ is well defined within a time interval of the order of $\tau$. We are now ready to compute the EM energy of atom A due to the presence of atom B at any time $T$ such that $0 \lesssim T \lesssim 2\pi\Gamma_A^{-1}$,

<figure><img src="image_17.png" /><figcaption>FIG. 1. Diagrammatic representation of the twelve time-ordered processes that contribute to $\langle W_A(T) \rangle$ at the lowest order in $W$. The time variable runs along the vertical.</figcaption></figure>

$$ \langle W_A(T) \rangle = -\langle \Psi(0) | \cup^{\dagger}(T) \mathbf{d}_A \cdot \mathbf{E}(\mathbf{R}_A) \cup(T) | \Psi(0) \rangle. \quad (1) $$

The above expression admits an expansion in powers of $W$, which can be developed out of the time-ordered exponential equation for $\cup(T)$, $\cup(T) = \cup_0(T)T \exp \int_0^T d\tau \cup_0^+(\tau) W \cup_0(\tau) dt$. At leading order, Eq. (1) contains a series of terms of fourth order in $W$ where an electric field operator creates or annihilates a photon at time $T$ at the position of atom A. They correspond to the twelve well-known time-ordered diagrams of Fig. 1 [6,15]. In the time-dependent approach, each diagram contributes to $\langle W_A(T) \rangle$ with two terms in which $W_A$ is flanked by two $\cup$ matrices which make the system evolve, in opposite time directions, from the initial state to two different states at time $T$, which differ from one another in the state of atom A and in the number of photons by one unit. In quasiresonant conditions, the greatest contribution to $\langle W_A(T) \rangle$ comes from diagram (a), in which two doubly resonant photons are exchanged one after the other. Doubly resonant photons are those emitted by one of the atoms in its upper level and absorbed by the other atom in its lower level, while for nonresonant photons the emission and absorption processes are inverted. Last, simply resonant photons are those emitted and absorbed by both atoms in either their upper or lower levels. In addition, the diagrams (b)–(f) of Fig. 1, which contain both doubly resonant and nonresonant photons, provide terms which make it possible to extend the frequency integrals of diagram (a) into the negative domain. Their contribution is indeed essential for establishing causality (cf. Ref. [23]). All the other contributions from these and from the rest of the diagrams are at the most of the order of $\Delta_{AB}/\omega_{A,B}$ times smaller and hence negligible. Putting everything together, transforming the sums over photon momenta into continuum integrals and imposing the causality condition $T > 2R/c$ with $\mathbf{R} = \mathbf{R}_B - \mathbf{R}_A$, we find at leading order,

$$ \begin{aligned} \langle W_A(T) \rangle &\approx \frac{1}{2\hbar^3} \int_{-\infty}^{\infty} \frac{\mathcal{V}k^2 dk}{(2\pi)^3} \int_{-\infty}^{\infty} \frac{\mathcal{V}k'^2 dk'}{(2\pi)^3} \int_0^{4\pi} d\Omega \int_0^{4\pi} d\Omega' \\ &\qquad \times \int_0^{\tau'} dt'' \langle \Psi(0) | \mathbf{d}_A \cdot \mathbf{E}_{\mathbf{k}'}^{(-)}(\mathbf{R}_A) \cup_0(T-t) \mathbf{d}_B \cdot \mathbf{E}_{\mathbf{k}'}^{(+)}(\mathbf{R}_B) \cup_0(t-t') \mathbf{d}_B \cdot \mathbf{E}_{\mathbf{k}}^{(-)}(\mathbf{R}_B) \\ &\qquad \times \cup_0(t'-t'') \mathbf{d}_A \cdot \mathbf{E}_{\mathbf{k}}^{(+)}(\mathbf{R}_A) \cup_0(t'') | \Psi(0) \rangle + [k \leftrightarrow k']^{\dagger}. \end{aligned} \quad (2) $$
The time integrals of the time-evolution operators in Eq. (2) yield a series of terms with poles along the real axis,

$$ \frac{c}{\Delta_{AB}(k-k_A)(k'-k_A)} - \frac{c \cos(\Delta_{AB}T)}{\Delta_{AB}(k-k_B)(k'-k_B)} - \frac{\cos[(\omega-\omega_A)T]}{(k-k_A)(k-k_B)(k'-k)} + \frac{\cos[(\omega'-\omega_A)T]}{(k'-k_A)(k'-k_B)(k'-k)}. \quad (3) $$

Further, the development of this expression contains terms in which both photons resonate either with the

transition of atom A or with the transition of atom B only. This is a direct consequence of energy conservation. Upon integration in frequencies, the former terms are time independent while the latter oscillate in time as ~cos $\Delta_{AB}T$. Important is the fact that only the first term in Eq. (3) arises in the stationary approach [6,9,15]. However, the integration in frequencies of the third and fourth terms provides additional time-independent contributions which are missing in the stationary approach. Last, replacing the time integrals in Eq. (2) with the result (3) and integrating in orientations and frequencies, we obtain [24]

$$ \begin{aligned} \langle W_A(T) \rangle = & \frac{\mathcal{U}_{ijpq}}{R^6} [\beta^{ij}\beta^{pq} - k_A^2 R^2 (\beta^{ij}\beta^{pq} + 2\alpha^{ij}\beta^{pq}) + k_A^4 R^4 \alpha^{ij}\alpha^{pq}] \cos(2k_A R) \\ & - \frac{2\mathcal{U}_{ijpq}}{R^6} k_A [\beta^{ij}\beta^{pq} - k_B^2 R^2 (\beta^{ij}\beta^{pq} + 2\alpha^{ij}\beta^{pq}) + k_B^4 R^4 \alpha^{ij}\alpha^{pq}] \cos(2k_B R + \Delta_{AB}T) \\ & - \frac{2\mathcal{U}_{ijpq}}{R^6} k_B [\beta^{ij}\beta^{pq} - k_B^2 R^2 \alpha^{ij}\beta^{pq}] \sin(2k_B R + \Delta_{AB}T) + \frac{\mathcal{U}_{ijpq}}{R^6} [1 + \dots + (k_{A,B}R)^4] \mathcal{O}(\Delta_{AB}/\omega_{A,B}) + \dots, \end{aligned} \quad (4) $$

where $\mathcal{U}_{ijpq} = \mu_A^4 \mu_B^4 \mu_B^2 / [(4\pi\epsilon_0)^2 \hbar \Delta_{AB}]$, $\mu^A = \langle A^- | d_A | A^+ \rangle$, $\mu^B = \langle B^- | d_B | B^+ \rangle$ and $\beta^{ij} = \delta^{ij} - 3R^i R^j/R^2$, $\alpha^{ij} = \delta^{ij} - R^i R^j/R^2$. It is worth stressing that the Heaviside function in Eq. (2), together with the time-order prescription, do not only guarantee causality, but also determine univocally the contours of integration over frequencies in the complex plane when taking the principal value [24]. The last term in Eq. (4) indicates the order of the leading corrections to the dominant doubly resonant photon exchange terms of Eq. (2) [25]. As anticipated, the time-independent terms of Eq. (4) oscillate only in space with frequency $2k_A$. On the contrary, the time-dependent terms oscillate in time with frequency $\Delta_{AB}$ and in space with frequency $2k_B$. Only for large integration times, $\delta T \gg |\Delta_{AB}^{-1}|$, their time average vanishes. In the short time limit, $T \to 2R/c$, $\langle W_A(T) \rangle$ vanishes identically at our order of approximation. This is a consequence of the fact that, in order to establish the interaction, it is necessary that the excitation be transferred actually to atom B. For $T > R/c$, the probability of excitation of atom B oscillates in time as $|\langle \Psi(T) | A_- \rangle \otimes |B_+ \rangle \otimes |0_+ \rangle|^2 \sim \sin^2[\Delta_{AB}(R/c - T)]/2$, being maximum for the first time at $T = R/c + \pi|\Delta_{AB}^{-1}|$. Correspondingly, $\langle W_A(T) \rangle$ becomes maximum for the first time at $T = 2R/c + \pi|\Delta_{AB}^{-1}|$. The lapse $R/c$ between these two times is the time for a photon to travel back from $R_B$ to $R_A$ after the excitation of atom B.

A long-standing debate exists in the literature concerning the spatial oscillations of the two-atom interaction in the retarded regime when one of the atoms is excited [14–21]. The existence of spatial oscillations is indeed supported by experiments [26,27]. According to our findings, for $k_{A,B}R \gg 1$ and $T > 2R/c$, the interaction oscillates both in time and in space as

$$ \begin{aligned} \langle W_A(T) \rangle &= \frac{\mathcal{U}_{ijpq}}{R^2} \alpha^{ij} \alpha^{pq} [k_A^4 \cos(2k_A R) - k_B^4 \cos(2k_B R + \Delta_{AB}T)] \\ &= \frac{-2\mathcal{U}_{ijpq}}{R^2} \alpha^{ij} \alpha^{pq} k_A^4 \sin[\Delta_{AB}(R/c - T/2)] \\ &\quad \times \sin[k_A(R + cT/2) + k_B(R - cT/2)]. \end{aligned} \quad (5) $$

From the last expression we read that, at fixed time, the interaction is modulated by long-range oscillations of frequency $\Delta_{AB}/c$, while short-range oscillations take place at frequency $k_A + k_B$. Also, as a function of time, the interaction is modulated by oscillations of frequency $\Delta_{AB}$. In Fig. 2 we plot the energy of the interaction between two alkali atoms, one of $^{87}$Rb which is excited to the state $5P_{1/2}$ and another one of $^{40}$K which is in its ground state, in the retarded regime.

In contrast to our result, the stationary approach of Power and Thirunamachandran in Ref. [14] predicts no oscillations for $\langle W_A \rangle$ in the far field. The key point in their calculation is the addition of small imaginary parts to the resonant frequency of atom A in such a way that poles get shifted off the real axis. They used the prescription that a positive (negative) imaginary part must be added for emitted (absorbed) photons in order to account for the finite line-width of the excited atom. In particular, for $\Delta_{AB} \ll \omega_{A,B}$ the dominant term in their stationary calculation is the first one in Eq. (3), but with the real poles shifted as $[\Delta_{AB}(k - k_A - i\eta/c)(k' - k_A + i\eta/c)]^{-1}$, $\eta \to 0^+$. After integrating in orientations an analogous equation to Eq. (2) [6,24], they must have obtained for the energy in the far field limit, $k_A R \gg 1$,
<figure><img src="image_18.png" /><figcaption>FIG. 2 (color online). Graphical representation of the interaction between a $^{87}$Rb atom in state 5P$_{1/2}$ ($k_A = 2\pi12578.95$ cm$^{-1}$) and a $^{40}$K atom in its ground state, for $k_A R \gg 1$. The black line corresponds to a snapshot of the interaction at time $3.0 \times 10^{-12}$ s, where the contributions of the D1 and D2 transition lines of the $^{40}$K atom ($k_B^1 = 2\pi12985.17$ cm$^{-1}$ and $k_B^2 = 2\pi13045.876$ cm$^{-1}$, respectively) add up approximately in phase. A long-range period $c\pi/|\Delta_{AB}|$, with $\Delta_{AB}/c = k_A - (k_B^1 + k_B^2)/2$, is identified. The red line corresponds to the time-independent result of the causal-adiabatic approximation. The average over dipole orientations has been taken.</figcaption></figure>

$$
- \frac{\mathcal{U}_{ijpq}}{4\pi^2 R^2} \alpha^{ij} \alpha^{pq} \int_{-\infty}^{+\infty} dk \int_{-\infty}^{+\infty} dk' k^2 k'^2 \\
\times \frac{e^{i(k+k')R} + e^{-i(k+k')R} - e^{i(k-k')R} - e^{-i(k-k')R}}{(k-k_A - i\eta/c)(k'-k_A + i\eta/c)}, \quad (6)
$$

with $\eta \to 0^+$. Since the pole in $k$ lies on the upper half of the complex plane and the pole in $k'$ lies on the lower half, the only nonvanishing contribution to the above integral comes from the term proportional to $e^{i(k-k')R}$. As the real part of the poles is in both cases $k_A$, taking the limit $\eta \to 0^+$, the exponent vanishes after evaluating the residues and the integral yields the nonoscillating result $(\mathcal{U}_{ijpq}/R^2)k_A^4\alpha^{ij}\alpha^{pq}$.

In the previous stationary calculation of McLone and Power [15] and Gomberoff et al. [16] the poles in Eq. (6) were not shifted. As a result, when taking the principal value of the integrals in Eq. (6) with $\eta=0$, the four exponentials in the numerator contribute as $\sim 2\cos 2k_A R + 2 = 4\cos^2 k_A R$ after adding up the residues, yielding the oscillating result $(\mathcal{U}_{ijpq}/R^2)k_A^4\alpha^{ij}\alpha^{pq}\cos^2 k_A R$. As mentioned after Eq. (3), some time-independent terms are missing in the stationary calculation, which explains the discrepancy of this result with the time-independent component of ours in Eq. (5). In fact, the correct time-independent result can be obtained by switching on adiabatically the interaction potential $W$, within the time-dependent causal approach. That is, by replacing $W_{A,B}(t)$ in Eq. (2) with $W_{A,B}(t)e^{i\eta}$, $\eta \to 0^+$, and extending the lower limit of the time integrals to $-\infty$ [24]. This results in an equation analogous to Eq. (6) but with both poles shifted along the

imaginary axis in the same direction an amount $i\eta/c$, which yields $(\mathcal{U}_{ijpq}/R^2)k_A^4\alpha^{ij}\alpha^{pq}\cos 2k_A R$ for $k_A R \gg 1$.

Recently, Safari and Karimpour have published a Letter [19] where they claim to obtain for $k_A R \gg 1$ the same oscillating behavior as Gomberoff, McLone, and Power [16]. However, a straightforward comparison of Eq. (19) of Ref. [19] and Eqs. (14), (26) of Ref. [16] reveals that this is indeed not the case. Whereas the result of the latter is the one outlined above, $\sim \cos^2 k_A R$, the authors of the former have found $\sim \cos 2k_A R$, despite the fact that both approaches are based on fourth order stationary perturbation theory. The origin of the discrepancy is in the algebraic manipulation inherited by the authors of Ref. [19] from Ref. [9]. In Appendix B of Ref. [9] the authors have tried to express the total contribution of the twelve diagrams of Fig. 1 as a single frequency integral whose integrand is a function of the ordinary polarizabilities of the two atoms. In doing so by means of Eq. (B2) of Ref. [9], the authors have replaced effectively the denominator of Eq. (6), which is a symmetric and separable function of $k$ and $k'$ for $\eta = 0$, by the expression $[\Delta_{AB}(k-k_A)]^{-1}[1/(k'-k) + 1/(k'+k)]$, which is neither symmetric nor separable. As a consequence, that replacement makes the frequency integrals depend arbitrarily on the order of integration. Next, integrating in $k'$ first and in $k$ later, one obtains $(\mathcal{U}_{ijpq}/R^2)k_A^4\alpha^{ij}\alpha^{pq}\cos 2k_A R$, which agrees with Eq. (19) of Ref. [19] for $k_A R \gg 1$, $\Delta_{AB} \ll \omega_{A,B}$, upon averaging in atomic orientations. Interestingly, this result equals the time-independent term of Eq. (5). However, this coincidence can only be accidental, since the above replacement and the subsequent prescription on the order of integration are neither connected to the time-dependent terms of Eq. (3), which cause the actual discrepancy with respect to the result of Refs. [15,16], nor to the causal-adiabatic approximation.

It is worth noting that while we have invoked the existence of finite lifetimes $\sim \Gamma_{A,B}^{-1}$ in order to impose physical constraints on the detuning $\Delta_{AB}$ and on the observation time $T$, no explicit reference to these quantities appear in our expression for $\langle W_A(T) \rangle$. As a matter of fact, only the emission through the exchange of resonant photons between the two atoms has been implicitly accounted for in our calculation of $\langle W_A(T) \rangle$. However, our calculation lacks the inclusion of the spontaneous emission of each atom into free space, whose rates are $k_{A,B}^3\mu_{A,B}^2/3\pi\epsilon_0\hbar$, respectively. The processes corresponding to the latter phenomenon are generally unimportant in comparison to those depicted in Fig. 1 since their leading contribution to $\langle W_A(T) \rangle$ is of order $\mathcal{O}(W^6) \sim \mu_A^3\mu_B^2$, $\mu_A^2\mu_B^4$—see Fig. 3. They might only be relevant for the case that the lifetimes are of the order of the temporal frequency of the interaction, $\Gamma_{A,B} \sim |\Delta_{AB}|$, but they cannot affect in any case the oscillatory behavior found here for the terms of order $\mu_A^2\mu_B^2$. This argument opposes the reasons
<figure><img src="image_19.png" /><figcaption>FIG. 3. Dominant diagrams at order $O(W^6)$, which incorporate in the calculation of $\langle W_A(T) \rangle$ the effect of photon emission into free space. In (i) and (ii) the photon $\gamma$ is emitted into free space from atom A, whereas in (iii) it is emitted from atom B.</figcaption></figure>

given in Ref. [14] to add imaginary shifts to the real poles at $O(W^4)$. We finalize by mentioning that Berman has shown in a recent publication [21] how to introduce spontaneous emission in a consistent manner in an adiabatic approximation. Assuming that atom A is excited adiabatically with $\tau$, $\delta T \gg |\Delta_{AB}^{-1}|$, he has obtained the correct time-independent result.

In this Letter we have shown that the van der Waals interaction between two dissimilar atoms, one of which is initially prepared in an excited state, presents generically oscillations both in time and in space. In quasiresonant conditions the interaction is dominated at all distances by the exchange of doubly resonant photons between the two atoms. It is modulated in space by long-range oscillations of frequency $\Delta_{AB}/c$, while short-range oscillations take place at frequency $k_A + k_B$. The time frequency is $\Delta_{AB}$, which determines the rate at which the excitation is transferred to atom B. In the retarded regime the interaction takes the form of Eq. (5). Only for large integration times, $\delta T \gg |\Delta_{AB}^{-1}|$, that expression reduces to a time-independent term which oscillates in space with frequency $2k_A$. The latter is also equivalent to the result of the causal-adiabatic approximations [21,24]. It does not agree, however, with the result of stationary perturbation theory [15,16].

We thank Paul Berman and Marie-Pascale Gorza for fruitful discussions. Financial support from ANR-10-IDEX-0001-02-PSL and ANR-13-BS04-0003-02 is gratefully acknowledged.

<footnote>* donaire@lkb.upmc.fr, mad37ster@gmail.com</footnote>

<footnote>[1] P. W. Milonni, *The Quantum Vacuum* (Academic Press, San Diego, 1994).</footnote>

<footnote>[2] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, *Photons and Atoms. Introduction to Quantum Electrodynamics* (Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, 2004).</footnote>

<footnote>[3] T. Juffmann, A. Milic, M. Müllneritsch, P. Asenbaum, A. Tsukernik, J. Tüxen, M. Mayor, O. Cheshnovsky, and M. Arndt, Nat. Nanotechnol. 7, 297 (2012).</footnote>

<footnote>[4] M. Saffman, T. G. Walker, and K. Mølmer, Rev. Mod. Phys. 82, 2313 (2010).</footnote>

<footnote>[5] J. Israelachvili, *Intermolecular and Surface Forces* (Academic Press, New York, 1992).</footnote>

<footnote>[6] D. P. Craig and T. Thirunamachandran, *Molecular Quantum Electrodynamics* (Dover, New York, 1998).</footnote>

<footnote>[7] F. London, Z. Phys. 63, 245 (1930).</footnote>

<footnote>[8] H. B. G. Casimir and D. Polder, Phys. Rev. 73, 360 (1948).</footnote>

<footnote>[9] H. Safari, S. Y. Buhmann, D.-G. Welsch, and H. T. Dung, Phys. Rev. A 74, 042101 (2006).</footnote>

<footnote>[10] J. M. Wylie and J. E. Sipe, Phys. Rev. A 30, 1185 (1984); 32, 2030 (1985).</footnote>

<footnote>[11] S. Scheel and S. Y. Buhmann, Acta Phys. Slovaca 58, 675 (2008).</footnote>

<footnote>[12] M.-P. Gorza and M. Ducloy, Eur. Phys. J. D 40, 343 (2006).</footnote>

<footnote>[13] M. Donaire, Phys. Rev. A 83, 022502 (2011); 85, 052518 (2012).</footnote>

<footnote>[14] E. A. Power and T. Thirunamachandran, Phys. Rev. A 51, 3660 (1995).</footnote>

<footnote>[15] R. R. McLone and E. A. Power, Proc. R. Soc. A 286, 573 (1965).</footnote>

<footnote>[16] L. Gomberoff, R. R. McLone, and E. A. Power, J. Chem. Phys. 44, 4148 (1966).</footnote>

<footnote>[17] Y. Sherkunov, Phys. Rev. A 72, 052703 (2005).</footnote>

<footnote>[18] H. R. Haakh, J. Schiefele, and C. Henkel, Int. J. Geom. Methods Mod. Phys. Conf. Series 14, 347 (2012).</footnote>

<footnote>[19] H. Safari and M. R. Karimpour, Phys. Rev. Lett. 114, 013201 (2015).</footnote>

<footnote>[20] L. Rizzuto, R. Passante, and F. Persico, Phys. Rev. A 70, 012107 (2004).</footnote>

<footnote>[21] P. R. Berman, Phys. Rev. A 91, 042127 (2015).</footnote>

<footnote>[22] J. J. Sakurai, *Advanced Quantum Mechanics* (Addison-Wesley, Reading, MA, 1994).</footnote>

<footnote>[23] P. R. Berman and B. Dubetsky, Phys. Rev. A 55, 4060 (1997).</footnote>

<footnote>[24] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.115.033201 for details about the evaluation of Eq. (2).</footnote>

<footnote>[25] Some of the diagrams of Fig. 1 discarded provide also terms with poles in both $k$ and $k'$. However, after integrated in frequencies, they yield contributions of the order of $\Delta_{AB}/(\omega_A + \omega_B)$ times smaller than the ones found here, and hence they are negligible. For instance, in diagram (g) this is due to its lack of resonance with the transition of atom B, since both photons are simply resonant.</footnote>

<footnote>[26] M. A. Wilson, P. Bushev, J. Eschner, F. Schmidt-Kaler, C. Becher, R. Blatt, and U. Dorner, Phys. Rev. Lett. 91, 213602 (2003).</footnote>

<footnote>[27] P. Bushev, A. Wilson, J. Eschner, C. Raab, F. Schmidt-Kaler, C. Becher, and R. Blatt, Phys. Rev. Lett. 92, 223602 (2004).</footnote>
### 2.3.2 Multiscattering formalism

In ref. [29], included hereafter, we studied the statistical properties of the Casimir-Polder interaction between a nanosphere and a dilute inhomogeneous medium. The inhomogeneous medium is itself modeled by a collection of nanospheres of size $a$ randomly distributed in a semi-infinite half space at a distance $z$ from the probe nanosphere. The nanospheres forming the inhomogeneous medium are distributed with a constant density $n$ and we have considered the dilute limit $na^3 \ll 1$.

The interaction energy between the probe sphere and the inhomogeneous medium is calculated using the multiscattering formalism presented before. If the inhomogeneous medium consists of $N$ nanospheres randomly distributed in a volume $V$ at a fixed density $n = N/V$, the interaction energy $U(z)$ between this medium and an additional probe sphere placed at a distance $z$ from the surface of the inhomogeneous medium is

$$U(z) = U^{(N+1)}(z) - U^{(N+1)}(z \to \infty) = U^{(N+1)}(z) - U^{(N)}(z) \quad (88)$$

where $U^{(N+1)}(z)$ is the total Casimir interaction energy between the $N$ spheres forming the inhomogeneous medium and the probe sphere.

The interaction energy between a collection of $N$ (or $N+1$) spheres can be readily calculated using the expression given in Eq. (36) as we chose to perform a calculation at zero temperature. We represent the scattering and translation operators in a basis of spherical waves. In this basis, the scattering operators of the spheres entering the matrix $S$ are simply Mie reflection coefficients. Closed expression for the translation operators entering the $W$ matrix can be found in ref. [30].

We have performed numerical calculations of the Casimir interaction energy $U(z)$ for different realizations of the inhomogeneous medium. In addition, using some approximations valid in the dilute limit $na^3 \ll 1$, we were able to derive an analytical expression for the probability distribution function obeyed by the Casimir interaction $U(z)$. Comparison between this analytical expression and histograms of the numerical calculations is found to be very good.
Statistical approach to Casimir-Polder potentials in heterogeneous media

Nicolas Cherroret, Romain Guérout, Astrid Lambrecht, and Serge Reynaud

Laboratoire Kastler Brossel, UPMC-Sorbonne Universités, CNRS, ENS-PSL Research University, Collège de France,
4 Place Jussieu, 75005 Paris, France

(Received 31 July 2015; published 26 October 2015)

We explore the statistical properties of the Casimir-Polder potential between a dielectric sphere and a three-dimensional heterogeneous medium, by means of extensive numerical simulations based on the scattering theory of Casimir forces. The simulations allow us to confirm recent predictions for the mean and standard deviation of the Casimir potential, and give us access to its full distribution function in the limit of a dilute distribution of heterogeneities. These predictions are compared with a simple statistical model based on a pairwise summation of the individual contributions of the constituting elements of the medium.

DOI: 10.1103/PhysRevA.92.042513

PACS number(s): 31.30.jh, 34.35.+a, 42.50.Ct, 42.25.Fx

I. INTRODUCTION

Materials placed in close vicinity to each other modify
the modes of the electromagnetic field. This results in a
change of the vacuum energy, which eventually manifests
itself as a net force known as the Casimir force [1,2].
The Casimir force has been the subject of a number of
experimental investigations at object separations ranging from
tens of nanometers to a few micrometers. Starting with the
experiments by Lamoreaux [3] and Mohideen [4], the Casimir
effect has experienced an enormous increase in experimental
activities in recent years [5–20].

Theoretical approaches to the Casimir force are usually
built on an effective medium description of the interacting
materials. Within such an approximation, the local details of
the materials' microstructure are neglected and the objects
are described by macroscopic, spatially independent dielectric
constants. While the effective medium description is in general
quite satisfactory for describing dense materials that indeed
look homogenous at the typical scales of the Casimir force,
this is not necessarily the case for strongly heterogeneous
("disordered") media that are made of many constituting
elements ("scatterers") well separated from one another.
Examples of such heterogeneous systems include nanoporous
materials [21], clouds of cold atoms [22], and, in a slightly
different context, corrugated surfaces [23,24].

From a theoretical viewpoint, interaction phenomena in-
volving strongly heterogeneous materials have been little
studied. Seminal works on that subject considered the thermal
Casimir interaction between slabs made of piled layers
separated from random distances (one-dimensional disor-
der) [25,26]. The question of disorder was also addressed
recently [27] in the context of the Casimir-Polder (CP)
interaction [28] between a sphere and a plate [29]. In a recent
work finally, the CP interaction between a dielectric sphere
(or an atom) and a three-dimensional disordered dielectric
material was also investigated [30]. This is the scenario we
consider in the present paper.

When a probe sphere or an atom interacts with a spatially
heterogeneous material such as a semi-infinite disordered
medium, the CP potential naturally fluctuates in space. In
other words, the Casimir interaction depends on the specific
statistical realization of the disorder. A shared conclusion of
Refs. [25–27,30] is that when the two objects are far enough

from each other, precisely when the distance between them
is large compared to the typical separation between two
heterogeneities, the value of the Casimir potential from a
realization to another is well captured by its configuration
average, which coincides with the prediction of the effective
medium prescription. In strong contrast, at smaller distances
fluctuations of the potential become larger than its mean,
which is consequently no longer representative. In practice,
this conclusion is crucial for measurements of quantum reflec-
tion [21,31–34], and more generally for any measurement of
the Casimir force involving heterogeneous materials.

In our previous work [30], we developed an exact math-
ematical treatment of the fluctuations of the CP interaction
between a dielectric sphere and a dilute disordered dielectric
medium, and applied it to the calculation of the mean value
of the CP potential and of its standard deviation. In this
paper, we consider the same geometry (recalled in Sec. II),
for which we perform extensive numerical simulations of the
CP potential. The results of these simulations confirm the
predictions of [30] (Sec. III), and additionally allow us to
compute the full probability distribution of the CP potential
which, for a given distribution of the scatterers, does not
depend on the microscopic properties of the latter. In a second
time (Sec. IV), we present a simple statistical model based on
a pairwise summation of the individual contributions of the
scatterers, and confront it with the simulations. Concluding
remarks are collected in Sec. V.

II. MEAN AND STANDARD DEVIATION OF THE
CASIMIR-POLDER POTENTIAL

We address the CP interaction between a probe dielectric
sphere (placed in a vacuum) of static polarizability α₀ (here
and in the following, polarizabilities are expressed in SI units
divided by ε₀) and a semi-infinite, three-dimensional disor-
dered medium consisting of a collection of many scatterers,
as illustrated in Fig. 1. We denote by z the distance between
the sphere and the surface of the disordered medium. For
definiteness, in this paper we restrict our discussion to the
retarded regime of the Casimir interaction where z much
exceeds the resonance wavelength λ₀ of the probe sphere (the
treatment of the opposite limit z ≪ λ₀ is analogous). Scatterers
are also modeled by dielectric spheres of size a and of static
polarizability αs. Throughout the paper, we assume that they
<figure><img src="image_20.png" /><figcaption>FIG. 1. (Color online) We consider the Casimir-Polder interaction between a dielectric sphere (placed in vacuum) and a semi-infinite disordered medium. The disordered medium consists of a collection of dielectric spheres (size a, density n) whose positions are uniformly distributed in space.</figcaption></figure>

are uniformly distributed in space with density n, and we consider the case of a dilute disordered medium, such that the average distance between the scattering spheres is larger than their size, na³ ≪ 1. This is the typical regime where the statistical fluctuations of the CP potential are the largest [30]. In the opposite limit na³ ~ 1 of scatterers very close to each other, the approach developed below does not apply but we expect the statistics of the Casimir-Polder potential to be approximately captured by a model where the atom interacts with a rough surface [35].

In [30], the question of fluctuations in the limit na³ ≪ 1 was tackled with the help of a statistical description of the disordered material, in which the CP potential U(z) becomes a random variable. Its mean, Ū(z), and its variance, δU²(z), were calculated from an exact treatment of radiation-matter interaction, based on the scattering approach to Casimir forces [36,37] combined with a diagrammatic description of radiation scattering off the disordered medium [38,39]. In the limit z ≫ λ₀, the following expression for the mean was found [30]:

$$ \bar{U}(z) = \frac{23}{60} n \alpha_s U^*(z), \quad (1) $$

where U*(z) = -3α₀ħc/(32π²z⁴) is the Casimir potential between the probe sphere and a perfect mirror. As announced, the result (1) coincides with the prediction of an effective medium description where the probe sphere interacts with an homogeneous surface of relative permittivity ε = 1 + nαs. The amplitude of fluctuations, quantified by the ratio γ of the standard deviation of the CP potential and its mean, was found to be (for z ≫ λ₀)

$$ \gamma = \frac{\sqrt{\delta U^2(z)}}{|\bar{U}(z)|} \approx \frac{a_1}{\sqrt{n}z^3}, \quad (2) $$

with a₁ ≃ 0.7. Equation (2) indicates that Ū(z), i.e., the prediction of the effective medium theory, is well representative of U(z) only when z ≫ n⁻¹/³. At smaller scales, γ becomes larger than unity and Ū(z) no longer provides a trustful estimation of the interaction.

### III. NUMERICAL SIMULATIONS

#### A. Methodology

We now propose to investigate the statistical properties of the CP potential from exact numerical simulations in the geometry of Fig. 1. For this purpose, we proceed as follows. We generate an ensemble of N dielectric spheres of radius a and frequency-dependent permittivity ε(ω), uniformly distributed in a cube of side L. This system constitutes a disordered medium of average density n = NL⁻³. An additional probe sphere is placed above this cube, at a distance z to the center of one face, as in Fig. 1. Denoting by U(N+1)(z) the total, internal Casimir energy between the N + 1 spheres [40], the CP interaction U(z) is by definition

$$ \begin{aligned} U(z) &= U^{(N+1)}(z) - U^{(N+1)}(z \to \infty) \\ &= U^{(N+1)}(z) - U^{(N)}. \end{aligned} \quad (3) $$

The strategy thus consists of calculating the interaction energy as the difference between the internal energies of N + 1 and N spheres. Within the scattering formalism [36,37], the total Casimir energy between N spheres is given by [40]

$$ U^{(N)} = \frac{\hbar}{2\pi} \int_{0}^{\infty} d\omega \log \det (MM_{\infty}^{-1}). \quad (4) $$

M is a block-square matrix of dimension N with the following structure:

$$ M = \begin{pmatrix} \mathbf{R}_1^{-1} & \mathbf{T}_{1 \to 2} & \dots & \mathbf{T}_{1 \to N} \\ \mathbf{T}_{2 \to 1} & \mathbf{R}_2^{-1} & \dots & \mathbf{T}_{2 \to N} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{T}_{N \to 1} & \mathbf{T}_{N \to 2} & \dots & \mathbf{R}_N^{-1} \end{pmatrix}. \quad (5) $$

The diagonal blocks of M are the inverse of the spheres' reflection operators Rᵢ. The (i, j) off-diagonal block of M contains the translation operator Tᵢⱼ, which relates an outgoing spherical wave centered on rᵢ to an incoming spherical wave centered on rⱼ [40]. Finally, M⁻¹ is the block-diagonal matrix diag(R₁, ..., Rₙ). For the simulations, we express the scattering and translation operators Rᵢ and Tᵢⱼ in a basis of spherical vector waves |ℓmP⟩, with ℓ > 1, -ℓ ≤ m ≤ ℓ, and P = {E, M}. In this basis, the matrix elements of Rᵢ are given by the standard Mie scattering amplitudes [41]. We compute these amplitudes without any approximation, taking into account the full multipole expansion. Finally, we evaluate the matrix elements of Tᵢⱼ using the formalism of Ref. [42].

#### B. Mean and standard deviation

Making use of the approach described in Sec. III A, we compute the CP potential U(z) between a dielectric sphere of radius a = 10 nm and a disordered medium consisting of N = 32 other identical spheres with the same radius a, uniformly distributed in a cube of side L = 12 µm. For these parameters, the disordered medium is dilute, na³ ≪ 1, and we effectively describing the geometry of a semi-infinite bulk system as long as z ≪ L. For definiteness, we give to all the spheres the permittivity of silicon: ε(ω) = 1 - [ε(0) - 1]ω₀²/(ω² - ω₀² + iγω), where ω₀ = 2πc/λ₀ with λ₀ = 295 nm, ε(0) = 11.6, and γ = 0.03ω₀ [43]. Figure 2 displays the absolute value
<figure><img src="image_21.png" /><figcaption>FIG. 2. (Color online) Absolute value of the mean CP potential between a dielectric sphere and a semi-infinite disordered medium, as a function $z$. Dots are the results of exact numerical simulations and the dashed line is the theoretical prediction (1).</figcaption></figure>

of the disorder-averaged CP potential $\bar{U}(z)$ computed with these parameters, for several values of $z$ (red dots). Each dot is obtained by generating thousands (typically between 3000 and 8000 depending on $z$) of disorder realizations, computing $U(z)$ for each of them and finally averaging the results. In Fig. 2, we also show the theoretical prediction (1) (dashed line), which is in very good agreement with the numerics. The small disagreement visible at large $z$ stems from deviations to the geometry of the semi-infinite medium: when $z$ becomes of the order of $L/2$, the probe sphere starts to be sensitive to the boundaries of the system, and a cross-over toward the sphere-cube geometry is expected. We also show in Fig. 3 the standard deviation of the CP potential relative to its mean, $\gamma$, as a function of $nz^3$. Red dots are the numerical results, and the dashed curve is the theoretical prediction (2). The agreement between theory and numerics is very good, up to small finite-size effects at large $z$.

### C. Probability distribution function

As was pointed out in [30], the fluctuations of $U(z)$ become significant at distances $z \lesssim n^{-1/3}$, when $\gamma$ becomes larger than unity, see Eq. (2) and Fig. 3. This suggests that at small distances, the mean $\bar{U}(z)$ is no longer representative of $U(z)$. In order to confirm this picture, we compute the

<figure><img src="image_22.png" /><figcaption>FIG. 3. (Color online) Relative fluctuations of the CP potential as a function of $nz^3$. Dots are the results of exact numerical simulations and the dashed line is the theoretical prediction (2).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>FIG. 4. (Color online) Scatter plot of the CP potential normalized to its mean value, $s = U/\bar{U}$, computed for 5500 disorder realizations, for $nz^3 = 0.5$.</figcaption></figure>

full probability distribution function $p(s)$ of the CP potential normalized to its mean, $s = U/\bar{U}$, by constructing histograms of the numerical data. In Fig. 4, we show as an example a scatter plot of the data obtained for 5500 disorder realizations, for $nz^3 = 0.5$. The associated histogram $p(s)$ is displayed in Fig. 5, together with the histograms corresponding to three other values of the parameter $nz^3$. A quick look at the distributions in Fig. 5 forthwith confirms the property already outlined by the analysis of $\gamma$: as the sphere gets closer to the disordered medium, the distribution function becomes more and more peaked around a value $s \ll 1$, corresponding to a CP potential much smaller than its mean. In other words, $U(z)$ is no longer a self-averaging quantity. Only when $nz^3 > 1$ does the maximum of the distribution approach $s = 1$. Such a phenomenon was previously observed in the context of the interaction between plates with one-dimensional disorder [25,26]. We see here that it is a quite general property, not restricted to one-dimensional systems.

<figure><img src="image_24.png" /><figcaption>FIG. 5. (Color online) Probability distribution function $p(s = U/\bar{U})$, for four increasing values of $nz^3$. Histograms (vertical lines) are the results of exact numerical simulations, and solid curves are the theoretical prediction, Eq. (15). The lower-right panel also displays as a dashed curve the Gaussian distribution expected in the limit of very large distances, Eq. (17).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>FIG. 6. (Color online) Numerical probability distribution $p(s)$ for two different sets of parameters of the scattering spheres, for $nz^3 = 1$ (histograms, vertical lines). No visible difference is seen between the two histograms. This is confirmed by the theoretical prediction (solid curve), Eq. (15), which depends on the single parameter $nz^3$.</figcaption></figure>

### D. Sensitivity to microscopic parameters

To conclude our numerical study, we address the question of the sensitivity of $p(s)$ with respect to changes in the microscopic properties of the scatterers. For this purpose, we perform additional numerical simulations involving scatterers with a different radius $a = 1$ nm and made of a different material with frequency-independent permittivity $\epsilon = 10$. We also set the radius of the probe sphere to $a = 1$ nm, but keep the same value of the permittivity of silicon as in the previous section. The distribution $p(s)$ obtained for these new parameters is shown in the right panel of Fig. 6, for $nz^3 = 1$. For comparison, we also display the distribution $p(s)$ computed with the parameters of the previous section. No visible difference is seen between the two histograms, which indicates that in the dilute limit, $p(s)$ is in fact a function of the parameter $nz^3$ only. In particular, the parameter $na^3$ is irrelevant. This could have been anticipated since in the limit $na^3 \ll 1$ of independent scatterers, $na^3$ enters both $U$ and $\bar{U}$ within the same prefactor, which thus cancels out when considering the ratio $s = U/\bar{U}$ (see Sec. IV for a general proof). This property is in particular fulfilled by the second moment of the distribution, $\gamma^2 = \delta U^2/\bar{U}^2$, see Eq. (2).

### IV. SIMPLE MODEL

We now develop a simplified statistical description of the CP interaction between a dielectric sphere and a disordered bulk medium, based on a pairwise summation (PWS) approximation [44]. This approximation describes the total CP interaction $U(z)$ as a sum of the pair interaction $\mathcal{E}$ between the probe sphere and each of the $N$ scatterers. It has to be distinguished from the perturbative expansion, as it can in principle be used for nonperturbative pair interactions. In the problem studied in this paper, however, the validity of the two approximations is a consequence of the same assumption of a dilute disorder ($na^3 \ll 1$).

#### A. Pairwise summation

As in the numerical simulations (Sec. III), we consider a situation where the distance $z$ much exceeds the sphere resonance wavelength $\lambda_0$ (in the opposite limit $z \ll \lambda_0$, the reasoning follows exactly the same lines). Consequently, the

<figure><img src="image_26.png" /><figcaption>FIG. 7. (Color online) Left: parametrization of the geometry for the statistical approach of Sec. IV A. Right: typical (rare) disorder realization contributing to the Lifshitz tail: the atom is above a large region free of scatterers, of size $R \gg z$.</figcaption></figure>

interaction potential between the probe sphere and a scatterer located at distance $r$ takes the simple form $\mathcal{E} = -\Gamma_7/r^7$, where $\Gamma_7$ is a constant characteristic of the microscopic properties of the two interacting objects. Let us consider a small spherical cap of volume $dV$ containing $dN$ scatterers, as illustrated in the left panel of Fig. 7. The elementary CP interaction between them and the probe sphere is $dU = \mathcal{E}dN = -\Gamma_7dN/r^7$. As in the simulations, we assume the positions of the scatterers to follow a uniform distribution, so the random variable $dN$ is Poisson distributed. The $m$th cumulant of $dN$, $K_m(dN)$, thus fulfills

$$K_m(dN) = K_1(dN) = ndV, \quad (6)$$

where $n = N/L^3$ is the average density of scatterers. Within the PWS approximation, the cumulant of the elementary CP potential $dU$ involving $dN$ scatterers reads

$$ \begin{aligned} K_m(dU) &= \left(-\frac{\Gamma_7}{r^7}\right)^m K_m(dN) \\ &= \left(-\frac{\Gamma_7}{r^7}\right)^m ndV \equiv dK_m(U), \end{aligned} \quad (7) $$

where we have used the property $K_m(aX) = a^m K_m(X)$ in the first equality, and Eq. (6) in the second. The cumulant of the total CP potential between the probe sphere and the $N$ scatterers is finally obtained by using the parametrization $dV = 2\pi r(r-z)dr$, see Fig. 7, and integrating over $r$ from $r = z$ to $\infty$. This yields

$$ K_m(U) = \frac{(-1)^m 2\pi n \Gamma_7^m}{(7m-3)(7m-2)z^{7m-3}}, \quad (8) $$

from which various statistical properties can be deduced, as we now discuss.

#### B. Mean and standard deviation

As a preliminary application of our statistical approach, we propose to rederive predictions (1) and (2), previously obtained from an exact treatment of the radiation-matter interaction. For simplicity we assume the scatterers to be spheres of static polarizability $\alpha_0$, identical to the probe sphere. The coefficient $\Gamma_7$ then describes the large-distance interaction between two identical spheres. It can be readily evaluated, for instance, from the Casimir-Polder law [1] for the energy between two
atoms [45]:

$$ \Gamma_7 = \frac{23hc\alpha_0^2}{(4\pi)^3}. \qquad (9) $$

The mean $\bar{U}(z)$ is by definition the first-order cumulant: $\bar{U}(z) = K_1(z) = -2\pi n\Gamma_7/(20z^3)$, where we have used Eq. (8) for $m=1$. Combining this expression with Eq. (9), we recover Eq. (1). The variance $\delta U^2(z)$ is, on the other hand, given by the second-order cumulant, $K_2(z) = 2\pi n\Gamma_7^2/(132z^{11})$. Taking the ratio with $K_1(z)$, we find

$$ \gamma = \frac{\sqrt{\delta U^2(z)}}{|\bar{U}(z)|} \approx \sqrt{\frac{50}{33\pi}} \frac{1}{\sqrt{n^{-3/2}}}, \qquad (10) $$

which is nothing but Eq. (2), with an analytic expression of the coefficient $\sqrt{50}/(33\pi) \approx 0.7$.

### C. Probability distribution function

Let us now derive the probability distribution function $p(s = U/\bar{U})$ that was studied numerically in Sec. III C. $p(s)$ is given by the inverse Laplace transform

$$ p(s) = \frac{1}{2\pi i} \int_{\delta-i\infty}^{\delta+i\infty} e^{st} e^{\varphi(t)} dt, \qquad (11) $$

where $\delta$ is greater than the real part of all singularities of $e^{\varphi(t)}$, $\varphi(t)$ is the cumulant generating function of $s$, and can be expressed as a power series of the cumulants (8):

$$ \varphi(t) = \sum_{m=0}^{\infty} \frac{(-t)^m K_m(U)}{m! \bar{U}^m}. \qquad (12) $$

Making use of Eq. (8) and of the relation $\bar{U}(z) = -2\pi n\Gamma_7/(20z^3)$ obtained above, we find

$$ \varphi(t) = \frac{2\pi n z^3}{6} [-1 + e^{-t} - 2\tau^{3/7} \gamma_{4/7}(\tau) + 3\tau^{2/7} \gamma_{5/7}(\tau)], \qquad (13) $$

where we have introduced $\tau = 20t/(2\pi n z^3)$ and where $\gamma_q(\tau) = \int_0^\tau x^{q-1}e^{-x}dx$ is the lower incomplete $\gamma$ function. Since $e^{\varphi(t)}$ has no singularities in the complex plane, we can set $\delta = 0$ in Eq. (11). Furthermore, we have the property $\varphi^*(t) = \varphi(t^*)$, such that after the substitution $t = ix$, Eq. (11) simplifies to

$$ p(s) = \frac{1}{\pi} \int_{0}^{\infty} \mathrm{Re}[e^{ixs} e^{\varphi(ix)}] dx. \qquad (14) $$

Inserting Eq. (13) into this relation, we finally obtain

$$ p(s) = \frac{n z^3}{10} \mathrm{Re} \int_0^\infty \exp \left\{ \frac{2\pi n z^3}{6} \left[ \frac{3}{10} s i \tau - 1 + e^{-i\tau} - 2(i\tau)^{3/7} \gamma_{4/7}(i\tau) + 3(i\tau)^{2/7} \gamma_{5/7}(i\tau) \right] \right\} d\tau. \qquad (15) $$

Distributions $p(s)$ obtained from numerical evaluation of Eq. (15) are displayed in Fig. 5 as solid curves on top of the numerical results of Sec. III. The agreement is excellent for all the values of $n z^3$. Furthermore, we notice that Eq. (15) confirms the conclusion drawn from the simulations in Sec. III D: $p(s)$ depends on the parameter $n z^3$ only, being

completely independent of the microscopic details of the probe sphere and of the scatterers.

### D. Asymptotics

Although the distribution $p(s)$, Eq. (15), has no evident analytic expression, several simple asymptotic limits can be readily examined.

Large-$s$ limit. We describe the limit of large $s$ by expanding the term inside the square brackets in Eq. (15) up to second order in $\tau \ll 1$. Since essentially the values of $\tau$ such that $\tau snz^3 \lesssim 1$ contribute to the integral, this expansion is a good approximation provided $1/(snz^3) \ll 1$. It results in a Gaussian integral which is straightforwardly performed to give

$$ p(s) \approx \frac{\sqrt{33n z^3}}{10} \exp \left[ - \frac{33\pi}{100} n z^3 (s-1)^2 \right]. \qquad (16) $$

Making use of Eq. (10), we rewrite Eq. (16) as

$$ \frac{1}{\sqrt{2\pi \delta U^2}} \exp \left[ - \frac{(U - \bar{U})^2}{2\delta U^2} \right]. \qquad (17) $$

At large $s$, $p(s)$ is thus simply a (normalized) Gaussian distribution of mean $\bar{U}$ and variance $\delta U^2$. The expansion used to derive Eq. (17) being valid as long as $s \gg (nz^3)^{-1}$, the Gaussian shape is a very good approximation of the whole distribution when $nz^3 \gg 1$, which is a direct consequence of the weakness of fluctuations at large distances. We expect this Gaussian distribution to be universal at large distances, as a consequence of the central-limit theorem [many scatterers contribute to $U(z)$ when $nz^3 \gg 1$], regardless the nature of the disordered medium. This conclusion is supported by similar predictions previously made in the context of the thermal Casimir effect in one-dimensional disordered media [25,26], as well as in recent studies of the CP interaction involving quasi-two-dimensional disordered metals [27]. For comparison, we show Eq. (17) in the lower-right panel of Fig. 5 as a dashed curve, for $nz^3 = 40$. Note that in the chosen geometry this limit is difficult to reach in the numerical simulations, because it requires the generation of a significant number of scattering spheres in order to satisfy the condition $nz^3 \gg 1$, while maintaining $z \lesssim L/2$ to avoid finite-size effects. When $nz^3 \ll 1$, Eq. (17) still holds but only in the very far tail of the distribution, $s \gg (nz^3)^{-1} \gg 1$. The physical reason for which we recover a Gaussian tail even for small values of $nz^3$ is the following. Large values of $s$ correspond to particular disorder realizations $\mathcal{D}$ for which the Casimir potential is very large, i.e., for which the density of scatterers $n_D$ underneath the probe atom is very high. Thus, for these specific disorder realizations the condition $n_D z^3 \gg 1$ is effectively fulfilled.

Moderate values of $s$. To describe the limit $s \ll (nz^3)^{-1}$, we expand the term inside the square brackets in Eq. (15) for large $\tau$:

$$ p(s) \approx \frac{n z^3}{10} \mathrm{Re} \int_0^\infty d\tau \exp \left\{ \frac{2\pi n z^3}{6} \left[ \frac{3}{10} i s \tau - 2\Gamma(4/7) \times (i\tau)^{3/7} + 3\Gamma(5/7)(i\tau)^{2/7} \right] \right\}. \qquad (18) $$
<figure><img src="image_27.png" /><figcaption>FIG. 8. (Color online) Distribution $p(s)$ for $nz^3 = 0.005$ (solid red curve), together with the approximate form (19) at moderate $s$ (dashed curve).</figcaption></figure>

We then perform a Wick rotation $\tau = ix$, and express both exponential terms $\exp(\alpha x^{3/7})$ and $\exp(\beta x^{2/7})$ as power series. This gives

$$p(s) = -\frac{nz^3}{10} \operatorname{Im} \int_0^\infty dx \exp\left(-\frac{\pi nz^3}{10}sx\right) \sum_{n,m=0}^\infty \frac{1}{n!m!} \\ \times \left[-\frac{2\pi nz^3}{3} \Gamma(4/7)(-x)^{3/7}\right]^n \left[\pi n z^3 \Gamma(5/7)(-x)^{2/7}\right]^m$$

If we additionally assume $s \gg nz^3$, the terms $(m,n) = (0,1)$ and $(m,n) = (1,0)$ give the leading contribution to $p(s)$ (the term $m=n=0$ is purely real and does not contribute). Keeping only these two terms and computing the remaining integral, we find

$$p(s) \approx \frac{2}{7} 10^{3/7} \frac{(\pi n z^3)^{4/7}}{s^{10/7}} \left[ 1 - \left( \frac{\pi n z^3 s}{10} \right)^{1/7} \right]. \quad (19)$$

Equation (19) holds for $nz^3 \ll s \ll (nz^3)^{-1}$. It is shown in Fig. 8 as a dashed curve for $nz^3 = 0.005$, together with the exact distribution calculated from Eq. (15) (solid red curve).

Small-$s$ limit. We finally consider the low-potential tail $s \ll nz^3$ of $p(s)$. In order to find an asymptotic expansion in that limit, we come back to Eq. (18) and apply the change of variables $\tau = (y/s)^7$. This yields

$$p(s) \approx \frac{7nz^3}{10s^7} \operatorname{Re} \int_0^\infty y^6 dy \exp \left\{ \frac{2\pi n z^3}{6s^6} \left[ \frac{3}{10} iy^7 - 2\Gamma(4/7)i^{3/7}s^3 y^3 + 3\Gamma(5/7)i^{2/7}s^4 y^2 \right] \right\}. \quad (20)$$

In the limit of very small $s$, we are thus led to evaluate

$$I = \operatorname{Re} \int_0^\infty dy y^6 \exp[\Lambda f(y)], \quad (21)$$

where $\Lambda = \pi n z^3/(3s^6) \gg 1$ and $f(y) = (3/10)iy^7 - 2\Gamma(4/7)i^{3/7}s^3 y^3 + 3\Gamma(5/7)i^{2/7}s^4 y^2$. Equation (21) naturally calls for the method of steepest descent. There are five saddle points, solutions of $f'(y) = 0$. Only one of them, denoted by $y_{SP}$, turns out to give a nonzero contribution to $p(s)$:

$$y_{SP} = \left[ \frac{60\Gamma(4/7)s^3}{21} \right]^{1/4} e^{-i\pi/14}. \quad (22)$$

Using Cauchy's theorem, we then deform the path of integration to a path coinciding with the path of steepest descent in the vicinity of $y_{SP}$. This is achieved by expanding $f(y)$ up to second order around $y = y_{SP}$ and performing the change of variables $x = (y - y_{SP}) \exp(-3i\pi/14)$, leaving us with a Gaussian integral whose evaluation leads to

$$p(s) \approx \alpha \frac{\sqrt{nz^3}}{s^{11/8}} \exp\left[-\beta \frac{nz^3}{s^{3/4}}\right], \quad (23)$$

with prefactors $\alpha = (20/7)^{3/8}\Gamma(4/7)^{7/8}$ and $\beta = 16\sqrt{2\pi}(5/7)^{3/4}\Gamma(4/7)^{7/4}/21$. The asymptotic form (23) is completely analogous to the so-called Lifshitz tail that describes the band edge of the density of states of disordered conductors in solid-state physics [46]. Physically, it can be understood from the following qualitative argument. Low values of $s$ are achieved for rare disorder realizations where the probe sphere stands above a large region free of scatterers, as illustrated in the right scheme of Fig. 7. Since the distribution of scatterers is Poissonian (Sec. IV), the distance between the events of this Poisson process follows an exponential distribution. Consequently, the probability to find a large region of size $R$ free of scatterers is $\propto \exp(-cnR^3)$, where $c$ is a numerical constant. In such a configuration, the Casimir potential felt by the atom can be estimated as (see Fig. 7)

$$U = -2\pi n \Gamma_7 \int_{\sqrt{z^2+R^2}}^{\infty} \frac{r(r-z)dr}{r^7} \propto \frac{-n\Gamma_7}{R^4}, \quad (24)$$

to leading order in $z/R \ll 1$. On the other hand, we have seen in Sec. IV B that the average Casimir potential is $\bar{U} \propto -n\Gamma_7/z^4$. Therefore, for the rare disorder realization displayed in Fig. 7, we have $U \propto (z/R)^4\bar{U}$, such that

$$\begin{aligned} p(s) &\sim \exp(-cnR^3) = \exp\left[-cn\left(\frac{z\bar{U}^{1/4}}{U^{1/4}}\right)^3\right] \\ &= \exp\left(-\frac{cnz^3}{s^{3/4}}\right), \end{aligned} \quad (25)$$

which is nothing but the asymptotic form (23).

## V. CONCLUSION

In this paper, we have developed a statistical description of Casimir-Polder potentials from both a numerical and an analytical perspective. This approach is well suited for describing the Casimir interaction between a simple dielectric object and a strongly heterogeneous medium made of a large number of independent constituting elements. It can be readily extended to other geometries and to heterogeneous media characterized by a more complex statistics involving, for instance, a nonuniform or polydisperse distribution of scatterers.

As a first extension of our work, it would be interesting to investigate deviations to the dilute limit $na^3 \ll 1$ where the scatterers can no longer be systematically considered independent. We expect these deviations to primarily affect the far tails of the distribution $p(s)$. A second open question concerns the change in the statistics of the Casimir-Polder potential when the host (homogeneous) medium has a dielectric constant differing from unity. This problem is more difficult
to treat since it now involves a surface, which implies multiple reflections inside and outside the medium. The presence of a host medium may therefore strongly affect the distribution $p(s)$ [in the obvious limit where the dielectric constant of the host medium goes to infinity, one recovers a perfectly reflecting interface and $p(s)$ should tend to the Dirac function $\delta(s - 1)$].

In practice, the distribution $p(s)$ could be experimentally accessed either by moving the sphere over a static disordered medium to record different disorder distributions, or by taking

advantage of a Brownian motion of the scatterers if the measurement process is fast enough. Indeed, in that case different disorder realizations can be obtained by detecting the Casimir force and then letting the scatterers move before carrying out the next measurement. If the measurement process is slow, the effect of the motion of the scatterers is to average the Casimir potential, giving him a value well approximated by Eq. (1) since Doppler shifts have a negligible effect at thermal velocities [47–49].

[1] H. B. G. Casimir, Proc. K. Ned. Akad. Wet. **51**, 793 (1948).

[2] H. B. G. Casimir and D. Polder, Phys. Rev. **73**, 360 (1948).

[3] S. K. Lamoreaux, Phys. Rev. Lett. **78**, 5 (1997).

[4] U. Mohideen and A. Roy, Phys. Rev. Lett. **81**, 4549 (1998).

[5] Th. Ederth, Phys. Rev. A **62**, 062104 (2000).

[6] H. B. Chan, V. A. Aksyuk, R. N. Kleiman, D. J. Bishop, and F. Capasso, Phys. Rev. Lett. **87**, 211801 (2001).

[7] F. Chen, U. Mohideen, G. L. Klimchitskaya, and V. M. Mostepanenko, Phys. Rev. Lett. **88**, 101801 (2002).

[8] R. S. Decca, D. López, E. Fischbach, and D. E. Krause, Phys. Rev. Lett. **91**, 050402 (2003).

[9] R. S. Decca, D. López, E. Fischbach, G. L. Klimchitskaya, D. E. Krause, and V. M. Mostepanenko, Phys. Rev. D **75**, 077101 (2007).

[10] P. J. van Zwol, G. Palasantzas, M. van de Schootbrugge, and J. Th. M. De Hosson, Appl. Phys. Lett. **92**, 054101 (2008).

[11] H. B. Chan, Y. Bao, J. Zou, R. A. Cirelli, F. Klemens, W. M. Mansfield, and C. S. Pai, Phys. Rev. Lett. **101**, 030401 (2008).

[12] G. Jourdan, A. Lambrecht, F. Comin, and J. Chevrier, Europhys. Lett. **85**, 31001 (2009).

[13] J. N. Munday, F. Capasso, and V. A. Parsegian, Nature **457**, 170 (2009).

[14] M. Masuda and M. Sasaki, Phys. Rev. Lett. **102**, 171101 (2009).

[15] S. de Man, K. Heeck, R. J. Wijngaarden, and D. Iannuzzi, Phys. Rev. Lett. **103**, 040402 (2009).

[16] Y. Bao, R. Guérout, J. Lussange, A. Lambrecht, R. A. Cirelli, F. Klemens, W. M. Mansfield, C. S. Pai, and H. B. Chan, Phys. Rev. Lett. **105**, 250402 (2010).

[17] G. Torricelli, I. Pirozhenko, S. Thornton, A. Lambrecht, and C. Binns, Europhys. Lett. **93**, 51001 (2011).

[18] A. O. Sushkov, W. J. Kim, D. A. R. Dalvit, and S. K. Lamoreaux, Nat. Phys. **7**, 230 (2011).

[19] D. Garcia-Sanchez, K. Y. Fong, H. Bhaskaran, S. Lamoreaux, and H. X. Tang, Phys. Rev. Lett. **109**, 027202 (2012).

[20] F. Intravaia, S. Koev, I. W. Jung, A. A. Talin, P. S. Davids, R. S. Decca, V. A. Aksyuk, D. A. R. Dalvit, and D. López, Nat. Commun. **4**, 2515 (2013).

[21] T. A. Pasquini, M. Saba, G.-B. Jo, Y. Shin, W. Ketterle, D. E. Pritchard, T. A. Savas, and N. Mulders, Phys. Rev. Lett. **97**, 093201 (2006).

[22] P. Schneeweiss, M. Gierling, G. Visanescu, D. P. Kern, T. E. Judd, A. Günther, and J. Fortágh, Nat. Nanotech. **7**, 515 (2012).

[23] C. Genet, A. Lambrecht, P. A. Maia Neto, and S. Reynaud, Europhys. Lett. **62**, 484 (2003).

[24] P. A. Maia Neto, A. Lambrecht, and S. Reynaud, Phys. Rev. A **72**, 012115 (2005).

[25] D. S. Dean, R. R. Horgan, A. Naji, and R. Podgornik, Phys. Rev. A **79**, 040101(R) (2009).

[26] D. S. Dean, R. R. Horgan, A. Naji, and R. Podgornik, Phys. Rev. E **81**, 051117 (2010).

[27] A. A. Allocca, J. H. Wilson, and V. Galitski, Phys. Rev. A **91**, 062512 (2015).

[28] F. Intravaia, C. Henkel, and M. Antezza, in *Casimir Physics*, edited by D. Dalvit *et al.*, Lecture Notes in Physics 834 (Springer-Verlag, Berlin, 2011), p. 345.

[29] A. Canaguier-Durand, A. Gérardin, R. Guérout, P. A. Maia Neto, V. V. Nesvizhevsky, A. Yu. Voronin, A. Lambrecht, and S. Reynaud, Phys. Rev. A **83**, 032508 (2011).

[30] N. Cherroret, R. Guérout, A. Lambrecht, and S. Reynaud, Eur. Phys. J. D **69**, 99 (2015).

[31] H. Friedrich, G. Jacoby, and C. G. Meister, Phys. Rev. A **65**, 032902 (2002).

[32] A. Y. Voronin, P. Froelich, and B. Zygelman, Phys. Rev. A **72**, 062903 (2005).

[33] G. Dufour, A. Gérardin, R. Guérout, A. Lambrecht, V. V. Nesvizhevsky, S. Reynaud, and A. Yu. Voronin, Phys. Rev. A **87**, 012901 (2013).

[34] G. Dufour, R. Guérout, A. Lambrecht, V. V. Nesvizhevsky, S. Reynaud, and A. Yu. Voronin, Phys. Rev. A **87**, 022506 (2013).

[35] G. A. Moreno, R. Messina, D. A. R. Dalvit, A. Lambrecht, P. A. Maia Neto, and S. Reynaud, Phys. Rev. Lett. **105**, 210401 (2010).

[36] A. Lambrecht, P. A. Maia Neto, and S. Reynaud, New J. Phys. **8**, 243 (2006).

[37] T. Emig, N. Graham, R. L. Jaffe, and M. Kardar, Phys. Rev. Lett. **99**, 170403 (2007).

[38] A. Lagendijk and B. A. van Tiggelen, Phys. Rep. **270**, 143 (1996).

[39] M. C. W. van Rossum and Th. M. Nieuwenhuizen, Rev. Mod. Phys. **71**, 313 (1999).

[40] S. J. Rahi, T. Emig, N. Graham, R. L. Jaffe, and M. Kardar, Phys. Rev. D **80**, 085021 (2009).

[41] A. Canaguier-Durand, P. A. Maia Neto, A. Lambrecht, and S. Reynaud, Phys. Rev. A **82**, 012511 (2010).

[42] R. C. Wittmann, IEEE Trans. Antennas Propag. **36**, 1078 (1988).

[43] E. D. Palik, *Handbook of Optical Constants of Solids* (Academic, New York, 1997).

[44] A.-F. Bitbol, A. Canaguier-Durand, A. Lambrecht, and S. Reynaud, Phys. Rev. B **87**, 045413 (2013).

[45] E. A. Power and T. Thirunamachandran, Phys. Rev. A **48**, 4761 (1993); **50**, 3929 (1994).

[46] I. M. Lifshitz, Adv. Phys. **13**, 483 (1964).

[47] T. L. Ferrell and R. H. Ritchie, Phys. Rev. A **21**, 1305 (1980).

[48] G. Barton, New J. Phys. **12**, 113045 (2010).

[49] G. V. Dedkov and A. A. Kyasov, Surf. Sci. **605**, 1077 (2011).
### 2.3.3 ForcaG experiment

In ref. [31], included hereafter, we have studied the effect of a material surface on the Wannier-Stark states for atoms trapped in an optical lattice. This work was done in the context of the ForcaG experiment at SYRTE which aims at measuring the Casimir-Polder and the gravitational force between atoms and a material surface using atomic interferometry techniques.

Wannier-Stark states are the solutions of a Hamiltonian consisting of a periodic potential plus a constant force. In the experiment, rubidium atoms are trapped in a vertical optical lattice. The periodic part of the Hamiltonian is then the optical potential $V_{opt}(z) = U(1 - \cos 2k_l z)/2$ created by the trapping laser and the constant force part of the Hamiltonian is the gravity force deriving from the gravitational potential $V_g(z) = -mgz$. Above, the quantity $U$ is the optical lattice depth. Depending on it, the Wannier-Stark states can be very localized on each lattice sites or spread out over multiple lattice sites.

In addition to the gravity and the trapping optical potential, the atoms are also subjected to the potential due to the presence of the mirror used to reflect the trapping laser. This surface potential consists of the Casimir-Polder potential and a short range part responsible for possible adsorption of atoms onto the surface which we model by a Lennard-Jones potential. Our work was to study the effects of this surface potential on the Wannier-Stark states of the atoms. The trapping laser has a wavelength of $\lambda_l = 532$ nm so that the lattice sites are located at a multiple of 266 nm from the mirror. The mirror is a Bragg mirror consisting of alternating layers of $SiO_2$ and $Ta_2O_5$ designed to be nearly transparent at 780 nm and 1064 nm while being reflecting at 532 nm. In order to calculate the Casimir-Polder potential, we calculated the reflection coefficient of this mirror using transfer matrix theory.

Data on the adsorption of rubidium atoms on the surface of the Bragg mirror are unknown. So, the parameters of the short range Lennard-Jones potential are unknown. We therefore use reasonable parameters for typical adsorption energies. The short range Lennard-Jones potential supports a number of surface bound states.

We found that the two least bound surface states lie in the energy range relevant to the Wannier-Stark states and strongly modify the Wannier-Stark spectrum. According to the LeRoy-Bernstein law, the position of those least bound states only depends on the depth $D$ of the Lennard-Jones potential (not its exact shape) and the $C_3$ coefficient of the non-retarded part of the Casimir-Polder potential. While the $C_3$ coefficient can be reliably determined from our calculations, the depth $D$ stays unknown. For that reason, we cannot draw any further conclusions concerning the surface-modified Wannier-Stark spectrum.

Nevertheless, the Wannier-Stark states localized far from the surface are basically insensitive to the unknown parameters of the Lennard-Jones potential. Using those states, it is possible to perform Raman interferometry and we have simulated a typical Raman spectrum. The effect of the Casimir-Polder potential between the atoms and the mirror is noticeable in the Raman spectrum and we concluded that the determination of the $C_3$ coefficient would be possible from a Raman interferometry experiment.
Surface-modified Wannier-Stark states in a one-dimensional optical lattice

A. Maury,\* M. Donaire, M.-P. Gorza,† A. Lambrecht, and R. Guérout

Laboratoire Kastler Brossel, UPMC-Sorbonne Universités, CNRS, ENS-PSL-Research University, Collège de France;
4 place Jussieu, F-75005 Paris, France

(Received 27 June 2016; published 2 November 2016)

We study the energy spectrum of atoms trapped in a vertical one-dimensional optical lattice in close proximity to a reflective surface. We propose an effective model to describe the interaction between the atoms and the surface at any distance. Our model includes the long-range Casimir-Polder potential together with a short-range Lennard-Jones potential, which are considered nonperturbatively with respect to the optical lattice potential. We find an intricate energy spectrum which contains a pair of loosely bound states localized close to the surface in addition to a surface-modified Wannier-Stark ladder at long distances. Atomic interferometry involving those loosely bound atom-surface states is proposed to probe the adsorption dynamics of atoms on mirrors.

DOI: 10.1103/PhysRevA.94.053602

I. INTRODUCTION

Trapping and manipulating cold neutral atoms in an optical lattice offers high control over the atomic locations and robust quantum coherence on the dynamics of the atomic states. These properties make of an optical lattice an ideal system for applications in metrology [1,2] and in precision measurements of the interactions between the atoms and the environment [3]. It is to the latter that the FORCA-G project applies [4]. In particular, the FORCA-G experiment aims at performing high-precision measurements of the electromagnetic and gravitational interactions between a neutral atom and a massive dielectric surface. Ultimately, it aims at establishing accurate constraints in the search of hypothetical deviations from the Newtonian law of gravitation at short length scales, the reason why an accurate knowledge of the electromagnetic interaction is also needed. It is on the electromagnetic interaction that we concentrate in this article.

In the setup of FORCA-G atoms of $^{87}$Rb are trapped in a vertical optical lattice by the potential generated by the standing waves of a laser source of wavelength $\lambda_I = 532$ nm, which reflect off a Bragg mirror (see Fig. 1). The optical potential takes the periodic form

$$V_{op}(z) = U(1 - \cos 2k_l z)/2, \quad (1)$$

where $k_l = 2\pi/\lambda_l$, $z$ is the vertical distance relative to the surface position and $U$ is the optical depth which depends on the laser intensity. In addition, the uniform Earth gravitational field creates a linear potential

$$V_g(z) = -mgz \quad (2)$$

with $m$ being the atomic mass and $g$ being the gravitational acceleration. Disregarding the atom-mirror interaction, the spectrum which results from the addition of the optical and gravitational potentials consists of a ladder of quasistable states referred to as Wannier-Stark (WS) states. The WS eigenstates are localized around the equilibrium points $z_n = n\lambda_l/2$, $n$ being an integer, and are uniformly distributed along the energy

spectrum at constant intervals $mg\lambda_l/2 = h\nu_B$. In this expression $\nu_B$ is referred to as Bloch frequency, and the degree of localization is determined by the relative optical depth with respect to the recoil energy, $U/(\hbar^2k_l^2/2m) = U/E_r$ (see Fig. 2).

In addition to $V_{op}(z)$ and $V_g(z)$, the neutral atoms interact with the surface through the mutual coupling of their charge fluctuations to the vacuum fluctuations of the electromagnetic field. This interaction is known generically as Casimir-Polder (CP) interaction [5,6]. At zero temperature its strength depends generally on the dielectric properties of the surface, the state of the atom, and the distance between them.

The *modus operandi* of FORCA-G consists of a sequence of pulses generated by Raman lasers and microwaves which are used to create an atomic interferometer. The pulses drive the atoms through a coherent superposition of low-lying Zeeman sublevels at different lattice sites [4]. The CP interaction induces a phase shift on the atoms which depends strongly on the distance of the atoms to the surface and slightly on the internal atomic states. The phase shift accumulated by the atomic wave function throughout the sequence of pulses is finally measured by atomic interferometry techniques. If the atoms are made to oscillate between lattice sites far from the surface, the CP-induced shift is additive. Therefore, once the phase shift associated to the passage through different WS levels, which is characteristic of the interferometer scheme, is subtracted, the remaining phase is the CP-induced shift we are interested in.

The latter applies to the case where the CP interaction is small with respect to the optical potential depth, so that it can be treated as a perturbation to the potential $V_{op}(z) + V_g(z)$ and hence to the original WS states. This takes place at separation distances of the order of microns, at which the perturbative development of Refs. [7,8] applies. On the contrary, at submicrometer distances and beyond the perturbative regime, it was already noticed in Ref. [7] that the CP corrections to the original WS energies diverge. This is especially relevant to the purposes of the FORCA-G project, as deviations from Newtonian gravity are expected to occur at submicrometer distances. Therefore, a precise knowledge of the CP interaction at this length scale as well as an accurate description of the spatial distribution of the atomic wave function are crucial in order to detect those gravitational effects. In Ref. [7] the authors apply a regularization scheme for the CP potential

<footnote>\*maury@lkb.upmc.fr</footnote>

<footnote>†on leave from Laboratoire de physique des lasers, Université Paris 13, Villetaneuse, France.</footnote>
<figure><img src="image_28.png" /><figcaption>FIG. 1. Scheme of the experimental device. Cold $^{87}$Rb atoms are trapped in a blue-detuned vertical optical lattice. An infrared laser at $\lambda = 1064$ nm ensures the transverse confinement. A pair of contrapropagating Raman lasers at $\lambda = 780$ nm (not shown) drives the transitions between lattice sites.</figcaption></figure>

based on the assumption that the minimum distance of the atom to the surface is determined by the atomic radius. However, it is found there that the resultant corrections strongly depend on this radius. Thus, nonreliable results were obtained.

It is the main purpose of the present article to develop a nonperturbative approach to this problem in order to determine accurately the energy spectrum and the profile of the atomic states at submicrometer distances. To this end, we model the short-range interaction between the atom and the surface by a Lennard-Jones potential which features the adsorption of the atoms on the surface. We find that, in addition to slightly modified WS states, the resultant spectrum contains a number of loosely bound atom-surface states whose properties depend critically on the parameters of the Lennard-Jones potential. Nonetheless, independent measurements can be performed to determine the unknowns of such potential.

The remainder of the article is organized as follows. In Sec. II, we present the features of the potential modeling the interaction between the atom and the surface. In Sec. III we show that the overall effect of the surface leads to a complex energy spectrum significantly departing from the usual Wannier-Stark states. We conclude by calculating a typical atomic interferometry spectrum obtained using stimulated Raman transitions between those surface-modified Wannier-Stark states.

## II. THE ATOM-SURFACE POTENTIAL

In addition to the optical potential described in the precedent section, the atoms interact with the mirror through

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /><figcaption>FIG. 2. Profile of the squared-norm of the wave function of the nth WS state for different values of the optical depth, $U = 3 E_r$ (upper figure) and $U = 9 E_r$ (lower figure).</figcaption></figure>

the electromagnetic field. Quite generally, this interaction is made of two distinct components, a short-range and a long-range potentials. The short-range potential results from the spatial overlap between the electronic clouds of the atoms and the surface at subnanometer distances. In turn, this potential depends on the precise profile of the electronic density distribution, which is difficult to determine both experimentally and theoretically. Hence, a parametrization scheme is required for the short-range potential. In contrast, the long-range potential originates from the mutual coupling of the charges within the atoms and the currents in the mirror to the fluctuating electromagnetic field. This is the so-called Casimir-Polder potential, which is computed in the electric dipole approximation at second order in stationary perturbation theory.

In the framework of the scattering theory [9], the Casimir-Polder potential between a flat mirror in the (xy) plane and an atom in the ground state separated by a distance $z$, at temperature $T$, is given by [10]

$$V_s^{\text{CP}}(z) = k_B T \sum_n \frac{\xi_n^2}{c^2} \frac{\alpha(i\xi_n)}{4\pi\epsilon_0} \int_0^\infty \frac{d^2\mathbf{k}}{\kappa} e^{-2\kappa z} \\ \times \left[ \rho^{\text{TE}} - \left(1 + \frac{2\kappa^2 c^2}{\xi_n^2}\right) \rho^{\text{TM}} \right] \quad (3)$$
with $k^2 = k_s^2 + k_i^2$, $\kappa = \sqrt{k^2 + \xi_n^2/c^2}$ and the sum runs over Matsubara frequencies $\xi_n = 2\pi nk_B T/\hbar$. In this equation $\rho^{\text{TE}}$ and $\rho^{\text{TM}}$ are the reflection coefficients of the mirror for the TE and TM polarizations, respectively, and $\alpha(i\xi)$ is the polarizability of a $^{87}$Rb atom in its ground state evaluated at imaginary frequencies [11],

$$ \alpha(i\xi) = \frac{2}{\hbar} \sum_j \frac{\omega_{j\xi} d_{j\xi}^2}{\omega_{j\xi}^2 + \xi^2}, \quad (4) $$

where $\omega_{j\xi} = \omega_j - \omega_g$ and $d_{j\xi}$ are, respectively, the transition frequency and the electric dipole matrix element between the states $|j\rangle$ and $|g\rangle$.

Concerning the optical properties of the mirror used in the FORCA-G experiment, its design is such that it is nearly transparent at 780 nm and 1064 nm while it is reflective at 532 nm. It is a Bragg mirror formed by alternating layers of SiO$_2$ and Ta$_2$O$_5$. Its reflection coefficients $\rho^{\text{TE}}$ and $\rho^{\text{TM}}$ are obtained using standard transfer matrix theory. Let us define first by $T_i$ the transfer matrix associated to the transmission through the interface between the layers $i$ and $i+1$, as well as to the propagation throughout the layer $i+1$ of width $d_{i+1}$. It relates the field on the left of the layer $i$ to the field on the right and reads [12]

$$ T_i = \frac{1}{\bar{t}_i} \begin{pmatrix} t_i \bar{t}_i - r_i \bar{r}_i & \bar{r}_i \\ -r_i & 1 \end{pmatrix} \begin{pmatrix} e^{ik_i d_{i+1}} & 0 \\ 0 & e^{-ik_i d_{i+1}} \end{pmatrix}. \quad (5) $$

In this equation, $r_i$ and $t_i$ are the Fresnel amplitudes from medium $i$ to medium $i+1$. The barred quantities are the reciprocal amplitudes from medium $i+1$ to medium $i$ and $k_z$ is the z component of the wave vector in medium $i+1$. The transfer matrix of the Bragg mirror, $\mathbb{T}$, is the product of the transfer matrices of all the layers $\mathbb{T} = \prod_i T_i$ and reads

$$ \mathbb{T} = \frac{1}{\bar{t}} \begin{pmatrix} \tau \bar{t} - \rho \bar{\rho} & \bar{\rho} \\ -\rho & 1 \end{pmatrix}, \quad (6) $$

from which the total reflection amplitude reads $\rho = -[\mathbb{T}]_{21}/[\mathbb{T}]_{22}$.

We show in Fig. 3 the Casimir-Polder potential calculated using Eq. (3) for a temperature $T = 300$ K. The potential is scaled with $z^3$, the third power of the atom-surface distance, in order to emphasize the nonretarded van der Waals regime characterized by its coefficient $C_3 \approx 3.28 a_0^3$ eV.

As for the short-range potential, we parametrize it using a 12-3 Lennard-Jones form,<fnref n="1" />

$$ V_s^{\text{LJ}}(z) = \frac{D}{3} \left[ \left( \frac{z_0}{z} \right)^{12} - 4 \left( \frac{z_0}{z} \right)^3 \right], \quad (7) $$

which is characterized by a well depth $D$ and an equilibrium distance $z_0$ which correspond to the energy and distance from the surface of an adsorbed atom, respectively. Continuity of the atom-surface potential demands that $V_s^{\text{LJ}}(z)$ and $V_s^{\text{CP}}(z)$ smoothly merge at some intermediate distance $z_m$. This implies

<figure><img src="image_31.png" /><figcaption>FIG. 3. The Casimir-Polder potential calculated between a Rubidium atom and a SiO$_2$-Ta$_2$O$_5$ Bragg mirror as a function of the distance $z$. The value of the van der Waals coefficient $C_3$ is indicated.</figcaption></figure>

that $D$ and $z_0$ are no longer independent but are related by the equation

$$ \frac{4}{3} D z_0^3 = C_3, \quad (8) $$

where $C_3$ is the van der Waals coefficient in the Casimir-Polder potential. With this condition between the parameters $D$ and $z_0$ in the Lennard-Jones potential, the matching distance $z_m$ is chosen where both potentials $V_s^{\text{LJ}}(z)$ and $V_s^{\text{CP}}(z)$ behave in $z^{-3}$ and leads to the total surface potential $V_s(z)$:

$$ V_s(z) = V_s^{\text{LJ}}(z)\Theta(z_m - z) + V_s^{\text{CP}}(z)\Theta(z - z_m), \quad (9) $$

where $\Theta(z)$ is the Heaviside function.

The form used for $V_s^{\text{LJ}}(z)$ is merely of a physisorption type and hence is expected to underestimate the adsorption energy. For instance, for an equilibrium distance $z_0 = 2.3$ Å we find $D \approx 30$ meV to be compared with a value of $\approx 350$ meV from a recent density functional theory calculation [13]. As a matter of fact, the parameters of the short-range potential carry the largest uncertainty in our calculation. An accurate determination of this part of the potential would require extensive ab initio calculations up to distances of the order of the nanometers which are beyond the scope of this work. Alternatively, the parameters $D$ and $z_0$ can be determined experimentally. Be that as it may, we will study in the next section the influence of our results upon the parameters of the Lennard-Jones model.

### III. SURFACE-MODIFIED WANNIER-STARK STATES

In the following and unless otherwise stated, we will refer to the distance $z$ to the surface in units of $\lambda_l/2 = 266$ nm and the energies in units of the recoil energy $E_r = \frac{\hbar^2 k_l^2}{2m} \approx 5.37 \times 10^{-30}$ J for a Rubidium atom. The surface-modified

<footnote n="1">¹In surface science, a 9-3 Lennard-Jones potential is also often used as it arises as pairwise summation of 12-6 Lennard-Jones atom-atom interactions.</footnote>
<figure><img src="image_32.png" /><figcaption>FIG. 4. Potential $V(z)$ in units of the recoil energy $E_r$ shown as the black curve. The dashed blue and red curves are, respectively, the surface potential $V_s(z)$ and the gravitational potential $-mgz$.</figcaption></figure>

Wannier-Stark states (SMWSSs) are solutions of the time-independent Schrödinger equation

$$ \frac{\hbar^2}{2m} \frac{d^2\psi_n(z)}{dz^2} + V(z)\psi_n(z) = E_n\psi_n(z), \quad (10) $$

$$ \text{with } V(z) = V_s(z) + V_g(z) + V_{op}(z). \quad (11) $$

In the situation where the mirror is above the atoms, the potential $V(z)$ is not bounded from below so that all states are rigorously Siegert states [14]. That corresponds to the situation where the atoms could ultimately “fall from the optical lattice.” Nevertheless, it has been shown in Ref. [8] that lifetimes are of the order of $10^{14}$ s for the first Bloch band, and hence they can be considered stable for any experimental realization. We show in Fig. 4 the potential $V(z)$ for an optical depth $U = 3E_r$.

At $z \approx 2$, the magnitudes of the gravitational and the Casimir-Polder potentials are similar. As a result, the very first optical well is strongly influenced by the surface to the point of becoming weakly bounding. Note that the minimum of the Lennard-Jones part of the surface potential has very different orders of magnitude, both in binding energy ($D \approx 10^9 E_r$) and in equilibrium distance ($z_0 = 2.3 \AA \approx 10^{-3} \lambda_{el}/2$), a reason why it does not appear in Fig. 4.

The SMWSSs $\psi_n(z)$ are conveniently characterized according to their mean distance to the surface $\langle z \rangle$:

$$ \langle z \rangle = \frac{\langle \psi_n | z | \psi_n \rangle}{\langle \psi_n | \psi_n \rangle}. \quad (12) $$

We show in Table I values of the mean distance $\langle z \rangle$ and the energy intervals for the first few SMWSSs calculated for an optical depth $U = 3 E_r$ and $z_0 = 2.3 \AA$, ordered according to an increasing value of $\langle z \rangle$ (the first excited Bloch band corresponds to energies greater than the optical depth $U$ and

<table><caption>TABLE I. SMWSSs for a lattice depth $U = 3 E_r$ ordered according to their mean distance to the surface $\langle z \rangle$. Energy intervals are given in units of $E_r$. Further analysis (see text) shows that surface-modified Wannier-Stark states begin at $n = 3$, while the first two states are atom-surface bound states. The last column refers to the energy intervals of an infinite potential surface (i.e., perfect surface).</caption><thead><tr><th>n</th><th>⟨z⟩</th><th>E<sub>n</sub> - E<sub>n-1</sub></th><th>Perfect surface</th></tr></thead><tbody><tr><td>1</td><td>0.799</td><td>E<sub>1</sub> = -0.0709</td><td></td></tr><tr><td>2</td><td>1.006</td><td>+1.9690</td><td></td></tr><tr><td>3</td><td>3.372</td><td>-0.5468</td><td></td></tr><tr><td>4</td><td>4.268</td><td>-0.1264</td><td>-0.1371</td></tr><tr><td>5</td><td>4.681</td><td>-0.0934</td><td>-0.0996</td></tr><tr><td>6</td><td>4.746</td><td>-0.0693</td><td>-0.0804</td></tr><tr><td>7</td><td>5.617</td><td>-0.0579</td><td>-0.0722</td></tr><tr><td>8</td><td>6.881</td><td>-0.0637</td><td>-0.0703</td></tr><tr><td>9</td><td>7.962</td><td>-0.0679</td><td>-0.0701</td></tr><tr><td>10</td><td>8.985</td><td>-0.0692</td><td>-0.0701</td></tr><tr><td>11</td><td>9.994</td><td>-0.0696</td><td>-0.0701</td></tr><tr><td>12</td><td>10.998</td><td>-0.0698</td><td>-0.0701</td></tr><tr><td>13</td><td>12.001</td><td>-0.0700</td><td>-0.0701</td></tr><tr><td>14</td><td>13.002</td><td>-0.0700</td><td>-0.0701</td></tr><tr><td>15</td><td>14.003</td><td>-0.0700</td><td>-0.0701</td></tr><tr><td>16</td><td>15.003</td><td>-0.0701</td><td>-0.0701</td></tr><tr><td>17</td><td>16.003</td><td>-0.0701</td><td>-0.0701</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

is therefore not trapped). The closest SMWSSs are modified very strongly by the presence of the surface, which reflects on the lack of regularity characteristic of a Wannier-Stark ladder. As $\langle z \rangle$ increases, though, we progressively recover a usual Wannier-Stark ladder spaced by the Bloch energy $h\nu_B$ and integer values of $\langle z \rangle$.

For the purpose of the FORCA-G experiment, we are mostly interested in the states closest to the surface.

We show in Fig. 5 the profile of the real wave functions corresponding to the first four SMWSSs according to Table I. The probability amplitudes of the first two states exhibit very

<figure><img src="image_33.png" /><figcaption>FIG. 5. Wave functions for the first four SMWSSs $\psi_n(z)$ according to Table I. As it is customary, the vertical offset of the wave functions correspond to their respective energies.</figcaption></figure>
rapid oscillations within the Lennard-Jones well, while they
are vanishingly small outside this well. On the other hand,
the states $n = 3$ and $n = 4$ already are well spread along the
optical potential as the ordinary Wannier-Stark states would.
The tail of their wave functions still shows some oscillations
caused by the Lennard-Jones potential.

A. Dependence upon the Lennard-Jones parameters

While modeling the short-range potential, our largest
uncertainty lies in the unknown shape of the potential well.
Although we have used a known analytical form which
correctly converges towards the Casimir-Polder potential, the
actual short-range potential may differ from the Lennard-Jones
form [13]. It is therefore crucial to study the dependence of our
results upon the free parameters in $V_s^{\text{LJ}}(z)$. Let us recall some
general features of the bound states of a 12 – 3 Lennard-Jones
potential. Having a finite depth $D$ and vanishing sufficiently
fast as $z \to \infty$, the potential $V_s^{\text{LJ}}(z)$ given by Eq. (7) possesses
a finite number of bound states. Those states represent
vibrational states for an atom bound to the surface and are
therefore indexed with an integer vibrational quantum number
$v$ starting with $v = 0$ for the ground state. When the total
number of bound states supported by a potential well is
unknown (e.g., due to uncertainties on the dissociation energy
$D$) it is customary to label the least bound states as $v = -1$,
the second least bounded states as $v = -2$, and so on. To a
very good approximation, the position of the few least bound
states depends only on the asymptotic behavior of the potential
as $z \to \infty$ and on a noninteger effective vibrational quantum
number at dissociation, $v_D$, which varies between 0 and -1
[15]. By decreasing continuously the depth of the potential
the states $v = -1$, $v = -2$ will be eventually expelled to the
continuum. From those considerations we see that, as far as
the few least bound states are concerned, the exact shape of the
potential energy well is not important. In our case, the effective
vibrational quantum number $v_D$ can be varied by simply
decreasing the depth $D$ of our 12 – 3 Lennard-Jones model.
Owing to Eq. (8), the dissociation energy $D$ is decreased
by increasing the equilibrium atom-surface distance $z_0$ as
$D(\text{eV}) \approx 0.36z_0^{-3}$ (Å).

We show in Fig. 6 the energies of the SMWSSs as a function of $z_0$ or, equivalently, as a function of decreasing dissociation energy $D$. In the first place, one sees that the states $n = 1$ and $n = 2$ have a very different behavior compared to all the others. The position of those states depends critically upon the dissociation energy $D$. As such, it is clear that the two SMWSS states $n = 1$ and $n = 2$ are basically the last two bound states $v = -2$ and $v = -1$, respectively, of the short-range potential. As the equilibrium distance $z_0$ increase, the energies of the states $n = 1$ and $n = 2$ increases, and they cross all the other states. Nonetheless there must be avoided crossings since all those states result from the diagonalization of the Hamiltonian operator.

On the other hand, the energies of the states starting from
the $n = 3$ are very much independent of the parameters used
in the short-range surface potential $V_s^{\text{LJ}}(z)$ except near an
avoided crossing with a bound atom-surface state. From Fig. 6
we conclude that the state $n = 3$ can be considered the first
surface-modified Wannier-Stark states. The coupling between

<figure><img src="image_34.png" /><figcaption>FIG. 6. Calculated energies of the SMWSSs as a function of the distance $z_0$ in the Lennard-Jones potential, at constant $C_3$ coefficient. The first four states depicted in Fig. 5 and calculated for $z_0 = 2.3$ Å are indicated by arrows.</figcaption></figure>

the few first SMWSSs and the atom-surface bound states
quickly vanishes as $n$ increases owing to the vastly different
mean atom-surface distance $\langle z \rangle$. This leads to negligible
avoided crossings between the state $n = 2$ and already the
state $n = 7$. Far from any avoided crossings, the SMWSSs are
still influenced by the surface. At $z_0 = 2.3$ Å it is shown in
Table I that the energy intervals between successive states are
not equal to the Bloch frequency for the first Wannier-Stark
states.

It is also illustrative to compare our results with those
obtained from the modeling of the short-range potential with
that of a perfect surface:

$$
V_s(z) = \begin{cases} +\infty & z < 0 \\ 0 & z > 0 \end{cases} \tag{13}
$$

In the first place, the repulsive part of the Lennard-Jones
potential plays the role of an infinite potential wall. However,
in the case of an infinite potential surface the wave function has
a different behavior at $z = 0$. In particular, the wave function
vanishes monotonically as $z \to 0$ [7,8], whereas it oscillates very
rapidly within a Lennard-Jones potential. Obviously, a major
drawback of an infinite potential surface is the total absence of
bound atom-surface states. Values of the corresponding energy
intervals can be found in Table I.

B. Simulated Raman spectrum

The experimental setup of the FORCA-G is detailed, e.g., in
Ref. [16]. In it, two counterpropagating Raman lasers operating
at $\lambda = 780$ nm drive coherent transitions between the ground
$|5^2S_{1/2}, F = 1, m_F = 0\rangle$ and excited $|5^2S_{1/2}, F = 2, m_F = 0\rangle$
hyperfine levels of trapped $^{87}$Rb atoms. Those transitions can
involve different SMWSSs with a probability proportional to
<figure><img src="image_35.png" /><figcaption>FIG. 7. Raman transition probabilities between an initial SMWSS $\psi_n(z)$ and a final state $\psi_m(z)$.</figcaption></figure>

the generator of translations along the z axis, $\langle\psi_n|e^{ik_{eff}z}|\psi_m\rangle$, with $k_{eff} \approx 4\pi/(780 \text{ nm})$.

We show in Fig. 7 the Raman transition probabilities between the states presented in Table I. The first two states, $\psi_1(z)$ and $\psi_2(z)$, which are the atom-surface bound states are only weakly coupled to the surface-modified Wannier-Stark states but strongly coupled to each other. We can see the smooth evolution of the SMWSSs towards “regular,” unmodified Wannier-Stark states whose transition probabilities become a function of $|n-m|$ only. For a lattice depth of 3 $E_r$, a given state $\psi_n(z)$ roughly couples to states up to $n \pm 6$.

With a low-density atomic cloud like in Ref. [16], some 10⁴ lattice sites are populated, and the Raman spectrum is dominated by transitions involving unmodified Wannier-

<figure><img src="image_36.png" /><figcaption>FIG. 9. Change in Raman transition frequencies $v_{n \to m}$, as a function of the van der Waals coefficient $C_3$, for selected states: (a) $|m-n|=1$ transitions; (b) $|m-n|=3$ transitions; (c) $|m-n|=5$ transitions.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>FIG. 8. Raman stick spectrum involving the states in Table I. Lines involving the atom-surface bound states $\psi_1(z)$ and $\psi_2(z)$ are not shown.</figcaption></figure>

Stark states. When the frequency difference between the two Raman lasers, $v_R = v_{R1} - v_{R2}$, is scanned around the rubidium hyperfine splitting $v_{HFS}$, this leads to a simple spectrum with lines at integer numbers of the Bloch frequency $v_B = h^{-1}mg\lambda_l/2 \approx 568.5$ Hz. One could imagine an experiment with a much more dense atomic sample with a size of a few microns where the contribution from the SMWSSs would be visible.

We show in Fig. 8 the simulated Raman stick spectrum (spectrum without line shapes) for the states listed in Table I. As we have shown in Fig. 6, the position of the atom-surface bound states $\psi_1(z)$ and $\psi_2(z)$ is largely unknown. Therefore, we do not show their contributions in the spectrum of Fig. 8. The energies of those atom-surface bound states will appear as additional lines in the Raman spectrum. As expected, the departure from the regular Wannier-Stark ladder generates many lines. Those lines have the tendency to bundle up around integer numbers of the Bloch frequency though. Recently a
relative sensibility of $4 \times 10^{-6}$ at 1 s on the measure of the Bloch frequency has been demonstrated using a Ramsey-type interferometry [17]. Such a sensibility would in principle allow to resolve the lines presented in Fig. 8.

### C. Determination of the Casimir-Polder potential

Up to now, the Casimir-Polder potential has been kept constant to its calculated value in Sec. II. The aim of the FORCA-G experiment is to determine this Casimir-Polder potential from a recorded Raman spectrum. Thus, we have to know how the Raman spectrum changes when one changes the Casimir-polder potential. For other references related to the use of Bloch oscillations in order to measure the coefficients in the Casimir-Polder potential, see, e.g., Refs. [18,19]. In the following, we focus on the van der Waals coefficient $C_3$.

We show in Fig. 9 the change in Raman transition frequencies $v_{n \to m} = h^{-1}(E_n - E_m)$ when the $C_3$ coefficient is allowed to vary from its nominal value of 3.28 $a_0^3$ eV calculated in Sec. II. We present in Fig. 9 selected transitions involving $|n-m|=1, 3$, and 5 and selected states $n \ge 6$. A precise analysis of the position of those lines with respect to integer values of the Bloch frequency will allow the van der Waals coefficient $C_3$ to be determined. In fact, from an absolute uncertainty of 20 mHz [17] on the determination of

Raman transition frequencies, we infer a relative uncertainty $\delta C_3/C_3$ on the van der Waals coefficient ranging from $10^{-2}$ to $10^{-4}$.

### IV. CONCLUSION

We have calculated the energies of atoms trapped in a one-dimensional vertical optical lattice taking into account the interaction between those atoms and the mirror used to realize the lattice. We have found that, in the range of energy of a few recoil energy $E_r$, loosely bound atom-mirror states appear as additional levels among an otherwise surface-modified Wannier-Stark ladder. The energies of those loosely bound atom-mirror states depend critically on the details of the adsorption atom-surface potential. Atomic interferometry involving those loosely bound atom-mirror states will shed light on the adsorption dynamics of rubidium atoms on mirrors. The close surface-modified Wannier-Stark states correspond to optically trapped atoms which nevertheless have a significant probability of being adsorbed by the mirror.

### ACKNOWLEDGMENTS

We thank Franck Pereira dos Santos and Peter Wolf for stimulating discussions. Financial support from Agence Nationale de la Recherche ANR-10-IDEX-0001-02-PSL and ANR-13-BS04-0003-02 is gratefully acknowledged.

[1] P. Cladé, E. de Mirandes, M. Cadoret, S. Guellati-Khélifa, C. Schwob, F. Nez, L. Julien, and F. Biraben, Determination of the Fine Structure Constant Based on Bloch Oscillations of Ultracold Atoms in a Vertical Optical Lattice, Phys. Rev. Lett. 96, 033001 (2006).

[2] M. Takamoto, F.-L. Hong, R. Higashi, and H. Katori, An optical lattice clock, Nature (London) 435, 321 (2005).

[3] F. Impens, C. Ccapa Tíra, R. O. Behunin, and P. A. Maia Neto, Dynamical local and nonlocal Casimir atomic phases, Phys. Rev. A 89, 022516 (2014).

[4] P. Wolf, P. Lemonde, A. Lambrecht, S. Bize, A. Landragin, and A. Clairon, From optical lattice clocks to the measurement of forces in the Casimir regime, Phys. Rev. A 75, 063608 (2007).

[5] H. B. G. Casimir and D. Polder, The influence of retardation on the London-van der Waals forces, Phys. Rev. 73, 360 (1948).

[6] E. M. Lifschitz, The theory of molecular attractive forces between solids, Sov. Phys. 2, 73 (1956).

[7] R. Messina, S. Pelisson, M.-C. Angonin, and P. Wolf, Atomic states in optical traps near a planar surface, Phys. Rev. A 83, 052111 (2011).

[8] S. Pelisson, R. Messina, M. C. Angonin, and P. Wolf, Lifetimes of atoms trapped in an optical lattice in proximity of a surface, Phys. Rev. A 88, 013411 (2013).

[9] A. Lambrecht, P. A. Maia Neto, and S. Reynaud, The Casimir effect within scattering theory, New J. Phys. 8, 243 (2006).

[10] G. Dufour, A. Gérardin, R. Guérout, A. Lambrecht, V. V. Nesvizhevsky, S. Reynaud, and A. Yu Voronin, Quantum reflection of antihydrogen from Casimir potential above matter slabs, Phys. Rev. A 87, 012901 (2013).

[11] A. Derevianko, S. G. Porsev, and J. F. Babb, Electric dipole polarizabilities at imaginary frequencies for hydrogen, the alkali-metal, alkaline-earth and noble gas atoms, At. Data Nucl. Data Tables 96, 323 (2010).

[12] G.-L. Ingold and A. Lambrecht, Casimir effect from a scattering approach, Am. J. Phys. 83, 156 (2015).

[13] J. A. Sedlacek, E. Kim, S. T. Rittenhouse, P. F. Weck, H. R. Sadeghpour, and J. P. Shaffer, Electric Field Cancellation on Quartz by Rb Adsorbate-Induced Negative Electron Affinity, Phys. Rev. Lett. 116, 133201 (2016).

[14] A. J. F. Siegert, On the derivation of the dispersion relation for nuclear reactions, Phys. Rev. 56, 750 (1939).

[15] R. J. LeRoy and R. B. Bernstein, Dissociation energy and long-range potential of diatomic molecules from vibrational spacings of higher levels, J. Chem. Phys. 52, 3869 (1970).

[16] Q. Beaufils, G. Tackmann, X. Wang, B. Pelle, S. Pelisson, P. Wolf, and F. Pereira dos Santos, Laser Controlled Tunneling in a Vertical Optical Lattice, Phys. Rev. Lett. 106, 213002 (2011).

[17] A. Hilico, C. Solaro, M.-K. Zhou, M. Lopez, and F. Pereira dos Santos, Contrast decay in a trapped-atom interferometer, Phys. Rev. A 91, 053616 (2015).

[18] I. Carusotto, L. Pitaevskii, S. Stringari, G. Modugno, and M. Inguscio, Sensitive Measurement of Forces at the Micron Scale using Bloch Oscillations of Ultracold Atoms, Phys. Rev. Lett. 95, 093202 (2005).

[19] F. Sorrentino, A. Alberti, G. Ferrari, V. V. Ivanov, N. Poli, M. Schioppo, and G. M. Tino, Quantum sensor for atom-surface interactions below $10 \mu m$, Phys. Rev. A 79, 013409 (2009).
## 3 Conclusion

We have presented our method of calculating the Casimir interaction energy between an arbitrary number of bodies based on the scattering theory. This method can be applied to any object whose scattering operator is known and do not require any regularization method leading automatically to a finite interaction energy. The comparison between theory and experiments in the case of metallic bodies is still not satisfactory regarding the role of the dissipation in the model of permittivity used in the calculations. We hope that further experiments will permit to shed some light on this puzzle.

In the recent years, our group has been involved in the GBAR experiment which takes place at CERN. This experiment aims at measuring the gravitational acceleration acting on antimatter as it falls. This has been the subject of our last three thesis and has shifted the focus of our research from Casimir physics. We hope that the field of Casimir physics will stay active in the coming years and that the remaining discrepancies between theory and experiments will be finally resolved.

## References

[1] Hendrick BG Casimir. On the attraction between two perfectly conducting plates. In *Proc. Kon. Ned. Akad. Wet.*, volume 51, page 793, 1948.

[2] Julian Schwinger, Lester L DeRaad Jr, and Kimball A Milton. Casimir effect in dielectrics. *Annals of Physics*, 115(1):1-23, 1978.

[3] J Dalibard, J Dupont-Roc, and C Cohen-Tannoudji. Vacuum fluctuations and radiation reaction: identification of their respective contributions. *Journal de Physique*, 43(11):1617-1638, 1982.

[4] PW Milonni. Casimir forces without the vacuum radiation field. *Physical Review A*, 25(3):1315, 1982.

[5] Satofumi Souma and Akira Suzuki. Local density of states and scattering matrix in quasi-one-dimensional systems. *Physical Review B*, 65(11):115307, 2002.

[6] Romain Pierrat, Philipp Ambichl, Sylvain Gigan, Alexander Haber, Rémi Carminati, and Stefan Rotter. Invariance property of wave scattering through disordered media. *Proceedings of the National Academy of Sciences*, 111(50):17765-17770, 2014.

[7] Felix T Smith. Lifetime matrix in collision theory. *Physical Review*, 118(1):349, 1960.

[8] Amos Egel, Lorenzo Pattelli, Giacomo Mazzamuto, Diederik S Wiersma, and Uli Lemmer. Celes: Cuda-accelerated simulation of electromagnetic scattering by large ensembles of spheres. *Journal of Quantitative Spectroscopy and Radiative Transfer*, 199:103-110, 2017.
[9] JJMR Gough and MR James. Quantum feedback networks: Hamiltonian formulation. *Communications in Mathematical Physics*, 287(3):1109-1132, 2009.

[10] Robert Cook, David I Schuster, Andrew N Cleland, and Kurt Jacobs. Input-output theory for superconducting and photonic circuits that contain weak retroreflections and other weak pseudocavities. *Physical Review A*, 98(1):013801, 2018.

[11] D Jackson John. *Classical electrodynamics*. John Wiley & Sons, Inc., New York, NY, 1999.

[12] EM Lifshitz. The theory of molecular attractive forces between solids. *Soviet Physics*, 2(1):73, 1956.

[13] RS Decca, D López, E Fischbach, GL Klimchitskaya, DE Krause, and VM Mostepanenko. Tests of new physics from precise measurements of the casimir pressure between two gold-coated plates. *Physical Review D*, 75(7):077101, 2007.

[14] Mingyue Liu, Jun Xu, GL Klimchitskaya, VM Mostepanenko, and U Mohideen. Examining the casimir puzzle with an upgraded afm-based technique and advanced surface cleaning. *Physical Review B*, 100(8):081406, 2019.

[15] F Intravaia and A Lambrecht. Surface plasmon modes and the casimir energy. *Physical review letters*, 94(11):110404, 2005.

[16] GL Klimchitskaya and VM Mostepanenko. Theory-experiment comparison for the casimir force between metallic test bodies: A spatially nonlocal dielectric response. *Physical Review A*, 105(1):012805, 2022.

[17] Romain Guérout, J Lussange, Ho Bun Chan, Astrid Lambrecht, and Serge Reynaud. Thermal casimir force between nanostructured surfaces. *Physical Review A*, 87(5):052514, 2013.

[18] MG Moharam, Drew A Pommet, Eric B Grann, and TK Gaylord. Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach. *JOSA A*, 12(5):1077-1086, 1995.

[19] Y Bao, Romain Guérout, J Lussange, Astrid Lambrecht, Raymond A Cirelli, F Klemens, William M Mansfield, CS Pai, and Ho Bun Chan. Casimir force on a surface with shallow nanoscale corrugations: Geometry and finite conductivity effects. *Physical review letters*, 105(25):250402, 2010.

[20] Ho Bun Chan, Y Bao, J Zou, RA Cirelli, F Klemens, WM Mansfield, and CS Pai. Measurement of the casimir force between a gold sphere and a silicon surface with nanoscale trench arrays. *Physical review letters*, 101(3):030401, 2008.

[21] Fujio Shimizu. Specular reflection of very slow metastable neon atoms from a solid surface. *Physical review letters*, 86(6):987, 2001.
[22] Hilmar Oberst, Yoshihisa Tashiro, Kazuko Shimizu, and Fujio Shimizu. Quantum reflection of he* on silicon. *Physical Review A*, 71(5):052901, 2005.

[23] Riccardo Messina, Diego AR Dalvit, Paulo A Maia Neto, Astrid Lambrecht, and Serge Reynaud. Dispersive interactions between atoms and nonplanar surfaces. *Physical Review A*, 80(2):022119, 2009.

[24] A Selmani, DR Salahub, and A Yelon. Interaction potential between a hydrogen atom and a silicon (100) surface. *Surface science*, 202(1-2):269-276, 1988.

[25] Manuel Donaire, Romain Guérout, and Astrid Lambrecht. Quasiresonant van der waals interaction between nonidentical atoms. *Physical Review Letters*, 115(3):033201, 2015.

[26] Ulrich D Jentschura and V Debierre. Long-range tails in van der waals interactions of excited-state and ground-state atoms. *Physical Review A*, 95(4):042506, 2017.

[27] Manuel Donaire. Two-atom interaction energies with one atom in an excited state: van der waals potentials versus level shifts. *Physical Review A*, 93(5):052706, 2016.

[28] M Donaire. Net force on an asymmetrically excited two-atom system from vacuum fluctuations. *Physical Review A*, 94(6):062701, 2016.

[29] Nicolas Cherroret, Romain Guérout, Astrid Lambrecht, and Serge Reynaud. Statistical approach to casimir-polder potentials in heterogeneous media. *Physical Review A*, 92(4):042513, 2015.

[30] Ronald C Wittmann. Spherical wave operators and the translation formulas. *IEEE Transactions on Antennas and Propagation*, 36(8):1078-1087, 1988.

[31] A Maury, M Donaire, M-P Gorza, A Lambrecht, and R Guérout. Surface-modified wannier-stark states in a one-dimensional optical lattice. *Physical Review A*, 94(5):053602, 2016.