# Modeling and testing screening mechanisms in the laboratory and in space

Valeri Vardanyan, Deaglan J Bartlett

► To cite this version:

Valeri Vardanyan, Deaglan J Bartlett. Modeling and testing screening mechanisms in the laboratory and in space. Universe, 2023, 9 (7), pp.340. 10.3390/universe9070340. hal-04127825

HAL Id: hal-04127825

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

Submitted on 10 Nov 2024

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
Review

# Modeling and Testing Screening Mechanisms in the Laboratory and in Space

Valeri Vardanyan<sup>1,2,*</sup> and Deaglan J. Bartlett<sup>3,*</sup>

1. Kavli Institute for the Physics and Mathematics of the Universe (WPI), UTIAS, The University of Tokyo, Chiba 277-8583, Japan

2. Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

3. Institut d'Astrophysique de Paris (IAP), CNRS & Sorbonne Université, UMR 7095, 98 bis bd Arago, F-75014 Paris, France

* Correspondence: valeri.vardanyan@ipmu.jp (V.V.); deaglan.bartlett@iap.fr (D.J.B.)

**Abstract:** The non-linear dynamics of scalar fields coupled to matter and gravity can lead to remarkable density-dependent screening effects. In this short review, we present the main classes of screening mechanisms, and discuss their tests in laboratory and astrophysical systems. We particularly focused on reviewing numerical and technical aspects involved in modeling the non-linear dynamics of screening and on tests using laboratory experiments and astrophysical systems, such as stars, galaxies, and dark matter halos.

**Keywords:** modified gravity; screening mechanisms; laboratory tests of gravity; astrophysical tests of gravity

## Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Summary of Theories</td><td>3</td></tr><tr><td>2.1. Thin-Shell Scenarios</td><td>3</td></tr><tr><td>2.2. Galileons</td><td>5</td></tr><tr><td>3. Numerical Methods and Technical Considerations</td><td>7</td></tr><tr><td>3.1. Relaxation Method</td><td>7</td></tr><tr><td>3.2. Finite Element Codes</td><td>10</td></tr><tr><td>3.3. Forces on Extended Objects</td><td>13</td></tr><tr><td>4. Laboratory Tests</td><td>14</td></tr><tr><td>4.1. Direct Force Measurements</td><td>15</td></tr><tr><td>4.2. Indirect Measurements</td><td>17</td></tr><tr><td>5. Astrophysical Scales</td><td>17</td></tr><tr><td>5.1. Stars</td><td>17</td></tr><tr><td>5.1.1. Hydrostatic Equilibrium—Stellar Structure and Evolution</td><td>17</td></tr><tr><td>5.1.2. Out of Equilibrium—Stellar Oscillations</td><td>19</td></tr><tr><td>5.2. Screening Maps</td><td>20</td></tr><tr><td>5.3. Galaxy Morphology</td><td>21</td></tr><tr><td>5.3.1. Thin-Shell Screened Theories</td><td>21</td></tr><tr><td>5.3.2. Vainshtein Screened Theories</td><td>22</td></tr><tr><td>5.4. Halo Properties</td><td>22</td></tr><tr><td>5.5. Splashback</td><td>23</td></tr></table>

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

check for
updates

**Citation:** Vardanyan, V.; Bartlett, D.J. Modeling and Testing Screening Mechanisms in the Laboratory and in Space. *Universe* **2023**, *9*, 340. https://doi.org/10.3390/universe9070340

Academic Editor: Lorenzo Iorio

Received: 5 April 2023

Revised: 30 May 2023

Accepted: 13 June 2023

Published: 20 July 2023

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

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
<table><tr><td>6. Conclusions</td><td>25</td></tr><tr><td>References</td><td>26</td></tr></table>

# 1. Introduction

The accelerated expansion of the Universe, discovered over two decades ago [1,2], still lacks a satisfactory theoretical explanation. Scalar fields have played a significant role in modeling cosmic acceleration and provide a relatively simple playground for establishing novel phenomenological features. Such scalar fields have been considered in two broad contexts, namely, dark energy (see e.g., Ref. [3]) and modifications of Einstein's general relativity (GR) (see e.g., Refs. [4–7]). The defining feature of modified gravity theories is the non-trivial coupling of scalar fields to gravity and matter. In modified gravity, therefore, one expects long-range “fifth” forces across a very broad range of scales, manifesting themselves from microscopically small scales to cosmologically large ones, characterized by the present-day expansion rate $H_0$.

In the cosmological regime, the phenomenology can be well-understood with linear perturbation theory methods. Particularly, the primary modifications are expressed as scale-and-time-dependent gravitational coupling and lensing potential (see e.g., Refs. [8,9]). Ongoing and upcoming large scale structure surveys are the primary testing ground for cosmological phenomenology of modified gravity. However, in order to have consistent dynamics on all the scales of interest, fifth forces should be screened away in regions much denser than the mean cosmological environments. Otherwise, such fifth forces would be ruled out on the grounds of local tests of gravity, which provide stringent constraints on additional interactions [10].

The inherently non-linear dynamics of non-trivially coupled scalar fields can lead to remarkably diverse effects depending on the environment in which they operate. Particularly, it has been realized that non-minimally coupled scalar fields with or without higher-derivative terms can lead to the so-called screening effects in dense environments. This environmental dependence allows the field to induce large deviations from GR on relatively poorly constrained cosmological structure formation and gravitational lensing scales, while having smaller effects in regimes where GR has been tested more precisely.

The most widely considered working scenario, known as *chameleon mechanism*, has been proposed and investigated in Refs. [11–13]. A related scenario is the so-called *symmetron screening*, introduced and studied in Refs. [14,15]. Both of these classes operate due to non-linearities in the potential sector. Non-linearities in the derivative sector can also lead to density-dependent screening, for example, the *Vainshtein mechanism* [16–18].

Due to their non-linear nature, the consistent modeling of screened fifth forces is a challenging task. Analytical progress is possible in certain situations, although such approaches are always problem-specific, and are hard to generalize. Numerical techniques have played a significant role in understanding screening phenomenology both in cosmological and astrophysical scenarios, and in laboratory experiments. In this review, we focused on summarizing the technical approaches, both analytical and numerical, to screen fifth force phenomenology. We refer the reader to already-existing excellent reviews covering theoretical and phenomenological aspects of screening mechanisms in greater detail. Among these references, Joyce et al. [6] provided a thorough review of modified gravity theories and their phenomenological aspects, Koyama [7] focused on cosmological phenomenology, Burrage and Sakstein [19,20] provided a detailed overview of chameleon and partly also symmetron tests, Sakstein [21] and Baker et al. [22] focused on astrophysical tests, and Brax et al. [23] focused on an effective field theory approach.

Instead of giving a detailed account of all the theoretical aspects of screening mechanisms, here, we focused on three representative classes. Namely, we will briefly introduce the chameleon, symmetron, and Vainshtein mechanisms in Section 2. In Section 3, we detail the state-of-the-art numerical methods used in studies of screened modified gravity. We particularly discuss the Finite-Difference and Finite-Element algorithms, as well as techniques for calculating the forces on extended objects. Section 4 is dedicated to laboratory tests. We
particularly discuss direct-force measurement experiments, such as Casimir and torsion balance tests, as well as indirect measurements, such as atomic interferometry and atomic spectroscopy experiments. In Section 5, we focused on a broad range of astrophysical tests. In particular, we reviewed the fifth force effects on stellar structure and evolution, galaxy morphology and dark matter halos structures. In this section, we also provide a review on a newly-emerging field of research which aims at producing screening maps of large-scale structures using constrained dark matter simulations. We conclude in Section 6.

## 2. Summary of Theories

### 2.1. Thin-Shell Scenarios

The simplest way to couple matter to scalar fields is via the universal conformal coupling. Chameleon [11–13] and symmetron [14,15] mechanisms, two of the most studied scenarios in the literature, rely on such conformal couplings. More concretely, we considered a generic theory described by the action,

$$S_{\varphi} = \int \sqrt{-g} \, d^4 x \left[ \frac{M_{\mathrm{Pl}}^2}{2} R - \frac{1}{2} \nabla^{\mu} \varphi \nabla_{\mu} \varphi - V(\varphi) \right] + S_M(\tilde{g}_{\mu\nu}, \Psi), \quad (1)$$

where $M_{\mathrm{Pl}}$ is the Planck mass, $R$ is the Ricci scalar corresponding to the Einstein-frame metric $g_{\mu\nu}$, and $S_M$ is the action for the matter fields. The real-valued scalar field $\varphi$ is assumed to be coupled to $g_{\mu\nu}$, while matter fields $\Psi$ couple to the Jordan frame metric $\tilde{g}_{\mu\nu}$. The two metrics are related by the following conformal transformation:

$$\tilde{g}_{\mu\nu} = A^2(\varphi)g_{\mu\nu}, \quad (2)$$

where $A(\varphi)$ is a generic function of the scalar field. The model is fully specified in terms of the coupling function $A(\varphi)$ and the potential $V(\varphi)$.

The scalar-field equation of motion is given by

$$\Box\varphi = V_{,\varphi} - \frac{\beta}{M_{\mathrm{Pl}}} A^4(\varphi) \tilde{T} \equiv V_{,\varphi}^{\mathrm{eff}}(\varphi), \quad \text{with } \beta(\varphi) \equiv M_{\mathrm{Pl}} \frac{d \log A}{d\varphi}, \quad (3)$$

where $\tilde{T}$ is the trace of the matter stress-energy tensor in the Jordan frame and we have introduced the effective potential $V^{\mathrm{eff}}(\varphi)$ for convenience. We note that the Einstein-frame stress-energy tensor is not conserved; however, it is often useful to introduce a matter density variable which is conserved in Einstein frame. Particularly, $\rho \equiv A^3\tilde{\rho}$ can be shown to be conserved with respect to the Einstein-frame metric. As a result, the effective potential takes the form

$$V^{\mathrm{eff}} = V(\varphi) + \rho A(\varphi). \quad (4)$$

Due to the non-trivial coupling, test particles feel an additional “fifth” force per unit mass given by

$$\vec{F}_{\varphi} = -\vec{\nabla} \log A(\varphi). \quad (5)$$

The qualitative realization of screening varies depending on the form of the coupling function $A(\varphi)$. In particular, the sum of runaway potentials and coupling functions in the chameleon mechanism leads to an environmentally-dependent effective mass for the scalar field. The larger mass in higher-density regions leads to shortened range of interactions. In contrast, the symmetron mechanism and its modifications rely on a suppression of the field expectation value in denser regions. The two mechanisms are schematically illustrated in Figure 1.

Prototypical examples of chameleons and symmetrons are given by

1. $V(\varphi) = \Lambda^4 \left(1 + \frac{\Lambda^n}{\varphi^n}\right)$, $A(\varphi) = e^{\beta\varphi/M_{\mathrm{Pl}}}$, with const $\Lambda, \beta$ and $n$ (chameleons),

2. $V(\varphi) = -\frac{1}{2}\mu^2\varphi^2 + \frac{1}{4}\lambda\varphi^4$, $A(\varphi) = 1 + \frac{1}{2}\frac{\varphi^2}{M^2}$, with const $\mu, \lambda$ and $M$ (symmetrons).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic illustration of chameleon (**left panel**) and symmetron (**right panel**) mechanisms. In high-density regions chameleons become more massive, while the symmetrons operate due to driving the field value to zero in such environments.</figcaption></figure>

The models described above and their variations have highly non-linear equations of motion. While we focused on methods for solving such equations in non-trivial scenarios, it is useful to recall that the most essential properties of solutions can already be demonstrated by a spherically symmetric, static top-hat dust configuration with radius $R$. The density profile is taken to be $\rho(r) = \rho\Theta(R - r)$, with $\Theta$ denoting a top-hat function. For concreteness, we will focus on symmetrons below. The exposition mainly follows Ref. [24]. In this one-dimensional configuration, the symmetron equation of motion is

$$ \varphi'' + \frac{2}{r}\varphi' = V_{,\varphi} + A_{,\varphi}\rho, \qquad (6) $$

where primes denote derivative with respect to the radial coordinate $r$. We are interested in finding the field profile satisfying the following boundary conditions

$$ \varphi'(r = 0) = 0, \quad \varphi(r \to \infty) = \varphi_0, \qquad (7) $$

where the vacuum expectation value in a $\rho = 0$ environment is given by $\varphi_0 = \mu / \sqrt{\lambda}$.

First, consider a sphere with density higher than the spontaneous symmetry breaking density $\rho_{SSB} \equiv \mu^2 M^2$. The field profile can be obtained analytically by approximating the system as a free scalar field with a potential centered at $\varphi = 0$. This leads to

$$ \varphi_{\text{in}}(r) = A \frac{R}{r} \sinh(\mu r \zeta), \qquad (8) $$

where $\zeta \equiv \sqrt{\frac{\rho}{\mu^2 M^2} - 1}$, and $A$ is an undetermined constant.

Outside the object, $r > R$, the field behaves as a free field with a quadratic potential centered at $\varphi = \varphi_0$, and the field profile is given by

$$ \varphi_{\text{out}}(r) = v + B \frac{R}{r} e^{-\sqrt{2}\mu r}, \qquad (9) $$

with $B$ being an undetermined constant. Matching the two solutions and their first derivatives at $r = R$ fixes the constants $A$ and $B$ as follows:

$$ B = -ve^{\sqrt{2}\mu R} \frac{\mu R \zeta \coth(\mu R \zeta) - 1}{\mu R \zeta \coth(\mu R \zeta) + \sqrt{2}\mu R}, \qquad (10) $$

$$ A = \frac{1}{\sinh \mu R \zeta} \left( v + B e^{-\sqrt{2}\mu R} \right). \qquad (11) $$
A limit of interest is the very strongly perturbing source with $\rho \gg \mu^2 M^2$. In this limit for the coefficient $B$ we have

$$B = -v \left( 1 - \frac{\tanh \sqrt{\alpha}}{\sqrt{\alpha}} \right), \qquad (12)$$

where the dimensionless parameter $\alpha \equiv \rho R^2/M^2$ controls the characteristic outer shell thickness within which the field profile is varied significantly, $\Delta R/R = 1/\alpha$.

It is instructive to consider large and small $\alpha$ limits. The coefficient $B$ takes the form

$$B(\alpha \ll 1) = -v \frac{\alpha}{3}, \quad B(\alpha \gg 1) = -v \left( 1 - \frac{1}{\sqrt{\alpha}} \right). \qquad (13)$$

The force between the sphere of mass $m_{\text{sphere}}$ and a test particle with mass $m_{\text{test}}$ is then given by

$$F \approx -\frac{v^2 m_{\text{test}} m_{\text{sphere}}}{4\pi M^4 r^2}, \quad \alpha \ll 1, \qquad (14)$$

$$F \approx -\frac{v^2 m_{\text{test}}(\lambda_{\text{sphere}} m_{\text{sphere}})}{4\pi M^4 r^2}, \quad \alpha \gg 1, \qquad (15)$$

where in the last result we have introduced the screening factor $\lambda_{\text{sphere}} \equiv 3/\alpha$. The essential assumption behind these results is that the test object does not alter the field profile of the sphere significantly. The important take-away is that, for large and dense objects, $\alpha \gg 1$, the fifth force on test particles is substantially suppressed due to the small screening factor $\lambda_{\text{sphere}} \ll 1$.

It is also instructive to compare the scalar force with the Newtonian gravitational interaction:

$$\frac{F_{\phi}}{F_{\text{N}}} = 2 \left( \frac{v M_{\text{Pl}}}{M^2} \right)^2, \quad \alpha \ll 1, \quad \frac{F_{\phi}}{F_{\text{N}}} = 2 \lambda_{\text{sphere}} \left( \frac{v M_{\text{Pl}}}{M^2} \right)^2, \quad \alpha \gg 1. \qquad (16)$$

The screening factor representation will be useful when discussing experimental configurations involving atoms and other test objects.

## 2.2. Galileons

Thus far, we have focused exclusively on *thin shell* screened modified gravity theories, where either the mass of the scalar field or the coupling of a test particle to the field is dependent on the environment, such that the fifth force is suppressed in regions with higher density or gravitational potential. We now consider the second main class of modified gravity theories: *kinetically screened theories*. This includes models exploiting, for example, the Vainshtein [16] or K-mouflage [25] mechanisms. Extensively studied models in this class are the Dvali-Gabadadze-Porrati (DGP) models [26]. In these cases, screening relies on higher order operators in the action, which may be beyond the validity of the relevant effective field theory. We do not dwell on this potential inconsistency in this section, and instead discuss the phenomenology which such models exhibit since future, consistent models may display similar behavior, and thus the results described here can be transferred to such models.

The canonical example of a Vainshtein-screened theory is the galileon [17]; scalar-tensor theories where the scalar, $\varphi$, obeys the galileon symmetry $\varphi \to \varphi + a + b_{\mu}x^{\mu}$, where $a$ and $b_{\mu}$ are constants. The cubic galileon has the action

$$S = \int d^4 x \sqrt{-g} \left[ \frac{R}{16\pi G} - \frac{1}{2} c_2 \nabla_{\mu} \varphi \nabla^{\mu} \varphi - \frac{c_3}{M^3} \Box \varphi \nabla_{\mu} \varphi \nabla^{\mu} \varphi - \mathcal{L}_{\text{m}} \right], \qquad (17)$$
for constants $c_2$ and $c_3$ and where $M^3 = M_{Pl}H_0^2$. The resulting equations of motion are second-order, and thus this is a special case of Horndeski theory [27]. In the mostly minus signature, the case $c_2 < 0$—dubbed the self-accelerating branch—has a non-canonical kinetic term and can thus self-accelerate, removing the necessity for a cosmological constant [28], whereas this cannot be achieved by the normal branch ($c_2 > 0$). Despite the negative $c_2$, the non-linear kinetic terms in Equation (17) admit ghost-free self-accelerating de Sitter solutions [17] and conditions on model parameters exist that prevent such instabilities arising in scalar and tensor perturbations [29]. Despite providing theoretically viable self-accelerating models, observations of the Integrated Sachs-Wolfe effect [30] make such a scenario unlikely for the cubic model, and measurements of the speed of gravitational waves with GW170817 [31,32] significantly constrain quartic and quintic galileons. Even if galileons cannot drive the Universe’s acceleration, the fifth-force due to galileons has interesting phenomenology, which we describe below.

In astrophysical and laboratory tests considered in this review, one can neglect terms suppressed by the Newtonian potential and its spatial derivatives and work in the quasi-static approximation to obtain the equation of motion [33],

$$ \nabla^2 \varphi + \frac{1}{3\beta_1 a^2 M^3} \left[ \left( \nabla^2 \varphi \right)^2 - \nabla_i \nabla_j \varphi \nabla^i \nabla^j \varphi \right] = \frac{M_{Pl}}{3\beta_2} 8\pi G_N a^2 \delta \rho_m, \quad (18) $$

where $i \in \{1, 2, 3\}$ and

$$ \beta_1 = \frac{1}{6c_3} \left[ -c_2 - \frac{4c_3}{M^3} (\ddot{\varphi} + 2H\dot{\varphi}) + \frac{16\pi G c_3^2}{M^6} \dot{\varphi}^4 \right], \quad \beta_2 = 2 \frac{M^3 M_{Pl}}{\dot{\varphi}^2} \beta_1. \quad (19) $$

Here, we are considering perturbations, $\varphi$, about a background, $\bar{\varphi}$. It is common to rewrite Equation (18) in terms of the cross-over scale, $r_c \equiv (\beta_1 a^2 M^3)^{-1}$ and coupling $\alpha \equiv M_{Pl}/(3\beta_2)$, instead of $\beta_1$ and $\beta_2$. Note that these parameters are function of time, but for low-redshift astrophysical tests the temporal evolution is typically neglected.

The acceleration experienced by a test object due to the fifth force is

$$ a_5 = -\alpha Q \nabla \varphi, \quad (20) $$

for

$$ Q = \int T^\mu_\mu d^3x, \quad (21) $$

where $T^{\mu\nu}$ is the stress-energy tensor of the test object [34]. For a non-relativistic object of mass $m$, $Q \approx m$ and thus the coupling to the field is the same as for gravity. In general, however, $Q < m$ since $T^\mu_\mu$ excludes the gravitational binding energy. The extreme case is for black holes, where $Q = 0$, which is a manifestation of the no hair theorem for galileons [35,36]. This is often phrased as a violation of the Strong Equivalence Principle (SEP), since the coupling is a function of an object’s composition. However, we note that the presence of an additional field means that this interpretation can be misleading; one would not say that the SEP is violated when an electrically charged and a neutral object take different trajectories in the presence of an electromagnetic field.

Despite the presence of non-standard kinetic terms in the equation of motion, for spherically symmetric systems, one can express the left hand side of Equation (18) as a total derivative, and thus derive the radial variation of $\varphi$ to be [37]

$$ \frac{d\varphi}{dr} = -\frac{GM\alpha}{r^2} \left( \left( \frac{r}{r_v} \right)^3 - \sqrt{\left( \frac{r}{r_v} \right)^6 + 4 \left( \frac{r}{r_v} \right)^3} \right), \quad (22) $$
where $M(r)$ is the mass (perturbation) enclosed at radius $r$, and the Vainshtein radius is defined to be

$$r_v = \frac{4}{3}GM(r)\alpha r_c^2, \qquad (23)$$

which is itself a function of $r$. In this case, we see that, well within the Vainshtein radius ($r \ll r_v$), the fifth force is suppressed by a factor $(r/r_v)^{3/2}$ relative to gravity, whereas far from a source the effect of the galileon is equivalent to increasing the gravitational constant by a factor $\Delta G/G_N = 2\alpha^2$. Similar behavior is observed for other variants of the galileon, but with different powers of $r/r_v$. A solar mass object has a Vainshtein radius of $r_v \sim O(100 \text{ pc})$ [38], allowing this model to evade solar system tests, as we are well within the regime where the fifth-force is suppressed.

Given the large Vainshtein radius of the Sun, one may suppose that the galileon field is screened in all regions of the Universe and thus any attempt to constrain such a model is futile, as there would not be any observational consequences with this level of suppression. There are two main reasons why this is incorrect.

First, we note that Equation (23) depends on the enclosed mass at a given radius, and is thus similar to Gauss' law. Therefore, within an extended mass distribution, the Vainshtein radius is smaller than for the corresponding point mass, allowing the galileon field to remain partially unscreened. This can lead to observational effects at distances $d \gtrsim 0.1 R_{200}$, where $R_{200}$ is the typical size of a dark matter halo [37]. Outside the mass distribution, the two fields are clearly equal.

Second, even if the field sourced by a given object is screened, it is important to remember that this object lives in an external environment. The galileon symmetry means that, if we have a solution for the galileon field $\varphi$, we can obtain another solution with $\partial_\mu \varphi \to \partial_\mu \varphi + b_\mu$, for constant $b_\mu$. This allows the galileon field sourced by large scale structure to remain unscreened, since this is approximately linear across the region of interest. This has been confirmed by cosmological simulations [39], which demonstrate that $\varphi$ has linear dynamics on $\gtrsim 10 \text{ Mpc}$ if $r_c \simeq 6 \text{ Gpc}$ [40–43].

### 3. Numerical Methods and Technical Considerations

#### 3.1. Relaxation Method

For many problems of interest, the density configurations can be assumed to be static and the scalar field equations of motion in e.g., Equation (3) or (18) can be considered as non-linear elliptical boundary value problems. Such an approximation is very well-motivated in the context of laboratory experiments as the experimental density configurations are kept static. The quasi-static assumption has also been widely considered in the astrophysical and cosmological contexts (see e.g., [44–46]) and has been verified in N-body simulations [47,48]. Qualitatively, the validity of this assumption can be understood by noticing that the typical timescale of the dynamics scales as $\sim \lambda_0/c$, with $\lambda_0$ denoting the vacuum Compton wavelength of the field, whereas the typical time-variation of matter in the cosmological context is given by $c/H_0$, which is several orders of magnitude larger than $\lambda_0/c$ in many relevant applications.

In order to introduce the relaxation algorithms which have been used in a number of analyses of screened modified gravity let us for concreteness focus on a 1-dimensional symmetron theory. The exposition below partly follows Ref. [49]. In numerical evaluation, it is useful to introduce dimensionless variables. We particularly introduce the normalized field variable $\chi \equiv \varphi/\varphi_0$ and consider the following equation:

$$\frac{d^2\chi}{dr^2} + \frac{2}{r}\frac{d\chi}{dr} = \frac{1}{2\lambda_0^2}\left[\left(\frac{\rho(r,t)}{\rho_{ssb}} - 1\right)\chi + \chi^3\right], \qquad (24)$$

where the vacuum Compton wavelength for symmetrons is given by

$$\lambda_0 = \frac{1}{\sqrt{2\mu}}. \qquad (25)$$
We impose the standard boundary conditions $\partial\chi/\partial r = 0$ at both $r = 0$ and $r \to \infty$.

We consider a discretization on a regular grid of size $h$ and employ a second-order discretization scheme. Higher order discretization schemes have been implemented and tested in, e.g., Ref. [50], without encountering significant performance differences in most of the problems of interest. The discretized equation can be written as

$$ \mathcal{L}[\chi_{i+1}, \chi_{i-1}; \chi_i] = 0, \qquad (26) $$

with

$$ \mathcal{L}[\chi_{i+1}, \chi_{i-1}; \chi_i] \equiv \mathcal{D}_K[\chi_{i+1}, \chi_{i-1}; \chi_i] - \mathcal{D}_P[\chi_i, \rho_i]. \qquad (27) $$

Here, we have introduced the following discretizations of the Laplace operator and the potential term,

$$ \mathcal{D}_K \equiv \frac{\chi_{i+1} + \chi_{i-1} - 2\chi_i}{h^2} + \frac{2}{r_i} \frac{\chi_{i+1} - \chi_{i-1}}{2h}, \quad \mathcal{D}_P = \frac{1}{\lambda_0^2} \left( \left( \frac{\rho_i}{\rho_{ssb}} - 1 \right) \chi_i + \chi_i^3 \right). \qquad (28) $$

In order to obtain a solution to the non-linear discretized problem, an iterative method should be employed. The particular structure of the operator $\mathcal{L}$ suggests the use of the Newton–Gauss–Seidel relaxation method. This has been a standard method used for obtaining the scalar field solutions in $N$-body simulations with additional fifth forces (see later).

The starting point of the iterative relaxation is the arbitrarily-chosen initial “guess” for $\chi$. At each iteration $n$ we use the current estimate $\chi_n(i)$ in order to obtain an improved estimate $\chi^{\text{new}}(i)$ according to

$$ \chi^{\text{new}}(i) = \chi_n(i) - \left. \frac{\mathcal{L}(\chi(i))}{\partial \mathcal{L}(\chi(i)) / \partial \chi(i)} \right|_{\chi(i)=\chi_n(i)}. \qquad (29) $$

In most of the practical problems, this iteration itself does not converge. Instead, one should typically use a fraction of $\chi^{\text{new}}(i)$ as a current estimate:

$$ \chi_{n+1}(i) = \omega\chi^{\text{new}} + (1 - \omega)\chi_n, \qquad (30) $$

where $0 < \omega \le 1$ is a constant hyper-parameter which should be tuned for each specific model separately. The optimal value of this parameter is problem-specific and is often chosen by empirically leveraging between being able to converge and the rate of change for $\chi_n$ between consecutive iterations.

Relaxation iterations are repeated multiple times until convergence is reached according to a pre-determined criterion. Two examples of convergence criteria are often employed in the literature. These are the global residual, defined as

$$ \mathcal{R}_1 \equiv \sqrt{\sum_i \mathcal{L}[\chi(i+1), \chi(i-1); \chi(i)]^2}, \qquad (31) $$

and the mesh-averaged change of the field profile between consecutive iterations

$$ \mathcal{R}_2 \equiv \sqrt{\sum_i (\chi^{\text{new}}(i) - \chi^{\text{old}}(i))^2}. \qquad (32) $$

One could, for example, choose $\mathcal{R}_2$ to reach a multiple of the numerical precision for the given implementation.

In order to validate the performance of the solver and assess the choice of the convergence criterion, one can use analytically solvable configurations. Let us consider two such examples. In order to construct the first analytical solution note that we can always plug a non-zero field profile in Equation (24) and reconstruct a unique density profile that serves
as a source for the mentioned field profile. As an example, we can choose $\chi \sim \tanh(r)$ and solve for the density configuration using Equation (24). The gray line in the left panel of Figure 2 shows the chosen field profile. The right panel of the same figure demonstrates the corresponding reconstruction of the density profile. The red dots in the left panel are the result of the numerical integration using the density profile from the right panel as an input source in our solver. As one can see, the numerical integration successfully matches the expected analytical field configuration.

<figure><img src="image_5.png" /><figcaption>Figure 2. Validation of the numerical solver. The gray line in the left panel represents an arbitrary non-zero function. The right panel depicts the reconstructed density profile obtained by inserting the function from the left panel into Equation (24). The red dots in the left panel show the numerical result for a problem with a density profile from the right panel. Figure adapted from [49].</figcaption></figure>

As our second example, we consider the configuration of two parallel plates with infinitely high density, separated by a vacuum gap [51]. Let the coordinate perpendicular to the plates be $z$ with the gap width being $\Delta z$ and the plate surfaces being placed at $-\Delta z/2$ and $+\Delta z/2$. The field equation of motion in this setup is given by

$$ \frac{d^2\chi}{d\hat{z}^2} = [(\tilde{\rho}(\hat{z}) - 1)\chi + \chi^3], \qquad (33) $$

where we have additionally defined $\hat{z} \equiv z/\sqrt{2}\lambda_0$ and $\tilde{\rho}(\hat{z}) \equiv \rho(\hat{z})/M^2\mu^2$ is assumed to be infinitely large inside the plates and zero in the gap. We can integrate this equation once in a $z$-interval where the density is constant. Choosing two subsequent intervals being $(0, \Delta\hat{z}/2)$ and $(\Delta\hat{z}/2, \infty)$, we can show that the value of the field on the plate surface $\chi_s$ is zero up to negligible corrections of the order of the ratio of the vacuum matter density to the plate density. Then, choosing an interval $(0, \hat{z})$ with $\hat{z} < \Delta\hat{z}/2$ we obtain

$$ \hat{z} = \frac{1}{\sqrt{1 - \frac{\chi_g^2}{2}}} \left[ \mathcal{F}\left(\pi/2, \sqrt{\frac{\chi_g^2}{2} - \chi_g^2}\right) - \mathcal{F}\left(\sin^{-1}\frac{\chi}{\chi_g}, \sqrt{\frac{\chi_g^2}{2} - \chi_g^2}\right) \right], \qquad (34) $$

where $\mathcal{F}$ is the elliptic integral of the first kind, and $\chi_g$ is the field value in the middle of the gap.

Fixing $\hat{z}$ to $\Delta\hat{z}/2$ and setting $\chi = 0$ as one of our boundary conditions, we can numerically solve for $\chi_g$. Having the latter we will then have $\chi$ as a function of $\hat{z}$ in the gap expressed in terms of the Jacobi elliptic function. The numerical relaxation solution of Equation (33) in the gap subject to boundary conditions $\chi(-\Delta\hat{z}/2) = 0 = \chi(+\Delta\hat{z}/2)$ is depicted in Figure 3 (red dots). This is in sub-percent-level agreement with the exact analytical solution obtained from Equation (34).
<figure><img src="image_6.png" /><figcaption>Figure 3. Analytical (gray line) and numerical (red) solutions of the symmetron field profile inside a gap between two parallel plates. Figure adapted from [49].</figcaption></figure>

Variations of the relaxation method have been implemented in N-Body simulation codes. The state-of-the-art codes are ECOSMOG [52], ISIS [53] and MG-GADGET [54]. See also Ref. [55] for a code-comparison study. In highly-demanding setups, such as N-body simulations, the scalar solvers can be further accelerated due to multigrid techniques. This relies on the observation that relaxation iterations efficiently reduce small-wavelength errors. Moving the problem to a coarser grid helps with efficiently reducing the longer-wavelength error modes in the original problem. Such multigrid techniques are implemented in the state-of-the-art modified gravity N-body codes. These codes additionally support adaptive mesh refinement.

## 3.2. Finite Element Codes

The finite element method (FEM) offers a very rich toolkit for solving partial differential equations. FEM-based methods are often employed in engineering tasks, while their use in astrophysics and cosmology has been limited. In this section we offer an introduction to the main FEM ideas and their applications to screening.

Let us consider the following general problem on a domain $\Omega$ with boundary $\partial\Omega$

$$ \nabla^2 u(x) = -f(x, u(x)), \quad x \in \Omega \qquad (35) $$

$$ u(x) = u_B, \quad x \in \partial\Omega, \qquad (36) $$

where $f(x, u)$ is a source which can generally depend non-linearly on both $x$ and $u$, and $u_B$ is the value of the solution at the boundary.

The first step in FEM is to reformulate the boundary-value problem as a variational problem. In order to achieve this, we multiply the equation by a well-defined “test” function $v(x)$:

$$ \int_{\Omega} v(x) \nabla^2 u(x) = - \int_{\Omega} fv, \qquad (37) $$

where $\int_{\Omega}$ is an integral over the domain $\Omega$. The test function $v(x)$ should satisfy certain properties, and we will particularly assume it vanishes on the entire boundary $\partial\Omega$. Next, by employing the Green’s theorem, the derivative order can be reduced, resulting in the following “weak-form” equation:

$$ \int_{\Omega} \nabla u(x) \nabla v(x) = \int_{\Omega} fv, \qquad (38) $$
which should be satisfied for all $v(x)$ in the relevant function space.

In practice, we seek a solution to a discretized problem, and instead of the functions $u(x)$ and $v(x)$ we consider the functions $u_h(x)$ and $v_h(x)$ belonging to a subspace of the original function space. We consider a basis expansion of both $u_h(x)$ and $v_h(x)$ as follows:

$$u_h(x) = \sum_i u_i e_i(x), \quad v_h(x) = \sum_i v_i e_i(x), \qquad (39)$$

where $u_i$ and $v_i$ are the values at vertices $\{P_i\}$ of the discretely triangulated grid of interest, and the basis functions can be chosen to satisfy $e_i(P_j) = \delta_{ij}$.

Inserting basis expansions into Equation (38), we arrive at

$$\sum_i \left( \int_{\Omega} \nabla e_i(x) \nabla v_j(x) \right) u_i = \int_{\Omega} f \left( x, \sum_i u_i e_i(x) \right) v_j. \qquad (40)$$

This system can be conveniently rewritten in a matrix form as

$$\mathbf{M}\mathbf{u} = \mathbf{b}, \qquad (41)$$

where

$$\mathbf{M}_{ij} = \int_{\Omega} \nabla e_i(x) \nabla v_j(x), \quad \mathbf{b}_j = \int_{\Omega} f v_j. \qquad (42)$$

Since $f$ is in general a non-linear function of unknown $u$ we need to rely on iterative non-linear solvers. Two general approaches are described below.

**Picard iteration:** The basic idea behind Picard iterations is to linearize the right-hand-side of Equation (41) around an initial guess $u_0$, and solve the resulting linear system of equations. Given an intermediate estimate $u_k$, the new estimate for $u$ is computed as follows. First, the non-linear function $f(x, u_{\text{new}})$ is replaced by its linear expansion around $u_k$:

$$f(x, u_{\text{new}}) \rightarrow f(x, u_k) - u_k \frac{\partial}{\partial u_k} f(x, u_k) + u_{\text{new}} \frac{\partial}{\partial u_k} f(x, u_k). \qquad (43)$$

Substituting this in Equation (41) and using Equation (39) again, we obtain

$$[\mathbf{M} + \mathbf{B}_k]\mathbf{u}_{\text{new}} = \mathbf{C}_k, \qquad (44)$$

where the matrix $\mathbf{B}_k$ and vector $\mathbf{C}_k$ are defined as

$$\mathbf{B}_k = - \int_{\Omega} \frac{\partial f(x, u_k)}{\partial u_k} e_i(x) v_j(x), \qquad (45)$$

$$\mathbf{C}_k = \int_{\Omega} \left( f(x, u_k) - u_k \frac{\partial}{\partial u_k} f(x, u_k) \right) v_i(x). \qquad (46)$$

Having solved for $\mathbf{u}_{\text{new}}$, the new estimate for the field profile is obtained as in Equation (30).

As a concrete example let us derive the corresponding linearized system for the power-law chameleon model described earlier. The equation of interest in dimensionless variables reads as

$$\alpha \nabla^2 \chi = -\chi^{-(n+1)} + \bar{\rho}, \qquad (47)$$

where $\alpha$ here is a dimensionless quantity determined in terms of the chameleon parameters, $\chi$ is a field variable normalized by the value minimizing the effective potential in the
ambient space, and $\tilde{\rho}$ is the matter density normalized by the density of the ambient space; see Ref. [56] for detailed definitions. The weak-form equation of motion is of the form

$$ \alpha \int_{\Omega} \nabla \chi(x) \nabla v_j(x) u_i = \int_{\Omega} \chi^{-(n+1)} v_j - \int_{\Omega} \tilde{\rho} v_j. \quad (48) $$

Following Equation (43), we expand the non-linear term around the k-th estimate of the solution, yielding

$$ \chi_{\text{new}}^{-(n+1)} \rightarrow (n+2)\chi_k^{-(n+1)} - (n+1)\chi_k^{-(n+2)}\chi_{\text{new}}. \quad (49) $$

In this case, the matrix $\mathbf{B}_k$ and vector $\mathbf{C}_k$ in Equation (44) take the form

$$ \mathbf{B}_k = \frac{n+1}{\alpha} \int_{\Omega} \chi_k^{-(n+2)} e_i(x) v_j(x), \quad (50) $$

$$ \mathbf{C}_k = \frac{n+2}{\alpha} \int_{\Omega} \chi_k^{-(n+1)} v_i(x) - \frac{1}{\alpha} \int_{\Omega} \tilde{\rho}(x) v_i(x). \quad (51) $$

Note that our notation here slightly differs from Ref. [56]. The resulting linear system can then be iterated over until convergence is achieved.

**Newton iteration:** In the Newton method, our goal is again to start from an initial guess and iteratively converge to the solution. Particularly, given the weak-form equation $\mathcal{L}(u, v_j) = 0$, we linearize the left-hand side around the current estimate for the solution, but now solve for the increment $\delta u \equiv u_{\text{new}} - u_k$. Equation (44) now changes into

$$ [\mathbf{M} + \mathbf{B}_k]\delta\mathbf{u} = -\mathbf{M}\mathbf{u}_k - \mathbf{D}_k, \quad (52) $$

where $\mathbf{M}$ and $\mathbf{B}$ are defined as in the Picard method, while $\mathbf{D}$ reads as

$$ \mathbf{D}_k = \int_{\Omega} f(x, u_k)v_i(x). \quad (53) $$

For the chameleon example discussed earlier, the vector $\mathbf{B}$ is given by Equation (50), and the quantity $\mathbf{D}$ coincides with $\mathbf{C}$ in Equation (51) but without the factor of $n+2$ in the first term. Once the solution for $\delta \mathbf{u}$ is found, the field estimate for the next iteration is obtained as

$$ u_{k+1} = u_k + \omega \delta \mathbf{u}. \quad (54) $$

Here, $0 < \omega \le 1$ is a hyperparameter analogous to the one in Equation (30), and controls the rate of iterative convergence.

In both Equations (44) and (52), the matrix $\mathbf{M}$ should only be computed once before the iterative evaluation starts. The matrix $\mathbf{M} + \mathbf{B}_k$, however, can in general be very large, especially for higher-dimensional problems. Inverting the linear equations with standard exact methods, therefore, is not always practical. In practice, additional iterative methods can be employed at each level of the Picard or Newton iterations in order to approximate the solutions to the linear systems of equations.

A class of powerful iterative methods for approximating the solutions of large linear systems of the form $\mathbf{A}\mathbf{x} = \mathbf{b}$ can be derived by treating the problem as a minimization problem for the function $f(x) = \frac{1}{2}\mathbf{x}^T\mathbf{A}\mathbf{x} - \mathbf{x}^T\mathbf{b}$. Such a problem can be efficiently solved with a conjugate gradient algorithm, which relies on starting with an initial guess $\mathbf{x}_0$, and updating the solution by minimizing the function $f$ along mutually conjugate directions $\{\mathbf{p}_i\}$ with $\mathbf{p}_i^T\mathbf{A}\mathbf{p}_j = 0$ if $i \ne j$. The latter construction ensures that minimizing solution along a direction $\mathbf{p}_k$ also minimizes the function along all the directions constructed in the previous iterations.

Currently, there are three publicly available, well-documented, and user-friendly FEM implementations for screening mechanisms in Python. The computational backbones of these are derived from the open-source FEM frameworks FEniCS [57,58] (implemented both
in C++ and Python) and SfePy [59] (implemented in Python). The φ-enics code [60], which is built on top of the FEniCS library, offers a 1-dimensional implementation of Vainshtein models. A more recent FEniCS-based code, SELCIE, has been presented in Ref. [56] and offers an implementation of chameleon theories in 2 and 3 spatial dimensions. The package has been employed for computing the chameleon field profiles in cosmic voids [61] and around halos [62].

The outer boundary conditions are set by the asymptotic behavior of the scalar field infinitely far away from matter sources. Since numerical implementations are limited to finite simulation boxes, the boundary conditions are typically imposed at distances sufficiently large compared to the Compton wavelength of the field. This approach is practical for many laboratory setups and isolated astrophysical configurations but can be computationally prohibitive if simulation domain is required to be very large. Adaptive mesh refinement techniques are often employed to reduce the computational cost in such situations. Alternatively, Ref. [63] presents a SfePy-based FEM package femtoscope, which offers a rigorous support for accurately implementing boundary conditions at infinity. The approach exploits a coordinate mapping from the infinite computational domain to a finite one, hence making the problem with an unbounded domain computationally tractable [64].

Earlier examples of using FEM solvers in the context of screening mechanisms can be found in Ref. [65] for chameleon theories and in Ref. [24] for symmetrons. Particularly, using FEniCS-based FEM solvers, Burrage et al. [65] simulated non-trivial density configurations in order to optimize for fifth force effects in laboratory experiments. Elder, Vardanyan, et al. [24] employed the commercially available FEM implementations in MATLAB alongside the custom-developed relaxation-based solvers in order to model symmetron fifth forces in Casimir experiments.

### 3.3. Forces on Extended Objects

Once the field profiles are computed, they can be used to calculate the fifth forces. So far we have only addressed the calculation of the fifth force on test particles which do not alter the ambient scalar field profile significantly. In realistic experiments, however, the relevant forces are often between extended, massive, and dense object. For example, in the Casimir experiments to be described later, the relevant force is either between two parallel plates or a plate and a large sphere. Calculation of the scalar force on extended objects is nuanced and we dedicate this subsection on introducing a method for systematically calculating the force. The discussion follows Ref. [34]; see also [24] for examples in the context of non-trivial laboratory configuration.

The starting point is the separation of the energy-momentum tensor $T_{\mu\nu}$ into the matter contribution $T_{\mu\nu}^{\text{matter}}$ and scalar field contribution $T_{\mu\nu}^{\phi}$. The quantity of interest is the 3-momentum in a spatial volume $V$ which is barely enclosing the surface of the extended object of interest and is therefore primarily dominated by the matter energy-momentum tensor:

$$P_i = \int_V d^3 x T_i^0. \qquad (55)$$

The force on the object can be computed from the time derivative of the 3-momentum. The latter can be written as:

$$\dot{P}_i = - \int_B d^2 \sigma_j T_i^j, \qquad (56)$$

where we have made use of the energy-momentum conservation, and the integral is over the surface $B$ of the volume $V$ with area element $d^2\sigma_i$.

As a specific example, let us consider a prototypical Casimir experiment with two infinitely-dense, parallel plates, separated by a vacuum gap. The left panel of Figure 4 depicts the configuration and the surface $B$ of a rectangular box. Our primary interest in
this example is the calculation of the Symmetron pressure on one of the plates, and we can make use of the analytical field profile derived using Equation (34).

<figure><img src="image_7.png" /><figcaption>Figure 4. Density configurations relevant to Casimir force experiments. The left configuration demonstrates the choice of the integration box used for estimating the scalar force on the plate. The right panel demonstrates the proximity-force calculations discussed in Section 4. Figure adapted from [24].</figcaption></figure>

Given the specific choice of the enclosing surface, the matter of the plate does not contribute to the integral in Equation (56). For the Symmetron field, we have

$$T_{\mu\nu}^{\varphi} = \partial_{\mu}\varphi\partial_{\nu}\varphi + \eta_{\mu\nu}\left(-\frac{1}{2}(\partial\varphi)^2 - V(\varphi)\right), \qquad (57)$$

From the symmetry of the problem, the only component of the force is in x-direction, and it reads as

$$\dot{P}_x = \frac{1}{2} \left[ \varphi'_{\text{in}} \left( \frac{L}{2} \right) \right]^2 - \frac{1}{2} \left[ \varphi'_{\text{out}} \left( \frac{L}{2} + W \right) \right]^2, \qquad (58)$$

where “in” and “out” denote the solutions inside and outside of the gap. The latter can be obtained by using the general form for the solution inside and assuming an infinitely wide gap. Primes denote derivatives with respect to the x-coordinate.

The fifth force per area A on the plate the follows:

$$\frac{F}{A} = \frac{\mu^4}{4\lambda} (\varphi_0^2(2 - \varphi_0^2) - 1), \qquad (59)$$

where a useful identity for Jacoby elliptic functions has been employed:

$$\left(\frac{d}{dx} \mathrm{cd}(x, m)\right)^2 = (1 - \mathrm{cd}(x, m)^2)(1 - m^2 \mathrm{cd}(x, m)^2). \qquad (60)$$

An alternative approach for calculating the pressure would have been to calculate the energy density stored in the scalar field configuration, and take the derivative with respect to the plate separation L. Direct computation shows agreement with the result presented in Equation (59).

## 4. Laboratory Tests

Laboratory experiments that are sensitive to the effects of weak forces and are capable of providing useful constraints on the screening mechanisms. We schematically divided such experiments into two broad categories. The experiments in the first category constrained the screening mechanisms due to the fifth force effects on test objects, such as atoms or neutrons. Such experiments are referred to as “indirect measurements” below, and include atomic interferometry [66–71], atomic spectroscopy [72,73], and bouncing
neutron [74] experiments. On the other hand, the experiments that try to directly measure the fifth forces between extended objects are referred to as “direct force measurements”, which include Casimir [75–79] and torsion balance experiments [80].

We made this distinction primarily because of the unique technical challenges as-
sociated with making robust and reliable theoretical predictions for those two classes of
experiments. While in the first class of experiments, the analytical approximations have
proven to be very useful and accurate enough, the direct force measurements often involve
very complex density geometries, and their analysis requires more careful numerical mod-
eling. Despite differences in modeling, it is important to combine the constraints derived
in different experiments since they often cover complementary regions of the parameter
space. Figure 5 represents a summary of constraints on symmetron models from laboratory
experiments. A similar summary for *n* = 1 chameleon models can be found in Ref. [81].

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Summary of laboratory constraints on symmetron parameters. Shaded regions and regions under the curves are ruled out by the respective experiment. Atomic interferometry, torsion balance and ultracold bouncing neutron results are reproduced from Refs. [51,70,74], respectively. The Casimir constraints should be interpreted as realistic forecast results. Figure adapted from [24].</figcaption></figure>

4.1. *Direct Force Measurements*

**Casimir Experiments:** The primary purpose of Casimir-force experiments is the mea-
surement of the force due to modified vacuum energy density in presence of conducting
extended bodies (see [82–84] for historical realizations of the experiment, and [76,85] for
modern and more precise experiments). Such experiments, however, can also be used
for constraining additional (classical) forces [75,76]. It is therefore natural to consider the
chameleon and symmetron fifth forces in the context of such experiments. Modern Casimir
experiments probe distance scales of ≈10 µm and are therefore sensitive to large effective
scalar-field masses.

The main theoretical complication lies in the difficulty of obtaining accurate estimates
of scalar field profiles and fifth forces in complex configurations of realistic experiments.
While the earlier Casimir force measurements were performed using two parallel plates
(see, e.g., Ref [77]), many modern setups can be approximated as a sphere placed in front of
the source plate [84–88] (see, however, Refs. [89,90] for modern experimental configurations
with planar geometry). The schematic setup of interest is depicted in the right panel of
Figure 4.

A careful treatment of systematic effects is essential to produce reliable constraints
on fifth-force models. For example, lateral forces and roughness are not captured by the
proximity force calculation [91,92] and the surfaces used have a finite conductivity [93],
both of which must be accounted for. Since the focus of this review is on the reconstruction
of the scalar fifth force profiles, we do not discuss these effects further here, but see e.g., Reynaud and Lambrecht [94] for a review.

Analytical progress can be made in the limiting cases of very small and large spheres. The exposition below focuses on symmetrons, and largely follows [24]. In the significantly simpler small-sphere limit, $R \ll \mu^{-1}$, we can treat the sphere as a test mass that does not significantly alter the ambient field profile sourced by the plate. The force on the sphere with mass $m_{sphere}$, in this case, can then be computed as

$$ \vec{F} = -\frac{\lambda_{sphere}m_{sphere}}{M^2}\varphi\vec{\nabla}\varphi, \qquad (61) $$

where $\lambda_{sphere}$ is the screening factor introduced in Section 2.1, and $\varphi$ is the analytically solvable field profile of the isolated plate (see the end of Section 3.1). As a result, we can obtain an accurate analytical force estimate when $R \ll \mu^{-1}$.

The situation is more subtle in the opposite limit of a very large sphere, $R \gg \mu^{-1}$. In this limit, the configuration can be approximated by a collection of parallel plates of varying gap openings by dividing the sphere into co-centric rings facing the plate. The exact (although not explicit) parallel-plate solutions can then be used to obtain the total force on the sphere. This approach is essentially the extension of the proximity force approximation widely employed for calculating the quantum Casimir forces (see e.g., [95,96]), and it has been used for studying Chameleon models in Casimir experiments [97].

Given the exact parallel-plate pressures $P$ for each of the co-centric rings, the total force can be expressed as

$$ F = -2\pi R^2 \int_{\pi/2}^{\pi} d\theta \sin\theta \cos\theta P(L), \qquad (62) $$

where $L = D + R + R \cos\theta$ is the x-distance between the plate and a ring (see Figure 4 for the notation).

We expect both the approximate approaches to fail in the intermediate regime in between of $R \gg \mu^{-1}$ and $R \ll \mu^{-1}$. For such intermediate radii, we should rely on exact numerical integration of the field profiles, which can be performed using any of the methods explained in Section 3 with sufficient care toward choosing an appropriately large box size so that the boundary effects do not alter the physical conclusions. Ref. [24] presents a detailed analysis of the plate-sphere configuration in symmetron models and finds that the proximity approximation provides an acceptable accuracy for large-enough spheres. The analysis also validates the results of the large and small sphere approximations. A similar analysis has been performed for chameleon models [81,97]. Interestingly, it has been demonstrated that, unlike the case of symmetrons, the proximity force approximation in the case of chameleons does not reproduce the numerical results even in the limit of large spheres. For chameleons, the approach based on the screening factor calculation for the sphere leads to more reliable analytical estimates [81]. This stresses once more the strong model-dependence of screening mechanisms.

**Torsion Balance:** One version of torsion-balance experiments [80] consists of two parallel cylinders. The bottom disc is turned uniformly while the upper one is a torsion pendulum and can rotate freely. The discs have regularly spaced holes with missing mass. In presence of a fifth force the torsion pendulum would experience a torque when its holes move across the holes of the lower rotating disc. The scalar field configuration is challenging to obtain in such a complex geometry. Particularly, the parallel-plane results can be used to obtain qualitatively robust results [98]. However, edge effects require a refined analytical modeling and numerical treatment.

Initial results have been obtained by numerically solving the field profile using the conjugate gradient minimization of the Hamiltonian of the system [99]. This direct numerical approach, however, can become prohibitive in certain parts of parameter space which require very fine grid sizes compared to the geometry of the problem. Refs. [51,100]
offered an approximate treatment of the system by introducing the so-called 1-dimensional plane-parallel approximation. The essential idea is to approximate the field on the pendulum cylinder surface by the corresponding surface field in a parallel-plate setup with the gap size corresponding to the distance of the point of interest and the nearest point on the bottom cylinder. Such an approach is similar to the proximity force approximation discussed earlier, and gives qualitatively correct results on the torque, although it quantitatively underestimates the numerical results.

## 4.2. Indirect Measurements

**Atomic Interferometry:** Atomic interferometry is a sensitive probe of gravitational free-fall acceleration. Such experiments are realized inside a vacuum chamber with a spherical source near its center. Atomic matter waves are beam-split and travel in free fall until their interference. The primary observable of interest is the accumulated phase difference between the split matter waves, which is linearly proportional to the acceleration of atomic clouds.

In the context of screening mechanisms, the acceleration is due to both gravity and the fifth force. The entire experiment is carried out inside a vacuum chamber with thick and dense walls. As a result, the field profile inside the chamber is independent of the ambient field profile outside. In order to induce fifth forces, a macroscopic source with e.g., spherical shape is placed near the center of the chamber. The total acceleration can be determined using the screening factors $\lambda_{atom}$ and $\lambda_{source}$ of the atom and the source and, for Chameleon theories, can be written as

$$a_{\varphi} + a_{\text{grav}} = \frac{G_{\text{N}}m_{\text{source}}}{r^2} \left( 1 + \lambda_{\text{atom}}\lambda_{\text{source}} \left( \frac{M_{\text{Pl}}}{M} \right)^2 \right). \qquad (63)$$

This expression includes both the gravitational acceleration $a_{grav}$ and the acceleration $a_{\varphi}$ due to fifth forces. A central feature of such configurations is that, while the source could be screened and have $\lambda_{source} \ll 1$, the atoms are typically not, $\lambda_{atom} = O(1)$, and the overall fifth-force acceleration is not doubly-suppressed.

Atomic interferometry experiments in the context of screened fifth force searches have been carried out by two independent groups [68,70,71]. Analytical modeling approximations have been extensively validated using exact numerical integration of the scalar field equations in the context of such experiments [101].

**Bouncing Neutrons and other methods:** Experiments including neutrons have been designed in order to measure the energy levels of quantum-mechanical particles in the gravitational field. In presence of an additional fifth force, the energy levels are perturbed and the parameters of the screening mechanisms can be strongly constrained. Such bouncing neutron experiments have been analyzed in the context of symmetrons in [74].

A logically similar method is to compute the shifts in atomic energy levels due to the fifth force on electrons [72,73]. The effect can be analytically described using the screening factor of the nucleus. In both the bouncing neutron and atomic spectroscopy methods the effect of the fifth force constitutes in perturbing the Hamiltonian of either the neutrons or electrons. Once the perturbed Hamiltonian is computed, the shifts in energy levels can be estimated by perturbation theory.

# 5. Astrophysical Scales

## 5.1. Stars

### 5.1.1. Hydrostatic Equilibrium—Stellar Structure and Evolution

In the presence of a fifth force, the structure and evolution of stars will be altered compared to GR, allowing stellar objects to be utilized as probes of such modifications. The gravitational collapse of a star is halted by its internal pressure such that, if gravity is stronger, then the required outward pressure must also increase, so stars are forced to burn their fuel at a faster rate. This shortens the lifetime of the star and will make them appear
brighter. The opposite would occur if the modification to gravity weakened the attractive force within the star.

To quantitatively study the impact of modifications to gravity, the only equation one needs to alter is the criterion for hydrostatic equilibrium [102,103]:

$$ \frac{1}{\rho} \frac{dP}{dr} = -\frac{GM}{r^2} - \frac{Y}{4} G \frac{d^2 M}{dr^2}, \qquad (64) $$

for pressure $P$, density $\rho$, and where $M$ is the mass enclosed at radius $r$. In standard gravity $Y = 0$, and is more generally given in terms of the $\alpha_B$, $\alpha_T$ and $\alpha_H$ parameters [104] as

$$ Y = \frac{4\alpha_H^2}{\alpha_H - \alpha_T - \alpha_B(1 + \alpha_T)}. \qquad (65) $$

The mere existence of stars already requires $Y > -2/3$, since below this value gravity becomes repulsive near the center of the star [103].

In the simplest models, one supposes that a star has a polytropic equation of state, $P = K\rho^{\frac{n+1}{n}}$, for positive constants $K$ and $n$. One can then define dimensionless variables in terms of the central density $\rho_c$ and pressure $P_c$ to be

$$ r = r_c \xi, \quad \rho = \rho_c \theta(\xi)^n, \quad P = P_c \theta(\xi)^{n+1}, \quad r_c^2 = \frac{(n+1)P_c}{4\pi G \rho_c^2}, \qquad (66) $$

and thus arrive at the modified Lame–Emden equation [102]

$$ \frac{1}{\xi^2} \frac{d}{d\xi} \left[ \left( 1 + \frac{n}{4} Y \xi^2 \theta^{n-1} \right) \xi^2 \frac{d\theta}{d\xi} + \frac{Y}{2} \xi^3 \theta^n \right] = -\theta^n, \qquad (67) $$

which has boundary conditions $\theta(0) = 1$, $\theta'(0) = 0$.

To make progress, one now needs to couple this equation with other equations describing stellar evolution. By using analytic approximations for the temperature profile and assuming a polytropic index $n$, Koyama and Sakstein [102] demonstrated that increasing $Y$ resulted in a reduction of the stellar luminosity at fixed mass, with stronger modifications for low mass stars since these are supported by gas pressure, rather than radiation pressure. This constrains $Y$ to be $\lesssim \mathcal{O}(1)$. Similarly, Sakstein [105,106] showed that the weakening of gravity within astrophysical bodies in scalar-tensor theories will increase the minimum mass a star must have to sustain hydrogen burning which, given the observed masses of Red Dwarf stars (which are greater than $0.08 M_{\odot}$ [107]), lead to similar constraints, ruling out $Y \gtrsim 1.6$. Modifications of gravity can also change the radius of brown dwarfs from their GR prediction ($0.1 R_{\odot}$) to arbitrarily large values if gravity is weaker, or down to $0.078 R_{\odot}$ if it is stronger.

To go beyond these semi-analytic models, one must solve all the relevant stellar structure equations simultaneously. This is commonly achieved by modifying the gravity model in the publicly available Modules for Experiments in Stellar Astrophysics (MESA) code [108–113].

After making such modifications, one can investigate the effect on different aspects of stellar evolution. For example, by considering the RGB phase, Chang and Hui [114] saw that, in Chameleon models, red giants can be significantly smaller (by tens of percent) and hotter (by hundred of Kelvins) than their GR counterparts. This effect is particularly pronounced in red giants since main sequence stars are much denser and thus have a suppressed scalar charge, whereas the outer envelope of red giants remains unscreened. If one considers $G^3$ galileons [104] exhibiting Vainstein screening, one finds that it is the main-sequence stars which are most strongly affected [102], and have a smaller effective temperature than GR stars at the same point in the Hertzsprung–Russell diagram. While such an effect could be degenerate with metallicity, modified gravity theories leave distinct evolutionary tracks in the Hertzsprung–Russell diagram, allowing these objects to constrain
such models. In the context of general modified gravity theories Saltas and Lopes [115] and Saltas and Christensen-Dalsgaard [116] have investigated the changes in the solar sound speed using analytical approximations and MESA simulations. These studies have concluded that precise helioseismic observations are able to improve the constraints on the fifth force coupling strength Y, and have derived an approximate constraint of $-10^{-3} < Y < 5 \times 10^{-4}$ at the $2\sigma$ confidence level.

### 5.1.2. Out of Equilibrium—Stellar Oscillations

Until now we have only considered stars in equilibrium. Moving beyond this, suppose one can consider small radius perturbations $\delta r$ in the star which obey the equation

$$ \frac{\partial^2 \delta r}{\partial t^2} = -\frac{1}{\rho} \frac{\partial P}{\partial r} - \frac{d\Phi}{dr}. \qquad (68) $$

In the absence of modified gravity, the period of small oscillations is found to be

$$ \Pi = \frac{2\pi}{\sqrt{4/3\pi G\rho(3\gamma - 4)}}. \qquad (69) $$

Hence, any enhancement in the strength of gravity ($G \to G + \Delta G$ with $\Delta G > 0$) would reduce the period [117]. The change in the period can be as large as 30% for $\mathcal{O}(1)$ couplings in chameleon screened theories [118]. This would be important for Cepheids—massive stars which pulsate with periods ranging from days to weeks when they pass through the instability strip—since the relation between their luminosity and period is used to determine their distance. This allows one to measure the strength of gravity in e.g., the Large Magellanic Cloud [119] by modeling the internal properties of Cepheids with MESA.

Moreover, since Cepheids are used in the construction of the distance ladder for local measurements of $H_0$, an incorrect distance estimation can bias the inferred value of $H_0$. Indeed, Desmond et al. [120] found that a fifth force with 5–30% the strength of gravity can reduce the Hubble tension. A similar effect can be obtained by instead calibrating to the TRGB [121]; if the hydrogen burning shell becomes unscreened, then the luminosity at the TRGB would be reduced [122] and thus distances based on this feature would be overestimated. Fifth-forces can therefore reduce the locally inferred value of $H_0$ thus reducing or alleviating the Hubble tension. Högas and Mörtsell [123] have studied the period-luminosity relation in the context of the symmetron theory. While this theory has been shown to worsen the Hubble tension, the distance ladder measurements impose strong constrains on the symmetron Compton wavelength.

Using the same notation as above, for Vainshtein screened theories, one can write $\delta r = r\chi(r) \exp(i\omega t)$ and linearise the pressure and potential to obtain a modified Linear Adiabatic Wave Equation [124]

$$ \frac{d}{dr} \left[ r^4 \left( \Gamma P_0 + Y \pi G r^2 \rho_0^2 \right) \frac{d\chi}{dr} \right] + r^3 \chi \frac{d}{dr} \left[ (3\Gamma - 4)P_0 \right] + r^4 \rho_0 \omega^2 W(r, Y) \chi = 0, \qquad (70) $$

where

$$ \Gamma = \left( \frac{\partial \log P}{\partial \log \rho} \right)_{\text{adiabatic}}, \quad W(r, Y) = 1 - \frac{Y \pi r^3 \rho_0 M}{2(M + \pi r^3 \rho_0)^2}, \qquad (71) $$

and a subscript '0' denotes an unperturbed quantity (see Sakstein [118] for the chameleon-screened version of this equation). For stars in GR, this results in unstable stars if $\Gamma < 4/3$, whereas in beyond Horndeski theory there is a second instability for sufficiently large Y, since this will result in $W < 0$ [124]. Using MESA, Sakstein et al. [124] showed that the modification to the period-luminosity relation for Cepheids also holds for these theories, and thus the distances inferred using this method can be compared to that from the TRGB to constrain Y.
## 5.2. Screening Maps

When determining whether a given object is screened or unscreened, one often computes properties of the local gravitational field: the potential, $\Phi$, for chameleons, the acceleration, $a$, for kinetically screened theories, or the local curvature (e.g., the Kretschmann scalar, $K$) in the presence of Vainshtein screening. These all have contributions from the object itself, but also from an object’s local environment. Particularly for astrophysical tests, it is therefore imperative to produce maps of these quantities to identify the screened and unscreened regions of the local Universe.

The Poisson equation for the potential, $\phi$, in the presence of a density perturbation, $\delta\rho$, is modified in $f(R)$ gravity to be

$$ \nabla^2 \phi = 4\pi G a^2 \delta\rho + \left( \frac{4\pi G a^2 \delta\rho}{3} + \frac{a^2}{6} \delta R \right) \equiv 4\pi G a^2 \delta\rho_{\text{eff}}, \quad (72) $$

where $\delta R$ is perturbation to Ricci scalar and $a$ is the scale factor. By integrating $a^2 \delta\rho_{\text{eff}}$ to obtain a dynamical mass, $M_D$, and comparing to true “lensing” mass, $M_L$, one can calculate a proxy for the degree of screening, $\Delta_m \equiv M_D/M_L - 1$ [125,126]. In this case, $\Delta_m \in [0, 1/3]$, where $\Delta_m = 0$ corresponds to a perfectly screened object, and $1/3$ is completely unscreened. Although this offers a useful metric for quantifying the degree of screening, we often do not have access to both the lensing and dynamical masses of objects, so must rely on alternative methods.

Cabré et al. [127] were the first to produce maps of the local gravitational potential for use in quantifying the degree of screening. For a set of galaxies of mass $M_i$ and virial radii $r_i$, they computed the external contribution to $\Phi$ as

$$ \Phi_{\text{ext}} = \sum_{d_i<r_i+\lambda_c} \frac{GM_i}{r_i}, \quad (73) $$

where $\lambda_c$ is the Compton wavelength of the field, and $d_i$ is the distance to the $i^{th}$ galaxy. One must introduce a cut-off based on the Compton wavelength, otherwise the above sum would diverge. This screening map covered $z < 0.05$ and the Sloan Digital Sky Survey (SDSS) footprint, since it was produced predominantly from a compilation of the SDSS group catalog [128], alongside other catalogs [129–132]. The approximations used were tested and calibrated against N-body simulations [133].

A more sophisticated analysis was performed by Desmond et al. [134], who extended these methods to produce full-sky maps up to ~200 Mpc, as well as maps of $a$ and $K$. Galaxies from the 2M++ galaxy catalog [132] were matched to ROCKSTAR [135] halos from the DARKSKY-400 simulation [136] using the abundance matching method of Lehmann et al. [137]. Assuming these halos have Navarro–Frenk–White (NFW) density profiles [138], $\Phi$ and $a$ were computed using all objects within $\lambda_c$ of a given point using standard formulae [139]. To compute the curvature, $K$ was approximated to be the sum of point mass contributions; this neglects the extended nature of halos and the non-linearity of this quantity, so should be treated as an order of magnitude approximation. To account for galaxies in halos which are too faint to be observed, correction factors for each quantity were applied by comparing to results from N-body simulations. To account for mass which does not reside in halos, a single constrained density field inferred using the Bayesian Origin Reconstruction from Galaxies (BORG) algorithm [140–144] applied to the 2M++ galaxy catalog was used.

Instead of computing screening proxies, one could directly compute the profiles of the dynamical fields which result in a modification to gravity, although one must do this separately for each theory one wishes to test. As part of the LOCal Universe Screening Test Suite (LOCUSTS) project, Shao et al. [145] solved for the scalar field for $n = 1$ Hu-Sawicki $f(R)$ gravity applying the ECOSMOG [52,146,147] code to a constrained N-body simulation
from the ELUCID project [148,149] (inferred using the SDSS DR7 galaxy catalog), using 20
$f_{R0}$ values in the range $10^{-7}$–$10^{-6}$.

Both the BORG and ELUCID density fields used to produce screening maps are con-
strained assuming no modification to gravity. Hence, these maps implicitly assume that
the low redshift Universe has matter fields which are similar for all reasonable modified
gravity and ΛCDM models on the scales of the reconstruction. The first constrained density
fields which include the effects of modified gravity have now been produced [150] from
CosmicFlows peculiar velocity data [151], through an extension of the ICECORE pack-
age [152]. By applying a Wiener filter [153] to obtain the linear density and velocity fields,
applying the reverse Zeldovich approximation [154], then adding fluctuations in poorly
constrained regions, they obtained estimates for the initial conditions at z = 49. They then
ran simulations using the MG-PICOLA [155] implementation of the COmoving Lagrangian
Acceleration (COLA) method [156] for nDGP and f(R) models, which solves for particle
trajectories in the frame of the reference given by Lagrangian perturbation theory. The
current implementation neglects the scale-dependence of the growth function in f(R) when
generating constrained realizations, hence they find that it is better to generate constrained
f(R) simulations from ΛCDM initial conditions than the ones they inferred. Future work
should be dedicated to producing fully-consistent constrained simulations and moving
away from the linear approximation required in this reconstruction [150]. Phenomeno-
logical approaches to incorporating screening effects in quasi-linear regimes of structure
formation could offer additional insights in understanding the large-scale environmental
effects on screening dynamics [157].

5.3. *Galaxy Morphology*

5.3.1. Thin-Shell Screened Theories

In thin-shell screened theories, where the degree of screening is determined by the
Newtonian potential, one can define a cut-off, Φ<sub>c</sub>, such that objects with potentials below
this value are unscreened and experience a fifth-force, and all other objects are screened.
If Φ<sub>c</sub>~10<sup>-6</sup>, then main sequence stars will always be screened, since this is the potential
at their surfaces, irrespective of their environment. However, if the local gravitational
potential is below this value, the surrounding gas and dark matter will be unscreened, and
hence feel an additional force. If there exists a gradient in the field responsible for the fifth
force, then the gas and dark matter will be accelerated whereas the stars will not, leading to
an offset between the center of mass of the stars and the remaining mass of the galaxy. This
offset will be stabilized by the induced gravitational potential, and this will also induce
a warp in the stellar disk. Thus, two morphological signatures are to be expected in such
a theory: (1) an offset parallel to the local gravitational field between the stars and gas of
galaxies in regions of low gravitational potential and (2) a warping of the stellar disk into a
characteristic 'U' shape [158].

Initial studies [159] of these phenomena, utilizing the screening maps of
Cabré et al. [127] (see Section 5.2) and combining data from the ALFALFA [160] and
SDSS [161] surveys, found no evidence for screened fifth forces, but the constraints were
not competitive with other astrophysical probes. More recently, galactic morphology has
played an increasingly important role in constraining such theories [162–164], with the
screening maps of Desmond et al. [134] being used to determine which objects are un-
screened and to predict the magnitude and direction of the offsets and warps. The most
recent constraints find that a fifth force with λ_c ≥ 0.3 Mpc must have a strength ΔG/G_N at
least as small as 0.8 at 1σ confidence [165]. n = 1 Hu-Sawicki f(R) with f_{R0} > 1.4 × 10⁻⁸ is
ruled out, resulting in practically all astrophysical objects being screened in the surviving
region of parameter space, and hence all astrophysically relevant models are disfavored.
The modeling of other astrophysical effects for these constraints has been shown to be
robust when calibrated against hydrodynamical ΛCDM simulations [166].
### 5.3.2. Vainshtein Screened Theories

One of the main astrophysical tests of the Vainshtein mechanism relies on the difference in coupling for non-relativistic matter ($Q = m$) and black holes ($Q = 0$) and was proposed by Hui and Nicolis [35]. This extreme difference means that, if the non-relativistic matter in a galaxy (stars, gas, dark matter) is falling down a gradient in the galileon field, then the supermassive black hole at its center lags behind as it does not experience the fifth force. This offset is stabilized by the induced gravitational potential, and thus one would expect to see an $O(kpc)$ offset between the black hole near a galaxy's center and its centroid, as measured by the stars.

This test was first performed by Sakstein et al. [38] by considering the black hole–galaxy offset in M87, utilizing the analysis of Asvathaman et al. [167]. Here, the galileon profile was determined analytically assuming spherical symmetry (Equation (22)), which was shown to be an excellent approximation for cosmological simulations (see e.g., Ref. [41]). As shown in Figure 6, $\alpha \sim O(1)$ was found to be excluded for cross-over scales $r_c \lesssim 1$ Gpc.

Bartlett et al. [168] extended these results by comparing the galaxy centroids to the X-ray and radio positions of active galactic nuclei (AGN) in a sample of 1916 galaxies compiled in [169]. They considered the galileon field sourced by large scale structure, utilizing a suite of constrained $N$-body simulations to produce maps of dark matter in the nearby Universe to predict both the magnitude and direction of the offset. By comparing to observations, the constraints on $\alpha$ were found to be tighter than those obtained from M87, with $\Delta G/G_N > 0.16$ ruled out at $1\sigma$ confidence (see Figure 6), although these constraints are applicable to larger $r_c$. We note that the constrained simulations were constructed assuming $\Lambda$CDM, and thus these constraints are only applicable for models which give rise to similar large scale structure as $\Lambda$CDM, and only affect objects on smaller scales. Furthermore, the Vainshtein mechanism was approximated as a hard cut-off, such that on large scales the galileon field traces the gravitational potential, but above some wavenumber the galileon field is set to zero. This is likely to be a less accurate approximation than that used by Sakstein et al. [38], although the advantage of using a statistical sample is that the non-galileon effects, which were ignored in Sakstein et al. [38], were empirically modelled and marginalized over.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Constraints on the galileon coupling parameter, $\alpha$, as a function of cross-over scale, $r_c$, from astrophysical tests. The shaded regions are excluded at $1\sigma$ confidence. For small $r_c$, the constraints are driven by Lunar Laser Ranging (LLR; [170,171]) and by considering the offset of the black hole in M87 [38]. For larger $r_c$, the constraints are derived from the offsets in a larger sample of galaxies, where the galileon field is sources by large scale structure [168]. Figure adapted from [168].</figcaption></figure>

### 5.4. Halo Properties

In GR, gravity determines the formation and growth of dark matter halos. Their abundance and internal structure can be modified in the presence of fifth forces, which might be partially unscreened especially for objects with smaller masses and on the outskirts
of halos. Since halos are strongly non-linear structures, N-body simulations have been the primary framework for deriving a quantitative understanding of these effects.

Modeling halo structure in modified gravity theories has been the focus of extensive past work, such as Refs. [172–179] for chameleon theories, Refs. [37,39,180–182] for Vainshtein mechanism, and Ref. [45] for symmetrons. Additionally, environmental dependences for both chameleon and Vainshtein mechanisms have been studied in Refs. [126,183].

Screening efficiency depends on the halo mass, and in the case of chameleon theories less massive halos tend to be unscreened. For Vainshtein, on the other hand, screening is less sensitive to the mass and practically all halos are screened [126,183]. As explained in Section 2.2, even though the inner regions of NFW halos in Vainshtein-screened theories are expected to be fully screened, one expects non-vanishing fifth force outside the virial radius. Simulations agree well with this theoretical prediction, with certain variations near the outskirts of the halos which is likely due to the breakdown of the NFW approximation [37,41,42].

As far as population-level mass distribution is concerned, it has been established that partially screened fifth forces induce significant modifications in the differential halo mass function $dn/dlogM_{halo}$, see e.g., [172,176,178]. While the exact results depend strongly on the model parameters, a general trend is that the distribution of very heavy halos is not affected. The main difference can be seen below a certain mass threshold, where halos are more abundant than in General Relativity. At even lower masses, halos are more abundant in GR, because such smaller halos are more efficiently absorbed by larger ones.

In ΛCDM, halos have universal density profiles described by the NFW parametrization:

$$ \rho(r)_{\text{dm}} = \frac{\rho_s}{\left(1 + \frac{r}{r_s}\right)^2}, \qquad (74) $$

where $\rho_s$ is a characteristic density, and $r_s$ is a characteristic scale for a halo. The compactness (concentration) of the profile is determined as the ratio of the virial radius $r_{200}$ and $r_s$; $c \equiv r_{200}/r_s$. The concentration and the characteristic density $\rho_s$ are connected via

$$ \rho_s = \frac{200}{3} \rho_c \frac{c^3}{g(c)}, \quad 4\pi r_s^3 \rho_s = \frac{M_{200}}{g(c)}, \qquad (75) $$

where $\rho_c$ is the critical density of the Universe, and $g(c) \equiv \log(1+c) - c/(1+c)$.

Earlier studies, such as Refs. [177,178], have examined the halo profiles in $f(R)$ N-body simulations and found that NFW is still an appropriate description, especially for larger halos. The profiles, however, can be more compact compared to GR. Detailed investigation of such an effect can be found in, e.g., Refs. [178,179], which found that halos with masses below a certain threshold are typically more concentrated. This, however, also depends on redshift, with the effect diminishing at relatively higher redshifts.

## 5.5. Splashback

In addition to modifying the halo concentrations and their mass function, fifth forces can also imprint non-trivial and potentially observable effects on the halo profiles. For larger, cluster-size halos this primarily concerns the outskirts of the density profiles, where fifth forces might be unscreened, hence significantly altering the dynamics of collapsing matter shells. In recent years, the *splashback* feature has gained importance as a physically motivated boundary of halos. This feature corresponds to the dividing boundary of single-streaming and multi-streaming regions and manifests itself as a steepening at outer regions [184]. The splashback location can, in principle, be inferred by measuring the radial distribution of subhalos, and assuming matter follows the same distribution [185–188]. Additionally, dedicated lensing measurements can also be used [189,190].

The position of the splashback feature is sensitive to any additional interactions and has been investigated in the context of $f(R)$ gravity and nDGP by Adhikari et al. [191].
Here, the authors used both semi-analytical and N-body approaches to predict the changes in splashback location. For symmetrons, splashback was investigated in Contigiani, Vardanyan, and Silvestri [50] by extending the self-similar collapse approach [192,193]. The latter approach, while not quantitatively very accurate, attains important insights into halo growth in presence of fifth forces. We will use this approach to demonstrate the changes in splashback location.

We will focus on the collapse in Einstein–de-Sitter (EdS) universe, and enforce the mass profile to satisfy a self-similar form:

$$M(r,t) = M(R,t)\mathcal{M}(r/R), \qquad (76)$$

where $R(t)$ is the turn-around radius at time $t$, and $M(r,t)$ is the mass within the radius $r$. The mass withing turn-around radius grows with scale factor as $M(R,t) \sim a(t)^s$, with accretion rate $s$, and in EdS is a power-law in cosmic time $t$.

In the absence of fifth force, the shell variables can be rescaled such that all the shells satisfy the same equation of motion,

$$\frac{d^2\lambda}{d\tau^2} = -\frac{\pi^2}{8} \frac{\tau^{2s/3}}{\lambda^2} \mathcal{M}\left(\frac{\lambda}{\Lambda(\tau)}\right), \qquad (77)$$

$$\mathcal{M}(y) = \frac{2s}{3} \int_{1}^{\infty} \frac{d\tau}{\tau^{1+2s/3}} \Theta\left(y - \frac{\lambda(\tau)}{\Lambda(\tau)}\right), \qquad (78)$$

where we have labeled each shell by its turn-around time $t_*$ and radius $r_*$, and have introduced $\lambda \equiv r/r_*$, $\tau \equiv t/t_*$, and $\Lambda = R(t)/r_*$. This system, complemented by initial conditions $\lambda(\tau)$ are $\lambda(\tau = 1) = 1$, $d\lambda/d\tau(\tau = 1) = 0$, can be iteratively solved until a converged mass-profile is obtained. The time-dependence of the density profile is determined in terms of the critical density $\rho_c(t)$ and $R(t)$. These equations can also be formulated for other 1-dimensional configurations [192], as well as tri-axial collapse [194]. The framework can also be extended to full $\Lambda$CDM backgrounds [195].

Once the self-similar halo profile is obtained as function of redshift, the symmetron field profiles can be computed with one of the methods introduced in Section 3. The extra force on outer shells can be added in Equation (77) in order to study the dynamics of the outer-most collapsing shells. It should be noted that, in the cosmological context presented here, the symmetron is fully screened at higher redshift. The unscreening redshift, $z_{ssb}$, is determined by the condition $\rho_c(z_{ssb}) = \rho_{ssb} = M^2\mu^2$. Figure 7 demonstrates the fifth-force profiles as a function of redshift, and the resulting effect on splashback location. From the right panel of the same figure it is clear that the innermost regions of the overdensity always stay screened. Substantial fifth forces start appearing in the outer regions when the $z \ll z_{ssb}$ condition is satisfied. At some point, a thin shell forms and the force gets concentrated in a narrow radial region. Due to this behavior, there exists an optimal value of $z_{ssb}$ which maximizes the change in splashback location due to fifth forces.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Splashback location in presence of symmetron fifth forces. The (**left**) panel depicts the phase-space distribution of shells in a halo, with colors denoting their turn-around redshifts. The dotted line depicts the distribution in the symmetron model. Splashback locations with and without symmetron force are marked by dotted and solid arrows, respectively. The (**right**) panel shows the symmetron force profile relative to the Newtonian force evaluated at different redshifts. Figure adapted from [50].</figcaption></figure>

**6. Conclusions**

The discovery of late-Universe cosmic acceleration motivated the extensive studies of alternatives to Einstein's theory of General Relativity. Such modified gravity theories generically involve additional scalar degrees of freedom. When such fields couple non-trivially to gravity and matter fields, they exert additional, "fifth" forces on matter particles. The presence of fifth forces leads to a remarkable variety of phenomenological effects, ranging from cosmologically large scales down to microscopic effects.

At large and linear scales, the novel effects have been relatively well-understood in terms of modified growth of structure—which in modified gravity can be scale- and time-dependent—and modified gravitational lensing. Significant progress has been made, first of all, by employing perturbation theory techniques in order to derive generic and largely model-independent observables, see e.g., [9,196]. Moreover, Effective Field Theory approaches have been very successful for systematically deriving unified prescriptions of linear perturbations [197–200]. Such unified approaches have been efficiently implemented in Einstein-Boltzmann numerical solvers [201–204].

On smaller scales, the dynamics of scalar fields exhibit highly non-linear behavior and lead to remarkable density-dependent screening effects. In this short review, we have summarized the main classes of screening mechanisms and have provided a broad summary of numerical and analytical methods developed in the context of non-linear dynamics. Unlike the linear regime, in screening regimes, scalar fields behave in a very model-dependent way. Even though, at the theoretical level, screening mechanisms can be categorized into distinct classes, such as chameleons, symmetrons, and Vainshtein mechanisms, and their precise phenomenological descriptions are given on a model-by-model basis. Even very similar mechanisms, such as chameleons and symmetrons, can often require different analysis approaches. An interesting and non-trivial example of this is the study of Casimir fifth forces in sphere-plane configurations, which closely resemble modern experimental setups. For example, the proximity-force approximation, explained in Section 4.1 can provide accurate description for the system for symmetrons, while fail for chameleons.

Numerical approaches are typically required in order to assess the validity of analytical approximations for a specific model. Providing an introductory summary of numerical algorithms is among the central topics of this review. While current numerical solvers have been designed with specific models in mind, most of the numerical approaches are applicable for a broad range of theories. This motivates the development of unified solver
packages. The finite element-based solvers introduced in Section 3.2, particularly, possess
very favorable dimensional scaling properties. A unified FEM-based code can be applied
not only to complex laboratory configurations, but also to astrophysical and cosmological
scenarios.

Identifying novel regimes where modified gravity can be tested is another promising avenue. Particularly, as outlined in Section 5.2, screening maps of the large-scale structure are useful for systematically exploring regions with discovery potential. To this end, constrained simulations are important in order to correctly account for environmental effects. In some of the existing studies, GR-based constrained simulations have been employed. Self-consistent modified gravity constrained simulations, such as the one presented in [150], will play an important role in producing precise screening maps.

Testing modified gravity theories in non-linear regimes will provide invaluable in-
sights about the nature of gravitational interactions. At the same time, the robust modeling
of non-linear phenomena is often very challenging and extensive cross-domain effort is
required. As the precision of both the laboratory and astrophysical tests is drastically
improving, approximate methods need to be replaced with more accurate approaches.

**Author Contributions:** All author contributed equally. All authors have read and agreed to the published version of the manuscript.

**Funding:** V.V. is supported by the WPI Research Center Initiative, MEXT, Japan and partly by JSPS KAKENHI grant 20K22348. D.J.B. acknowledges support by the Simons Collaboration on "Learning the Universe".

**Data Availability Statement:** Figures 1–3 are adapted from Vardanyan [49], Figures 4 and 5 are adapted from Elder, Vardanyan, et al. [24], Figure 6 is adapted from Bartlett, Desmond, and Ferreira [168]. Figure 7 is adapted from Contigiani, Vardanyan, and Silvestri [50].

**Conflicts of Interest:** The authors declare no conflict of interest.

References

1. Perlmutter, S. Measurements of Ω and Λ from 42 high redshift supernovae. *Astrophys. J.* **1999**, *517*, 565–586. [CrossRef]

2. Riess, A.G.; Filippenko, A.V.; Challis, P.; Clocchiatti, A.; Diercks, A.; Garnavich, P.M.; Gillil, R.L.; Hogan, C.J.; Jha, S.; Kirshner, R.P.; et al. Observational evidence from supernovae for an accelerating universe and a cosmological constant. *Astron. J.* **1998**, *116*, 1009–1038. [CrossRef]

3. Copeland, E.J.; Sami, M.; Tsujikawa, S. Dynamics of dark energy. *Int. J. Mod. Phys. D* **2006**, *15*, 1753–1936. [CrossRef]

4. Silvestri, A.; Trodden, M. Approaches to Understanding Cosmic Acceleration. *Rept. Prog. Phys.* **2009**, *72*, 096901. [CrossRef]

5. Clifton, T.; Ferreira, P.G.; Padilla, A.; Skordis, C. Modified Gravity and Cosmology. *Phys. Rept.* **2012**, *513*, 1–189. [CrossRef]

6. Joyce, A.; Jain, B.; Khoury, J.; Trodden, M. Beyond the Cosmological Standard Model. *Phys. Rept.* **2015**, *568*, 1–98. [CrossRef]

7. Koyama, K. Cosmological Tests of Modified Gravity. *Rept. Prog. Phys.* **2016**, *79*, 046902. [CrossRef]

8. Ade, P.A.; Aghanim, N.; Arnaud, M.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Bay, A.J.; Barreiro, R.B.; Bartolo, N.; Battaner, E.; et al. Planck 2015 results. XIV. Dark energy and modified gravity. *Astron. Astrophys.* **2016**, *594*, A14. [CrossRef]

9. Pogosian, L.; Silvestri, A. What can cosmology tell us about gravity? Constraining Horndeski gravity with Σ and µ. *Phys. Rev. 2016*, *D94*, 104014. [CrossRef]

10. Will, C.M. *Theory and Experiment in Gravitational Physics*; Cambridge University Press: Cambridge, UK, 1993.

11. Khoury, J.; Weltman, A. Chameleon cosmology. *Phys. Rev. D* **2004**, *69*, 044026. [CrossRef]

12. Khoury, J.; Weltman, A. Chameleon fields: Awaiting surprises for tests of gravity in space. *Phys. Rev. Lett.* **2004**, *93*, 171104. [CrossRef] [PubMed]

13. Brax, P.; van de Bruck, C.; Davis, A.C.; Khoury, J.; Weltman, A. Detecting dark energy in orbit: The cosmological chameleon. *Phys. Rev. D* **2004**, *70*, 123518. [CrossRef]

14. Hinterbichler, K.; Khoury, J. Symmetron Fields: Screening Long-Range Forces Through Local Symmetry Restoration. *Phys. Rev. Lett.* **2010**, *104*, 231301. [CrossRef] [PubMed]

15. Hinterbichler, K.; Khoury, J.; Levy, A.; Matas, A. Symmetron Cosmology. *Phys. Rev. D* **2011**, *84*, 103521. [CrossRef]

16. Vainshtein, A.I. To the problem of nonvanishing gravitation mass. *Phys. Lett. B* **1972**, *39*, 393–394. [CrossRef]

17. Nicolis, A.; Rattazzi, R.; Trincherini, E. Galileon as a local modification of gravity. *Phys. Rev. D* **2009**, *79*, 064036. [CrossRef]

18. Babichev, E.; Deffayet, C. An introduction to the Vainshtein mechanism. *Class. Quant. Grav.* **2013**, *30*, 184001. [CrossRef]

19. Burrage, C.; Sakstein, J. A Compendium of Chameleon Constraints. *J. Cosmol. Astropart. Phys.* **2016**, *11*, 045. [CrossRef]

20. Burrage, C.; Sakstein, J. Tests of Chameleon Gravity. *Living Rev. Relativ.* **2018**, *21*, 1. [CrossRef]

21. Sakstein, J. Astrophysical tests of screened modified gravity. *Int. J. Mod. Phys. D* **2018**, *27*, 1848008. [CrossRef]
22. Baker, T.; Barreira, A.; Desmond, H.; Ferreira, P.; Jain, B.; Koyama, K.; Li, B.; Lombriser, L.; Nicola, A.; Sakstein, J.; et al. Novel Probes Project: Tests of gravity on astrophysical scales. *Rev. Mod. Phys.* **2021**, *93*, 015003. [CrossRef]

23. Brax, P.; Casas, S.; Desmond, H.; Elder, B. Testing Screened Modified Gravity. *Universe* **2021**, *8*, 11. [CrossRef]

24. Elder, B.; Vardanyan, V.; Akrami, Y.; Brax, P.; Davis, A.C.; Decca, R.S. Classical symmetron force in Casimir experiments. *Phys. Rev. D* **2020**, *101*, 064065. [CrossRef]

25. Babichev, E.; Deffayet, C.; Ziour, R. k-MOUFLAGE Gravity. *Int. J. Mod. Phys. D* **2009**, *18*, 2147–2154. [CrossRef]

26. Dvali, G.; Gabadadze, G.; Porrati, M. 4D gravity on a brane in 5D Minkowski space. *Phys. Lett. B* **2000**, *485*, 208–214. [CrossRef]

27. Horndeski, G.W. Second-Order Scalar-Tensor Field Equations in a Four-Dimensional Space. *Int. J. Theor. Phys.* **1974**, *10*, 363–384. [CrossRef]

28. Babichev, E.; Esposito-Farèse, G. Time-dependent spherically symmetric covariant Galileons. *Phys. Rev. D* **2013**, *87*, 044032. [CrossRef]

29. de Felice, A.; Tsujikawa, S. Cosmology of a Covariant Galileon Field. *Phys. Rev. Lett.* **2010**, *105*, 111301. [CrossRef]

30. Renk, J.; Zumalacárregui, M.; Montanari, F.; Barreira, A. Galileon gravity in light of ISW, CMB, BAO and $H_0$ data. *J. Cosmol. Astropart. Phys.* **2017**, *2017*, 020. [CrossRef]

31. LIGO Scientific Collaboration; Virgo Collaboration. GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral. *Phys. Rev. Lett.* **2017**, *119*, 161101. [CrossRef]

32. Ezquiaga, J.M.; Zumalacárregui, M. Dark Energy After GW170817: Dead Ends and the Road Ahead. *Phys. Rev. Lett.* **2017**, *119*, 251304. [CrossRef] [PubMed]

33. Barreira, A.; Li, B.; Hellwing, W.A.; Baugh, C.M.; Pascoli, S. Nonlinear structure formation in the cubic Galileon gravity model. *J. Cosmol. Astropart. Phys.* **2013**, *2013*, 027. [CrossRef]

34. Hui, L.; Nicolis, A.; Stubbs, C.W. Equivalence principle implications of modified gravity models. *Phys. Rev. D* **2009**, *80*, 104002. [CrossRef]

35. Hui, L.; Nicolis, A. Proposal for an Observational Test of the Vainshtein Mechanism. *Phys. Rev. Lett.* **2012**, *109*, 051304. [CrossRef] [PubMed]

36. Hui, L.; Nicolis, A. No-Hair Theorem for the Galileon. *Phys. Rev. Lett.* **2013**, *110*, 241104. [CrossRef]

37. Schmidt, F.; Hu, W.; Lima, M. Spherical collapse and the halo model in braneworld gravity. *Phys. Rev. D* **2010**, *81*, 063005. [CrossRef]

38. Sakstein, J.; Jain, B.; Heyl, J.S.; Hui, L. Tests of Gravity Theories Using Supermassive Black Holes. *Astrophys. J. Lett.* **2017**, *844*, L14. [CrossRef]

39. Falck, B.; Koyama, K.; Zhao, G.B.; Li, B. The Vainshtein Mechanism in the Cosmic Web. *J. Cosmol. Astropart. Phys.* **2014**, *07*, 058. [CrossRef]

40. Cardoso, A.; Koyama, K.; Seahra, S.S.; Silva, F.P. Cosmological perturbations in the DGP braneworld: Numeric solution. *Phys. Rev. D* **2008**, *77*, 083512. [CrossRef]

41. Schmidt, F. Self-consistent cosmological simulations of DGP braneworld gravity. *Phys. Rev. D* **2009**, *80*, 043001. [CrossRef]

42. Khoury, J.; Wyman, M. N-body simulations of DGP and degravitation theories. *Phys. Rev. D* **2009**, *80*, 064023. [CrossRef]

43. Chan, K.C.; Scoccimarro, R. Large-scale structure in brane-induced gravity. II. Numerical simulations. *Phys. Rev. D* **2009**, *80*, 104005. [CrossRef]

44. Davis, A.C.; Li, B.; Mota, D.F.; Winther, H.A. Structure Formation in the Symmetron Model. *Astrophys. J.* **2012**, *748*, 61. [CrossRef]

45. Clampitt, J.; Jain, B.; Khoury, J. Halo Scale Predictions of Symmetron Modified Gravity. *J. Cosmol. Astropart. Phys.* **2012**, *1201*, 030. [CrossRef]

46. Brax, P.; Davis, A.C.; Li, B.; Winther, H.A.; Zhao, G.B. Systematic Simulations of Modified Gravity: Symmetron and Dilaton Models. *J. Cosmol. Astropart. Phys.* **2012**, *1210*, 002. [CrossRef]

47. Llinares, C.; Mota, D.F. Cosmological simulations of screened modified gravity out of the static approximation: effects on matter distribution. *Phys. Rev.* **2014**, *D89*, 084023. [CrossRef]

48. Noller, J.; von Braun-Bates, F.; Ferreira, P.G. Relativistic scalar fields and the quasistatic approximation in theories of modified gravity. *Phys. Rev.* **2014**, *D89*, 023521. [CrossRef]

49. Vardanyan, V. Aspects of Cosmic Acceleration. Ph.D. Thesis, Leiden University, Leiden, The Netherlands, 2019.

50. Contigiani, O.; Vardanyan, V.; Silvestri, A. Splashback radius in symmetron gravity. *Phys. Rev. D* **2019**, *99*, 064030. [CrossRef]

51. Upadhye, A. Symmetron dark energy in laboratory experiments. *Phys. Rev. Lett.* **2013**, *110*, 031301. [CrossRef]

52. Li, B.; Zhao, G.B.; Teyssier, R.; Koyama, K. ECOSMOG: an Efficient COde for Simulating MOdified Gravity. *J. Cosmol. Astropart. Phys.* **2012**, *2012*, 051. [CrossRef]

53. Llinares, C.; Mota, D.F.; Winther, H.A. ISIS: A new N-body cosmological code with scalar fields based on RAMSES. Code presentation and application to the shapes of clusters. *Astron. Astrophys.* **2014**, *562*, A78. [CrossRef]

54. Puchwein, E.; Baldi, M.; Springel, V. Modified Gravity-GADGET: A new code for cosmological hydrodynamical simulations of modified gravity models. *Mon. Not. R. Astron. Soc.* **2013**, *436*, 348. [CrossRef]

55. Winther, H.A.; Schmidt, F.; Barreira, A.; Arnold, C.; Bose, S.; Llinares, C.; Baldi, M.; Falck, B.; Hellwing, W.A.; Koyama, K.; et al. Modified Gravity N-body Code Comparison Project. *Mon. Not. R. Astron. Soc.* **2015**, *454*, 4208–4234. [CrossRef]

56. Briddon, C.; Burrage, C.; Moss, A.; Tamosiunas, A. SELCIE: a tool for investigating the chameleon field of arbitrary sources. *J. Cosmol. Astropart. Phys.* **2021**, *2021*, 043. [CrossRef]
57. Alnaes, M.S.; Blechta, J.; Hake, J.; Johansson, A.; Kehlet, B.; Logg, A.; Richardson, C.; Ring, J.; Rognes, M.E.; Wells, G.N. The FEniCS Project Version 1.5. *Arch. Numer. Softw.* **2015**, *3*, 9–23. [CrossRef]

58. Logg, A.; Mardal, K.A.; Wells, G. *Automated Solution of Differential Equations by the Finite Element Method*; Springer: Berlin, Germany, 2012. [CrossRef]

59. Cimrman, R.; Lukeš, V.; Rohan, E. Multiscale finite element calculations in Python using SfePy. *arXiv* **2018**, arXiv:1810.00674. https://doi.org/10.48550/arXiv.1810.00674.

60. Braden, J.; Burrage, C.; Elder, B.; Saadeh, D. *qenics: Vainshtein screening with the finite element method*. *J. Cosmol. Astropart. Phys.* **2021**, *2021*, 010. [CrossRef]

61. Tamosiunas, A.; Briddon, C.; Burrage, C.; Cutforth, A.; Moss, A.; Vincent, T. Chameleon screening in cosmic voids. *J. Cosmol. Astropart. Phys.* **2022**, *2022*, 056. [CrossRef]

62. Tamosiunas, A.; Briddon, C.; Burrage, C.; Cui, W.; Moss, A. Chameleon screening depends on the shape and structure of NFW halos. *J. Cosmol. Astropart. Phys.* **2022**, *2022*, 047. [CrossRef]

63. Lévy, H.; Bergé, J.; Uzan, J.P. Solving nonlinear Klein-Gordon equations on unbounded domains via the finite element method. *Phys. Rev. D* **2022**, *106*, 124021. [CrossRef]

64. Grosch, C.E.; Orszag, S.A. Numerical Solution of Problems in Unbounded Regions: Coordinate Transforms. *J. Comput. Phys.* **1977**, *25*, 273–295. [CrossRef]

65. Burrage, C.; Copeland, E.J.; Moss, A.; Stevenson, J.A. The shape dependence of chameleon screening. *J. Cosmol. Astropart. Phys.* **2018**, *2018*, 056. [CrossRef]

66. Burrage, C.; Copeland, E.J.; Hinds, E.A. Probing Dark Energy with Atom Interferometry. *J. Cosmol. Astropart. Phys.* **2015**, *2015*, 042. [CrossRef]

67. Burrage, C.; Copeland, E.J. Using Atom Interferometry to Detect Dark Energy. *Contemp. Phys.* **2016**, *57*, 164–176. [CrossRef]

68. Hamilton, P.; Jaffe, M.; Haslinger, P.; Simmons, Q.; Müller, H.; Khoury, J. Atom-interferometry constraints on dark energy. *Science* **2015**, *349*, 849–851. [CrossRef]

69. Burrage, C.; Kuribayashi-Coleman, A.; Stevenson, J.; Thrussell, B. Constraining symmetron fields with atom interferometry. *J. Cosmol. Astropart. Phys.* **2016**, *1612*, 041. [CrossRef]

70. Jaffe, M.; Haslinger, P.; Xu, V.; Hamilton, P.; Upadhye, A.; Elder, B.; Khoury, J.; Müller, H. Testing sub-gravitational forces on atoms from a miniature, in-vacuum source mass. *Nat. Phys.* **2017**, *13*, 938. [CrossRef]

71. Sabulsky, D.O.; Dutta, I.; Hinds, E.A.; Elder, B.; Burrage, C.; Copeland, E.J. Experiment to detect dark energy forces using atom interferometry. *Phys. Rev. Lett.* **2019**, *123*, 061102. [CrossRef]

72. Brax, P.; Burrage, C. Atomic Precision Tests and Light Scalar Couplings. *Phys. Rev. D* **2011**, *83*, 035020. [CrossRef]

73. Brax, P.; Davis, A.C.; Elder, B. Screened scalar fields in hydrogen and muonium. *Phys. Rev. D* **2023**, *107*, 044008. [CrossRef]

74. Cronenberg, G.; Brax, P.; Filter, H.; Geltenbort, P.; Jenke, T.; Pignol, G.; Pitschmann, M.; Thalhammer, M.; Abele, H. Acoustic Rabi oscillations between gravitational quantum states and impact on symmetron dark energy. *Nat. Phys.* **2018**, *14*, 1022–1026. [CrossRef]

75. Decca, R.S.; Lopez, D.; Chan, H.B.; Fischbach, E.; Krause, D.E.; Jamell, C.R. Constraining new forces in the Casimir regime using the isoelectronic technique. *Phys. Rev. Lett.* **2005**, *94*, 240401. [CrossRef]

76. Decca, R.S.; Lopez, D.; Fischbach, E.; Krause, D.E. Measurement of the Casimir Force between Dissimilar Metals. *Phys. Rev. Lett.* **2003**, *91*, 050402. [CrossRef]

77. Almasi, A.; Brax, P.; Iannuzzi, D.; Sedmik, R.I.P. Force sensor for chameleon and Casimir force experiments with parallel-plate configuration. *Phys. Rev.* **2015**, *D91*, 102002. [CrossRef]

78. Sedmik, R.; Brax, P. Status Report and first Light from Cannex: Casimir Force Measurements between flat parallel Plates. In *Journal of Physics: Conference Series, Proceedings of the 10th International Conference on Precision Physics of Simple Atomic Systems (PSAS 2018)*, Vienna, Austria, 14–18 May 2018; IOP Publishing Ltd.: Bristol, UK, 2018; Volume 1138, p. 012014. [CrossRef]

79. Brax, P.; van de Bruck, C.; Davis, A.C.; Shaw, D.J.; Iannuzzi, D. Tuning the Mass of Chameleon Fields in Casimir Force Experiments. *Phys. Rev. Lett.* **2010**, *104*, 241101. [CrossRef] [PubMed]

80. Adelberger, E.G. Sub-millimeter tests of the gravitational inverse square law. In *Proceedings of the 2nd Meeting on CPT and Lorentz Symmetry*, Bloomington, IN, USA, 15–18 August 2002; pp. 9–15. [CrossRef]

81. Brax, P.; Davis, A.C.; Elder, B. Casimir Tests of Scalar-Tensor Theories. *Phys. Rev. D* **2022**, *107*, 084025. [CrossRef]

82. Sparnaay, M.J. Measurements of attractive forces between flat plates. *Physica* **1958**, *24*, 751–764. [CrossRef]

83. van Blokland, P.H.G.M.; Overbeek, J.T.G. van der Waals forces between objects covered with a chromium layer. *J. Chem. Soc. Faraday Trans. 1* **1978**, *74*, 2637–2651. [CrossRef]

84. Lamoreaux, S.K. Demonstration of the Casimir force in the 0.6 to 6 micrometers range. *Phys. Rev. Lett.* **1997**, *78*, 5–8; Erratum in *Phys. Rev. Lett.* **1998**, *81*, 5475. [CrossRef]

85. Decca, R.S.; Lopez, D.; Fischbach, E.; Klimchitskaya, G.L.; Krause, D.E.; Mostepanenko, V.M. Precise comparison of theory and new experiment for the Casimir force leads to stronger constraints on thermal quantum effects and long-range interactions. *Annals Phys.* **2005**, *318*, 37–80. [CrossRef]

86. Mohideen, U.; Roy, A. Precision measurement of the Casimir force from 0.1 to 0.9 micrometers. *Phys. Rev. Lett.* **1998**, *81*, 4549–4552. [CrossRef]
87. Harris, B.W.; Chen, F.; Mohideen, U. Precision measurement of the Casimir force using gold surfaces. *Phys. Rev. A* **2000**, *62*, 052109. [CrossRef]

88. Chen, Y.J.; Tham, W.K.; Krause, D.E.; Lopez, D.; Fischbach, E.; Decca, R.S. Stronger Limits on Hypothetical Yukawa Interactions in the 30-8000 nm Range. *Phys. Rev. Lett.* **2016**, *116*, 221102. [CrossRef] [PubMed]

89. Bressi, G.; Carugno, G.; Onofrio, R.; Ruoso, G. Measurement of the Casimir force between parallel metallic surfaces. *Phys. Rev. Lett.* **2002**, *88*, 041804. [CrossRef] [PubMed]

90. Zou, J.; Marcet, Z.; Rodriguez, A.W.; Reid, M.T.H.; McCauley, A.P.; Kravchenko, I.I.; Lu, T.; Bao, Y.; Johnson, S.G.; Chan, H.B. Casimir forces on a silicon micromechanical chip. *Nat. Commun.* **2013**, *4*, 1845. [CrossRef]

91. Neto, P.A.M.; Lambrecht, A.; Reynaud, S. Roughness correction to the Casimir force: Beyond the Proximity Force Approximation. *EPL (Europhys. Lett.)* **2005**, *69*, 924–930. [CrossRef]

92. Rodrigues, R.B.; Neto, P.A.M.; Lambrecht, A.; Reynaud, S. Lateral Casimir Force beyond the Proximity-Force Approximation. *Phys. Rev. Lett.* **2006**, *96*, 100402. [CrossRef]

93. Lambrecht, A.; Reynaud, S. Casimir force between metallic mirrors. *Eur. Phys. J. D* **2000**, *8*, 309–318. [CrossRef]

94. Reynaud, S.; Lambrecht, A. 407Casimir forces and vacuum energy. In *Quantum Optics and Nanophotonics*; Oxford University Press: Oxford, UK, 2017. [CrossRef]

95. Dalvit, D.A.R.; Onofrio, R. On the use of the proximity force approximation for deriving limits to short-range gravitational-like interactions from sphere-plane Casimir force experiments. *Phys. Rev. D* **2009**, *80*, 064025. [CrossRef]

96. Krause, D.E.; Decca, R.S.; Lopez, D.; Fischbach, E. Experimental Investigation of the Casimir Force beyond the Proximity-Force Approximation. *Phys. Rev. Lett.* **2007**, *98*, 050403. [CrossRef]

97. Brax, P.; van de Bruck, C.; Davis, A.C.; Mota, D.F.; Shaw, D.J. Detecting chameleons through Casimir force measurements. *Phys. Rev.* **2007**, D76, 124034. [CrossRef]

98. Brax, P.; Davis, A.C. Casimir, Gravitational and Neutron Tests of Dark Energy. *Phys. Rev. D* **2015**, *91*, 063503. [CrossRef]

99. Upadhye, A.; Gubser, S.S.; Khoury, J. Unveiling chameleons in tests of gravitational inverse-square law. *Phys. Rev. D* **2006**, *74*, 104024. [CrossRef]

100. Upadhye, A. Dark energy fifth forces in torsion pendulum experiments. *Phys. Rev. D* **2012**, *86*, 102003. [CrossRef]

101. Elder, B.; Khoury, J.; Haslinger, P.; Jaffe, M.; Müller, H.; Hamilton, P. Chameleon Dark Energy and Atom Interferometry. *Phys. Rev. D* **2016**, *94*, 044051. [CrossRef]

102. Koyama, K.; Sakstein, J. Astrophysical probes of the Vainshtein mechanism: Stars and galaxies. *Phys. Rev. D* **2015**, *91*, 124066. [CrossRef]

103. Saito, R.; Yamauchi, D.; Mizuno, S.; Gleyzes, J.; Langlois, D. Modified gravity inside astrophysical bodies. *J. Cosmol. Astropart. Phys.* **2015**, *2015*, 008. [CrossRef]

104. Gleyzes, J.; Langlois, D.; Piazza, F.; Vernizzi, F. Exploring gravitational theories beyond Horndeski. *J. Cosmol. Astropart. Phys.* **2015**, *2015*, 018. [CrossRef]

105. Sakstein, J. Hydrogen Burning in Low Mass Stars Constrains Scalar-Tensor Theories of Gravity. *Phys. Rev. Lett.* **2015**, *115*, 201101,

106. Sakstein, J. Testing gravity using dwarf stars. *Phys. Rev. D* **2015**, *92*, 124045. [CrossRef]

107. Ségransan, D.; Delfosse, X.; Forveille, T.; Beuzit, J.L.; Udry, S.; Perrier, C.; Mayor, M. Accurate masses of very low mass stars. III. 16 new or improved masses. *Astron. Astrophys.* **2000**, *364*, 665–673.

108. Paxton, B.; Bildsten, L.; Dotter, A.; Herwig, F.; Lesaffre, P.; Timmes, F. Modules for Experiments in Stellar Astrophysics (MESA). *Astrophys. J. Suppl. Ser.* **2011**, *192*, 3. [CrossRef]

109. Paxton, B.; Cantiello, M.; Arras, P.; Bildsten, L.; Brown, E.F.; Dotter, A.; Mankovich, C.; Montgomery, M.H.; Stello, D.; Timmes, F.X.; et al. Modules for Experiments in Stellar Astrophysics (MESA): Planets, Oscillations, Rotation, and Massive Stars. *Astrophys. J. Suppl. Ser.* **2013**, *208*, 4. [CrossRef]

110. Paxton, B.; Marchant, P.; Schwab, J.; Bauer, E.B.; Bildsten, L.; Cantiello, M.; Dessart, L.; Farmer, R.; Hu, H.; Langer, N.; et al. Modules for Experiments in Stellar Astrophysics (MESA): Binaries, Pulsations, and Explosions. *Astrophys. J. Suppl. Ser.* **2015**, *220*, 15. [CrossRef]

111. Paxton, B.; Schwab, J.; Bauer, E.B.; Bildsten, L.; Blinnikov, S.; Duffell, P.; Farmer, R.; Goldberg, J.A.; Marchant, P.; Sorokina, E.; et al. Modules for Experiments in Stellar Astrophysics (MESA): Convective Boundaries, Element Diffusion, and Massive Star Explosions. *Astrophys. J. Suppl. Ser.* **2018**, *234*, 34. [CrossRef]

112. Paxton, B.; Smolec, R.; Schwab, J.; Gautschy, A.; Bildsten, L.; Cantiello, M.; Dotter, A.; Farmer, R.; Goldberg, J.A.; Jermyn, A.S.; et al. Modules for Experiments in Stellar Astrophysics (MESA): Pulsating Variable Stars, Rotation, Convective Boundaries, and Energy Conservation. *Astrophys. J. Suppl. Ser.* **2019**, *243*, 10. [CrossRef]

113. Jermyn, A.S.; Bauer, E.B.; Schwab, J.; Farmer, R.; Ball, W.H.; Bellinger, E.P.; Dotter, A.; Joyce, M.; Marchant, P.; Mombarg, J.S.G.; et al. Modules for Experiments in Stellar Astrophysics (MESA): Time-Dependent Convection, Energy Conservation, Automatic Differentiation, and Infrastructure. *arXiv* **2022**, arXiv:2208.03651. [CrossRef]

114. Chang, P.; Hui, L. Stellar Structure and Tests of Modified Gravity. *Astrophys. J.* **2011**, *732*, 25. [CrossRef]

115. Saltas, I.D.; Lopes, I. Obtaining Precision Constraints on Modified Gravity with Helioseismology. *Phys. Rev. Lett.* **2019**, *123*, 091103. [CrossRef]

116. Saltas, I.D.; Christensen-Dalsgaard, J. Searching for dark energy with the Sun. *Astron. Astrophys.* **2022**, *667*, A115. [CrossRef]
117. Jain, B.; Vikram, V.; Sakstein, J. Astrophysical Tests of Modified Gravity: Constraints from Distance Indicators in the Nearby Universe. Astrophys. J. **2013**, 779, 39. [CrossRef]

118. Sakstein, J. Stellar oscillations in modified gravity. Phys. Rev. D **2013**, 88, 124013. [CrossRef]

119. Desmond, H.; Sakstein, J.; Jain, B. Five percent measurement of the gravitational constant in the Large Magellanic Cloud. Phys. Rev. D **2021**, 103, 024028. [CrossRef]

120. Desmond, H.; Jain, B.; Sakstein, J. Local resolution of the Hubble tension: The impact of screened fifth forces on the cosmic distance ladder. Phys. Rev. D **2019**, 100, 043537. [CrossRef]

121. Desmond, H.; Sakstein, J. Screened fifth forces lower the TRGB-calibrated Hubble constant too. Phys. Rev. D **2020**, 102, 023007. [CrossRef]

122. Sakstein, J.; Desmond, H.; Jain, B. Screened fifth forces mediated by dark matter-baryon interactions: Theory and astrophysical probes. Phys. Rev. D **2019**, 100, 104035. [CrossRef]

123. Högås, M.; Mörtsell, E. Impact of symmetron screening on the Hubble tension: New constraints using cosmic distance ladder data. Phys. Rev. D **2023**, 108, 024007. [CrossRef]

124. Sakstein, J.; Kenna-Allison, M.; Koyama, K. Stellar pulsations in beyond Horndeski gravity theories. J. Cosmol. Astropart. Phys. **2017**, 2017, 007. [CrossRef]

125. Zhao, G.B.; Li, B.; Koyama, K. Testing Gravity Using the Environmental Dependence of Dark Matter Halos. Phys. Rev. Lett. **2011**, 107, 071303. [CrossRef]

126. Schmidt, F. Dynamical masses in modified gravity. Phys. Rev. D **2010**, 81, 103002. [CrossRef]

127. Cabré, A.; Vikram, V.; Zhao, G.B.; Jain, B.; Koyama, K. Astrophysical tests of gravity: a screening map of the nearby universe. J. Cosmol. Astropart. Phys. **2012**, 2012, 034. [CrossRef]

128. Yang, X.; Mo, H.J.; van den Bosch, F.C.; Pasquali, A.; Li, C.; Barden, M. Galaxy Groups in the SDSS DR4. I. The Catalog and Basic Properties. Astrophys. J. **2007**, 671, 153–170. [CrossRef]

129. Abell, G.O.; Corwin, H.G., Jr.; Olowin, R.P. A Catalog of Rich Clusters of Galaxies. Astrophys. J. Suppl. Ser. **1989**, 70, 1. [CrossRef]

130. Ebeling, H.; Voges, W.; Bohringer, H.; Edge, A.C.; Huchra, J.P.; Briel, U.G. Properties of the X-ray-brightest Abell-type clusters of galaxies (XBACs) from ROSAT All-Sky Survey data—I. The sample. Mon. Not. R. Astron. Soc. **1996**, 281, 799–829. [CrossRef]

131. Karachentsev, I.D.; Karachentseva, V.E.; Huchtmeier, W.K.; Makarov, D.I. A Catalog of Neighboring Galaxies. Astron. J. **2004**, 127, 2031–2068. [CrossRef]

132. Lavaux, G.; Hudson, M.J. The 2M++ galaxy redshift catalogue. Mon. Not. R. Astron. Soc. **2011**, 416, 2840–2856. [CrossRef]

133. Zhao, G.B.; Li, B.; Koyama, K. N-body simulations for f(R) gravity using a self-adaptive particle-mesh code. Phys. Rev. D **2011**, 83, 044007. [CrossRef]

134. Desmond, H.; Ferreira, P.G.; Lavaux, G.; Jasche, J. Reconstructing the gravitational field of the local Universe. Mon. Not. R. Astron. Soc. **2018**, 474, 3152–3161. [CrossRef]

135. Behroozi, P.S.; Wechsler, R.H.; Wu, H.Y. The ROCKSTAR Phase-space Temporal Halo Finder and the Velocity Offsets of Cluster Cores. Astrophys. J. **2013**, 762, 109. [CrossRef]

136. Skillman, S.W.; Warren, M.S.; Turk, M.J.; Wechsler, R.H.; Holz, D.E.; Sutter, P.M. Dark Sky Simulations: Early Data Release. arXiv **2014**, arXiv:1407.2600.

137. Lehmann, B.V.; Mao, Y.Y.; Becker, M.R.; Skillman, S.W.; Wechsler, R.H. The Concentration Dependence of the Galaxy-Halo Connection: Modeling Assembly Bias with Abundance Matching. Astrophys. J. **2017**, 834, 37. [CrossRef]

138. Navarro, J.F.; Frenk, C.S.; White, S.D.M. A Universal Density Profile from Hierarchical Clustering. Astrophys. J. **1997**, 490, 493–508. [CrossRef]

139. Cole, S.; Lacey, C. The structure of dark matter haloes in hierarchical clustering models. Mon. Not. R. Astron. Soc. **1996**, 281, 716. [CrossRef]

140. Jasche, J.; Wandelt, B.D. Bayesian inference from photometric redshift surveys. Mon. Not. R. Astron. Soc. **2012**, 425, 1042–1056. [CrossRef]

141. Jasche, J.; Wandelt, B.D. Bayesian physical reconstruction of initial conditions from large-scale structure surveys. Mon. Not. R. Astron. Soc. **2013**, 432, 894–913. [CrossRef]

142. Jasche, J.; Kitaura, F.S.; Wandelt, B.D.; Enßlin, T.A. Bayesian power-spectrum inference for large-scale structure data. Mon. Not. R. Astron. Soc. **2010**, 406, 60–85. [CrossRef]

143. Jasche, J.; Leclercq, F.; Wandelt, B.D. Past and present cosmic structure in the SDSS DR7 main sample. J. Cosmol. Astropart. Phys. **2015**, 2015, 036. [CrossRef]

144. Lavaux, G.; Jasche, J. Unmasking the masked Universe: the 2M++ catalogue through Bayesian eyes. Mon. Not. R. Astron. Soc. **2016**, 455, 3169–3179. [CrossRef]

145. Shao, S.; Li, B.; Cautun, M.; Wang, H.; Wang, J. Screening maps of the local Universe I - Methodology. Mon. Not. R. Astron. Soc. **2019**, 489, 4912–4925. [CrossRef]

146. Li, B.; Zhao, G.B.; Koyama, K. Exploring Vainshtein mechanism on adaptively refined meshes. J. Cosmol. Astropart. Phys. **2013**, 2013, 023. [CrossRef]

147. Li, B.; Barreira, A.; Baugh, C.M.; Hellwing, W.A.; Koyama, K.; Pascoli, S.; Zhao, G.B. Simulating the quartic Galileon gravity model on adaptively refined meshes. J. Cosmol. Astropart. Phys. **2013**, 2013, 012. [CrossRef]
148. Wang, H.; Mo, H.J.; Yang, X.; Jing, Y.P.; Lin, W.P. ELUCID—Exploring the Local Universe with the Reconstructed Initial Density Field. I. Hamiltonian Markov Chain Monte Carlo Method with Particle Mesh Dynamics. *Astrophys. J.* **2014**, *794*, 94. [CrossRef]

149. Wang, H.; Mo, H.J.; Yang, X.; Zhang, Y.; Shi, J.; Jing, Y.P.; Liu, C.; Li, S.; Kang, X.; Gao, Y. ELUCID—Exploring the Local Universe with Reconstructed Initial Density Field III: Constrained Simulation in the SDSS Volume. *Astrophys. J.* **2016**, *831*, 164. [CrossRef]

150. Naidoo, K.; Hellwing, W.; Bilicki, M.; Libeskind, N.; Pfeifer, S.; Hoffman, Y. Constrained simulations of the local Universe with Modified Gravity. *arXiv* **2022**, arXiv:2209.14386.

151. Kourkchi, E.; Tully, R.B.; Eftekharzadeh, S.; Llop, J.; Courtois, H.M.; Guinet, D.; Dupuy, A.; Neill, J.D.; Seibert, M.; Andrews, M.; et al. Cosmicflows-4: The Catalog of ~10,000 Tully-Fisher Distances. *Astrophys. J.* **2020**, *902*, 145. [CrossRef]

152. Doumler, T.; Gottlöber, S.; Hoffman, Y.; Courtois, H. Reconstructing cosmological initial conditions from galaxy peculiar velocities—III. Constrained simulations. *Mon. Not. R. Astron. Soc.* **2013**, *430*, 912–923. [CrossRef]

153. Zaroubi, S.; Hoffman, Y.; Fisher, K.B.; Lahav, O. Wiener Reconstruction of the Large-Scale Structure. *Astrophys. J.* **1995**, *449*, 446. [CrossRef]

154. Doumler, T.; Hoffman, Y.; Courtois, H.; Gottlöber, S. Reconstructing cosmological initial conditions from galaxy peculiar velocities—I. Reverse Zeldovich Approximation. *Mon. Not. R. Astron. Soc.* **2013**, *430*, 888–901. [CrossRef]

155. Winther, H.A.; Koyama, K.; Manera, M.; Wright, B.S.; Zhao, G.B. COLA with scale-dependent growth: applications to screened modified gravity models. *J. Cosmol. Astropart. Phys.* **2017**, *2017*, 006. [CrossRef]

156. Tassev, S.; Zaldarriaga, M.; Eisenstein, D.J. Solving large scale structure in ten easy steps with COLA. *J. Cosmol. Astropart. Phys.* **2013**, *2013*, 036. [CrossRef]

157. Fasiello, M.; Vlah, Z. Screening in perturbative approaches to LSS. *Phys. Lett. B* **2017**, *773*, 236–241. [CrossRef]

158. Jain, B.; VanderPlas, J. Tests of modified gravity with dwarf galaxies. *J. Cosmol. Astropart. Phys.* **2011**, *2011*, 032. [CrossRef]

159. Vikram, V.; Cabré, A.; Jain, B.; VanderPlas, J.T. Astrophysical tests of modified gravity: the morphology and kinematics of dwarf galaxies. *J. Cosmol. Astropart. Phys.* **2013**, *2013*, 020. [CrossRef]

160. Giovanelli, R.; Haynes, M.P.; Kent, B.R.; Perillat, P.; Saintonge, A.; Brosch, N.; Catinella, B.; Hoffman, G.L.; Stierwalt, S.; Spekkens, K.; et al. The Arecibo Legacy Fast ALFA Survey. 1. Science goals, survey design and strategy. *Astron. J.* **2005**, *130*, 2598–2612. [CrossRef]

161. Abazajian, K.N.; Adelman-McCarthy, J.K.; Agüeros, M.A.; Allam, S.S.; Prieto, C.A.; An, D.; Anderson, K.S.; Anderson, S.F.; Annis, J.; Bahcall, N.A.; et al. The Seventh Data Release of the Sloan Digital Sky Survey. *Astrophys. J. Suppl.* **2009**, *182*, 543–558. [CrossRef]

162. Desmond, H.; Ferreira, P.G.; Lavaux, G.; Jasche, J. Fifth force constraints from the separation of galaxy mass components. *Phys. Rev. D* **2018**, *98*, 064015. [CrossRef]

163. Desmond, H.; Ferreira, P.G.; Lavaux, G.; Jasche, J. Fifth force constraints from galaxy warps. *Phys. Rev. D* **2018**, *98*, 083010. [CrossRef]

164. Desmond, H.; Ferreira, P.G.; Lavaux, G.; Jasche, J. The Fifth Force in the Local Cosmic Web. *Mon. Not. R. Astron. Soc.* **2019**, *483*, L64–L68. [CrossRef]

165. Desmond, H.; Ferreira, P.G. Galaxy morphology rules out astrophysically relevant Hu-Sawicki f (R) gravity. *Phys. Rev. D* **2020**, *102*, 104060. [CrossRef]

166. Bartlett, D.J.; Desmond, H.; Ferreira, P.G. Calibrating galaxy formation effects in galactic tests of fundamental physics. *Phys. Rev. D* **2021**, *103*, 123502. [CrossRef]

167. Asvathaman, A.; Heyl, J.S.; Hui, L. Eötvös experiments with supermassive black holes. *Mon. Not. R. Astron. Soc.* **2017**, *465*, 3261–3266. [CrossRef]

168. Bartlett, D.J.; Desmond, H.; Ferreira, P.G. Constraints on Galileons from the positions of supermassive black holes. *Phys. Rev. D* **2021**, *103*, 023523. [CrossRef]

169. Bartlett, D.J.; Desmond, H.; Devriendt, J.; Ferreira, P.G.; Slyz, A. Spatially offset black holes in the Horizon-AGN simulation and comparison to observations. *Mon. Not. R. Astron. Soc.* **2021**, *500*, 4639–4657. [CrossRef]

170. Khoury, J. Les Houches Lectures on Physics Beyond the Standard Model of Cosmology. *arXiv* **2013**, arXiv:1312.2006.

171. Murphy, T.W., Jr.; Adelberger, E.G.; Battat, J.B.R.; Hoyle, C.D.; Johnson, N.H.; McMillan, R.J.; Stubbs, C.W.; Swanson, H.E. APOLLO: millimeter lunar laser ranging. *Class. Quantum Gravity* **2012**, *29*, 184005. [CrossRef]

172. Schmidt, F.; Lima, M.V.; Oyaizu, H.; Hu, W. Non-linear Evolution of f(R) Cosmologies III: Halo Statistics. *Phys. Rev. D* **2009**, *79*, 083518. [CrossRef]

173. Li, Y.; Hu, W. Chameleon Halo Modeling in f(R) Gravity. *Phys. Rev. D* **2011**, *84*, 084033. [CrossRef]

174. Pourhasan, R.; Afshordi, N.; Mann, R.B.; Davis, A.C. Chameleon Gravity, Electrostatics, and Kinematics in the Outer Galaxy. *J. Cosmol. Astropart. Phys.* **2011**, *12*, 005. [CrossRef]

175. Lombriser, L.; Koyama, K.; Zhao, G.B.; Li, B. Chameleon f(R) gravity in the virialized cluster. *Phys. Rev. D* **2012**, *85*, 124054. [CrossRef]

176. Lombriser, L.; Li, B.; Koyama, K.; Zhao, G.B. Modeling halo mass functions in chameleon f(R) gravity. *Phys. Rev. D* **2013**, *87*, 123511. [CrossRef]

177. Lombriser, L.; Koyama, K.; Li, B. Halo modelling in chameleon theories. *J. Cosmol. Astropart. Phys.* **2014**, *03*, 021. [CrossRef]

178. Shi, D.; Li, B.; Han, J.; Gao, L.; Hellwing, W.A. Exploring the liminality: properties of haloes and subhaloes in borderline f(R) gravity. *Mon. Not. R. Astron. Soc.* **2015**, *452*, 3179–3191. [CrossRef]
179. Mitchell, M.A.; Arnold, C.; He, J.H.; Li, B. A general framework to test gravity using galaxy clusters II: A universal model for the halo concentration in $f(R)$ gravity. *Mon. Not. R. Astron. Soc.* **2019**, *487*, 1410–1425. [CrossRef]

180. Barreira, A.; Li, B.; Hellwing, W.A.; Lombriser, L.; Baugh, C.M.; Pascoli, S. Halo model and halo properties in Galileon gravity cosmologies. *J. Cosmol. Astropart. Phys.* **2014**, *2014*, 029. [CrossRef]

181. Mitchell, M.A.; Hernández-Aguayo, C.; Arnold, C.; Li, B. A general framework to test gravity using galaxy clusters IV: cluster and halo properties in DGP gravity. *Mon. Not. R. Astron. Soc.* **2021**, *508*, 4140–4156. [CrossRef]

182. Pizzuti, L.; Saltas, I.D.; Umetsu, K.; Sartoris, B. Probing vainsthein-screening gravity with galaxy clusters using internal kinematics and strong and weak lensing. *Mon. Not. Roy. Astron. Soc.* **2022**, *512*, 4280–4290. [CrossRef]

183. Falck, B.; Koyama, K.; Zhao, G.B. Cosmic Web and Environmental Dependence of Screening: Vainshtein vs. Chameleon. *J. Cosmol. Astropart. Phys.* **2015**, *2015*, 049. [CrossRef]

184. Diemer, B.; Kravtsov, A.V. Dependence of the outer density profiles of halos on their mass accretion rate. *Astrophys. J.* **2014**, *789*, 1. [CrossRef]

185. More, S.; Miyatake, H.; Takada, M.; Diemer, B.; Kravtsov, A.V.; Dalal, N.K.; More, A.; Murata, R.; Melbaum, R.; Rozo, E.; et al. Detection of the Splashback Radius and Halo Assembly bias of Massive Galaxy Clusters. *Astrophys. J.* **2016**, *825*, 39. [CrossRef]

186. Baxter, E.; Chang, C.; Jain, B.; Adhikari, S.; Dalal, N.; Kravtsov, A.; More, S.; Rozo, E.; Rykoff, E.; Sheth, R.K. The Halo Boundary of Galaxy Clusters in the SDSS. *Astrophys. J.* **2017**, *841*, 18. [CrossRef]

187. Chang, C.; Baxter, E.; Jain, B.; Sánchez, C.; Adhikari, S.; Varga, T.N.; Fang, Y.; Rozo, E.; Rykoff, E.S.; Kravtsov, A.; et al. The Splashback Feature around DES Galaxy Clusters: Galaxy Density and Weak Lensing Profiles. *Astrophys. J.* **2018**, *864*, 83. [CrossRef]

188. Shin, T.H.; Adhikari, S.; Baxter, E.J.; Chang, C.; Jain, B.; Battaglia, N.; Bleem, L.; Bocquet, S.; DeRose, J.; Gruen, D.; et al. Measurement of the Splashback Feature around SZ-selected Galaxy Clusters with DES, SPT and ACT. *Mon. Not. R. Astron. Soc.* **2018**, *487*, 2900–2918. [CrossRef]

189. Umetsu, K.; Diemer, B. Lensing Constraints on the Mass Profile Shape and the Splashback Radius of Galaxy Clusters. *Astrophys. J.* **2017**, *836*, 231. [CrossRef]

190. Contigiani, O.; Hoekstra, H.; Bahé, Y.M. Weak lensing constraints on splashback around massive clusters. *Mon. Not. R. Astron. Soc.* **2018**, *485*, 408–415. [CrossRef]

191. Adhikari, S.; Sakstein, J.; Jain, B.; Dalal, N.; Li, B. Splashback in galaxy clusters as a probe of cosmic expansion and gravity. *J. Cosmol. Astropart. Phys.* **2018**, *11*, 033. [CrossRef]

192. Fillmore, J.A.; Goldreich, P. Self-similar gravitational collapse in an expanding universe. *Astrophys. J.* **1984**, *281*, 1–8. [CrossRef]

193. Bertschinger, E. Self-similar secondary infall and accretion in an Einstein-de Sitter universe. *Astrophys. J. Suppl. Ser.* **1985**, *58*, 39–65. [CrossRef]

194. Lithwick, Y.; Dalal, N. Self-Similar Solutions of Triaxial Dark Matter Halos. *Astrophys. J.* **2011**, *734*, 100. [CrossRef]

195. Shi, X. The outer profile of dark matter haloes: an analytical approach. *Mon. Not. R. Astron. Soc.* **2016**, *459*, 3711–3720. [CrossRef]

196. Amendola, L.; Appleby, S.; Avgoustidis, A.; Bacon, D.; Baker, T.; Baldi, M.; Bartolo, N.; Blanchard, A.; Bonvin, C.; Borgani, S.; et al. Cosmology and fundamental physics with the Euclid satellite. *Living Rev. Relativ.* **2018**, *21*, 2. [CrossRef]

197. Gubitosi, G.; Piazza, F.; Vernizzi, F. The Effective Field Theory of Dark Energy. *J. Cosmol. Astropart. Phys.* **2013**, *2013*, 032. [CrossRef]

198. Gleyzes, J.; Langlois, D.; Piazza, F.; Vernizzi, F. Essential Building Blocks of Dark Energy. *J. Cosmol. Astropart. Phys.* **2013**, *2013*, 025. [CrossRef]

199. Creminelli, P.; Gleyzes, J.; Hui, L.; Simonović, M.; Vernizzi, F. Single-Field Consistency Relations of Large Scale Structure. Part III: Test of the Equivalence Principle. *J. Cosmol. Astropart. Phys.* **2014**, *2014*, 009. [CrossRef]

200. Bellini, E.; Sawicki, I. Maximal freedom at minimum cost: linear large-scale structure in general modifications of gravity. *J. Cosmol. Astropart. Phys.* **2014**, *2014*, 050. [CrossRef]

201. Hu, B.; Raveri, M.; Frusciante, N.; Silvestri, A. Effective Field Theory of Cosmic Acceleration: an implementation in CAMB. *Phys. Rev. D* **2014**, *89*, 103530. [CrossRef]

202. Raveri, M.; Hu, B.; Frusciante, N.; Silvestri, A. Effective Field Theory of Cosmic Acceleration: constraining dark energy with CMB data. *Phys. Rev. D* **2014**, *90*, 043513. [CrossRef]

203. Zumalacárregui, M.; Bellini, E.; Sawicki, I.; Lesgourgues, J.; Ferreira, P.G. hi_class: Horndeski in the Cosmic Linear Anisotropy Solving System. *J. Cosmol. Astropart. Phys.* **2017**, *2017*, 019. [CrossRef]

204. Bellini, E.; Sawicki, I.; Zumalacárregui, M. hi_class: Background Evolution, Initial Conditions and Approximation Schemes. *J. Cosmol. Astropart. Phys.* **2020**, *2020*, 008. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.