Modeling and testing screening mechanisms in the laboratory and in space

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 focus on reviewing numerical and technical aspects involved in modeling the non-linear dynamics of screening. In this review, we focus on tests using laboratory experiments and astrophysical systems, such as stars, galaxies and dark matter halos.

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 focus on a broad range of astrophysical tests. We particularly review 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 newly emerging field of research which aims at producing screening maps of large-scale structure using constrained dark matter simulations. We conclude in Section 6.

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][12][13] and symmetron [14,15] mechanisms, two of the most studied scenarios in the literature, rely on such conformal couplings. More concretely, we consider a generic theory described by the action where M Pl is the Planck mass, R is the Ricci scalar corresponding to the Einstein-frame metric g µν , and S M is the action for the matter fields. The real-valued scalar field φ is assumed to be coupled to g µν , while matter fields Ψ couple to the Jordan frame metric g µν . The two metrics are related by the following conformal transformation where A(φ) is a generic function of the scalar field. The model is fully specified in terms of the coupling function A(φ) and the potential V(φ).
The scalar-field equation of motion is given by where T is the trace of the matter stress-energy tensor in the Jordan frame and we have introduced the effective potential V eff (φ) 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, ρ ≡ A 3 ρ can be shown to be conserved with respect to the Einstein-frame metric. As a result, the effective potential takes the form Due to the non-trivial coupling, test particles feel an additional "fifth" force per unit mass given by The qualitative realization of screening varies depending on the form of the coupling function A(φ). Particularly, the sum of runaway potentials and coupling functions in chameleon mechanism leads to environmentally-dependent effective mass for the scalar field. The larger mass in higher-density regions leads to shortened range of interactions. In contrast, the cymmetron 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.
The models described above and their variations have highly non-linear equations of motion. While we will focus 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 ρ(r) = ρΘ(R − r), with Θ 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 where primes denote derivative with respect to the radial coordinate r. We are interested in finding the field profile satisfying the following boundary conditions where the vacuum expectation value in a ρ = 0 environment is given by φ 0 = µ/ √ λ. First, consider a sphere with density higher than the spontaneous symmetry breaking density ρ SSB ≡ µ 2 M 2 . The field profile can be obtained analytically by approximating the system as a free scalar field with a potential centered at φ = 0. This leads to where ζ ≡ ρ µ 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 φ = φ 0 , and the field profile is given by 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 A limit of interest is the very strongly perturbing source with ρ ≫ µ 2 M 2 . In this limit for the coefficient B we have where the dimensionless parameter α ≡ ρR 2 /M 2 controls the characteristic outer shell thickness within which the field profile is varied significantly, ∆R/R = 1/α. It is instructive to consider large and small α limits. The coefficient B takes the form The force between the sphere of mass m sphere and a test particle with mass m test is then given by where in the last result we have introduced the screening factor λ sphere ≡ 3/α. 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, α ≫ 1, the fifth force on test particles is substantially suppressed due to the small screening factor λ sphere ≪ 1. It is also instructive to compare the scalar force with the Newtonian gravitational interaction: The screening factor representation will be useful when discussing experimental configurations involving atoms and other test objects.

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 behaviour, and thus the results described here can be transferred to such models.
The canonical example of a Vainshtein-screened theory is the galileon [17]: scalartensor theories where the scalar, φ, obeys the galileon symmetry φ → φ + a + b µ x µ , where a and b µ are constants. The cubic galileon has the action for constants c 2 and c 3 and where M 3 = M Pl H 2 0 . 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 Eq. (17) admit ghost-free self-accelerating de Sitter solutions [17], and conditions on model parameters exist which 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 quasistatic approximation to obtain the equation of motion [33] where i ∈ {1, 2, 3} and Here we are considering perturbations, φ, about a background,φ. It is common to rewrite Eq. (18) in terms of the cross-over scale, r C ≡ β 1 a 2 M 3 −1 and coupling α ≡ M Pl /(3β 2 ), instead of β 1 and β 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 for where T µν is the stress-energy tensor of the test object [34]. For a non-relativistic object of mass m, Q ≈ m and thus the coupling to the field is the same as for gravity. In general, however, Q < m since T µ µ 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 Eq. (18) as a total derivative, and thus derive the radial variation of φ to be [37] dφ where M(r) is the mass (perturbation) enclosed at radius r, and the Vainshtein radius is defined to be which is itself a function of r. In this case, we see that well within the Vainshtein radius (r ≪ 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 ∆G/G N = 2α 2 . Similar behaviour 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 ∼ O(100 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 Eq. (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 ≳ 0.1R 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 this object lives in an external environment. The galileon symmetry means that, if we have a solution for the galileon field φ, we can obtain another solution with ∂ µ φ → ∂ µ φ + b µ , for constant b µ . 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 φ has linear dynamics on ≳ 10 Mpc if r C ≃ 6 Gpc [40][41][42][43].

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. Eq. (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][45][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 ∼ λ 0 /c, with λ 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 λ 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 χ ≡ φ/φ 0 and consider the following equation: where the vacuum Compton wavelength for symmetrons is given by We impose the standard boundary conditions ∂χ/∂r = 0 at both r = 0 and r → ∞. 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. [50], without encountering significant performance differences in most of the problems of interest. The discretized equation can be written as Here we have introduced the following dicretizations of the Laplace operator and the potential term In order to obtain a solution to the non-linear discretized problem an iterative method should be employed. The particular structure of the operator 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 χ. At each iteration n we use the current estimate χ n (i) in order to obtain an improved estimate χ new (i) according to In most of the practical problems, this iteration itself does not converge. Instead, one should typically use a fraction of χ new (i) as a current estimate: where 0 < ω ⩽ 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 χ 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 and the mesh-averaged change of the field profile between consecutive iterations One could, for example, choose R 2 to reach a multiple of the numerical precision for the given implementation. In order to validate the performance of the solver, and asses 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 Eq. (24) and reconstruct a unique density profile that serves as a source for the mentioned field profile. As an example, we can choose χ ∼ tanh(r) and solve for the density configuration using Eq. (24). The gray line in the left panel of Fig. 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. 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 ∆z and the plate surfaces being placed at −∆z/2 and +∆z/2. The field equation of motion in this setup is given by where we have additionally definedẑ ≡ z/ √ 2λ 0 andρ(ẑ) ≡ ρ(ẑ)/M 2 µ 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, ∆ẑ/2) and (∆ẑ/2, ∞) we can show that the value of the field on the plate surface χ 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,ẑ) withẑ < ∆ẑ/2 we obtain where F is the elliptic integral of the first kind, and χ g is the field value in the middle of the gap. Fixingẑ to ∆ẑ/2 and setting χ = 0 as one of our boundary conditions we can numerically solve for χ g . Having the latter we will then have χ as a function ofẑ in the gap expressed in terms of the Jacobi elliptic function. The numerical relaxation solution of Eq. (33) in the gap subject to boundary conditions χ(−∆ẑ/2) = 0 = χ(+∆ẑ/2) is depicted in Fig. 3 (red dots). This is in sub-percent-level agreement with the exact analytical solution obtained from Eq. (34).
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 in 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.

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 is 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 Ω with boundary ∂Ω 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): where Ω is an integral over the domain Ω. The test function v(x) should satisfy certain properties, and we will particularly assume it vanishes on the entire boundary ∂Ω. Next, by employing the Green's theorem, the derivative order can be reduced, resulting in the following "weak-form" equation: 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 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 ) = δ ij . Inserting basis expansions into Eq. (38) we arrive at This system can be conveniently rewritten in a matrix form as where 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-handside of Eq. (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 new ) is replaced by its linear expansion around u k : Substituting this in Eq. (41) and using Eq. (39) again we obtain where the matrix B k and vector C k are defined as Having solved for u new , the new estimate for the field profile is obtained as in Eq. (30). As a concrete example let us derive the corresponding linearized system for the powerlaw chameleon model described earlier. The equation of interest in dimensionless variables reads as where α here is a dimensionless quantity determined in terms of the chameleon parameters, χ is a field variable normalized by the value minimizing the effective potential in the ambient space, andρ 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 Following Eq. (43), we expand the non-linear term around the k-th estimate of the solution, yielding In this case, the matrix B k and vector C k in Eq. (44) take the form 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 L u, v j = 0, we linearize the left-hand side around the current estimate for the solution, but now solve for the increment where M and B are defined as in the Picard method, while D reads as For the chameleon example discussed earlier, the vector B is given by Eq. (50), and the quantity D coincides with C in Eq. (51) but without the factor of n + 2 in the first term. Once the solution for δu is found, the field estimate for the next iteration is obtained as Here 0 < ω ⩽ 1 is a hyperparameter analogous to the one in Eq. (30), and controls the rate of iterative convergence. In both Eqs. (44) and (52) the matrix M should only be computed once before the iterative evaluation starts. The matrix M + 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 iteration 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 Ax = b can be derived by treating the problem as a minimization problem for the function Such a problem can be efficiently solved with a conjugate gradient algorithm, which relies on starting with an initial guess x 0 , and updating the solution by minimizing the function f along mutually conjugate directions The latter construction ensures that minimizing solution along a direction 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 behaviour 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.

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 µν into the matter contribution T matter µν and scalar field contribution T φ µν . 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: The force on the object can be computed from the time derivative of the 3-momentum. The latter can be written as:Ṗ 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 σ 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 Fig. 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 Eq. (34). Given the specific choice of the enclosing surface, the matter of the plate does not contribute to the integral in Eq. (56). For the Symmetron field, we have From the symmetry of the problem, the only component of the force is in x-direction, and it reads asṖ 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: where a useful identity for Jacoby elliptic functions has been employed 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 Eq. (59).

Laboratory tests
Laboratory experiments that are sensitive to the effects of weak forces are capable of providing useful constraints on the screening mechanisms. We will schematically divide such experiments into two broad categories. The experiments in the first category constrain 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][67][68][69][70][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][76][77][78][79] and torsion balance experiments [80]. We make this distinction primarily because of the unique technical challenges associated 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 modeling. 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. Fig. 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].

Direct force measurements
Casimir Experiments: The primary purpose of Casimir-force experiments is the measurement of the force due to modified vacuum energy density in presence of conducting extended bodies (see [82][83][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][85][86][87][88] (see, however, [89,90] for modern experimental configurations with planar geometry). The schematic setup of interest is depicted in the right panel of Fig. 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 ≪ µ −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 where λ sphere is the screening factor introduced in Subsection 2.1, and φ is the analytically solvable field profile of the isolated plate (see the end of Subsection 3.1). As a result we can obtain an accurate analytical force estimate when R ≪ µ −1 .
The situation is more subtle in the opposite limit of a very large sphere, R ≫ µ −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 where L = D + R + R cos θ is the x-distance between the plate and a ring (see Fig.4 for the notation). We expect both the approximate approaches to fail in the intermediate regime in between of R ≫ µ −1 and R ≪ µ −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 pen-dulum 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.

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 λ atom and λ source of the atom and the source, and for Chameleon theories can be written as This expression includes both the gravitational acceleration a grav and the acceleration a φ due to fifth forces. A central feature of such configurations is that, while the source could be screened and have λ source ≪ 1, the atoms are typically not, λ 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.

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 utilised 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] for pressure P, density ρ and where M is the mass enclosed at radius r. In standard gravity Υ = 0, and is more generally given in terms of the α B , α T and α H parameters [104] as The mere existence of stars already requires Υ > −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ρ n+1 n , for positive constants K and n. One can then define dimensionless variables in terms of the central density ρ c and pressure P c to be and thus arrive at the modified Lane-Emden equation [102] 1 which has boundary conditions θ(0) = 1, θ ′ (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 Υ 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 Υ to be ≲ 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.08M ⊙ [107]), lead to similar constraints, ruling out Υ ≳ 1.6. Modifications of gravity can also change the radius of brown dwarfs from their GR prediction (0.1R ⊙ ) to arbitrarily large values if gravity is weaker, or down to 0.078R ⊙ 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][109][110][111][112][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 Υ, and have derived an approximate constraint of −10 −3 < Υ < 5 × 10 −4 at the 2σ confidence level.

Out of equilibrium -stellar oscillations
Until now we have only considered stars in equilibrium. Moving beyond this, suppose one can consider small radius perturbations δr in the star which obey the equation In the absence of modified gravity, the period of small oscillations is found to be . (69) Hence any enhancement in the strength of gravity (G → G + ∆G with ∆G > 0) would reduce the period [117]. The change in the period can be as large as 30% for 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 become 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ögås 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 δr = rχ(r) exp(iωt) and linearise the pressure and potential to obtain a modified Linear Adiabatic Wave Equation [124] d dr where and a subscript '0' denotes an unperturbed quantity (see Sakstein [118] for the chameleonscreened version of this equation). For stars in GR, this results in unstable stars if Γ < 4/3, whereas in beyond Horndeski theory there is a second instability for sufficiently large Υ, 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 Υ.

Screening maps
When determining whether a given object is screened or unscreened, one often computes properties of the local gravitational field: the potential, Φ, 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, ϕ, in the presence of a density perturbation, δρ, is modified in f (R) gravity to be where δR is perturbation to Ricci scalar and a is the scale factor. By integrating a 2 δρ 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, ∆ m ≡ M D /M L − 1 [125,126]. In this case, ∆ m ∈ [0, 1/3], where ∆ 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 compute the external contribution to Φ to be where λ 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][130][131][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], Φ and a were computed using all objects within λ 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][141][142][143][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 constrained 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 package [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 obtain estimates for the initial conditions at z = 49. They then run 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 realisations, hence they find that it is better to generate constrained f (R) simulations from ΛCDM initial conditions than the ones they infer. Future work should be dedicated to producing fully-consistent constrained simulations and moving away from the linear approximation required in this reconstruction [150]. Phenomenological 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].

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, Φ c , such that objects with potentials below this value are unscreened and experience a fifth-force, and all other objects are screened. If Φ c ∼ 10 −6 , 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 centre of mass of the stars and the remaining mass of the galaxy. This offset will be stabilised 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 utilising the screening maps of Cabré et al. [127] (see Subsection 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][163][164], with the screening maps of Desmond et al. [134] being used to determine which objects are unscreened 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 −8 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 disfavoured. The modeling of other astrophysical effects for these constraints has been shown to be robust when calibrated against hydrodynamical ΛCDM simulations [166].  Figure 6. Constraints on the galileon coupling parameter, α, as a function of cross-over scale, r C , from astrophysical tests. The shaded regions are excluded at 1σ 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].

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 centre lags behind as it does not experience the fifth force. This offset is stabilised by the induced gravitational potential, and thus one would expect to see an O(kpc) offset between the black hole near a galaxy's centre and its centroid, as measured by the stars.
This test was first performed by Sakstein et al. [38] by considering the black holegalaxy offset in M87, utilising the analysis of Asvathaman et al. [167]. Here the galileon profile was determined analytically assuming spherical symmetry (Eq. (22)), which was shown to be an excellent approximation for cosmological simulations (see e.g. Ref. [41]). As shown in Fig. (6), α ∼ O(1) was found to be excluded for cross-over scales r C ≲ 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, utilising 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 α were found to be tighter than those obtained from M87, with ∆G/G N > 0.16 ruled out at 1σ confidence (see Fig. (6)), although these constraints are applicable to larger r C . We note that the constrained simulations were constructed assuming ΛCDM, and thus these constraints are only applicable for models which give rise to similar large scale structure as Λ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.

Halo properties
In GR, gravity determines the formation and growth of dark matter halos. Their abundance and internal structure can be modified in 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.
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 Subsection 2.2, even though the inner regions of NFW halos in Vainstein-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 variation 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: where ρ 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 ≡ r 200 /r s . The concentration and the characteristic density ρ s are connected via , 4πr 3 s ρ s = where ρ c is the critical density of the Universe, and g(c) ≡ 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.

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][186][187][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 have used both semi-analytical and N-body approaches to predict the changes in splashback location. For symmetrons, splashback has been investigated in Contigiani, Vardanyan, and Silvestri [50] by extending the self-similar collapse approach [192,193]. The latter approach, while not quantitatively very accurate, allows gaining 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: 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) ∼ a(t) s , with accretion rate s, and in EdS is a power-law in cosmic time t.
In absence of fifth force the shell variables can be rescaled such that all the shells satisfy the same equation of motion where we have labeled each shell by its turn-around time t * and radius r * , and have introduced λ ≡ r/r * , τ ≡ t/t * , and Λ = R(t)/r * . This system, complemented by initial conditions λ(τ) are λ(τ = 1) = 1, dλ/dτ(τ = 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 ρ 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 Λ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 Eq. (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 ρ c (z ssb ) = ρ ssb = M 2 µ 2 . Fig. 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 ≪ 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 behaviour there exists an optimal value of z ssb which maximizes the change in splashback location due to fifth forces.

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 nontrivially to gravity and matter fields, they exert additional, "fifth" forces on matter particles. The presence of fifth forces leads to 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][198][199][200]. Such unified approaches have been efficiently implemented in Einstein-Boltzmann numerical solvers [201][202][203][204].
On smaller scales the dynamics of scalar fields exhibit highly non-linear behaviour, 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, their precise phenomenological descriptions are given at 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 Subsection 4.1 can provide accurate description for the system for symmetrons, while fail for chameleons.
Numerical approaches are typically required in order to asses 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 for developing unified solver packages. The Finite-Elements-Based solvers introduced in Subsection 3.2, particularly, posses very favourable 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 Subsection 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 insights about the nature of gravitational interactions. At the same time, 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.  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.