# Equivalence principle violation in nonminimally coupled gravity and
constraints from Lunar Laser Ranging

Riccardo March<EMAIL_ADDRESS>Istituto per le Applicazioni del Calcolo,
CNR, Via dei Taurini 19, 00185 Roma, Italy  INFN - Laboratori Nazionali di
Frascati (LNF), Via E. Fermi 40, Frascati 00044 Roma, Italy  Orfeu Bertolami
<EMAIL_ADDRESS>Departamento de Física e Astronomia,
Faculdade de Ciências da Universidade do Porto, Rua do Campo Alegre 687,
4169-007 , Porto, Portugal  Centro de Física das Universidades do Minho e do
Porto, Rua do Campo Alegre 687, 4169-007, Porto, Portugal  Marco Muccino
<EMAIL_ADDRESS>INFN - Laboratori Nazionali di Frascati (LNF), Via
E. Fermi 40, Frascati 00044 Roma, Italy  Simone Dell’Agnello
<EMAIL_ADDRESS>INFN - Laboratori Nazionali di Frascati (LNF),
Via E. Fermi 40, Frascati 00044 Roma, Italy

(today)

###### Abstract

We analyze the dynamics of the Sun-Earth-Moon system in the context of a
particular class of theories of gravity where curvature and matter are
nonminimally coupled (NMC). These theories can potentially violate the
Equivalence Principle as they give origin to a fifth force and a extra non-
Newtonian force that may imply that Earth and Moon fall differently towards
the Sun. We show, through a detailed analysis, that consistency with the bound
on Weak Equivalence Principle arising from 48 years of Lunar Laser Ranging
data, for a range of parameters of the NMC gravity theory, can be achieved via
the implementation of a suitable screening mechanism.

##  I Introduction

General Relativity (GR) can account for astrophysical and cosmological
phenomena such as the flattening of the rotation curves of galaxies and the
accelerated expansion of the Universe provided about $95\%$ of the content of
the Universe is composed of dark energy and dark matter. In principle, this
somewhat puzzling situation can be circumvented without these dark components
in the context of some alternative theories of gravity. Some alternatives
include, $f(R)$ gravity Capoz-1 ; Carroll ; Capoz-2 ; DeFTs , in which the
linear Ricci curvature scalar, $R$, in the Einstein-Hilbert action is replaced
by a more general function $f(R)$ and theories where a nonminimally coupling
(NMC) between curvature and matter is introduced BBHL . In the latter, the
Einstein-Hilbert action is replaced by two functions of curvature, $f^{1}(R)$
and $f^{2}(R)$. The function $f^{1}(R)$ is analogous to $f(R)$ gravity theory,
while the function $f^{2}(R)$ multiplies the matter Lagrangian density, which
couples nontrivially geometry and matter BBHL . This theoretical route has
been extensively examined in what concerns dark matter drkmattgal , dark
energy curraccel , inflation c , energy density fluctuations d , gravitational
waves e , cosmic virial theorem f  and black holes g . This modified theory
has also been examined through the Newton-Schrodinger approach h ; i .

In a previous work MPBD  the case of functions $f^{1}(R),f^{2}(R)$ analytic at
$R=0$ was considered, and implications of the NMC model were examined via the
perturbations to the perihelion precession by using data from observations of
Mercury’s orbit.

It turns out that NMC gravity modifies gravity as it introduces both a Yukawa
type fifth force and an extra force which depends on the spatial gradient of
the Ricci scalar. While the Yukawa force is typical also of $f(R)$ gravity,
the existence of the extra force is specific of NMC gravity BBHL ; BLP , as
the nonminimal coupling induces a non-vanishing covariant derivative of the
energy-momentum tensor. On its hand, the fifth force can give rise to static
solutions even though in the absence of pressure i .

Constraints to the NMC gravity model with analytic $f^{1},f^{2}$ functions
have been computed using the results of a geophysical experiment in Ref. MBMBD
. The idea was to consider deviations from Newton’s inverse square law in the
ocean Zum . It was found that the presence of the extra force in a fluid such
as seawater imposes more stringent constraints on the NMC gravity model than
the observation of both Mercury’s perihelion precession and lunar geodetic
precession. Hence for the NMC gravity model, Solar System constraints are
weaker than geophysical constraints.

In the present paper we look for meaningful Solar System constraints to NMC
gravity, for a function $f^{2}(R)$ which contains a term proportional to
$R^{\alpha}$, with $\alpha<0$, so that $f^{2}(R)$ is not analytic at $R=0$.
The resulting model has been used in Ref. drkmattgal  to predict the
flattening of the galaxy rotation curves, and to predict the current
accelerated expansion of the Universe curraccel .

In Ref. MPBD  the method based on the $1/c$ expansion was used to study the
NMC gravity model. However, since the $f^{2}(R)$ function is not analytic, a
different nonlinear approach has to be employed. It turns out that in Solar
System the above NMC model exhibits a screening mechanism, which is a version
of the so called chameleon mechanism KW  adapted to the NMC gravity and used
to obtain Solar System constraints MBMGDeA .

In Ref. MBMGDeA  the complications of the NMC where considered. These
generalise the nonlinear computations of Ref. HS  for the chameleon mechanism
for the gravitational field of the Sun and the corresponding calculation for
the the case of $f(R)$ gravity HS . Spherical symmetry was considered and the
constraints arising from the Cassini measurement of PPN parameter $\gamma$
Cassini  were used to constrain the parameters of the NMC gravity model. It
was shown that the chameleon solution in NMC gravity turns out to be close to
GR inside a _screening radius_ $r_{s}$ that has to be large enough, so that
$r_{s}$ either lies inside the solar convection zone, close to the top of the
zone, or it is larger. Deviations from GR are sourced by the fraction of solar
mass, including solar atmosphere, contained in the region with radii,
$r>r_{s}$, so that if $r_{s}$ lies in the convection zone then such deviations
are essentially sourced by a _thin shell_ of mass in the convection zone. This
is a typical feature of the chameleon mechanism KW .

In the present work we consider the nonlinear computations made in Ref.
MBMGDeA  and evaluate the contribution of the interactions in the NMC gravity
in the interior of the Sun, Earth and Moon. It is shown that in order to
satisfy Solar System constraints for the NMC gravity, the solution of the
chameleon mechanism must remain close to the minimizer of an effective
potential $V_{\rm eff}$ in most of the interior of massive astronomical bodies
(Sun, Earth and Moon), so that GR is approximately satisfied KW ; HS . More
specifically, for each astronomical body, the solution has to be close to the
minimizer of $V_{\rm eff}$ inside a critical radius, the screening radius,
which must be determined. If the screening radius is close to the radius of
the astronomical body, then the thin shell condition is satisfied and
deviations from GR are screened.

The paper is organized as follows. In Section II the NMC gravity model is
presented. In Section III we consider the field equations inside and around
the astronomical bodies (Sun, Earth, Moon) for the chameleon mechanism. In
Section IV we compute the solutions for the interior of the astronomical
bodies. In Section V the gravitational field of the astronomical bodies is
evaluated using the method of the images for a system of spheres. In Section
VI the dynamics of the continuous bodies is considered in order to compute the
fifth force and the extra force inside the bodies and to evaluate the yielding
jump in the pressure. Both the fifth force and the extra force are shown to be
negligible inside the screening radii. In Section VII the acceleration of
Earth and Moon due to the fifth force is computed. In Section VIII the
acceleration of Earth and Moon due to the extra NMC force is computed. In
Section IX the potential violation of the Equivalence Principle is quantified
and in section X the bound on the Weak Equivalence Principle (WEP) arsing from
the Lunar Laser Ranging (LLR) data is used to constrain the parameters of the
NMC gravity model. It is shown that the screening mechanism is successful in
ensuring that the WEP can be respected for a suitable range of model
parameters. Finally, our conclusions are presented in Section XI. Appendices A
to E contain the details of the calculations needed to obtain the results of
Section X.

##  II Nonminimally coupled gravity

The action functional of NMC gravity theory here considered is of the form
BBHL ,

$S=\int\left[\frac{1}{2}f^{1}(R)+[1+f^{2}(R)]\mathcal{L}_{m}\right]\sqrt{-g}\,d^{4}x,$
(1)

where $f^{i}(R)$ (with $i=1,2$) are functions of the Ricci curvature scalar
$R$, $\mathcal{L}_{m}$ is the Lagrangian density of matter, and $g$ is the
metric determinant. The Einstein-Hilbert action of GR is recovered by taking

$f^{1}(R)=\frac{c^{4}}{8\pi G}R,\qquad f^{2}(R)=0,$ (2)

where $G$ is Newton’s gravitational constant. We work in the Jordan frame
throughout this paper.

The first variation of the action functional with respect to the metric
$g_{\mu\nu}$ yields the field equations:

$\displaystyle\left(f^{1}_{R}+2f^{2}_{R}\mathcal{L}_{m}\right)R_{\mu\nu}-\frac{1}{2}f^{1}g_{\mu\nu}$
(3)
$\displaystyle=\left(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square\right)\left(f^{1}_{R}+2f^{2}_{R}\mathcal{L}_{m}\right)+\left(1+f^{2}\right)T_{\mu\nu},$

where $f^{i}_{R}\equiv df^{i}/dR$. The trace of the field equations is given
by

$\displaystyle\left(f^{1}_{R}+2f^{2}_{R}\mathcal{L}_{m}\right)R-2f^{1}+3\square
f^{1}_{R}+6\square\left(f^{2}_{R}\mathcal{L}_{m}\right)$
$\displaystyle=\left(1+f^{2}\right)T,$ (4)

where $T$ is the trace of the energy-momentum tensor $T_{\mu\nu}$.

In NMC gravity the energy-momentum tensor of matter is not covariantly
conserved multiscalar ; Sotiriou1 : applying the Bianchi identities to Eq.
(3), it follows

$\nabla_{\mu}T^{\mu\nu}=\frac{f^{2}_{R}}{1+f^{2}}(g^{\mu\nu}\mathcal{L}_{m}-T^{\mu\nu})\nabla_{\mu}R,$
(5)

a property which gives rise to an extra force which is added to the fifth
force which is typical of $f(R)$ gravity theory. We will find that this extra
force has a negligible effect on the motion of Earth and Moon for values of
NMC gravity parameters of astrophysical and cosmological interest, while such
a force is expected to have important effects at the galactic scale.

###  II.1 Metric and energy-momentum tensors

We use the following notation for indices of tensors: Greek letters denote
space-time indices ranging from 0 to 3, whereas Latin letters denote spatial
indices ranging from 1 to 3. Cartesian three-vectors are indicated by boldface
type and scalar product is indicated by a dot. The signature of the metric
tensor is $(-,+,+,+)$.

The Sun is modeled as a static spherically symmetric distribution of matter,
while the Earth and Moon are modeled as orbiting spherically symmetric bodies.
The metric tensor which describes the spacetime in the Sun-Earth-Moon system
is given, in the Newtonian gauge, by

$\displaystyle ds^{2}=$ $\displaystyle-$
$\displaystyle\left[1-2\Phi(\mathbf{x},t)+2\Psi(\mathbf{x},t)\right]c^{2}dt^{2}$
(6) $\displaystyle+$
$\displaystyle\left[1+2\Phi(\mathbf{x},t)\right]\delta_{ij}dx^{i}dx^{j},$

where the potentials $\Phi$ and $\Psi$ are perturbations of the Minkowski
metric of order ${\cal O}(1/c^{2})$.

The Sun is considered as a perfect fluid in hydrostatic equilibrium, while the
Earth and Moon are approximately described as continuous bodies in a
hydrostatic state of stress, i.e., the normal stresses are equal to the
pressure and shear stresses are neglected Turcotte . In this approximation the
components of the energy-momentum tensor, to the relevant order for our
computations, for all the astronomical bodies are given by (Ref. Wi , Chapter
4.1):

$\displaystyle T^{00}$ $\displaystyle=$ $\displaystyle\rho c^{2}+{\cal
O}\left(1\right),$ (7) $\displaystyle T^{0i}$ $\displaystyle=$
$\displaystyle\rho cv^{i}+{\cal O}\left(\frac{1}{c}\right),$ (8)
$\displaystyle T^{ij}$ $\displaystyle=$ $\displaystyle\rho
v^{i}v^{j}+p\delta_{ij}+{\cal O}\left(\frac{1}{c^{2}}\right),$ (9)

where matter is characterized by density $\rho$, velocity field $v^{i}$, and
pressure $p$. The trace of the energy-momentum tensor is

$T=-\rho c^{2}+{\cal O}\left(1\right).$ (10)

In the present paper we use $\mathcal{L}_{m}=-\rho c^{2}+{\cal O}(1)$ for the
Lagrangian density of matter BLP .

###  II.2 Choice of functions $f^{1}(R)$ and $f^{2}(R)$

We choose the following functions:

$f^{1}(R)=\frac{c^{4}}{8\pi G}R,\qquad f^{2}(R)=qR^{\alpha},\quad\alpha<0,$
(11)

where the function $f^{1}(R)$ corresponds to GR and $q,\alpha$ are real
numbers that have to be considered as parameters of the NMC model of gravity.

The functions (11) have been used in Ref. drkmattgal  to model the rotation
curves of galaxies, and in Ref. curraccel  to model the current accelerated
expansion of the Universe.

##  III Approximation of the field equations

We approximate the field equations (3) and (II) taking into account that the
metric potentials $\Phi$ and $\Psi$ are small perturbations of the Minkowski
metric, so that we neglect the higher order terms that include products of
potentials or their derivatives, and cross-products between their derivatives
and the potentials. Moreover, velocities of bodies are negligible with respect
to $c$. By computing the Ricci tensor and the Ricci curvature scalar $R$ from
the metric (6), it then follows that the functions $\Phi$ and $\Psi$ satisfy
the following equations:

$\displaystyle\nabla^{2}(\Phi+\Psi)=-\frac{R}{2},$ (12)
$\displaystyle\nabla^{2}\Psi=\frac{1}{2}\left(R_{00}-\frac{R}{2}\right),$ (13)

where $\nabla^{2}$ denotes Laplace operator in flat three-dimensional space.
We introduce the scalar field $\eta$ which is a function of curvature $R$ also
explicitly depending on spacetime coordinates $(\mathbf{x},t)$ through mass
density:

$\eta=\eta(\mathbf{x},t,R)=f^{1}_{R}-2f^{2}_{R}\,\rho(\mathbf{x},t)c^{2}.$
(14)

Using the metric (6) and Eqs. (12-13), the time-time component of the field
equations (3) at leading order is given by

$\displaystyle\eta(-\nabla^{2}\Phi+\nabla^{2}\Psi)+\frac{f^{1}}{2}(1-2\Phi+2\Psi)$
(15) $\displaystyle=$
$\displaystyle\nabla_{0}\nabla_{0}\eta+(1-2\Phi+2\Psi)\square\eta+(1+f^{2})\rho
c^{2},$

and the trace (II) of the field equations becomes

$\eta R-2f^{1}+3\square\eta=(1+f^{2})T.$ (16)

Eliminating $\square\eta$ in Eq. (15) by means of the trace equation, and
$\nabla^{2}\Phi$ by means of Eq. (12), we obtain

$\displaystyle\eta\left(2\nabla^{2}\Psi+\frac{5}{6}R\right)-\frac{f^{1}}{6}(1-2\Phi+2\Psi)$
$\displaystyle=$ $\displaystyle\nabla_{0}\nabla_{0}\eta+\frac{2}{3}\eta
R(\Phi-\Psi)+\frac{2}{3}(1+f^{2})\rho c^{2}(1+\Phi-\Psi).$

We require the functions $f^{1}(R)$ and $f^{2}(R)$ to satisfy the following
conditions:

$\left|\frac{8\pi G}{c^{4}}\,\frac{f^{1}}{R}-1\right|\ll
1,\qquad\left|f^{2}\right|\ll 1,$ (18)

and the following condition on the derivatives of $f^{1}$ and $f^{2}$ with
respect to $R$,

$\left|\frac{8\pi G}{c^{4}}\eta-1\right|\ll 1.$ (19)

The conditions (18) mean that the Lagrangian density in Eq. (1) is a small
perturbation of the Lagrangian of GR. While the first of conditions (18) is
trivial for the choice (11) of function $f^{1}$, the other two conditions will
be verified a posteriori. Using such conditions Eq. (III) is approximately
given by:

$\displaystyle\nabla^{2}\Psi+\frac{R}{3}$ $\displaystyle=$
$\displaystyle\frac{8\pi
G}{3c^{2}}\rho+\frac{1}{3}\left(\frac{R}{2}+\frac{8\pi
G}{c^{2}}\rho\right)(\Phi-\Psi)$ (20) $\displaystyle+$
$\displaystyle\frac{4\pi G}{c^{4}}\nabla_{0}\nabla_{0}\eta,$

from which, keeping terms of order ${\cal O}(1/c^{2})$, we find

$\nabla^{2}\Psi=\frac{1}{3}\left(\frac{8\pi G}{c^{2}}\rho-R\right),$ (21)

and, using Eq. (12),

$\nabla^{2}\Phi=-\frac{4\pi G}{c^{2}}\rho+\frac{1}{6}\left(\frac{8\pi
G}{c^{2}}\rho-R\right).$ (22)

###  III.1 Equation for the scalar field $\eta$

Equations (21) and (22) have to be completed with an equation for the scalar
field $\eta$. Neglecting cross-products between the potentials $\Phi,\Psi$ and
their derivatives we have at leading order

$\displaystyle\square\eta$ $\displaystyle=$
$\displaystyle(1-2\Phi)\nabla^{2}\eta-\frac{1}{c^{2}}(1+2\Phi-2\Psi)\frac{\partial^{2}\eta}{\partial
t^{2}}+\frac{\partial\Psi}{\partial x_{i}}\frac{\partial\eta}{\partial x^{i}}$
(23) $\displaystyle\approx$ $\displaystyle\nabla^{2}\eta.$

Using conditions (18-19) the trace equation (16) is approximately given by

$\nabla^{2}\eta=\frac{c^{4}}{24\pi G}\,R-\frac{1}{3}\rho c^{2}.$ (24)

By introducing a potential function $V=V(\mathbf{x},t,\eta)$ and an effective
potential $V_{\rm eff}$ as in Refs. KW ; HS ,

$\frac{\partial V}{\partial\eta}=\frac{c^{4}}{24\pi
G}\,\omega(\eta,\rho),\qquad V_{\rm eff}=V-\frac{1}{3}\rho c^{2}\eta,$ (25)

where the function $\omega(\eta,\rho)$ is obtained by solving the equation
(14) with respect to $R$, the equation for the scalar field $\eta$ becomes

$\nabla^{2}\eta=\frac{\partial V_{\rm eff}}{\partial\eta}.$ (26)

Note that for the choice (11) of functions $f^{1}(R),f^{2}(R)$ the function
$\omega(\eta,\rho)$ exists and it is unique. The effective potential has an
extremum which corresponds to the GR solution

$R=\omega(\eta,\rho)=\frac{8\pi G}{c^{2}}\rho,$ (27)

and we require that such an extremum is a minimum KW ; HS , which yields the
condition

$\frac{\partial^{2}V_{\rm eff}}{\partial\eta^{2}}=\frac{c^{4}}{24\pi
G}\,\frac{1}{\eta_{R}}\geq 0,$ (28)

with $\eta_{R}=f^{1}_{RR}-2f^{2}_{RR}\rho c^{2}$ and $R=\omega(\eta,\rho)$,
the double subscript in $f^{i}_{RR}$ denoting second derivative with respect
to $R$. For the choice (11) of $f^{1},f^{2}$ such a minimum condition requires

$\alpha(\alpha-1)q[\omega(\eta,\rho)]^{\alpha-2}\leq 0,$ (29)

from which, for $q\neq 0$ and $R=\omega(\eta,\rho)>0$, it follows

$\alpha<0\implies q<0,$ (30)

which is an application of a general stability condition against Dolgov-
Kawasaki instability in NMC gravity found in Refs. Faraoni ; BeSeq .

At the minimum of the effective potential $V_{\rm eff}$ we set

$\frac{\partial^{2}V_{\rm
eff}}{\partial\eta^{2}}=\frac{1}{\lambda^{2}}>0,\quad\mbox{with}\quad\omega(\eta,\rho)=\frac{8\pi
G}{c^{2}}\rho,$ (31)

where $\lambda=\lambda(\rho)>0$ has dimension of length and depends on mass
density. For the choice (11) of functions of curvature we have

$\lambda^{2}=6q\alpha(1-\alpha)\left(\frac{8\pi
G}{c^{2}}\rho\right)^{\alpha-1},$ (32)

and the function $\lambda(\rho)$ decreases as density $\rho$ increases. In the
next sections we consider only the choice (11) of functions
$f^{1}(R),f^{2}(R)$, and we compute an analytic approximate solution of the
equation for the scalar field $\eta$.

##  IV Solution in the interior of bodies and in interplanetary space

Equation (26) for the scalar field $\eta$ is typical of chameleon theories of
gravity KW ; HS . The difference with respect to other chameleon theories such
as $f(R)$ gravity consists in the explicit dependence of $\partial
V/\partial\eta$ on $\rho(\mathbf{x},t)$ due to the nonminimal coupling (see
the discussion in Ref. MBMGDeA ).

In order to satisfy the stringent bounds from Solar System experiments on
modified gravity, a chameleon theory requires the solution $\eta$ to remain
close to the minimizer of the effective potential $V_{\rm eff}$ in most of the
interior of massive astronomical bodies such as the Sun, Earth and Moon, so
that GR is approximately satisfied KW ; HS . More precisely, in each body
$\eta$ has to be close to the minimizer of $V_{\rm eff}$ inside a critical
radius, called the screening radius, that has to be determined. If the
screening radius is close to the radius of the astronomical body for each
body, then the thin shell condition is satisfied and deviations from GR are
screened. In the following we denote by $r_{S},r_{E}$ and $r_{M}$ the
screening radii of the Sun, Earth and Moon, respectively.

For each body we compute a solution inside the screening radius by using the
approximation of spherical symmetry around the center. In Sec. V.5 we will
match the interior solution with the solution outside of the screening radii
and we will show that the approximation of interior spherical symmetry can be
used, provided that the distances between the astronomical bodies are much
greater than the radii, a condition which is satisfied for the Sun-Earth-Moon
system.

For a spherically symmetric solution the finiteness of $\nabla^{2}\eta$
imposes the boundary condition

$\frac{d\eta}{dr}=0,\quad\mbox{at }r=0,$ (33)

at the center of each astronomical body, where the variable $r$ is distance
from the center.

###  IV.1 Solution in the Sun’s interior

The Sun is modeled as a static spherically symmetric distribution of matter
with density $\rho_{S}=\rho_{S}(r)$ where $r$ is distance from the center. A
model of mass density profile for the Sun has been used in Ref. MBMGDeA  and
the parts of the model that will be used in the sequel are reported in Section
.1 of Appendix A.

Since the effective potential has an extremum which corresponds to the GR
solution, $\partial V_{\rm eff}/\partial\eta=0$, then expression (27) of
curvature yields an exact solution of the equation (24) only if
$\nabla^{2}\eta=0$. Under spherical symmetry, the only harmonic function which
satisfies the boundary condition (33) is a constant, however, one can check
that the solution $\eta=constant$ implies that density $\rho_{S}(r)$ must also
be constant, which is not the case for the Sun’s interior MBMGDeA .

Though the GR solution is not an exact solution of Eq. (26), it is an an
approximate solution if the following consistency condition is satisfied for
$r<r_{S}$ HS ; MBMGDeA :

$\left|\nabla^{2}\left(\eta(r,R=8\pi
G\rho_{S}(r)/c^{2})\right)\right|\ll\frac{1}{3}\rho_{S}c^{2},$ (34)

with $r_{S}<R_{\odot}$, where $R_{\odot}$ is the Sun’s radius. Computing the
Laplacian of $\eta$ according to (14), the consistency condition has been
evaluated in Ref. MBMGDeA  and reads

$\lambda^{2}\left|\frac{\alpha}{1-\alpha}\nabla^{2}\rho_{S}-\frac{\alpha+1}{\rho_{S}}\left(\frac{d\rho_{S}}{dr}\right)^{2}\right|\ll\rho_{S}.$
(35)

Then, using definition (14) of $\eta$, the GR expression of curvature (27),
and formula (32) of $\lambda^{2}$, we have the following approximate solution:

$\displaystyle\eta$ $\displaystyle\approx$ $\displaystyle\frac{c^{4}}{8\pi
G}-2\alpha q\left(\frac{8\pi G}{c^{2}}\rho_{S}\right)^{\alpha-1}\rho_{S}c^{2}$
(36) $\displaystyle=$ $\displaystyle\frac{c^{4}}{8\pi
G}-\frac{\lambda^{2}(\rho_{S})}{3(1-\alpha)}\rho_{S}c^{2},$

which holds for $r<r_{S}$. The boundary condition (33) is satisfied provided
that the Sun’s density model has the property $d\rho_{S}/dr=0$ at the center.

Eventually, by using the Sun density profile in Section .1 of Appendix A, we
find $\lambda(\rho_{S})\ll R_{\odot}$, particularly $\lambda(\rho_{S})$ turns
out to be a negligible quantity for $|\alpha|$ not too close to zero.

###  IV.2 Solution in the Earth’s interior

Differently from the Sun, the density profile of the Earth’s interior is
conveniently modeled by resorting to density discontinuities, detected by
seismology, such as the Mohorovičić discontinuity, or Moho, at the boundary
between the crust and the mantle. On the two sides of a density discontinuity
the values of $\eta$ that minimize $V_{\rm eff}$ are different, however, since
the function $\eta$ and its gradient have to be continuous in order to
guarantee the existence of $\nabla^{2}\eta$, then at the discontinuity the
solution $\eta$ has to interpolate between the two minimizers so that it
cannot be close to both minimizers. Hence, the solution locally deviates from
GR and one has to check if such a deviation is small enough in such a way that
screening takes place anyway.

The Earth is modeled as a spherically symmetric distribution of matter with
density $\rho_{E}=\rho_{E}(r)$, where $r$ is distance from the center, and
axial rotation of Earth is neglected. An average Earth model PREM  is
considered, and the planet is divided into four homogeneous regions separated
by spherical surfaces of density discontinuities: ocean layer, crust, mantle
and core. Numerical values of density and radii of discontinuity surfaces are
reported in Section .2 of Appendix A.

Since, for $r<r_{E}$, the solution $\eta$ has to remain close to the minimizer
of $V_{\rm eff}$, then inside each region of constant density the derivative
of the potential $V$ is approximated by

$\frac{\partial V}{\partial\eta}(\eta,\rho_{E})\approx\frac{\partial
V}{\partial\eta}(\eta_{\rm
min},\rho_{E})+\frac{\partial^{2}V}{\partial\eta^{2}}(\eta_{\rm
min},\rho_{E})(\eta-\eta_{\rm min}),$ (37)

where $\eta_{\rm min}$ is the minimizer of $V_{\rm eff}$ in the considered
region. Since

$\frac{\partial V}{\partial\eta}(\eta_{\rm
min},\rho_{E})=\frac{1}{3}\rho_{E}c^{2},$ (38)

Eq. (26) for $\eta$ becomes

$\nabla^{2}\eta\approx\frac{1}{\lambda^{2}(\rho_{E})}(\eta-\eta_{\rm min}).$
(39)

Since in each region the value of density $\rho_{E}$ is constant, using (32)
$\lambda^{2}(\rho_{E})$ is also constant, then Eq. (39) admits a closed-form
solution that, for a given screening radius $r_{E}$ of Earth, is completely
determined by conditions of continuity of $\eta$ and its radial derivative at
the discontinuity surfaces, and by the following condition at $r_{E}$ given in
Ref. TaTsu :

$\frac{1}{\lambda^{2}(\rho_{E})}\left[\eta(r_{E})-\eta_{\rm
min}\right]=-\frac{1}{3}\rho_{E}c^{2},$ (40)

where $r_{E}<R_{\oplus}$, with $R_{\oplus}$ the Earth’s radius, and $\rho_{E}$
is evaluated at the uppermost layer of the Earth’s interior model, i.e., the
ocean layer.

Note that the Sun’s interior is modeled by means of a continuously varying
density profile (see Section .1 of Appendix A), hence Eq. (39) in the case of
Sun does not admit a closed-form solution, so that we have used a different
approximation of the solution $\eta$ given by Eq. (36). Conversely, the
presence of discontinuity surfaces in the Earth’s interior prevents us from
using a consistency condition of the type (35) where derivatives of density
are involved, so that we resorted to a piecewise constant density profile
which permits us to compute an analytic solution that can be proved to be
close to the minimizer of the effective potential, except at the density
discontinuities.

In the sequel $R_{E,c}$ is the outer radius of the crust, $R_{E,m}$ is the
outher radius of the mantle, and $R_{E,n}$ is the radius of the core
(nucleus). The mass densities of ocean, crust, mantle and core are denoted by
$\rho_{E,w},\rho_{E,c},\rho_{E,m}$ and $\rho_{E,n}$, respectively. The
corresponding values of $\lambda(\rho_{E})$ are denoted by
$\lambda_{E,w},\lambda_{E,c},\lambda_{E,m}$ and $\lambda_{E,n}$. Analogously,
the values of $\eta$ minimizing the effective potential in the various layers
are denoted by $\eta_{E,w},\eta_{E,c},\eta_{E,m}$ and $\eta_{E,n}$, and they
are given by

$\eta_{\rm min}=\frac{c^{4}}{8\pi
G}-\frac{\lambda^{2}(\rho_{E})}{3(1-\alpha)}\rho_{E}c^{2}.$ (41)

By using the numerical values of density of the various Earth’s layers, we
find in Section .2 of Appendix A $\lambda(\rho_{E})\ll R_{\oplus}$,
particularly $\lambda(\rho_{E})$ turns out to be a completely negligible
quantity for $|\alpha|$ not too close to zero.

The expression of the analytic solution of Eq. (39) in the various Earth’s
layers is cumbersome, nevertheless, it admits a manageable approximation, that
guarantees continuity of $\eta$ and approximate continuity of its derivative,
and it is given in the following.

Ocean layer. The ocean and seas cover $70.8\%$ of the surface of the Earth, so
that we approximate the uppermost layer with seawater. The approximate
solution for $R_{E,c}\leq r\leq r_{E}$ is given by

$\displaystyle\eta(r)\approx\eta_{E,w}-\frac{1}{r}\left[\frac{1}{3}r_{E}\lambda^{2}_{E,w}\rho_{E,w}c^{2}\exp\left(\frac{r-r_{E}}{\lambda_{E,w}}\right)\right.$
$\displaystyle-$
$\displaystyle\left.\frac{\lambda_{E,w}/\lambda_{E,c}}{1+\lambda_{E,w}/\lambda_{E,c}}R_{E,c}(\eta_{E,c}-\eta_{E,w})\exp\left(\frac{R_{E,c}-r}{\lambda_{E,w}}\right)\right].$

If $r_{E}$ is close to $R_{\oplus}$ then the thin shell condition for Earth is
satisfied, and $r_{E}-R_{E,c}\gg\lambda_{E,w}$. It turns out that the
difference $|\eta(r)-\eta_{E,w}|$ between $\eta$ and the minimizer is
exponentially suppressed for $r_{E}-r\gg\lambda_{E,w}$ and
$r-R_{E,c}\gg\lambda_{E,w}$, hence in most of the ocean layer, due to the
smallness of $\lambda_{E,w}$. The value of $\eta$ at the boundary between
seawater and the oceanic crust is

$\eta(R_{E,c})\approx\eta_{E,w}+\frac{\lambda_{E,w}/\lambda_{E,c}}{1+\lambda_{E,w}/\lambda_{E,c}}(\eta_{E,c}-\eta_{E,w}),$
(43)

and we see that the solution $\eta$ interpolates between the minimizers
$\eta_{E,w}$ and $\eta_{E,c}$, hence deviating from GR in a thin shell of
thickness of order $\lambda_{E,w}$. In Secs. VI.1 and VI.3 we will prove that
the resulting perturbation of the Newtonian gravitational force, hence of
hydrostatic equilibrium, is negligible.

Crust. The approximate solution for $R_{E,m}\leq r\leq R_{E,c}$ is given by

$\displaystyle\eta(r)\approx\eta_{E,c}-\frac{1}{r}\left[\frac{R_{E,c}(\eta_{E,c}-\eta_{E,w})}{1+\lambda_{E,w}/\lambda_{E,c}}\exp\left(\frac{r-R_{E,c}}{\lambda_{E,c}}\right)\right.$
$\displaystyle-$
$\displaystyle\left.\frac{\lambda_{E,c}/\lambda_{E,m}}{1+\lambda_{E,c}/\lambda_{E,m}}R_{E,m}(\eta_{E,m}-\eta_{E,c})\exp\left(\frac{R_{E,m}-r}{\lambda_{E,c}}\right)\right].$

Again, it turns out that the difference between $\eta$ and the minimizer
$\eta_{E,c}$ is exponentially suppressed in most of the crust. The value of
$\eta$ at Moho, the discontinuity between the crust and the mantle, is

$\eta(R_{E,m})\approx\eta_{E,c}+\frac{\lambda_{E,c}/\lambda_{E,m}}{1+\lambda_{E,c}/\lambda_{E,m}}(\eta_{E,m}-\eta_{E,c}),$
(45)

a formula analogous to Eq. (43). In the crust the solution $\eta$ deviates
from GR in thin shells of thickness of order $\lambda_{E,c}$ adjacent to the
upper boundary at $r=R_{E,c}$ and to the lower boundary at $r=R_{E,m}$,
respectively. The resulting perturbation will turn out to be again negligible,
so that screening takes place.

Mantle. The approximate solution for $R_{E,n}\leq r\leq R_{E,m}$ is given by

$\displaystyle\eta(r)\approx\eta_{E,m}-\frac{1}{r}\left[\frac{R_{E,m}(\eta_{E,m}-\eta_{E,c})}{1+\lambda_{E,c}/\lambda_{E,m}}\exp\left(\frac{r-R_{E,m}}{\lambda_{E,m}}\right)\right.$
$\displaystyle-\left.\frac{\lambda_{E,m}/\lambda_{E,n}}{1+\lambda_{E,m}/\lambda_{E,n}}R_{E,n}(\eta_{E,n}-\eta_{E,m})\exp\left(\frac{R_{E,n}-r}{\lambda_{E,m}}\right)\right].$

The properties of the solution in the mantle are analogous to the properties
in the crust.

Core. The approximate solution for $0<r\leq R_{E,n}$ is given by

$\displaystyle\eta(r)$ $\displaystyle\approx$
$\displaystyle\eta_{E,n}-2\frac{R_{E,n}}{r}\frac{(\eta_{E,n}-\eta_{E,m})}{1+\lambda_{E,m}/\lambda_{E,n}}\exp\left(-\frac{R_{E,n}}{\lambda_{E,n}}\right)$
(47) $\displaystyle\times$
$\displaystyle\sinh\left(\frac{r}{\lambda_{E,n}}\right).$

The difference $|\eta(r)-\eta_{E,n}|$ between $\eta$ and the minimizer is
exponentially suppressed for $R_{E,n}-r\gg\lambda_{E,n}$, hence in the whole
core except in a thin shell of thickness of order $\lambda_{E,n}$ adjacent to
the boundary of the core. The approximate solution in the core satisfies the
boundary condition (33).

Using expression (32) of $\lambda^{2}$ we see that, for given values of NMC
gravity parameters $\alpha$ and $q$, the approximate solution in the Earth’s
interior is completely determined in all the layers except the ocean layer
where it depends on the screening radius $r_{E}$, hence it is determined
everywhere once the screening radius is determined.

###  IV.3 Solution in the Moon’s interior

The model of the lunar interior is analogous to the Earth’s model. The Moon is
modeled as a spherically symmetric distribution of matter with density
$\rho_{M}=\rho_{M}(r)$, where $r$ is distance from the center, and the
satellite is divided into three homogeneous regions separated by spherical
surfaces of density discontinuities: crust, mantle and core. Numerical values
of density and radii of discontinuity surfaces are reported in Section .3 of
Appendix A. The length $\lambda(\rho_{M})$ turns out to be again a completely
negligible quantity.

In the sequel $R_{M}$ is the Moon’s radius, $R_{M,m}$ is the outher radius of
the mantle, and $R_{M,n}$ is the radius of the core. The mass densities of
crust, mantle and core are denoted by $\rho_{M,c},\rho_{M,m}$ and
$\rho_{M,n}$, respectively. The corresponding values of $\lambda(\rho_{M})$
are denoted by $\lambda_{M,c},\lambda_{M,m}$ and $\lambda_{M,n}$. Analogously,
the values of $\eta$ minimizing the effective potential in the various layers
are denoted by $\eta_{M,c},\eta_{M,m}$ and $\eta_{M,n}$.

The approximate solution for $\eta$ is analogous to the one found for Earth.

Crust. The approximate solution for $R_{M,m}\leq r\leq r_{M}$ is given by

$\displaystyle\eta(r)\approx\eta_{M,c}-\frac{1}{r}\left[\frac{1}{3}r_{M}\lambda^{2}_{M,c}\rho_{M,c}c^{2}\exp\left(\frac{r-r_{M}}{\lambda_{M,c}}\right)\right.$
$\displaystyle-$
$\displaystyle\left.\frac{\lambda_{M,c}/\lambda_{M,m}}{1+\lambda_{M,c}/\lambda_{M,m}}R_{M,m}(\eta_{M,m}-\eta_{M,c})\exp\left(\frac{R_{M,m}-r}{\lambda_{M,c}}\right)\right].$

The properties of the solution are analogous to the ones in the ocean layer of
Earth. The value of $\eta$ at the lunar Moho is

$\eta(R_{M,m})\approx\eta_{M,c}+\frac{\lambda_{M,c}/\lambda_{M,m}}{1+\lambda_{M,c}/\lambda_{M,m}}(\eta_{M,m}-\eta_{M,c}).$
(49)

Mantle. The approximate solution for $R_{M,n}\leq r\leq R_{M,m}$ is given by

$\displaystyle\eta(r)\approx\eta_{M,m}-\frac{1}{r}\left[\frac{R_{M,m}(\eta_{M,m}-\eta_{M,c})}{1+\lambda_{M,c}/\lambda_{M,m}}\exp\left(\frac{r-R_{M,m}}{\lambda_{M,m}}\right)\right.$
$\displaystyle-\left.\frac{\lambda_{M,m}/\lambda_{M,n}}{1+\lambda_{M,m}/\lambda_{M,n}}R_{M,n}(\eta_{M,n}-\eta_{M,m})\exp\left(\frac{R_{M,n}-r}{\lambda_{M,m}}\right)\right].$

Core. The approximate solution for $0<r\leq R_{M,n}$ is given by

$\displaystyle\eta(r)$ $\displaystyle\approx$
$\displaystyle\eta_{M,n}-2\frac{R_{M,n}}{r}\frac{(\eta_{M,n}-\eta_{M,m})}{1+\lambda_{M,m}/\lambda_{M,n}}\exp\left(-\frac{R_{M,n}}{\lambda_{M,n}}\right)$
(51) $\displaystyle\times$
$\displaystyle\sinh\left(\frac{r}{\lambda_{M,n}}\right).$

The approximate solution in the core satisfies the boundary condition (33).
Again, the approximate solution in the lunar interior is completely determined
once the Moon’s screening radius $r_{M}$ is determined.

###  IV.4 Solution in the outskirts of the Solar System

We assume that in the solar neighborhood of the Galaxy the field $\eta$ is
close to the minimizer of the effective potential $V_{\rm eff}$, so that the
spacetime curvature $R$ is approximately given by the GR solution MBMGDeA .
This assumption implies that the Milky Way is screened within a distance of
about $8\,{\rm kpc}$ from its center, where the Solar System is approximately
located. Such a screening condition may impose additional constraints on the
NMC gravity model whose assessment requires the solution for the gravitational
field of the Milky Way, possibly taking also into account the effect of the
other galaxies in the local group. That will be the subject of a future paper.

The galactic mass density $\rho_{g}$ in the solar neighborhood of the Milky
Way is $\rho_{g}\approx 6.9\times 10^{-24}\,{\rm g}/{\rm cm}^{3}$ HoFly , so
that we have $R\approx R_{g}$, with $R_{g}$ the GR solution $R_{g}=8\pi
G\rho_{g}/c^{2}$. The minimizer $\eta_{g}$ of the effective potential $V_{\rm
eff}(\eta,\rho)$, corresponding to $\rho(\mathbf{x},t)=\rho_{g}$, is given by

$\eta_{g}=\frac{c^{4}}{8\pi G}-2\alpha q\left(\frac{8\pi
G}{c^{2}}\right)^{\alpha-1}\rho_{g}^{\alpha}c^{2}.$ (52)

We denote by $r_{g}$ a distance from the Sun’s center such that mass density
is dominated by the galactic density component at points $\mathbf{x}$ such
that $|\mathbf{x}-\mathbf{x}_{S}|>r_{g}$, where $\mathbf{x}_{S}$ is the
position vector of Sun’s center. We choose $r_{g}$ at the heliopause, the
boundary between the solar wind and the interstellar medium MBMGDeA ,
corresponding to a heliocentric radial distance of about $120\mbox{ AU
}=2.58\times 10^{4}R_{\odot}$, where $R_{\odot}$ is the Sun’s radius.

Since $\eta$ approximately minimizes $V_{\rm eff}$ in the solar neighborhood
of the Galaxy, then Eq. (26) becomes

$\nabla^{2}\eta\approx\frac{1}{\lambda_{g}^{2}}(\eta-\eta_{g}),\qquad\mbox{for
}|\mathbf{x}-\mathbf{x}_{S}|>r_{g},$ (53)

where $\lambda_{g}=\lambda(\rho_{g})$. The length $\lambda(\rho)$ increases as
density decreases, which is a typical property of the chameleon mechanism KW ,
so that the length $\lambda_{g}$ is an upper bound for $\lambda$ in the Solar
System. The computations in the present paper will be made under the condition
$\lambda_{g}\gg r_{g}$ which will permit us to find analytic estimates of the
results. We will find that such a condition is satisfied when the constraint
from LLR measurements is saturated. Then we have

$\nabla^{2}\eta=0\qquad\mbox{for
}r_{g}<|\mathbf{x}-\mathbf{x}_{S}|\ll\lambda_{g},$ (54)

and, for $|\mathbf{x}-\mathbf{x}_{S}|$ large enough, we assume

$\eta(\mathbf{x},t)\approx\eta_{g}.$ (55)

Under our assumptions the solution $\eta$ is a harmonic function in the
outskirts of the Solar System.

###  IV.5 Solution in interplanetary space

In interplanetary space, where mass density $\rho$ is small and gradually
approaches the galactic density $\rho_{g}$ as $|\mathbf{x}-\mathbf{x}_{S}|$
approaches $r_{g}$, we proceed as in Ref. MBMGDeA  and we expand the
derivative of the potential $V$ around the minimizer $\eta_{g}$:

$\frac{\partial V}{\partial\eta}(\eta,\rho)\approx\frac{\partial
V}{\partial\eta}(\eta_{g},\rho)+\frac{\partial^{2}V}{\partial\eta^{2}}(\eta_{g},\rho)(\eta-\eta_{g}).$
(56)

Solving the expression (14) of $\eta$,

$\eta=\frac{c^{4}}{8\pi G}-2\alpha qR^{\alpha-1}\rho c^{2},$ (57)

with respect to curvature $R=\omega(\eta,\rho)$ we find

$\omega(\eta,\rho)=\left(\frac{16\pi G}{c^{2}}\alpha
q\rho\right)^{1/(1-\alpha)}\left(1-\frac{8\pi
G}{c^{4}}\eta\right)^{1/(\alpha-1)},$ (58)

from which, using the first equality in Eq. (25), we obtain the property

$\frac{\partial
V}{\partial\eta}(\eta,\rho)=\left(\frac{\rho}{\rho_{g}}\right)^{1/(1-\alpha)}\frac{\partial
V}{\partial\eta}(\eta,\rho_{g}).$ (59)

Taking now into account that at density $\rho_{g}$ (in the solar vicinity of
the Galaxy) the field $\eta$ approximately minimizes the effective potential
$V_{\rm eff}$, so that

$\frac{\partial
V}{\partial\eta}(\eta_{g},\rho_{g})\approx\frac{1}{3}\,\rho_{g}c^{2},$ (60)

we can compute the approximation (56) of the derivative of the potential:

$\displaystyle\frac{\partial V}{\partial\eta}(\eta,\rho)$
$\displaystyle\approx$
$\displaystyle\left(\frac{\rho}{\rho_{g}}\right)^{1/(1-\alpha)}\left[\frac{\partial
V}{\partial\eta}(\eta_{g},\rho_{g})+\frac{\partial^{2}V}{\partial\eta^{2}}(\eta_{g},\rho_{g})(\eta-\eta_{g})\right]\approx\left(\frac{\rho}{\rho_{g}}\right)^{1/(1-\alpha)}\left[\frac{1}{3}\,\rho_{g}c^{2}+\frac{1}{\lambda_{g}^{2}}(\eta-\eta_{g})\right]$
(61) $\displaystyle\approx$
$\displaystyle\frac{1}{3}\,\rho_{g}c^{2}\left(\frac{\rho}{\rho_{g}}\right)^{1/(1-\alpha)},$

where, taking into account that $\lambda_{g}$ is assumed large with respect to
$r_{g}$, we have used the inequality

$\frac{1}{\lambda_{g}^{2}}(\eta-\eta_{g})\ll\frac{1}{3}\rho_{g}c^{2},$ (62)

that will be verified a posteriori. Then we have

$\frac{\partial V_{\rm eff}}{\partial\eta}=\frac{\partial
V}{\partial\eta}-\frac{1}{3}\,\rho
c^{2}\approx-\frac{1}{3}\,\widetilde{\rho}c^{2},$ (63)

with

$\widetilde{\rho}=\rho\left[1-\left(\frac{\rho_{g}}{\rho}\right)^{-\alpha/(1-\alpha)}\right].$
(64)

Then, in interplanetary space we look for the field $\eta$ which solves the
Poisson equation

$\nabla^{2}\eta=-\frac{1}{3}\,\widetilde{\rho}c^{2}.$ (65)

As $|\mathbf{x}-\mathbf{x}_{S}|$ tends to $r_{g}$, density $\rho$ tends to
$\rho_{g}$ and the Poisson equation turns to Laplace equation (54) which holds
for $|\mathbf{x}-\mathbf{x}_{S}|\geq r_{g}$.

##  V Gravitational field of astronomical bodies

In the thin shells inside the astronomical bodies, the function $\eta$ has to
interpolate between the solution inside the screening radii computed in the
previous section and the solution in interplanetary space. By adapting to NMC
gravity the chameleon mechanism developed in Ref. KW , we require the
interpolating function $\eta$ to satisfy the condition (see also Ref. MBMGDeA
,HS )

$0<\frac{\partial V}{\partial\eta}(\eta,\rho)\ll\frac{1}{3}\,\rho c^{2},$ (66)

inside the thin shells and in the vicinity of the astronomical bodies, so that
the equation for $\eta$ becomes the Poisson equation

$\nabla^{2}\eta=-\frac{1}{3}\,\rho c^{2}.$ (67)

Inequality (66) will be verified a posteriori in Appendix E for $\alpha<-1/2$.
Since we have

$\alpha<0\implies 0<-\frac{\alpha}{1-\alpha}<1,$ (68)

for $|\alpha|$ not too small, the second term inside the square bracket in Eq.
(64) is negligible in comparison to 1 in the thin shells inside the
astronomical bodies, where mass density is much larger than the galactic
density $\rho_{g}$, so that $\widetilde{\rho}\approx\rho$ inside the thin
shells and the Poisson equation,

$\nabla^{2}\eta=-\frac{1}{3}\,\widetilde{\rho}c^{2},$ (69)

is valid with a good approximation everywhere outside of the screening radii.
In the next sections we compute an approximate solution of this equation.

###  V.1 Boundary conditions

We rewrite the Poisson equation (69) in the form

$\nabla^{2}\left(\eta-\eta_{g}\right)=-\frac{1}{3}\,\widetilde{\rho}c^{2},$
(70)

and we solve this equation in the unbounded domain $\Omega$ outside of the
screening radii, with the boundary condition (55) at large distance from the
Sun’s center. We denote by $\partial\Omega$ the boundary of $\Omega$ and by
$\mathbf{x}_{S},\mathbf{x}_{E},\mathbf{x}_{M}$ the centers of the Sun, Earth
and Moon, respectively. Then the boundary $\partial\Omega$ consists of the
union of the three spherical surfaces with centers at
$\mathbf{x}_{S},\mathbf{x}_{E},\mathbf{x}_{M}$ and radii given by
$r_{S},r_{E},r_{M}$, respectively.

We have to match the solution in $\Omega$ with the interior solution computed
in Sec. IV. Using the value (41) of $\eta$ minimizing the effective potential
inside the Earth, and the solution (IV.2) for $\eta$ in the Earth’s ocean
layer, the value of $\eta$ at the Earth’s screening surface is given by

$\eta\approx\frac{c^{4}}{8\pi
G}-\frac{2-\alpha}{3(1-\alpha)}\lambda^{2}_{E,w}\rho_{E,w}c^{2}.$ (71)

Using the expression (52) of the minimizer $\eta_{g}$ in the Galaxy, and
formula (32) for $\lambda(\rho)$, we have

$\eta_{g}=\frac{c^{4}}{8\pi
G}-\frac{\lambda^{2}_{g}}{3(1-\alpha)}\rho_{g}c^{2}.$ (72)

In the sequel we assume $|\alpha|$ not too close to zero in such a way that,
using again the formula for $\lambda(\rho)$, we have

$\frac{\lambda^{2}_{E,w}\rho_{E,w}}{\lambda^{2}_{g}\rho_{g}}=\left(\frac{\rho_{E,w}}{\rho_{g}}\right)^{\alpha}\ll
1,$ (73)

from which it follows that at the Earth’s screening surface the boundary
condition $\eta-\eta_{g}\approx c^{4}/(8\pi G)-\eta_{g}$ holds. Note that the
expression (57) of $\eta$, for $R>0$ and $\alpha<0$, implies $\eta\neq
c^{4}/(8\pi G)$, so that the above boundary condition is proposed as an
approximate condition in the sense of the inequality $|(8\pi
G/c^{4})\eta-1|\ll 1$. Analogous considerations can be applied for the
boundary conditions at the screening surfaces of the Sun and Moon.

Then we have the following Dirichlet boundary conditions for the Poisson
equation (70) in the set $\Omega$:

$\begin{cases}\eta(\mathbf{x},t)-\eta_{g}\approx\frac{c^{4}}{8\pi
G}-\eta_{g}\qquad\mbox{for $\mathbf{x}\in\partial\Omega$ and any $t$},\\\
\eta(\mathbf{x},t)-\eta_{g}\approx 0\quad\mbox{for
$r_{g}<|\mathbf{x}-\mathbf{x}_{S}|$ large and any $t$}.\end{cases}$ (74)

Moreover, the existence of $\nabla^{2}\eta$ requires the continuity of the
first partial derivatives of $\eta$ across the screening surfaces. The normal
derivative of $\eta$ at the Earth’s screening surface, computed from the
interior by using solution (IV.2), is given by

$\frac{d\eta}{dr}(r_{E})\approx-\frac{1}{3}\lambda_{E,w}\rho_{E,w}c^{2},$ (75)

and an analogous result holds for the Moon. Moreover, using the approximate
solution (36) in the Sun’s interior, the normal derivative of $\eta$ at the
Sun’s screening surface is given by

$\frac{d\eta}{dr}(r_{S})\approx-\frac{\alpha\lambda^{2}(\rho_{S}(r_{S}))}{3(1-\alpha)}\frac{d\rho_{S}}{dr}(r_{S})c^{2}.$
(76)

Then, taking into account the smallness of $\lambda(\rho)$ in the interior of
the bodies, we neglect the derivatives of $\eta$ from the interior of
$\partial\Omega$ and we impose the Neumann boundary condition:

$\nabla\eta(\mathbf{x})\cdot\hat{\mathbf{n}}(\mathbf{x})\approx
0\qquad\mbox{for $\mathbf{x}\in\partial\Omega$ and any $t$},$ (77)

where $\hat{\mathbf{n}}(x)$ denotes a normal unit vector to the surface
$\partial\Omega$ at point $\mathbf{x}$. For given screening radii
$r_{S},r_{E},r_{M}$ the Dirichlet condition uniquely determines the solution
of the Poisson equation, while the Neumann condition will be used to find the
screening radii.

###  V.2 Solution by means of Green’s function

In order to compute effects in the Sun-Earth-Moon system it is enough to
compute the solution $\eta(\mathbf{x},t)$ for $|\mathbf{x}-\mathbf{x}_{S}|$ of
the order of Earth’s distance from the Sun. For a given screening radii and
time instant $t$, the solution of the Dirichlet problem for the Poisson
equation is represented by means of the Green’s function
$G(\mathbf{x},\mathbf{x}_{0})$ Haberm :

$\displaystyle\eta(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\eta_{g}-\frac{c^{2}}{3}\int_{\Omega}\widetilde{\rho}(\mathbf{x}_{0},t)G(\mathbf{x},\mathbf{x}_{0})d^{3}x_{0}$
$\displaystyle+$ $\displaystyle\left(\frac{c^{4}}{8\pi
G}-\eta_{g}\right)\int_{\partial\Omega}\nabla_{\mathbf{x}_{0}}G(\mathbf{x},\mathbf{x}_{0})\cdot\hat{\mathbf{n}}(\mathbf{x}_{0})dS_{0},$

provided that, for $|\mathbf{x}-\mathbf{x}_{S}|$ of the order of Earth’s
distance from the Sun, the following inequality is satisfied

$|\eta(\mathbf{x},t)|\gg|I(\eta,G,R^{\ast})|,$ (79)

where $I(\eta,G,R^{\ast})$ is the following integral evaluated on a sphere
$S(R^{\ast})$ of large enough radius $R^{\ast}>r_{g}$ and center in
$\mathbf{x}_{S}$,

$I(\eta,G,R^{\ast})=\int_{S(R^{\ast})}\left[\left(\eta-\eta_{g}\right)\nabla
G-G\nabla\eta\right]\cdot\hat{\mathbf{n}}\,dS.$ (80)

In Green’s representation formula (V.2) the normal unit vector
$\hat{\mathbf{n}}$ to $\partial\Omega$ points towards the interior of bodies.
In the next section we compute an analytical approximation of the Green’s
function.

###  V.3 Method of images for a system of spheres

If the set $\partial\Omega$ would consist of a single sphere, then the Green’s
function could be obtained by using the method of images. Since in our case
the set $\partial\Omega$ consists of three spheres, which are the three
screening surfaces of Sun, Earth and Moon, respectively, then we apply to our
problem the extension of the method of images to a system of spheres that has
been proposed in Ref. MetzLa . This is an iterative method that involves an
infinite series of images so that it yields a representation of Green’s
function by means of an infinite series:

$G(\mathbf{x},\mathbf{x}_{0})=G^{(0)}(\mathbf{x},\mathbf{x}_{0})+G^{(1)}(\mathbf{x},\mathbf{x}_{0})+G^{(2)}(\mathbf{x},\mathbf{x}_{0})+\cdots$
(81)

The first terms of the series are obtained as follows (see MetzLa  for further
details). The zeroth-order term is

$G^{(0)}(\mathbf{x},\mathbf{x}_{0})=-\frac{1}{4\pi}\,\frac{1}{|\mathbf{x}-\mathbf{x}_{0}|}\,,$
(82)

with $\mathbf{x}_{0}\in\Omega$, so that $\mathbf{x}_{0}$ lies outside of the
three screening spheres.

The first order term $G^{(1)}(\mathbf{x},\mathbf{x}_{0})$ involves three image
points: one image inside each screening sphere. In each sphere the image point
is obtained by applying the usual method of electrostatics to the unitary
source at $\mathbf{x}_{0}$:

$\displaystyle
G^{(1)}(\mathbf{x},\mathbf{x}_{0})=\frac{1}{4\pi}\left(\frac{r_{S}}{|\mathbf{x}_{S}-\mathbf{x}_{0}|\cdot|\mathbf{x}-\widetilde{\mathbf{x}}_{S}|}\right.$
$\displaystyle+$
$\displaystyle\left.\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{0}|\cdot|\mathbf{x}-\widetilde{\mathbf{x}}_{E}|}+\frac{r_{M}}{|\mathbf{x}_{M}-\mathbf{x}_{0}|\cdot|\mathbf{x}-\widetilde{\mathbf{x}}_{M}|}\right),$

where $\widetilde{\mathbf{x}}_{S}$ is the image of $\mathbf{x}_{0}$ inside the
screening sphere of the Sun, which is given by

$\widetilde{\mathbf{x}}_{S}=\mathbf{x}_{S}+r_{S}^{2}\,\frac{\mathbf{x}_{0}-\mathbf{x}_{S}}{|\mathbf{x}_{0}-\mathbf{x}_{S}|^{2}}\,,$
(84)

and the image points $\widetilde{\mathbf{x}}_{E}$ and
$\widetilde{\mathbf{x}}_{M}$, inside the screening spheres of Earth and Moon,
are obtained by replacing in the expression of $\widetilde{\mathbf{x}}_{S}$
the subscript $S$ with $E$ and $M$, respectively.

The second order term $G^{(2)}(\mathbf{x},\mathbf{x}_{0})$ involves six image
points: two images inside each screening sphere. The image points in each
sphere are obtained by iterating the procedure used for the first order term:
consider the three sources at points
$\widetilde{\mathbf{x}}_{S},\widetilde{\mathbf{x}}_{E},\widetilde{\mathbf{x}}_{M}$
with charges

$-\frac{r_{S}}{|\mathbf{x}_{S}-\mathbf{x}_{0}|},\quad-\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{0}|},\quad-\frac{r_{M}}{|\mathbf{x}_{M}-\mathbf{x}_{0}|},$
(85)

respectively. Then, the six image points are obtained by applying the method
of electrostatics to the above sources and they are given by:

  * (i)

the images inside the Sun’s screening sphere of the sources at
$\widetilde{\mathbf{x}}_{E}$ and $\widetilde{\mathbf{x}}_{M}$;

  * (ii)

the images inside the Earth’s screening sphere of the sources at
$\widetilde{\mathbf{x}}_{S}$ and $\widetilde{\mathbf{x}}_{M}$;

  * (iii)

the images inside the Moon’s screening sphere of the sources at
$\widetilde{\mathbf{x}}_{S}$ and $\widetilde{\mathbf{x}}_{E}$.

The resulting expression of $G^{(2)}(\mathbf{x},\mathbf{x}_{0})$ is given in
Appendix B. This procedure is repeated iteratively giving rise to an infinite
series of images and terms in the Green’s function. The convergence of the
series is discussed in Ref. MetzLa . If the distances between the spheres are
much greater than the radii, a condition which is satisfied for the Sun-Earth-
Moon system, then the contribution of the higher order terms decreases
quickly, and we find that the Green’s function up to the second order will
suffice while higher order terms will not be necessary. Eventually, we observe
that notwithstanding the ratio of Sun’s radius and the Earth-Moon distance is
not small, such a ratio never appears in the computations.

###  V.4 Solution for the scalar field $\eta$

In this section we give the solution $\eta$ in $\Omega$ by using the Green’s
function up to the first order, while the second order terms are given in
Appendix B. Then the function
$G^{(0)}(\mathbf{x},\mathbf{x}_{0})+G^{(1)}(\mathbf{x},\mathbf{x}_{0})$ is
substituted in Green’s representation formula (V.2), and the volume integral
over $\Omega$ and the surface integral over $\partial\Omega$ are computed
according to the following scheme:

  * (i)

we have $\widetilde{\rho}\approx\rho$ inside the thin shells, and density
becomes immediately much smaller outside the thin shells in solar atmosphere,
in terrestrial atmosphere and outside of the Moon; then in the volume integral
over $\Omega$ the contribution outside of the thin shells can be safely
neglected, so that we have to add the integrals over the three shells; these
integrals are then evaluated in closed-form by resorting to spherical
coordinates for each of the three astronomical bodies;

  * (ii)

the surface integrals over the three screening spheres that constitute
$\partial\Omega$ are evaluated by using Gauss theorem.

The result of the computation is the following. For each astronomical body we
introduce the effective mass which is a function of the screening radius. The
effective mass of the Sun is

$M_{\odot,{\rm eff}}(r_{S})=4\pi\int_{r_{S}}^{R_{\odot}}\rho_{S}(r)r^{2}dr\,,$
(86)

and the effective masses $M_{\oplus,{\rm eff}}(r_{E})$ and $M_{M,{\rm
eff}}(r_{M})$ of Earth and Moon are defined analogously. The solution $\eta$
is given by

$\eta(\mathbf{x},t)=\eta_{S}(\mathbf{x},t)+\eta_{E}(\mathbf{x},t)+\eta_{M}(\mathbf{x},t)+\eta_{g},$
(87)

where $\eta_{S},\eta_{E}$ and $\eta_{M}$ are the contributions from the thin
shells of Sun, Earth and Moon, respectively. The contribution from the Earth’s
shell is

$\displaystyle\eta_{E}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\left(\frac{c^{4}}{8\pi
G}-\eta_{g}\right)\frac{r_{E}}{|\mathbf{x}-\mathbf{x}_{E}|}+\mathcal{I}_{E}(\mathbf{x},t)+\mathcal{J}_{E}(\mathbf{x},t)$
(88) $\displaystyle-$
$\displaystyle\frac{c^{2}}{3}\,\frac{r_{E}}{|\mathbf{x}-\mathbf{x}_{E}|}\int_{r_{E}}^{R_{\oplus}}\rho_{E}(r)r\,dr\,,$

where the term $\mathcal{I}_{E}$ is given by

$\displaystyle\mathcal{I}_{E}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\frac{c^{2}}{3|\mathbf{x}-\mathbf{x}_{E}|}\int_{r_{E}}^{|\mathbf{x}-\mathbf{x}_{E}|}\rho_{E}(r)r^{2}dr$
$\displaystyle+$
$\displaystyle\frac{c^{2}}{3}\int_{|\mathbf{x}-\mathbf{x}_{E}|}^{R_{\oplus}}\rho_{E}(r)r\,dr\qquad\mbox{for
}|\mathbf{x}-\mathbf{x}_{E}|<R_{\oplus},$

and

$\mathcal{I}_{E}(\mathbf{x},t)=\frac{c^{2}}{12\pi}\,\frac{M_{\oplus,{\rm
eff}}(r_{E})}{|\mathbf{x}-\mathbf{x}_{E}|}\qquad\mbox{for
}|\mathbf{x}-\mathbf{x}_{E}|\geq R_{\oplus},$ (90)

and the term $\mathcal{J}_{E}=\mathcal{J}_{E,S}+\mathcal{J}_{E,M}$ is given by

$\displaystyle\mathcal{J}_{E,S}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle-\frac{c^{2}}{12\pi}\,\frac{r_{S}M_{\oplus,{\rm
eff}}(r_{E})}{\left||\mathbf{x}-\mathbf{x}_{S}|(\mathbf{x}_{S}-\mathbf{x}_{E})+r_{S}^{2}\hat{\mathbf{n}}_{S}\right|},$
(91) $\displaystyle\mathcal{J}_{E,M}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle-\frac{c^{2}}{12\pi}\,\frac{r_{M}M_{\oplus,{\rm
eff}}(r_{E})}{\left||\mathbf{x}-\mathbf{x}_{M}|(\mathbf{x}_{M}-\mathbf{x}_{E})+r_{M}^{2}\hat{\mathbf{n}}_{M}\right|},$

where $\hat{\mathbf{n}}_{S}$ and $\hat{\mathbf{n}}_{M}$ are the outward unit
normal vectors to the screening surfaces of the Sun and Moon, respectively.

The function $\eta_{E}(\mathbf{x},t)$ depends on time through the centers
$\mathbf{x}_{S}(t),\mathbf{x}_{E}(t),\mathbf{x}_{M}(t)$ of the bodies, which
vary with time along the respective orbits.

The meaning of the terms in the expression (88) of $\eta_{E}$ is the
following: the first term is the surface integral over Earth’s screening
surface of the term in $G^{(1)}(\mathbf{x},\mathbf{x}_{0})$ corresponding to
the image of $\mathbf{x}_{0}$ inside Earth’s screening sphere (the other three
surface integrals vanish); $\mathcal{I}_{E}$ is the volume integral over the
Earth’s shell of the zeroth-order term of Green’s function;
$\mathcal{J}_{E,S}$ and $\mathcal{J}_{E,M}$ are the volume integrals of the
terms in $G^{(1)}(\mathbf{x},\mathbf{x}_{0})$ corresponding to the images of
$\mathbf{x}_{0}$ inside the screening spheres of the Sun and the Moon,
respectively; the last term is the volume integral corresponding to the image
of $\mathbf{x}_{0}$ inside Earth’s screening sphere.

The solution $\eta$ is formally symmetric with respect to the three bodies:

  * (i)

the contribution $\eta_{S}$ from the Sun’s shell is obtained by replacing in
Eqs. (88-90) the subscript $E$ with $S$, the radius $R_{\oplus}$ with
$R_{\odot}$, and the effective mass $M_{\oplus,{\rm eff}}$ with $M_{\odot,{\rm
eff}}$; then

$\mathcal{J}_{E}=\mathcal{J}_{E,S}+\mathcal{J}_{E,M}\,\mbox{is replaced
with}\,\mathcal{J}_{S}=\mathcal{J}_{S,E}+\mathcal{J}_{S,M},$ (92)

in the first of Eqs. (91) the subscripts $E$ and $S$ are exchanged and
$M_{\oplus,{\rm eff}}$ is replaced with $M_{\odot,{\rm eff}}$, and in the
second of Eqs. (91) the subscript $E$ is replaced with $S$ and $M_{\oplus,{\rm
eff}}$ with $M_{\odot,{\rm eff}}$;

  * (ii)

the contribution $\eta_{M}$ from the Moon’s shell is obtained with analogous
changes.

We observe that the terms $\mathcal{J}_{S},\mathcal{J}_{E}$ and
$\mathcal{J}_{M}$ originate from the term $G^{(1)}(\mathbf{x},\mathbf{x}_{0})$
in the Green’s function which involves the first order images.

Eventually, since the various terms in $\eta(\mathbf{x},t)$ decrease as
$1/|\mathbf{x}-\mathbf{x}_{S}|$, $1/|\mathbf{x}-\mathbf{x}_{E}|$ and
$1/|\mathbf{x}-\mathbf{x}_{M}|$ in interplanetary space, and the integral
$I(\eta,G,R^{\ast})$ decreases as $1/R^{\ast}$, then inequality (79) is
satisfied for $|\mathbf{x}-\mathbf{x}_{S}|$ of the order of Earth’s distance
from the Sun and $R^{\ast}>r_{g}$ large enough.

In Appendix C we report the error in the verification of the Dirichlet
condition on $\partial\Omega$ when the solution $\eta$ is computed by using
the Green’s function up to the second order.

###  V.5 Determination of the screening radii

Now we impose the Neumann boundary condition (77) on $\partial\Omega$ and we
use such a condition to find integral equations that determine the three
screening radii $r_{S},r_{E},r_{M}$.

Let us consider Earth’s screening sphere where we have to compute the scalar
product $\nabla\eta\cdot\hat{\mathbf{n}}_{E}$. We compute the contribution to
the scalar product given by the leading terms, the other ones being negligible
because suppressed by factors involving the small radius-to-distance ratios
(and their powers) of the astronomical bodies.

Let us first consider the contribution from the solar term $\eta_{S}$. Since
$\mathcal{I}_{S}$ is the volume integral over the Sun’s shell of the zeroth-
order term (82) of Green’s function, and $\mathcal{J}_{S,E}$ is the volume
integral of the term in $G^{(1)}(\mathbf{x},\mathbf{x}_{0})$ corresponding to
the image of $\mathbf{x}_{0}$ inside the Earth’s screening sphere, then, by
the properties of image points we have
$\mathcal{I}_{S}(\mathbf{x},t)+\mathcal{J}_{S,E}(\mathbf{x},t)=0$ on the
screening sphere, so that the vector
$\nabla\left(\mathcal{I}_{S}+\mathcal{J}_{S,E}\right)$ is orthogonal to the
sphere. By means of a Taylor approximation we have

$\displaystyle\nabla\left(\mathcal{I}_{S}+\mathcal{J}_{S,E}\right)\cdot\hat{\mathbf{n}}_{E}=\frac{c^{2}}{12\pi}\left\\{1-3\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|}\hat{\mathbf{n}}_{ES}\cdot\hat{\mathbf{n}}_{E}\right.$
(93) $\displaystyle-$
$\displaystyle\left.\frac{5}{2}\left[1-3\left(\hat{\mathbf{n}}_{ES}\cdot\hat{\mathbf{n}}_{E}\right)^{2}\right]\left(\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|}\right)^{2}\right\\}\frac{M_{\odot,{\rm
eff}}(r_{S})}{r_{E}|\mathbf{x}_{E}-\mathbf{x}_{S}|}$ $\displaystyle+$
$\displaystyle{\cal
O}\left(\left(\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|}\right)^{3}\right),$

where
$\hat{\mathbf{n}}_{ES}=(\mathbf{x}_{E}-\mathbf{x}_{S})/|\mathbf{x}_{E}-\mathbf{x}_{S}|$,
from which, since

$\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|}<10^{-4},\qquad\frac{r_{E}^{2}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|^{2}}<10^{-8},$
(94)

it follows that the scalar product has the constant dominant term

$\nabla\left(\mathcal{I}_{S}+\mathcal{J}_{S,E}\right)\cdot\hat{\mathbf{n}}_{E}\approx\frac{c^{2}}{12\pi}\,\frac{M_{\odot,{\rm
eff}}(r_{S})}{r_{E}|\mathbf{x}_{E}-\mathbf{x}_{S}|},$ (95)

plus a variable part on the sphere suppressed by the small geometric factor
$r_{E}/|\mathbf{x}_{E}-\mathbf{x}_{S}|$ and its powers.

Let us now consider the integral term in $\eta_{S}$,

$-\frac{c^{2}}{3}\,\frac{r_{S}}{|\mathbf{x}-\mathbf{x}_{S}|}\int_{r_{S}}^{R_{\odot}}\rho_{S}(r)r\,dr\,,$
(96)

obtained by replacing $E$ with $S$ in the last term of the solution (88) for
$\eta_{E}$. By using the method of images it is shown in Appendix C that such
an integral term cancels on the Earth’s screening sphere with a term in
$\eta_{S}$ resulting from the second order Green’s function
$G^{(2)}(\mathbf{x},\mathbf{x}_{0})$. Hence the gradient of the sum of these
two terms is orthogonal to the sphere and the contribution to the scalar
product $\nabla\eta_{S}\cdot\hat{\mathbf{n}}_{E}$, by means of a computation
analogous to Eq. (93), has the constant dominant term

$-\frac{c^{2}}{3}\,\frac{r_{S}}{r_{E}|\mathbf{x}_{E}-\mathbf{x}_{S}|}\int_{r_{S}}^{R_{\odot}}\rho_{S}(r)r\,dr\,,$
(97)

plus a small variable part on the sphere. The surface term in $\eta_{S}$,

$\left(\frac{c^{4}}{8\pi
G}-\eta_{g}\right)\frac{r_{S}}{|\mathbf{x}-\mathbf{x}_{S}|},$ (98)

obtained by replacing $E$ with $S$ in the first term of Eq. (88), gives rise
to an analogous cancellation discussed in Appendix C. Then, arguing as before,
we find a contribution to the scalar product with the constant dominant term

$\frac{r_{S}}{r_{E}|\mathbf{x}_{E}-\mathbf{x}_{S}|}\left(\frac{c^{4}}{8\pi
G}-\eta_{g}\right).$ (99)

Further contributions from $\eta_{S}$ turn out to be negligible. The
contribution from the lunar term $\eta_{M}$ is obtained by replacing $S$ with
$M$ in the previous expressions. Eventually, the contribution from the
terrestrial term $\eta_{E}$ is given by

$\nabla\eta_{E}\cdot\hat{\mathbf{n}}_{E}\approx\frac{c^{2}}{3r_{E}}\int_{r_{E}}^{R_{\oplus}}\rho_{E}(r)r\,dr-\frac{1}{r_{E}}\left(\frac{c^{4}}{8\pi
G}-\eta_{g}\right).$ (100)

By substituting all these contributions in the scalar product
$\nabla\eta\cdot\hat{\mathbf{n}}_{E}$, and imposing the Neumann boundary
condition (77) on Earth’s screening sphere, we obtain the following integral
equation:

$\displaystyle\frac{c^{2}}{3}\int_{r_{E}}^{R_{\oplus}}\rho_{E}(r)r\,dr-\left(1-\frac{r_{S}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|}-\frac{r_{M}}{|\mathbf{x}_{E}-\mathbf{x}_{M}|}\right)\left(\frac{c^{4}}{8\pi
G}-\eta_{g}\right)$ $\displaystyle+$
$\displaystyle\frac{c^{2}}{3}\left[\frac{1}{|\mathbf{x}_{E}-\mathbf{x}_{S}|}\left(\frac{M_{\odot,{\rm
eff}}(r_{S})}{4\pi}-r_{S}\int_{r_{S}}^{R_{\odot}}\rho_{S}(r)r\,dr\right)+\frac{1}{|\mathbf{x}_{E}-\mathbf{x}_{M}|}\left(\frac{M_{M,{\rm
eff}}(r_{M})}{4\pi}-r_{M}\int_{r_{M}}^{R_{M}}\rho_{M}(r)r\,dr\right)\right]=0.$

Then, repeating the computation on the screening spheres of the Sun and Moon,
we obtain a system of integral equations. Such equations are obtained from Eq.
(V.5) by exchanging the subscripts $S,E,M$. Using formula (52) for the
minimizer $\eta_{g}$, the solution to the resulting system of integral
equations determines the three screening radii $r_{E},r_{S},r_{M}$ for given
values of the NMC gravity parameters $\alpha$ and $q$. The system of equations
generalizes the integral equation found in Ref. MBMGDeA  to a system of three
gravitationally interacting extended bodies.

Neglecting factors involving radius-to-distance ratios, the following
relations follow from Eq. (V.5) and the other integral equations, as a first
approximation:

$\int_{r_{S}}^{R_{\odot}}\rho_{S}(r)r\,dr\approx\int_{r_{E}}^{R_{\oplus}}\rho_{E}(r)r\,dr\approx\int_{r_{M}}^{R_{M}}\rho_{M}(r)r\,dr\,.$
(102)

If we also assume that all bodies have a thin shell, so that we have
$R_{\oplus}-r_{E}\ll R_{\oplus}$ for the Earth, and analogous inequalities for
the Sun and Moon, then the above approximate relations are equivalent to the
following relations between effective masses:

$\frac{M_{\odot,{\rm eff}}(r_{S})}{R_{\odot}}\approx\frac{M_{\oplus,{\rm
eff}}(r_{E})}{R_{\oplus}}\approx\frac{M_{M,{\rm eff}}(r_{M})}{R_{M}}\,.$ (103)

Eventually, we observe that the results obtained in this section show that the
normal derivative of $\eta$ on the screening spheres has a main part which is
constant on each sphere, and a much smaller variable part which can be
neglected thanks to inequalities of type (94). Since the interior solution
computed in Sec. IV is spherically symmetric inside each body, then it has a
constant normal derivative on each screening sphere, so that matching the
solution outside of the screening radii with the interior solution yields a
consistent approximation, as it was anticipated at the beginning of Sec. IV.

###  V.6 Verification of inequalities

The solution for $\eta$ has been computed by assuming inequalities (18-19),
(62) and (66), necessary in order to find an analytic approximation of the
solution, that have to be verified a posteriori. In this section we show that
the computed solution satisfies inequalities (18-19) and (62), while
inequality (66) is verified in Appendix E.

Let us first consider inequality (19) and the solution inside the screening
radii. In the case of Earth, using the value (41) of $\eta$ minimizing the
effective potential, the solution for $\eta$ found in Sec. IV.2 satisfies

$\left|\frac{8\pi
G}{c^{4}}\eta-1\right|\sim\frac{G}{c^{2}}\lambda^{2}(\rho_{E})\rho_{E}\sim\left(\frac{\lambda(\rho_{E})}{R_{\oplus}}\right)^{2}\frac{GM_{\oplus}}{c^{2}R_{\oplus}}\ll
1.$ (104)

Analogous results can be found for the Moon and the Sun. Then, using now the
solution outside of the screening radii found in Sec. V.4, for
$|\mathbf{x}-\mathbf{x}_{S}|<r_{g}$ the expression $|(8\pi G/c^{4})\eta-1|$ is
bounded by a sum of terms each of which is bounded by a quantity of type

$\mbox{either}\quad\frac{GM_{\rm eff}}{c^{2}R}\ll
1,\quad\mbox{or}\quad\frac{GM_{\rm eff}}{c^{2}d}\ll 1,$ (105)

where $M_{\rm eff}$ is an effective mass, $R$ is the radius of a body, and $d$
is a distance between the astronomical bodies. Moreover, using formula (52)
for the minimizer $\eta_{g}$ and the integral equations (V.5), we have $|(8\pi
G/c^{4})\eta_{g}-1|\ll 1$. Since $\eta$ is a harmonic function for
$|\mathbf{x}-\mathbf{x}_{S}|>r_{g}$ by Eq. (54), and $\eta\approx\eta_{g}$ at
large distance from the Sun by the boundary condition (55), then by the
maximum principle for harmonic functions $\eta$ satisfies the desired
inequality also for $|\mathbf{x}-\mathbf{x}_{S}|>r_{g}$. Hence the computed
solution $\eta$ satisfies inequality (19) everywhere.

We now consider the second of inequalities (18). We set $R_{\rm GR}=8\pi
G\rho/c^{2}$ and we observe that $R\approx R_{\rm GR}$ both inside the
screening radii and in the solar neighborhood of the Galaxy, $R\ll R_{\rm GR}$
in the thin shells (see Appendix E), and $R$ interpolates between such values
in interplanetary space, so that $R/R_{\rm GR}\lesssim 1$ everywhere. Then,
using formula (57), we have

$\left|f^{2}(R)\right|=|q|R^{\alpha}=\frac{1}{2|\alpha|}\,\frac{R}{R_{\rm
GR}}\left|\frac{8\pi G}{c^{4}}\eta-1\right|\ll 1,$ (106)

for $|\alpha|$ not too close to zero, which is the case of interest for
applications to astrophysics drkmattgal  and cosmology curraccel . Then the
second of inequalities (18) is satisfied, the first being trivial.

Eventually we verify inequality (62) under the assumption that the thin shell
condition is satisfied for all the astronomical bodies, a condition that will
follow from the constraint from LLR measurements. First we observe that the
leading terms of the solution $\eta$ in interplanetary space are given by

$\eta-\eta_{g}\approx\frac{c^{2}}{12\pi}\left(\frac{M_{\odot,{\rm
eff}}}{|\mathbf{x}-\mathbf{x_{S}}|}+\frac{M_{\oplus,{\rm
eff}}}{|\mathbf{x}-\mathbf{x_{E}}|}+\frac{M_{M,{\rm
eff}}}{|\mathbf{x}-\mathbf{x_{M}}|}\right),$ (107)

where we have dropped the dependence of the effective masses on the screening
radii for simplicity. Then, using the integral equations that determine the
screening radii under the approximation (102), the thin shell condition and
the definition (32) of $\lambda_{g}$, we find

$M_{\oplus,{\rm eff}}\approx 4\pi
R_{\oplus}\int_{r_{E}}^{R_{\oplus}}\rho_{E}(r)rdr\approx\frac{4\pi}{1-\alpha}R_{\oplus}\lambda_{g}^{2}\rho_{g}.$
(108)

Using the approximate relations (103) between effective masses, we have

$\displaystyle\eta-\eta_{g}\approx\frac{c^{2}\lambda_{g}^{2}\rho_{g}}{3(1-\alpha)}\left(\frac{R_{\odot}}{|\mathbf{x}-\mathbf{x_{S}}|}+\frac{R_{\oplus}}{|\mathbf{x}-\mathbf{x_{E}}|}+\frac{R_{M}}{|\mathbf{x}-\mathbf{x_{M}}|}\right),$

from which, for $\alpha<0$, inequality (62) is satisfied if

$\frac{R_{\odot}}{|\mathbf{x}-\mathbf{x_{S}}|}+\frac{R_{\oplus}}{|\mathbf{x}-\mathbf{x_{E}}|}+\frac{R_{M}}{|\mathbf{x}-\mathbf{x_{M}}|}\ll
1,$ (110)

which is satisfied in interplanetary space at large enough distance from the
astronomical bodies.

###  V.7 Solution for the potentials $\Phi$ and $\Psi$

We have assumed the field $\eta$ is close to the minimizer of the effective
potential $V_{\rm eff}$ in the solar neighbourhood of the Galaxy for
$|\mathbf{x}-\mathbf{x}_{S}|>r_{g}$, so that GR is approximately satisfied.

In what follows we denote by $U$ the Newtonian potential of the mass
distribution with density $\rho$,

$U(\mathbf{x},t)=G\int\frac{\rho(\mathbf{y},t)}{|\mathbf{x}-\mathbf{y}|}d^{3}y,$
(111)

which satisfies the Poisson equation $\nabla^{2}U=-4\pi G\rho$.

Using Eqs. (21-22) it follows that the potential $\Psi$ of the metric is
related to the deviation from GR, then we impose the following boundary
conditions in the Galaxy at large distance $|\mathbf{x}-\mathbf{x}_{S}|$ from
the Sun’s center, where GR is satisfied by our assumptions:

$\Phi(\mathbf{x},t)\approx\frac{1}{c^{2}}U(\mathbf{x},t),\qquad\Psi(\mathbf{x},t)\approx
0.$ (112)

Combining equations (21) and (24) for $\Psi$ and $\eta$ we have

$\nabla^{2}\left(\Psi+\frac{8\pi G}{c^{4}}\,\eta\right)=0.$ (113)

Using now the second of boundary conditions (112) for $\Psi$ and the boundary
condition (55) for $\eta$, for all points $\mathbf{x}$ on a sphere with center
in $\mathbf{x}_{S}$ and large enough radius $R^{\ast}>r_{g}$, and for any time
$t$, we have

$\Psi(\mathbf{x},t)+\frac{8\pi G}{c^{4}}\,\eta(\mathbf{x},t)\approx\frac{8\pi
G}{c^{4}}\,\eta_{g}.$ (114)

Hence, by the maximum principle for harmonic functions, it follows that the
harmonic function (with respect to the spatial variables $\mathbf{x}$)

$\Psi(\mathbf{x},t)+\frac{8\pi G}{c^{4}}\,\eta(\mathbf{x},t)$ (115)

is constant inside the sphere of radius $R^{\ast}$, so that the solution for
the potential $\Psi$ is given by

$\Psi(\mathbf{x},t)=-\frac{8\pi
G}{c^{4}}\left[\eta(\mathbf{x},t)-\eta_{g}\right].$ (116)

The solution for $\Psi$ then follows immediately from the solution for $\eta$
found in the previous sections. Combining now equations (21) and (22) for
$\Psi$ and $\Phi$ we have

$\nabla^{2}\left(\Phi-\frac{U}{c^{2}}-\frac{1}{2}\,\Psi\right)=0.$ (117)

Then, applying both boundary conditions (112) and using again the maximum
principle for harmonic functions, the solution $\Phi$ of this equation is
given by

$\Phi(\mathbf{x},t)=\frac{1}{c^{2}}\,U(\mathbf{x},t)+\frac{1}{2}\,\Psi(\mathbf{x},t).$
(118)

The solutions found for $\Phi$ and $\Psi$ define the space-time metric (6).

##  VI Dynamics of continuous bodies

We consider the motion of Earth and Moon in the gravitational field of the
Sun. The equations describing the dynamics of the system are obtained by
taking the covariant divergence of the energy-momentum tensor and applying
Bianchi identities to the gravitational field equations (see Ref. BBHL ), as
given by Eq. (5) that we repeat for convenience:

$\nabla_{\mu}T^{\mu\nu}=\frac{f^{2}_{R}}{1+f^{2}}(g^{\mu\nu}\mathcal{L}_{m}-T^{\mu\nu})\nabla_{\mu}R.$
(119)

In the following computation we neglect mass density of solar and terrestrial
atmospheres, and of solar wind in interplanetary space, so that density
$\rho(\mathbf{x},t)$ has a compact support consisting of the three spheres of
radii $R_{\odot},R_{\oplus}$ and $R_{M}$. We assume densities assigned inside
the bodies according to the profiles given in Appendix A, and we assume all
the astronomical bodies in hydrostatic equilibrium. According to Sec. II.1 the
Sun is considered as a perfect fluid, while the Earth and Moon are
approximately described as continuous bodies in a hydrostatic state of stress,
so that inside the bodies the equations describing the dynamics of continuous
media formally coincide with the equations of hydrodynamics of a perfect
fluid.

We begin by computing the 0th component of Eq. (119): using the components of
the energy-momentum tensor given by Eqs. (7) and (8), the left-hand side of
this equation yields

$\nabla_{\mu}T^{\mu 0}=c\frac{\partial\rho}{\partial
t}+c\frac{\partial}{\partial x^{i}}(\rho v^{i})+{\cal
O}\left(\frac{1}{c}\right).$ (120)

Now we compute the right-hand side of Eq. (119). First we observe that inside
the screening surfaces we have $R={\cal O}(1/c^{2})$. Then, the integral
equations (V.5) imply

$|q|=\rho_{g}^{-\alpha}\cdot{\cal O}\left(\frac{1}{c^{2-2\alpha}}\right),$
(121)

being $\alpha<0$, from which, using formula (58) for curvature
$R=\omega(\eta,\rho)$ and the property $|(8\pi G/c^{4})\eta-1|={\cal
O}(1/c^{2})$ found in Sec. V.6, we have $R={\cal O}(1/c^{2})$ also in the thin
shells of the astronomical bodies. Using now property (121) and the definition
$f^{2}(R)=qR^{\alpha}$, it follows

$f^{2}(R)=\rho_{g}^{|\alpha|}\cdot{\cal O}\left(\frac{1}{c^{2}}\right),\quad
f^{2}_{R}=\rho_{g}^{|\alpha|}\cdot{\cal O}(1),$ (122)

from which the evaluation of the right-hand side of Eq. (119) yields

$\frac{f^{2}_{R}}{1+f^{2}}(g^{\mu 0}\mathcal{L}_{m}-T^{\mu 0})\frac{\partial
R}{\partial x^{\mu}}={\cal O}\left(\frac{1}{c}\right).$ (123)

Using Eq. (120) and neglecting terms of order ${\cal O}(1/c^{2})$, the
continuity equation then follows in the nonrelativistic limit as usual:

$\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x^{i}}(\rho
v^{i})=0.$ (124)

The NMC term on the right-hand side of Eq. (119) gives a distinctive
contribution to the spatial part of this equation that now we compute. Using
the components (7)-(9) of the energy-momentum tensor, for $i=1,2,3$ the left-
hand side yields

$\displaystyle\nabla_{\mu}T^{\mu i}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial t}(\rho v^{i})+\frac{\partial}{\partial
x^{j}}(\rho v^{i}v^{j})-\rho\frac{\partial U}{\partial
x_{i}}+\frac{c^{2}}{2}\rho\frac{\partial\Psi}{\partial x_{i}}$ (125)
$\displaystyle+$ $\displaystyle\frac{\partial p}{\partial x_{i}}+{\cal
O}\left(\frac{1}{c^{2}}\right).$

Using now the continuity equation Eq. (124), at order $O(1)$, we get

$\nabla_{\mu}T^{\mu i}=\rho\frac{dv^{i}}{dt}+\frac{\partial p}{\partial
x_{i}}-\rho\frac{\partial U}{\partial
x_{i}}+\frac{c^{2}}{2}\rho\frac{\partial\Psi}{\partial x_{i}}+{\cal
O}\left(\frac{1}{c^{2}}\right),$ (126)

where $d/dt=\partial/\partial t+v^{i}\partial/\partial x^{i}$ is the material
derivative of continuum mechanics.

For $i=1,2,3$, taking into account that $|f^{2}|\ll 1$, the right-hand side of
Eq. (119) yields

$\frac{f^{2}_{R}}{1+f^{2}}(g^{\mu i}\mathcal{L}_{m}-T^{\mu i})\frac{\partial
R}{\partial x^{\mu}}=-c^{2}f^{2}_{R}\,\rho\frac{\partial R}{\partial
x_{i}}+{\cal O}\left(\frac{1}{c^{2}}\right).$ (127)

Combining equations (126) and (127), and neglecting terms of order ${\cal
O}(1/c^{2})$, we obtain the equations of NMC dynamics of continuous bodies in
hydrostatic state of stress and in the nonrelativistic limit,

$\rho\frac{d\mathbf{v}}{dt}=\rho\nabla U-\nabla p-\frac{1}{2}\rho
c^{2}\nabla\Psi-c^{2}f^{2}_{R}\,\rho\nabla R,$ (128)

where the vector notation has been used. These equations are the Eulerian
equations of Newtonian hydrodynamics with the presence of two additional
terms:

  * (i)

a fifth force density proportional to the gradient of the metric potential
$\Psi$;

  * (ii)

an extra force density proportional to the product of $f^{2}_{R}$ by the
gradient of curvature $R$.

The extra force density in (ii) has been extensively discussed in Ref. BBHL ,
and for relativistic perfect fluids in Ref. BLP . While the fifth force is
typical of $f(R)$ gravity theory, the extra force is specific of NMC gravity
for the choie Eq. (11).

###  VI.1 Fifth force inside the screening spheres

The fifth force density is given by

$-\frac{1}{2}\rho c^{2}\nabla\Psi=\frac{4\pi G}{c^{2}}\rho\nabla\eta,$ (129)

where the expression (116) of $\Psi$ in terms of the function $\eta$ has been
used. In the interior of the screening spheres the magnitude of fifth force is
expected to be largest at the surfaces of density discontinuity where the
largest deviations from GR take place. Let us compute such a force at the
Moho, the crust-mantle discontinuity in the Earth’s interior.

The radial derivative of $\eta$ is continuous, and using the solution for
$\eta$ inside the Earth’s screening radius, found in Sec. IV.2, we have

$\frac{d\eta}{dr}(R_{E,m})\approx\frac{\eta_{E,c}-\eta_{E,m}}{\lambda_{E,c}+\lambda_{E,m}}.$
(130)

The quantity $\rho_{E}d\eta/dr$ is discontinuous at the Moho and the maximum
value of its magnitude is $\rho_{E,m}|d\eta/dr(R_{E,m})|$. Then, using the
value (41) of $\eta$ minimizing the effective potential $V_{\rm eff}$ and
formula (32) for $\lambda$, the magnitude of the fifth force density is
estimated by

$\displaystyle\left|F_{\rm f}\right|$ $\displaystyle\approx$
$\displaystyle\frac{4}{3}\,\frac{\pi
G}{1+|\alpha|}\rho_{E,m}\frac{|\lambda^{2}_{E,c}\rho_{E,c}-\lambda^{2}_{E,m}\rho_{E,m}|}{\lambda_{E,c}+\lambda_{E,m}}$
(131) $\displaystyle\leq$
$\displaystyle\frac{8\pi}{3}G\lambda_{E,c}\rho_{E,m}\rho_{E,c}.$

Now the magnitude of the Newtonian force density at the Moho is

$\left|F_{N}\right|=\rho_{E,m}\left|\frac{dU}{dr}(R_{E,m})\right|=\frac{4\pi}{3}GR_{E,m}\rho_{E,m}\overline{\rho}_{E,c},$
(132)

where $\overline{\rho}_{E,c}$ is the average Earth density below the crust.
Since we have $\lambda(\rho_{E})\ll R_{E,m}$ with $\lambda(\rho_{E})$
completely negligible ($\lambda(\rho_{E})\lesssim 10^{-7}$ m for $\alpha=-1$
according to Section .2 of Appendix A), it then follows $\left|F_{\rm
f}\right|\ll\left|F_{N}\right|$. The same behavior is found at the other
density discontinuities, moreover, far from discontinuities the fifth force is
further decreased by exponential suppression with decay constant
$1/\lambda(\rho_{E})$ in the various layers. Analogous results hold for the
Moon, while for the Sun see Ref. MBMGDeA . Eventually, the perturbation of the
Newtonian gravitational force, hence of hydrostatic equilibrium, is completely
negligible inside the screening radii, confirming the effectiveness of the
screening mechanism.

###  VI.2 Jump conditions for the pressure

The expression (58) for curvature $R=\omega(\eta,\rho)$, which we rewrite in
the form

$R=\left(\frac{c^{2}/(8\pi G)-\eta/c^{2}}{2\alpha
q\rho}\right)^{\frac{1}{\alpha-1}},$ (133)

shows that $R$ is discontinuous at surfaces across which mass density is
discontinuous, such as the external surfaces of Earth and Moon or the
Mohorovičić discontinuity. Since the extra force is proportional to $\nabla
R$, then such a force is concentrated at the surfaces of density discontinuity
and the concentration gives rise to a jump of pressure $p$ across these
surfaces that we now compute. We will find that such a jump of pressure is
undetectable for interesting values of parameter $\alpha$.

In the following the astronomical body considered is either the Earth or the
Moon. Let $R_{\rm d}$ be the radius of a discontinuity surface, and let
$\rho^{+}$ and $\rho^{-}$ be the values of density on the two sides of the
discontinuity. We adopt the method of Ref. Obrien-Sy  and we introduce a
boundary layer across which density changes continuously, then we proceed to
the limit in which the thickness of the layer tends to zero. Let $\varepsilon$
be a positive small parameter and let $r$ be the distance from the center of
the body, then we introduce a family of mass density functions
$\rho_{\varepsilon}$ defined by

$\begin{cases}\rho_{\varepsilon}(r)=\rho(r)\quad\mbox{for }r<R_{\rm
d}-\varepsilon\mbox{ and }r>R_{\rm d}+\varepsilon,\\\
\rho_{\varepsilon}(R_{\rm
d}+\varepsilon)=\rho^{+},\quad\rho_{\varepsilon}(R_{\rm
d}-\varepsilon)=\rho^{-},\end{cases}$ (134)

where $\rho_{\varepsilon}(r)$ is continuously differentiable in the interval
$(R_{\rm d}-\varepsilon,R_{\rm d}+\varepsilon)$. Moreover, we denote by means
of the subscript $\varepsilon$ the functions which are solution of the field
equations and equations of motion for assigned density $\rho_{\varepsilon}$,
and we assume that such functions converge as $\varepsilon\to 0$ to the
solution of the problem with density $\rho$.

Let now
$\hat{\mathbf{n}}=(\mathbf{x}-\mathbf{x}_{o})/|\mathbf{x}-\mathbf{x}_{o}|$,
where $\mathbf{x}_{o}$ is the position vector of the center of the body, then
replacing $\rho$ with $\rho_{\varepsilon}$ in Eqs. (128), taking the scalar
product of the resulting equations by the unit vector $\hat{\mathbf{n}}$, and
integrating radially over the boundary layer of thickness $2\varepsilon$, at
the time instant $t$ we have

$\displaystyle\Delta p_{\varepsilon}$ $\displaystyle=$
$\displaystyle-c^{2}\int_{R_{\rm d}-\varepsilon}^{R_{\rm
d}+\varepsilon}\rho_{\varepsilon}f^{2}_{R_{\varepsilon}}\frac{\partial
R_{\varepsilon}}{\partial r}dr$ $\displaystyle+$ $\displaystyle\int_{R_{\rm
d}-\varepsilon}^{R_{\rm
d}+\varepsilon}\rho_{\varepsilon}\left[-\frac{d\mathbf{v}_{\varepsilon}}{dt}\cdot\hat{\mathbf{n}}+\frac{\partial
U_{\varepsilon}}{\partial
r}-\frac{c^{2}}{2}\frac{\partial\Psi_{\varepsilon}}{\partial r}\right]dr,$

where $\Delta p_{\varepsilon}=p_{\varepsilon}(R_{\rm
d}+\varepsilon,\theta,\varphi,t)-p_{\varepsilon}(R_{\rm
d}-\varepsilon,\theta,\varphi,t)$ and $(r,\theta,\varphi)$ are spherical
coordinates on the discontinuity surface. Taking the limit as
$\varepsilon\rightarrow 0$, and observing that the derivatives in the
integrand in the second row of the above equation are bounded above and below
for all $\varepsilon$, we find

$\Delta p=-c^{2}\lim_{\varepsilon\rightarrow 0}\int_{R_{\rm
d}-\varepsilon}^{R_{\rm
d}+\varepsilon}\rho_{\varepsilon}f^{2}_{R_{\varepsilon}}\frac{\partial
R_{\varepsilon}}{\partial r}dr,$ (136)

where $\Delta p=p(R_{\rm d}^{+},\theta,\varphi,t)-p(R_{\rm
d}^{-},\theta,\varphi,t)$ is the pressure jump. In order to compute the
pressure jump we now compute the above limit of the integral. Using
$f^{2}(R_{\varepsilon})=qR_{\varepsilon}^{\alpha}$ and Eq. (133), the
integrand is given by

$-\left(2|\alpha|\right)^{\frac{\alpha}{1+|\alpha|}}|q|^{\frac{1}{1+|\alpha|}}\rho_{\varepsilon}\frac{\partial}{\partial
r}\left[\rho_{\varepsilon}^{\frac{\alpha}{1+|\alpha|}}\left(\frac{c^{2}}{8\pi
G}-\frac{\eta_{\varepsilon}}{c^{2}}\right)^{\frac{\alpha}{\alpha-1}}\right],$
(137)

which can be decomposed into the sum of a bounded term whose contribution
vanishes in the limit as $\varepsilon\rightarrow 0$, and an unbounded term
which in the limit gives rise to the pressure jump. The result of the
computation is

$\Delta p=\mathcal{C}(\alpha,q)\left(\frac{c^{2}}{8\pi
G}-\frac{\eta}{c^{2}}\right)^{\frac{\alpha}{\alpha-1}}\left({\rho^{-}}^{\frac{1}{1+|\alpha|}}-{\rho^{+}}^{\frac{1}{1+|\alpha|}}\right),$
(138)

where the coefficient $\mathcal{C}(\alpha,q)$ is given by

$\mathcal{C}(\alpha,q)=c^{2}2^{\frac{\alpha}{1-\alpha}}\left(\alpha
q\right)^{\frac{1}{1+|\alpha|}},$ (139)

and the function $\eta=\eta(R_{\rm d},\theta,\varphi,t)$ is evaluated on the
discontinuity surface at the point with angular coordinates
$(\theta,\varphi)$. The jump $\Delta p$ does not depend on the angular
coordinates on the discontinuity surfaces which are located inside the
screening spheres, while the jump varies with position on the external
surfaces of Earth and Moon.

We see that the jump of pressure depends both on the solution $\eta$ and,
explicitly, on the NMC gravity parameters $\alpha$ and $q$ (take into account
that in Sec. III.1 we have found that $\alpha<0$ implies $q<0$).

###  VI.3 Extra force inside the screening spheres

We consider the perturbation of the hydrostatic equilibrium due to extra force
in the interior of astronomical bodies. Let us apply the jump condition for
pressure at the mantle-core discontinuity inside the Moon. Separating the
motion of the center of mass from the internal condition of hydrostatic
equilibrium in Eqs. (128), integrating radially and taking into account the
pressure jump, we find

$\displaystyle p(R_{\rm d}^{+})-p(0)-\Delta p=\int_{0}^{R_{\rm
d}}\rho_{M}\frac{\partial U_{M}}{\partial r}dr$ (140)
$\displaystyle-\frac{c^{2}}{2}\int_{0}^{R_{\rm
d}}\rho_{M}\frac{\partial\Psi_{M}}{\partial r}dr-c^{2}\int_{0}^{R_{\rm
d}}\rho_{M}f^{2}_{R}\left[\frac{\partial R}{\partial r}\right]dr,$

where $R_{d}=R_{M,n}$ is the radius of the mantle-core interface, $p(0)$ is
pressure at the Moon’s center, $U_{M}$ is the Newtonian potential sourced by
Moon’s mass density, $\Psi_{M}$ is the potential inside the Moon’s screening
radius, and $[\partial R/\partial r]$ denotes the derivative $\partial
R/\partial r$ outside of the discontinuity surface.

In order to compute the pressure jump $\Delta p$ we use Eq. (138) and the
solution (51) for $\eta$ inside the Moon’s core:

$\displaystyle\frac{c^{2}}{8\pi
G}-\frac{1}{c^{2}}\eta(R_{M,n})=\frac{1}{3(1+|\alpha|)}\,\frac{\lambda_{M,m}\lambda_{M,n}}{\lambda_{M,m}+\lambda_{M,n}}$
$\displaystyle\times$
$\displaystyle\left(\lambda_{M,m}\rho_{M,m}+\lambda_{M,n}\rho_{M,n}\right)<\frac{2}{3(1+|\alpha|)}\,\lambda_{M,m}^{2}\rho_{M,m},$

from which we obtain the following estimate for the pressure jump:

$|\Delta p|<\frac{16}{3}\,\pi G\lambda_{M,m}^{2}\rho_{M,n}^{2}.$ (142)

For the Newtonian term we have

$\int_{0}^{R_{\rm d}}\rho_{M}\frac{\partial U_{M}}{\partial
r}dr=-\frac{2}{3}\,\pi GR_{M,n}^{2}\rho_{M,n}^{2},$ (143)

from which, being $\lambda_{M,m}\ll R_{M,n}$ with $\lambda_{M,m}$ completely
negligible ($\lambda(\rho_{M})\lesssim 10^{-7}$ m for $\alpha=-1$ according to
Section .3 of Appendix A), it then follows

$|\Delta p|\ll\left|\int_{0}^{R_{\rm d}}\rho_{M}\frac{\partial U_{M}}{\partial
r}dr\right|.$ (144)

The volume contribution of the extra force inside Moon’s core is given by

$-c^{2}\int_{0}^{R_{\rm d}}\rho_{M}f^{2}_{R}\left[\frac{\partial R}{\partial
r}\right]dr=-c^{2}\rho_{M,n}\left[f^{2}(R_{d}^{-})-f^{2}(0)\right],$ (145)

Using expression (133) of curvature $R$, the solution (51) for $\eta$ inside
Moon’s core, taking into account that $\lambda_{M,n}\ll R_{d}$ implies
$\eta(0)\approx\eta_{M,n}$, we have

$|f^{2}(R_{d}^{-})|>|f^{2}(0)|,$ (146)

then we compute the following estimate:

$\displaystyle c^{2}\left|\int_{0}^{R_{\rm
d}}\rho_{M}f^{2}_{R}\left[\frac{\partial R}{\partial
r}\right]dr\right|<2c^{2}\rho_{M,n}|f^{2}(R_{d}^{-})|$
$\displaystyle<\frac{16}{3}2^{\frac{1}{\alpha-1}}\pi
G\left(\frac{\rho_{M,n}}{\rho_{M,m}}\right)^{\frac{1}{1+|\alpha|}}\frac{\lambda_{M,m}^{2}}{|\alpha|}\rho^{2}_{M,m}.$
(147)

Since $\lambda_{M,m}^{2}/|\alpha|\ll R_{M,n}^{2}$ for values of $\alpha$ of
interest for application to astrophysics and cosmology (see Refs. drkmattgal
,curraccel  and Section .3 of Appendix A), then for such values the volume
contribution of the extra force turns out to be completely negligible in
comparison with the Newtonian term (143).

The same behavior is found in other regions of the Moon’s interior and in the
Earth’s interior, while for the Sun see Ref. MBMGDeA . The effectiveness of
the screening mechanism is then confirmed by taking into account also the
extra force.

###  VI.4 Motion of centers of mass

In the following we denote by ${\cal V}_{E}(t)$ and ${\cal V}_{M}(t)$ the
regions of space occupied by Earth and Moon, respectively, at the time instant
$t$. By using the continuity equation and Reynolds transport theorem of
continuum mechanics, we have

$M_{\oplus}\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}=\int_{{\cal
V}_{E}(t)}\frac{d\mathbf{v}}{dt}\rho_{E}(\mathbf{x},t)d^{3}x,$ (148)

where $M_{\oplus}$ is the mass of Earth, and an analogous equation holds for
the Moon. By substituting the expression of $\rho d\mathbf{v}/dt$ given by Eq.
(128) we obtain the various contributions to the acceleration of Earth and
Moon. Because of the presence of discontinuities, the integrals involving the
pressure gradient and the extra force are computed as follows. We consider the
Earth motion. The results of Section VI.2 show that the extra force makes
pressure discontinuous across the external surface of Earth, and the pressure
jump,

$\Delta p=p_{\rm atm}-p(R_{\oplus}^{-},\theta,\varphi,t),$ (149)

where $p_{\rm atm}$ is atmospheric pressure at sea level, is given by formula
(138) with $\rho^{-}=\rho_{E,w}$ and $\rho^{+}=\rho_{\rm atm}$, the
atmospheric density at sea level.

The contribution of the extra force to the integral on the right-hand side of
Eq. (148) is then given by

$-c^{2}\int_{{\cal V}_{E}(t)}\rho_{E}f^{2}_{R}\left[\nabla
R\right]d^{3}x-\int_{\partial{\cal
V}_{E}(t)}p(R_{\oplus}^{-},\theta,\varphi,t)\hat{\mathbf{n}}_{E}d\sigma,$
(150)

where $[\nabla R]$ denotes the vector function $\nabla R$ outside of the
discontinuity surfaces, and we have taken into account that the contributions
of pressure jumps across the discontinuity surfaces in the Earth’s interior
vanish because of spherical symmetry of curvature $R$, assuming such surfaces
contained in the screening sphere. The surface integral of
$p(R_{\oplus}^{-},\theta,\varphi,t)$ is computed by using the expression (138)
of $\Delta p$ taking into account that the integral of $p_{\rm atm}$ vanishes
because of spherical symmetry, assuming $p_{\rm atm}$ uniform on the Earth
surface.

Then, substituting the expression of $\rho d\mathbf{v}/dt$ given by Eq. (128)
into the integral in Eq. (148), and using Eq. (150), we obtain

$\displaystyle M_{\oplus}\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}$ $\displaystyle=$
$\displaystyle\int_{{\cal V}_{E}(t)}\rho_{E}\nabla
U\,d^{3}x-\frac{c^{2}}{2}\int_{{\cal V}_{E}(t)}\rho_{E}\nabla\Psi\,d^{3}x$
$\displaystyle-$ $\displaystyle\int_{\partial{\cal
V}_{E}(t)}p^{-}\hat{\mathbf{n}}_{E}d\sigma-c^{2}\int_{{\cal
V}_{E}(t)}\rho_{E}f^{2}_{R}\left[\nabla R\right]d^{3}x,$

where $p^{-}=p(R_{\oplus}^{-},\theta,\varphi,t)$. The first integral is the
contribution of Newtonian gravity to the acceleration of Earth, the second
integral is the contribution of fifth force, the surface integral is the
contribution of the extra pressure on Earth’s surface, and the last integral
is the volume contribution of the extra force. An analogous formula holds for
the Moon.

In the next section we evaluate the integral which gives the contribution of
the fifth force.

##  VII Acceleration of Earth and Moon due to the fifth force

The fifth force contribution to the acceleration of Earth is given by

$M_{\oplus}\left(\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}\right)_{\rm
f}=-\frac{c^{2}}{2}\int_{{\cal
V}_{E}(t)}\rho_{E}(\mathbf{x},t)\nabla\Psi(\mathbf{x},t)d^{3}x,$ (152)

and an analogous integral holds for the Moon. In Sec. VI.1 we have found that
the contribution to the fifth force from the interior of the screening sphere
is negligible due to the smallness of $\lambda(\rho_{E})$, so that the
contribution to the integral over ${\cal V}_{E}(t)$ only comes from the thin
shell defined by $r_{E}<|\mathbf{x}-\mathbf{x}_{E}|<R_{\oplus}$.

The integral over the thin shell is evaluated by using the expression (116) of
the potential $\Psi$ in terms of $\eta$, which gives $\nabla\Psi=-(8\pi
G/c^{4})\nabla\eta$, and using the solution for the function $\eta$ computed
by means of the method of images.

Let us first consider the contribution from the solar term $\eta_{S}$. Using
spherical coordinates we find the values of the following integrals:

$\frac{4\pi G}{c^{2}}\int_{{\cal
V}_{E}(t)}\rho_{E}\nabla\mathcal{I}_{S}d^{3}x=\frac{G}{3}M_{\odot,{\rm
eff}}M_{\oplus,{\rm
eff}}\frac{\mathbf{x}_{S}-\mathbf{x}_{E}}{|\mathbf{x}_{S}-\mathbf{x}_{E}|^{3}}\,,$
(153)

where we have dropped the dependence of effective masses on the screening
radii, and

$\int_{{\cal V}_{E}(t)}\rho_{E}\nabla\mathcal{J}_{S,E}\,d^{3}x=0.$ (154)

Let us now consider the contribution from the following terms of $\eta_{S}$:

$\frac{r_{S}}{|\mathbf{x}-\mathbf{x}_{S}|}\left(\frac{c^{4}}{8\pi
G}-\eta_{g}-\frac{c^{2}}{3}\int_{r_{S}}^{R_{\odot}}\rho_{S}(r)r\,dr\right).$
(155)

Using the integral equation on the Sun’s screening sphere, which is obtained
from Eq. (V.5) by exchanging $E$ with $S$, we have

$\frac{c^{2}}{3}\int_{r_{S}}^{R_{\odot}}\rho_{S}(r)r\,dr=\frac{c^{4}}{8\pi
G}-\eta_{g}+{\cal O}\left(\frac{R}{d}\right),$ (156)

where ${\cal O}(R/d)$ denotes terms multiplied by a factor of type $R/d$,
where $R$ is a radius and $d$ is a distance between the astronomical bodies.
Hence, by taking the gradient of the expression (155) and integrating over the
Earth’s thin shell, we find contributions to the fifth force multiplied by
factors $R/d$.

The contribution from the lunar term $\eta_{M}$ is obtained by replacing $S$
with $M$ in the previous expressions.

For the contribution from the terrestrial term $\eta_{E}$ we find

$\int_{{\cal
V}_{E}(t)}\rho_{E}(\mathbf{x},t)\nabla\mathcal{I}_{E}(\mathbf{x},t)d^{3}x=0,$
(157) $\displaystyle\left(\frac{c^{4}}{8\pi
G}-\eta_{g}-\frac{c^{2}}{3}\int_{r_{E}}^{R_{\oplus}}\rho_{E}(r)r\,dr\right)\times$
(158) $\displaystyle\times$ $\displaystyle\int_{{\cal
V}_{E}(t)}\rho_{E}(\mathbf{x},t)\nabla\frac{r_{E}}{|\mathbf{x}-\mathbf{x}_{E}|}d^{3}x=0,$

because of spherical symmetry.

All other contributions from $\eta$ resulting from Green’s function
$G(\mathbf{x},\mathbf{x}_{0})$ up to second order are discussed in Appendix D
where it is argued that they either vanish or are multiplied by factors $R/d$.
Moreover, it turns out that all contributions to the fifth force of the order
of $R/d$ cancel each other. Since in the Sun-Earth-Moon system the ratios
$R/d$ are small, then corrections of the order of $(R/d)^{2}$ have not been
computed because exceedingly small to give rise to observable effects.

Combining all these results, we obtain the acceleration of the Earth due to
fifth force and computed by using the Green’s function up to second order:

$\displaystyle M_{\oplus}\left(\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}\right)_{\rm
f}$ $\displaystyle=$ $\displaystyle\frac{G}{3}M_{\oplus,{\rm
eff}}\left[M_{\odot,{\rm
eff}}\frac{\mathbf{x}_{S}-\mathbf{x}_{E}}{|\mathbf{x}_{S}-\mathbf{x}_{E}|^{3}}\right.$
(159) $\displaystyle+$ $\displaystyle\left.M_{M,{\rm
eff}}\frac{\mathbf{x}_{M}-\mathbf{x}_{E}}{|\mathbf{x}_{M}-\mathbf{x}_{E}|^{3}}\right].$

The acceleration of the Moon due to fifth force is analogous:

$\displaystyle M_{M}\left(\frac{d^{2}\mathbf{x}_{M}}{dt^{2}}\right)_{\rm f}$
$\displaystyle=$ $\displaystyle\frac{G}{3}M_{M,{\rm eff}}\left[M_{\odot,{\rm
eff}}\frac{\mathbf{x}_{S}-\mathbf{x}_{M}}{|\mathbf{x}_{S}-\mathbf{x}_{M}|^{3}}\right.$
(160) $\displaystyle+$ $\displaystyle\left.M_{\oplus,{\rm
eff}}\frac{\mathbf{x}_{E}-\mathbf{x}_{M}}{|\mathbf{x}_{E}-\mathbf{x}_{M}|^{3}}\right].$

In the case of astronomical bodies with uniform mass density the above
expressions coincide with the accelerations of Earth and Moon found in KW  for
chameleon gravity and in Capoz-Tsu  for $f(R)$ gravity theory. In the more
realistic case of bodies with varying density such expressions give different
results.

Since the accelerations of the Earth and Moon depend on the effective masses
of the bodies, which are the masses of the respective thin shells and these
depend on the internal structure of the bodies through density and size, then
a violation of the Weak Equivalence Principle takes place. Such a violation
can be looked for in the Earth-Moon system by means of LLR measurements.

##  VIII Acceleration of Earth and Moon due to the extra force

The acceleration due to extra force takes contributions from the extra
pressure on the surface of the moving astronomical body and from the volume
part inside the body. In the sequel we compute the acceleration of Earth,
while the acceleration of the Moon is obtained by an analogous computation.
Using Eqs. (VI.4) and (138), the contribution of the surface integral of
pressure to the acceleration is given by

$\displaystyle M_{\oplus}\left(\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}\right)_{\rm
e,s}=-\int_{\partial{\cal
V}_{E}(t)}p(R_{\oplus}^{-},\theta,\varphi,t)\hat{\mathbf{n}}_{E}d\sigma$
$\displaystyle=$
$\displaystyle\mathcal{C}(\alpha,q)\rho_{E,w}^{\frac{1}{1+|\alpha|}}\int_{\partial{\cal
V}_{E}(t)}\left(\frac{c^{2}}{8\pi
G}-\frac{\eta}{c^{2}}\right)^{\frac{\alpha}{\alpha-1}}\hat{\mathbf{n}}_{E}d\sigma,$

where we have neglected the atmospheric density with respect to the density of
seawater. Using again Eq. (VI.4) the volume part of the extra force, taking
into account that the contribution from the interior of Earth’s screening
sphere vanishes because of spherical symmetry, is given by

$\displaystyle M_{\oplus}\left(\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}\right)_{\rm
e,v}=-c^{2}\int_{{\cal V}_{E}(t)}\rho_{E}f^{2}_{R}\left[\nabla R\right]d^{3}x$
$\displaystyle=$ $\displaystyle-c^{2}\int_{\Omega_{E}}\rho_{E}f^{2}_{R}\nabla
R\,d^{3}x=-c^{2}\rho_{E,w}\int_{\Omega_{E}}\nabla f^{2}\,d^{3}x,$

where $\Omega_{E}=\\{r_{E}<|\mathbf{x}-\mathbf{x}_{E}|<R_{\oplus}\\}$ is the
thin shell of Earth, which is assumed to lie inside seawater. Using now Gauss
theorem, the formula $f^{2}(R)=qR^{\alpha}$, the expression (133) of
curvature, and taking into account that $R$ is constant on the screening
surface, we have

$\displaystyle M_{\oplus}\left(\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}\right)_{\rm
e,v}=-c^{2}\rho_{E,w}\int_{\partial{\cal
V}_{E}(t)}f^{2}(R)\hat{\mathbf{n}}_{E}d\sigma$ $\displaystyle=$
$\displaystyle-\frac{1}{\alpha}\mathcal{C}(\alpha,q)\rho_{E,w}^{\frac{1}{1+|\alpha|}}\int_{\partial{\cal
V}_{E}(t)}\left(\frac{c^{2}}{8\pi
G}-\frac{\eta}{c^{2}}\right)^{\frac{\alpha}{\alpha-1}}\hat{\mathbf{n}}_{E}d\sigma,$

so that we have the following relation between the volume and surface
contributions to Earth’s acceleration:

$\left(\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}\right)_{\rm
e,v}=-\frac{1}{\alpha}\left(\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}\right)_{\rm
e,s}.$ (164)

For the function $\eta$ inside the integral (VIII) we use the solution
computed by means of the Green’s function up to the second order: note that
such a solution was computed by neglecting curvature $R$ inside the thin shell
according to inequality (66), so that with this procedure we compute a first
correction to $R=0$. This correction satisfies itself inequality (66) for
$\alpha<-1/2$ as it is shown in Appendix E.

In order to approximate the surface integral (VIII) we introduce the small
radius-to-distance ratios $\varepsilon_{1},\varepsilon_{2},\dots$, where

$\varepsilon_{1}=\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|},\qquad\varepsilon_{2}=\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{M}|},$
(165)

the other geometric ratios are defined analogously, and the thin-shell
parameter of Earth,

$\delta_{E}=\frac{\Delta
R_{\oplus}}{R_{\oplus}}=\frac{R_{\oplus}-r_{E}}{R_{\oplus}},$ (166)

so that the thin-shell condition for Earth reads $\delta_{E}\ll 1$. If such a
condition is satisfied, then Earth’s screening surface lies inside seawater
and, at first order in $\delta_{E}$, we have the relations

$M_{\oplus,{\rm eff}}(r_{E})\approx 4\pi\rho_{E,w}R_{\oplus}^{3}\delta_{E},$
(167)

and

$\int_{r_{E}}^{R_{\oplus}}\rho_{E}(r)rdr-\frac{M_{\oplus,{\rm
eff}}(r_{E})}{4\pi R_{\oplus}}\approx\frac{M_{\oplus,{\rm eff}}(r_{E})}{8\pi
R_{\oplus}}\,\delta_{E}.$ (168)

Now we observe that the Dirichlet condition (74) requires the function
$(c^{2}/(8\pi G)-\eta/c^{2})$ to vanish on Earth’s screening surface, so that
such a function, when evaluated on the external surface of Earth, is
infinitesimal with respect to $\delta_{E}$. Such an infinitesimal property has
to be satisfied by the approximation of function $\eta$ in order to have a
valid approximation of the surface integral (VIII). Then, if we consider the
solution $\eta$ computed with the Green’s function up to the first order, we
have argued in Sec. V.5 that the sums $\mathcal{I}_{S}+\mathcal{J}_{S,E}$ and
$\mathcal{I}_{M}+\mathcal{J}_{M,E}$ vanish on Earth’s screening surface, so
that they are infinitesimal with respect to $\delta_{E}$ on $\partial{\cal
V}_{E}$. Moreover, using Eqs. (87)-(88), the radial terms of $(c^{2}/(8\pi
G)-(\eta_{E}+\eta_{g})/c^{2})$, which depend on $|\mathbf{x}-\mathbf{x}_{E}|$,
also vanish on the screening surface. The other terms in the first order
solution do not vanish on the screening surface, so that the solution computed
with the Green’s function up to the second order has to be used. By using the
method of images for a system of spheres, we show in Appendix C that the
remaining terms in the first order solution cancel on the screening surface
with terms of the second order solution, so that the sum of all such terms is
infinitesimal with respect to $\delta_{E}$ on $\partial{\cal V}_{E}$.

Using these results we make a second order Taylor approximation on the
external surface of Earth, with respect to $\delta_{E}$, and
$\varepsilon_{1},\varepsilon_{2},\dots$, of the terms in the function
$(c^{2}/(8\pi G)-\eta/c^{2})$ which vanish on Earth’s screening surface. By
the above arguments such a Taylor approximation has the overall multiplicative
factor $\delta_{E}$. By using the method of images one can check that the
remaining terms in the second order solution which do not vanish on the
screening surface, and which are not considered in the Taylor approximation,
cancel with terms of the third order solution giving rise to smaller
corrections, and this procedure can be iterated at higher orders. The leading
terms in the Taylor approximation are computed in Appendix D and they are the
following:

$\displaystyle\frac{c^{2}}{8\pi
G}-\frac{1}{c^{2}}\eta(R_{\oplus},\theta,\varphi,t)\approx\frac{\delta_{E}}{12\pi}\left[\frac{M_{\oplus,{\rm
eff}}}{2R_{\oplus}}+3\times\right.$ $\displaystyle\times$
$\displaystyle\left.\left(\varepsilon_{1}\frac{M_{\odot,{\rm
eff}}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|}\hat{\mathbf{n}}_{ES}+\varepsilon_{2}\frac{M_{M,{\rm
eff}}}{|\mathbf{x}_{E}-\mathbf{x}_{M}|}\hat{\mathbf{n}}_{EM}\right)\cdot\hat{\mathbf{n}}_{E}+\dots\right],$

where $\theta,\varphi$ are angular coordinates on Earth’s surface,
$\hat{\mathbf{n}}_{E}=\hat{\mathbf{n}}_{E}(R_{\oplus},\theta,\varphi)$, the
term multiplied by $\varepsilon_{1}$ comes from the Taylor approximation of
$\mathcal{I}_{S}+\mathcal{J}_{S,E}$, the term multiplied by $\varepsilon_{2}$
comes from the approximation of $\mathcal{I}_{M}+\mathcal{J}_{M,E}$, and the
dots represent all terms that involve the other geometric ratios. We will find
that these further terms give a negligible contribution to the extra force.

The further first order Taylor approximation with respect to
$\varepsilon_{1},\varepsilon_{2},\dots$ then follows:

$\displaystyle\left[\frac{c^{2}}{8\pi
G}-\frac{1}{c^{2}}\eta(R_{\oplus},\theta,\varphi,t)\right]^{\frac{\alpha}{\alpha-1}}\approx\left(\frac{\delta_{E}}{12\pi}\right)^{\frac{|\alpha|}{1+|\alpha|}}$
$\displaystyle\times$ $\displaystyle\left[\left(\frac{M_{\oplus,{\rm
eff}}}{2R_{\oplus}}\right)^{\frac{|\alpha|}{1+|\alpha|}}+\frac{3|\alpha|}{1+|\alpha|}\left(\frac{M_{\oplus,{\rm
eff}}}{2R_{\oplus}}\right)^{-\frac{1}{1+|\alpha|}}\right.$
$\displaystyle\times$
$\displaystyle\left.\left(\varepsilon_{1}\frac{M_{\odot,{\rm
eff}}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|}\hat{\mathbf{n}}_{ES}+\varepsilon_{2}\frac{M_{M,{\rm
eff}}}{|\mathbf{x}_{E}-\mathbf{x}_{M}|}\hat{\mathbf{n}}_{EM}\right)\cdot\hat{\mathbf{n}}_{E}+\dots\right].$

Now, using the integral equation (V.5), substituting the expression (52) of
$\eta_{g}$ into the integral equation, neglecting factors involving radius-to-
distance ratios in the integral equation and using the relation (168), the
term $\alpha q$ inside the coefficient $\mathcal{C}(\alpha,q)$ is approximated
by

$\alpha q\approx\left(\frac{8\pi
G}{c^{2}}\right)^{1+|\alpha|}\rho_{g}^{-\alpha}\frac{M_{\oplus,{\rm
eff}}}{24\pi R_{\oplus}}.$ (171)

Substituting the Taylor approximation (VIII) and the above expression of
$\alpha q$ inside the contributions (VIII) and (VIII) to the Earth
acceleration, using the relation (167) for the Earth’s effective mass,
executing the surface integrals, and adding the surface and volume parts, we
obtain the total expression of the extra force on Earth. We find that a
computation identical to the one executed for the fifth force shows that the
terms represented by dots in the Taylor approximation (VIII) give
contributions to the extra force of the order of $(R/d)^{2}$, so that they
have been neglected. The final expression of the extra force is given by

$\displaystyle M_{\oplus}\left(\frac{d^{2}\mathbf{x}_{E}}{dt^{2}}\right)_{\rm
e}\approx\left(\frac{2R_{\oplus}}{\Delta
R_{\oplus}}\right)^{\frac{1}{1+|\alpha|}}\frac{G}{3}\,\frac{r_{E}}{R_{\oplus}}\left(\frac{\rho_{g}}{\rho_{E,w}}\right)^{\frac{|\alpha|}{1+|\alpha|}}$
$\displaystyle\times$ $\displaystyle M_{\oplus,{\rm eff}}\left[M_{\odot,{\rm
eff}}\frac{\mathbf{x}_{E}-\mathbf{x}_{S}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|^{3}}+M_{M,{\rm
eff}}\frac{\mathbf{x}_{E}-\mathbf{x}_{M}}{|\mathbf{x}_{E}-\mathbf{x}_{M}|^{3}}\right].$

Comparing with Eq. (159), we see that the extra force and the fifth force are
parallel vectors which point in opposite directions. The ratio of magnitude of
the two vectors is given by

$Q_{\oplus}\approx\left(\frac{2R_{\oplus}}{\Delta
R_{\oplus}}\right)^{\frac{1}{1+|\alpha|}}\frac{r_{E}}{R_{\oplus}}\left(\frac{\rho_{g}}{\rho_{E,w}}\right)^{\frac{|\alpha|}{1+|\alpha|}}.$
(173)

The contribution of the extra force to the Moon’s acceleration is analogous:

$\displaystyle M_{M}\left(\frac{d^{2}\mathbf{x}_{M}}{dt^{2}}\right)_{\rm
e}\approx\left(\frac{2R_{M}}{\Delta
R_{M}}\right)^{\frac{1}{1+|\alpha|}}\frac{G}{3}\,\frac{r_{M}}{R_{M}}\left(\frac{\rho_{g}}{\rho_{M,c}}\right)^{\frac{|\alpha|}{1+|\alpha|}}$
$\displaystyle\times$ $\displaystyle M_{M,{\rm eff}}\left[M_{\odot,{\rm
eff}}\frac{\mathbf{x}_{M}-\mathbf{x}_{S}}{|\mathbf{x}_{M}-\mathbf{x}_{S}|^{3}}+M_{\oplus,{\rm
eff}}\frac{\mathbf{x}_{M}-\mathbf{x}_{E}}{|\mathbf{x}_{M}-\mathbf{x}_{E}|^{3}}\right],$

and the corresponding ratio of magnitude of this vector and the fifth force on
the Moon is given by

$Q_{M}\approx\left(\frac{2R_{M}}{\Delta
R_{M}}\right)^{\frac{1}{1+|\alpha|}}\frac{r_{M}}{R_{M}}\left(\frac{\rho_{g}}{\rho_{M,c}}\right)^{\frac{|\alpha|}{1+|\alpha|}}.$
(175)

If the thin-shell condition for the Moon is satisfied, then the Moon’s
screening surface lies inside the Moon’s crust so that the crustal density
$\rho_{M,c}$ appears in the above formulae.

##  IX Equivalence principle violation

The results of the previous sections imply that the Earth and Moon fall toward
the Sun with different accelerations, unless the astronomical bodies are
completely screened (the effective masses vanish), hence a violation of the
universality of free fall (UFF) takes place.

Using the expressions of the acceleration of Earth and Moon due to fifth force
and extra force, computed in Secs. VII and VIII, the leading terms of the
relative Earth-Moon acceleration are given by

$\displaystyle\mathbf{a}_{M}-\mathbf{a}_{\oplus}=-GM^{\ast}\frac{\mathbf{x}_{M}-\mathbf{x}_{E}}{|\mathbf{x}_{M}-\mathbf{x}_{E}|^{3}}+\Delta_{\rm
ESM}GM_{\odot}\frac{\mathbf{x}_{S}-\mathbf{x}_{E}}{|\mathbf{x}_{S}-\mathbf{x}_{E}|^{3}}$
$\displaystyle+\left(1+\Delta_{\rm
tidal}\right)GM_{\odot}\left(\frac{\mathbf{x}_{E}-\mathbf{x}_{S}}{|\mathbf{x}_{E}-\mathbf{x}_{S}|^{3}}-\frac{\mathbf{x}_{M}-\mathbf{x}_{S}}{|\mathbf{x}_{M}-\mathbf{x}_{S}|^{3}}\right),$

where

$\displaystyle\Delta_{\rm ESM}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\frac{M_{\odot,{\rm
eff}}}{M_{\odot}}\left[(1-Q_{M})\frac{M_{M,{\rm
eff}}}{M_{M}}-(1-Q_{\oplus})\frac{M_{\oplus,{\rm eff}}}{M_{\oplus}}\right],$
$\displaystyle M^{\ast}$ $\displaystyle=$ $\displaystyle M_{\oplus}+M_{M}$
$\displaystyle+$ $\displaystyle\frac{1}{3}\frac{M_{\oplus,{\rm
eff}}}{M_{\oplus}}\frac{M_{M,{\rm
eff}}}{M_{M}}\left[M_{\oplus}(1-Q_{M})+M_{M}(1-Q_{\oplus})\right],$
$\displaystyle\Delta_{\rm tidal}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(1-Q_{M})\frac{M_{M,{\rm
eff}}}{M_{M}}\frac{M_{\odot,{\rm eff}}}{M_{\odot}}.$ (179)

The meaning of the terms in the expression of
$\mathbf{a}_{M}-\mathbf{a}_{\oplus}$ is the following Wi :

  * (i)

the first term is the relative acceleration due to the gravitational
attraction between the Earth and Moon;

  * (ii)

the second term, which can be written in the form,

$\Delta_{\rm ESM}\,\mathbf{g}_{S},$ (180)

$\mathbf{g}_{S}$ being the Newtonian acceleration of Earth due to the Sun, is
the UFF violation-related difference between the Earth and the Moon
accelerations toward the Sun, hence, in the framework of NMC gravity, this
term gives rise to a violation of the Weak Equivalence Principle (WEP);

  * (iii)

the third term is the solar tidal perturbation of the Moon’s orbit,
$\Delta_{\rm tidal}$ being the NMC gravity correction to the Newtonian
perturbation.

The size of the UFF violation is represented by the parameter $\Delta_{\rm
ESM}$, where ESM stands for Sun, Earth and Moon, since, by definition of the
effective mass, such a parameter depends on the composition (density) and size
of all the three astronomical bodies Visw-Fienga . Particularly, the WEP
violation depends on the composition of the Sun, in addition to the more usual
dependence on the compositions of the Earth and the Moon.

If the astronomical bodies are completely screened (screening radii equal to
radii of the bodies), then, neglecting the contribution to density of the
solar and terrestrial atmospheres, the effective masses vanish and
$\Delta_{\rm ESM}=0$, so that WEP is satisfied.

If we denote $d_{\rm EM}$ the Earth-Moon distance, then, the UFF acceleration
(180) gives rise to a polarization of the Moon’s orbit in the direction of the
Sun with a periodic perturbation $\Delta d_{\rm EM}$ of the form Nordtvedt-1

$\Delta d_{\rm EM}=S\Delta_{\rm ESM}\cos D,$ (181)

where $D$ is the synodic phase measured from the new Moon, and $S$ is a
scaling factor whose theoretical computation gives $S=2.9\times 10^{12}$ cm
Nordtvedt-2 .

The Newtonian solar tidal perturbation of the Moon’s orbit also produces a
periodic perturbation of $d_{\rm EM}$ that has a component at the synodic
frequency; nevertheless, such a perturbation can be accounted for to a very
small uncertainty Williams09 , so that the corresponding NMC gravity
correction can be neglected for $|\Delta_{\rm tidal}|\ll 1$.

In the next section constraints on the NMC gravity model will be obtained by
means of the test of WEP performed in Ref. Visw-Fienga  using lunar LLR data.

##  X LLR constraints on NMC gravity parameters

In Ref. Visw-Fienga  the authors give a general constraint in terms of
difference between the Earth and the Moon accelerations toward the Sun,
without assuming metric theories or other types of modified gravity theories.
In order to test UFF violations, a supplementary acceleration of the form
(180) is introduced in the geocentric equation of motion of the Moon. The
parameter $\Delta_{\rm ESM}$ is estimated in the LLR adjustment together with
a set of parameters of the lunar ephemerides listed in Visw-Fienga .

Then, using the expression (LABEL:Delta-ESM), specific of our framework, the
estimate of $\Delta_{\rm ESM}$ obtained in Visw-Fienga  can be directly
translated into a constraint on NMC gravity parameters. The result on the WEP
violation parameter in Visw-Fienga , based on 48 years of LLR data, is given
by

$\Delta_{\rm ESM}=\left(3.8\pm 7.1\right)\times 10^{-14}.$ (182)

In our framework $\Delta_{\rm ESM}$ is a function of $\alpha$ and of the
screening radii of the Sun, Earth and Moon,

$\Delta_{\rm ESM}=\Delta_{\rm ESM}(\alpha,r_{S},r_{E},r_{M}),$ (183)

which are determined by Eqs. (V.5) as functions of the parameters $\alpha$ and
$q$ of the NMC gravity model:

$r_{S}=r_{S}(\alpha,q),\quad r_{E}=r_{E}(\alpha,q),\quad
r_{M}=r_{M}(\alpha,q).$ (184)

Substituting such functions in the constraint (182) we obtain the set of the
admissible values of parameters $\alpha$ and $q$.

In order to avoid either too small or too large numbers we replace parameter
$q$ with the following rescaled, dimensionless parameter:

$\widetilde{q}=qR_{g}^{\alpha},$ (185)

with $R_{g}=8\pi G\rho_{g}/c^{2}$. With this substitution, the function
$f^{2}(R)$ can be written in the form

$f^{2}(R)=\widetilde{q}\left(\frac{R}{R_{g}}\right)^{\alpha}.$ (186)

We represent the constraints from LLR by means of a two-dimensional exclusion
plot in the plane with coordinates $\alpha,q$. The resulting admissible region
in the plane is restricted by means of the condition $\lambda_{g}\gg r_{g}$
introduced in Section IV.4 which, expressed in terms of parameters
$\alpha,\widetilde{q}$ and using Eq. (32) with $\rho=\rho_{g}$, becomes

$\left[\frac{3}{4\pi}\alpha(1-\alpha)\frac{c^{2}}{G\rho_{g}}\widetilde{q}\right]^{1/2}>10^{2}r_{g},$
(187)

where we have required $\lambda_{g}>10^{2}r_{g}$.

The LLR constraints are graphically reported in Figures 2 and 2: admissible
regions for parameters are plotted in white, while the excluded regions are
plotted in grey.

Figure 1: LLR constraints on the parameter quarter plane
$|\alpha|,|\widetilde{q}|$ for $0>\alpha>-10$. The solid black line yields the
upper bound on $|\widetilde{q}|$ from LLR data, the dashed black line yields
the lower bound on $|\widetilde{q}|$ from inequality
$\lambda_{g}>10^{2}r_{g}$. The dashed red line represents the condition Sun’s
screening radius $r_{S}$ equal to the radius $r_{p}$ at the base of the solar
photosphere.

Figure 2: LLR constraints on the parameter quarter plane
$|\alpha|,|\widetilde{q}|$ for $0>\alpha>-1$. The solid black line yields the
upper bound on $|\widetilde{q}|$ from LLR data. The meaning of the dashed
black line and of the dashed red line is the same as in Fig. 2.

Fig. 2 shows the admissible region for values of parameter $\alpha$ in the
range $(-10,0)$; since both $\alpha$ and $\widetilde{q}$ are negative the
admissible region is plotted in the quarter plane with coordinates
$(|\alpha|,|\widetilde{q}|)$. Fig. 2 shows the admissible region for values of
$\alpha$ in the range $(-1,0)$. The portion of the admissible region which
lies above the dashed red line $r_{S}=r_{p}$ corresponds to values of $r_{S}$
such that the Sun’s screening sphere lies in the solar convection zone.

In Ref. MBMGDeA  it has been found that the Cassini measurement of PPN
parameter $\gamma$ constrains the parameters $\alpha,\widetilde{q}$ to be of
the order $|\widetilde{q}|<10^{-12}$ for $-1>\alpha>-10$. Hence, the
constraints from WEP violation and LLR data provide tighter bounds on model
parameters compared to bounds from Cassini measurement.

Figures 4 and 4 show the thin-shell parameter of Earth $\delta_{E}$ and the
corresponding parameters $\delta_{S}$ and $\delta_{M}$ of Sun and Moon, which
are defined analogously to formula (166), plotted versus $|\widetilde{q}|$ for
a fixed value of $\alpha$. Fig. 4 shows the plot for $\alpha=-1$ and Fig. 4
for $\alpha=-10$.

Figure 3: Plot of thin-shell parameters $\delta_{S},\delta_{M}$ and
$\delta_{E}$ versus $|\widetilde{q}|$ for $\alpha=-1$. The vertical solid
black line on the right corresponds to the upper bound on $|\widetilde{q}|$
for $\alpha=-1$. The vertical dashed black line corresponds to the condition
$\lambda_{g}=10^{2}r_{g}$ and the red one corresponds to $r_{S}=r_{p}$. The
horizontal dashed and dotted lines correspond to $r_{S}=r_{p}$ and
$r_{S}=r_{\rm conv}$, respectively.

Figure 4: Plot of thin-shell parameters $\delta_{S},\delta_{M}$ and
$\delta_{E}$ versus $|\widetilde{q}|$ for $\alpha=-10$. The vertical solid
black line on the right corresponds to the upper bound on $|\widetilde{q}|$
for $\alpha=-10$. The meaning of the vertical dashed black and red lines, and
of the horizontal dashed and dotted lines, is the same as in Fig 4.

The excluded regions are colored in grey. The figures show that increasing
$|\widetilde{q}|$ the Sun’s screening radius enters into the convection zone.

##  XI Conclusions

In this work we have examined the NMC theories of gravity specified by Eq.
(11). These theories for certain values of $q$ and $\alpha$ are potentially
inconsistent with the Equivalence Principle given that the fifth force and the
extra non-Newtonian force they give rise could imply that Earth and Moon fall
differently towards the Sun. We have shown that the screening chameleon-type
mechanism, discussed for the first time in the context of these theories in
Ref. MBMGDeA , can be used to fully shield any violation of the Equivalnce
Principle provided the astronomical bodies are completely screened, that is,
their screening radii are equal to the radii of the bodies as discussed in
section IX. Furthermore, we have shown in section X for which range of values
of the parameters $\alpha$ and $q$, compatibility with 48 years of LLR data
can be ensured. Hence, the screening mechanism here discussed does allow for
compatibility with data even of versions of the NMC gravity theories that are
potentially problematic. We find that constraints from WEP compatible with the
LLR data yield tighter bounds on the model parameters compared than the
corresponding bounds from the Cassini measurement of the parametrized post-
Newtonian parameter, $\gamma$. Indeed, WEP and LLR data, compatible with lunar
ephemerides, do provide stringent bounds on the parameters $q$ and $\alpha$,
more specifically, as depicted in Figs. (2) and (2), for $|\alpha|$ and the
rescaled value $|\widetilde{q}|$ given by Eq. (185). Compatibility with data
is found for $|\widetilde{q}|<10^{-13}$ and for $0.2<|\alpha|<10$. As for the
robustness of the screening mechanism, the thin-shell conditions were tested
for different model parameters for Earth, Moon and Sun, as shown in Figs. (4)
and (4). As can one inspect, consistency is found for the range of model
parameters that are compatible with data.

As a final remark, we can state, on general grounds, that our results provide
yet another confirmation that the NMC gravity model under study remains a
workable alternative to address certain issues for which GR is not fully
adequate. Furthermore, given that the particular subclass of NMC model under
study requires the development of various techniques for the implementation of
the screening mechanism, its relevance goes beyond the specificities of the
model as the technical issues addressed in our work can, in principle, be
useful in any model of gravity.

## Acknowledgments

The work of R.M. is partially supported, and the work of M.M. and S.DA is
fully supported, by INFN (Istituto Nazionale di Fisica Nucleare, Italy), as
part of the MoonLIGHT-2 experiment in the framework of the research activities
of the Commissione Scientifica Nazionale n. 2 (CSN2). S.DA and M.M. would like
to thank ASI (Agenzia Spaziale Italiana) for the support granted to them with
ASI-INFN Agreement n. 2019-15-HH.0.

O.B. acknowledges the partial support from Fundação para a Ciência e a
Tecnologia (Portugal) through the research project CERN/FIS-PAR/0027/2021.

## Appendix A: interior density profiles

We report on models of mass density profiles for the Sun, Earth and Moon which
have been used in order to find analytical estimates of the constraints on the
NMC gravity parameters at a suitable order of magnitude.

###  .1 Sun density profile

The complete model of mass density of the Sun is reported in Ref. MBMGDeA . In
this appendix we report the density profiles of the convection zone and
photosphere since the screening radius $r_{S}$ that saturates the LLR bound
lies in such zones. For the details see Ref. MBMGDeA .

The radius of the Sun is $R_{\odot}=6.9634\times 10^{5}\,{\rm km}$, the solar
atmosphere begins below the spherical surface of radius $R_{\odot}$ and center
in the origin, at a depth of about 500 km, and extends outward from the Sun.
Then $r_{p}=R_{\odot}-500\,{\rm km}$ is the radius at the base of the
photosphere. Matter in the Sun is modeled as a perfect gas in hydrostatic
equilibrium.

Convection zone. We use a polytrope model with an effective polytropic index
$n_{c}=2.33$ Mullan :

$\rho_{S}(r)=K_{c}\left[T_{S}(r)\right]^{n_{c}},$ (188)

with $K_{c}=3.44\times 10^{-16}$, and the radius $r$ varying in the range

$r_{\rm conv}\leq r<r_{p},\qquad r_{\rm conv}=5.3185\times 10^{5}\,\mbox{km}.$
(189)

The temperature profile is approximated by Mullan

$T_{S}(r)=\frac{GM_{\odot}}{C_{p}r}-T_{0},$ (190)

with $M_{\odot}=1.989\times 10^{33}\,{\rm g}$, $C_{p}=2.95\times 10^{8}\,{\rm
erg\,g^{-1}\,K^{-1}}$ is an averaged value of the specific heat at constant
pressure, and $T_{0}=6.461\times 10^{6}\,{\rm K}$.

Photosphere. For the density in the photosphere we use the following model,
adapted from DeJager , for $r_{p}\leq r<R_{\odot}$:

$\rho_{S}(r)=\frac{\mu
m_{p}p_{m}}{k_{B}[T_{m}+A(R_{\odot}-r)^{2}])}\exp\left(\frac{R_{\odot}-r}{H_{p}}\right),$
(191)

where $\mu=1.26$, $m_{p}=1.66\times 10^{-24}\,{\rm g}$ is the proton mass,
$k_{B}=1.3806\times 10^{-16}\,{\rm erg\,K^{-1}}$ is the Boltzmann constant,
$H_{p}=117\,{\rm km}$, $A=8.8\times 10^{-3}\,{\rm km^{-2}K}$, $T_{m}=4.4\times
10^{3}\,{\rm K}$ is the temperature minimum at the top of the photosphere,
$p_{m}$ is pressure corresponding at the temperature minimum, such that
$p_{m}/k_{B}=1.2\times 10^{19}\,{\rm K\,cm^{-3}}$.

The Sun density profile permits us to compute the values of the length
$\lambda(\rho_{S})$. By using the definition (32) of $\lambda$ we have

$\lambda(\rho_{S})=\lambda_{g}\left(\frac{\rho_{S}}{\rho_{g}}\right)^{\frac{\alpha-1}{2}},$
(192)

and by the discussion in Sec. IV.4 we require $\lambda_{g}\gg r_{g}$. For
instance, if $\alpha=-1$ and $\lambda_{g}=10^{3}r_{g}\sim 10^{5}$ AU, at the
bottom of the solar convection zone where $\rho_{S}\approx 1.65\times
10^{-1}\,{\rm g}/{\rm cm}^{3}$, we have $\lambda(\rho_{S})\approx
10^{-6}\,{\rm m}$. At the top of the convection zone where $\rho_{S}\approx
2.73\times 10^{-7}\,{\rm g}/{\rm cm}^{3}$, we have $\lambda(\rho_{S})\approx
7\times 10^{-1}\,{\rm m}$.

###  .2 Earth density profile

We consider an average Earth model in the sense of PREM . The Earth is divided
into four regions with constant mass density separated by spherical surfaces
of density discontinuities: ocean layer, crust, mantle and core. The following
numerical values of radii and densities are taken from PREM ; Stacey . The
radius of the Earth is $R_{\oplus}=6371$ km.

Ocean layer. Since the ocean and seas cover $70.8\%$ of the surface of the
Earth, the uppermost layer of the average Earth model consists of seawater
with a depth of 3 km. The layer corresponds to distances $r$ from the Earth’s
center $R_{E,c}<r\leq R_{\oplus}$:

$R_{E,c}=6368\,\,{\rm km},\qquad\rho_{E,w}=1.02\,\,{\rm g\,cm^{-3}}.$ (193)

Crust. The layer corresponds to $R_{E,m}<r\leq R_{E,c}$:

$R_{E,m}=6346.6\,\,{\rm km},\qquad\rho_{E,c}=2.7\,\,{\rm g\,cm^{-3}}.$ (194)

At radius $R_{E,m}$ there is the Mohorovičić discontinuity.

Mantle. The layer corresponds to $R_{E,n}<r\leq R_{E,m}$:

$R_{E,n}=3480\,\,{\rm km},\qquad\rho_{E,m}=4.5\,\,{\rm g\,cm^{-3}}.$ (195)

At radius $R_{E,n}$ there is the Gutenberg discontinuity.

Core. The layer corresponds to $0<r\leq R_{E,n}$:

$\rho_{E,n}=11\,\,{\rm g\,cm^{-3}}.$ (196)

Using Eq. (192) we may compute $\lambda(\rho_{E})$ in the various Earth’s
layers. We have $\lambda(\rho_{E})\leq\lambda(\rho_{E,w})$, then, using for
instance $\alpha=-1$ and $\lambda_{g}=10^{5}$ AU, we find
$\lambda(\rho_{E,w})\approx 10^{-7}\,{\rm m}$, hence a completely negligible
quantity.

###  .3 Moon density profile

The Moon is divided into three regions with constant mass density separated by
spherical surfaces of density discontinuities: crust, mantle and core. The
following numerical values of radii and densities are taken from Garcia ;
Wieczor ; Matsumoto . The radius of the Moon is $R_{M}=1737$ km.

Crust. The layer corresponds to distances $r$ from the Moon’s center
$R_{M,m}<r\leq R_{M}$ Wieczor :

$R_{M,m}=1697\,\,{\rm km},\qquad\rho_{M,c}=2.55\,\,{\rm g\,cm^{-3}}.$ (197)

At radius $R_{M,m}$ there is the lunar Moho.

Mantle. The layer corresponds to $R_{M,n}<r\leq R_{M,m}$:

$\rho_{M,m}=3.4\,\,{\rm g\,cm^{-3}}.$ (198)

Core. The layer corresponds to $0<r\leq R_{M,n}$ Garcia ; Matsumoto :

$R_{M,n}=330\,\mbox{--}\,400\,\,{\rm
km},\qquad\rho_{M,n}=3.9\,\mbox{--}\,5.5\,\,{\rm g\,cm^{-3}}.$ (199)

Considerable uncertainty is connected with the radius and physical state of a
metallic core.

In the Moon’s layers, using again $\alpha=-1$ and $\lambda_{g}=10^{5}$ AU, we
have $\lambda(\rho_{M})\lesssim 10^{-7}\,{\rm m}$.

## Appendix B: second order Green’s function

First we give the six second order images inside the screening spheres. The
position vectors of the two image points inside the Sun’s screening sphere are
given by

$\displaystyle\widehat{\mathbf{x}}_{SE}=\mathbf{x}_{S}+r_{S}^{2}\frac{\widetilde{\mathbf{x}}_{E}-\mathbf{x}_{S}}{|\widetilde{\mathbf{x}}_{E}-\mathbf{x}_{S}|^{2}},$
$\displaystyle\widehat{\mathbf{x}}_{SM}=\mathbf{x}_{S}+r_{S}^{2}\frac{\widetilde{\mathbf{x}}_{M}-\mathbf{x}_{S}}{|\widetilde{\mathbf{x}}_{M}-\mathbf{x}_{S}|^{2}},$
(200)

where $\widehat{\mathbf{x}}_{SE},\widehat{\mathbf{x}}_{SM}$ are the images of
the sources at $\widetilde{\mathbf{x}}_{E},\widetilde{\mathbf{x}}_{M}$,
respectively. The other four second order images are obtained by changing the
subscripts $S,E,M$.

The second order term of Green’s function contains six terms and its
expression is then given by

$\displaystyle G^{(2)}(\mathbf{x},\mathbf{x}_{0})$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\left[\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{0}|}\,\frac{r_{S}}{|\widetilde{\mathbf{x}}_{E}-\mathbf{x}_{S}|}\,\frac{1}{|\mathbf{x}-\widehat{\mathbf{x}}_{SE}|}\right.$
$\displaystyle+$
$\displaystyle\frac{r_{M}}{|\mathbf{x}_{M}-\mathbf{x}_{0}|}\,\frac{r_{S}}{|\widetilde{\mathbf{x}}_{M}-\mathbf{x}_{S}|}\,\frac{1}{|\mathbf{x}-\widehat{\mathbf{x}}_{SM}|}$
$\displaystyle+$
$\displaystyle\frac{r_{S}}{|\mathbf{x}_{S}-\mathbf{x}_{0}|}\,\frac{r_{E}}{|\widetilde{\mathbf{x}}_{S}-\mathbf{x}_{E}|}\,\frac{1}{|\mathbf{x}-\widehat{\mathbf{x}}_{ES}|}$
$\displaystyle+$
$\displaystyle\frac{r_{M}}{|\mathbf{x}_{M}-\mathbf{x}_{0}|}\,\frac{r_{E}}{|\widetilde{\mathbf{x}}_{M}-\mathbf{x}_{E}|}\,\frac{1}{|\mathbf{x}-\widehat{\mathbf{x}}_{EM}|}$
$\displaystyle+$
$\displaystyle\frac{r_{S}}{|\mathbf{x}_{S}-\mathbf{x}_{0}|}\,\frac{r_{M}}{|\widetilde{\mathbf{x}}_{S}-\mathbf{x}_{M}|}\,\frac{1}{|\mathbf{x}-\widehat{\mathbf{x}}_{MS}|}$
$\displaystyle+$
$\displaystyle\left.\frac{r_{E}}{|\mathbf{x}_{E}-\mathbf{x}_{0}|}\,\frac{r_{M}}{|\widetilde{\mathbf{x}}_{E}-\mathbf{x}_{M}|}\,\frac{1}{|\mathbf{x}-\widehat{\mathbf{x}}_{ME}|}\right].$

We now give the contribution $\eta_{E}^{(2)}(\mathbf{x},t)$ from the Earth’s
thin shell to the solution $\eta(\mathbf{x},t)$ obtained by adding the second
order term of Green’s function. The following expression has to be added to
formula (88). We drop the dependence on $(\mathbf{x},t)$ for simplicity:

$\displaystyle\eta_{E}^{(2)}$ $\displaystyle=$
$\displaystyle-\left(\frac{c^{4}}{8\pi
G}-\eta_{g}\right)r_{E}\left[\frac{r_{S}}{\left||\mathbf{x}-\mathbf{x}_{S}|(\mathbf{x}_{S}-\mathbf{x}_{E})+r_{S}^{2}\hat{\mathbf{n}}_{S}\right|}+\frac{r_{M}}{\left||\mathbf{x}-\mathbf{x}_{M}|(\mathbf{x}_{M}-\mathbf{x}_{E})+r_{M}^{2}\hat{\mathbf{n}}_{M}\right|}\right]$
(202) $\displaystyle+$
$\displaystyle\mathcal{F}_{E,ES}+\mathcal{F}_{E,SE}+\mathcal{F}_{E,EM}+\mathcal{F}_{E,ME}+\mathcal{F}_{E,SM}+\mathcal{F}_{E,MS},$
$\mathcal{F}_{E,ES}+\mathcal{F}_{E,SE}=\frac{c^{2}}{3}\,r_{E}r_{S}\left\\{\frac{\int_{r_{E}}^{R_{\oplus}}\rho_{E}(r)r\,dr}{\left||\mathbf{x}-\mathbf{x}_{S}|(\mathbf{x}_{S}-\mathbf{x}_{E})+r_{S}^{2}\hat{\mathbf{n}}_{S}\right|}+\frac{M_{\oplus,{\rm
eff}}(r_{E})/4\pi}{\left|r_{E}^{2}(\mathbf{x}_{E}-\mathbf{x}_{S})+|\mathbf{x}-\mathbf{x}_{E}|\left(|\mathbf{x}_{E}-\mathbf{x}_{S}|^{2}-r_{S}^{2}\right)\hat{\mathbf{n}}_{E}\right|}\right\\},$
(203) $\displaystyle\mathcal{F}_{E,SM}+\mathcal{F}_{E,MS}$ $\displaystyle=$
$\displaystyle\frac{c^{2}}{12\pi}\,r_{M}r_{S}\left(A_{SM}+A_{MS}\right)M_{\oplus,{\rm
eff}}(r_{E}),$ $\displaystyle A_{SM}$ $\displaystyle=$