# On the impact of perturbative counterterms on black holes

Jesse Daas<EMAIL_ADDRESS>Cristobal Laporte
<EMAIL_ADDRESS>Frank Saueressig<EMAIL_ADDRESS>Institute for Mathematics, Astrophysics and Particle Physics (IMAPP)
Radboud University, Heyendaalseweg 135, 6525 AJ Nijmegen,The Netherlands

###### Abstract

We determine the corrections to the Schwarzschild geometry arising from
including the Goroff-Sagnotti counterterm in the gravitational dynamics. We
find that static, asymptotically flat, and spherically symmetric geometries
are completely characterized by their asymptotic mass and the coupling
associated with the counterterm. The latter induces distinct corrections at
sixth order of the parameterized post-Newtonian expansion. The resulting
spacetime geometries still exhibit an event horizon. In the parameter space
accessible to numerical integration, the horizon area is smaller than its
Schwarzschild counterpart, leading to an increase in the Hawking temperature.
Corrections to the shadow size can be determined analytically and are used to
give a first bound on the new coupling. While it is difficult to access the
geometry inside of the event horizon, our analysis also provides evidence that
the counterterm could resolve the curvature singularity appearing in the
Schwarzschild geometry.

##  I Introduction

The detection of gravitational waves from binary black hole mergers [1, 2] and
shadow imaging of supermassive black holes [3, 4, 5] have put a spotlight on
black hole research. Most of the physics insights gained in this way rely on
simulations which model the black hole spacetimes by the vacuum solutions of
classical general relativity, i.e., the Schwarzschild spacetime and its
rotating extension given by the Kerr solution. While this is perfectly
reasonable at the classical level, the tension between general relativity and
quantum mechanics is a strong pointer that this is not the final picture. On
general grounds any theory of quantum gravity is expected to modify the
classical black hole spacetimes in one way or another, e.g., when aiming at
resolving the spacetime singularities featured by the classical solutions [6]
or when dealing with the information paradox [7].

This letter takes an important step towards understanding black holes beyond
general relativity. Concretely, we determine the corrections to the
Schwarzschild geometry arising from the two-loop counterterm found by Goroff
and Sagnotti [8, 9, 10]. We show that the admissible spherically symmetric,
asymptotically flat and static geometries have “no hair” in the sense that
they are completely characterized by their asymptotic mass and the value of
the new coupling. The latter modifies the geometry at sixth order in the
parameterized post-Newtonian (PPN) expansion. Moreover, we show that the event
horizon persists and that the solution in the exterior can be accurately
approximated by a simple analytic form. We also argue for singularity
resolution, provided that the external solution connects to our local analysis
at the core. A more detailed analysis will be provided in [11].

##  II General Motivation and Setup

Constructing alternatives to the black hole spacetimes known from general
relativity has, by now, evolved into a rich subject [12, 13, 14, 15, 16, 17,
18, 19, 20]. From a phenomenological perspective, one may construct models
geared towards resolving deficits intrinsic to the black holes of general
relativity. This motivates, e.g., the study of regular black hole spacetimes
[21, 22, 23, 24] with the Hayward solution [25] being one of the most
prominent examples. Alternatively, one may study black holes based on
classical extensions of general relativity. Prototypical examples are the
black holes arising in scalar-tensor theories [26, 27], recently reviewed in
[28], where the scalar field is used to model a generic fifth force.

From the perspective of a quantum theory, one can regard the Einstein-Hilbert
action as the leading term in an effective field theory for the gravitational
interactions [29]. The gravitational dynamics then receives corrections in the
form of local higher-derivative terms and non-local contributions encoded in
form factors.111The effect of non-local corrections on the Schwarzschild
geometry has been considered in [30]. In the context of non-local, ghost-free
gravity it has been argued that non-locality resolves the curvature
singularities appearing in the Schwarzschild solution [31, 32, 33].
Restricting to the local part, one-loop corrections lead to new terms in the
dynamics which are quadratic in the spacetime curvature and one naturally
arrives at quadratic gravity [34, 35]. These contributions vanish when the
vacuum Einstein equations are imposed. From the perspective of quantizing
general relativity, this is the statement that the one-loop counterterms are
either topological or can be absorbed by a field redefinition [36]. For the
spacetime geometry, this entails that solutions of classical general
relativity are also solutions of quadratic gravity. While the solution space
of the latter is significantly richer [37, 38, 39, 40, 41, 42, 43, 44, 45, 46,
47, 48], any confirmation of vacuum properties of general relativity can also
be interpreted as a confirmation of quadratic gravity [47].

In order to arrive at a non-vanishing correction to the Schwarzschild
geometry, one needs to go beyond the terms quadratic in the spacetime
curvature. One attempt in this direction is Einstein-cubic gravity [49] where
the relative coefficients among the terms cubic in the spacetime curvature are
engineered in such a way that the resulting static and spherically symmetric
spacetime geometry is characterized by a single metric function of the radial
coordinate [50, 51]. This feature is non-generic though.

Taking inspiration from the perturbative quantization of general relativity,
one may then ask what is the “simplest” local modification, coming with the
lowest number of spacetime derivatives, which expels the Schwarzschild
geometry as a solution of the equations of motion. The unique answer to this
question is the Goroff-Sagnotti counterterm, coming with six spacetime
derivatives. This correction exhibits several outstanding properties. Firstly,
its contributions to the gravitational dynamics does not vanish when the
vacuum Einstein equations are imposed. Hence it gives rise to genuine
modifications of the black hole solutions known from general relativity.
Secondly, since one is dealing with a counterterm in the quantization of
general relativity, it is expected to appear in the effective gravitational
dynamics in a rather universal way. Thirdly, we verified that the counterterm
does not give rise to additional degrees of freedom, at least as a
perturbation around a flat background. Thus the Yukawa-type interactions
associated with the massive ghost-modes appearing in quadratic gravity are
absent in this case. These properties provide a strong motivation for studying
static, spherically symmetric spacetimes with this particular contribution
added to the dynamics.

Following the path motivated above, our starting point is the Einstein-Hilbert
(EH) action with the cosmological constant set to zero,

$S^{\rm EH}[g]=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}R\,,$ (1)

where $G$ is Newton’s constant and $g\equiv\det(g_{\mu\nu})$ is the
determinant of the spacetime metric $g_{\mu\nu}$. The covariant derivative,
Ricci scalar, Ricci tensor and Weyl tensor constructed from $g_{\mu\nu}$ are
denoted by $D_{\mu}$, $R$, $R_{\mu\nu}$ and $C_{\mu\nu\rho\sigma}$,
respectively. The dynamics (1) is complemented by the Goroff-Sagnotti (GS)
term,

$S^{\rm GS}[g]=\frac{G\,\lambda}{16\pi}\int
d^{4}x\sqrt{-g}\left(C_{\mu\nu}^{\phantom{\mu\nu}\kappa\gamma}\,C_{\kappa\gamma}^{\phantom{\kappa\gamma}\rho\sigma}\,C_{\rho\sigma}^{\phantom{\rho\sigma}\mu\nu}\right)\,,$
(2)

where $\lambda$ is an a priori undetermined, dimensionless coupling constant.
For later purposes, it is convenient to introduce
$\tilde{\lambda}\equiv\lambda G^{2}$ as a short-hand. The perturbative
quantization of (1) necessitates the inclusion of (2) in order to absorb the
infinities appearing in the quantization procedure at two loop-level [8, 9,
10].

The equations of motion are obtained from the variation principle

$H_{\mu\nu}\equiv\frac{1}{\sqrt{-g}}\,\frac{\delta}{\delta
g^{\mu\nu}}\left(S^{\rm EH}+S^{\rm GS}\right)=0\,.$ (3)

The explicit form of $H_{\mu\nu}$ is readily obtained using computer algebra
packages like the xAct package for Mathematica [52]. The result is lengthy and
little illuminating, so we refrain from giving it explicitly. The fact that
$H_{\mu\nu}$ is obtained from a diffeomorphism invariant action functional
entails that the components appearing in (3) are related by
$D^{\mu}H_{\mu\nu}=0$.

We are interested in static, spherically symmetric solutions of the system
(3). Using Schwarzschild coordinates, the most general form of the line-
element compatible with these symmetry requirements is [53]

$ds^{2}=-h(\bar{r})dt^{2}+f(\bar{r})^{-1}d\bar{r}^{2}+\bar{r}^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)\,.$ (4)

where the radial coordinate $\bar{r}$ is defined geometrically via the area of
the two-spheres coordinatized by $\theta,\phi$. For the Schwarzschild
geometry, the two metric functions $h(\bar{r})$ and $f(\bar{r})$ agree and are
given by

$h(\bar{r})=f(\bar{r})=1-\frac{2GM}{\bar{r}}\,.$ (5)

Here $M$ denotes the asymptotic mass of the geometry. Birkhoff’s theorem
guarantees that this is the only free parameter appearing in spherically
symmetric vacuum solutions of general relativity. Due to the contributions of
the Goroff-Sagnotti term in (3), the geometry (5) is not a solution of the
modified equations of motion though.

In order to study the solutions of the system (3) efficiently, we follow a
strategy highly successful in the context of quadratic gravity [40, 41, 42,
43, 44], and recast the line-element (4) into conformal-to-Kundt form

$ds^{2}=\Omega^{2}(r)\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}-2dudr+\mathcal{H}(r)du^{2}\right)\,.$ (6)

The coordinate transformation relating (4) and (6) is

$\bar{r}=\Omega(r)\,,\qquad t=u-\int\frac{dr}{\mathcal{H}(r)}\,,$ (7)

and the metric functions obey

$h=-\Omega^{2}(r)\,\mathcal{H}(r)\,,\quad
f=-\left(\frac{\Omega^{\prime}(r)}{\Omega(r)}\right)^{2}\,\mathcal{H}(r)\,.$
(8)

A direct consequence of eqs. (7) and (8) is that the center of the geometry,
$\bar{r}=0$, is mapped to $\Omega(r)=0$. Horizons, satisfying
$f(\bar{r}_{h})=h(\bar{r}_{h})=0$, are encoded in the roots
$\mathcal{H}(r_{h})=0$ with $\Omega(r_{h})=\bar{r}_{h}$ being finite.

The main advantage of the conformal-to-Kundt description is that it does not
contain any explicit $r$-dependence. As a consequence, the equations of motion
for the metric functions $\Omega(r)$ and $\mathcal{H}(r)$, resulting from
substituting the line-element (6) into (3), are an _autonomous_ system of
coupled, non-linear differential equations

$H_{rr}=0\,,\quad H_{ru}=0\,,\quad H_{\theta\theta}=0\,,$ (9)

which contain fourth derivatives of $\Omega$ and $\mathcal{H}$ with respect to
$r$. Selecting two independent equations from the set (9), one can determine
the resulting geometries based on the _necessary_ conditions

$\displaystyle\mathcal{H}^{\prime\prime\prime}$
$\displaystyle=\tfrac{\tilde{\lambda}\left(3\mathcal{H}^{\prime}\Omega^{\prime}+\Omega(\mathcal{H}^{\prime\prime}-1)\right)\left(2+\mathcal{H}^{\prime\prime}\right)^{2}-54\mathcal{H}\Omega^{3}(\Omega^{\prime})^{2}-18\mathcal{H}^{\prime}\Omega^{\prime}\Omega^{4}-18\Omega^{5}}{3\,\tilde{\lambda}\,\Omega\,\mathcal{H}^{\prime}\left(2+\mathcal{H}^{\prime\prime}\right)},$
$\displaystyle\Omega^{\prime\prime}$
$\displaystyle=\tfrac{\tilde{\lambda}\left(2+\mathcal{H}^{\prime\prime}\right)^{3}-18\Omega^{4}(2+\mathcal{H}^{\prime\prime})-108\Omega^{3}\Omega^{\prime}\mathcal{H}^{\prime}}{108\mathcal{H}\Omega^{3}}\,.$
(10b)

Here the primes denote derivatives with respect to $r$ and we solved for the
highest derivative terms. Once a solution of the system is found, either in
the form of a power series or numerical integration, one then checks that the
full system (9) does not lead to any further restrictions.

The system (10) shows that the degree of the differential equation is
determined by the Goroff-Sagnotti term. In particular (10) changes from a
third order to a first order equation when $\tilde{\lambda}$ is set to zero.
Thus the limit $\tilde{\lambda}\rightarrow 0$ is highly non-trivial and it is
not guaranteed that all solutions of the system with $\tilde{\lambda}\not=0$
smoothly connect to the Schwarzschild geometry when the counterterm is
switched off.

##  III Black hole spacetimes

While the conformal-to-Kundt prescription of the geometry leads to a
significant simplification of the equations of motion, an analytic solution is
currently still out of reach. A first step towards understanding the impact of
the correction term then consists in determining local solutions of the
system, which can be obtained as a series expansion at an arbitrary point
$r_{0}$ through the Frobenius method. Building on the conformal-to-Kundt
coordinates, the generic form of this expansion is222Our analysis covers
solutions which admit a series expansion at the expansion point only. Tracking
non-analytic terms is beyond the present work.

$\begin{split}\Omega(r)=&(r-r_{0})^{n}\sum_{i=0}a_{i}\,(r-r_{0})^{i}\,,\\\
\mathcal{H}(r)=&(r-r_{0})^{p}\sum_{i=0}b_{i}\,(r-r_{0})^{i}\,.\end{split}$
(11)

Here $(n,p)$ are the scaling exponents of the solution and we have
$a_{0}\not=0$ and $b_{0}\not=0$ by definition.

Substituting these expansions into the system (10) and matching terms with
equal powers of $(r-r_{0})^{i}$ leads to a hierarchy of equations which
determines the admissible values for $(n,p)$ as solutions of the incidental
polynomial. The free coefficients $a_{i},b_{i}$ parameterizing the solution
space typically appear at low index numbers and the hierarchy of equations
then expresses the higher-order expansion parameters in terms of these free
parameters.

###  III.1 Corrections in the asymptotically flat region

We start by studying the corrections in the asymptotically flat regions. For
this purpose, we use the Schwarzschild coordinates and expand the metric
functions in terms of inverse powers of $\bar{r}$:333In the case of quadratic
gravity, this expansion fails to capture the contributions of the massive
graviton excitations which are non-analytic in the asymptotically flat region
[54]. In the present framework, we find these modes are absent. Thus a
polynomial expansion is expected to work.

$h(\bar{r})=\sum_{n=0}h_{n}\,\bar{r}^{-n}\,,\quad
f(\bar{r})=\sum_{n=0}f_{n}\,\bar{r}^{-n}\,,$ (12)

with $h_{0}=f_{0}=1$. Substituting this expansion into the corresponding
equations of motion and comparing powers of $1/\bar{r}$ yields

$\begin{split}h(\bar{r})=&1-\frac{\bar{r}_{h}}{\bar{r}}+\frac{6\,\tilde{\lambda}\,\bar{r}_{h}^{2}}{\bar{r}^{6}}-\frac{4\,\tilde{\lambda}\,\bar{r}_{h}^{3}}{\bar{r}^{7}}+O\left(\frac{\tilde{\lambda}^{2}\,\bar{r}_{h}^{3}}{\bar{r}^{11}}\right)\,,\\\
f(\bar{r})=&1-\frac{\bar{r}_{h}}{\bar{r}}+\frac{18\,\tilde{\lambda}\,\bar{r}_{h}^{2}}{\bar{r}^{6}}-\frac{16\,\tilde{\lambda}\,\bar{r}_{h}^{3}}{\bar{r}^{7}}+O\left(\frac{\tilde{\lambda}^{2}\,\bar{r}_{h}^{3}}{\bar{r}^{11}}\right)\,,\end{split}$
(13)

where we abbreviated $\bar{r}_{h}\equiv 2GM$. For $\tilde{\lambda}=0$ this
expression reduces to the Schwarzschild geometry (5) and we verified that the
expansion using conformal-to-Kundt coordinates gives the same result.

Eq. (13) is remarkable for three reasons. Firstly, the solution is determined
completely in terms of the asymptotic mass $M$ and the new coupling
$\tilde{\lambda}$. Thus the solution exhibits a “no-hair”-type property in the
sense that it is fixed completely once the asymptotic mass and the coupling
$\tilde{\lambda}$ are given. Secondly, the leading corrections occur at
$O(\bar{r}^{-6})$ only. Converting the line-element to isotropic coordinates
shows that this corresponds to corrections at the _sixth order in the PPN
expansion_. The lower-order PPN parameters agree with general relativity.
Thirdly, the corrections break the degeneracy in the metric functions so that
for $\tilde{\lambda}\not=0$ one has $h(\bar{r})\not=f(\bar{r})$ with _a fixed
ratio between the leading correction terms_. This characteristic fingerprint
invites probing (13) with experiments that are sensitive to a difference in
the metric functions.

###  III.2 Admissible local solutions

The Frobenius analysis is readily extended beyond the asymptotically flat
region. The most efficient way to proceed is to analyze the system (10) based
on the expansion (11) and subsequently check for potential additional
constraints coming from the full system (9). The independent solution classes
identified in this way are summarized in Table 1.

$(n,p)^{\phantom{\infty}}$ | expansion point | free parameters | comments
---|---|---|---
$(0,0)^{\phantom{\infty}}$ . | regular point |  $a_{0},a_{1},b_{0},b_{1},b_{2};r_{0}$ |
$(-1,2)^{0}_{\phantom{\infty}}$ |  asymptotically flat region | $a_{0},b_{0},b_{1}=-2M$ |  $a_{0}=-1,b_{0}=-1$ are fixed by
$(r\rightarrow 0,r_{0}=0)$ |  asymptotic flatness and canonical time coordinate
$(-1,2)^{\infty}$ | center of the geometry |  | complex solution space owed to non-linear
$(r\rightarrow\infty,r_{0}=0)$ |  | nature of the system (10)
$(0,1)^{\phantom{\infty}}$ . | horizon | $a_{0},b_{0},b_{1};r_{0}=r_{h}$ |
$(0,2)^{\phantom{\infty}}$ . | double horizon | $a_{0};r_{0}=r_{h}$ | may not allow to fix the coordinate freedom completely
Table 1:  Classification of the admissible local scaling behaviors arising
from solving the system (10) based on the ansatz (13).

We highlight the most important properties of the results. The _asymptotically
flat region_ is covered by the solutions $(-1,2)^{0}$ coming with three free
parameters. The condition that the geometry approaches Minkowski space fixes
$a_{0}$ and $b_{0}$. This is the no-hair property highlighted in eq. (13). At
a _generic point_ , both $\Omega(r)$ and $\mathcal{H}(r)$ are finite by
definition. This is reflected in the $(0,0)$ class. Keeping $r_{0}$ fixed,
this class has five free parameters. Two of these may be fixed by a canonical
normalization of $t,r$ in the asymptotically flat region. The fact that one is
then left with more free parameters than in the $(-1,2)^{0}$ case is not a
contradiction per se: not all solutions found in the interior need to connect
to the asymptotically flat region. Our analysis also shows that there is a
consistent family of solutions $(0,1)$. This ensures that the equations of
motion actually admit geometries exhibiting an event horizon.

The most intriguing results are obtained in the core region described by the
$(-1,2)^{\infty}$ solution. In this case the analysis is very involved since
one has to track a significant number of branches due to the non-linear
appearance of the coefficients $a_{i},b_{i}$ in the hierarchy. Our analysis
suggests that there is one non-trivial, consistent solution. This solution has
the remarkable property that all sixteen Carminati–McLenaghan invariants [55],
including the Kretschmann scalar $R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$
and the square of the Weyl tensor, remain finite as $\bar{r}\rightarrow 0$.
Thus the class does not admit spacetime singularities. This is in stark
contrast to a similar analysis in quadratic gravity where the higher-
derivative terms lead to an amplification of the curvature singularity for
asymptotically flat solutions [43, 47].

###  III.3 Numerical integration

The Frobenius method gives local solutions at the core, the horizon, and the
asymptotically flat region. In order to arrive at the global picture, we then
integrate the system (10) numerically, placing initial conditions in the
asymptotically flat region. Our analysis focuses on the case $\lambda>0$.
While $\lambda<0$ can be plausible as well, the numerical integration of the
equations of motion is much more involved in this case and we have not been
able to connect the local solutions in this way.

Figure 1: Black hole solution obtained from integrating the system (10)
numerically with initial conditions $M=10$, $\lambda=0.1$ set at $r_{\rm
init}=-0.02$. The main diagram shows the metric functions $\mathcal{H}(r)<0$
(left-side scale) and $\Omega(r)>0$ (right-side scale) describing the black
hole in the conformal-to-Kundt coordinates. Its horizon is located at
$r_{h}^{\rm GS}=-0.0500000423$. The solution closely follows the Schwarzschild
solution (superimposed as black dashed lines) up to the horizon. Inset:
comparison between the metric function $f(\bar{r})$ of the numerical solution
(orange line) and the approximation (13), superimposed as a dashed-dotted
line. The difference between the exact solution and the analytic approximation
is found to be of the order $10^{-6}$. The metric function $h(\bar{r})$
exhibits the same features.

A prototypical solution is shown in Fig. 1. Starting from the asymptotically
flat region, the numerical integration of (10), reveals that the geometries
exhibit an event horizon where $\mathcal{H}(r)|_{r=r^{\rm GS}_{h}}=0$.
Moreover, the comparison between the numerical integration of (10) and the
analytic approximation (13) shows that the latter is reliable up to the
horizon, at least in the parametric regime where corrections of the order
$\tilde{\lambda}^{2}/M^{8}$ are negligible. Thus the asymptotically flat
solutions connect continuously to the local solutions at the event horizon and
we obtain a comprehensive picture of the geometry in the outside region of the
black hole.

##  IV Thermodynamics and shadow imaging

The series expansion (13) gives a good description of the geometry in the
outer region of the black hole. Thus we take it as a starting point for
determining the position of the event horizon and its properties. All analytic
results have been corroborated by integrating the exact equations numerically.

_Horizons and black hole thermodynamics._ Solving the horizon condition
$f(\bar{r}_{h}^{\rm GS};\lambda)=h(\bar{r}_{h}^{\rm GS};\lambda)=0$
perturbatively in $\lambda$ yields

$\bar{r}_{h}^{\rm
GS}(\tilde{\lambda})-\bar{r}_{h}=-\frac{\tilde{\lambda}}{4\,G^{3}M^{3}}+O\left(\tilde{\lambda}^{2}\right)\,.$
(14)

Thus, for $\tilde{\lambda}>0$, the corrections due to the Goroff-Sagnotti term
lead to black holes which have a smaller horizon area than their Schwarzschild
counterparts with the same asymptotic mass $M$.

The horizon also allows to associate a horizon temperature,

$T\equiv\frac{\kappa}{2\pi}\,.$ (15)

The surface gravity $\kappa$ is defined via
$\xi^{\mu}D_{\mu}\xi^{\nu}=\kappa\xi^{\nu}$ and $\xi^{\mu}=(1,0,0,0)$ is the
Killing vector associated with time translations. A brief computation based on
the expansion (13) yields

$T^{\rm GS}(\tilde{\lambda})-T^{\rm
SS}=\frac{\tilde{\lambda}}{16\pi\,G^{5}M^{5}}+O\left(\tilde{\lambda}^{2}\right)\,.$
(16)

Thus, $\tilde{\lambda}>0$ leads to an increase in the horizon temperature.
This is illustrated in Fig. 2 where the analytic result (16) is also confirmed
by numerical integration.

Figure 2:  Dependence of the horizon temperature (15) on $\lambda$ for an
asymptotic mass $M=10$ ($G=1$). The dashed horizontal line gives the
Schwarzschild result, the dashed-dotted line is the result (16) and we have
superimposed the result from integrating the equations of motion numerically
as the gray line. The spikes constitute an artifact of our numerical
integration procedure.

_Bounding $\lambda$ from shadow observations._ A powerful tool for probing the
spacetime geometry in the vicinity of a black hole horizon is taking an image
resolving the central brightness depleted region (shadow imaging). The size of
the shadow is determined by the photon ring located at coordinate radius
$\bar{r}_{\rm ph}$ where light rays can undergo an unstable circular orbit.
For the geometry (4) one has [56, 47]

$\bar{r}_{\rm ph}=2\,h(\bar{r}_{\rm
ph})\left(\left.\frac{dh(\bar{r})}{d\bar{r}}\right|_{\bar{r}=\bar{r}_{\rm
ph}}\right)^{-1}\,,\;\bar{r}_{\rm sh}=\frac{\bar{r}_{\rm
ph}}{\sqrt{h(\bar{r}_{\rm ph})}}\,.$ (17)

Here the first relation is an implicit equation for $\bar{r}_{\rm ph}$. The
second identity determines the shadow size $\bar{r}_{\rm sh}$, taking into
account that the bending of light close to the black hole causes the photon
ring to appear bigger than it really is. The angular radius of the photon ring
on the sky is then given by $\theta_{\rm sh}=\bar{r}_{\rm sh}/D$ with $D$ the
distance to the object. For the $M_{87}^{*}$ black hole, one of two Event
Horizon Telescope candidates, observations give $\theta_{\rm sh}=42\pm 3$
$\mu$as [3].

For the Schwarzschild case the relations (17) evaluate to

$\bar{r}_{\rm ph}^{\rm SS}=3M\,,\quad\bar{r}_{\rm sh}^{\rm
SS}=3\sqrt{3}\,M\,.$ (18)

We then determine the impact of $\lambda$ on the shadow size, based on eq.
(13). Defining $\zeta\equiv\frac{\bar{r}_{\rm sh}^{\rm
GS}(\tilde{\lambda})-\bar{r}_{\rm sh}^{\rm SS}}{\bar{r}_{\rm sh}^{\rm SS}},$
one obtains

$\zeta=-\frac{320}{729}\,\frac{\tilde{\lambda}}{\bar{r}_{h}^{4}}+\mathcal{O}\big{(}\frac{\tilde{\lambda}^{2}}{M^{8}}\big{)}\,.$
(19)

The sign indicates that for $\tilde{\lambda}>0$ the shadow size is decreased.
This is consistent with the black hole being more compact. The Schwarzschild
radius of $M^{*}_{87}$ is estimated as $\bar{r}_{h}\approx 2\times 10^{13}$ m.
Asking that the correction term changes the angular radius by at most 10$\%$
then bounds the dimensionless parameter $|\lambda|<5.4\times 10^{191}$.

##  V Summary and outlook

The need for theoretically well-motivated black hole spacetimes going beyond
the framework of classical general relativity can hardly be overstated. In
this letter, we contributed to this theme by determining the imprints of the
Goroff-Sagnotti counterterm on the Schwarzschild geometry. Eq. (13) shows that
this leads to characteristic correction terms which are excluded by the
uniqueness theorems of general relativity and modify the geometry at sixth
order in the PPN. Finding conclusive evidence for these extra terms would be a
major step towards testing fundamental ideas about quantum gravity, since any
quantum theory of the gravitational interactions has to address the role of
the perturbative singularity canceled by the Goroff-Sagnotti term at some
point. This includes, in particular, the gravitational asymptotic safety
program [57, 58, 59, 60, 61], where the value of the new coupling can be
predicted from first principle computations [62]. It would be very exciting to
explore to which extend gravitational wave signals emitted in a binary black
hole mergers are sensitive to our modifications.

###### Acknowledgements.

We thank A. Bonanno, A. Held, B. Knorr, J. Oliva, J. Podolsky, A. Ribes, S.
Silveravalle, M. Wondrak and J. Zanelli for insightful discussions. The work
of C.L. is supported by the scholarship Becas Chile ANID-PCHA/2020-72210073.
This publication is part of the Dutch Black Hole Consortium with project
number NWA.1292.19.202 of the research programme NWA which is (partly)
financed by the Dutch Research Council (NWO).

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2023] R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3, Phys. Rev. X 13, 011048 (2023), arXiv:2111.03634 [astro-ph.HE] .
  * Akiyama  _et al._ [2019a] K. Akiyama _et al._ (Event Horizon Telescope), First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole, Astrophys. J. Lett. 875, L1 (2019a), arXiv:1906.11238 [astro-ph.GA] .
  * Akiyama  _et al._ [2019b] K. Akiyama _et al._ (Event Horizon Telescope), First M87 Event Horizon Telescope Results. VI. The Shadow and Mass of the Central Black Hole, Astrophys. J. Lett. 875, L6 (2019b), arXiv:1906.11243 [astro-ph.GA] .
  * Akiyama  _et al._ [2022] K. Akiyama _et al._ (Event Horizon Telescope), First Sagittarius A* Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole in the Center of the Milky Way, Astrophys. J. Lett. 930, L12 (2022), arXiv:2311.08680 [astro-ph.HE] .
  * Bambi [2023] C. Bambi,  _Regular Black Holes_, Springer Series in Astrophysics and Cosmology (Springer Singapore, 2023) arXiv:2307.13249 [gr-qc] .
  * Chen  _et al._ [2015] P. Chen, Y. C. Ong, and D.-h. Yeom, Black Hole Remnants and the Information Loss Paradox, Phys. Rept. 603, 1 (2015), arXiv:1412.8366 [gr-qc] .
  * Goroff and Sagnotti [1985] M. H. Goroff and A. Sagnotti, Quantum Gravity at Two Loops, Phys. Lett. B 160, 81 (1985).
  * Goroff and Sagnotti [1986] M. H. Goroff and A. Sagnotti, The Ultraviolet Behavior of Einstein Gravity, Nucl. Phys. B 266, 709 (1986).
  * van de Ven [1992] A. E. M. van de Ven, Two loop quantum gravity, Nucl. Phys. B 378, 309 (1992).
  * [11] J. Daas, C. Laporte, B. Knorr, and F. Saueressig, in preparation .
  * Svarc  _et al._ [2020] R. Svarc, J. Podolsky, and O. Hruska, Kundt spacetimes in the Einstein-Gauss-Bonnet theory, Phys. Rev. D 102, 084012 (2020), arXiv:2007.06648 [gr-qc] .
  * Knorr and Platania [2022] B. Knorr and A. Platania, Sifting quantum black holes through the principle of least action, Phys. Rev. D 106, L021901 (2022), arXiv:2202.01216 [hep-th] .
  * Bueno  _et al._ [2023] P. Bueno, P. A. Cano, R. A. Hennigar, M. Lu, and J. Moreno, Generalized quasi-topological gravities: the whole shebang, Class. Quant. Grav. 40, 015004 (2023), arXiv:2203.05589 [hep-th] .
  * Bonanno  _et al._ [2022] A. Bonanno, S. Silveravalle, and A. Zuccotti, Nonsymmetric wormholes and localized big rip singularities in Einstein-Weyl gravity, Phys. Rev. D 105, 124059 (2022), arXiv:2204.04966 [gr-qc] .
  * Silveravalle and Zuccotti [2023] S. Silveravalle and A. Zuccotti, Phase diagram of Einstein-Weyl gravity, Phys. Rev. D 107, 064029 (2023), arXiv:2210.13877 [gr-qc] .
  * Held and Lim [2021] A. Held and H. Lim, Nonlinear dynamics of quadratic gravity in spherical symmetry, Phys. Rev. D 104, 084075 (2021), arXiv:2104.04010 [gr-qc] .
  * Aguilar-Gutierrez  _et al._ [2023] S. E. Aguilar-Gutierrez, P. Bueno, P. A. Cano, R. A. Hennigar, and Q. Llorens, Aspects of higher-curvature gravities with covariant derivatives, (2023), arXiv:2310.09333 [hep-th] .
  * Pawlowski and Tränkle [2023] J. M. Pawlowski and J. Tränkle, Effective action and black hole solutions in asymptotically safe quantum gravity, (2023), arXiv:2309.17043 [hep-th] .
  * Álvarez  _et al._ [2023] E. Álvarez, J. Anero, and E. Velasco-Aja, The fate of horizons under quantum corrections, JCAP 02, 029, arXiv:2207.12721 [hep-th] .
  * Bardeen [1968] J. Bardeen, Non-singular general-relativistic gravitational collapse, in proceedings of The International Conference GR5, Tbilisi, U.S.S.R. (1968) , 174 (1968).
  * Dymnikova [1992] I. Dymnikova, Vacuum nonsingular black hole, Gen. Rel. Grav. 24, 235 (1992).
  * Bonanno and Reuter [2000] A. Bonanno and M. Reuter, Renormalization group improved black hole space-times, Phys. Rev. D 62, 043008 (2000), arXiv:hep-th/0002196 .
  * Frolov [2016] V. P. Frolov, Notes on nonsingular models of black holes, Phys. Rev. D 94, 104056 (2016), arXiv:1609.01758 [gr-qc] .
  * Hayward [2006] S. A. Hayward, Formation and evaporation of regular black holes, Phys. Rev. Lett. 96, 031103 (2006), arXiv:gr-qc/0506126 .
  * Sotiriou and Faraoni [2012] T. P. Sotiriou and V. Faraoni, Black holes in scalar-tensor gravity, Phys. Rev. Lett. 108, 081103 (2012), arXiv:1109.6324 [gr-qc] .
  * Hui and Nicolis [2013] L. Hui and A. Nicolis, No-Hair Theorem for the Galileon, Phys. Rev. Lett. 110, 241104 (2013), arXiv:1202.1296 [hep-th] .
  * Herdeiro and Radu [2015] C. A. R. Herdeiro and E. Radu, Asymptotically flat black holes with scalar hair: a review, Int. J. Mod. Phys. D 24, 1542014 (2015), arXiv:1504.08209 [gr-qc] .
  * Donoghue  _et al._ [2017] J. F. Donoghue, M. M. Ivanov, and A. Shkerin, EPFL Lectures on General Relativity as a Quantum Field Theory, (2017), arXiv:1702.00319 [hep-th] .
  * Calmet  _et al._ [2022] X. Calmet, R. Casadio, S. D. H. Hsu, and F. Kuipers, Quantum Hair from Gravity, Phys. Rev. Lett. 128, 111301 (2022), arXiv:2110.09386 [hep-th] .
  * Koshelev  _et al._ [2018] A. S. Koshelev, J. a. Marto, and A. Mazumdar, Schwarzschild $1/r$-singularity is not permissible in ghost free quadratic curvature infinite derivative gravity, Phys. Rev. D 98, 064023 (2018), arXiv:1803.00309 [gr-qc] .
  * Buoninfante  _et al._ [2018] L. Buoninfante, A. S. Koshelev, G. Lambiase, J. a. Marto, and A. Mazumdar, Conformally-flat, non-singular static metric in infinite derivative gravity, JCAP 06, 014, arXiv:1804.08195 [gr-qc] .
  * Boos [2021] J. Boos, Non-singular ”Gauss” black hole from non-locality: a simple model with a de Sitter core, mass gap, and no inner horizon, (2021), arXiv:2104.00555 [gr-qc] .
  * Stelle [1977] K. S. Stelle, Renormalization of Higher Derivative Quantum Gravity, Phys. Rev. D 16, 953 (1977).
  * Stelle [1978] K. S. Stelle, Classical Gravity with Higher Derivatives, Gen. Rel. Grav. 9, 353 (1978).
  * ’t Hooft and Veltman [1974] G. ’t Hooft and M. J. G. Veltman, One loop divergencies in the theory of gravitation, Ann. Inst. H. Poincare Phys. Theor. A 20, 69 (1974).
  * Lü  _et al._ [2015] H. Lü, A. Perkins, C. N. Pope, and K. S. Stelle, Black holes in D = 4 higher-derivative gravity, Int. J. Mod. Phys. A 30, 1545016 (2015).
  * Lu  _et al._ [2015] H. Lu, A. Perkins, C. N. Pope, and K. S. Stelle, Black Holes in Higher-Derivative Gravity, Phys. Rev. Lett. 114, 171601 (2015), arXiv:1502.01028 [hep-th] .
  * Lü  _et al._ [2015] H. Lü, A. Perkins, C. N. Pope, and K. S. Stelle, Spherically Symmetric Solutions in Higher-Derivative Gravity, Phys. Rev. D 92, 124019 (2015), arXiv:1508.00010 [hep-th] .
  * Pravda  _et al._ [2017] V. Pravda, A. Pravdova, J. Podolsky, and R. Svarc, Exact solutions to quadratic gravity, Phys. Rev. D 95, 084025 (2017), arXiv:1606.02646 [gr-qc] .
  * Svarc  _et al._ [2018] R. Svarc, J. Podolsky, V. Pravda, and A. Pravdova, Exact black holes in quadratic gravity with any cosmological constant, Phys. Rev. Lett. 121, 231104 (2018), arXiv:1806.09516 [gr-qc] .
  * Podolsky  _et al._ [2018] J. Podolsky, R. Svarc, V. Pravda, and A. Pravdova, Explicit black hole solutions in higher-derivative gravity, Phys. Rev. D 98, 021502 (2018), arXiv:1806.08209 [gr-qc] .
  * Podolský  _et al._ [2020] J. Podolský, R. Švarc, V. Pravda, and A. Pravdova, Black holes and other exact spherical solutions in Quadratic Gravity, Phys. Rev. D 101, 024027 (2020), arXiv:1907.00046 [gr-qc] .
  * Pravda  _et al._ [2021] V. Pravda, A. Pravdova, J. Podolsky, and R. Svarc, Black holes and other spherical solutions in quadratic gravity with a cosmological constant, Phys. Rev. D 103, 064049 (2021), arXiv:2012.08551 [gr-qc] .
  * Bonanno and Silveravalle [2019] A. Bonanno and S. Silveravalle, Characterizing black hole metrics in quadratic gravity, Phys. Rev. D 99, 101501 (2019), arXiv:1903.08759 [gr-qc] .
  * Bonanno and Silveravalle [2021] A. Bonanno and S. Silveravalle, The gravitational field of a star in quadratic gravity, JCAP 08, 050, arXiv:2106.00558 [gr-qc] .
  * Daas  _et al._ [2023] J. Daas, K. Kuijpers, F. Saueressig, M. F. Wondrak, and H. Falcke, Probing quadratic gravity with the Event Horizon Telescope, Astron. Astrophys. 673, A53 (2023), arXiv:2204.08480 [gr-qc] .
  * Silveravalle [2022] S. Silveravalle, The gravitational field of isolated objects in quadratic gravity, Nuovo Cim. C 45, 153 (2022), arXiv:2202.00999 [gr-qc] .
  * Bueno and Cano [2016a] P. Bueno and P. A. Cano, Einsteinian cubic gravity, Phys. Rev. D 94, 104005 (2016a), arXiv:1607.06463 [hep-th] .
  * Hennigar and Mann [2017] R. A. Hennigar and R. B. Mann, Black holes in Einsteinian cubic gravity, Phys. Rev. D 95, 064055 (2017), arXiv:1610.06675 [hep-th] .
  * Bueno and Cano [2016b] P. Bueno and P. A. Cano, Four-dimensional black holes in Einsteinian cubic gravity, Phys. Rev. D 94, 124051 (2016b), arXiv:1610.08019 [hep-th] .
  * Martin-Garcia and et. al. [2013] J. M. Martin-Garcia, et. al., http://xact.es/ (2002-2013).
  * Hawking and Ellis [2023] S. W. Hawking and G. F. R. Ellis,  _The Large Scale Structure of Space-Time_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2023).
  * Saueressig  _et al._ [2021] F. Saueressig, M. Galis, J. Daas, and A. Khosravi, Asymptotically flat black hole solutions in quadratic gravity, Int. J. Mod. Phys. D 30, 2142015 (2021).
  * Carminati and McLenaghan [1991] J. Carminati and R. G. McLenaghan, Algebraic invariants of the Riemann tensor in a four-dimensional Lorentzian space, J. Math. Phys. 32, 3135 (1991).
  * Psaltis  _et al._ [2020] D. Psaltis _et al._ (Event Horizon Telescope), Gravitational Test Beyond the First Post-Newtonian Order with the Shadow of the M87 Black Hole, Phys. Rev. Lett. 125, 141104 (2020), arXiv:2010.01055 [gr-qc] .
  * Percacci [2017] R. Percacci,  _An Introduction to Covariant Quantum Gravity and Asymptotic Safety_, 100 Years of General Relativity, Vol. 3 (World Scientific, 2017).
  * Reuter and Saueressig [2019] M. Reuter and F. Saueressig,  _Quantum Gravity and the Functional Renormalization Group: The Road towards Asymptotic Safety_ (Cambridge University Press, 2019).
  * Saueressig [2023] F. Saueressig, The Functional Renormalization Group in Quantum Gravity, (2023), arXiv:2302.14152 [hep-th] .
  * Eichhorn and Schiffer [2022] A. Eichhorn and M. Schiffer, Asymptotic safety of gravity with matter, (2022), arXiv:2212.07456 [hep-th] .
  * Pawlowski and Reichert [2023] J. M. Pawlowski and M. Reichert, Quantum Gravity from dynamical metric fluctuations, (2023), arXiv:2309.10785 [hep-th] .
  * Gies  _et al._ [2016] H. Gies, B. Knorr, S. Lippoldt, and F. Saueressig, Gravitational Two-Loop Counterterm Is Asymptotically Safe, Phys. Rev. Lett. 116, 211302 (2016), arXiv:1601.01800 [hep-th] .


# Dispersive interaction between two atoms in Proca Quantum Electrodynamics

Gabriel Camacho de Pinho, Carlos Augusto Domingues Zarro, Carlos Farina
Instituto de Fisica, Universidade Federal do Rio de Janeiro, Av. Athos da
Silveira Ramos, 149, Centro de Tecnologia - bloco A, Cidade Universitária, Rio
de Janeiro - RJ - CEP: 21941-909, Brasil  Reinaldo de Melo e Souza  Instituto
de Física, Universidade Federal Fluminense, Av. Litorânea, s/n, Niterói, RJ,
24210-346, Brasil  Maurício Hippert  Illinois Center for Advanced Studies of
the Universe & Department of Physics, University of Illinois Urbana-Champaign,
Urbana, IL 61801, USA  Instituto de Física, Universidade do Estado do Rio de
Janeiro, Rua São Francisco Xavier, 524, Rio de Janeiro, RJ, 20550-013, Brasil

###### Abstract

We analyze the influence of a massive photon in the dispersive interaction
between two atoms in their fundamental states. We work in the context of Proca
Quantum Electrodynamics. The photon mass not only introduces a new length
scale but also gives rise to a longitudinal polarization for the
electromagnetic field. We obtain explicitly the interaction energy between the
atoms for any distance regime and consider several particular cases. We show
that, for a given interatomic distance, the greater the photon mass the better
it is the non-retarded approximation.

##  I Introduction

Neutral but polarizable atoms with no permanent multipoles placed in vacuum
will still interact with each other. This can be understood from the
Heisenberg uncertainty principle which enforces dipole fluctuation even in
states where the expectation value of the dipole operator vanishes. These
forces are the so-called dispersive interactions [1, 2, 3, 4, 5]. They are
responsible for several phenomena not only in Physics, but also in Chemistry,
Biology and many other areas of science. For instance, they explain why noble
gases condensate. Besides, since dispersive forces strongly depend on the
polarizabilities of the atoms/molecules involved, which in turn scale with the
volumes of such atoms/molecules, these forces explain why the temperatures of
condensation of noble gases are greater for the larger atoms (with Helium
having the smallest and so on). By the same argument, dispersive forces may
become dominant in macromolecules interactions. Increase attention are being
given to these forces since they are responsible for the cohesion between
different layers of the van der Waals heterostructures, materials obtained by
stacking different 2d materials and which display fascinating properties[6].
There is still an intense research going also from the quantum field theory
perspective[7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. In Biology, they are
responsible for the softness of our skin and explain the remarkable adhesion
of geckos on the walls of our houses[19]. We can go even further with more
bizarre examples, like the important role of dispersive forces in the
generation of the high electrostatic potentials ocurring during the storms,
among others [20, 21]. Differently from other kinds of intermolecular forces,
dispersion forces are always present due to the ubiquity of Heisenberg
uncertainty principle.

We can think of dispersive forces as a consequence of the exchange of virtual
photons between the atoms. In Maxwell electrodynamics in vacuum, photons are
massless, and this is responsible for some important features of these forces,
such as the long range character of the interaction. However, under certain
circumstances, photons may acquire an effective mass due to shielding
mechanism. One example arises in colloidal systems. When two atoms interact
within an electrolyte solution, the potential generated by the atoms ionize
the solution and when the potential is not very strong - that is if the atoms
are not very close to each other - we may treat the interaction within the
linearized Poisson-Boltzmann equation[22]. In this scenario the Maxwell field
in the solution behaves as a free Proca effective electromagnetic field, with
the photon mass being given by the inverse of the Debye wavelength. Similar
mechanisms arise in metals, or inside waveguides. Indeed, in many situations,
a gas of photons between conducting plates may be thought of as a massive bi-
dimensional photon gas [23]. Such an effective photon mass affects not only
real, but also virtual photons, and thus will have consequences for dispersive
interactions. In fact, a recent work [24] has also shown that the dispersive
interaction between two atoms inside a plane capacitor is screened due to non-
additive effects, in a way that can be effectively described by endowing
photons with a mass inversely proportional to the separation between capacitor
plates. Also, electric as well as magnetic condensates may lead to effectively
massive gauge fields and can be used to investigate superconductivity,
confinement in QCD, and an holographic duality scenario to explain the metal-
insulator transition in condensed matter [25, 26, 27, 28].

Even in vacuum, a tiny, yet non-vanishing, photon mass is also present in
scenarios of physics beyond the standard model [29, 30]. Moreover, these
Yukawa-type corrections appear in axion physics and modifications of Newtonian
gravity. Casimir effect experiments can be used to impose bounds on these
Yukawa parameters[31, 32, 33]. However small it may be, a finite mass for the
photon changes a lot our physical picture of the world. Indeed, as pointed out
in [34], the above bounds may not be correct depending upon the microscopical
origin of the mass. For example, if it appears from the Higgs mechanism, it is
possible for large-scale fields to be effectively Maxwellian, and in this case
astrophysical observations will be insensitive to the mass of the photon.

Many different observations place very stringent bounds for the photon mass.
The most rigorous restrictions come from astrophysical observations, which put
an upper bound of $10^{-27}$ eV to the mass of the photon [35]. Particle
properties, such as the anomalous magnetic moment of the electron, also place
strict limits on its mass [36]. Interestingly, some mass bound for photons can
be obtained using Schuman ressonances [37]. In Ref. [38] another terrestrial
experiment was proposed, in which the longitudinal mode of photons would be
measured to restrict the mass of the photon. Two comprehensive reviews,
written by the same authors almost 40 years apart, can be found in [39, 40].

Ultimately, dispersive interactions akin to van der Waals forces are to be
expected in any theory with Abelian vector bosons, as is the case in some
models of nuclear physics. In that context, vector mesons have a considerable
mass in vacuum, which are not protected by the global symmetries linked to
these particles. In particular, short-range repulsion between nucleons due to
the exchange of massive omega vector mesons plays an important role in the
structure of both nuclei and neutron stars [41, 42, 43]. In fact, the
interplay between attractive interactions mediated by scalar mesons and
repulsive interactions, mediated by such vector mesons, are thought to be the
mechanism behind nuclear saturation [44, 45]. The omega meson couples to
baryon charge, while mesons in general — composed of a valence quark-antiquark
pair — should display fluctuations of the dipole moment of baryon number
density. Therefore, omega meson exchange should lead to dispersive forces
between mesons. Recently, the coupling of omega mesons to pions and sigmas has
been proposed as a mechanism to explain the speed of sound peak expected in
ultra-dense nuclear matter [46]. Dispersive van der Waals forces within the
context of Quantum Chromodynamics have been previously explored in the
literature [47, 48].

The simplest description of massive vector mesons is given by Proca
electrodynamics, in which gauge invariance is explicitly lost [49]. While, in
some applications, a more fundamental description may be provided by other
formulations, which might preserve gauge symmetry, such as the Higgs mechanism
or the Stueckelberg action, Proca electrodynamics can be motivated as a
simplified version of these formulations, obtained, for instance, by treating
extra fields (e.g., the Higgs field) as static, uniform backgrounds. The
investigation of dispersive forces mediated by massive photons also raises
interesting conceptual questions. For instance, what is the nature of the
limit in which the photon mass is taken to zero? The disappearance of the
longitudinal polarization could raise the possibility of a discontinuous
transition to Maxwell theory, whereas a continuous transition seems to always
be found. Yet, the precise dependence of van der Waal forces on the mediator
mass near this limit will be revealing of how sensitive these forces are to
small photon mass in scenarios beyond the standard model. In the context of
the Casimir effect this discussion has revealed surprising subtleties. It was
initially predicted that the zero mass limit would be singular due to the
longitudinal modes[50], which was subsequently disproved[51, 52] by showing
that in the zero mass limit these modes decouple from matter and the plate
becomes transparent to them. In this way only transversal modes contribute to
the Casimir attraction between plates assuring that this limit is well
behaved. Furthermore, understanding dispersive interactions mediated by
massive photons may also furnish new insights into Maxwell electrodynamics.
For example, a photon mass breaks the degeneracy between the group, phase and
signal propagation velocities, present in standard electrodynamics. This will
enable us to verify which one plays the dominant role in the physics of
dispersive forces.

In this work we employ Proca electrodynamics to investigate the effects of a
finite photon mass on the dispersive interaction between two atoms. This
subject has recently got renewed interest in the search for signatures of
physics beyond standard models[53]. In this paper we employ a Hamiltonian
originally proposed by P. Milonni[1] and show that it is possible to analyze
the Proca dispersive interaction within first order perturbation theory.
Afterwards we analyze the short distance limit, known as the non-retarded
regime. We study this limit by two complementary approaches. In section II we
evaluate the interaction energy without quantizing the electromagnetic field,
while in section III we re-obtain the results of section II by taking the
appropriate limit of the full quantum electrodynamics treatment. By doing so,
we are able to address the question of when it is necessary to quantize the
electromagnetic showing that it leads to a weaker condition for the
interatomic distance. In other words, a photonic mass extends the range of
validity of the non-retarded regime. We leave section IV for our final
remarks.

##  II London-Proca interaction

Let us consider two atoms, $A$ and $B$, held at positions $\mathbf{r_{A}}$ and
$\mathbf{r_{B}}$, respectively and for convenience let us define
$\bm{R}=\bm{r}_{B}-\bm{r}_{A}$. Throughout this paper we shall assume the
magnitude of the distance $R$ to be much larger than the typical size of the
atoms, enabling us to employ the electric dipole approximation. Hence,
dispersive forces arise from correlations between the fluctuating dipoles in
each atom. In the London regime, we consider the field generated by the atomic
dipoles as electrostatic. This is a good approximation when the atoms are
close to each other, and in the next section we employ a full quantum
electrodynamics calculation in order to establish the precise condition $R$
must satify. In the dipole approximation, the interaction Hamiltonian
describing the coupling between the atoms is given by

$H_{int}=-\mathbf{d}_{A}\cdot\mathbf{E}_{B}(\mathbf{r_{A}})\,,$ (1)

where $\mathbf{d}_{A}$ denotes the electric dipole operator of atom $A$ and
$\mathbf{E}_{B}(\mathbf{r_{A}})$ stands for the electrostatic field created by
the dipole $B$ at the position of atom $A$. The electrostatic field
$\mathbf{E}_{B}$ depends upon the dipole operator of atom $B$,
$\mathbf{d}_{B}$, and is given by (see appendix A)

$\displaystyle\mathbf{E}_{B}(\mathbf{r}_{A})$ $\displaystyle=$
$\displaystyle-\frac{e^{-\mu R}}{R^{3}}(\mu R+1)\mathbf{d}_{B}+$ (2)
$\displaystyle+$ $\displaystyle\frac{e^{-\mu R}}{R^{5}}(\mu^{2}R^{2}+3\mu
R+3)(\mathbf{d}_{B}\cdot\mathbf{R}))\mathbf{R}\,,$ (4)

where $\mu\equiv mc\,(\hbar)^{-1}>0$, with $m$ denoting the photon mass in
Proca electrodynamics. Substituting the last equation into the Hamiltonian (1)
we obtain

$\displaystyle H_{int}=\frac{e^{-\mu R}}{R^{3}}(\mu
R+1)\mathbf{d}_{A}\cdot\mathbf{d}_{B}+$ (5) $\displaystyle-$
$\displaystyle\frac{e^{-\mu R}}{R^{5}}(\mu^{2}R^{2}+3\mu
R+3)(\mathbf{d}_{A}\cdot\mathbf{R})(\mathbf{d}_{B}\cdot\mathbf{R})\,.$ (7)

Notice that this Hamiltonian is symmetric under the exchange $A\leftrightarrow
B$ as it should - the interaction does not depend on which atom we choose as
the source of the electric field in Eq. (1). Denoting the ground state of atom
$j$ by $|0_{j}\rangle$, the fact that the atom does not possess permanent
electric dipole means that $\langle 0_{j}|\mathbf{d}_{j}|0_{j}\rangle=0$.
Therefore, perturbation theory applied to the above Hamiltonian vanishes at
first order. At second order we have

$U_{\rm NR}=-\sum_{r\neq 0,s\neq 0}\frac{\langle
0_{A}0_{B}|H_{int}|rs\rangle\langle
rs|H_{int}|0_{A}0_{B}\rangle}{\hbar(\omega_{r0}+\omega_{s0})}\,$ (8)

where we denote the excited states of atoms $A$ and $B$ by $|r\rangle$ and
$|s\rangle$, respectively. We assume that the atomic ground state is non-
degenerate, as usual. Choosing the $Z$-axis parallel to $\mathbf{R}$ we may
rewrite Eq. (7) to obtain

$\displaystyle H_{int}=\frac{e^{-\mu R}}{R^{3}}\Bigg{[}(\mu
R+1)(d^{x}_{A}d^{x}_{B}+d^{y}_{A}d^{y}_{B})$ (9) $\displaystyle-$
$\displaystyle(\mu^{2}R^{2}+2\mu R+2)d^{z}_{A}d^{z}_{B}\Bigg{]}\,.$ (11)

Substituting Eq. (11) into Eq. (8) we obtain

$U_{\rm NR}=-\Lambda\frac{e^{-2\mu R}}{9R^{6}}p(\mu R)\,,$ (12)

where

$p(\xi)=6+12\xi+10\xi^{2}+4\xi^{3}+\xi^{4}$ (13)

and $\Lambda>0$ is a constant depending only on the atomic internal structure,
given by

$\Lambda=\sum_{r\neq 0,s\neq
0}\frac{|\mathbf{d}_{A}^{0r}|^{2}|\mathbf{d}_{B}^{0r}|^{2}}{\hbar(\omega_{r0}+\omega_{s0})}\,.$
(14)

In the last expression we assume an isotropic response of atoms, which enables
us to make the replacement $\langle 0|d_{m}|I\rangle\langle
I|d_{n}|0\rangle=|\mathbf{d}^{0I}|^{2}\delta_{mn}/3$ for each atom. We
emphasize that, in principle, $\Lambda$ depends on $\mu$, since the photon
mass also affects the fields inside the atoms. Nonetheless, assuming $\mu a\ll
1$, where $a$ is the atom size, we can consider the atom structure to be
independent of $\mu$, and suppose that $\mu$ affects only the intermolecular
interaction. We shall assume this to be the case henceforth.

As expected, if $\mu=0$ we reobtain London’s formula, hereafter denoted by
$U_{\rm London}$. For $\mu R\ll 1$, the leading corrections to the London
potential are given by

$U=U_{\rm London}\left(1-\frac{1}{3}(\mu
R)^{2}+\mathcal{O}(\mu^{4}R^{4})\right)\,,$ (15)

and hence we see that a small $\mu$ weakens the interaction, as could be
expected, since the field produced by each fluctuating dipole becomes
exponentially damped with distance. Notice that the linear term is not
present, which could be anticipated by Eq.(4), since the Proca correction for
the dipole electric field is already of order $\mu^{2}$, as the reader may
verify. In the regime $\mu R\gg 1$, on the other hand, the interaction between
atoms is exponentially suppressed. In Fig.1 we plot the interaction energy
normalized by the usual London interaction as a function of the dimensionless
parameter $\mu R$.

Figure 1: Interaction energy in the non-retarded regime divided by the London-
van der Waals potential as a function of the dimensionless radius parameter
$\mu R$.

Equation (12) reproduces exactly the interaction energy of two atoms
interacting inside a colloidal system[22], justifying the use of Proca
electrodyamics as an effective theory for these systems. Within this effective
scenario the intraatomic field is Maxwellian and thus $\Lambda$ is not
dependent on $\mu$. This means that the ratio $U_{\rm NR}/U_{\rm London}$
depends only on the colloidal properties which defines $\mu$ and not on the
atoms partaking in the interaction.

Before closing this section, we follow [54] (see section 7.6) and rewrite
$\Lambda$ in the form

$\Lambda=\frac{9}{4\pi}\int_{-\infty}^{\infty}\alpha_{A}(i\chi)\alpha_{B}(i\chi)d\chi\,,$
(16)

where $\alpha_{j}$ denotes the dynamical polarizability of atom $j$, defined
by

$\alpha_{j}(\omega)=\frac{2}{3}\sum_{I\neq
0}\frac{\omega_{I0}|\mathbf{d}_{j}^{0I}|^{2}}{\hbar(\omega_{I0}^{2}-\omega^{2})}\,.$
(17)

Since the intermolecular interaction is due to the correlation between induced
dipoles, the polarizability is a key concept to the physical understanding of
the intermolecular interaction. Equation (17) illustrates the dispersive
character of the interaction, expressed by the dependence of the atomic
response $\alpha_{j}$ on the frequency of the external perturbation. As we
shall see in the next section, the polarizability is the most convenient
starting point for analyzing the interatomic interaction within quantum
electrodynamics.

##  III Quantum field theory approach

###  III.1 The complete potential

When retardation effects in the electromagnetic interaction are relevant, we
must include into the description of our physical system the mediator of the
interaction between the atoms, namely, the electromagnetic field. And since
dispersive forces have a quantum nature, we must quantize this field [1]. In
this section we obtain the dispersive interaction between two atoms in the
framework of Proca QED for any distance regime. For simplicity, we consider
the atoms inside a perfectly conducting cube of volume $V$, to be taken to
infinity after performing the calculation. The electric field operator in
Proca QED is given by[55]

$\displaystyle\mathbf{E}_{0}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{k}\lambda}[\mathbf{E}^{(+)}_{0,\mathbf{k}\lambda}+\mathbf{E}^{(-)}_{0,\mathbf{k}\lambda}]$
(18)

with

$\mathbf{E}^{(+)}_{0,\mathbf{k}\lambda}=i\sum_{\mathbf{k}\lambda}\left(\frac{2\pi\hbar\omega_{k}}{V}\right)^{1/2}f_{\lambda}\bm{{\epsilon}}_{\mathbf{k}\lambda}a_{\mathbf{k}\lambda}e^{i(\mathbf{r}\cdot\mathbf{k}-\omega_{k}t)}=\mathbf{E}^{(-)\dagger}_{0,\mathbf{k}\lambda}\,,$
(19)

where $a_{\mathbf{k}\lambda}$ and $a^{\dagger}_{\mathbf{k}\lambda}$ are the
usual annihilation and creation operators for the mode with wavevector
$\bm{k}$ and polarization $\lambda$, $\omega_{k}=c\sqrt{k^{2}+\mu^{2}}$,
$\bm{\epsilon}_{\mathbf{k}\lambda}$ are the three unitary polarization
vectors, with $\lambda=1,2$ denoting transversal polarizations, while
$\lambda=3$ corresponds to the longitudinal one, not present in Maxwell
electrodynamics. Choosing $\bm{\epsilon}_{\mathbf{k}\lambda}$ to be normalized
requires the introduction of the factor

$\displaystyle f_{\lambda}=\left\\{\begin{matrix}1&&\lambda=1,2\\\
\frac{c\mu}{\omega_{k}}&&\lambda=3\,,\end{matrix}\right.$ (20)

The interaction between the atoms can be described by the interaction
Hamiltonian introduced by P. Milonni [1](see specially his Eq. (3.73)), namely

$H_{int}=-\frac{1}{2}\sum_{\mathbf{k}\lambda}\alpha_{A}(\omega_{k})[\mathbf{E}_{0,\mathbf{k}\lambda}(\mathbf{r}_{A})+\mathbf{E}_{B,\mathbf{k}\lambda}(\mathbf{r}_{A})]^{2}\,,$
(21)

where $\mathbf{E}_{0,\mathbf{k}\lambda}$ is given in Eq. (18) and
$\mathbf{E}_{B,\mathbf{k}\lambda}$ is the electromagnetic field created by the
dipole induced in atom $B$ by the vacuum field, given by (see appendix A)

$\displaystyle\mathbf{E}_{B}^{(\pm)}(\mathbf{r}_{A},\omega)$ $\displaystyle=$
$\displaystyle\frac{e^{\mp ikR}}{R^{5}}(-k^{2}R^{2}\pm
3ikR+3)(\mathbf{d}_{B}^{(\pm)}(\omega_{k})\cdot\mathbf{R})\mathbf{R}\,$ (22)
$\displaystyle+$ $\displaystyle\frac{e^{\mp
ikR}}{R^{3}}(R^{2}\omega_{k}^{2}\mp ikR-1)\mathbf{d}_{B}^{(\pm)}(\omega_{k}),$
(24)

where the superscript $(\pm)$ follows the same convention as decomposition
(18). The dipole operator $\mathbf{d}_{B}^{(\pm)}(\omega_{k})$ inhabits the
field Hilbert space and is given by

$\mathbf{d}_{B}^{(\pm)}(\omega_{k})=\alpha_{B}(\omega_{k})\mathbf{E}^{(\pm)}_{0,\mathbf{k}\lambda}(\mathbf{r}_{B}).$
(25)

Physically, Hamiltonian (21) can be interpreted in the following way: it is
the energy of a dipole induced in the atom $A$ due to the field acting in it.
This field is the superposition of vacuum electromagnetic field with the
electric dipole field generated by atom $B$. This latter dipole is also not
permanent but instead induced by the vacuum field according to Eq.(25). The
interaction energy comes from a first order perturbation calculation, this
being the convenience of employing this Hamiltonian. The dominant order for
the interaction comes from the terms involving $\alpha_{A}\alpha_{B}$ which
are given by

$\displaystyle H_{int}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{\mathbf{k}\lambda}\alpha_{A}(\omega_{k})[\mathbf{E}_{0,\mathbf{k}\lambda}(\mathbf{r}_{A})\cdot\mathbf{E}_{B,\mathbf{k}\lambda}(\mathbf{r}_{A})+$
(26) $\displaystyle+$
$\displaystyle\mathbf{E}_{B,\mathbf{k}\lambda}(\mathbf{r}_{A})\cdot\mathbf{E}_{0,\mathbf{k}\lambda}(\mathbf{r}_{A})]\,.$
(28)

Its expected value yields the interaction energy

$\displaystyle
U_{\text{total}}\\!\\!=\\!\\!-\sum_{\mathbf{k}\lambda}\alpha_{A}(\omega_{k})\alpha_{B}(\omega_{k})g_{ij}(\mathbf{R})$
(29) (30)
$\displaystyle\times\bra{0}\\{E^{i}_{0,\mathbf{k}\lambda}(\mathbf{r}_{A}),E^{j}_{0,\mathbf{k}\lambda}(\mathbf{r}_{B})\\}\ket{0},$
(31)

where

$\displaystyle g_{ij}(\mathbf{R})$ $\displaystyle=$
$\displaystyle\text{Re}\,\Bigg{\\{}\frac{e^{-ikR}}{R^{3}}(R^{2}\omega^{2}_{k}-ikR-1)\delta_{ij}+$
(32) $\displaystyle+$
$\displaystyle\frac{e^{-ikR}}{R^{5}}(-k^{2}R^{2}+3ikR+3)R_{i}R_{j}\Bigg{\\}}\,.$
(34)

From Eqs. (18) and (19) we obtain

$\displaystyle U_{\text{total}}$ $\displaystyle=$
$\displaystyle-\text{Re}\sum_{\mathbf{k}\lambda}\left(\frac{2\pi\hbar\omega_{k}}{V}\right)\alpha_{A}(\omega_{k})\alpha_{B}(\omega_{k})$
(37) $\displaystyle
g_{ij}(\mathbf{R})f_{\lambda}^{2}(\bm{\epsilon}_{\mathbf{k}\lambda})_{i}(\bm{\epsilon}_{\mathbf{k}\lambda})_{j}\,e^{i\mathbf{k}\cdot\mathbf{R}}\,.$

Substituting Eqs. (20) and (34) into the previous equation, we obtain, after
some straightforward algebra,

$\displaystyle U_{\text{total}}$ $\displaystyle=$
$\displaystyle-\frac{\hbar}{4\pi^{2}R^{3}}\text{Re}\int
d^{3}k\,\alpha_{A}(\omega_{k})\alpha_{B}(\omega_{k})\omega_{k}e^{-ikR}e^{i\mathbf{k}\cdot\mathbf{R}}$
(38) $\displaystyle\times$
$\displaystyle\left\\{\left[\frac{2\omega_{k}^{2}R^{2}}{c^{2}}-2ik-2-k^{2}R^{2}\sin^{2}\theta_{k}\right.\right.$
(40) $\displaystyle+$
$\displaystyle\left.\left.3ikR\sin^{2}\theta_{k}+3\sin^{2}\theta_{k}\right]\right.$
(42) $\displaystyle+$
$\displaystyle\left.\frac{c^{2}\mu^{2}}{\omega_{k}^{2}}\left[\frac{\omega_{k}^{2}R^{2}}{c^{2}}-ikR-1-k^{2}R^{2}\cos^{2}\theta_{k}+\right.\right.$
(44) $\displaystyle+$
$\displaystyle\left.\left.3ikR\cos^{2}\theta_{k}+\cos^{2}\theta_{k}\right]\right\\},$
(46)

where we have defined $\theta_{\mathbf{k}}$ as the angle between
$\mathbf{\hat{k}}$ and $\mathbf{R}$, and took the continuum limit
$V\rightarrow\infty$. Defining $x\equiv kR$ and performing the angular
integrals, the interaction energy reads

$\displaystyle U_{\text{total}}=\frac{-\hbar c}{\pi
R^{7}}\int_{0}^{\infty}\,\frac{dx\,x\,\alpha_{A}\quantity(\omega_{k}(x))\alpha_{B}\quantity(\omega_{k}(x))}{2\sqrt{x^{2}+\quantity(\mu
R)^{2}}}\times$ $\displaystyle\left\\{\left(2x^{4}+2x^{2}\left(2\quantity(\mu
R)^{2}-5\right)+3\left(\quantity(\mu R)^{4}+2\right)\right)\sin(2x)\right.$
$\displaystyle+\left.4x\left(x^{2}-3\right)\cos(2x)\right\\},$ (47)

At this point it is important to make a self-consistent check of our results
and reobtain as a particular case the result with usual QED. In fact, by
taking the $\mu\to 0$ limit we re-obtain the well known result for the
interaction energy within Maxwell QED (see Eq. 3.85 of [1]).

Integrating by residues in the complex $k$ plane, we may transform the
integration in Eq. (III.1) into an imaginary $k$ one, thus exchanging
oscillatory terms by evanescent ones and recasting that equation in a form
amenable to numerical integration (see details in appendix B). By expressing
our results in terms of the imaginary frequency $\omega=i\frac{c}{R}\zeta$
(with $\zeta\in\mathbb{R}$) we obtain

$\displaystyle U_{\text{total}}\\!=\\!\frac{-\hbar c}{2\pi
R^{7}}\int_{0}^{\infty}d\zeta$ $\displaystyle\Bigg{\\{}$ $\displaystyle
e^{-2\sqrt{\zeta^{2}+(\mu R)^{2}}}f\left(\sqrt{\zeta^{2}+(\mu
R)^{2}}\right)\times$ (48) $\displaystyle\times$
$\displaystyle\alpha_{A}\left(i\frac{c}{R}\zeta\right)\alpha_{B}\left(i\frac{c}{R}\zeta\right)\Bigg{\\}}\,,$
(50)

where

$f(\xi)\\!=\\!6+3\quantity(\mu R)^{4}+12\xi+(10\\!-\\!4\quantity(\mu
R)^{2})\xi^{2}+4\xi^{3}+2\xi^{4}.$ (51)

The mass of the photon has a crucial role in the calculation, since it
introduces a cut due to the square root in the dispersion relation
$\omega_{k}=c\sqrt{k^{2}+\mu^{2}}$.

###  III.2 Qualitative discussion

A striking feature of Eq. (48) concerns the relevant physical scales. Pure
dimensional analysis shows the existence of three timescales, corresponding to
three typical frequencies: $c/R$, $c\mu$, and the dominant atomic transition
frequency $\omega_{0}$. Nonetheless, an inspection of the integrand of Eq.
(48) reveals that for the dispersive interaction only two scales are involved:
$\omega_{0}$, which is present in the polarizability functions and

$\omega_{F}\sim c\sqrt{\frac{1}{R^{2}}+\mu^{2}}\,.$ (52)

which is contained in the field functions $f$ and the decaying exponential.
Let us now employ these scales in order to determine the values of $\zeta$
which dominate the integral in Eq. (50). Equation (17) shows that the atoms
become transparent ($\alpha_{A,B}\approx 0$) at large imaginary frequencies,
$(c/R)\,\zeta\gg\omega_{0}$, suppressing large values of $\zeta$. The field
modes are also suppressed for high frequencies - this is the main advantadge
of the rotation in the complex plane discussed in Appendix B. This yields to
an exponential attenuation, cutting off $(c/R)\,\zeta\gg\omega_{F}$ in Eq.
(50). These frequencies correspond to large wavevectors $k\gg R^{-1}$ which
oscillate a large number of times between atoms $A$ and $B$, leading to a
decorrelation between the field at positions $\bm{r}_{A}$ and $\bm{r}_{B}$
after superposing a continuum of modes. The interplay between large-frequency
transparency and large-wavenumber decorrelation leads to two opposite
asymptotic regimes, which we discuss below.

These frequencies correspond to large wavevectors $k\gg R^{-1}$ which
oscillate a large number of times between atoms $A$ and $B$, which, upon
superposing a continuum of modes, leads to decorrelation between the field at
position $\bm{r}_{A}$ and the atomic dipole at position $\bm{r}_{B}$. The
interplay between large-frequency transparency and large-wavenumber
decorrelation leads to two opposite asymptotic regimes, which we discuss
below.

When $\omega_{0}\gg\omega_{F}$, the contribution of the high frequency in
integrand (48) is suppressed by decorrelation in the response of the field
before atomic high-frequency transparency sets in. Therefore, interatomic
dipole-dipole correlations are limited by the retardation of the
electromagnetic field. Moreover, from Eq. (17), one may approximate
$\alpha(\omega)\approx\alpha(0)$, where the latter denotes the static
polarizability. Physically this is due to the fact that in this regime the
atom easily follows the much slower field oscillations. In other words, in
this regime the correlation between the atomic dipoles are limited by the
delay in electromagnetic interaction. This case is denominated the retarded
limit, known as the Casimir-Polder regime in the case of Maxwell
electrodynamics. Note that in Proca QED this regime is possible only for small
masses $\mu\ll\omega_{0}/c$.

In the opposite limit, $\omega_{0}\ll\omega_{F}$, large values of $\zeta$ are
suppressed due to high-frequency transparency before field-decorrelation or
retardation effects become relevant. We can therefore approximate the
imaginary frequency to be $\zeta\approx 0$ for the field — that is, everywhere
in the integrand of Eq.(50) except in the polarizabilities. In this case we
may treat the electromagnetic field in the electrostatic approximation, thus
recovering the results of section II, as we shall explicitly confirm in
sectionIII.4. From Eq. (52) we see that the non-retarded regime applies to a
wider range of frequencies in Proca than in Maxwell electrodynamics, since the
condition for $R$ is weaker. In particular, if $\mu c\gg\omega_{0}$ we are
always in the non-retarded limit, regardless of distance. Naturally, however,
for $R\gg\mu^{-1}$ the interaction is exponentially suppressed and thus very
weak.

By defining the phase velocity
$v_{p}(k)\equiv\omega(k)/k=c\sqrt{(1+\mu^{2}/k^{2})}$, the non-retarded regime
can be written as $\omega_{0}R\ll v(k_{\rm field})$ where $k_{\rm field}=1/R$
(with an analogous expression for the retarded regime). In this way, the
condition defining each regime in Proca QED is the same as in the Maxwellian
case as long as in the latter we substitute $c$ by the phase velocity of the
photon evaluated at the wavenumber $1/R$. Nonetheless there are now remarkable
differences. In Maxwell, $t_{\gamma}=R/c$ is the time it takes to a photon to
travel the distance separating the atoms and corresponds to the retardation
time of the interaction. Therefore, the non-retarded condition means that
$t_{\gamma}$ must be much smaller than the atomic typical time scale
$1/\omega_{0}$, representing the dipole’s fluctuation time. No such simple
interpretation holds in Proca, since due to dispersion there is not a single
retardation time describing the interaction and the phase velocity can not be
interpreted as the velocity of propagation of the electromagnetic wave.
Indeed, phase velocity is superluminal in Proca electrodynamics. This poses no
paradox since we are dealing only with virtual photons and not with any
information propagation. The key concept here is that the photon mass
increases the photon frequency for every wavelength and thus retardation
effects are less relevant in Proca electrodynamics than in the Maxwellian
case, since it is the slower dynamics which is more relevant to limit the
dipole-dipole correlation. We now turn to the detailed calculation which
supports the above statements.

###  III.3 Retarded regime ($\omega_{f}\ll\omega_{0}$)

In this regime, we may substitute the atomic polarizabilities given in Eq.
(17) by their static values $\alpha(\omega=0)$. With this approximation, the
integral (50) takes the form

$\displaystyle U_{\text{Ret}}$ $\displaystyle=$ $\displaystyle\frac{-\hbar
c}{2\pi R^{7}}\alpha_{A}(0)\alpha_{B}(0)$ (53) $\displaystyle\times$
$\displaystyle\int_{0}^{\infty}d\zeta\,e^{-2\sqrt{\zeta^{2}+(\mu
R)^{2}}}f\left(\sqrt{\zeta^{2}+(\mu R)^{2}}\right)\,.$ (55)

After performing the integration we are left with

$\displaystyle U_{\text{Ret}}$ $\displaystyle=-\frac{\hbar
c\alpha_{A}\alpha_{B}}{4\pi R^{7}}\left\\{8\quantity(\mu
R)^{2}\left(\quantity(\mu R)^{2}+3\right)K_{0}(2\mu R)+\right.$ (56)
$\displaystyle+$ $\displaystyle\left(2\quantity(\mu R)^{5}+31\quantity(\mu
R)^{3}+24(\mu R)\right)K_{1}(2\mu R)+$ (58) $\displaystyle+$
$\displaystyle\left.22\quantity(\mu R)^{2}K_{2}(2\mu R)\right\\},$ (60)

where $K_{\nu}(x)$ are the modified Bessel functions of the second kind.

As $\mu R\rightarrow 0$, the term $(\mu R)K_{0}(2\mu R)$ in Eq. (56) becomes
subdominant, while, using
$K_{\alpha}(x)\sim\frac{\Gamma(\alpha)}{2}\left(\frac{2}{x}\right)^{\alpha}$,
valid for $\alpha>0$, we see that the term inside the brackets becomes
$\sim(24+22)/2=23$, thus re-obtaining the famous result by Casimir and Polder
[56] hereafter denoted $U_{\rm CP}$. Working out the next term in the small
mass limit $\mu R\ll 1$ we obtain

$U=U_{\rm CP}\left(1-\frac{15\quantity(\mu
R)^{2}}{46}\right)+\mathcal{O}[\quantity(\mu R)^{4}].$ (61)

Figure 2: Interaction potential normalized by the Casimir-Polder potential
$U_{\text{CP}}$ as a function of the dimensionless radius parameter $\mu R$.
Figure 3: Complete interaction energy $U_{\text{total}}$, normalized by the
retarded interaction $U_{R}$, as a function of $\omega_{0}R/c$, for identical
two-level atoms with transition frequency $\omega_{0}R$ for different values
of the photon mass.

As we had obtained for the non-retarded regime, the first correction is of
order $(\mu R)^{2}$ and weakens the interaction. This is due to the fact that
the mass increases the frequency for every wavevector, thus weakening the
correlation between the atoms. This behavior is illustrated in Fig. 2 where we
plot the Casimir-Polder retarded energy normalized by the Casimir-Polder
result in terms of $\mu R$.

In the opposite limit, $\mu R\gg 1$, we have an exponential decay given by the
asymptotic form of the modified Bessel functions,

$U=-\frac{\hbar
c\alpha_{A}(0)\alpha_{B}(0)\mu^{9/2}}{4\pi^{1/2}R^{5/2}}e^{-2\mu
R}\left\\{1+\mathcal{O}\left(\frac{1}{\mu R}\right)\right\\}$ (62)

The previous result shows that, for a given distance between the atoms, their
interaction potential decreases exponentially with the mass of the photon, as
expected, since the quantum fluctuations of the field are drastically
suppressed as the photon mass increases. Note also the fractional power law
dependence with the distance between the atoms. This fact resembles somehow
the correction to the Coulomb interaction which appears in the so called
Uehling potential (see [57] for more detail). In the last case, the electron
loop in the vacuum polarization Feynman diagram introduces the electron mass
scale ($m_{e}$) into the problem and leads, for large $R$ ($R\gg 1/m_{e}$), to
a correction to the Coulomb potential which is suppressed expeonentially with
$m_{e}R$ and which falls with the distance with the same fractional power law
that we have obtained for a massive photon.

In Fig. 3 we plot the exact interaction energy given by Eq.(50) normalized by
$U_{\rm Ret}$ written in Eq.(56) for identical two level atoms as a function
of $\omega_{0}R$, where $\omega_{0}$ denotes the atomic transition frequency.
This is done by assuming that only a single term is present in the
polarizability defined in Eq.(17). Notice that the Casimir-Polder
approximation overestimates the interaction energy, which is expected since
the static polarizability is higher than the dynamical one far from resonance
(which is always the case in the retarded regime). Another striking feature is
that the higher the transition frequency is, for a fixed distance, the worse
the retarded approximation becomes, along the lines discussed in the previous
subsection.

###  III.4 Non-retarded regime ($\omega_{f}\gg\omega_{0}$)

Here we take the opposite limit and, as discussed in section III.2, we take
the zero frequency limit for the field degrees of freedom, setting $\zeta=0$
in Eq. (50), thus obtaining for interaction energy in the non-retarded regime
the expression

$\displaystyle U_{\text{NR}}(R)\\!\\!=\\!\\!\frac{-\hbar ce^{-2\mu R}}{2\pi
R^{7}}p(\mu R)I\,,$ (63)

where $p$ is defined in Eq. (13) and

$\displaystyle I$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\zeta\,\alpha_{A}\left(\frac{i\zeta
c}{R}\right)\alpha_{B}\left(\frac{i\zeta c}{R}\right)$ (64) $\displaystyle=$
$\displaystyle\frac{4}{9\hbar^{2}}\sum_{r,s\neq
0}|\bm{d}_{A}^{0r}|^{2}|\bm{d}_{B}^{0r}|^{2}\int_{0}^{\infty}\frac{d\zeta\,\omega_{r0}\omega_{s0}}{\left(\omega_{r0}^{2}+\frac{R^{2}\zeta^{2}}{c^{2}}\right)\left(\omega_{s0}^{2}+\frac{R^{2}\zeta^{2}}{c^{2}}\right)}$
(66) $\displaystyle=$ $\displaystyle\frac{2\pi R}{9\hbar c}\Lambda\,,$ (68)

where $\Lambda$ is defined in Eq. (14).

Substituting the result in Eq. (68) into Eq. (63) we reobtain Eq. (12) of
section II. Therefore, in Proca we see that when the retardation of the
electromagnetic field can be neglected we do not need to quantize the
electromagnetic field, as it is also true in Maxwell. Nonetheless, in Proca we
have an extra subtlety worth mentioning. In section II we did not quantize the
classical degrees of freedom for the field, so that the quantum fluctuations
came entirely from the atomic dipole fluctuations. In particular, since we are
in the electrostatic regime ($\omega=0$), the field produced by the
fluctuating dipoles is evanescent. In section III on the other hand, our
Hamiltonian (28) contained only field fluctuations with the atomic dipoles
being induced by them, as in Eq. (25). Here we have only propagating modes for
the fields, which requires frequencies greater than $\mu c$, in marked
contrast with the previous calculation. Nevertheless, as we have demonstrated
in this section, the results of the electrostatic treatment of section II are
contained in the full quantum electrodynamics analysis. This is physically
expected since we can choose whether we consider the interatomic interaction
as fluctuating dipoles which generates electric field or the other way around,
with the vacuum electromagnetic fluctuations playing the leading role. We
could as well start from a democratic Hamiltonian which involved fluctuations
of both degrees of freedom but then we would have to perform a tedious 4th
order perturbation theory calculation [53].

Figure 4: Complete interaction energy $U_{\text{total}}$, normalized by the
non-retarded regime interaction $U_{NR}$, as a function of $\omega_{0}R/c$ for
different values of the photon mass. In this analysis, atoms $A$ and $B$ are
considered to be identical two-level systems of energy gap $\omega_{0}$.

The ratio between the full result from Eq. (50) for identical two-level atoms
and the interaction energy in the non-retarded approximation is shown as a
function of $\mu R$ in Fig. 4. We note that, as $\mu c/\omega_{0}$ increases,
the non-retarded regime extends to larger values of $\omega_{0}R/c$, but the
interaction energy itself becomes more strongly suppressed for $R\gtrsim
c/\omega_{0}$.

The attentive reader may have noticed the similarity between Figs. 1 and 2,
which show the ratio between the massive and massless interaction energies,
$U(\mu R)/U(\mu R=0)$, evaluated in the non-retarded and the retarded
approximations, respectively. As shown in Fig. 5, the curves for this ratio as
a function of the dimensionless parameter $\mu R$ are indeed remarkably
similar for both regimes, even though the two approximations correspond to
opposite limits. In the particular case where $\mu R\ll 1$, this can be
readily verified by comparing Eqs.( 15) and (61), since the coefficients of
the quadratic correction in $\mu R$ differ by less than 1% iiiWe note that the
asymptotic expressions for $\mu R\gg 1$ have different behaviors for the two
approximations. However, this difference only becomes pronounced for small
values of $U(\mu R)/U(\mu R=0)$, which explains why they are not visible in
Fig. 5.. This happens because for a given interatomic distance, the ratio
$U_{\rm total}(\mu R)/U_{\rm total}(\mu R=0)$, evaluated from Eq. (50) is very
nearly independent of the transition frequency $\omega_{0}$. Indeed, by
varying $\omega_{0}$ between $10^{-4}c/R$ and $10^{4}c/R$ for every distance
we obtain the grey band depicted in Fig. 5 connecting the retarded curve
(which is the limite for $\omega_{0}\rightarrow\infty)$ with the non-retarded
curve (corresponding to the $\omega_{0}\rightarrow 0$ limit).

Figure 5: Ratio between the interaction energies for massive and massless QED,
as a function of the dimensionless radius parameter $\mu R$ for two regimes:
the non-retarded (blue dashed line) and the retarded (purple dot-dashed line)
approximations. The gray band visible on the inset displays the curves $U_{\rm
total}(\mu R)/U_{\rm total}(\mu R=0)$ for $\omega_{0}$ in the range
$10^{-4}c/R-10^{4}c/R$.

##  IV Conclusions and final remarks

In this paper we have derived the dispersive interaction energy between two
polarizable atoms, with no permanent dipole moments, within Proca
electrodynamics for any distance regime. We emphasize that we employed a very
convenient Hamiltonian, introduced in Ref.[1], which allows for a first order
perturbative calculation.

We have found that, despite the inclusion of an extra energy scale
corresponding to the photon mass, dispersive interactions within Proca
electrodynamics are still marked by the existence of no more than two
asymptotic regimes, corresponding to the London and Casimir-Polder limits of
Maxwellian electrodynamics. A timescale consisting of interatomic distance
over phase velocity may be used to distinguish between these two regimes, but,
in contrast to Maxwell’s theory, the nonlinear relationship between frequency
and wavenumber in the massive case prevents its unambiguous physical
interpretation as a retardation time. This is due to the fact that the phase
velocity is both frequency-dependent and superluminal in the massive theory.
We have computed closed-form analytical expressions that approximate the
dispersive interaction energy in the retarded and non-retarded regimes and
presented a detailed discussion on their differences with respect to the case
of a massless photon. Our results reveal an enhancement of the non-retarded
regime, which extends over a broader range of interatomic distances in the
presence of a photon mass.

We have also shown that the effects of the photon mass weaken the interaction
energy. More strikingly, we have shown that, for a given interatomic distance,
this weakening is very similar for both the retarded and the non-retarded
regimes. This means that the mass correction is nearly independent of the
atomic transition frequency. As a perspective, it would be interesting to
explore if this still holds beyond the atomic scenario, analyzing for instance
the Casimir-Proca force between dielectric media.

Our results should be relevant not only to beyond the standard model physics
scenarios, but also for any theories or effective descriptions featuring
massive vector bosons or an effective photon mass. An effective photon mass
may emerge, for instance, inside electrolyte solutions, metals or waveguides
[23, 24]. Other massive vector bosons can be found in effective descriptions
of nuclear matter, such as relativistic mean-field models of nuclear
interactions [44, 45, 41, 42, 43]. In the nuclear context, the possible
polarization of mesons into a baryon-numer dipole could lead to new couplings
with vector mesons, with potentially interesting phenomenological consequences
[46].

We hope our findings will inform future investigations on potential observable
effects of a finite, albeit small, photon mass [53]. Because bounds on the
mass of the photon are rather small, the detection of a photon mass via van
der Waals interactions would be far from trivial. Nonetheless, possible
implications of the finite-mass modifications in the condensation of gases and
other physical phenomena could potentially prove more revealing of the photon
mass then the van der Waals interaction energy in itself. It would be
interesting to explore this possibility in future work.

###### Acknowledgements.

The authors would like to thank Dr. Pedro C. Malta for the careful reading of
this manuscript. This study was financed in part by the Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code
001. C.A.D.Z. are partially supported by Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq) under the grant 310703/2021-2. C.A.D.Z. is
also partially supporteprocessação Carlos Chagas Filho de Amparo à Pesquisa do
Estado do Rio de Janeiro (Faperj) under Grant E-26/201.447/2021 (Programa
Jovem Cientista do Nosso Estado). M.H. was supported in part by the National
Science Foundation (NSF) within the framework of the MUSES Collaboration,
under grant number OAC-2103680, and by Universidade Estadual do rio de
Janeiro, within the Programa de Apoio à Docência (PAPD).

##  Appendix A The Proca electric dipole field

In this appendix we obtain the Proca field generated by an electric dipole,
following an alternative approach to the one presented in Ref.[59]. We work in
the Fourier frequency domain, but keep the dependence in spatial variables.
With this choice the potential quadrivector $(\phi,\bm{A})$ satisfies in Proca
electrodynamics the equation

$\left(\bm{\nabla}^{2}+\frac{\omega^{2}}{c^{2}}-\mu^{2}\right)A^{\mu}(\bm{r})=-\frac{4\pi}{c}j^{\mu}(\bm{r})\,,$
(69)

which can be readily solved with the aid of the Green function satisfying

$\left(\bm{\nabla}^{2}+\frac{\omega^{2}}{c^{2}}-\mu^{2}\right)G(\bm{r},\bm{r}^{\prime},\omega)=-4\pi\delta(\bm{r}-\bm{r}^{\prime})\,,$
(70)

whose solution is given by

$G(\mathbf{r},\mathbf{r^{\prime}},\omega)=\frac{e^{\mp
ik|\mathbf{r}-\mathbf{r^{\prime}}|}}{|\mathbf{r}-\mathbf{r^{\prime}}|}\,,$
(71)

where we have definied $k:=\sqrt{\frac{\omega^{2}}{c^{2}}-\mu^{2}}$ and the
minus (plus) sign refers to the solution with (positive) negative frequency.
For the purposes of this paper, we need only propagating frequencies
$\omega>\mu c$ for these are the only ones present in the vacuum
electromagnetic field. The Proca field produced by a dipole follows from
solving Eq. (69) with

$j^{\mu}=(-c\bm{d}(\omega)\cdot\bm{\nabla}\delta(\bm{r}-\bm{r}^{\prime}),-i\omega\bm{d}(\omega)\delta(\bm{r}-\bm{r}^{\prime}))\,,$
(72)

which corresponds to a current 4-vector of a point dipole[60] $\bm{d}(\omega)$
at position $\bm{r}^{\prime}$. Hence, the potentials are given by

$\displaystyle\phi(\bm{r},\omega)$ $\displaystyle=$
$\displaystyle\bm{d}(\omega)\cdot\bm{\nabla}^{\prime}G(\bm{r},\bm{r}^{\prime},\omega)\,,$
(73) $\displaystyle\bm{A}(\bm{r},\omega)$ $\displaystyle=$
$\displaystyle-\frac{i\omega}{c}\bm{d}(\omega)G(\bm{r},\bm{r}^{\prime},\omega)\,.$
(74)

Hence, the Proca electric field for a dipole is given by

$E_{i}(\mathbf{r},\omega)=d_{j}(\omega)\left(\partial_{i}\partial_{j}+\delta_{ij}\omega^{2}\right)G(\mathbf{r},\mathbf{r^{\prime}},\omega)\,,$
(75)

where we have employed $\partial_{j}^{\prime}=-\partial_{j}$, which is valid
since $G$ \- given in Eq. (71) - is a function only of
$\bm{r}-\bm{r}^{\prime}$. By evaluating the derivatives present in the
previous equation we obtain Eq. (24) of the main paper. Finally, the
electrostatic case is obtained by taking the static limit $\omega\rightarrow
0$, which is equivalent to take $k=\mp\mu$. The choice of the sign is made to
comply with the boundary condition $G\rightarrow 0$ for $r\rightarrow\infty$.
With this choice, Eq. (24) furnishes Eq. (4) of section II.

##  Appendix B The integral of Sec. III

In this appendix we will explicitly calculate integrals of the form

$I=A\int_{0}^{\infty}\frac{dx}{\sqrt{x^{2}+y^{2}}}\left[P_{s}(x)\sin(2x)+P_{c}(x)\cos(2x)\right],$
(76)

where $A$ is some real constant, $P_{s}(x)$ is a real odd function and
$P_{c}(x)$ is a real even function with no branch cuts. Those integrals are
clearly divergent on the UV regime, so we will need to perform it on the
extended complex plane.

The first step is simply to use the parity of the integrand to write the above
integral as

$I=\frac{A}{2}\text{Re}\int_{-\infty}^{\infty}\frac{dx}{\sqrt{x^{2}+y^{2}}}\left[P_{c}(x)-iP_{s}(x)\right]e^{2ix}.$
(77)

We then move to the complex plane, defining a complex variable $z$ such that
$\real z=x$, so that integral (77) lies on the real $z$ line.

Figure 6: Contour in the complex plane used to compute Eq. (77).

The above integrand has two branch cuts on the complex plane (one running from
$x=iy$ to $x=i\infty$ and the other running from $x=-iy$ to $x=-\infty$ ), as
shown in Fig. 6. To apply the residue theorem, we avoid these branch cuts, by
employing the contour shown in Fig. 6. We choose to close the contour on the
upper half plane, so that the segments $C_{R}^{\pm}$ on this figure evaluate
to zero at complex infinity. The only parts of the contour that are non-zero
are $I$, $u^{+}$ and $u^{-}$, and the residue theorem implies that

$\displaystyle\int_{I}g(z)\,dz+\int_{u_{+}}g(z)\,dz+\int_{u_{-}}g(z)\,dz=0,$
(78)

where

$g(z)=\frac{\tilde{P}(z)\,e^{2iz}}{\sqrt{z^{2}+y^{2}}}$ (79)

and

$\tilde{P}(z)=P_{c}(z)-iP_{s}(z).$ (80)

The branch cut is chosen such that the square root has positive imaginary part
on the right side (the $u^{-}$ line) and negative imaginary part on the left
side (the $u^{-}$ line). Hence,

$\displaystyle\int_{u_{+}}g(z)dz=\frac{A}{2}\int_{\infty}^{y}\frac{(idu)}{i\sqrt{u^{2}-y^{2}}}\tilde{P}(iu)e^{-2u}$
$\displaystyle=-\frac{A}{2}\int_{y}^{{\infty}}\frac{du}{\sqrt{u^{2}-y^{2}}}\tilde{P}(iu)\,e^{-2u}$
(81)

and

$\displaystyle\int_{u_{-}}g(z)dz=\frac{A}{2}\int_{y}^{\infty}\frac{(idu)}{-i\sqrt{u^{2}-y^{2}}}\tilde{P}(iu)e^{-2u}$
$\displaystyle=-\frac{A}{2}\int_{y}^{{\infty}}\frac{du}{\sqrt{u^{2}-y^{2}}}\tilde{P}(iu)e^{-2u}$
(82) $\displaystyle=\int_{u_{+}}g(z)dz$ (83)

So that integral $I$ becomes:

$I=-A\,\text{Re}\int_{y}^{{\infty}}\frac{du}{\sqrt{u^{2}-y^{2}}}\tilde{P}(iu)\,e^{-2u}.$
(84)

## References

  * Milonni [1994] P. W. Milonni,  _The Quantum Vacuum, an introduction to quantum electrodynamics. Academic Press Inc._ (Academic Press, 1994).
  * Israelachvili [2011] J. N. Israelachvili,  _Intermolecular and Surface Forces,_ (Academic Press, 2011).
  * Buhmann [2013a] S. Y. Buhmann,  _Dispersion Forces I: Macroscopic quantum electrodynamics and ground-state Casimir, Casimir–Polder and van der Waals forces_ , Vol. 247 (Springer, 2013).
  * Buhmann [2013b] S. Y. Buhmann,  _Dispersion Forces II: Many-Body Effects, Excited Atoms, Finite Temperature and Quantum Friction_ , Vol. 248 (Springer, 2013).
  * Milton [2001] K. A. Milton,  _The Casimir Effect: Physical Manifestations of Zero-Point Energy_ (World Scientific, 2001).
  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Castillo-López  _et al._ [2024] S. G. Castillo-López, R. Esquivel-Sirvent, G. Pirruccio, and C. Villarreal, Physics 6, 394 (2024).
  * Castillo-López  _et al._ [2022] S. G. Castillo-López, R. Esquivel-Sirvent, G. Pirruccio, and C. Villarreal, Scientific Reports 12, 2905 (2022).
  * Fosco  _et al._ [2024] C. D. Fosco, F. C. Lombardo, and F. D. Mazzitelli, Physics 6, 290 (2024).
  * Fosco and Mazzitelli [2020] C. D. Fosco and F. D. Mazzitelli, Physical Review D 101, 045012 (2020).
  * Fosco  _et al._ [2016] C. D. Fosco, F. C. Lombardo, and F. D. Mazzitelli, Physical Review D 93, 125015 (2016).
  * Romaniega  _et al._ [2023] C. Romaniega, J. Munoz-Castaneda, and I. Cavero-Pelaez, Physical Review D 107, 025002 (2023).
  * Cavero-Peláez  _et al._ [2021] I. Cavero-Peláez, J. M. Munoz-Castaneda, and C. Romaniega, Physical Review D 103, 045005 (2021).
  * Shajesh  _et al._ [2016] K. V. Shajesh, I. Brevik, I. Cavero-Peláez, and P. Parashar, Physical Review D 94, 065003 (2016).
  * Shajesh  _et al._ [2017] K. Shajesh, P. Parashar, I. Cavero-Peláez, J. Kocik, and I. Brevik, Physical Review D 96, 105010 (2017).
  * Li  _et al._ [2022] Y. Li, K. A. Milton, P. Parashar, G. Kennedy, N. Pourtolami, and X. Guo, Physical Review D 106, 036002 (2022).
  * Milton  _et al._ [2017] K. A. Milton, P. Kalauni, P. Parashar, and Y. Li, Physical Review D 96, 085007 (2017).
  * Milton  _et al._ [2019] K. A. Milton, P. Kalauni, P. Parashar, and Y. Li, Physical Review D 99, 045013 (2019).
  * Autumn  _et al._ [2002] K. Autumn, M. Sitti, Y. A. Liang, A. M. Peattie, W. R. Hansen, S. Sponberg, T. W. Kenny, R. Fearing, J. N. Israelachvili, and R. J. Full, Proceedings of the National Academy of Sciences 99, 12252 (2002).
  * Lamoreaux [2007] S. K. Lamoreaux, Physics Today 60, 40 (2007).
  * Scheeres  _et al._ [2010] D. J. Scheeres, C. M. Hartzell, P. Sánchez, and M. Swift, Icarus 210, 968 (2010).
  * Mahanty and Ninham [1976] J. Mahanty and B. W. Ninham,  _Dispersion Forces_ (Academic Press, 1976).
  * Chiao and Boyce [1999] R. Y. Chiao and J. Boyce, Physical Review A 60, 4114 (1999).
  * de Melo e Souza  _et al._ [2015] R. de Melo e Souza, W. J. M. Kort-Kamp, F. S. S. Rosa, and C. Farina, Phys.Rev. A 91, 052708 (2015).
  * Grigorio  _et al._ [2012] L. S. Grigorio, M. S. Guimaraes, R. Rougemont, C. Wotzasek, and C. A. D. Zarro, Phys. Rev. D 86, 027705 (2012), arXiv:1202.3798 [hep-th] .
  * Guimaraes  _et al._ [2013] M. S. Guimaraes, R. Rougemont, C. Wotzasek, and C. A. D. Zarro, Phys. Lett. B 723, 422 (2013), arXiv:1209.3073 [hep-th] .
  * Reinosa [2024] U. Reinosa (2024) arXiv:2404.06118 [hep-ph] .
  * Rougemont  _et al._ [2015] R. Rougemont, J. Noronha, C. A. D. Zarro, C. Wotzasek, M. S. Guimaraes, and D. R. Granado, JHEP 07, 070, arXiv:1505.02442 [hep-th] .
  * Pospelov [2009] M. Pospelov, Phys. Rev. D 80, 095002 (2009), arXiv:0811.1030 [hep-ph] .
  * Spallicci  _et al._ [2021] A. D. A. M. Spallicci, J. A. Helayël-Neto, M. López-Corredoira, and S. Capozziello, Eur. Phys. J. C 81, 4 (2021), arXiv:2011.12608 [astro-ph.CO] .
  * Bordag  _et al._ [1998] M. Bordag, B. Geyer, G. Klimchitskaya, and V. Mostepanenko, Physical Review D 58, 075003 (1998).
  * Fischbach  _et al._ [2001] E. Fischbach, D. Krause, V. Mostepanenko, and M. Novello, Physical Review D 64, 075010 (2001).
  * Mostepanenko [2016] V. Mostepanenko, International Journal of Modern Physics A 31, 1641020 (2016).
  * Adelberger  _et al._ [2007] E. Adelberger, G. Dvali, and A. Gruzinov, Physical Review Letters 98, 010402 (2007).
  * Chibisov [1976] G. V. Chibisov, Soviet Physics Uspekhi 19, 624 (1976).
  * Accioly  _et al._ [2010] A. Accioly, J. Helayël-Neto, and E. Scatena, Physical Review D 82, 065026 (2010).
  * Malta and Helayël-Neto [2022] P. Malta and J. Helayël-Neto, Physical Review D 106, 116014 (2022).
  * Malta [2023] P. Malta, Physical Review D 107, 115018 (2023).
  * Goldhaber and Nieto [1971] A. S. Goldhaber and M. M. Nieto, Reviews of Modern Physics 43, 277 (1971).
  * Goldhaber and Nieto [2010] A. S. Goldhaber and M. M. Nieto, Reviews of Modern Physics 82, 939 (2010).
  * Glendenning [1997] N. K. Glendenning,  _Compact Stars: Nuclear Physics, Particle Physics and General Relativity_, Astronomy and Astrophysics Library (Springer New York, 1997).
  * Horowitz and Piekarewicz [2001] C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001), arXiv:astro-ph/0010227 .
  * Chen and Piekarewicz [2014] W.-C. Chen and J. Piekarewicz, Phys. Rev. C 90, 044305 (2014), arXiv:1408.4159 [nucl-th] .
  * Walecka [1974] J. D. Walecka, Annals Phys. 83, 491 (1974).
  * Boguta [1983] J. Boguta, Phys. Lett. B 120, 34 (1983).
  * Pisarski [2021] R. D. Pisarski, Phys. Rev. D 103, L071504 (2021), arXiv:2101.05813 [nucl-th] .
  * Willey [1978] R. Willey, Physical Review D 18, 270 (1978).
  * Fujii and Kharzeev [1999] H. Fujii and D. Kharzeev, Phys. Rev. D 60, 114039 (1999), arXiv:hep-ph/9903495 .
  * Proca [1936] A. Proca, Journal de Physique et le Radium 7, 347 (1936).
  * Davies and Unwin [1981] P. Davies and S. Unwin, Physical Letters B 98, 274 (1981).
  * Barton and Dombey [1984] G. Barton and N. Dombey, Nature 311, 336 (1984).
  * Barton and Dombey [1985] G. Barton and N. Dombey, Annals of Physics 162, 231 (1985).
  * Mattioli  _et al._ [2019] L. Mattioli, A. Frassino, and O. Panella, Physical Review D 100, 116023 (2019).
  * Craig and Thirunamachandran [1998] D. P. Craig and T. Thirunamachandran,  _Molecular quantum electrodynamics: an introduction to radiation-molecule interactions_ (Courier Corporation, 1998).
  * Greiner and Reinhardt [1996] W. Greiner and J. Reinhardt,  _Field quantization_ (Springer Science & Business Media, 1996).
  * Casimir and Polder [1948] H. B. G. Casimir and D. Polder, Physical Review 73, 360 (1948).
  * Peskin and Schroeder [1995] M. E. Peskin and D. V. Schroeder,  _An Introduction to quantum field theory_ (Addison-Wesley, Reading, USA, 1995).
  * [58] We note that the asymptotic expressions for $\mu R\gg 1$ have different behaviors for the two approximations. However, this difference only becomes pronounced for small values of $U(\mu R)/U(\mu R=0)$, which explains why they are not visible in Fig. 5.
  * Dragulin  _et al._ [2008] P. Dragulin, P. Leung, _et al._ , Physical Review E 78, 026605 (2008).
  * Pitombo  _et al._ [2021] R. S. Pitombo, M. Vasconcellos, C. Farina, and R. de Melo e Souza, European Journal of Physics 42, 025202 (2021).