# What to expect from scalar-tensor space geodesy

Hugo Lévy, Joël Bergé, Jean-Philippe Uzan

► To cite this version:

Hugo Lévy, Joël Bergé, Jean-Philippe Uzan. What to expect from scalar-tensor space geodesy. Physical Review D, 2024, 109 (8), pp.084009. 10.1103/PhysRevD.109.084009. hal-04256905

HAL Id: hal-04256905

https://hal.science/hal-04256905v1

Submitted on 17 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
What to expect from scalar-tensor space geodesy

Hugo Lévy,<fnref n="1" />,²,\* Joël Bergé,<fnref n="1" /> and Jean-Philippe Uzan²

¹DPHY, ONERA, Université Paris Saclay F-92322 Châtillon - France

²Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France

Scalar-tensor theories with screening mechanisms come with non-linearities that make it difficult to study setups of complex geometry without resorting to numerical simulations. In this article, we use the *femtoscope* code that we introduced in a previous work in order to compute the fifth force arising in the chameleon model in the Earth orbit. We go beyond published works by introducing a departure from spherical symmetry — embodied by a mountain on an otherwise spherical Earth — as well as by implementing several atmospheric models, and quantify their combined effect on the chameleon field. Building on the numerical results thus obtained, we address the question of the detectability of a putative chameleon fifth force by means of space geodesy techniques and, for the first time, quantitatively assess the back-reaction created by the screening of a satellite itself. We find that although the fifth force has a supposedly measurable effect on the dynamics of an orbiting spacecraft, the imprecise knowledge of the mass distribution inside the Earth greatly curtails the constraining power of such space missions. Finally, we show how this degeneracy can be lifted when several measurements are performed at different altitudes.

I. Introduction

Scalar fields appear in most of the extensions beyond the standard models. Theories involving extra dimensions, from Kaluza-Klein theories up to string theories in the low energy limit, predict the existence of a light spin-0 particle. Scalar fields are also key ingredients in cosmology phenomenology, in particular for the dark sector and inflation. Coupling the scalar field to matter<fnref n="1" /> automatically gives rise to a so-called *fifth force*, resulting in deviations from general relativity (GR) in gravitational phenomena. Evading the Solar system tests of GR and laboratory experiments [2] comes at the price of introducing non-linearities in the model which enable *screening mechanisms* (e.g. Damour-Polyakov [3], chameleon [4, 5], K-mouflage [6, 7], or Vainshtein [8, 9]).

Although screening mechanisms are precisely designed to recover GR — and thus in the weak field regime, Newtonian gravity — at Solar system scales, they leave nonetheless a small imprint which we can attempt to measure. Tests can be performed in a very wide range of length scales, from laboratory experiments [10–12], to spacecraft in orbit around the Earth [13–15] or traveling through the Solar system [16], planetary motion [17–19], and to astrophysical tests [20–22] (see Refs. [1, 23] and references therein for a more comprehensive review). Here, we are interested in the category of space-based experiments, which have long been expected to provide new constraints in the case of the chameleon model [5]. Several space missions were successfully launched in the past decades: MICROSCOPE [13, 24] for testing the weak equivalence principle (see Refs. [25, 26] for how constraints on the chameleon model could be derived from

those data), Gravity Probe A and B [27], LAGEOS 1 and 2, LARES 1 and 2 [14].

Beside these space missions specifically tailored for fundamental physics, artificial satellites have also given rise to space geodesy. Initially, space geodesy primarily focused on measuring the Earth's shape and size, but technological advancements have propelled it into a realm of unprecedented accuracy and multifaceted applications. Cutting edge instruments onboard satellites allow for the implementation of complementary geodetic techniques such as laser and Doppler ranging, Global Navigation Satellite Systems, gravimetry (e.g. GOCE, CHAMP, GRACE-FO satellite missions), etc. The determination of the Earth's figure (mass distribution) constitutes an inverse problem: given the data $d_{\text{obs}}$ collected by the various satellite missions and a model describing the laws of gravitation $M$ with forward map $F_M$, the goal is to determine the model parameters $p$ such that the residual $d_{\text{obs}} - F_M(p)$ is minimized (in some specific sense, e.g. least-squares or probabilistic approaches). In space geodesy, this inverse problem is solved with the central assumption that the governing equation is Newton's law of gravity (and $p$ would represent the distribution of mass) [28].

The goal of the present article is to assess the pertinence of orbitography techniques to test screened scalar-tensor theories, illustrated with the chameleon model, and to characterize the *best site* in the Solar system to perform such tests. This is a follow-up to our previous article Ref. [29] where we laid the foundations in terms of numerical simulations. There, we saw that the unconstrained region of the chameleon parameter space (see Fig. 3 of Ref. [30]) corresponds to a situation where the Earth is screened, i.e. where the chameleonic force is sourced only by its outer layers. This mere observation suggests that the local landform — specifically any local deviation from spherical symmetry — can leave a significant imprint on the chameleon profile. Consequently, if the chameleon's effects differ sufficiently from Newton-

<footnote>\* hugo.levy@onera.fr</footnote>

<footnote n="1">¹ From a quantum mechanical perspective, the introduction of a scalar field in the gravity sector *always* generate interactions between this scalar and matter fields [1].</footnote>
nian gravity, it should leave a distinctive signature on the Earth's gravity.

Mountains and craters are typical examples of as-
phericities that can be sensed through space geodesy.
Relative to the size of a planet, a mountain represents
a spiky feature. Several works bring to light the paral-
lel between chameleon (and symmetron) gravity in the
screened regime and electrostatics: the behavior of the
scalar field is roughly the same as the behavior of the
electrostatic potential for a perfect conductor<fnref n="2" /> [31–33].
Taking the analogy a step further, Ref. [31] mentions the
“lightning rod effect” in electromagnetism, exhibited by
needle-like conductors around which the electric field (∞
gradient of the potential) is enhanced. In the case of the
chameleon, the counterpart of the electric field would be
the fifth force (∞ gradient of the scalar field) — mak-
ing the mountain an interesting case study. Neverthe-
less as Ref. [33] underlines, while this analogy provides
valuable qualitative insights, numerical computations re-
main essential to establish a quantitative connection with
real-world observations and experimental data. In that
respect, we aim to address the long-standing question of
how much an atmosphere smooths out the mountain’s
contribution to the fifth force in space. More generally,
existing work accounting for the atmosphere [4, 5, 34–
37] are, in our opinion, not extensive enough: the models
are not realistic (one layer of constant density) and con-
clusions are drawn on qualitative arguments that can be
misleading (see e.g. the introduction of Ref. [38]). We
shall also pay attention to the influence of a spacecraft on
the background field, and evaluate how this perturbation
impacts the overall fifth force that it experiences.

The article is organized as follows. In Sec. II, we briefly recall the main equations describing both Newtonian and chameleon gravity, and give precise meaning to physical models outlined above, namely the modeling of the mountainous planet together with its atmosphere. In this setup, the total gravitational potential is computed numerically using femtoscope, a code that was specifically designed to solve these equations with asymptotic boundary conditions [29]. It allows for the computation of both the Newtonian potential and the chameleon field in space. The numerical results are presented and discussed in Sec. III. We explore a vast region of the chameleon parameter space and ascertain the influence of an atmosphere in several scenarios, making this a quite comprehensive study compared to what has been done in previous work. Finally, Sec. IV takes us back to space geodesy as we compare the dynamics of a spacecraft with and without a fifth force acting on it as it orbits the mountainous planet. We address the issue of being

able to discriminate between the two in the presence of
model uncertainties, and further suggest ways to break
this source of degeneracy. These analyses pave the way
to the design of orbitography experiments in the Solar
system and their subtle interpretation. We conclude in
Sec. V.

II. Model & Numerical techniques

A. General equations

1. Newtonian gravity

It is well known that, in the weak-field regime and when the sources are moving very slowly compared to the speed of light, GR reduces to Newtonian gravity which is described by the *Newtonian potential* Φ with dimension [*L*<fnref n="2" />·*T*⁻<fnref n="2" />]. For a static configuration, we define it as

$$
\Phi(\mathbf{x}) = -G \int_{\mathbb{R}^3} \frac{\rho(\mathbf{x}')}{\|\mathbf{x} - \mathbf{x}'\|} d^3\mathbf{x}', \quad (1)
$$

where *G* is the Newtonian gravitational constant and
*ρ* = *ρ*(*x*) is the matter density function which depends
on position *x*. Assuming that the *weak equivalence principle*
holds perfectly (Ref. [13] shows that it holds at
less than 10<sup>-15</sup>) and from a classical mechanics perspec-
tive, the gravitational acceleration undergone by a point-
like particle is simply **a**<sub>Φ</sub> = −∇Φ. Eq. (1) provides a
straightforward way of computing the Newtonian poten-
tial by evaluating some three-dimensional integral (see
e.g. Ref. [39]). However, it may be more convenient
from a numerical standpoint to solve the following Pois-
son’s equation

$$
\Delta \Phi = 4\pi G \rho, \quad (2)
$$

implied by the definition of Φ. Indeed, on the one hand
one has to evaluate the integral appearing in Eq. (1) for
each point **x** where the Newtonian potential is sought,
whereas on the other hand solving the partial differential
equation (PDE) (2) provides an approximation of Φ over
the whole numerical domain.

Assuming that the mass density vanishes as one moves
away from the source of gravity, the gravitational acceler-
ation **a**<sub>Φ</sub> = −∇Φ is expected to decay to zero at infinity.
The essential boundary condition is therefore defined at
infinity and a very common choice for the constant of
integration is

$$
\Phi(\mathbf{x}) \underset{\|\mathbf{x}\| \to +\infty}{\longrightarrow} 0. \qquad (3)
$$

2. Chameleon gravity

In the Newtonian limit, the chameleon field $\phi$ is gov-
erned by a nonlinear Klein-Gordon equation which takes

<footnote n="2">² Indeed, it can be shown that the equation of motion of the chameleon field in the quasi-static Newtonian limit with thin-shell can be well-approximated by the electrostatic potential equation. Then, same differential equations lead to same solutions.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Mapping from the chameleon parameter space in the plane $n = 1$ to the dimensionless parameter $\alpha$ appearing in Eq. (7), where $M = M_{Pl}/\beta$. The gray lines represent the iso-values of the $\alpha$ parameter covered in this study, ranging from $10^{-5}$ to $10^{-28}$. The orange horizontal dashed line corresponds to $\Lambda = \Lambda_{DE} = 2.4 \times 10^{-3}$ eV, the dark energy scale.</figcaption></figure>

the form

$$ \Delta\phi = \frac{dV_{\text{eff}}}{d\phi} = \frac{\beta}{M_{\text{Pl}}} \rho - \frac{n\Lambda^{n+4}}{\phi^{n+1}}, \quad (4) $$

where $M_{Pl} \equiv 1/\sqrt{8\pi G}$ is the reduced Planck mass and $V_{\text{eff}}$ is the so-called *effective potential* of the scalar field. The model further has three parameters — $\beta$ a positive dimensionless constant which encodes the coupling of the scalar field to matter, $\Lambda$ a mass scale and $n$ a natural number. The 3-acceleration experienced by a point-like particle induced by its coupling to the chameleon field is proportional to the gradient of the scalar field and takes the form

$$ \mathbf{a}_{\phi} = -\frac{\beta}{M_{\text{Pl}}} \nabla \phi. \quad (5) $$

If we assume that the density uniformly decays to some vacuum density $\rho_{\text{vac}}$ far away from the source, then the chameleonic acceleration is expected to decay to zero at infinity, just as in the Newtonian gravity case discussed above. Equating the r.h.s. of Eq. (4) to zero and solving for $\phi$ yields the following asymptotic boundary condition:

$$ \phi(\mathbf{x}) \underset{\|\mathbf{x}\| \to +\infty}{\longrightarrow} \left( M_{\text{Pl}} \frac{n\Lambda^{n+4}}{\beta\rho_{\text{vac}}} \right)^{\frac{1}{n+1}} \equiv \phi_{\text{vac}}. \quad (6) $$

In Ref. [29], we introduced *femtoscope* — a PYTHON numerical tool based on the finite element method which enables us to solve Eq. (4) on spatially unbounded domains. We perform the same nondimensionalization as in

Refs. [29, 40] by introducing (i) $\rho_0$ a characteristic density of the problem, (ii) $\phi_0 \equiv (nM_{Pl}\Lambda^{n+4}/\beta\rho_0)^{1/(n+1)}$ the expectation value of the chameleon field in an ambient medium of density $\rho_0$ and (iii) $L_0$ a characteristic length scale of the system under study. Denoting the new dimensionless quantities with a tilde, trivial algebra leads to

$$ \alpha \tilde{\Delta} \tilde{\phi} = \tilde{\rho} - \tilde{\phi}^{-(n+1)}, \\ \text{with } \alpha \equiv \left( \frac{M_{Pl}\Lambda}{\beta L_0^2 \rho_0} \right) \left( \frac{n M_{Pl} \Lambda^3}{\beta \rho_0} \right)^{1/(n+1)}. \quad (7) $$

The mapping $(\beta, \Lambda) \mapsto \alpha$ for $n = 1$ is illustrated in Fig. 1. Note that Eq. (7) now only depends on two parameters, $\alpha$ and $n$, instead of the three initial ones, which allows for a more efficient numerical exploration of the chameleon parameter space<fnref n="3" />. The chameleonic acceleration (5) then scales as

$$ \mathbf{a}_{\phi} \propto \Lambda^{\frac{n+4}{n+1}} \beta^{\frac{n}{n+1}} \tilde{\nabla} \tilde{\phi}. \quad (8) $$

We denote $a_0$ the multiplicative constant appearing in front of the dimensionless gradient, which reads

$$ a_0 [\text{m/s}^2] = (\Lambda [\text{eV}] \times \varepsilon [\text{J/eV}])^{\frac{n+4}{n+1}} \frac{\beta^{\frac{n}{n+1}}}{M_{\text{Pl}} L_0} \left[ \frac{n M_{\text{Pl}}}{\rho_0 (\hbar c)^3} \right]^{\frac{1}{n+1}}. \quad (9) $$

In Eq. (9), physical quantities are expressed in SI units unless specified using square brackets and $\varepsilon \sim 1.6022 \times 10^{-19}$ J/eV is the conversion factor from electron-volts to joules. As a rule of thumb, the smaller $\alpha$, the more screened the setup. All physical results issued in this article are evaluated with $L_0 = R_\oplus = 6371$ km (the Earth radius) and $\rho_0 = 1$ kg/m<fnref n="3" />.

The Newtonian potential and the chameleon field do not have the same physical dimension. In order to be able to compare these two quantities, we define a new field

$$ \Psi = \frac{\beta}{M_{\text{Pl}}} \phi \quad (10) $$

which can be expressed in $m^2/s^2$. We refer to $\Psi$ as the *chameleon potential* since it plays the same role as $\Phi$. The total gravitational acceleration undergone by a point-like particle will simply be $-\nabla(\Phi + \Psi)$. Furthermore, the term ‘fifth-force’ will be used loosely throughout this article. Most occurrences of it should be taken as a synonym for ‘chameleon acceleration’, i.e. a quantity homogeneous to an acceleration and not a force per say. Finally, we will often refer to the ‘screened regime’ or to the ‘thin-shell of a body’ in this article. These notions

<footnote n="3">³ Naturally, the mapping $(\beta, \Lambda, n) \mapsto (\alpha, n)$ described above is not bijective.</footnote>
can be given precise meanings now that we have introduced the main notations. A macroscopic body is said to be *screened* when the chameleon field reaches the value that minimizes its effective potential $V_{\text{eff}}$ deep inside the body. In that case, the field remains essentially frozen in that body except in a (usually) thin surface layer, which is referred to as the thin-shell.

## B. Physical models

### 1. Mountains

At first order and seen from afar, planetary-mass objects have a rounded, ellipsoidal shape due to their self-gravity and rotation. It is only when we take a closer look at such bodies in the Solar System that smaller, more complex features become visible: mountains, ridges, craters, volcanoes, etc. This rich variety of topographies results in perturbations (with respect to the spherically symmetric case) in the gravitational field which, in the case of the Earth, can be measured by geodetic satellites. With a view to understand how 5th-forces affect Newtonian gravity in the vicinity of these topographical features, it is desirable to first work with a simple toy-model. We thus consider a spherical body together with a single, axisymmetric mountain on top of it as depicted in Fig. 2. It is mainly described by two dimensionless parameters:

- $h_m$, the height of the mountain divided by the radius $R_{\text{body}}$ of the spherical body (which is unitary on Fig. 2);

- $\theta_m$, the mountain's half-angle.

Note that these two parameters are deliberately exaggerated on Fig. 2 for better visualization, and are clearly not representative of any realistic mountain in the Solar system — see Table I. All numerical computations presented in this article were performed with $h_m = 10^{-2}$ and $\theta_m = 10^{-2}$ rad comparatively. The resulting setup is itself axisymmetric which means FEM computations can be performed in two dimensions rather than three, greatly reducing computational complexity.

For the model to be complete, we further need to specify the density function $\rho(\mathbf{x})$ inside and outside the body. For the sake of simplicity, we assign a constant density to the body $\rho_{\text{body}}$. The body may or may not be surrounded by an atmosphere. In either case, the density outside the body depends solely on the radial distance from the center $r$ and always goes down to a constant vacuum value $\rho_{\text{vac}}$. For all FEM computations, we set

$$ \tilde{\rho}_{\text{body}} = \frac{\rho_{\text{body}}}{\rho_0} = 10^3 \quad \text{and} \quad \tilde{\rho}_{\text{vac}} = \frac{\rho_{\text{vac}}}{\rho_0} = 10^{-15}. $$

Additionally, we will work most of the time with the dimensionless variable $\tilde{r} = r/L_0$, and set $L_0 = R_{\text{body}}$. The

<figure><img src="image_3.png" /><figcaption>FIG. 2. Mountain visualization and notations. The Cartesian frame $(O, x, y, z)$ is centered at the geometric center of the sphere devoid of mountain. The actual mountain profile used in numerical computations is drawn using B-splines in polar coordinates so as to form a smooth manifold.</figcaption></figure>

various fields involved in this study (chameleon potential, Newtonian potential, together with their gradient) will be probed at fixed discrete values of $\tilde{r}$ for the sake of consistency. We made the choice to show results for $\tilde{r} \in \{1.059, 1.111, 1.314, 4.645, 6.617\}$, which for the case of the Earth corresponds roughly to peculiar orbits: the International Space Station, MICROSCOPE, a Medium Earth orbit, Galileo and geostationary satellites, respectively.

### 2. Atmospheres

Some Solar system bodies are surrounded by an atmospheric layer — a gas envelope held in place by the gravity of the body. This slight over-density with respect to the case with no-atmosphere is expected to have an influence on the chameleon field profile and, therefore, on the 5th-force in space [4, 5, 37]. However, works that take account of the atmosphere often model it as an additional shell of matter with constant density satisfying $\rho_{\text{body}} > \rho_{\text{atm}} > \rho_{\text{vac}}$, or at best as a constant piecewise function [4, 34–37]. It is actually difficult to be more precise than this using analytical techniques only. Here, we go a step further by taking advantage of femtoscope to analyze the chameleon field profile in more realistic atmospheric setups. To avoid confusion, the requirement that the atmosphere must have a thin-shell stipulated in Ref. [4] only holds in the case of non-universal coupling, wherein unacceptably large violations of the weak equivalence principle would be observed in ground based experiments. Here, we work on the assumption of a universal coupling (characterized by a single dimensionless constant $\beta$) and so there is no particular reason for imposing this condition$^4$.

<footnote>$^4$ The fact remains that, even in the case of a universal coupling, deviations from the inverse square law can be suppressed by the</footnote>
<table><caption>TABLE I. List of some peculiar mountains in the Solar system<sup>a</sup>.</caption><thead><tr><th rowspan="2">Site</th><th rowspan="2">Body density [kg/m<sup>3</sup>]</th><th colspan="2">Atmosphere</th><th colspan="2">height (base to peak)</th><th rowspan="2">&theta;<sub>m</sub> [rad]</th></tr><tr><th>density [kg/m<sup>3</sup>]</th><th>thickness [km]</th><th>[km]</th><th>h<sub>m</sub></th></tr></thead><tbody><tr><td>Earth<br>Mount Everest</td><td>2.6 &times; 10<sup>3</sup><br>(Earth crust)</td><td>1.2 (sea level)</td><td>~ 100</td><td>4.6</td><td>7.2 &times; 10<sup>-4</sup></td><td>~ 10<sup>-3</sup></td></tr><tr><td>Earth<br>Mauna Kea</td><td>2.6 &times; 10<sup>3</sup><br>(Earth crust)</td><td>1.2 (sea level)</td><td>~ 100</td><td>10.2</td><td>1.6 &times; 10<sup>-3</sup></td><td>~ 10<sup>-2</sup></td></tr><tr><td>Mars<br>Mons Olympus</td><td>2582<br>(Mars crust)</td><td>2 &times; 10<sup>-2</sup> (max.)</td><td>~ 10</td><td>21.9</td><td>6.5 &times; 10<sup>-3</sup></td><td>~ 9 &times; 10<sup>-2</sup></td></tr><tr><td>Moon<br>Mons Huygens</td><td>2550<br>(Moon crust)</td><td>no atmosphere</td><td></td><td>5.5</td><td>3.2 &times; 10<sup>-3</sup></td><td>~ 6 &times; 10<sup>-2</sup></td></tr><tr><td>Io<br>Boösaule Montes</td><td>3500<br>(mean density)</td><td>&lt; 10<sup>-6</sup></td><td>—</td><td>18.2</td><td>10<sup>-2</sup></td><td>~ 1.5 &times; 10<sup>-2</sup></td></tr><tr><td>Vesta<br>Rheasilvia central peak</td><td>2800<br>(crust estimate)</td><td>no atmosphere</td><td></td><td>25</td><td>10<sup>-2</sup></td><td>~ 0.4</td></tr></tbody></table>

<sup>a</sup> Mainly based on https://en.wikipedia.org/wiki/List_of_tallest_mountains_in_the_Solar_System, last visited: August 22<sup>th</sup>, 2023

Three atmospheric density profiles are considered in this study: *Earth-like*, *Tenuous* and *Dense*. The Earth-like model is built from the 1976 version of the U.S. Standard Atmosphere model [41], commonly known as the US76 model<sup>5</sup>. It provides an estimate of the Earth atmospheric density $\rho_{US}$ as a continuous function of the altitude, up to $R_{atm} \sim 36 \times 10^3$ km. Because we want the minimum dimensionless density in the numerical domain to be exactly $\tilde{\rho}_{vac} = 10^{-15}$, we apply the following transformation on the original data:

$$ \log \tilde{\rho}_{\text{Earth-like}} = \log \tilde{\rho}_{\text{US}} + k \left[ \log \tilde{\rho}_{\text{US}} - \log(\min \tilde{\rho}_{\text{US}}) \right] \\ \text{with } k = \frac{\log(\tilde{\rho}_{\text{US}}/\tilde{\rho}_{\text{vac}})}{\log(\max \tilde{\rho}_{\text{US}}/\min \tilde{\rho}_{\text{US}})} $$

for $r < R_{atm}$, which is nothing but an affine transformation on the logarithmic densities. Beyond $R_{atm}$, we set $\tilde{\rho}_{Earth-like} = \tilde{\rho}_{vac}$. The other two models — Tenuous and Dense — are purely empirical in the sense that they are not based on actual atmospheric data. Both are constructed via the expression

$$ \log \tilde{\rho}(r) = \begin{cases} A \exp \left[ \frac{(r - R_{\text{atm}})^2}{\sigma^2} \right] + B & \text{if } r < R_{\text{atm}}, \\ \log \tilde{\rho}_{\text{vac}} & \text{otherwise} \end{cases} $$

where the parameters $(A, B, \sigma)$ are adjusted by hand to obtain either a very tenuous, thin atmosphere or a very dense, thick one instead. The resulting density profiles are depicted in Fig. 3.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Atmospheric profiles investigated in this study.</figcaption></figure>

### C. Decomposition of scalar fields into spherical harmonics

In geophysics and physical geodesy, the Earth gravitational potential is conveniently modeled as a spherical harmonics expansion [42]. Any well-behaved function $f: \mathbb{R}^3 \to \mathbb{R}$ may be decomposed as

$$ f(r, \mathbf{n}) = \sum_{l=0}^{+\infty} \sum_{m=-l}^{+l} f_{lm}(r) Y_{lm}(\mathbf{n}), \quad (11) $$

where $r, \mathbf{n} = (\theta, \varphi)$ refer to spherical coordinates, $Y_{lm}$ is the real spherical harmonic function of degree $l$ and order $m$ (see Ref. [43] for its definition), and $f_{lm}$ are the spherical harmonic coefficients that only depend on the radial coordinate — they are referred to as the bare coefficients in this article. There are several normalization conventions for an unequivocal definition of spherical harmonic functions. In this study, we stick to the orthonormalized

<footnote>atmosphere.</footnote>

<footnote n="5">5 Data downloaded from http://www.braeunig.us/space/atmos.htm, (especially for the density between 1000 km - 36000 km altitude). Last visited: June 1<sup>st</sup>, 2022.</footnote>
<table><caption>TABLE II. Notations for the spherical harmonic coefficients.</caption><thead><tr><th></th><th>Bare coefficients</th><th>Rescaled coefficients</th></tr></thead><tbody><tr><td>Newtonian potential</td><td>Φ<sub>lm</sub>(r)</td><td>y<sup>N</sup><sub>lm</sub></td></tr><tr><td>Chameleon potential</td><td>Ψ<sub>lm</sub>(r)</td><td>y<sup>C</sup><sub>lm</sub>(r)</td></tr></tbody></table>

convention for which

$$ \int_{S^2} Y_{lm}(\mathbf{n})Y_{l'm'}(\mathbf{n}) d^2\Omega = \delta_{ll'} \delta_{mm'}, \quad (12) $$

where $S^2$ is the unit 2-sphere, $d\Omega$ is the differential surface $\sin(\theta)d\theta d\phi$ and $\delta_{ij}$ is the Kronecker delta function. The notations used to refer to the spherical harmonic coefficients of the Newtonian potential $\Phi$ and the chameleon potential $\Psi$ are gathered in Table II.

### 1. Rescaled coefficients

The bare spherical harmonic coefficients of the Newtonian potential $\Phi_{lm}$ further exhibit a scaling property. Let us denote by $\mu_{body} \equiv GM_{body}$ the standard gravitational parameter of the central body of mean radius $R_{body}$ and mass $M_{body}$. Then, the rescaled coefficients

$$ y_{lm}^{N} = \frac{r}{\mu_{body}} \left( \frac{r}{R_{body}} \right)^{l} \Phi_{lm}(r) \quad (13) $$

can be shown to be independent of $r$ [44],<fnref n="6" /> owing to the specific form of the Newtonian potential (1). Such rescaled coefficients are thus *universal* to the body under consideration. Similarly to Eq. (13), we denote by $y_{lm}^{C}(r)$ the rescaled coefficients of the chameleon potential which, for their part, have no particular reason to be independent of the radial distance. In that sense, Ref. [44] shows the explicit dependence of such coefficients with respect to $r$ in the case of a Yukawa interaction.

This relation can also serve as a means of checking the numerical results obtained for the Newtonian potential. This test is performed in Appendix B.

### 2. Recovery of the coefficients

We use the software SHTools [43] to compute the spherical harmonic coefficients of the scalar fields of interest. The PYTHON package pyshtools comes with the routine SHGrid.expand which calculates the coefficients

by means of some numerical quadrature<fnref n="7" />. The only detail worth mentioning is the fact that this routine outputs separate variables for the cosine $C_{lm}$ and sine $S_{lm}$ coefficients (sometimes referred to as the Stokes coefficients). The conversion from $(C_{lm}, S_{lm})$ to bare coefficients is outlined in Appendix A — Eq. (A5).

## D. Numerical techniques

### 1. Using femtoscope to solve linear and nonlinear PDEs with asymptotic boundary conditions

As mentioned earlier, femtoscope is a ready-to-use PYTHON program which plays a central role in this study as it enables us to compute both the Newtonian potential and the chameleon field by solving Eqs. (2) and (4) respectively. It is based on the finite element method — building on top of the open-source package Sfepy [45] — and further implements techniques to deal with nonlinearities and asymptotic boundary conditions (3, 6).

The proper treatment of these asymptotic boundary conditions is of noticeable importance in this study. Indeed, it is tempting to simply truncate the numerical domain at a fixed radius and apply a homogeneous Dirichlet boundary condition on the artificial border resulting from that process. This procedure has several flaws:

1. For the error that arise therefrom to be small, the domain must be sufficiently large, which translates to higher computational cost.

2. Selecting the size of that domain is a blind experiment in the sense that the dependence of the error on the truncation radius is not easily accessible without additional tricks.

3. It wantonly imposes spherical symmetry on the solution as we approach the artificial boundary. This is particularly undesirable in this study where we are interested in the small deviations from spherical symmetry introduced by the presence of the mountain.

This latter point is illustrated on Fig. 4 where it can be seen that, as we approach the artificial boundary, the truncation method (labeled ‘FEM bounded’, dash-dotted pink line) exhibits a poor approximation.

Instead, we employ a technique based on the splitting of the numerical domain $\Omega$ into two subdomains $\Omega_{int}$ and $\Omega_{ext}$ such that $\bar{\Omega} = \bar{\Omega}_{int} \cup \bar{\Omega}_{ext}$. $\Omega_{int}$ is the bounded, interior domain, while $\Omega_{ext}$ is the unbounded, exterior domain. An inversion transform is then applied to $\Omega_{ext}$, resulting in a bounded domain $\bar{\Omega}_{ext}$ (called the *inversed exterior domain*) which can be meshed on a

<footnote n="6">⁶ The numerical values of $\mu_{body}$ and $R_{body}$ could in theory be chosen arbitrarily. However the numerical values of the rescaled coefficients are tied to this choice.</footnote>

<footnote n="7">⁷ In this study, we use a $N \times 2N$ Driscoll and Healy sampled grid.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Orthoradial profiles of the dimensionless Newtonian potential $\delta\tilde{\Phi}$ sourced by the mountain at three different altitudes, corresponding to $\tilde{r}_1 = 1.059$, $\tilde{r}_2 = 4.645$, and $\tilde{r}_3 = 6.617$ (top, middle and bottom panels respectively). The black dots together with their error bar represent the benchmark solution, obtained through the computation of the integral Eq. (1) with scipy's pltquad routine. The pink dash-dotted line is obtained by solving Poisson's equation (2) with an homogeneous Dirichlet boundary condition applied at $\tilde{r} = \tilde{R}_c = 7$ while the green solid line is the solution provided by femtoscope with asymptotic boundary condition. Finally, the gray dash-dotted line is an analytical approximation where the mountain is replaced by a point-mass, whose location and mass were fitted to provide a good match with respect to the benchmark: $m_{\text{mountain}}/M_{\text{body}} = 2.23 \times 10^{-7}$ and $z/R_{\text{body}} = 2.22 \times 10^{-3}$.</figcaption></figure>

computer. There are many possible numerical implementations based on this method, see e.g. Refs. [46–48]. In this study, we make use of the so-called *virtual connection of DOFs* described in our previous work [29].

## 2. Numerical challenges and verification

There are several inconspicuous challenges associated with the numerical computation of the field profiles in the setup described in Sec. II B. To start with, let us stress the fact that we are looking for small deviations from spherical symmetry, owing to the presence of a very localized over-density at the pole that we here call a mountain. Quantitatively speaking, a back-of-the-envelope calculation shows that — at a fixed altitude $h$ — the relative variation of the Newtonian potential $\Phi(R_{\text{body}} + h, \theta)$ along the latitudes with respect to its mean value at this altitude is no larger than a few $10^{-6}$. The higher we go, the smaller this ratio, which means our numerical approximations have to be correct up to at least seven significant digits to be deemed good. This mere order-of-magnitude calculation raises an additional concern: how do we actually check that the numerical approximations we obtain are compliant with the required levels of precision?

The Poisson's equation (2) governing the Newtonian potential being linear, it is possible to apply the superposition principle, where the total field is simply the mountain's contribution on top of a spherically symmetric background: $\Phi_{\text{tot}}(r, \theta) = \delta\Phi(r, \theta) + \Phi_0(r)$. Turning to the chameleon field, the nonlinearity in the r.h.s. of the Klein-Gordon equation (4) prevents us from following the same path. Even if one were to decompose the chameleon field as $\phi_{\text{tot}}(r, \theta) = \delta\phi(r, \theta) + \phi_0(r)$, the term $(\phi_0 + \delta\phi)^{-(n+1)}$ becomes linearizable only under the assumption that $\delta\phi \ll \phi_0$ everywhere. Unfortunately, this assumption has no reason to hold in all scenarios, owing to the very nature of the screening mechanism. Indeed, it is far from being valid in the case where the mountain itself becomes screened, which turns out to be the most interesting case given the current constraints on the chameleon field [30]. For lack of a better workaround, we abandoned perturbation-based techniques and put our efforts into solving for the full field. It is therefore necessary to compare the FEM approximation obtained with femtoscope against some benchmark. Failing to have an analytical solution for the Newtonian potential of a mountain, we can still resort to the numerical integration of Eq. (1). In this respect, we use scipy's pltquad routine [49] to evaluate the integral with an estimated relative error of a few $10^{-9}$. This semi-analytical approach constitutes our benchmark and is depicted by the black dots together with their error bar in Fig. 4. Note that while it takes only a few seconds to evaluate the potential at a single point with this method, it is not conceivable to construct a full map of the field in this way. Rather, this semi-analytical computation should be employed sparingly to assess the error of the FEM computations.

In contrast, the chameleon field does not enjoy a similar integral representation which in turns means that we cannot easily define a benchmark profile. Nonetheless, we came up with the following strategies:
<figure><img src="image_6.png" /><figcaption>FIG. 5. Evolution of the dimensionless chameleon profile $\tilde{\phi}$ (top) and the associated residual (bottom) after different numbers of iterations of Newton's method (2, 30 and 120). These two quantities are displayed as a function of $\theta$ at fixed altitude $\tilde{r} = 1.059$. The residual becomes stationary and thus no longer decreases after a sufficient number of iterations has been reached.</figcaption></figure>

- Select the set of FEM-related parameters (number and distribution of DOFs, order of the base functions, etc.) so that the FEM-approximation of the Newtonian potential matches the benchmark and use those parameters for the FEM computation of the chameleon field. The light green curves on Fig. 4 correspond to such FEM-approximations (using the aforementioned ‘virtual connection of DOFs’ method) and show that it is indeed possible to reach a high level of accuracy as they stay within the error bars of the benchmarks.

- It is also good practice to refer to pre-established *FEM convergence curves*, which are simple charts relating the error to the number of DOFs — see e.g. Fig. 1 of Ref. [50]. We can then construct our meshes in an enlightened way, ensuring they are fine enough to meet the stated accuracy.

- Evaluate the *strong residual*, which can be done by inputting the FEM approximation obtained for the

chameleon field into its equation of motion (7) — schematically: Residual = $\alpha\tilde{\Delta}\tilde{\phi} - \tilde{\rho} + \tilde{\phi}^{-(n+1)}$. The closer the quantity is to zero, the better the numerical approximation. In order to make this criterion more quantitative, we can monitor (i) the strong residual's decrease throughout the Newton's iterations (see Fig. 5) and especially how small the final residual is compared to the initial one, and (ii) its size relative to the size of each term in it: the final residual should be at least a few orders of magnitude smaller than the dominant terms. This criterion is assessed on all 2D numerical computations of the chameleon field discussed in this article. As an example, Fig. 17 in Appendix C demonstrates that this criterion is indeed met on three distinct numerical solutions, at three altitudes.

Yet, formulating criteria based on the strong residual alone is not entirely satisfactory as it is an absolute quantity. Consequently, there is a priori no simple connection between the relative error committed on the approximation and the strong residual, since the latter quantity is dependent on the PDE's parameters (value of $\alpha^8$, density model, etc.). Computing a reduction factor, that is by how much the strong residual has decreased over the Newton's iterations, is not sufficient either as it depends on how well the initial guess has been chosen (see discussion in the next paragraph). One idea to break this deadlock is to compare our numerical approximations with the chameleon radial profile around a ball. Indeed, the spherically symmetric case is much more under control as we have analytical approximations at our disposal (see e.g. Refs. [4, 51]) and the Klein-Gordon equation boils down to a one-dimensional ordinary differential equation (ODE) which can be solved numerically with a much higher density of DOFs and higher-order finite elements. In terms of residual, the numerical solutions are actually better than their analytical counterparts (see e.g. Tab. II from Ref. [29]), which is why we propose to use 1D numerical solutions as a benchmark for the spherically symmetric case. Because the addition of a mountain on top of the spherical planet is not expected to have a huge impact on the field's strength outside it, we can check that the evolution of the field along the outgoing radial direction follows that of the benchmark. We provide a quantitative way of assessing that statement in Appendix C, which is applied for all the numerical solutions discussed in this article. Finally, the orthoradial variations of the field at fixed altitudes seems more difficult to verify. As a rough check, we can set $h_m = 0$ and verify that this leads to $\partial_{\theta}\phi = 0$. In practice, we do not expect this equality to hold exactly so we rather make sure that the amplitude $\max_{\theta}\phi(r, \theta) - \min_{\theta}\phi(r, \theta)$

<footnote n="8">⁸ In particular, we observed in Ref. [29] that the 2-norm of the strong residual was increasing with $\alpha$, all other things being equal.</footnote>
is much smaller in the case $h_m = 0$ compared to the case $h_m = 0.01$. Doing this sanity check on a handful of cases (doing it on all cases would have been too costly) consistently shows that the two quantities differ by at least two orders of magnitude, so that we can confidently state that the orthoradial profiles showed later on do not originate from numerical noise.

Ultimately, the most critical point in this FEM computation is the convergence of the Newton's iterations. Whether or not the method converges depends on a lot of factors. Unfortunately, there are no miracle techniques to address convergence issues but rather *recipes* and good practices which we concisely report here. Perhaps the most important one is to start from a *good* initial guess, i.e. an initial approximation that is as close as possible to the true solution. In most cases, we use a pre-computed 1D radial profile of the field to this end. Other common practices are to refine the meshes where the field is expected to vary quickly (large gradient) — that is near the transition between the inside and the outside of the body, and near the area representing spatial infinity in the inverse exterior domain $\hat{\Omega}_{\text{ext}}$ — or to tweak the relaxation parameter [52]. Additionally in the particular case of the chameleon field entering the so-called screened regime, we can get rid of the region of the mesh $r < R_{\text{screened}}$ where the field is screened (i.e. constant) and apply a Dirichlet boundary condition at $r = R_{\text{screened}}$. When all the above failed, we resorted to so-called *ramping* [53, 54] or *numerical continuation* methods [52, 55]. For example in some cases, we would gradually vary the $\alpha$ parameter from Eq. (7) from a value where the solution is known to the desired value which is problematic convergence-wise, using the solution at each intermediate step as an initial guess of the next one. In spite of all these additional tricks, some combinations of $\{\alpha, \text{atmosphere model}\}$ resisted all our attempts and were thus discarded from this study. As a closing remark, let us emphasize the fact that we made use of many widely spread techniques in the literature for nonlinear FEM problems (see e.g. Ref. [52] chapter 4), both for implementation and verification purposes. While we are unable to quantify the relative error made on each solution obtained in this study, we grant them a sufficiently high level of confidence that the orders of magnitude discussed hereafter are correct, leaving the physical conclusions unchanged.

In total, we ran FEM computations for four different density profiles outside the main body — the constant vacuum value $\tilde{\rho}_{\text{vac}} = 10^{-15}$ as well as the three atmospheric models depicted in Fig. 3 — and for $\log_{10} \alpha \in \{-5, \dots, -28\}$. This amounts to nearly a hundred problems to solve on meshes with roughly $10^6 \mathbb{P}_2$-triangles. The computations were performed on an ONERA's computing platform equipped with Broadwell and Cascade Lake nodes.

### III. Modified gravity around and above a mountain

In this section, we present and analyze simulations results. We start off with the atmosphere-free case before discussing the influence of each atmospheric models. Due to the parameters' degeneracy mentioned in Sec. II A 2, we decided to fix $\Lambda = \Lambda_{\text{DE}}$ for all numerical results and figures presented in the following. One can refer to Fig. 1 to get a better grasp of the $(\beta, \Lambda, n) \mapsto (\alpha, n)$ mapping.

#### A. Simulation of an atmosphere-free planet

##### 1. Gravitational potential profiles and spherical harmonics decomposition

The total gravitational potential is the sum of the Newtonian potential $\Phi$ and the chameleon potential $\Psi$ as defined in Sec. II A. These are the direct results of FEM computations, i.e. *femtoscope's* outputs. As this raw data can sometimes be noisy, we had recourse to smoothing splines notably for post-process operations requiring the evaluation of the fields outside mesh data points like the computation of spherical harmonic coefficients using SHTools [43]. Note that the azimuthal symmetry of our setup imposes that the only nonzero coefficients are the ones for which the order $m$ is equal to zero.

In Fig. 6, we represent the potential profiles as a function of the colatitude $\theta$ — the radial coordinate being fixed at $\tilde{r} = 1.059$ — (left column) and their associated spherical harmonic coefficients for degrees $l \in \{1, \dots, 100\}$ (right column). The top row corresponds to the specific case of the Newtonian potential while the following four rows correspond to chameleon potentials for $\log_{10} \alpha \in \{-4, -6, -15, -25\}$ respectively. The Newtonian potential $\Phi$ is by far the largest contribution to the total potential: roughly $-10^7 \text{ m}^2/\text{s}^2$ compared to $0.2 \text{ m}^2/\text{s}^2$ for the chameleon potential in the $\alpha = 10^{-25}$ case (last row of the same figure). The variation of the potential with respect to $\theta$ around this mean value has the same kind of shape where both ends of the curves have a slope that goes down to zero for symmetry reasons. Note that the potential is always smaller at $\theta = 0$ than at $\theta = \pi$. This is because the mass excess that the mountain represents is located at $\theta = 0$, forming a deeper potential well.

While all four potential profiles share this apparently common trend, the spherical harmonic coefficients displayed on the right column reveal important differences and two types of spectrum emerge. On the one hand, the Newtonian potential and the chameleon potential for $\alpha = 10^{-4}$ have a similar, monotonically decreasing spectrum. This is due to the fact that here, the chameleon field is unscreened which means that all the mass of the main body contributes to the field just like in the Newtonian case. On the other hand, as soon as $\alpha < 10^{-5}$, the chameleon field enters the screened regime, changing the shape of the spectra. We recall that the smaller $\alpha$,
the more screened the setup. These spectra all have a maximum for $l > 1$. This distinctive feature of screened chameleon potentials, which could not be seen by eye on the left-hand-side curves, is nevertheless small in front of the Newtonian potential's coefficients $\Phi_{l0}$.

As $\alpha$ decreases, the chameleon potential mean value increases. Indeed, we have $\Psi = K\tilde{\phi}$ where $K \propto \alpha^{-(n+1)/(n+2)}$ at fixed $\Lambda$. As a result, the spherical harmonic coefficients also get amplified as $\alpha$ decreases, leading to a more disturbed gravitational potential.

Fig. 18 in Appendix D further shows how the spectra evolve as the altitude is increased for the Newtonian potential and the chameleon potential ($\alpha = 10^{-25}$).

## 2. Newtonian gravity and fifth-forces

Once the gravitational potential is known, the actual gravitational acceleration is easily derived by computing its gradient. It is convenient to decompose the acceleration vector $\mathbf{a}$ onto the unit vectors ($\mathbf{e}_r, \mathbf{e}_\theta$) (there is no component of the acceleration on $\mathbf{e}_\varphi$ due to rotational invariance) such that

$$ \mathbf{a} = a_r \mathbf{e}_r + a_\theta \mathbf{e}_\theta. $$

In practice, the dimensionless gradient is computed numerically and then multiplied by the relevant coefficient $a_0$ with units m/s² — see Eq. (9). Fig. 7 gives an overview of both Newtonian acceleration (top panel) and fifth-forces for $\log_{10} \alpha \in \{-15, -27, -28\}$. Specifically, we represent the component $a_r$ (purple curve) and $a_\theta$ (crimson curve) as a function of $\theta$ while the altitude is held fixed at $\tilde{r} = 1.059$.

An important point to discuss here is the limit $\alpha \to 0$. On the one hand, we have seen that for chameleon gravity, $a_0$ is proportional to $\alpha^{-(n+1)/(n+2)}$ at fixed $\Lambda$, and consequently

$$ a_0 \xrightarrow{\alpha \to 0^+} +\infty \text{ with constraint } \Lambda = \Lambda_{DE}. $$

This gives the impression that one can make the fifth-force as large as desired simply by taking an ever-decreasing value of $\alpha$. On the other hand, we know that in the limit $\alpha = 0$, the chameleon field profile is trivially given by $\tilde{\phi} = \tilde{\rho}^{-1/(n+1)}$ (take $\alpha = 0$ in Eq. 7). Yet for altitudes higher than the mountain's height, our models are such that $\partial_\theta \tilde{\rho} \equiv 0$ so that $\tilde{a}_\theta$ is expected to vanish for sufficiently small values of $\alpha$. In front of this apparent paradox, we raise two points:

1. Taking the limit $\alpha \to 0$ at fixed $\Lambda$ coerces $\beta \to +\infty$. Yet, a glimpse at the chameleon constraints plot (see e.g. Fig 3 from Ref. [23]) reveals that chameleon models with $\beta > 10^{14}$ are ruled-out by precision atomic tests. In our case, this corresponds to $\alpha < 10^{-37}$, which is out of the range of $\alpha$ values covered in this study.

<table><caption>TABLE III. Assessment of the maximum fifth-force at $\tilde{r} = 1.059$</caption><thead><tr><th></th><th>arg max<sub>α</sub></th><th>value</th></tr></thead><tbody><tr><td>$\tilde{a}_r$</td><td>$10^{-24}$</td><td>$1.47 \times 10^8$</td></tr><tr><td>$a_r$</td><td>$10^{-25}$</td><td>$1.46 \times 10^{-7}$ [m/s<sup>2</sup>]</td></tr><tr><td>$\tilde{a}_\theta$</td><td>$10^{-25}$</td><td>$1.09 \times 10^6$</td></tr><tr><td>$a_\theta$</td><td>$10^{-25}$</td><td>$1.40 \times 10^{-9}$ [m/s<sup>2</sup>]</td></tr><tr><td>$||\tilde{\mathbf{a}}||$</td><td>$10^{-24}$</td><td>$4.58 \times 10^9$</td></tr><tr><td>$||\mathbf{a}||$</td><td>$10^{-25}$</td><td>$4.50 \times 10^{-6}$ [m/s<sup>2</sup>]</td></tr></tbody></table>

2. Another argument that does not involve referring to current chameleon constraints can be made on the basis of Figs. 7 and 8. On Fig. 8, we decompose $a_\theta$ into the product $a_0(\alpha)$ times $\tilde{a}_\theta = \partial_\theta \tilde{\phi}/\tilde{r}$. The two terms of this product both depend on $\alpha$: while $a_0$ is simply a power law of $\alpha$, $\tilde{a}_\theta$ clearly exhibits the phenomenon aforementioned, namely that the dimensionless gradient — after reaching a peak for $\alpha = 10^{-25}$ — vanishes for $\alpha < 10^{-28}$. When multiplied together, these two quantities result in $a_\theta$ which is scattered in log-scale on the bottom panel of this figure. We recognize the power law behavior $a_\theta \propto \alpha^{-n/(n+2)}$ in the range $[10^{-10}, 10^{-21}]$ where $\tilde{a}_\theta$ is roughly constant, followed by a sharp decline due to the vanishing of $\tilde{a}_\theta$. This explains why on Fig. 7, the transition from $\alpha = 10^{-27}$ to $\alpha = 10^{-28}$ completely destroys $\nabla\Psi$. There only remains numerical noise, whose amplitude has no genuine physical meaning.

We also conducted the same analysis as performed in Fig. 8 for the radial component of the acceleration vector as well as for its norm (still for $\tilde{r} = 1.059$). The results are reported in Table III. In terms of dimensionless quantities, the orthoradial component $\tilde{a}_\theta$ is maximum for $\alpha = 10^{-25}$ while the radial component and the norm of the gradient are both maximized for $\alpha = 10^{-25}$. When these quantities are expressed in units homogeneous to an acceleration (m/s²), $\alpha = 10^{-25}$ is again the argument that maximizes them all. This corresponds to the traditional parameters $(\beta, \Lambda, n) \sim (10^6, \Lambda_{DE}, 1)$.

Besides, we can discuss the direction of the acceleration vector by computing the ratio $a_\theta/a_r$ for both Newtonian and chameleon gravity. We conclude that the Newtonian part of the total acceleration vector is very radial, with $\max_\theta (a_\theta/a_r) \le 3 \times 10^{-5}$, whereas the chameleon acceleration has a more significant orthoradial component since $\max_\theta (a_\theta/a_r) \le 10^{-2}$ for $\log_{10} \alpha = -25$. The physical interpretation for this discrepancy is that in the screened regime, the chameleon acceleration is sourced only by a thin outer layer of the planet which is commonly referred to as the *thin-shell* [4].

Finally, one may be surprised by the fact that maximum fifth forces are obtained for values of the parameters $(\beta, \Lambda, n)$ which belong to the thin-shell regime; as
<figure><img src="image_7.png" /><figcaption>FIG. 6. Newtonian and chameleon potential profiles (left column) together with their spherical harmonic coefficients up to degree 100 (right column) computed at $\tilde{r} = 1.059$. The top row corresponds to the Newtonian case while the four remaining rows are associated with chameleon potentials with $\log_{10} \alpha \in \{-4, -6, -15, -25\}$ from top to bottom. The monopole ($\Phi_{00}$, $\Psi_{00}$) is deliberately excluded from the bar graphs because it is only dependent on the field's mean value. All quantities are expressed in $m^2/s^2$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7. Gravitational field $\mathbf{a} = a_r\mathbf{e}_r + a_\theta\mathbf{e}_\theta$ in Newtonian gravity (top) and in the chameleon model for the set of parameters $\{\alpha \in \{10^{-15}, 10^{-27}, 10^{-28}\}, n = 1, \Lambda = \Lambda_{\text{DE}}\}$ at $\tilde{r} = 1.059$. The orthoradial acceleration $a_\theta$ is depicted by the red curve (left axis) while the radial acceleration $a_r$ is depicted by the purple curve (right axis).</figcaption></figure>

this appears to contradict the usual rule of thumb that fifth forces should be suppressed in this regime. There is actually no contradiction, provided we clearly define the context. Indeed, the total fifth force acting on a given macroscopic body can be computed via the integration of the gradient of the field on its whole volume — as done later in Sec. IV A 2 for instance. It is true that, if

<figure><img src="image_9.png" /><figcaption>FIG. 8. Study of the chameleon orthoradial acceleration $a_\theta = a_0 \times \tilde{a}_\theta$ with respect to $\alpha$ at fixed $\Lambda$ and fixed altitude $\tilde{r} = 1.059$. The top panel features each term separately, $\tilde{a}_\theta$ in magenta dots (dimensionless) and $a_0$ as the blue curve (in m/s²). The bottom panel is simply the product of these two terms $a_\theta$ (in m/s²). Finally, the red vertical dotted lines correspond to values of $\alpha$ for which the FEM computation was deemed unsatisfactory (failure of Newton's method to converge or unacceptably large residuals).</figcaption></figure>

the macroscopic object at stake has a thin-shell, the integral of the gradient of the field vanishes everywhere but in that thin-shell, greatly reducing the overall fifth force experienced by that body. Here however, the situation is radically different: we are interested in the fifth force experienced by a point-mass (which by essence, cannot possess a thin-shell) sourced by a mountainous planet. In this framework, we merely observe that increasing the value of the coupling constant $\beta$ in a given range, while keeping $n$ and $\Lambda$ fixed, results in greater fifth forces. Incidentally, increasing $\beta$ while keeping $n$ and $\Lambda$ fixed means decreasing $\alpha$ (see Fig. 1) and results in a more screened body. This phenomenon was already observed in Figs. 14 and 15 of our previous work [29] and in Fig. 7 of Ref. [56] for instance. This can also be understood in the framework of the analytical approximation of the chameleon fifth force for spherical objects. Taking Eq. (2.64) from Ref. [51] reads

$$a_{\phi} = 3\beta^2 \frac{\Delta R}{R} \frac{GM_{\text{ball}}}{r^2} (1 + m_{\phi}) e^{-m_{\phi}(r-R_{\text{ball}})}.$$
In this expression, $\Delta R/R \propto \beta^{-1}$ and $m_{\phi} \propto \beta^{\frac{n+2}{2(n+1)}}$ so that, at fixed $r > R_{\text{ball}}$, the function: $\beta \mapsto a_{\phi}$ is increasing on the interval $[0, \beta^*[\text{and decreasing on } ]\beta^*, +\infty[, \text{ for a certain parameter } \beta^* > 0$. We also see that $a_{\phi} \to 0$ when $\beta \to \infty$, because of the exponential term. This is exactly the phenomenology that we observe on numerical simulations: beyond a certain value of $\beta \sim 10^8$ (i.e. below a certain value of $\alpha$, which is around $10^{-28}$), the fifth force vanishes — see Figs. 7 and 8.

## B. Adding an atmosphere

Here we study the influence of adding an atmosphere to the density model. To the best of our knowledge, only a handful of studies deal with the influence of the atmosphere [4, 5, 34–36]. In this section, we address simple questions: how is the fifth-force mitigated by the presence of an atmosphere? Can the mountain still be somehow seen in the field profile? How does all of this depend on the atmospheric model?

Part of the answer can be unveiled by first studying a simplified version of the setup. More precisely, we got rid of any orthoradial dependence in the density distribution function — which amounts to taking the mountain out of our model to end up with a purely radial setup. This simplification allows us to perform computationally inexpensive 1D FEM simulations with femtoscope and still get valuable insight into how the chameleon field behaves in the presence of an atmosphere. We ran computations for all atmospheric models outlined in Sec. II B 2 and for all values of $\alpha \in \{10^{-5}, \dots, 10^{-29}\}$. Part of this simulation campaign has been compiled into Fig. 9, where the vast range of $\alpha$ values explored has been boiled down to only four distinct values for the sake of clarity and conciseness. On each sub-panel, the grey dash-dotted line is associated with the fully screened profile obtained in the limit $\alpha = 0$ which is given by $\tilde{\phi}(\tilde{r}) = \tilde{\rho}(\tilde{r})^{-1/(n+1)}$. Contrary to the previous atmosphere-free case, where the radial density would have been a mere Heaviside step function, the atmospheric density function smoothly interpolates between $\tilde{\rho}_{\text{atm}}(\tilde{r} = 1)$ to $10^{-15}$ such that the asymptotic profile’s gradient $\partial_r \tilde{\phi}(\alpha \to 0)$ is not identically zero.

It is only when we put these 1D chameleon profiles into perspective with the full 2D simulation’s results that a clear understanding of the influence of the atmosphere emerges. Starting from $\alpha = 10^{-5}$ and gradually decreasing the value of this parameter, we witness the succession of several regimes:

1. For the larger values of $\alpha$, the planet is not fully screened, i.e. there is still a thin-shell. This can be seen on the first column of Fig. 9 ($\alpha = 10^{-7}$) where the chameleon kicks in (i.e. departs from limit profile $\tilde{\phi}(\alpha \to 0)$) before $\tilde{r} = 1$ (which corresponds to the transition between the planet and the atmosphere). This regime is particularly visible

on sub-panels a) to d). The impact of the atmosphere on the fifth-force at higher altitudes is then minor — see Tab. IV thereafter where we compare the amplitude of the fifth-force with and without atmosphere at $\tilde{r} \in \{1.059, 1.314\}$.

2. At some point when decreasing $\alpha$, the lowest part of the atmosphere becomes screened itself. This is especially visible on sub-panels e) to h). We provide a zoomed-in view of this very region in order to be able to compare the fraction of the atmosphere that is screened against the relative size of the mountain $\tilde{h}_m = 0.01$. As soon as the screened area overflows the mountain, i.e. everything below $\tilde{r} = 1.01$ is screened, the imprint of the mountain of the chameleon field is definitely lost at higher altitudes. In other words, the orthoradial acceleration vanishes, giving way to numerical noise. This is why some entries of Table IV are set to N/A. When it comes to radial component of the fifth-force, it is hardly modified compared to the scenario without atmosphere.

3. For even smaller values of $\alpha$, the screening eventually reaches the probed region at high altitude. This is particularly clear in sub-panels i) to l), where the chameleon field profile is getting closer to the limit profile (grey dash-dotted line). Here, the orthoradial acceleration remains drowned in the numerical noise while the radial acceleration is fully dictated by the density profile. This is why in some cases, $a_r$ can even become larger with an atmosphere than without (see entries of Table IV greater than unity).

Once we know that these three regimes exist regardless of the specific form of the atmospheric profile (as long as density decreases with altitude), we can start to be more quantitative by

- specifying where the transition between each regime occurs for the atmospheric models at stake;

- computing the attenuation factor on the fifth-force for the different density models.

These quantitative results are reported in Table IV where each entry is a pair $(\mu_r, \mu_\theta)$ defined as

$$ \mu_r = \max_{\theta} a_r^{\text{with-atm}}(\tilde{r}, \theta) / \max_{\theta} a_r^{\text{no-atm}}(\tilde{r}, \theta) \\ \mu_\theta = \max_{\theta} a_\theta^{\text{with-atm}}(\tilde{r}, \theta) / \max_{\theta} a_\theta^{\text{no-atm}}(\tilde{r}, \theta) \quad (14) $$

at a specific radial coordinate $\tilde{r}$. We refer to these coefficients as the attenuation factors, which are of course dependent on the atmospheric model as well the altitude at which they are computed.

The take home message from this study of atmospheric models is that the presence of an atmosphere, as tenuous as it may be, prevents access to the biggest fifth-force attainable without atmosphere. Indeed, we saw earlier
<figure><img src="image_10.png" /><figcaption>FIG. 9. Radial profiles of the chameleon field for $log_{10} \alpha \in \{-7, -10, -18, -23\}$ (columns) and for all three atmospheric models defined in Sec. II B 2, namely *Tenuous*, *Earth-like* and *Dense* (rows). On each sub-panel, the grey dash-dotted line corresponds to *fully screened* profile, that is obtained in the limit $\alpha = 0$ and is given by $\phi(\alpha \to 0) = \bar{\rho}^{-1/(n+1)}$. The radial chameleon profile is depicted by the bi-color solid line, where the transition from the darker color to the lighter one occurs at the chosen interface ($\tilde{r} = 7$) between the interior domain and the kelvin-inversed exterior domain (see Ref. [29] for more details).</figcaption></figure>

<table><caption>TABLE IV. Attenuation coefficients of the radial and orthoradial component of the chameleon acceleration with an atmosphere compared to the atmosphere-free case. The first number of each pair corresponds to the radial part and is computed as $\max_{\theta} a_r^{\text{with-atm}}(\tilde{r}, \theta) / \max_{\theta} a_r^{\text{no-atm}}(\tilde{r}, \theta)$. Similarly, the second figure of each pair is the orthoradial attenuation factor and is defined by $\max_{\theta} a_{\theta}^{\text{with-atm}}(\tilde{r}, \theta) / \max_{\theta} a_{\theta}^{\text{no-atm}}(\tilde{r}, \theta)$. These attenuation factors are computed for $\tilde{r} \in \{1.059, 1.314\}$.</caption><thead><tr><th rowspan="2"></th><th colspan="4">r̃ = 1.059</th><th colspan="4">r̃ = 1.314</th></tr><tr><th>α = 10<sup>-6</sup></th><th>α = 10<sup>-10</sup></th><th>α = 10<sup>-15</sup></th><th>α = 10<sup>-20</sup></th><th>α = 10<sup>-6</sup></th><th>α = 10<sup>-10</sup></th><th>α = 10<sup>-15</sup></th><th>α = 10<sup>-20</sup></th></tr></thead><tbody><tr><td>Tenuous</td><td>(1.00 - 1.00)</td><td>(1.00 - 0.89)</td><td>(1.00 - 0.15)</td><td>(1.03 - N/A)</td><td>(1.00 - 0.99)</td><td>(1.00 - 0.72)</td><td>(1.01 - 0.11)</td><td>(1.02 - N/A)</td></tr><tr><td>Earth-like</td><td>(1.00 - 1.00)</td><td>(1.00 - 0.76)</td><td>(1.01 - N/A)</td><td>(0.07 - N/A)</td><td>(1.00 - 0.99)</td><td>(1.00 - 0.61)</td><td>(1.02 - N/A)</td><td>(1.03 - N/A)</td></tr><tr><td>Dense</td><td>(1.00 - 0.99)</td><td>(7×10<sup>-7</sup> - N/A)</td><td>(7×10<sup>-7</sup> - N/A)</td><td>(6×10<sup>-7</sup> - N/A)</td><td>(1.00 - 0.98)</td><td>(1.12 - N/A)</td><td>(6×10<sup>-5</sup> - N/A)</td><td>(6×10<sup>-5</sup> - N/A)</td></tr></tbody></table>

that, around $\tilde{r} = 1.059$, the fifth-force was reaching its maximum value for $\alpha$ in the order of $10^{-25}$. Yet in all three atmospheric models under study, the screening of the atmosphere at low altitude occurs for much bigger values of $\alpha$, putting a lower threshold on the maximum accessible fifth-force. When put into perspective with current bounds on $n = 1$ chameleon theory, these results show that the largest part of the unconstrained region maps to a screened atmosphere in the LEO altitude range. All other things remaining equal, the radial component of the fifth force can be recovered by going higher up in altitude, where the atmospheric density is lower.

## IV. Influence on spacecraft trajectory

In this section, we shift our focus to how geodesics get modified in the presence of a putative chameleon fifth-force with respect to the purely Newtonian case. We want to ascertain the effects of the fifth force in a rather quantitative way: is the deviation from Newtonian dynamics large enough to be detected by current satellite technology? Is it possible to discriminate the presence of a fifth force from the imperfect knowledge of the model at stake or small perturbations of the initial conditions? When does a satellite in Low Earth Orbit (LEO) become screened? Besides, we will refrain from commenting too
much on secular drifts that can arise between modified gravity and Newtonian gravity. That is because in any realistic scenario — where many additional forces of different nature come into play —, it would be merely impossible to discriminate the fifth force from such forces. We thus keep our analysis *local*, by focusing our attention on the dynamics at the passage over the mountain.

## A. Screening of the fifth force by the spacecraft

### 1. Existing criteria

We stress that modeling a spacecraft by a material point (in the framework of chameleon gravity) roughly amounts to making the hypothesis that it does not possess a thin-shell. Ref. [4] derives an analytical criterion for a typical satellite in low Earth orbit not to have a thin-shell (see Eq. 80 of this reference). Applying this criterion with the density values employed in our study (except for that of the satellite itself which is set at $8 \times 10^3 \text{ kg/m}^3$) leads straightforwardly to the requirement that $\beta \lesssim 2 \times 10^2 \Leftrightarrow \alpha \gtrsim 3 \times 10^{-20}$ (for $\Lambda = \Lambda_{\text{DE}}$ and $n = 1$). The following orbit propagation results being performed with ($n = 1, \Lambda = \Lambda_{\text{DE}}, \alpha = 10^{-25}$), the satellite would be partially screened according to this criterion and the chameleon effects would thus be smaller than presented.

Ref. [25] derives another criterion based on numerical simulations claiming that the satellite will be fully screened when the thickness of its walls is larger than $100\lambda_{c,\text{wall}}$, where $\lambda_{c,\text{wall}}$ refers to the Compton wavelength in the wall. However, this criterion must be taken with a grain of salt as it was derived for a density contrast $\rho_{\text{vac}}/\rho_{\text{wall}} = 10^{-3}$, far from a realistic setup. Still, applying this second criterion for a wall of thickness 10 cm leads to the fact that the satellite will not have a thin-shell if $\beta \lesssim 2 \times 10^{-2} \Leftrightarrow \alpha \gtrsim 4 \times 10^{-14}$.

### 2. Computation of the field with femtoscope

Although quite qualitative, these two criteria provide us with a comprehension of how the parameters in our model affect the screening of the satellite. For instance, increasing the overall density of the satellite (other things being equal) results in more screening. Ideally, one would compute the scalar field profile sourced by the Earth and the spacecraft all at once — which would avoid having to rely on such criteria and provide a definitive answer. The problem then becomes a numerical one, because the simulation should accommodate a thousand-kilometer-size object (a planet) together with a meter-size object (a satellite). We create a mesh using the Gmsh software [57] that captures both scales (whose ratio is equal or less than $10^{-6}$) thanks to h-adaptivity — a technique that adjusts the mesh resolution by refining or coarsening elements to focus computational resources where they are

most needed. The setup is as follows: we place a cylindrical object centered at coordinates $(\tilde{x}_{\text{Sat}}, \tilde{z}_{\text{Sat}}) = (0, 1.1)$ whose axis is aligned with the z-axis (Fig. 2). The diameter and height of the cylinder are set equal to $L_{\text{Sat}}$ and we denote by $\rho_{\text{Sat}}$ its density. In order to get an order of magnitude of a satellite mean density, we take the example of a CubeSat<fnref n="9" /> whose density is around $\sim 10^3 \text{ kg/m}^3$. We then compute the chameleon field map in the $(x, z)$-plane for various combinations of $\rho_{\text{Sat}}$, $L_{\text{Sat}}$ and $\alpha$. The global acceleration undergone by the cylindrical satellite $\mathbf{a}_{\text{cham}}^{\text{tot}}$ is obtained by integrating the gradient of the scalar field over its whole volume. Under the assumption that the satellite is made of a material of constant density, one gets

$$ \mathbf{a}_{\text{cham}}^{\text{tot}} = -\frac{1}{V} \int_V \nabla \Psi \, dV = -\frac{\beta}{V M_{\text{Pl}}} \int_V \nabla \phi \, dV, \quad (15) $$

where $V = \pi L_{\text{Sat}}^3/4$ is the volume of the cylinder. Now because the setup admits $O_{xz}$ and $O_{yz}$ as planes of symmetry, $a_z^{\text{tot}} = \mathbf{a}_{\text{cham}} \cdot \mathbf{e}_z$ is the only non-zero component of the acceleration vector. Setting $x_{\text{max}} = L_{\text{Sat}}/2$, $z_{\pm} = z_{\text{Sat}} \pm L_{\text{Sat}}/2$, the calculation thus simplifies to

$$ \begin{aligned} a_z^{\text{tot}} &= -\frac{1}{V} \int_0^{2\pi} \int_0^{x_{\text{max}}} \int_{z_{-}}^{z_{+}} x \partial_z \Psi \, dz \, dx \, d\theta \\ &= \left(\frac{2}{L_{\text{Sat}}}\right)^3 \int_0^{x_{\text{max}}} x [\Psi(x, z_{-}) - \Psi(x, z_{+})] \, dx. \end{aligned} \quad (16) $$

The resulting 1D integral can easily be computed using any numerical integration routine.

Fig. 10 shows the chameleon potential $\Psi$ (top row) together with the elementary acceleration $a_z = -\partial_z \Psi$ (bottom row) along the axis of the cylinder that passes through the Earth. On panel a), we recognize the customary chameleon field profile of a screened ball, perturbed nearby $z = 1.1 R_⊕$ by the presence of the satellite. When we zoom-in, we see the potential well imputed to the satellite in panel b). This localized variation of the chameleon field results in a large gradient in absolute value (bigger than anywhere else in the numerical domain). However big the field's gradient may be, looking at panel d) with naked eyes could lead us to believe that it is an odd function with respect to $z = z_{\text{Sat}}$. If that turned out to be the case, then performing the integration (15) would result in a net zero acceleration and the satellite's trajectory would coincide with GR geodesic (in the absence of any non-gravitational perturbation).

We tackle this issue by computing $a_z^{\text{tot}}$ using Eq. (16) for several physical parameters ($\rho_{\text{Sat}}$, $L_{\text{Sat}}$) and several chameleon parameters $\alpha$. From there, the whole question is to determine how the total chameleon acceleration undergone by the satellite compares against that of a

<footnote n="9">⁹ CubeSats have a form factor of 10 cm cubes and have a mass of no more than 2 kg.</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 10. Chameleon potential $\Psi$ (top row) and acceleration $a_z$ (bottom row) along the z-axis. The panels b) and d) are a zoomed version around $z_{Sat} = 1.1 R_⊕$ of panels a) and c) respectively. The parameters used to produce this figure are: $\rho_⊕ = 10^3$ kg/m$^3$, $\rho_{vac} = 10^{-15}$ kg/m$^3$, $\rho_{Sat} = 10^3$ kg/m$^3$, $L_{Sat} = 2 \times 10^{-6} R_⊕ \sim 12.7$ m, $z_{Sat} = 1.1 R_⊕ \sim 7 \times 10^3$ km, $\alpha = 10^{-15}$, $n = 1$, $\beta = 0.24$, $\Lambda = \Lambda_{DE}$. On panels b) and d), the dotted line is centered at $z = z_{Sat}$ while the dashed lines represent the extent of the satellite.</figcaption></figure>

<table><caption>TABLE V. Total chameleon acceleration undergone by a satellite (extended object) $|a_z^{tot}|$ compared to that of a point-like particle $|a_z|$. The accelerations are expressed in m/s$^2$. Each of the three cases corresponds to three different satellites: Case 1 is a benchmark, Case 2 represents a 10 times denser satellite, Case 3 represents a 4 times smaller satellite. As long as the satellite is not screened, $|a_z^{tot}| \simeq |a_z|$. When the satellite is screened (which occurs at a different $\alpha$ depending on the satellite's characteristics), $|a_z^{tot}|$ drops down to nearly zero. Note that Fig. 10 corresponds to Case 1 with $\alpha = 10^{-15}$.</caption><thead><tr><th rowspan="2">α</th><th colspan="3">Case 1 — benchmark</th><th colspan="3">Case 2 — denser</th><th colspan="3">Case 3 — smaller</th></tr><tr><th colspan="3">(ρ<sub>Sat</sub> = 10<sup>3</sup> kg/m<sup>3</sup>, L<sub>Sat</sub> = 2 × 10<sup>-6</sup> R<sub>⊕</sub>)</th><th colspan="3">(ρ<sub>Sat</sub> = 10<sup>4</sup> kg/m<sup>3</sup>, L<sub>Sat</sub> = 2 × 10<sup>-6</sup> R<sub>⊕</sub>)</th><th colspan="3">(ρ<sub>Sat</sub> = 10<sup>3</sup> kg/m<sup>3</sup>, L<sub>Sat</sub> = 5 × 10<sup>-7</sup> R<sub>⊕</sub>)</th></tr><tr><th>|a<sub>z</sub><sup>tot</sup>|</th><th>10<sup>-14</sup></th><th>10<sup>-15</sup></th><th>10<sup>-16</sup></th><th>10<sup>-14</sup></th><th>10<sup>-15</sup></th><th>10<sup>-16</sup></th><th>10<sup>-14</sup></th><th>10<sup>-15</sup></th><th>10<sup>-16</sup></th></tr></thead><tbody><tr><td>|a<sub>z</sub>|</td><td>7.25 × 10<sup>-12</sup></td><td>1.56 × 10<sup>-11</sup></td><td>~ 0</td><td>7.25 × 10<sup>-12</sup></td><td>~ 0</td><td>~ 0</td><td>7.25 × 10<sup>-12</sup></td><td>1.56 × 10<sup>-11</sup></td><td>3.36 × 10<sup>-11</sup></td></tr><tr><td>|a<sub>z</sub>|</td><td>7.25 × 10<sup>-12</sup></td><td>1.56 × 10<sup>-11</sup></td><td>3.37 × 10<sup>-11</sup></td><td>7.25 × 10<sup>-12</sup></td><td>1.56 × 10<sup>-11</sup></td><td>3.37 × 10<sup>-11</sup></td><td>7.25 × 10<sup>-12</sup></td><td>1.56 × 10<sup>-11</sup></td><td>3.37 × 10<sup>-11</sup></td></tr></tbody></table>

point-like particle not affecting the background field. The results set out in Table V provide some answers. We consider three cases which correspond to three satellites with distinct characteristics, namely different length scale and density. For each case, we vary $\alpha \in \{10^{-14}, 10^{-15}, 10^{-16}\}$ and compute the total chameleon acceleration undergone by the satellite (extended object) $|a_z^{tot}|$ and that of a point-like particle $|a_z|$. Surprisingly, the outcome of this experiment is binary:

- When the satellite is unscreened — that is when the scalar field does not reach the value that minimizes the effective potential inside the cylinder $\phi_{Sat}$ — we find that the total chameleon acceleration it undergoes is equal to that of a test particle placed at $z_{Sat}$. This is a remarkable fact, which we did not antic-

ipate by simply looking at Eq. (16) and we thus provide an attempt to explain this phenomenon in Appendix E. In other words, the satellite *feels* the fifth force sourced by the Earth as if it did not perturb the field at all. Consequently, it behaves as a point-like particle and will follow the geodesics of the Jordan frame metric $\tilde{g}_{\mu\nu} = \exp(2\beta\phi/M_{Pl}) g_{\mu\nu}$, where $g_{\mu\nu}$ refers to the Einstein frame metric.

- When the satellite is screened, the integral of the field's gradient over the volume occupied by the satellite vanishes almost completely — essentially because there, the gradient is null. The satellite only *feels* the Newtonian part of the gravitational force and thus follows the geodesics of the Einstein frame metric $g_{\mu\nu}$.
Of course, there actually exists an intermediate case where the satellite would only be partially screened, i.e. where the field would indeed reach $\phi_{Sat}$ deep inside the cylinder while still having some space to vary in its outermost regions. In this specific case, the ratio $|a_z^{tot}|/|a_z|$ lies somewhere between 0 and 1. However, the results reported in Table V suggest that the transition from the unscreened case and the fully screened case does not cover a wide region of the chameleon parameter space. Indeed, taking the Case 1 as an example, the transition occurs between $\alpha = 10^{-15}$ and $\alpha = 10^{-16}$ — refer to Fig. 1 to get a better idea of the narrowness of this region in the chameleon parameter space.

### 3. Discussion

We can check that the reported results are in accordance with the qualitative predictions made by the first two criteria discussed earlier. They both predict that increasing the density and/or the length of the satellite should make it more likely to be screened. This is in agreement with our findings: (i) going from Case 1 to Case 2 shows the effect of an increase by one order of magnitude of the satellite's density, (ii) going from Case 3 to Case 1 illustrates the effect of increasing the satellite's overall size. A follow-up question is whether it is possible to find a distribution of mass inside the satellite $\rho_{Sat}(x, z)$ such that $|a_z^{tot}| > |a_z|$. The simple tests we performed so far — for instance, setting different densities for the upper and lower halves of the satellite — all resulted in $|a_z^{tot}| \le |a_z|$. The question remains open. Additionally, dealing with an extended object means that new rotational degrees of liberty can enter the scene and it would be interesting to look at similar optimization process in order to find the maximum torque (note that Refs. [31, 33] mention this effect and highlight the fact that it can stand out from Newtonian gravity).

Although the satellite model implemented in this section is very simple, this study shows that it is possible for a realistic satellite not to be screened in parts of the chameleon parameter space. This has implications for space-based tests of gravity. For instance, in chameleon models where the scalar field does not couple universally to all matter fields, violations of the weak equivalence principle are not necessarily suppressed by the satellite walls or the experimental setup (as opposed to what was claimed in Ref. [26]). Another example (which holds for a universal coupling constant $\beta$) is that of an accelerometer with a screened test mass onboard an unscreened satellite: the accelerometer would measure a force akin to a bias.

## B. Orbital dynamics of an artificial satellite

Sec. IV A made it clear that in practice, a satellite orbiting some planetary body in the framework of

chameleon gravity cannot be treated as a point-like particle in the entire parameter space. We have highlighted that there is a narrow transition zone beyond which the satellite becomes fully screened and the net fifth force acting on it vanishes almost completely. In what follows however, we make the assumption that we can treat the satellite as a point-like particle. This is justified by at least two reasons:

1. This is a valid approximation in parts of the parameter space (see Table V).

2. One can always, at least at the thought experiment stage, make the satellite smaller or less dense so that it is not subject to screening.

That being said, in all the orbit propagation results presented in the following, we choose the chameleon parameters that produce the strongest fifth force at $\tilde{r} = 1.059$ (which represents an altitude of approximately 376 km) in the absence of atmosphere: ($n = 1, \Lambda = \Lambda_{DE}, \beta = 1.1 \times 10^6$). Notice that this point of the parameter space is already constrained by atom interferometry, see e.g. Ref. [58].

Suppose a point-like particle is placed in a gravitational potential $U$. The equations of motion in spherical coordinates are

$$ 
\begin{aligned}
\ddot{r} - r(\dot{\theta}^2 + \dot{\varphi}^2 \sin^2 \theta) &= -\partial_r U \\
r\left(\ddot{\theta} + 2\frac{\dot{r}}{r}\dot{\theta} - \dot{\varphi}\sin\theta\cos\theta\right) &= -\frac{1}{r}\partial_{\theta}U, \quad (17) \\
r\sin\theta\left(\ddot{\varphi} + 2\frac{\cos\theta}{\sin\theta}\dot{\theta}\dot{\varphi} + 2\frac{\dot{r}}{r}\dot{\varphi}\right) &= -\frac{1}{r\sin\theta}\partial_{\varphi}U
\end{aligned}
$$

where dots refer to time derivatives. The massic energy is given by

$$ \mathcal{E} = \frac{1}{2} (\dot{r}^2 + r^2 \dot{\theta}^2 + r^2 \sin^2 \theta \dot{\varphi}^2) + U \quad (18) $$

and it is conserved along the trajectory, i.e. $\dot{\mathcal{E}} \equiv 0$. Our setup being axisymmetric, we can get rid of the $\varphi$-dependence. Then, note that Eq. (17) implies that the angular momentum $L \equiv r^2 \dot{\theta}$ satisfies

$$ \dot{L} = -\partial_{\theta}U. \quad (19) $$

The problem at stake is a perturbed Kepler problem (the mountain and fifth force contributions are small compared to the central force), whose total gravitational potential $U$ can therefore be decomposed into

$$ U = -\mu/r + \delta U, $$

where $\mu \equiv GM_{body}$ is the standard gravitational parameter of the main body (note that $\mu$ does not encompass the mass contained in the mountain itself). The perturbation $\delta U$ is the sum of the Newtonian potential of the mountain $\delta\Phi$ and the chameleon potential $\Psi$ of the whole system. With this in mind, it is also more appropriate to
decompose the motion into a Keplerian part — that we assume to be circular — and a perturbed part, reading

$$
\begin{align}
r &= a + \delta r & \theta &= \theta_0 + \omega t + \delta \theta \nonumber \\
\dot{r} &= \dot{\delta} r & \dot{\theta} &= \omega + \dot{\delta} \theta \nonumber \\
\ddot{r} &= \ddot{\delta} r & \ddot{\theta} &= \ddot{\delta} \theta \tag{20} \\
L &= L_0 + \delta L & \dot{L} &= \dot{\delta} L \nonumber
\end{align}
$$

In the above, $a$ is the radius of the circular orbit and $\omega$ is the Keplerian pulsation, satisfying $\omega^2 = \mu/a^3$. $L_0$ is the initial angular momentum with $L_0^2 = \mu a$ and $\theta_0$ is the initial co-latitude for a Keplerian motion. This lets us rewrite the equations of motion (17) as

$$ \ddot{r} = L^2/r^3 - \partial_r U, \quad \dot{\delta}\theta = L/r^2 - \omega, \quad \dot{\delta}L = -\partial_\theta U; \quad (21) $$

while the energy conservation reads

$$ (\dot{\delta}r)^2 + (L/r)^2 + 2U - 2\mathcal{E} = 0. \quad (22) $$

Note that at the 0th-order, Eq. (22) boils down to the usual energy conservation in a circular Keplerian orbit

$$ (a\dot{\theta})^2 - \frac{2\mu}{a} - 2\mathcal{E} = 0. $$

### C. Numerical integration with energy conservation

The state vector that we wish to propagate over time is $\mathbf{X} = (\delta r, \dot{\delta} r, \delta \theta, \dot{\delta} \theta) \in \mathbb{R}^4$. It is governed by the ordinary differential equation (ODE) $\dot{\mathbf{X}} = F(t, \mathbf{X})$, where $F: \mathbb{R} \times \mathbb{R}^4 \to \mathbb{R}^4$ is given by Eq. (21). Note that while energy conservation (22) is derivable from the ODE itself, there is no a priori reason for it to hold on the numerical approximation. For one thing, the energy might fluctuate on short time scales depending on the numerical integrator employed, leading to an increase or decrease over longer time scales. Additionally, the r.h.s. of ODE (21) is obtained through FEM computation and is hence noisy, meaning that even so-called energy-preserving integrators would exhibit the energy-drift phenomenon.

Appending the energy conservation (22) to the ODE defines an over-determined differential-algebraic system of equations (DAE). One convenient way to preserve first integrals such as energy conservation when numerically integrating dynamical systems is to resort to projection techniques. The idea behind this class of techniques is to slightly perturb the state after each solver's step so that the energy remains constant. This technique is described in Ref. [59]. As for the implementation, one can easily modify any existing general purpose ODE solver to perform this projection. We provide a minimally modified version of scipy's Runge-Kutta solvers that was used for the numerical integration as supplementary material.

The simulations presented below are performed at $\tilde{r} = 1.059$, which corresponds to an altitude of roughly 376 km. The Newtonian potential and its gradient are

evaluated using the point-mass approximation introduced in Sec. II D 2 as it is hardly distinguishable from the semi-analytical solution. As for the chameleon field, we have the freedom to select an operating point in the parameter space. We choose ($n = 1, \Lambda = \Lambda_{DE}, \alpha = 10^{-25}, \beta = 1.1 \times 10^6$) which we have identified as the point that concurrently results in the strongest fifth force and the greatest field's strength (in the atmosphere-free scenario, see Sec. III A). The experiment performed in Sec. IV A indicates that any medium to large size satellite would presumably be screened in this case. Consequently, the point-like approximation we adopt can be understood as a best case scenario, i.e. an upper bound on the maximum fifth force. Indeed, escaping the screened regime comes at the prize of restricting the allowed range of parameter $\alpha$ to $\alpha > \alpha_{screened}$, which limits the maximum fifth force — see Sec. III A. In terms of initial conditions, the point-like particle is set in a Keplerian motion so that the initial state vector reads $\mathbf{X}(t=0) = \mathbf{0} \in \mathbb{R}^4$, with $\theta_0 = \pi$.

### D. Results and discussion

Here we present and discuss the orbit propagation results. For the sake of clarity and concision, we denote by $\mathbf{X}^{\text{New}}$ and $\mathbf{X}^{\text{Cham}}$ the state vectors in the purely Newtonian case and in the modified gravity (i.e. the sum of chameleon and Newtonian gravity) respectively.

#### 1. Results of the simulations

The evolution with respect to time of the main quantities of interest are presented in Fig. 11. The time spans 10 hours which encompasses roughly three full orbits. The first row of this figure shows how, in a purely Newtonian setting, the presence of the mountain breaks the Keplerian, circular motion. Some elements, such as $\delta L$ are very correlated to the passage of the point-mass above the mountain (denoted by the vertical light-gray dashed lines on each panel). The angular momentum is indeed roughly constant along the trajectory, except nearby $\theta \sim 0$ where it peaks very sharply ($L_0$ being negative, this corresponds to an increase in the absolute value of $L$). On the other hand, some other elements are irreversibly imprinted by the mountain after the first passage above it, see e.g. $\delta r, \dot{\delta} r$ or $\delta \theta$, leaving traces on the longer term. The physical intuition for this is that, although the gravity field is symmetric with respect to $\theta = 0$, the dynamics is not. Indeed, in the $(\theta > 0)$-plane, the system acquires non-zero velocity $\dot{\delta} \theta$ which leads $\delta \theta$ to slowly drift away from zero initial state. Right after the passage of the mountain — that is when $\theta$ becomes negative — the mountain's gravity acts as a restoring force, which has the immediate effect of slowing down $\delta \theta$. But it is already too late: in the meantime, the altitude has been disturbed ($\delta r \neq 0$) and $\theta$ continues on its run (at
<figure><img src="image_12.png" /><figcaption>FIG. 11. Orbit propagation over three Keplerian periods. The first row shows the evolution of the state vector $\mathbf{X}^{\text{New}} = (\delta r^{\text{New}}, \dot{\delta} r^{\text{New}}, \delta\theta^{\text{New}}, \delta L^{\text{New}})$ and $\dot{\delta}\theta^{\text{New}}$ with respect to time, where the dynamics is purely Newtonian. The second row lays emphasis on the orbital dynamics in modified gravity by showing $\mathbf{X}^{\text{Cham}} - \mathbf{X}^{\text{New}}$ and $\dot{\delta}\theta^{\text{Cham}} - \dot{\delta}\theta^{\text{New}}$. The vertical light-gray dashed lines correspond to the instants at which the point-like particle passes over the mountain, at $\theta = 0$, in the purely Newtonian case.</figcaption></figure>

an increased $\omega + \dot{\delta}\theta$ pace), so that the restoring force at $-\theta_* < 0$ is not equal to the force that disturbed the Keplerian motion at $\theta_*$. Once the symmetry is broken, the orbit can no longer be circular — it has a non-zero osculating eccentricity — which is why $(\delta r, \delta r', \delta\theta, \delta\theta')$ exhibit an oscillatory behavior at approximately the Keplerian frequency. We dedicate Appendix F to prove this point in a more rigorous way.

In the second row of Fig. 11, we illustrate what we call the *anomaly* $\mathbf{X}^{\text{Cham}} - \mathbf{X}^{\text{New}}$, that is simply the difference between the geodesic in modified gravity and in Newtonian gravity — for the same set of initial conditions. Surprisingly, apart from $\delta\theta^{\text{Cham}} - \delta\theta^{\text{New}}$ which undergoes a steady decline, the other elements of the anomaly seem to be periodic and remain around zero. In Fig. 12, we show the slow drift of the distance anomaly between the two trajectories, that is $\|\mathbf{r}^{\text{Cham}} - \mathbf{r}^{\text{New}}\|$. This steady increase of the distance has a mean slope of $\sim 2 \text{ m/h}$, but the rate of increase is maximized at each passage of the point-mass above the mountain where it exceeds $4 \text{ m/h}$.

All these orders of magnitude relating to the anomaly should be put into perspective with the current level of precision with which we are able to determine a satellite's position and other orbital elements. This process goes under the name 'Precise Orbit Determination' (POD) and involves analyzing various observational data, often obtained from ground-based tracking stations or satellite-based instruments — see e.g. Refs. [60–62] for the implementation of these techniques and the reachable orders of magnitude in terms of precision. One of the main space geodetic techniques is Satellite Laser Ranging (SLR) which measures the time it takes for a laser beam to travel from the ground station to a retro-reflector on the satellite and back again, providing unambiguous range measurements to millimeter precision

<figure><img src="image_13.png" /><figcaption>FIG. 12. Distance anomaly as a function of time. The various passages above the mountain, depicted by the vertical light-gray dashed lines, correspond to the most rapid increase in this distance.</figcaption></figure>

[63, 64]. This technique is also placed at the service of fundamental physics; in this respect let us mention the recent launch of the LARES 2 satellite [14] to test GR.

Therefore, with no uncertainty on the model and initial conditions, the anomaly caused by the fifth force is of the order of a meter (see the leftmost column of Fig. 11), which is around three orders of magnitude larger than the best attainable precision. At this stage of the discussion, it would seem easy to detect the fifth force.
## 2. The GRACE-FO scenario

The GRACE-FO<fnref n="10" /> mission, currently in operation, aims at monitoring the Earth's gravitational field. It uses a pair of satellites flying on the same orbital path, approximately 220 km apart. As they orbit the Earth, the spacecraft are affected by the uneven gravity field caused by the uneven distribution of mass inside the planet — e.g. the presence of a mountain, which produces a slightly stronger gravitational pull. As a result, the distance between the two satellites varies continuously over time. This distance variation is measured down to the micron level thanks to a microwave ranging system<fnref n="11" /> [66]. Ultimately, the changes in the distance between the satellites are used to monitor the time variations of the Earth gravity field due to mass changes (ice melting, droughts, floods, etc.).

Given the extreme level of precision GRACE-FO is able to reach in terms of ranging, we investigate whether or not fifth force effects would end up being in its sensitivity range. To do so, we simulate a pair of satellites following each other by duplicating the trajectory and shifting it in time by a few minutes, mimicking the real mission configuration. We can then reconstruct the change in inter-spacecraft distance with respect to time $d(t)$. An example of such a curve is given in Fig. 13 (red solid line or salmon dashed line, the two being indistinguishable by eye), where roughly three orbits have been completed. The passage above the mountain can easily be spotted on the curve by the little spikes they spawn. They can be understood fairly intuitively: approaching the mountain's latitude, the leading satellite starts feeling a slightly stronger gravity relative to the trailing one and is pulled slightly ahead, increasing the distance between the two satellites. When the first satellite has eventually passed on the other side of the mountain (that is $\theta < 0$), it is slowed down while the trailing satellite is accelerating, resulting in a decrease in the inter-satellite distance overall. Long after the occurrence of this short-term event, this distance continues to vary in sinusoidal fashion. This is due to the fact that, as discussed above and brought to light in Fig. 11, the orbit is no longer circular after the passage of the mountain and therefore the velocity varies along an orbit.

The difference between the modified gravity case $d_{\text{Cham}}$ and the Newtonian case $d_{\text{New}}$ can hardly be seen on those curves. It is depicted by the solid blue line on Fig. 13 (again called the 'anomaly'). Choosing the same set of initial conditions for both models ensures that the anomaly is null at time $t = 0$. The anomaly exhibits three maxima — at a level of a few centimeters — cor-

<figure><img src="image_14.png" /><figcaption>FIG. 13. Left y-axis: inter-spacecraft distance with respect to time in Newtonian gravity (solid red line) and in modified gravity (dashed salmon line). Right y-axis: the blue curve corresponds to the anomaly, that is the difference between the two models. The initial time-delay between the two satellites is set to 100 seconds.</figcaption></figure>

responding to the passage of the pair of satellites above the mountain.

In view of these results, we may believe that chameleon-like fifth forces should be detectable with our current space technology as the anomaly is ~ $10^4$ times larger than the sensitivity threshold of GRACE-FO. That would be true under the (unrealistic) assumptions that:

1. the initial conditions are perfectly known, that is there is no uncertainty in our initial state $\mathbf{X}_0$ prior to the propagation;

2. the density model of the main body (the Earth) is perfectly known.

Neither of the two hypotheses can be fulfilled in practice. In the two forthcoming sections, we tackle these points and strongly mitigate our previous statement in regard to fifth force detectability in space.

## 3. Perturbation of initial conditions

Here, we investigate whether a slight modification of the initial state vector $\mathbf{X}_0 \leftarrow \mathbf{X}_0 + \delta\mathbf{X}_0$ could account for the anomaly that we unveiled in Fig. 13. For that purpose, we employ a Nelder-Mead optimizer where our objective function is

$$g: \mathbb{R}^4 \to \mathbb{R}_+ \qquad (23)$$

$$\delta\mathbf{X}_0 \mapsto \|d_{\text{Cham}}(\mathbf{t}, \mathbf{X}_0) - d_{\text{New}}(\mathbf{t}, \mathbf{X}_0 + \delta\mathbf{X}_0)\|_2,$$

<footnote n="10">¹⁰ Gravity Recovery and Climate Experiment Follow-On.</footnote>

<footnote n="11">¹¹ GRACE-FO also employs laser-ranging interferometry for a more accurate inter-satellite ranging which can improve the separation distance measurement by a factor of more than 20 relative to the GRACE mission [65].</footnote>
where $\|\cdot\|_2$ is the two-norm in $\mathbb{R}^4$ and $\mathbf{t} = [t_0, t_1, \dots, t_N]$ is the discrete time vector with $t_N \sim 10$ h and $N = 10^4$. We find an optimum at

$$ \delta\mathbf{X}_0^{\text{opt}} = \begin{bmatrix} -6.04 \times 10^{-1} & \text{m} \\ +6.16 \times 10^{-7} & \text{m/s} \\ -7.49 \times 10^{-7} & \text{rad} \\ +1.98 \times 10^3 & \text{m}^2/\text{s} \end{bmatrix}, \quad (24) $$

with a residual smaller than 5 mm<fnref n="12" />. This is an extremely good fit given the characteristic length of the problem (several hundred kilometers). The conclusion to be drawn from this is clear: on this specific inter-satellite distance tracking example, an extra chameleonic acceleration cannot be distinguished from a small perturbation of the initial state vector. A brief analysis indicates that the parameter that has the biggest weight in Eq. (24) is $\delta r_0$. Now the question is how this small perturbation compares to the precision with which we have access to the initial state. As discussed previously in Sec. IV D 1, it turns out that the initial radial distance $r_0$ could in principle be determined with at most centimetric precision which is smaller than the 60 cm perturbation found in Eq. (24). Although this does not constitute a rigorous proof, this brief study tends to indicate that the ‘unknown initial state’ hypothesis can be ruled out.

### 4. Perturbation of the mass distribution

Nonetheless, the knowledge of initial conditions is not the only potential source of degeneracy. Indeed, the mass distribution inside the main body — the very source of gravity — is perhaps the most important degree of freedom to have knowledge of. In that perspective, can the fifth force effects on a satellite be interpreted in the framework of Newtonian gravity as a slightly altered density model? In order to answer that question, we continue in the same spirit as in Sec. IV D 3 by constructing an optimization problem. We saw earlier on, notably in Fig. 4, that the Newtonian potential of the mountain could very well be approximated by a point-mass. We can thus try and perturb the density model — and consequently the Newtonian potential — by adding a point-mass somewhere along the z-axis (see Fig. 2), as we do not wish to break the azimuthal symmetry. This simple model has only two parameters:

- $m_*$ the mass of the point-mass;

- $z_*$ the z coordinate of the point-mass.

<table><caption>TABLE VI. Best fit parameters of the approximation of the chameleon acceleration by a point-mass in Newtonian gravity.</caption><thead><tr><th></th><th>r&#773; = 1.059</th><th>r&#773; = 1.111</th></tr></thead><tbody><tr><td>m&#773;<sub>*</sub></td><td>6.6 &times; 10<sup>-6</sup></td><td>1.8 &times; 10<sup>-6</sup></td></tr><tr><td>z&#773;<sub>*</sub></td><td>1.03</td><td>1.06</td></tr><tr><td>f(m&#773;<sub>*</sub>, z&#773;<sub>*</sub>)/f(0, 1)</td><td>3.1 &times; 10<sup>-2</sup></td><td>5.6 &times; 10<sup>-2</sup></td></tr><tr><td>||&part;&theta;(&Phi;<sub>*</sub> &minus; &Psi;)||<sub>L<sup>2</sup></sub>/||&part;&theta;(&Phi;<sub>*</sub> + &Psi;)||<sub>L<sup>2</sup></sub></td><td>0.14</td><td>0.22</td></tr><tr><td>||&part;r(&Phi;<sub>*</sub> &minus; &Psi;)||<sub>L<sup>2</sup></sub>/||&part;r(&Phi;<sub>*</sub> + &Psi;)||<sub>L<sup>2</sup></sub></td><td>0.07</td><td>0.08</td></tr></tbody></table>

The goal is then to find the pair $(m_*, z_*)$ for which Newtonian gravity best mimics the modified gravity case. Precisely, our objective function is

$$ \begin{align} f: \mathbb{R}^2 \to \mathbb{R}_+ \tag{25} \\ (m_*, z_*) \mapsto \int_0^\pi (\partial_\theta \Phi_* - \partial_\theta \Psi)^2 \, d\theta \nonumber \\ \phantom{(m_*, z_*) \mapsto \int_0^\pi (\partial_\theta \Phi_* - \partial_\theta \Psi)^2 \, d\theta} + \int_0^\pi (\partial_r \Phi_* - \partial_r \Psi)^2 \, d\theta \nonumber \end{align} $$

where $\Phi_*$ is the Newtonian potential created by the extra point-mass and the integral is carried out at fixed $\tilde{r}$. We denote by $(m_r^*, z_r^*)$ the pair that minimizes the function $f$ at radius $\tilde{r}$. Using $\|\cdot\|_{L^2}$ to denote the $L^2$-norm over the space of square-integrable function on $[0, \pi]$, one can rewrite

$$ f(m_*, z_*) = \|\partial_\theta(\Phi_* - \Psi)\|_{L^2}^2 + \|\partial_r(\Phi_* - \Psi)\|_{L^2}^2. $$

Basically, we aim at approximating both the radial and orthoradial parts of the chameleon acceleration at the same time. This optimization problem being low-dimensional, we can dispense with a sophisticated optimization algorithm and do a full exploration of the parameter space instead (see Fig. 14). Note that our point-mass model cannot reproduce the chameleon monopole (which is, in other words, a constant radial acceleration offset). Therefore, we removed it by hand before proceeding to the optimization phase. This offset is tiny: $\sim 1.4 \times 10^{-7} \text{ m/s}^2$ which corresponds to relative change of the mean density of the main body of only $\sim 9 \times 10^{-8}$. In comparison, let us mention that the Earth mass is known with a relative uncertainty of $10^{-4}$.

The results for $\tilde{r} \in \{1.059, 1.111\}$ are reported in Table VI, where $\tilde{m}_r^* = m_r^*/M_{\text{mountain}}$ and $\tilde{z}_r^* = z_r^*/R_{\text{body}}$.

To further assert the quality of the fit in quantitative terms, we compute the ratio $f(m_*, z_*)/f(0, 1)$ (where $\{m_* = 0, z_* = 0\}$ corresponds to flat profiles) as well as relative errors in $L^2$-norm. Several comments must be made:

- With only a simplistic model (i.e. a single point-mass has been added to the pre-existing model, contributing to the global Newtonian potential), we manage to approximate the fifth-force profile at a

<footnote n="12">¹² The last entry of vector $\delta\mathbf{X}_0^{\text{opt}}$ in Eq. 24 corresponds to the perturbation in the initial angular momentum and may attract attention due to the fact it is orders of magnitude bigger than the other entries. To provide a benchmark, the unperturbed initial angular momentum is $L_0 \simeq 2.2 \times 10^{10} \text{ m}^2/\text{s}$.</footnote>
<figure><img src="image_15.png" /></figure>

FIG. 14. a) Contour plot of the objective function (in log scale) in the ($m_*, z_*$)-plane$^a$. b) Contour plot of the anomaly in the ($m_*, z_*$)-plane. The white cross and blue circle are located at the objective function's minimum when $\tilde{r} = 1.059$ and $\tilde{r} = 1.111$ respectively. The area left in plain white is not physically accessible as it corresponds to a 'negative extra mass' in vacuum. Warning: looking at the anomaly (bottom panel), one might expect the special case $m_* = 0$ to reduce to the case displayed on Fig. 13 and exhibit an anomaly of a few centimeters. The difference lies in the fact that here, the chameleon radial acceleration offset is artificially reproduced by slightly increasing the main body's mass.

given altitude with remarkable accuracy (see the various figures in Table VI).

- This approximation is good enough to almost reproduce the dynamics of the satellite's orbit over the mountain. In fact, we can repeat for instance the same exercise as we did in Fig. 13 and compute the so-called *anomaly*, i.e. the difference between the “modified gravity without extra mass” case and the “Newtonian gravity with extra mass” case. We find it to be no greater than 15 µm. This is more than a thousand times smaller than the anomaly

computed in Fig. 13. This invites us to moderate the statements made earlier, since we are approaching here the precision limits of the GRACE-FO's LRI system.

- The objection could be made that the characteristics of the point-mass associated with the best fit do not correspond to any physical reality. Indeed, taking the second column of Table VI with entries ($\tilde{m}_* = 6.6 \times 10^{-6}$, $\tilde{z}_* = 1.03$) suggests that there would be a ~2 × 10$^{12}$ kg mass at an altitude of 186 km (or equivalently, 123 km above the mountain's top) — which is obviously absurd! Nevertheless, it can be seen on the top panel of Fig. 14 — which represents the cost function (25) in the ($m_*, z_*$)-plane — that lowering a bit $z_*$ from the optimum (depicted by the white cross) while maintaining $m_*$ constant has only a slight effect on the cost function. For this reason, the dynamics is not much affected by a shift of $z_*$ towards the planet. As a matter of fact, setting $z_* = 1$ (i.e. bringing the extra mass at the planet's surface) leads to an anomaly bounded below 40 µm. The bottom panel of that same figure is intended to illustrate this phenomenon, and the strong correlation between the cost function and the anomaly is visible to the naked eye.

We can even place this extra mass at the same location as the point-mass Newtonian approximation of the mountain itself (see caption of Fig. 4) without any major change in the dynamics. This can therefore be interpreted as slightly increasing the mountain's density, by roughly 10$^{-3}$%. Such a slight deviation could equivalently be attributed to the fact that the gravitational constant $G$ is only known with some certainty with four significant digits [67, 68].

These orders of magnitude on the density must be put into perspective with our current knowledge of the Earth inner density, with all the attendant uncertainties. Despite advancements in geophysical techniques, our knowledge of mass distribution is still imperfect, for simple reasons:

1. The planet's interior is out of reach. As a matter of fact, the deepest human-made hole ever dug is only 12.3 km deep (less than 0.2% of the Earth radius).

2. The density is not uniform (even at fixed depth), which means extrapolation is not a valid procedure unless strong assumptions are made about the Earth's composition and structure.

3. As it happens, we also have to rely on indirect measurements, ranging from gravitational anomalies and magnetic anomalies [69–71] to seismic analysis [72–74]. All these techniques are in turn limited in both resolution and accuracy.

On this latter point, we stress that one should be careful when trying to put constraints on a given modified grav-

<footnote>$^a$ Note that negative mass (which would correspond to an *extrusion* for $\tilde{z}_* < 1 + h_m$) cannot represent the chameleon acceleration as well as positive mass.</footnote>
ity model, using a model of the Earth that comes from gravitational measurements in the first place. Indeed, the inversion of a gravity map into, say, a density map is model dependent (and unless contraindicated, would have been performed in a Newtonian framework). See Ref. [44], where this topic is discussed at length. In this regard, let us mention a recent work [75] that proposes to use the preliminary reference Earth model (PREM) [72], which is a radial seismic model, to constrain some alternative theories to GR.

## 5. Breaking the degeneracy

We have seen with two simple examples that drawing a distinction between a fifth force and model uncertainties (of different natures) is no easy task. These uncertainties spearhead degeneracies, which we partially address here.

In Sec. IV D 3, we provided the relevant orders of magnitude of the perturbation of a satellite initial state vector necessary to alone mimic a fifth force influence. The perturbation on the initial altitude was then put in comparison against the available level of precision for LEO satellites. It turned out POD techniques are good enough to the relatively large perturbation found in Eq. (24). Yet one must bear in mind that this was done on a very specific test-case and the conclusion may not generalize to others.

In Sec. IV D 4, we looked at how to distinguish a chameleon acceleration on top of Newtonian dynamics from a slight change of the density model in a purely Newtonian framework. We showed that it was possible to imperceptibly tweak the mass distribution in the mountain and in the planet to reproduce the chameleon acceleration profile at a given altitude (see Table VI). This naturally raises the question of whether such a fit works at different altitudes, or rather how well. Elements of response can be found in Table VI and Fig. 14. The first rows of Table VI bring out the fact that inferring the mountain's characteristics at two orbital radii leads to two clashing physical realities: strikingly, the extra mass inferred at $\tilde{r} = 1.059$ is almost 4 times greater than it appears at $\tilde{r} = 1.111$. In order to sharpen the analysis, we represent in Fig. 14 by a white cross and a blue circle the cost function's minimum at $\tilde{r} = 1.059$ and $\tilde{r} = 1.111$ respectively, while the contour plots are performed for $\tilde{r} = 1.059$. Following the notations introduced above, $(m_*^{1.059}, z_*^{1.059})$ and $(m_*^{1.111}, z_*^{1.111})$ are coordinates of the white cross and the blue circle respectively. Similarly, $f_{1.059}$ and $f_{1.111}$ refer to the cost functions at the two altitudes. At $(m_*^{1.111}, z_*^{1.111})$, we see that the cost function $f_{1.059}$ is much above its minimum and, in turn, corresponds to a large anomaly. Quantitatively speaking, we have

$$ \frac{f_{1.059}(m_*^{1.111}, z_*^{1.111})}{f_{1.059}(m_*^{1.059}, z_*^{1.059})} \approx 7.8 \times 10^2, $$

which is a big ratio and reflects that $(m_*^{1.111}, z_*^{1.111})$ does

not produce a good fit of the chameleon acceleration profile at $\tilde{r} = 1.059$. On the other hand, changing our perspective to $f_{1.111}$, we have

$$ \frac{f_{1.111}(m_*^{1.059}, z_*^{1.059})}{f_{1.111}(m_*^{1.111}, z_*^{1.111})} \approx 13, $$

meaning that $(m_*^{1.059}, z_*^{1.059})$ is better tolerated by $f_{1.111}$ and $f_{1.059}$ than $(m_*^{1.111}, z_*^{1.111})$.

Fig. 15 provide more visual insights into these tensions. As in panel b) of Fig. 14, we computed the anomaly as a function of the pair $(m_*, z_*)$, in a scenario where the two GRACE-FO-like satellites orbit at $\tilde{r} = 1.059$ (associated with blue colors) and in another scenario where they orbit at $\tilde{r} = 1.111$ (associated with orange colors). On panel a), we display two contours corresponding to anomaly thresholds of $3 \times 10^{-6}$ m and $7 \times 10^{-6}$ m, for both altitudes. The less the blue contours overlap with the orange ones, the greater the tension. Panels b) and c) complement the figure by representing the anomaly along the dotted lines visible on panel a) which pass through the minimal anomaly for each altitude. Ideally, the performances showcased by the GRACE-FO laser-link technology would allow for an exclusion of any $(m_*, z_*)$ pair mapping to an anomaly greater than a micrometer, revealing the incompatibility between the two density models.

In conclusion of this section, the use of different altitudes in the analysis is a first step toward breaking the degeneracy. This idea was already put forward in Ref. [44] where the authors study the impact of a Yukawa potential on the spherical harmonic coefficients of the Earth. Precisely, the rescaled coefficients $y_{lm}$ (introduced in Sec. II C 1) become dependent on the altitude meaning for instance that measurements of the $J_2$ zonal term at GOCE and GRACE altitudes could provide a test of the model<fnref n="13" />. On the whole, the difficulty lies in being able to find a set of several physical measures that would be in tension one with another when adding a fifth force to the play:

- The greater the tension, the tighter the potential constraints on the modified gravity model.

- The more measurements we have, the greater the likelihood of ending-up with a significant tension.

We tried to bring out such a tension with the GRACE-FO setup deployed at two different altitudes (see Table VI and Fig. 14). Nevertheless, this does not look practical in actual experiments. Even though we manage to create a small tension in the anomaly, one must bear in mind that our fitting model consisting of a single point-mass remains overly simplistic. More complex (and viable) density models may relieve this tension very well. In

<footnote n="13">¹³ Note that this would also be true for the chameleon model as $y_{lm}^C$ depends on the radial coordinate $r$.</footnote>
<figure><img src="image_16.png" /></figure>

FIG. 15. Tensions in the inferred mountain's characteristics. The blue elements refer to the altitude $h_1 = 376$ km ($\tilde{r} = 1.059$) and the orange elements refer to the altitude $h_2 = 707$ km ($\tilde{r} = 1.111$). The anomaly observed in the inter-satellite distance at $h_1$ (resp. $h_2$) is best explained in the framework of Newtonian gravity by an extra point-mass with characteristics ($m_*^{1.059}$, $z_*^{1.059}$) (resp. ($m_*^{1.111}$, $z_*^{1.111}$)) depicted by the black cross (resp. gray circle) in panel a)ª, which minimizes the anomaly down to $1.3 \times 10^{-6}$ m (resp. $2.6 \times 10^{-7}$ m). In panel a), the darker contours map to an anomaly below $3 \times 10^6$ m while the lighter ones map to an anomaly below $7 \times 10^{-6}$ m. The gray shaded area is not physically accessible as it corresponds to a 'negative extra mass' in vacuum. Panel b) (resp. c)) represent the anomaly along $\tilde{m}_*$ at $\tilde{z}_*^{\tilde{r}}$ (resp. $\tilde{z}_*$ at $\tilde{m}_*$). The anomaly peaks (maxima) visible on panel c) correspond to the horizontal feature seen previously in Fig. 14-b) around $\tilde{z}_* = 1.06$.

short, the ultimate goal of constraining the chameleon model with space-based geodesy is impeded by:

- the uncertainty on the source of gravity, that is the density model;

- the uncertainty on the measurements themselves;

- all the other forces acting on a satellite, ranging from atmospheric drag and solar radiation pressure, to third-body perturbation, which also come with error bars in our models. Note that these perturbing forces are also a nuisance for geodesy, hence the use of accelerometers on board satellites [76, 77].

In Appendix G, we further look at orbital periods at two different altitudes and compute their difference, in a Newtonian framework and in a modified gravity framework.

V. Discussion and conclusion

Effect of a mountain

This article investigated the testability of chameleon gravity by space geodesy experiments, with a focus on the influence of the local landform and the atmosphere. The motivations were twofold. First, viable regions of the chameleon parameter space all map to a screened Earth, that is only a thin-shell contributes to the fifth force. Therefore, it seemed important to study departures from spherical symmetry, hereby embodied by a mountain. Second, while published works sometimes account for the atmosphere in their study, the models implemented are simplistic (often one layer of constant density) and the determination of whether it has a thin-shell is based on

<footnote>ª Note that the pairs ($\tilde{m}_*^{\tilde{r}}$, $z_*^{\tilde{r}}$) are different from the ones reported in Table VI for two reasons: (i) here, we minimize the anomaly in the inter-satellite distance between the “modified gravity without extra mass” case and the “Newtonian gravity with extra mass” case, which is different from minimizing the objective function *f* given by Eq. (25); and (ii) we allowed ourselves to modify the initial inter-satellite distance for the two considered altitudes in order to better showcase the tension.</footnote>
rather qualitative arguments. Addressing such questions
is not possible by means of analytical techniques alone
due to the complexity of the physical models we wished
to study, and to the nonlinear nature of the chameleon
equation of motion. We thus resorted to numerical sim-
ulations — performed with the code *femtoscope* — to
conduct this work.

We obtained the chameleon contribution to the total gravitational potential of a mountainous planet, scanning through an extended region of the parameter space. As already pointed out in Ref. [29], the unscreened regime shares similarities with pure Newtonian gravity in that, in both cases, the fields are sourced by the entire mass of the main body. Consequently, the chameleon potential in the unscreened regime is roughly the same as the Newtonian potential up to an affine transform (and the same goes for the accelerations). As we enter the screened regime however, the multipole expansion of the chameleon field starts to depart from that of the Newtonian potential, revealing a distinct signature. In terms of acceleration, the chameleon acceleration vector is a bit more directed towards the mountain compared to the Newtonian acceleration. Their norm ratio remains small though, bounded from above by ~$10^{-6}$ at the equivalent of LEO altitudes in the atmosphere-free case<fnref n="14" />.

Effect of an atmosphere

Based on our study of three distinct atmospheric den-
sity profiles, we found that the addition of an extra layer
of air surrounding the main body can mitigate the ef-
fect of the fifth force. We showed that there exists a
threshold on the value of the parameter $\alpha$. Above this
threshold, the atmosphere acts as an attenuator, effec-
tively reducing the chameleon acceleration by a certain
amount compared to the case without atmosphere. Be-
low this cutoff, the effect of the atmosphere is more dras-
tic: any non-radial dependence of the scalar field vanishes
— the mountain is plainly invisible. For even smaller val-
ues of $\alpha$, the atmosphere itself becomes screened, and the
chameleon field is thereupon fully determined by the at-
mospheric density profile. As we saw, it is even possible
in this case to enhance the radial fifth force at given al-
titude with respect to the atmosphere-free case. This
study represents a step forward with respect to previous
work discussing the influence of the atmosphere. More-
over, this clearly gives the edge to bodies devoid of at-

mosphere when it comes to select a Solar System site for
testing this screened scalar-tensor model.

Space geodesy experiments

Our knowledge of the geopotential comes to a large extent from spaceborne geodesy. From this standpoint, we thus investigated whether constraints could be put on modified gravity models using satellites in orbit. For that purpose, we performed orbit propagations, with and without the putative fifth force, and studied the resulting anomaly<fnref n="15" /> on several observables (such as the variations of the distance between two satellites following each other as in the GRACE and GRACE-FO setups). While the anomalies we find are technically well within the detection range of current on-board and ground-based space-technology, we showed that uncertainties in the model for the distribution of matter are large enough to allow for degeneracies.

We laid emphasis on the fact that one way to distin-
guish a chameleon acceleration from a slight change of
the Earth density model in a purely Newtonian frame-
work is to rely on experiments performed at (at least)
two different altitudes. Indeed, in the regime where the
Earth (or any other planetary body) is screened, the
chameleon acceleration does not decrease as $r^{-2}$ like the
Newtonian acceleration (this is particularly stressed in
Refs. [44, 80]). If the chameleon field *actually exists*, then
inferring density models under the assumption of Newto-
nian gravity at several altitudes should result in tensions
between those models. Of course, these tensions should
be accounted for in a probabilistic way, which is beyond
the scope of this article. Conversely, if it were not for
all the other perturbing forces that greatly complexify
the model, this method could be used to put constraints
on the chameleon model, and more generally on massive
scalar-tensor theories.

Back-reaction of a satellite on the scalar field

We also took into account the back-reaction of an ob-
ject as small as a satellite in orbit on the scalar field.
For the first time, we went beyond the various approx-
imations found in the literature and computed the *full
solution* of the {Earth + Satellite} system using *femto-
scope*. This involves taking advantage of the h-adaptivity
technique granted by FEM. We could then compute the
overall fifth force acting on an object with a simple ge-
ometry and characteristics close to that of a real space-
craft (length-scale and density). Surprisingly, as long as

<footnote n="14">¹⁴ It is insightful to compare this ratio with the ratio of the Solar radiation pressure over the Newtonian acceleration, which is around $10^{-8}$ [78]. Despite being so tiny, the Solar radiation pressure perturbing acceleration, when integrated over many orbits, is enough to cause significant drifts of orbital elements [79]. What makes the chameleon acceleration difficult to distinguish from the Newtonian acceleration is the fact that they are both sourced by the same body.</footnote>

<footnote n="15">¹⁵ The term ‘anomaly’ is used to refer to the difference for a given observable between the {Newtonian gravity} case and the {Newtonian gravity + fifth force} case.</footnote>
the satellite is not screened and despite the background field being disturbed, the global chameleon acceleration undergone by the satellite is the same as the one a point-particle (not disturbing the background field) would experience. We provide mathematical insights into why this is the case in Appendix E. In the screened regime however, the net fifth force vanishes to zero. The transition between those two regimes occurs over a narrow band in the chameleon parameter space.

## Outlook

While we focused on fifth force searches, other venues exist to test scalar-tensor theories. For instance, in any such theory involving a conformal coupling of the scalar field $\phi$ to matter fields in the Einstein frame, the gravitational redshift effect has a $\phi$-dependence (see e.g. Ref. [4] for the chameleon's contribution to this effect). We will take a deeper look at this effect in an upcoming article. Most importantly, we will tackle the question of what is actually measurable, and with which precision. As a complement to fifth force searches where we look for dynamical effects whose amplitude depends inherently of the field's gradient, gravitational redshift (or equivalently, gravitational time-dilation) can be measured in a static configuration and is sensitive to the field's strength. Whether clocks are put into orbit (as envisaged in the ACES mission [81]) or left on Earth, they have become so precise<fnref n="16" /> that their constraining power (i.e. the possibility to use this technology to rule out modified gravity models) has to be quantified. The bound given in Ref. [4] has to be revisited, given two decades elapsed since the writing of this article.

## Acknowledgments

HL thanks Pablo Richard for fruitful discussions about numerical computations. We acknowledge the financial support of CNES through the APR program ("MICRO-SCOPE" project).

### A. Conversion of cosine and sine coefficients to bare coefficients

Let $f: S^2 \to \mathbb{R}$ be a real-valued function on the unit sphere and $L \in \mathbb{N}^*$ a maximum spherical harmonic degree. The truncated spherical harmonic expansion, which defines an approximation $f_{\text{trunc}} \simeq f$, may be written as

$$f_{\text{trunc}}(\mathbf{n}) = \sum_{m=0}^{L} \sum_{l=m}^{L} \left[ C_{lm} \bar{P}_{lm}(\cos \theta) \cos(m\varphi) + S_{lm} \bar{P}_{lm}(\cos \theta) \sin(m\varphi) \right] \quad (\text{A1})$$

with $\mathbf{n} = (\theta, \varphi)$ and $\bar{P}_{lm}$ the normalized associated Legendre functions<fnref n="17" /> which relate to their unnormalized counterparts $P_{lm}$ via

$$\bar{P}_{lm}(x) = \sqrt{\frac{(2 - \delta_{m0})(2l + 1)(l - m)!}{4\pi(l + m)!}} P_{lm}(x). \quad (\text{A2})$$

We want to convert the cosine and sine coefficients ($C_{lm}, S_{lm}$) into the bare coefficients $f_{lm}$ that appear in the usual expansion

$$f_{\text{trunc}}(\mathbf{n}) = \sum_{l=0}^{L} \sum_{m=-l}^{+l} f_{lm} Y_{lm}(\mathbf{n}). \quad (\text{A3})$$

In order to express ($C_{lm}, S_{lm}$) as a function of $f_{lm}$, we start from Eq. (A3) and interchange the order of summations over $l$ and $m$ to obtain

<footnote n="16">¹⁶ High-precision clocks, such as optical lattice clocks, currently achieve astonishing levels of accuracy with a fractional frequency uncertainty of approximately $10^{-19}$ to $10^{-20}$ [82, 83].</footnote>

<footnote n="17">¹⁷ In this respect, the definition of normalized associated Legendre functions is consistent with the definition of orthonormalized spherical harmonic functions. See Table 1 from Ref. [43].</footnote>
$$
\begin{align*}
f_{\text{trunc}}(\mathbf{n}) &= \sum_{l=0}^{L} \sum_{m=-l}^{+l} f_{lm} Y_{lm}(\mathbf{n}) \\
&= \sum_{l=0}^{L} \left[ \sum_{m=0}^{+l} f_{lm} \bar{P}_{lm}(\cos \theta) \cos(m\varphi) + \sum_{m=-l}^{-l} f_{lm} \bar{P}_{l|m|}(\cos \theta) \sin(|m|\varphi) \right] \\
&= \sum_{l=0}^{L} \left[ \sum_{m=0}^{+l} f_{lm} \bar{P}_{lm}(\cos \theta) \cos(m\varphi) + \sum_{m'=1}^{+l} f_{l,-m'} \bar{P}_{lm'}(\cos \theta) \sin(m'\varphi) \right] \tag{A4} \\
&= \sum_{l=0}^{L} \left[ \sum_{m=0}^{+l} C_{lm} \bar{P}_{lm}(\cos \theta) \cos(m\varphi) + \sum_{m=0}^{+l} S_{lm} \bar{P}_{lm}(\cos \theta) \sin(m\varphi) \right] \\
&= \sum_{m=0}^{L} \sum_{l=m}^{L} \left[ C_{lm} \bar{P}_{lm}(\cos \theta) \cos(m\varphi) + S_{lm} \bar{P}_{lm}(\cos \theta) \sin(m\varphi) \right].
\end{align*}
$$

The above computation is consistent if we set

$$
\begin{align}
C_{lm} &= f_{lm} \quad \text{if } m \ge 0, \nonumber \\
S_{lm} &= \begin{cases} 0 & \text{if } m = 0 \\ f_{l,-m} & \text{if } m \ge 1 \end{cases}. \tag{A5}
\end{align}
$$

B. Verification of the scaling relation for the
spherical harmonic coefficients of the
Newtonian potential

The Newtonian potential defined by Eq. (1) is special
in that its bare spherical harmonic coefficients Φ<sub>lm</sub>(r) can
be rescaled according to Eq. (13) which yields altitude-
independent coefficients. We denote these rescaled co-
efficients y<sub>lm</sub><sup>N</sup>. This peculiar property can be used as
an additional means of test ascertaining the quality of
our numerical approximations. Indeed, from our numer-
ical Φ(r, θ) maps of the Newtonian potential, we can
compute the rescaled coefficients y<sub>l0</sub><sup>N</sup> at several altitudes
and check whether or not they actually depend on the
altitude. Fig. 16 shows the result of this process for
~r ∈ {1.059, 1.111, 1.314} and l ∈ {1, ..., 10}. At first
sight, the scaling relation seems consistent with the nu-
merical data at low degree. It is however more difficult
to verify it at higher altitude and for higher degrees as
the rescaling process involves multiplying the bare coef-
ficients by ~r<sup>l+1</sup> which quickly blows up to infinity. The
bare coefficients being themselves plagued with numer-
ical errors — they are derived from spherical harmonic
decomposition algorithm on top of FEM computations —
we clearly do not expect this relation to perfectly hold in
this regime.

<figure><img src="image_17.png" /><figcaption>FIG. 16. Verification of the scaling relation between bare spherical harmonic coefficients Φ<sub>l0</sub>(~r) and dimensionless coefficients y<sub>l0</sub><sup>N</sup>(~r) obtained numerically. The rescaled coefficients should in principle be independent of the altitude at which they are computed.</figcaption></figure>

C. Additional checks on 2D numerical computations

In this Appendix, we present two additional checks
that were performed on all FEM computations of the
chameleon field done in this article. We mainly elabo-
rate on the ideas introduced in Sec. II D 2.
# 1. Check of the radial evolution of the chameleon field

The {Earth + mountain} system constitutes a small departure from spherical symmetry. Therefore, the behavior of the chameleon field in the outgoing radial direction should be close to that of the {Earth} system, which in turn is spherically symmetrical, and so purely radial. From a numerical viewpoint, such radial profiles are much easier to obtain than a less symmetrical case. Indeed in the former case, the Klein-Gordon equation (7) boils down to a simple ODE

$$ \alpha \frac{d}{d\tilde{r}} \left( \tilde{r}^2 \frac{d\tilde{\phi}}{d\tilde{r}} \right) = \tilde{r}^2 \tilde{\rho} - \tilde{r}^2 \tilde{\phi}^{-(n+1)}, \quad (C1) $$

where numerical resources (density of DOFs, order of the finite elements) can be increased without blowing up the time complexity of the algorithm. As a result, we can obtain benchmark solutions at relatively low cost, for all the cases discussed in this study treated as purely radial (i.e. without mountain, all other physical parameters being equal). We denote $\tilde{\phi}_{1D}(\tilde{r})$ such benchmarks, and $\tilde{\phi}_{2D}(\tilde{r}, \theta)$ the 2D field profiles presented throughout the article. We then implement the following metric:

$$ \text{for } \tilde{r} > 1 + h_m, \quad \Gamma_{\tilde{r}} = \frac{\min_{\theta \in [0, \pi]} \left| \tilde{\phi}_{1D}(\tilde{r}) - \tilde{\phi}_{2D}(\tilde{r}, \theta) \right|}{\left| \tilde{\phi}_{1D}(\tilde{r}) \right|} \quad (C2) $$

Note that this metric has the advantage of being relative, as opposed to the absolute criteria discussed in Sec. II D 2.

Tab. VII includes $\Gamma_{\tilde{r}}$ for radial coordinates $\tilde{r} \in \{1.059, 1.111, 1.314, 4.645, 6.617\}$ and for all $(\alpha, \text{atmospheric scenario})$ pairs considered in this work. Although there is no physical motivation for having $\Gamma_{\tilde{r}} \equiv 0$ systematically, the fact that it remains below one part in a thousand in the vast majority of cases reflects the good agreement between the radial benchmark and $\tilde{\phi}_{2D}$. By way of comparison, applying the same metric on the Newtonian potential yields $\Gamma_{\tilde{r}} \sim 10^{-7}$. Evaluating this metric at different altitudes is also a way to make sure that none of the 2D solutions behaves unexpectedly in the radial direction.

# 2. Check of the strong residual amplitude with respect to each term

As mentioned in Sec. II D 2, the strong residual alone does not provide much insight into how good the numerical approximation is at the end of Newton's iterations. However, it is meaningful to compare locally the size of the residual against the size of each term in it, namely

$$ \left| \alpha \tilde{\Delta} \tilde{\phi} \right|, \quad \left| \tilde{\rho} \right|, \quad \left| \tilde{\phi}^{-(n+1)} \right|. \quad (C3) $$

A numerical approximation deemed acceptable must be such that the residual should be at least a few orders of magnitude smaller than the dominant term in (C3).

This criterion was assessed for five specific values of $\tilde{r} \in \{1.059, 1.111, 1.314, 4.645, 6.617\}$ on all numerical approximations discussed in this study. In Fig. 17, we show several examples of scatter plots that allowed us to do this monitoring. Each sub-panel corresponds to a given altitude and a given $(\alpha, \text{atmospheric scenario})$ pair — both randomly chosen —, and depicts the absolute value of the residual (black dots) vs terms appearing in (C3) (pastel-colored squares) as functions of $\theta$. We see that the residual remains well-below the dominant term in absolute values.

# D. Spherical harmonic coefficients at different altitudes

For the sake of comprehensiveness, we provide histograms of the spherical harmonic coefficients of both the Newtonian potential $\Phi$ and the chameleon potential $\Psi$ (up to degree $l = 200$) at three altitudes in Fig. 18 (see Sec. III A 1 in the main text). The specific shapes of both potential decompositions hold at all three altitudes, although they get squashed toward lower degrees the higher we go. The oscillations that we observe at high degrees for $\tilde{r} \in \{1.111, 1.314\}$ are not deemed physical but can be rather attributed to numerical noise.

# E. Further insights into the fifth force experienced by an unscreened satellite

In Sec. IV A, we computed the total chameleon acceleration undergone by a satellite in orbit. We found that, as long as the satellite was not screened, the resulting force (computed numerically by integrating the gradient of the scalar field over the whole volume occupied by the satellite) was equal to that acting on an equal mass point-like particle. In other words, the back-reaction of the satellite on the scalar field, in the unscreened regime, is such that there is no self-force perturbing the Jordan frame geodesics. In this Appendix, we provide an explanation for this phenomenon observed through numerical simulations, based on several approximations that can be justified in the {Earth + Satellite} system.

# 1. The case of Newtonian gravity

Given two massive bodies labeled by the subscripts $i \in \{1, 2\}$, the total gravitational force acting on the second body is

$$ \mathbf{F}_2 = - \int_{V_2} \nabla \Phi_N(\mathbf{x}) \, d\mathbf{m}(\mathbf{x}). $$
<table><caption>TABLE VII. Metric $\Gamma_{\tilde{r}}$ for $\tilde{r} \in \{1.059, 1.111, 1.314, 4.645, 6.617\}$ and for all $(\alpha, \text{atmospheric scenario})$ pairs considered in this work.</caption><thead><tr><th></th><th>α</th><th>τ̃ = 1.059</th><th>τ̃ = 1.111</th><th>τ̃ = 1.314</th><th>τ̃ = 4.645</th><th>τ̃ = 6.617</th></tr></thead><tbody><tr><td rowspan="14">No-Atmosphere</td><td>10<sup>-4</sup></td><td>0</td><td>0</td><td>0</td><td>8.7 × 10<sup>-7</sup></td><td>8.0 × 10<sup>-7</sup></td></tr><tr><td>10<sup>-5</sup></td><td>0</td><td>0</td><td>0</td><td>7.9 × 10<sup>-7</sup></td><td>7.4 × 10<sup>-7</sup></td></tr><tr><td>10<sup>-6</sup></td><td>2.6 × 10<sup>-3</sup></td><td>1.4 × 10<sup>-3</sup></td><td>5.0 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-5</sup></td><td>2.7 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-10</sup></td><td>0</td><td>0</td><td>0</td><td>2.5 × 10<sup>-6</sup></td><td>1.7 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-11</sup></td><td>6.9 × 10<sup>-6</sup></td><td>1.7 × 10<sup>-6</sup></td><td>1.7 × 10<sup>-6</sup></td><td>1.7 × 10<sup>-6</sup></td><td>1.7 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-12</sup></td><td>6.9 × 10<sup>-6</sup></td><td>6.9 × 10<sup>-6</sup></td><td>6.9 × 10<sup>-6</sup></td><td>6.8 × 10<sup>-6</sup></td><td>6.9 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-14</sup></td><td>3.9 × 10<sup>-5</sup></td><td>3.9 × 10<sup>-5</sup></td><td>3.9 × 10<sup>-5</sup></td><td>3.8 × 10<sup>-5</sup></td><td>3.8 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-15</sup></td><td>1.5 × 10<sup>-5</sup></td><td>1.4 × 10<sup>-5</sup></td><td>1.4 × 10<sup>-5</sup></td><td>1.3 × 10<sup>-5</sup></td><td>1.3 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-16</sup></td><td>1.3 × 10<sup>-5</sup></td><td>7.5 × 10<sup>-6</sup></td><td>3.3 × 10<sup>-6</sup></td><td>9.7 × 10<sup>-7</sup></td><td>9.0 × 10<sup>-7</sup></td></tr><tr><td>10<sup>-18</sup></td><td>1.7 × 10<sup>-4</sup></td><td>1.9 × 10<sup>-4</sup></td><td>2.1 × 10<sup>-4</sup></td><td>2.0 × 10<sup>-4</sup></td><td>2.0 × 10<sup>-4</sup></td></tr><tr><td>10<sup>-20</sup></td><td>0</td><td>0</td><td>1.1 × 10<sup>-4</sup></td><td>9.0 × 10<sup>-6</sup></td><td>5.1 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-21</sup></td><td>0</td><td>0</td><td>5.1 × 10<sup>-5</sup></td><td>2.9 × 10<sup>-6</sup></td><td>1.3 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-23</sup></td><td>0</td><td>0</td><td>1.7 × 10<sup>-5</sup></td><td>5.7 × 10<sup>-9</sup></td><td>4.5 × 10<sup>-9</sup></td></tr><tr><td>10<sup>-24</sup></td><td>0</td><td>0</td><td>1.1 × 10<sup>-6</sup></td><td>4.8 × 10<sup>-10</sup></td><td>4.3 × 10<sup>-10</sup></td></tr><tr><td>10<sup>-25</sup></td><td>0</td><td>0</td><td>0</td><td>4.5 × 10<sup>-11</sup></td><td>4.9 × 10<sup>-11</sup></td></tr><tr><td>10<sup>-26</sup></td><td>0</td><td>0</td><td>1.1 × 10<sup>-12</sup></td><td>4.7 × 10<sup>-12</sup></td><td>4.7 × 10<sup>-12</sup></td></tr><tr><td>10<sup>-27</sup></td><td>0</td><td>0</td><td>0</td><td>4.0 × 10<sup>-13</sup></td><td>5.8 × 10<sup>-13</sup></td></tr><tr><td>10<sup>-28</sup></td><td>5.7 × 10<sup>-14</sup></td><td>0</td><td>0</td><td>2.8 × 10<sup>-14</sup></td><td>4.3 × 10<sup>-14</sup></td></tr><tr><td rowspan="10">Tenuous</td><td>10<sup>-6</sup></td><td>4.5 × 10<sup>-6</sup></td><td>4.3 × 10<sup>-6</sup></td><td>4.1 × 10<sup>-6</sup></td><td>4.0 × 10<sup>-6</sup></td><td>3.6 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-10</sup></td><td>4.7 × 10<sup>-6</sup></td><td>4.7 × 10<sup>-6</sup></td><td>4.7 × 10<sup>-6</sup></td><td>4.8 × 10<sup>-6</sup></td><td>4.4 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-11</sup></td><td>0</td><td>0</td><td>0</td><td>6.5 × 10<sup>-6</sup></td><td>6.2 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-12</sup></td><td>1.7 × 10<sup>-3</sup></td><td>9.4 × 10<sup>-4</sup></td><td>3.4 × 10<sup>-4</sup></td><td>4.0 × 10<sup>-5</sup></td><td>3.0 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-14</sup></td><td>1.0 × 10<sup>-3</sup></td><td>5.7 × 10<sup>-4</sup></td><td>2.4 × 10<sup>-4</sup></td><td>6.0 × 10<sup>-5</sup></td><td>5.4 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-15</sup></td><td>2.4 × 10<sup>-4</sup></td><td>1.5 × 10<sup>-4</sup></td><td>7.4 × 10<sup>-5</sup></td><td>2.6 × 10<sup>-5</sup></td><td>2.3 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-17</sup></td><td>1.9 × 10<sup>-5</sup></td><td>1.3 × 10<sup>-5</sup></td><td>8.4 × 10<sup>-6</sup></td><td>5.3 × 10<sup>-6</sup></td><td>4.7 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-20</sup></td><td>1.9 × 10<sup>-6</sup></td><td>2.0 × 10<sup>-6</sup></td><td>2.1 × 10<sup>-6</sup></td><td>2.5 × 10<sup>-6</sup></td><td>2.3 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-5</sup></td><td>7.5 × 10<sup>-6</sup></td><td>6.9 × 10<sup>-6</sup></td><td>5.7 × 10<sup>-6</sup></td><td>4.3 × 10<sup>-6</sup></td><td>3.9 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-6</sup></td><td>4.5 × 10<sup>-6</sup></td><td>4.3 × 10<sup>-6</sup></td><td>4.1 × 10<sup>-6</sup></td><td>4.0 × 10<sup>-6</sup></td><td>3.6 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-8</sup></td><td>2.6 × 10<sup>-3</sup></td><td>1.4 × 10<sup>-3</sup></td><td>5.0 × 10<sup>-4</sup></td><td>4.7 × 10<sup>-5</sup></td><td>3.1 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-10</sup></td><td>3.8 × 10<sup>-7</sup></td><td>2.8 × 10<sup>-6</sup></td><td>4.3 × 10<sup>-6</sup></td><td>4.7 × 10<sup>-6</sup></td><td>4.4 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-12</sup></td><td>1.1 × 10<sup>-2</sup></td><td>6.2 × 10<sup>-3</sup></td><td>2.2 × 10<sup>-3</sup></td><td>2.2 × 10<sup>-4</sup></td><td>1.5 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-14</sup></td><td>2.5 × 10<sup>-3</sup></td><td>1.4 × 10<sup>-3</sup></td><td>5.5 × 10<sup>-4</sup></td><td>8.8 × 10<sup>-5</sup></td><td>7.2 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-15</sup></td><td>4.9 × 10<sup>-5</sup></td><td>3.6 × 10<sup>-5</sup></td><td>2.5 × 10<sup>-5</sup></td><td>1.9 × 10<sup>-5</sup></td><td>1.9 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-17</sup></td><td>1.1 × 10<sup>-5</sup></td><td>8.3 × 10<sup>-6</sup></td><td>6.7 × 10<sup>-6</sup></td><td>5.7 × 10<sup>-6</sup></td><td>5.3 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-20</sup></td><td>0</td><td>9.7 × 10<sup>-7</sup></td><td>2.1 × 10<sup>-6</sup></td><td>2.4 × 10<sup>-6</sup></td><td>2.2 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-23</sup></td><td>0</td><td>7.7 × 10<sup>-5</sup></td><td>1.7 × 10<sup>-5</sup></td><td>7.3 × 10<sup>-4</sup></td><td>1.4 × 10<sup>-3</sup></td></tr><tr><td rowspan="10">Earth-like</td><td>10<sup>-6</sup></td><td>4.3 × 10<sup>-6</sup></td><td>4.3 × 10<sup>-6</sup></td><td>4.1 × 10<sup>-6</sup></td><td>4.0 × 10<sup>-6</sup></td><td>3.6 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-8</sup></td><td>1.9 × 10<sup>-5</sup></td><td>1.1 × 10<sup>-5</sup></td><td>5.9 × 10<sup>-6</sup></td><td>4.1 × 10<sup>-6</sup></td><td>3.7 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-10</sup></td><td>0</td><td>0</td><td>4.8 × 10<sup>-2</sup></td><td>3.2 × 10<sup>-3</sup></td><td>2.1 × 10<sup>-3</sup></td></tr><tr><td>10<sup>-15</sup></td><td>0</td><td>0</td><td>0</td><td>2.8 × 10<sup>-5</sup></td><td>2.7 × 10<sup>-5</sup></td></tr><tr><td>10<sup>-18</sup></td><td>0</td><td>0</td><td>0</td><td>3.3 × 10<sup>-6</sup></td><td>2.9 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-20</sup></td><td>0</td><td>0</td><td>0</td><td>1.7 × 10<sup>-6</sup></td><td>1.6 × 10<sup>-6</sup></td></tr><tr><td>10<sup>-21</sup></td><td>0</td><td>0</td><td>0</td><td>5.3 × 10<sup>-7</sup></td><td>6.3 × 10<sup>-7</sup></td></tr><tr><td>10<sup>-25</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>FIG. 17. Representation of the strong residual (black circles) and the various terms of the dimensionless Klein-Gordon equation (7) (pastel-colored squares) in absolute values as a function of $\theta \in [0, \pi]$. Each column corresponds to a given radial coordinate $\tilde{r} \in \{1.059, 1.314, 4.645\}$ whereas each row corresponds to a given pair $(\alpha, \text{atmospheric scenario})$. In all cases, the absolute value of the strong residual remains at least several orders of magnitude below the dominant term of the Klein-Gordon equation, which is in line with the criterion set out in Sec. II D 2. The splitting of the curves associated with each term is due to the fact that we use second-order finite elements.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 18. Spherical harmonic coefficients of the Newtonian potential (top row) and of the chameleon potential for $\alpha = 10^{-25}$ (bottom row). The spectra are computed at three different altitudes, namely $\tilde{r} \in \{1.059, 1.111, 1.314\}$.</figcaption></figure>
In the above expression, $V_2$ is the volume occupied by the body 2 and $\Phi_N$ is the total Newtonian potential created by the two bodies. Thanks to the linearity of the Poisson equation governing the Newtonian potential, one can apply the superposition principle $\Phi_N = \Phi_1 + \Phi_2$, where $\Phi_i$ is the potential sourced by the body $i$ alone. The $\nabla$ operator and the integral being linear, we get

$$ \mathbf{F}_2 = - \int_{V_2} \nabla \Phi_1(\mathbf{x}) \, dm(\mathbf{x}) - \int_{V_2} \nabla \Phi_2(\mathbf{x}) \, dm(\mathbf{x}). $$

Physically speaking, the first integral represents the force exerted by 1 on 2 while the second integral is the force exerted by 2 on 2, which must be zero according to Newton's third law. This can be mathematically proven fairly easily given that

$$ \Phi_2(\mathbf{x}) = -G \int_{V_2} \frac{dm(\mathbf{x}')}{\|\mathbf{x} - \mathbf{x}'\|} \\ \text{and } \nabla (\|\mathbf{x} - \mathbf{x}'\|^{-1}) = - \frac{\mathbf{x} - \mathbf{x}'}{\|\mathbf{x} - \mathbf{x}'\|^3}. $$

We thus get

$$ \begin{aligned} \int_{V_2} \nabla \Phi_2(\mathbf{x}) dm(\mathbf{x}) &= G \int_{V_2} \left( \int_{V_2} \frac{\mathbf{x} - \mathbf{x}'}{\|\mathbf{x} - \mathbf{x}'\|^3} dm(\mathbf{x}') \right) dm(\mathbf{x}) \\ &= \frac{G}{2} \int_{V_2} \int_{V_2} \frac{\mathbf{x} - \mathbf{x}'}{\|\mathbf{x} - \mathbf{x}'\|^3} dm(\mathbf{x}') dm(\mathbf{x}) \\ &\quad - \frac{G}{2} \int_{V_2} \int_{V_2} \frac{\mathbf{x}' - \mathbf{x}}{\|\mathbf{x}' - \mathbf{x}\|^3} dm(\mathbf{x}) dm(\mathbf{x}') \\ &= 0 \end{aligned} $$

In conclusion, despite disturbing the overall Newtonian potential, the body 2 experiences the force sourced by the body 1 only. Furthermore, if the body 2 is small enough that $\nabla\Phi_1$ is approximately constant over $V_2$, we recover the point-mass approximation, i.e. $\mathbf{F}_2 \simeq -m_2\nabla\Phi_1(\mathbf{x}_2)$.

## 2. The case of the chameleon field in the unscreened regime

The above demonstration relies mainly on the superposition principle, which is lost in the case of the chameleon field because of the nonlinear nature of the Klein-Gordon equation governing the scalar field Eq. (4). Nonetheless, let $\phi_{\text{tot}} := \phi_{\oplus} + \delta\phi$ be the chameleon field of the {Earth + Satellite} system, where $\phi_{\oplus}$ is the background field of the Earth alone. Working with the dimensionless version of the Klein-Gordon equation (7), we have by definition

$$ \left\{ \begin{aligned} \alpha \Delta \phi_{\text{tot}} &= (\rho_{\oplus} + \rho_{\text{Sat}} + \rho_{\text{vac}})(\mathbf{x}) - \phi_{\text{tot}}^{-(n+1)} \\ \alpha \Delta \phi_{\oplus} &= (\rho_{\oplus} + \rho_{\text{vac}})(\mathbf{x}) - \phi_{\oplus}^{-(n+1)} \end{aligned} \right. . $$

In the unscreened case, $\delta\phi$ can indeed represent a small perturbation with respect to the background field $\phi_{\oplus}$ —

see e.g. the case illustrated in Fig. 10. Then, the nonlinear term can be approximated as

$$ \phi_{\text{tot}}^{-(n+1)} \approx \phi_{\oplus}^{-(n+1)} - (n+1)\phi_{\oplus}^{-(n+2)}\delta\phi, $$

so that

$$ \alpha \Delta \delta \phi \approx \rho_{\text{Sat}}(\mathbf{x}) + (n+1) \phi_{\oplus}^{-(n+1)} \frac{\delta \phi}{\phi_{\oplus}}. \quad (\text{E1}) $$

The r.h.s. of Eq. (E1) can be further simplified if we assume that, at the satellite's altitude, $\phi_{\oplus}$ is close to its asymptotic value in vacuum, that is $\phi_{\oplus}(\mathbf{x}_{\text{Sat}}) \sim \rho^{-1/(n+1)}$. Then we have, depending on whether $\mathbf{x} \in V_{\text{Sat}}$,

- Inside the satellite: $\rho_{\text{Sat}}(\mathbf{x}) \neq 0$ and so

$$ \frac{\phi_{\oplus}^{-(n+2)} \delta\phi}{\rho_{\text{Sat}}(\mathbf{x})} \sim \frac{\delta\phi}{\phi_{\oplus}} \frac{\rho_{\text{vac}}}{\rho_{\text{Sat}}} \ll 1. $$

Consequently, Eq. (E1) can be legitimately approximated by a Poisson equation inside the satellite.

- Outside the satellite: $\rho_{\text{Sat}}(\mathbf{x}) = 0$. We still have $\delta\phi/\phi_{\oplus}^{n+2} \ll 1$ and $\delta\phi \to 0$ as one moves away from the satellite (while $\phi_{\oplus}^{-(n+2)}$ remains bounded) so that we essentially recover a Laplace equation.

In brief, we showed that, under some assumptions, $\delta\phi$ obeys a Poisson equation inside the satellite, and a Laplace equation outside the satellite. The Newtonian potential sourced by the satellite (denoted by $\Phi_2$ in the previous discussion) is governed by the same partial differential equation. Yet, same equations have the same solutions, which means that $\delta\phi$ has a role similar to $\Phi_2$. Therefore, following the demonstration made in the case of the Newtonian potential above, we get

$$ \begin{aligned} \mathbf{F}_{\text{Sat}}^{\text{th}} &= -\frac{\beta}{M_{\text{Pl}}} \int_{V_{\text{Sat}}} \nabla \phi_{\text{tot}} dm(\mathbf{x}) \\ &\approx -\frac{\beta}{M_{\text{Pl}}} \int_{V_{\text{Sat}}} \nabla \phi_{\oplus} dm(\mathbf{x}), \end{aligned} \quad (\text{E2}) $$

qed.

## F. Mathematical proof of the absence of symmetry in the orbital dynamics

*Context & Notations.* In Sec. IV D, we have seen on simulation results that, although the gravity field is exactly symmetric with respect to the line $\theta = 0$, the dynamics of a point-mass in orbit is not. Let us translate this statement into mathematical terms. Let $\{\mathbf{X}_*(t), t > 0\}$ be a trajectory in phase space that is a solution of the ODE of interest, i.e. $\forall t > 0$, $\dot{\mathbf{X}}_*(t) = F(t, \mathbf{X}_*(t))$. At some point, the particle will pass over the mountain so that we can define $t_m$, the time at which $\theta(t_m) = 0$ for the
<figure><img src="image_20.png" /><figcaption>FIG. 19. Visual support for the proof.</figcaption></figure>

first time. Demanding that the trajectory is symmetric with respect to $\theta = 0$ actually means

$$ \forall t \in [0, t_m], \quad \mathbf{X}_*(t) = \mathbf{X}_*(2t_m - t). \quad (\text{F1}) $$

The state vector $\mathbf{X}_*(t)$ has components

$$ \begin{aligned} \mathbf{X}_*(t) &= [\delta r(t), \dot{\delta} r(t), \delta \theta(t), \delta L(t)] \\ &= [x_1(t), x_2(t), x_3(t), x_4(t)]. \end{aligned} $$

We furthermore recall that the vector field $F: (s, \mathbf{Y}) \in \mathbb{R} \times \mathbb{R}^4 \mapsto (F_1, F_2, F_3, F_4) \in \mathbb{R}^4$ is given by

$$ \begin{aligned} F_1 &= y_2 \\ F_2 &= \frac{(L_0 + y_4)^2}{(a + y_1)^3} + g(a + y_1, \theta_0 + \omega s + y_3) \\ F_3 &= \frac{L_0 + y_4}{(a + y_1)^2} - \omega \\ F_4 &= h(a + y_1, \theta_0 + \omega s + y_3). \end{aligned} \quad (\text{F2}) $$

In the above, functions $g$ and $h$ refer to the gravitational potential partial derivatives $-\partial_r U$ and $-\partial_\theta U$ respectively.

**Theorem.** The perturbed Keplerian problem with vector field (F2) is not symmetric on both sides of the mountain in the general case.

**Proof by contradiction.** Let us suppose that (F1) holds and derive a (necessary) condition on function $g$ and $h$. Because $x_1$ is continuous, there exists $t_p \in [0, t_m]$ such that $x_1$ is monotonous over $[t_p, t_m]$. Letting $t_q = 2t_m - t_p$, $x_1$ is also monotonous over $[t_m, t_q]$ due to the symmetry (F1). We further set

$$ \begin{cases} A = x_1(t_p) = x_1(t_q) \\ B = x_1(t_m) \end{cases}. $$

Notations introduced so far are shown in Fig. 19. Let us assume for now that $A \neq B$ so that $x_1$ is actually strictly

monotonous over $[t_p, t_m]$ and $[t_m, t_q]$ respectively. Then $V := [\min(A, B), \max(A, B)]$ is not a degenerate interval and we can define

$$ \begin{array}{l} x_1^p: [t_p, t_m] \to V \\ t \mapsto x_1(t) \end{array} \qquad \begin{array}{l} x_1^q: [t_m, t_q] \to V \\ t \mapsto x_1(t) \end{array} $$

together with their respective inverse

$$ \begin{array}{l} z_1^p: V \to [t_p, t_m] \\ u \mapsto z_1^p(u) \end{array} \qquad \begin{array}{l} z_1^q: V \to [t_m, t_q]. \\ u \mapsto z_1^q(u) \end{array} $$

We will make use of the following property on the inverse functions

$$ \forall u \in V, \quad z_1^p(u) = 2t_m - z_1^q(u). \quad (\text{F3}) $$

Indeed, for $u \in V$, there exist two unique times $t_\alpha \in [t_p, t_m]$ and $t_\beta \in [t_m, t_q]$ such that $u = x_1^p(t_\alpha) = x_1^q(t_\beta)$. Reciprocally, $t_\alpha = z_1^p(u)$ and $t_\beta = z_1^q(u)$. Then

$$ \begin{aligned} z_1^p(u) &= z_1^p(x_1^q(t_\beta)) = z_1^p(x_1(t_\beta)) = z_1^q(x_1(2t_m - t_\beta)) \\ &= z_1^p(x_1^p(2t_m - t_\beta)) = 2t_m - t_\beta \\ &= 2t_m - z_1^q(u). \end{aligned} $$

We then compute the integral

$$ I := \int_{t_p}^{t_q} \frac{dx_1}{ds}(s) x_2(s) ds $$

by two different ways. On the one hand,

$$ I = \int_{t_p}^{t_q} \left| \frac{dx_1}{ds}(s) \right|^2 ds \quad (\text{F4}) $$

because $\dot{x}_1(s) = x_2(s)$ along the trajectory. On the other hand,

$$ I = \int_{t_p}^{t_m} \frac{dx_1}{ds}(s) x_2(s) ds + \int_{t_m}^{t_q} \frac{dx_1}{ds}(s) x_2(s) ds, $$

from which we can make the changes of variable $u = x_1^p(s)$ in the first integral and $u = x_1^q(s)$ in the second one, yielding

$$ \begin{aligned} I &= \int_{A}^{B} x_2(z_1^p(u)) du + \int_{B}^{A} x_2(z_1^q(u)) du \\ &= \int_{A}^{B} x_2(z_1^p(u)) du + \int_{B}^{A} x_2(2t_m - z_1^p(u)) du \\ &= 0 \quad \text{because of symmetry (F1).} \end{aligned} $$

From Eq. (F4), we immediately deduce that $\dot{x}_1 \equiv 0$ on $[t_p, t_q]$. The fact that $x_1$ is constant contradicts our previous assumption that $A \neq B$. Therefore, $A$ has to be equal to $B$. Put in perspective with the fact that $x_1$ is monotonous on $[t_p, t_m]$ and on $[t_m, t_q]$, we have
• $x_1$ is constant over $[t_p, t_q]$. Let $H$ be this constant and define the radial distance $R := a + H$.

• $\dot{x}_2 \equiv 0$ on $[t_p, t_q]$ as well.

The final stage of this demonstration follows from the specific form of the vector field $F$. Let $s \in [t_p, t_q]$. For convenience, we recall that $\theta = \theta_0 + \omega s + x_3(s)$ and we denote by $\partial_1, \partial_2$ the partial derivatives of a two-variable function with respect to the first and second variable respectively. Taking the derivative with respect to $s$ in the second equation of the ODE system $\dot{\mathbf{X}}_*(s) = F(s, \mathbf{X}_*(s))$ yields

$$ \frac{d}{ds}\{\dot{x}_2(s)\} = \frac{d}{ds}\left\{\frac{[L_0 + x_4(s)]^2}{R^3} + g(R, \theta)\right\} = 0 \\ \iff 2\dot{x}_4(s)\frac{L_0 + x_4(s)}{R^3} + [\omega + \dot{x}_3(s)]\partial_2 g(R, \theta) = 0. $$

In this last equation, we can substitute $\dot{x}_4$ and $\dot{x}_3$ by the r.h.s. of the ODE, yielding

$$ [L_0 + x_4(s)] [2h(R, \theta) + R\partial_2 g(R, \theta)] = 0. \quad (F5) $$

The total angular momentum $L(s) = L_0 + x_4(s) = R^2\dot{\theta}(s)$ cannot be zero, because otherwise the point-mass would be frozen it time. Moreover, the interval $J := \{\theta_0 + \omega s + x_3(s), s \in [t_p, t_q]\}$ is not a singleton because again, the point-mass cannot remain frozen in time. Therefore, Eq. (F5) leads straightforwardly to

$$ \forall \theta \in J, \quad 2h(R, \theta) + R\partial_2 g(R, \theta) = 0. \quad (F6) $$

Replacing functions $g$ and $h$ by their definition in relation to the gravitational potential $U$, we arrive at the final conclusion that $\forall \theta \in J$:

$$ \begin{array}{l} \partial_{\theta} [2U(R, \theta) + R\partial_{r}U(R, \theta)] = 0 \\ \text{i.e.} \quad 2U(R, \theta) + R\partial_{r}U(R, \theta) = C^{st} \end{array} \quad (F7) $$

Condition (F7) is very restrictive on the form of admissible gravitational potentials and one can check that the potential of the {sphere + mountain} system that we have been using throughout this article does not satisfy this criterion.

**Conclusion.** We found a necessary condition on the gravitational potential (F7) for the dynamics of a point-mass in orbit to be symmetric with respect to $\theta = 0$. As a remark, the potential created by a perfect sphere trivially satisfies the criterion.

## G. Orbital periods

In classical central force problems, the period $T$ of a satellite in circular orbit around a planet can be expressed as a function of the distance $r$ to the planet's center and the acceleration $a$ it undergoes: $T = 2\pi\sqrt{r/a}$. A direct consequence of this formula is that the addition of the chameleon acceleration to the Newtonian one will slightly modify the orbital period. In Sec. IV D 5, we laid emphasis on the fact that the use of different altitudes was one possible way of circumventing the issue of model uncertainties. We can thus examine the difference in orbital period for the two gravity models, namely

$$ \Delta T^{\text{New}} = 2\pi \left( \sqrt{\frac{R_{\oplus} + h_1}{a_{\text{New}}(r_1)}} - \sqrt{\frac{R_{\oplus} + h_2}{a_{\text{New}}(r_2)}} \right), $$

$$ \Delta T^{\text{cham}} = 2\pi \left( \sqrt{\frac{R_{\oplus} + h_1}{(a_{\text{New}} + a_{\text{cham}})(r_1)}} - \sqrt{\frac{R_{\oplus} + h_2}{(a_{\text{New}} + a_{\text{cham}})(r_2)}} \right), $$

with $r = R_{\oplus} + h$ and $a_{\text{New}}(r_1) = \mu_{\oplus}/r^2$.

Fig. 20 illustrates the difference $\mathcal{D} = \Delta T^{\text{New}} - \Delta T^{\text{cham}}$ (expressed in seconds) in the $(h_1, h_2)$-plane, for several values of the parameter $\alpha$. Equivalently, one can parameterize deviation from Newtonian gravity as $\Delta T^{\text{New}} = \Delta T^{\text{cham}}(1 + \epsilon)$. Then we have $\epsilon = \Delta T^{\text{New}}\mathcal{D}$. In order to remain consistent with the rest of this article, we have fixed $\Lambda = \Lambda_{\text{DE}}$ and $n = 1$. The setup being spherically symmetric, the numerical computation of the chameleon fifth force can be performed with 1D finite elements. We can see that the orbital period anomaly $\mathcal{D}$ can be greater than $10^{-5}$ s. This effect scales linearly with the number of completed orbits: after 100 000 orbits, one can expect the anomaly to be of the order of a second — which would take approximately 20 years for a satellite orbiting 1000 km above the Earth surface (ignoring all the perturbing forces otherwise present in a realistic scenario).

[1] P. Brax, S. Casas, H. Desmond, and B. Elder, Universe 8, 10.3390/universe8010011 (2022).

[2] C. M. Will, Living Reviews in Relativity 17, 4 (2014).

[3] P. Brax, C. van de Bruck, A.-C. Davis, J. Khoury, and A. Weltman, Phys. Rev. D 70, 123518 (2004).

[4] J. Khoury and A. Weltman, Phys. Rev. D 69, 044026 (2004).

[5] J. Khoury and A. Weltman, Phys. Rev. Lett. 93, 171104 (2004).

[6] E. Babichev, C. Deffayet, and R. Ziour, International Journal of Modern Physics D 18, 2147 (2009), https://doi.org/10.1142/S0218271809016107.
<figure><img src="image_21.png" /><figcaption>FIG. 20. $\mathcal{D} = \Delta T^{\text{New}} - \Delta T^{\text{cham}}$ [s]. We set $2 \times 10^4 \text{ km} \ge h_2 \ge h_1 \ge 5 \times 10^2 \text{ km}$ and the black solid line corresponds to $h_1 = h_2$. For $-\log_{10} \alpha \in \{12, 15, 21\}$, the anomaly $\mathcal{D}$ is maximal when the two considered altitudes are far apart. The case $\alpha = 10^{-22}$ exhibits the transition between the latter regime and a new regime ($\alpha \le 10^{-23}$) where $h_1$ has to remain low to produce a notable anomaly. This second regime arises because in the limit $\alpha \to 0$, the chameleon field quickly (i.e. at low altitude) gets close to its asymptotic value. The fifth force thus vanishes at higher altitudes, and so does the orbital period anomaly $\mathcal{D}$.</figcaption></figure>

[7] P. Brax, C. Burrage, and A.-C. Davis, *Journal of Cosmology and Astroparticle Physics* **2013** (01), 020.

[8] A. Vainshtein, *Physics Letters B* **39**, 393 (1972).

[9] A. Nicolis, R. Rattazzi, and E. Trincherini, *Phys. Rev. D* **79**, 064036 (2009).

[10] M. Jaffe, P. Haslinger, V. Xu, P. Hamilton, A. Upadhye, B. Elder, J. Khoury, and H. Müller, *Nature Physics* **13**, 938 (2017).

[11] P. Brax, C. van de Bruck, A.-C. Davis, D. F. Mota, and D. Shaw, *Phys. Rev. D* **76**, 124034 (2007).

[12] A. Upadhye, *Phys. Rev. D* **86**, 102003 (2012).

[13] P. Touboul, G. Métris, M. Rodrigues, J. Bergé, A. Robert, Q. Baghi, Y. André, J. Bedouet, D. Boulanger, S. Bremer, P. Carle, R. Chhun, B. Christophe, V. Cipolla, T. Damour, P. Danto, L. Demange, H. Dittus, O. Dhuicque, P. Fayet, B. Foulon, P.-Y. Guidotti, D. Hagedorn, E. Hardy, P.-A. Huynh, P. Kayser, S. Lala, C. Lämmerzahl, V. Lebat, F. m. c. Liorzou, M. List, F. Löffler, I. Panet, M. Pernot-Borràs, L. Perraud, S. Pires, B. Pouilloux, P. Prieur, A. Rebray, S. Reynaud, B. Rievers, H. Selig, L. Serron, T. Sumner, N. Tanguy, P. Torresi, and P. Visser (MICROSCOPE Collaboration), *Phys. Rev. Lett.* **129**, 121102 (2022).

[14] I. Ciufolini, A. Paolozzi, E. C. Pavlis, J. C. Ries, R. Matzner, C. Paris, E. Ortore, V. Gurzadyan, and R. Penrose, *The European Physical Journal C* **83**, 87 (2023).

[15] I. Ciufolini, A. Paolozzi, E. C. Pavlis, R. Koenig, J. Ries, V. Gurzadyan, R. Matzner, R. Penrose, G. Sindoni, C. Paris, H. Khachatryan, and S. Mirzoyan, *The European Physical Journal C* **76**, 120 (2016).

[16] B. Bertotti, L. Iess, and P. Tortora, *Nature* **425**, 374 (2003).

[17] J. G. Williams, S. G. Turyshev, and D. H. Boggs, *Phys. Rev. Lett.* **93**, 261101 (2004).

[18] J. G. Williams, S. G. Turyshev, and D. H. Boggs, *Classical and Quantum Gravity* **29**, 184004 (2012).

[19] G. Dvali, A. Gruzinov, and M. Zaldarriaga, *Phys. Rev. D* **68**, 024012 (2003).

[20] B. Jain, V. Vikram, and J. Sakstein, *The Astrophysical Journal* **779**, 39 (2013).

[21] V. Vikram, J. Sakstein, C. Davis, and A. Neil, *Phys. Rev. D* **97**, 104055 (2018).

[22] H. Wilcox, D. Bacon, R. C. Nichol, P. J. Rooney, A. Terukina, A. K. Romer, K. Koyama, G.-B. Zhao, R. Hood, R. G. Mann, M. Hilton, M. Manolopoulou, M. Sahlén, C. A. Collins, A. R. Liddle, J. A. Mayers, N. Mehrtens, C. J. Miller, J. P. Stott, and P. T. P. Viana, *Monthly Notices of the Royal Astronomical Society* **452**, 1171 (2015), https://academic.oup.com/mnras/article-pdf/452/2/1171/18754194/stv1366.pdf.

[23] C. Burrage and J. Sakstein, *Journal of Cosmology and Astroparticle Physics* **2016** (11), 045.

[24] J. Bergé, *Reports on Progress in Physics* **86**, 066901 (2023).

[25] M. Pernot-Borràs, J. Bergé, P. Brax, and J.-P. Uzan,
Phys. Rev. D **100**, 084006 (2019).

[26] M. Pernot-Borràs, J. Bergé, P. Brax, J.-P. Uzan, G. Métris, M. Rodrigues, and P. Touboul, Phys. Rev. D **103**, 064070 (2021).

[27] C. W. F. Everitt, D. B. DeBra, B. W. Parkinson, J. P. Turneaure, J. W. Conklin, M. I. Heifetz, G. M. Keiser, A. S. Silbergleit, T. Holmes, J. Kolodziejczak, M. Al-Meshari, J. C. Mester, B. Muhlfelder, V. G. Solomonik, K. Stahl, P. W. Worden, W. Bencze, S. Buchman, B. Clarke, A. Al-Jadaan, H. Al-Jibreen, J. Li, J. A. Lipa, J. M. Lockhart, B. Al-Suwaidan, M. Taber, and S. Wang, Phys. Rev. Lett. **106**, 221101 (2011).

[28] H. Moritz, Classical physical geodesy, in *Handbook of Geomathematics*, edited by W. Freeden, M. Z. Nashed, and T. Sonar (Springer Berlin Heidelberg, Berlin, Heidelberg, 2015) pp. 253–289.

[29] H. Lévy, J. Bergé, and J.-P. Uzan, Phys. Rev. D **106**, 124021 (2022).

[30] P. Yin, R. Li, C. Yin, X. Xu, X. Bian, H. Xie, C.-K. Duan, P. Huang, J.-h. He, and J. Du, Nature Physics **18**, 1181 (2022).

[31] K. Jones-Smith and F. Ferrer, Phys. Rev. Lett. **108**, 221101 (2012).

[32] R. Pourhasan, N. Afshordi, R. Mann, and A. Davis, Journal of Cosmology and Astroparticle Physics **2011** (12), 005.

[33] L. Ogden, K. Brown, H. Mathur, and K. Rovelli, Phys. Rev. D **96**, 124029 (2017).

[34] T. P. Waterhouse, An introduction to chameleon gravity (2006), arXiv:astro-ph/0611816 [astro-ph].

[35] T. Katsuragawa, T. Nakamura, T. Ikeda, and S. Capozziello, Phys. Rev. D **99**, 124050 (2019).

[36] J.-A. Gu and W.-T. Lin, Solar-system constraints on f(r) chameleon gravity (2011), arXiv:1108.1782 [gr-qc].

[37] D. F. Mota and D. J. Shaw, Phys. Rev. D **75**, 063501 (2007).

[38] L. Kraiselburd, S. J. Landau, M. Salgado, D. Sudarsky, and H. Vucetich, Phys. Rev. D **97**, 104044 (2018).

[39] T. Fukushima, Monthly Notices of the Royal Astronomical Society **462**, 2138 (2016).

[40] C. Briddon, C. Burrage, A. Moss, and A. Tamosiunas, Journal of Cosmology and Astroparticle Physics **2021** (12), 043.

[41] *U.S. Standard Atmosphere, 1976*, Tech. Rep. (U.S. Government Printing Office, Washington, D.C., 1976).

[42] N. K. Pavlis, S. A. Holmes, S. C. Kenyon, and J. K. Factor, Journal of Geophysical Research: Solid Earth **117** (2012).

[43] M. A. Wieczorek and M. Meschede, Geochemistry, Geophysics, Geosystems **19**, 2574 (2018).

[44] J. Bergé, P. Brax, M. Pernot-Borràs, and J.-P. Uzan, Classical and Quantum Gravity **35**, 234001 (2018).

[45] R. Cimrman, V. Lukeš, and E. Rohan, Advances in Computational Mathematics **45**, 1897 (2019).

[46] H.-S. Oh, B. Jang, and Y. Jou, Numerical Methods for Partial Differential Equations **19**, 301 (2003).

[47] Boulmezaoud, Tahar Zamène, ESAIM: M2AN **39**, 109 (2005).

[48] M. S. Nabizadeh, R. Ramamoorthi, and A. Chern, ACM Trans. Graph. **40**, 10.1145/3450626.3459809 (2021).

[49] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J.

Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods **17**, 261 (2020).

[50] H. Lévy, in *57th Rencontres de Moriond 2023 [Gravitation]* (2023).

[51] M. Pernot-Borràs, *Testing gravity in space : towards a realistic treatment of chameleon gravity in the MICROSCOPE mission*, phdthesis, Sorbonne Université (2020).

[52] H. P. Langtangen, *Computational Partial Differential Equations* (Springer Berlin Heidelberg, 2003).

[53] W. Frei, Nonlinearity ramping for improving convergence of nonlinear problems, COMSOL Blog (2013), accessed: May 26th, 2023.

[54] W. Frei, Load ramping of nonlinear problems, COMSOL Blog (2013), accessed: May 26th, 2023.

[55] E. L. Allgower and K. Georg, *Numerical Continuation Methods: An Introduction* (Springer Berlin Heidelberg, Berlin, Heidelberg, 1990).

[56] C. Burrage, E. J. Copeland, and J. A. Stevenson, Phys. Rev. D **91**, 065030 (2015).

[57] C. Geuzaine and J.-F. Remacle, *International Journal for Numerical Methods in Engineering* **79**, 1309 (2009).

[58] C. Burrage and J. Sakstein, *Living Reviews in Relativity* **21**, 1 (2018).

[59] E. Hairer, G. Wanner, and C. Lubich, Conservation of first integrals and methods on manifolds, in *Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations* (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 97–142.

[60] B. Tapley, B. Schutz, and G. Born, *Statistical Orbit Determination*, edited by B. D. Tapley, B. E. Schutz, and G. H. Born (Academic Press, Burlington, 2004).

[61] P. Schreiner, R. König, K. H. Neumayer, and A. Reinhold, *Advances in Space Research* **72**, 47 (2023).

[62] Y. Lou, X. Dai, X. Gong, C. Li, Y. Qing, Y. Liu, Y. Peng, and S. Gu, *Satellite Navigation* **3**, 15 (2022).

[63] D. Arnold, O. Montenbruck, S. Hackel, and K. Sośnica, Journal of Geodesy **93**, 2315 (2019).

[64] M. Pearlman, D. Arnold, M. Davis, F. Barlier, R. Biancale, V. Vasiliev, I. Ciufolini, A. Paolozzi, E. C. Pavlis, K. Sośnica, and M. Bloßfeld, *Journal of Geodesy* **93**, 2181 (2019).

[65] K. Abich, A. Abramovici, B. Amparan, A. Baatzsch, B. B. Okihiro, D. C. Barr, M. P. Bize, C. Bogan, C. Braxmaier, M. J. Burke, K. C. Clark, C. Dahl, K. Dahl, K. Danzmann, M. A. Davis, G. de Vine, J. A. Dickson, S. Dubovitsky, A. Eckardt, T. Ester, G. F. Barranco, R. Flatscher, F. Flechtner, W. M. Folkner, S. Francis, M. S. Gilbert, F. Gilles, M. Gohlke, N. Grossard, B. Guenther, P. Hager, J. Hauden, F. Heine, G. Heinzel, M. Herding, M. Hinz, J. Howell, M. Katsumura, M. Kaufer, W. Klipstein, A. Koch, M. Kruger, K. Larsen, A. Lebeda, A. Lebeda, T. Leikert, C. C. Liebe, J. Liu, L. Lobmeyer, C. Mahrdt, T. Mangoldt, K. McKenzie, M. Misfeldt, P. R. Morton, V. Müller, A. T. Murray, D. J. Nguyen, K. Nicklaus, R. Pierce, J. A. Ravich, G. Reavis, J. Reiche, J. Sanjuan, D. Schütze, C. Seiter, D. Shaddock, B. Sheard, M. Sileo, R. Spero, G. Spiers, G. Stede, M. Stephens, A. Sutton, J. Trinh, K. Voss, D. Wang, R. T. Wang, B. Ware, H. Wegener, S. Windisch, C. Woodruff, B. Zender, and M. Zimmer-
mann, Phys. Rev. Lett. 123, 031101 (2019).

[66] F. W. Landerer, F. M. Flechtner, H. Save, F. H. Webb, T. Bandikova, W. I. Bertiger, S. V. Bettadpur, S. H. Byun, C. Dahle, H. Dobslaw, E. Fahne-stock, N. Harvey, Z. Kang, G. L. H. Kruizinga, B. D. Loomis, C. McCullough, M. Murböck, P. Nagel, M. Paik, N. Pie, S. Poole, D. Strekalov, M. E. Tamisiea, F. Wang, M. M. Watkins, H.-Y. Wen, D. N. Wiese, and D.-N. Yuan, Geophysical Research Letters **47**, https://doi.org/10.1029/2020GL088306 (2020).

[67] E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, Rev. Mod. Phys. **93**, 025010 (2021).

[68] Q. Li, C. Xue, J.-P. Liu, J.-F. Wu, S.-Q. Yang, C.-G. Shao, L.-D. Quan, W.-H. Tan, L.-C. Tu, Q. Liu, H. Xu, L.-X. Liu, Q.-L. Wang, Z.-K. Hu, Z.-B. Zhou, P.-S. Luo, S.-C. Wu, V. Milyukov, and J. Luo, Nature **560**, 582 (2018).

[69] A. K. Chouhan, P. Choudhury, and S. K. Pal, Journal of Earth System Science **129**, 64 (2020).

[70] Y. Zhao, J. Zhang, J. Li, S. Liu, P. Miao, Y. Shi, and E. Zhao, Measurement Science and Technology **32**, 042002 (2021), publisher: IOP Publishing.

[71] G. Gabriel, D. Vogel, R. Scheibe, H. Lindner, R. Pucher, T. Wonik, and C. M. Krawczyk, Geophysical Journal International **184**, 1113 (2011), https://academic.oup.com/gji/article-pdf/184/3/1113/17362259/184-3-1113.pdf.

[72] A. M. Dziewonski and D. L. Anderson, Physics of the Earth and Planetary Interiors **25**, 297 (1981).

[73] R. Butler and S. Tsuboi, Physics of the Earth and Planetary Interiors **321**, 106802 (2021).

[74] R. S. Matoza and D. C. Roman, Bulletin of Volcanology **84**, 86 (2022).

[75] A. Kozak and A. Wojnar, Earthquakes as probing tools for gravity theories (2023), arXiv:2308.01784 [gr-qc].

[76] B. Christophe, D. Boulanger, B. Foulon, P.-A. Huynh, V. Lebat, F. Liorzou, and E. Perrot, Acta Astronautica **117**, 1 (2015).

[77] P. Touboul, B. Foulon, B. Christophe, and J. P. Marque, in *Geodesy for Planet Earth*, edited by S. Kenyon, M. C. Pacino, and U. Marti (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012) pp. 215–221.

[78] O. Montenbruck and E. Gill, *Satellite Orbits* (Springer Berlin Heidelberg, Berlin, Heidelberg, 2000).

[79] P. Musen, Journal of Geophysical Research (1896-1977) **65**, 1391 (1960).

[80] A. Schärer, R. Angélil, R. Bondarescu, P. Jetzer, and A. Lundgren, Phys. Rev. D **90**, 123005 (2014).

[81] E. Savalle, C. Guerlin, P. Delva, F. Meynadier, C. le Poncin-Lafitte, and P. Wolf, Classical and Quantum Gravity **36**, 10.1088/1361-6382/ab4f25 (2019).

[82] S. M. Brewer, J.-S. Chen, A. M. Hankin, E. R. Clements, C. W. Chou, D. J. Wineland, D. B. Hume, and D. R. Leibrandt, Phys. Rev. Lett. **123**, 033201 (2019).

[83] T. Bothwell, C. J. Kennedy, A. Aeppli, D. Kedar, J. M. Robinson, E. Oelker, A. Staron, and J. Ye, Nature **602**, 420 (2022).