Testing Screened Modified Gravity

Long range scalar fields with a coupling to matter appear to violate known bounds on gravitation in the solar system and the laboratory. This is evaded thanks to screening mechanisms. In this short review, we shall present the various screening mechanisms from an effective field theory point of view. We then investigate how they can and will be tested in the laboratory and on astrophysical and cosmological scales.

These non-linearities lead to the screening mechanisms that we review here. We do so irrespective of the origin and phenomenology of these scalar fields, be it dark matter-or dark energy-related, and present the screening mechanisms as a natural consequence of the use of effective field theory methods to describe the dynamics of scalar fields at all scales in the Universe. Given the ubiquity of new light degrees of freedom in modified gravity models-and the empirical necessity for screening-screened scalars represent one of the most promising avenues for physics beyond ΛCDM.
There are a number of excellent existing reviews on screening and modified gravity [36][37][38][39][40][41][42][43][44][45]. In [36] the emphasis is mostly on chameleons, in particular the inverse power law model, and symmetrons. K-mouflage is reviewed in [37] together with Galileons as an example of models characterised by the Vainshtein screening mechanism. A very comprehensive review on screening and modified gravity can be found in [39] where the screening mechanisms are classified into non-derivative and derivative, up to second order, mechanisms for the first time. There are subsequent more specialised reviews such as [42] on the chameleon mechanism, and [43] with an emphasis on laboratory tests. Astrophysical applications and consequences are thoroughly reviewed in [40,41] whilst more theoretical issues related to the construction of scalar-tensor theories of the degenerate type (DHOST) are presented in [45]. Finally a whole book [44] is dedicated to various approaches to modified gravity. In this review, we present the various screening mechanisms in a synthetic way based on an effective field theory approach. We then review and update results on the main probes of screening from the laboratory to astrophysics and then cosmology with future experiments in mind. Some topics covered here have not been reviewed before and range from neutron quantum bouncers to a comparison between matter spectra of Brans-Dicke and K-mouflage models.
We begin with a theoretical overview (Sec. II) before discussing tests of screening on laboratory (Sec. III), astrophysical (Sec. IV) and cosmological (Sec. V) scales. Sec. VI summarises and discusses the complementarity (and differences) between these classes of test.

II. SCREENING LIGHT SCALARS
A. Coupling scalars to matter Screening is most easily described using perturbations around a background configuration. The background could be the cosmology of the Universe on large scales, or the solar system. The perturbation is provided by a matter over density. This could be a planet in the solar system, a test mass in the laboratory or matter fluctuations on cosmological scales. We will simplify the analysis by postulating only a single scalar φ, although the analysis is straightforwardly generalised to multiple scalars. The scalar's background configuration is denoted by φ 0 and the induced perturbation of the scalar field due to the perturbation by an over density will be denoted by ϕ ≡ φ − φ 0 . At lowest order in the perturbation and considering a setting where space-time is locally flat (i.e. Minkowski), the Lagrangian describing the dynamics of the scalar field coupled to matter is simply [37,46] at the second order in the scalar perturbation ϕ and the matter perturbation δT µν . The latter is the perturbed energymomentum tensor of matter compared to the background. In this Lagrangian, matter is minimally coupled to the perturbed Jordan metric δg µν and the Jordan frame energy-momentum tensor is therefore conserved ∂ µ δT µν = 0. The expansion of the Lagrangian starts at second order as the background satisfies the equations of motion of the system. Notice that we restrict ourselves to situations where the Lorentz invariance is preserved locally. For instance in the laboratory we assume that the time variation of the background field is much slower than the ones of experiments performed on Earth. There are three crucial ingredients in this Lagrangian. The first is m(φ 0 ), i.e. the mass of the scalar field. The second is the wave function normalisation Z(φ 0 ). The third is the composite metric δg µν , which is not the local metric of space-time but the leading 2-tensor mediating the interactions between the scalar field and matter. This composite metric can be expanded as At leading order, the first term is the dominant one and corresponds to a conformal coupling of the scalar to matter with the dimensionless coupling constant β(φ 0 ) 1 . One can also introduce a term in second derivatives of ϕ which depends on a dimensionful coupling of dimension minus three. Finally going to higher order, there are also terms proportional to the first order derivatives of ϕ squared and a coupling constant of dimension minus four. These two terms can be seen as disformal interactions [47]. The equations of motion for ϕ are given by where δT ≡ δT µ ν and we have used the conservation of matter ∂ µ δT µν = 0. This equation will allow us to describe the different screening mechanisms.

B. Modified gravity
Let us now specialise the Klein-Gordon equation to experimental or observational cases where δT 00 = ρ is a static matter over density locally and the background is static too. This corresponds to a typical experimental situation where over densities are the test masses of a gravitational experiment. In this case, we can focus on the case where φ 0 can be considered to be locally constant. As a result we have The kinetic terms are modified by the tensor When the over densities are static, the disformal term in K µν which depends on the matter energy momentum tensor does not contribute and we have K ij (φ 0 ) Z(φ 0 )δ ij leading to the modified Yukawa equation whereβ(φ 0 ) = β(φ0) Z(φ0) . For nearly massless fields we can neglect the mass term within the Compton wavelength of size m −1 (φ 0 ) which is assumed to be much larger than the experimental setting. In this case the Yukawa equation becomes a Poisson equation ∆ϕ =β (φ 0 ) m Pl ρ.
As a result the scalar field behaves like the Newtonian potential and matter interacts with the effective Newtonian potential 2 i.e. gravity is modified with an increase of Newton's constant by Notice that the scalar field does not couple to photons as δT = 0, hence matter particles are deviated with a larger Newtonian interaction than photons, G eff ≥ G N . As a result, the modification of G N into G eff is not just a global rescaling and gravitational physics in genuinely modified. This appears for instance in the Shapiro effect (the time delay of photons in the presence of a massive object) as measured by the Cassini probe around the Sun. When the mass of the scalar field cannot be neglected, the effective Newton constant becomes distance-dependent: where r is the distance to the object sourcing the field. This equation allows us to classify the screening mechanisms.
C. The non-derivative screening mechanisms: chameleon and Damour-Polyakov The first mechanism corresponds to an environment-dependent mass m(φ 0 ). If the mass increases sharply inside dense matter, the scalar field emitted by any mass element deep inside a compact object is strongly Yukawa suppressed by the exponential term e −m φ (φ0)r , where r is the distance from the mass element. This implies that only a thin shell of mass ∆M at the surface of the object sources a scalar for surrounding objects to interact with. As a result the coupling of the scalar field to this dense object becomes where M is the mass of the object. As long as ∆M/M 1, the effects of the scalar field are suppressed. This is the chameleon mechanism [48][49][50][51].
The second mechanism appears almost tautological. If in dense matter the coupling β(φ 0 ) = 0, all small matter elements deep inside a dense object will not couple to the scalar field. As a result and similarly to the chameleon mechanism, only a thin shell over which the scalar profile varies at the surface of the objects interacts with other compact bodies. Hence the scalar force is also heavily suppressed. This is the Damour-Polyakov mechanism [52].
In fact this classification can be systematised and rendered more quantitative using the effective field theory approach that we have advocated. Using Eq. (7), we get Let us first consider the case of a normalised scalar field with Z(φ 0 ) = 1. The scalar field is screened when its response to the presence of an over density is suppressed compared to the Newtonian case. This requires that |ϕ| 2m Pl |Φ N | ≤ β(φ out ) (13) where ϕ = φ in −φ 0 is the variation of the scalar field inside the dense object. Here Φ N is the Newtonian potential at the surface of the object. This is the quantitative criterion for the chameleon and Damour-Polyakov mechanisms [48,53]. In particular, in objects which are sufficiently dense, the field φ in nearly vanishes and ϕ −φ 0 only depends on the environment. As a result for such dense objects, screening occurs when |Φ N | ≥ φ0 2m Pl β(φout) which depends only on the environment. Chameleon and Damour-Polyakov screenings occur for objects with a large enough surface Newtonian potential. In fact it turns out that for a screened object labelled by A is the scalar charge of this object 3 , i.e. its coupling to matter. The screening criterion (13) simply requires that the scalar charge of an object is less than the coupling of a test particle β(φ 0 ). 3 One can also introduce the screening factor λ A = β A β(φ 0 ) whereby screening occurs when λ A ≤ 1. The screening factor is also related to the mass of the thin shell ∆M A as ∆M A The third case in fact covers two mechanisms. If locally in a region of space the normalisation factor Z(φ 0 ) 1 (15) then obviously the effective couplingβ(φ 0 ) 1 and gravitational tests can be evaded. Notice that we define screening as reducing the effective coupling. This case covers the K-mouflage 4 and Vainshtein mechanisms.
The normalisation factor is a constant at leading order. Going beyond leading order, i.e. including higher order operators in the effective field theory, Z(φ 0 ) can be expanded in a power series where r c is a cross over scale and has the dimension of length and M is an energy scale. The scale Λ plays the role of the strong coupling scale of the models. The functions a, b and c are assumed to be smooth and of order unity.

K-mouflage
The K-mouflage screening mechanism [46,54,55] is at play when Z(φ 0 ) ≥ 1 and the term in (∂ϕ) 2 /Λ 4 dominates in (16), i.e. | ∇ϕ| ≥ Λ 2 (17) and therefore the Newtonian potential must satisfy Hence K-mouflage screening occurs where the gravitational acceleration a N = − ∇Φ N is large enough. Let us consider two typical situations. First the Newtonian potential of a point-like source of mass M has a gradient satisfying (18) inside a radius R K The scalar field is screened inside the K-mouflage radius R K . Another interesting example is given by the large scale structures of the Universe where the Newtonian potential is sources by over densities δρ compared to the background energy densityρ. In this case, screening takes place for wave-numbers k such that where δ ≡ δρ ρ . In particular for models motivated by dark energy Λ 4 m 2 Pl H 2 0 screening occurs on scales such that k/H 0 > ∼ β(φ 0 )δ, i.e. large scale structures such as galaxy clusters are not screened as they satisfy k/H 0 < ∼ 1 [56,57].

Vainshtein
The Vainshtein mechanism [58,59] follows the same pattern as K-mouflage. The main difference is that now the dominant term in Z(φ 0 ), i.e. (16), is given by the ϕ term. This implies that i.e. screening occurs in regions where the spatial curvature is large enough. Taking once again a point source of mass M , the Vainshtein mechanism is at play on scales smaller than the Vainshtein radius. 5 Notice the power 1/3 compared to the 1/2 in the K-mouflage case. Similarly on large scales where the density contrast is δ, the scalar field is screened for wave numbers such that where Ω m is the matter fraction of the Universe when r c = 1/H 0 . The Vainshtein mechanism is stronger than K-mouflage and screens all structures reaching the non-linear regime δ > ∼ 1 as long as β(φ 0 ) > ∼ 1. Finally let us consider the case where the term in 2 ϕ dominates in (16). This corresponds to 6 For a point source the transition happens at the radius As expected the power is now 1/5 which can be obtained by power counting. This case is particularly relevant as this corresponds to massive gravity and the original investigation by Vainshtein. In the massive gravity case [60,61] where m G is the graviton mass. In all these cases, screening occurs in a regime where one would expect the effective field theory to fail, i.e. when certain higher order operators start dominating. Contrary to expectation, this is not always beyond the effective field theory regime. Indeed, scalar field theories with derivative interactions satisfy non-renormalisation theorems which guarantee that these higher order terms are not corrected by quantum effects [62,63]. Hence the classical regime where some higher order operators dominate can be trusted. This is in general not the case for non-derivative interaction potentials, which are corrected by quantum effects. As a result, the K-mouflage and Vainshtein mechanisms appear more robust that the chameleon and Damour-Polyakov ones under radiative effects. E. Screening criteria: the Newtonian potential and its derivatives Finally let us notice that the screening mechanisms can be classified by inequalities of the type where C is a dimensionful constant and k = 0 for chameleons, k = 1 for K-mouflage and k = 2 for the Vainshtein mechanism. This implies that it is the Newtonian potential, acceleration and spacetime curvature respectively that govern objects' degrees of screening in these models. The case k = 4 appears for massive gravity. Of course, if higher order terms in the expansion of Z(φ 0 ) in powers of derivatives were dominant, larger values of k could also be relevant. As we have seen, from an effective field theory point of view, the powers k = 0, 1, 2 are the only ones to be considered. The case of massive gravity k = 4 only matters as the other cases k ≤ 2 are forbidden due to the diffeomorphism invariance of the theory, see the discussion in section II G 1. 5 The equation (21) should be understood as integrated over a ball of radius r. The left hand side is proportional to the point mass and the right hand side to the volume of the ball 6 This inequality can be understood as ∆Φ N ≥ where the integration volume is taken as a ball of radius r and ∆ −1 (r) = − 1 4πr .

F. Disformally induced charge
Let us now come back to a situation where the time dependence of the background is crucial. For future observational purposes, black holes are particularly important as the waves emitted during their collisions could carry much information about fundamental physics in previously untested regimes. For scalar fields mediating new interactions, this seems to be a perfect new playground. In most scalar field theories, no-hair theorems prevent the existence of a coupling between black holes and a scalar field, implying that black holes have no scalar charge (see Sec. IV B for observational consequences of this). However, these theorems are only valid in static configurations; in a time-dependent background the black hole can be surrounded by a non-trivial scalar cloud.
Let us consider a canonically normalised and massless scalar field in a cosmological background. As before we assume that locally Lorentz invariance is respected on the time scales under investigation. The Klein-Gordon equation becomes in the presence of a static overdensity ∆ϕ =γ(φ 0 )ρ.
As a result, we see that a scalar charge is induced by the cosmological evolution of the background [64,65] where γ 1 = dγ dφ and γ 2 = d 2 γ dφ 2 . This is particularly relevant to black holes solutions with a linear dependence in timė φ 0 = q. In this case the induced charge is strictly constant which could lead to interesting phenomena in binary systems.
G. Examples of screened models

Massive gravity
The first description of screening in a gravitational context was given by Vainshtein and can be easily described using the Fierz-Pauli [60] modification of General Relativity (GR). In GR and in the presence of matter represented by the energy-momentum tensor T µν , the response of the weak gravitational field h µν = g µν − η µν is given in momentum space by 7 where two features are important. The first is that 1/p 2 = 1/p λ p λ is characteristic of the propagation of a massless graviton. The second is the 1/2 factor which follows from the existence of two propagating modes. When the graviton becomes massive, the following mass term is added The tensorial structure of the Fierz-Pauli mass term guarantees the absence of ghosts in a flat background. The response to matter becomes The factor in 1/(p 2 + m 2 G ) is the propagator of a massive field of mass m G . More surprising is the change 1/2 → 1/3 in the tensorial expression. In particular, in the limit m G → 0 one does not recover the massless case of GR. This is the famous vDVZ (van Dam-Veltman-Zakharov) discontinuity [66,67]. Its origin can be unraveled as follows. Writing 7 As the background metric is the Minkowskian one, the use of Fourier modes is legitimate. and corresponding to a scalar satisfying we find that (33) is satisfied provided that Hence we have decomposed the massive graviton into a helicity two parth µν and a scalar part ϕ coupled to matter with a scalar charge β = 1/ √ 6. These are three of the five polarisations of a massive graviton. Notice that the scalar polarisation is always present however small the mass m G , i.e. the massless limit is discontinuous as the number of propagating degrees of freedom is not continuous. As it stands, massive gravity with such a large coupling and a mass experimentally constrained to be m G ≤ 10 −22 eV would be excluded by solar system tests. This is not the case thanks to the Vainshtein mechanism.
Indeed non-linear interactions must be included as GR is not a linear theory. At the next order one expects terms in h 3 leading to Lagrangian interactions of the type [61] where Λ 5 = m Pl m 4 G . The structure in ϕ follows from the symmetry ϕ → ϕ + λ µ x µ which can be absorbed in to . The Klein-Gordon equation is modified by terms in (( ϕ) 2 ). As a result, the normalisation factor is dominated by Z(φ 0 ) ∼ 2 ϕ/Λ 5 as mentioned in the previous section. This leads to the Vainshtein mechanism inside R MG which allows massive gravity to evade gravitational tests in the solar system for instance.

Cubic Galileon models
The cubic Galileon models [59] provide an example of Vainshtein mechanism with the 1/5 power instead of the 1/3. They are defined by the Lagrangian The normalisation factor for the kinetic terms involves φ as expected. These theories are amongst the very few Horndeski models which do not lead to gravitational waves with a speed differing from the speed of light. Unfortunately as theories of self-accelerating dark energy, i.e. models where the acceleration is not due to a cosmological constant, they suffer from an anomalously large Integrated-Sachs-Wolfe (ISW) effect in the Cosmic Microwave Background (CMB). See section II H for more details.

Quartic K-Mouflage
The simplest example of K-mouflage model is provided by the Lagrangian [64] which is associated to a normalisation factor containing a term in (∂φ) 2 . These models pass the standard tests of gravity in the solar system but need to be modified to account for the very small periastron anomaly of the Moon orbiting around the Earth. See section II I for more details.

Ratra-Peebles and f(R) chameleons
Chameleons belong to a type of scalar-tensor theories [68] specified entirely by two function of the field. The first one is the interaction potential V (φ) and the second one is the coupling function A(φ). The dynamics are driven by the effective potential [48,49] V eff (φ) = V (φ) + (A(φ) − 1)ρ (42) where ρ is the conserved matter density. When the effective potential has a minimum φ 0 ≡ φ(ρ), its second derivative defines the mass of the chameleon Cosmologically the chameleon minimum of the effective potential is an attractor when m(ρ) H, i.e. the mass is greater than the Hubble rate [50]. This is usually guaranteed once the screening of the solar system has been taken into account, see section II I. A typical example of chameleon theory is provided by [48,49] associated to a constant coupling constant β. More generally the coupling becomes density dependent as Chameleons with n = 1 are extremely well constrained by laboratory experiments, see section III E. Surprisingly models of modified gravity defined by the Lagrangian [34] which is a function of the Ricci scalar R can be transformed into a scalar-tensor setting. First of all the the field equations of f (R) gravity can be obtained after a variation of the Lagrangian (46) with respect to the metric g µν and they read where f R ≡ df (R)/dR and is the d'Alembertian operator. These equations naturally reduce to Einstein's field equations when f (R) = R. This theory can be mapped to a scalar field theory via where β = 1 √ 6 . The coupling function is given by the exponential leading to the same coupling to matter as massive gravity. Contrary to the massive gravity case, f (R) models evade solar system tests of gravity thanks to the chameleon mechanism when the potential is appropriately chosen. A popular model has been proposed by Hu-Sawicki [69] and reads which involves two parameters, the exponent (positive definite) n > 0 and the normalisation f R0 which is constrained to be |f R0 | < ∼ 10 −6 by the requirement that the solar system is screened [69] (see Sec. II I). On large scales, structures are screened for which for the n = 1 Hu-Sawicki model, where χ is the "self-screening parameter". This follows directly from the fact that where δφ is the variation of the scalaron due to a structure in the present Universe. Assessing the inequality in (52) -or equivalently requiring that the scalar charge Q = |δφ| 2βm Pl Φ N must be less that β = 1/ √ 6 -gives a useful criterion for identifying unscreened objects (see Sec. IV).

f (R) and Brans-Dicke
The f (R) models can be written as a scalar-tensor theory of the Brans-Dicke type. The first step is to replace the f (R) Lagrangian density by which reduces to the original model by solving for λ. Then an auxiliary field ψ ≡ df (λ) dλ (55) can be introduced, together with the potential V (ψ) = m 2 Pl (f (λ(ψ)) − λ(ψ)ψ)/2, which corresponds to the Legendre transform of the function f (λ). After replacing back into the original action, one recovers a scalar field action for ψ in the Jordan frame that reads This theory corresponds to the well known Generalized Jordan-Fierz-Brans-Dicke [70] theories with ω BD = 0. When the ω BD parameter is non-vanishing and a constant, this reduces to the popular Jordan-Brans-Dicke theory. Exact solutions of these theory have been tested against observations of the Solar System [33,71] and the Cassini mission sets the constraint ω BD > 40, 000, so that JBD has to be very close to GR. This bound is a reformulation of (88), see Sec. II I for more details. After going to the Einstein frame the theory must a scalar-tensor with the Chameleon or Damour-Polyakov mechanisms in order to evade the gravitational tests in the solar system.

The symmetron
The symmetron [72] is another scalar-tensor theory with a Higgs-like potential and a non-linear coupling function where the quadratic term is meant to be small compared to unity. The coupling is given by which vanishes at the minimum of the effective potential when ρ > µ 2 M 2 . This realises the Damour-Polyakov mechanism.
7. Beyond 4d: Dvali-Gabadadze-Porrati gravity The Dvali-Gabadadze-Porrati (DGP) gravity model [73] is a popular theory of modified gravity that postulates the existence of an extra fifth-dimensional Minkowski space, in which a brane of 3+1 dimensions is embedded. Its solutions are known to have two branches, one which is self-accelerating (sDGP), but is plagued with ghost instabilities [74] and another branch, the so-called normal branch (nDGP) which is non-self-accelerating, and has better stability properties. At the nonlinear level, the fifth-force is screened by the effect of the Vainshtein mechanism and therefore can still pass solar system constraints. This model can be written as a pure scalar-field model and in the following we will use the notations of [75] to describe the model and its cosmology. The action is given by where L matter is the matter Lagrangian, R is the Ricci scalar built from the bulk metric g ab and M 4 and M 5 are the Planck scales on the brane and in thebulk, respectively. The metric g µν is on the brane, R its Ricci scalar, and K = g µν K µν is the trace of extrinsic curvature, K µν . Finally, σ is the tension or bare vacuum energy on the brane. These two different mass scales give rise to a characteristic scale that can be written as For scales larger than r 5 , the 5 dimensional physics contributes to the dynamics, while for scales smaller than r 5 , gravity is 4 dimensional and reduces to GR. The reader can find the complete set of field equations in [75]. After solving the Friedmann equations, the effective equation of state of this model is given by where κ is the 3-dimensional spatial curvature. During the self-accelerating phase w eff → −1 in (62), therefore emulating a cosmological constant.

H. Horndeski theory and beyond
For four dimensional scalar-tensor theories used so far, the action defining the system in the Einstein frame can be expressed as where φ is the scalar field, V (φ) its potential and it couples to the matter fields ψ (i) m through the Jordan frame metric g µν , which is related to the metric g µν as The disformal factor term in B 2 (φ, X) leads to the derivative interactions in (2). In the previous discussions, see Sec. II G 4, we focused on the conformal parameter A(φ) chosen to be X-independent where X = −(∂φ) 2 /2Λ 2 and Λ is a given scale. Other choices are possible which will dot be detailed here, in particular in the case of DHOST theories for which the dependence of A(φ, X) is crucial [45]. As can be expected, (63) can be generalized to account for all possible theories of a scalar field coupled to matter and the metric tensor. When only second order equations of motion are considered, this theory is called the Horndeski theory. Its action can be written as where the four Lagrangian terms corresponds to different combinations of 4 functions G 2,3,4,5 of the scalar field and its kinetic energy χ = −∂ µ ∂ µ φ/2, the Ricci scalar and the Einstein tensor G µν and are given by After the gravitational wave event GW170817 ( [76,77], and as already anticipated in [? ], the propagation of gravitational waves is practically equal to the speed of light, implying that a large part of Horndeski theory with cosmological effects, is ruled out, leaving mostly only models of type L 2 and Cubic Galileons (Horndeski with Lagrangians up to L 3 ) as the surviving class of models [78][79][80]. They are the ones that will be dealt with in this review and can be linked most directly to the screening mechanisms described here. When going beyond the Horndeski framework [81], the Vainshtein mechanism can break within massive sources [82]. This phenomenology was studied further in [83], and may be used to constrain such theories as described in Sec. IV A.

I. Solar system tests
Screening mechanisms have been primarily designed with solar system tests in mind. Indeed light scalar fields coupled to matter should naturally violate the gravitational tests in the solar system as long as the range of the scalar interaction, i.e. the fifth force, is large enough and the coupling to matter is strong enough. The first and most prominent of these tests is provided by the Cassini probe [33] and constrains the effective coupling between matter and the scalar to be as long as the range of the scalar force exceeds several astronomical units and β 2 eff corresponds to the strength of the fifth force acting on the satellite. As we have mentioned this translates into the effective bound where φ 0 is the value of the scalar field in the interplanetary medium of the solar system. Here we have assumed that the Cassini satellite is not screened and the Sun is screened. As a result, the scalar charges are respectively the background one β(φ 0 ) for the satellite and β for the Sun. In the case of the K-mouflage and Vainshtein mechanisms, the scalar charges of the Sun and the satellite are equal and the Cassini bound can be achieved thanks to a large Z(φ 0 ) factor. As an example, for cubic Galileon models the ratio between the fifth force and the Newtonian force behaves like where β(φ 0 ) = β and R V is the Vainshtein radius. For cosmological models where L ∼ H −1 0 , the Vainshtein radius of the Sun is around 0.1 kpc. As a result for the planets of the solar system r/r V 1 and the fifth force is negligible. K-mouflage models of cosmological interest with Λ Λ DE ∼ 10 −3 eV lead to the same type of phenomenology with a K-mouflage radius of the Sun larger than 1000 a.u. and therefore no fifth force effects in the solar system. For chameleon-like models the Cassini constraint becomes where we have assumed that Z(φ 0 ) = 1 and β(φ 0 ) 1. This is a stringent bound which translates into for the values of the scalar in the solar system. Indeed we have assumed that in dense bodies such as the Sun or planets, the scalar field vanishes. We have also used the Newtonian potential of the Sun Φ N 10 −6 .
In fact chameleon-screened theories are constrained even more strongly by the Lunar Ranging experiment (LLR) [84,85]. This experiment constrains the Eötvos parameter for two bodies falling in the presence of a third one C. The accelerations a A,B are towards C and due to C. For bodies such as the Earth A = ⊕, the moon B = moon and the Sun C = , a non-vanishing value of the Eötvos parameter would correspond to a violation of the strong equivalence principle, i.e. a violation of the equivalence principle for bodies with a non-negligible gravitational self-energy. Such a violation is inherent to chameleon-screened models. Indeed, screened bodies have a scalar charge β A which is dependent on the Newtonian potential of the body β A ∝ Φ −1 A implying a strong dependence on the nature of the objects. As the strength of the gravitational interaction between two screened bodies is given by as long as the two objects are closer than the background Compton wavelength m −1 (φ 0 ), the Eötvos parameter becomes In the case of the LLR experiment we have β A φ 0 /2β(φ 0 )Φ A m Pl and therefore β moon β ⊕ . Using Φ N moon 10 −11 and Φ ⊕ 10 −9 we find that the LLR constraint This becomes for the scalar charge of the Earth [49] β ⊕ < ∼ 10 −6 (76) which is stronger than the Cassini bound, i.e. we must impose that This corresponds to the energy scale of particle accelerators such as the Large Hadron Collider (LHC). This bound leads to relevant constraint on the parameter space of popular models. Let us first consider the n = 1 inverse power law chameleon model with This model combines the screening property of inverse power law chameleon and the cosmological constant term Λ 4 DE leading to the acceleration of the expansion of the Universe. The mass of the scalar is given by implying that in the solar system m 0 > ∼ 10 6 H 0 . Now as long as the chameleon sits at the minimum of its effective potential we have m cosmo ( ρcosmo ρ G ) 1/2 m 0 where ρ cosmo is the cosmological matter density and ρ G 10 6 ρ cosmo is the one in the Milky Way. As a results we have the constraints on the cosmological mass of the chameleon [86] [87] m cosmo > ∼ 10 3 H 0 .
As the Hubble rate is smaller than the cosmological mass, the minimum of the effective potential is a tracking solution for the cosmological evolution of the field. This bound (80) is generic for chameleon-screened models with an effect on the dynamics of the Universe on large scale. In the context of the Hu-Sawicki model, and as m cosmo /H 0 ∝ f −1/2 R0 , the solar system tests imply typically that f R0 < ∼ 10 −6 [69]. For models with the Damour-Polyakov mechanism such as the symmetron, and if ρ G ≤ µ 2 M 2 , the field value in the solar system is close to φ 0 µ √ λ . The mass of the scalar is also of order µ implying that the range of the symmetron is very short unless µ < ∼ 10 −18 eV. In this case the LLR bound applies and leads to which implies that the symmetron models must be effective field theory below the grand unification scale. Models with derivative screening mechanisms such as K-mouflage and Vainshtein do not violate the strong equivalence principle but lead to a variation of the periastron of objects such as the Moon [88]. Indeed, the interaction potential induced by a screening object does not vary as 1/r anymore. As a result Bertrand's theorem 8 is violated and the planetary trajectories are not closed anymore. For K-mouflage models defined by a Lagrangian L = Λ 4 K(X) where X = −(∂φ) 2 /2Λ 4 and Λ Λ DE , the periastron is given by [89] where x = r/R K is the reduced radius (R K is the K-mouflage radius). For the Moon, the LLR experiment implies that δθ ≤ 2 · 10 −11 which constrains the function K(X) and its derivatives K (X) and K (X). A typical example of models passing the solar system tests is given by K(X) = −1 + X + K (X − X arctan( X X )) with X 1 and K > ∼ 10 3 . In these models, the screening effect is obtained as K ∼ K 1 as long as |X| > ∼ |X |. For cubic Galileons, the constraint from the periastron of the Moon reduces to a bound on the suppression scale [88,90] The lower bound corresponds to Galileon models with an effect on cosmological scales. Finally, models with the K-mouflage or the Vainshtein screening properties have another important characteristic. In the Jordan frame where particles inside a body couple to gravity minimally, the Newton constant is affected by the conformal coupling function A(φ), i.e.
For chameleon-screened objects, the difference between the Jordan and Einstein values of the Newton constant is irrelevant as deep inside screened objects φ is constant and A(φ) can be normalised to be unity. This is what happens for symmetrons or inverse power law chameleons for instance. For models with derivative screening criteria, i.e. K-mouflage or Vainshtein, the local value of the field can be decomposed as where t 0 is the present time. Here φ( x) is the value of the field due to the local and static distribution of matter whilst the correction term depends on time and follows from the contamination of the local values of the field by the large scale and cosmological variations of the field. In short, regions screened by the K-mouflage or Vainshtein mechanisms are not shielded from the cosmological time evolution of matter. As a result, the Newton constant in the Jordan frame becomes time dependent with a drift [91] d ln where we have taken the scalar to be coupled conformally with a constant strength β. The LLR experiment has given a bound in the solar system [84] | i.e. Newton's constant must vary on timescales larger than the age of the Universe. This can be satisfied by Kmouflage or Vainshtein models provided β ≤ 0.1 as long asφ ∼ m Pl H 0 , i.e. the scalar field varies by an order of magnitude around the Planck scale in one Hubble time [89].

III. TESTING SCREENING IN THE LABORATORY
Light scalar fields have a long range and could induce new physical effects in laboratory experiments. We will consider some typical experiments which constrain screened models in a complementary way to the astrophysical and cosmological observations discussed below. In what follows, the bounds on the screened models will mostly follow from the classical interaction between matter and the scalar field. A light scalar field on short enough scales could lead to quantum effects. As a rule, if the mass of the scalar in the laboratory environment is smaller than the inverse size of the experiment, the scalar can considered to be massless. Quantum effects of the Casimir type imply that two metallic plates separated by a distance d will then interact and attract according to as long as the coupling between the scalar and matter is large enough 9 . In the Casimir or Eötwash context, this would mean that the usual quantum effects due to electromagnetism would be be increased by a factor of 3/2. Such a large effect is excluded and therefore the scalar field cannot be considered as massless on the scales of these typical experiments. In the following we will consider the case where the scalar is screened on the scales of the experiments, i.e. its typical mass is larger than the inverse of the size of the experimental set up. In this regime where quantum effects can be neglected, the classical effects induced by the scalars are due to the non-trivial scalar profile and its non-vanishing gradient 10 . In the following, we will mostly focus on the classical case and the resulting constraints.

A. Casimir interaction and Eötwash experiment
We now turn to the Casimir effect [93], associated with the classical field between two metallic plates separated by a distance d. The classical pressure due to the scalar field with a non-trivial profile between the plates is attractive and with a magnitude given by [? ] where A is the surface area of the plates and V eff is the effective potential. This is the difference between the potential energy in vacuum (i.e., without the plates) where the field takes the constant value φ 0 and in the vacuum chamber halfway between the plates. In general the field acquires a bubble-like profile between the plates and φ(0) is where the field is maximal. The density inside the plates is much larger than between the plates, so the field value inside the plates is zero to a very good approximation. For a massive scalar field of mass m with a coupling strength β, the resulting pressure between two plates separated by distance d is given by which makes explicit the Yukawa suppression of the interaction between the two plates. In the screened case the situation can be very different. Let us first focus on the symmetron case. As long as µ > ∼ d −1 , the value φ(0) is very close to the vacuum value φ(0) µ √ λ implying that F φ /A 0, i.e. the Casimir effect does not probe efficiently symmetrons with large masses compared to the inverse distance between the plates. On the other hand when µ < ∼ d −1 , the field essentially vanishes in the plates and between the plates [94]. As a result the classical pressure due to the scalar becomes Notice that is this regime, the symmetron decouples from matter inside the experiment as β(φ(0)) = 0. We will see how this compares to the quantum effects in section III F. We can now turn to the chameleon case where we assume that the density between the plates vanishes and is infinite in the plates. This simplified the expression of the pressure which becomes [95] where B(. , .) is Euler's B function. In the chameleon case, the pressure is a power law depending on 2n/(n + 2) which can be very flat in d and therefore dominates of the photon Casimir pressure at large distance. Quantum effects can also be taken into account when the chameleon's mass is small enough, see section III F. The most stringent experimental constraint on the intrinsic value of the Casimir pressure has been obtained with a distance d = 746 nm between two parallel plates and reads | ∆F φ A | ≤ 0.35 mPa [96]. The plate density is of the order of ρ plate = 10 g cm −3 . The constraints deduced from the Casimir experiment can be seen in section III E. It should be noted that realistic experiments sometimes employ a plate-and-sphere configuration, which can have an O(1) modification to (92) [97].
The Eöt-Wash experiment [98] is similar to a Casimir experiment and involves two rotating plates separated by a distance d. Each plate is drilled with holes of radii r h spaced regularly on a circle. The gravitational and scalar interactions vary in time as the two plates rotate, hence inducing a torque between the plates. This effect can be understood by evaluating the potential energy of the configuration. The potential energy is obtained by calculating the amount of work required to approach one plate from infinity [35,99]. Defining by A(θ) the surface area of the two plates which face each other at any given time, a good approximation to energy is simply the work of the force between the plates corresponding to the amount of surface area in common between the two plates. The torque is then obtained as the derivative of the potential energy of the configuration with respect to the rotation angle θ and is given by where a θ = dA dθ depends on the experiment and is a well-known quantity. As can be seen, the torque is a direct consequence of the classical pressure between two plates.
For a Yukawa interaction and upon using the previous expression (89) for the classical pressure, we find that the torque is given by which is exponentially suppressed with the separation between the two plates d. Let us now consider the symmetron and chameleon cases. In the symmetron case, the classical pressure is non-vanishing only when d < ∼ µ −1 implying that Hence the torque increases linearly before saturating at a maximal value. For chameleons, three cases must be distinguished. First when n > 2, the torque is insensitive to the long range behaviour of the chameleon field in the absence of the plates and we have which decreases with the distance. In the case n < 2, the torque is sensitive to the Yukawa suppression of the scalar field at distances larger that d ∼ m −1 0 , where m 0 is the mass in the vacuum between the plates. This becomes for d < ∼ d and essentially vanishes for larger distances. In the case n = 2, a logarithmic behaviour appears.
The 2006 Eöt-Wash experiment [100] gives the bound for a separation between the plates of d = 55µm is where Λ T = 0.35Λ DE [35] and Λ DE = 2.4 meV. We must also need to modify the torque calculated previously in order to take into account the effects of a thin electrostatic shielding sheet of width d s = 10µm between the plates in the Eöt-Wash experiment. This reduces the observed torque which becomes T obs = e −mcds T . As a result we have that Surprisingly, the Eötwash experiment tests the dark energy scale in the laboratory as Λ T ≈ Λ DE .

B. Quantum bouncer
Neutrons behave quantum mechanically in the terrestrial gravitational field. The quantised energy levels of the neutrons have been observed in Rabi oscillation experiments [101]. Typically a neutron is prepared in its ground state by selecting the width of its wave function using a cache, then a perturbation induced either mechanically or magnetically makes the neutron state jump from the ground state to one of its excited levels. Then the ground state is again selected by another cache. The missing neutrons are then compared with the probabilities of oscillations from the ground state to an excited level. This allows one to detect the first few excited states and measure their energy levels. Now if a new force complements gravity, the energy levels will be perturbed. Such perturbations have been investigated and typically the bounds are now at the 10 −14 eV level.
The wave function of the neutron satisfies the Schrödinger equation where m N is the neutron's mass and the potential over a horizontal plate is where φ(z) is the vertical profile of the scalar field. We put the mirror at z ≤ 0. The contribution due to the scalar field is which depends on the model. In the absence of any scalar field, the wavefunctions are Airy functions where c k is a normalisation constant, z 0 = (h 2 /2m 2 N g) 1/3 , − k are the zeros of the Airy function. Typically k = {2.338, 4.088, 5.521, 6.787, 7.944, 9.023 . . . } for the first levels k = 1 . . . . At the first order of perturbation theory, the energy levels are where E is the averaged value of the perturbed potential in the excited states.
Let us see what this entails for chameleon models [102,103]. In this case the perturbation depends on where the profile of the chameleon over the plate is given by Using this form of the correction to the potential energy, i.e. power laws, and the fact that the corrections to the energy levels are linear in β, one can deduce useful constraints on the parameters of the model. So far we have assumed that the neutrons are not screened. When they are screened, the correction to the energy levels are easily obtained by replacing β → λβ where λ is the corresponding screening factor. In the case of symmetrons, the correction to the potential energy depends on whilst the symmetron profile is given by [104] where we assume that the plate is completely screened. The averaged values of δV are constrained by which leads to strong constraints on symmetron models. See section III E.

C. Atomic interferometry
Atomic interferometry experiments are capable of very precisely measuring the acceleration of an atom in free fall [105,106]. By placing a source mass in the vicinity of the atom and performing several measurements with the source mass in different positions, the force F source between the atom and the source mass can be isolated. That force is a sum of both the Newtonian gravitational force and any heretofore undiscovered interactions: As such, atom interferometry is a sensitive probe of new theories that predict a classical fifth force F φ . In experiments such as [107,108] the source is a ball of matter and the extra acceleration 95cm is the radius of the ball and d B = 0.88cm is the distance to the interferometer. The whole set up is embedded inside a cavity of radius R c = 6.1cm. Scalar fields, of the type considered in this review, generically predict such a force. The fifth force is of the form where A(φ) is the coupling function to matter. In essence, the source mass induces a nonzero field gradient producing a fifth force, allowing atom interferometry to test scalar field theories. The fifth force depends on the scalar charge q A of the considered object A, i.e. on the way an object interact with the scalar field. In screened theories, it is often written as the product q A = β A m A of the mass m A of the objects and the reduced scalar charge β A . The reduced scalar charge can be factorised as β A = λ A β(φ 0 ) where β(φ 0 ) is a the coupling of a point-particle to the scalar field in the background environment characterised by the scalar field value φ 0 . The screening factor λ A takes a numerical value between 0 and 1 and in general depends on the strength and form of the scalar-matter coupling function A(φ), the size, mass, and geometry of the object, as well as the ambient scalar field value φ 0 . For a spherical object, the screening factor of object A is given by when the object is screened, otherwise λ A = β(φ 0 ). Here φ A is the value of the scalar field deep inside the body A, Φ A the Newtonian potential at its surface and φ 0 is the ambient field value far away from the object. In terms of the screening factors, the force between two bodies A, B is where m c is the effective mass of the scalar particle's fluctuations. In screened theories, the screening factors of macroscopic objects are typically tiny, necessitating new ways to test gravity in order to probe the screened regime of these theories. Atom interferometry fits the bill perfectly [109,110], as small objects like atomic nuclei are typically unscreened. Consequently, screened theories predict large deviations from Newtonian gravity inside those experiments. Furthermore, the experiment is performed in a chamber where the mass m 0 = m(φ 0 ) of the scalar particles is small, and distance scales of order ∼ cm are probed. The strongest bounds are achieved when the source mass is small, approximately the size of a marble, and placed inside the vacuum chamber, as a metal vacuum chamber wall between the bodies would screen the interaction. Within the approximations that led to Eq. (114) one only needs to determine the ambient field value φ 0 inside the vacuum chamber. This quantity depends on the precise theory in question, but some general observations may be made. First, in a region with uniform density ρ, the field will roll to minimise its effective potential V eff (φ) given by (42) for a value φ(ρ). In a dense region like the vacuum chamber walls, φ(ρ) is small, while in the rarefied region inside the vacuum chamber φ(ρ) is large. The field thus starts at a small value φ min,wall near the walls and rolls towards a large value φ min,vac near the center. However, the field will only reach φ min,vac if the vacuum chamber is sufficiently large. The energy of the scalar field depends upon both potential energy V (φ) and gradient energy ( ∇φ) 2 . A field configuration that rolls quickly to the minimum has relatively little potential energy but a great deal of gradient energy, and vice-versa. The ground state classical field configuration is the one that minimises the energy, and hence is a balance between the potential and field gradients. If the vacuum chamber is small, then the minimum energy configuration balances these two quantities by rolling to a value such that the mass of the scalar field is proportional to the size R of the vacuum chamber [49,53] If the vacuum chamber is large, though, then there is plenty of room for the field to roll to the minimum of the effective potential. The condition for this to occur is As such, the field inside the vacuum chamber is It should be noted that in practical experiments, where there can be significant deviations from the approximations used here, i.e. non-spherical source masses and an irregularly shaped vacuum chambers, numerical techniques have been used to solve the scalar field's equation of motion in three dimensions. This enables the experiments to take advantage of source masses that boost the sensitivity of the experiment to fifth forces by some 20% [108]. More exotic shapes have been shown to boost the sensitivity even further, by up to a factor ∼ 3 [111].
Atom interferometry experiments of this type, with an in-vacuum source mass, have now been performed by two separate groups [107,112,113]. In these experiments, the acceleration between an atom and a marble-sized source mass have been constrained to a < ∼ 50 nm/s 2 at a distance of r < ∼ 1 cm. These experiments have placed strong bounds on the parameters of chameleon and symmetron [114] modified gravity, as will be detailed in section III E.

D. Atomic spectroscopy
In the previous section we saw that the scalar field mediates a new force, Eq. (114), between extended spherical objects. This same force law acts between atomic nuclei and their electrons, resulting in a shift of the atomic energy levels. Consequently, precision atomic spectroscopy is capable of testing the modified gravity models under consideration in this review.
The simplest system to consider is hydrogen, consisting of a single electron orbiting a single proton. The force law of Eq. (114) perturbs the electron's Hamiltonian where λ p , m p are the screening factor and mass of the proton, and we have assumed that the scalar field's Compton wavelength m c is much larger than the size of the atom. The electron is pointlike, and is therefore unscreened. 11 The perturbation to the electron's energy levels are computed via the first order perturbation theory result where ψ n are the unperturbed electron's eigenstates. This was first computed for a generic scalar field coupled to matter with a strength β(φ) = m Pl /M [117], using measurements of the hydrogen 1s-2s transition [118][119][120] to rule out However, that study did not account for the screening behavior exhibited by chameleon and symmetron theories. That analysis was recently extended to include screened theories [121], resulting in the bound that is illustrated in Fig. 1.

E. Combined laboratory constraints
Combined bounds on theory parameters deriving from the experimental techniques detailed in this section are plotted in Figs 1 and 2. Chameleons and symmetrons have similar phenomenology, and hence are constrained by similar experiments. Theories exhibiting Vainshtein screening, however, are more difficult to constrain with local tests, as the presence of the Earth or Sun nearby suppresses the fifth force. Such effects were considered in [90] and only restricted to planar configurations where the effects of the Earth are minimised.
The chameleon has a linear coupling to matter, often expressed in terms of a parameter M = m Pl /β. Smaller M corresponds to a stronger coupling. Experimental bounds on the theory are dominated by three tests. At sufficiently small M , the coupling to matter is so strong that collider bounds rule out a wide region of parameter space. At large M > ∼ m Pl , the coupling is sufficiently weak that even macroscopic objects are unscreened, so torsion balances are capable of testing the theory. In the intermediate range the strongest constraints come from atom interferometry. One could also consider chameleon models with n = 1. In general, larger values of n result in more efficient screening effects, hence the plots on constraints would look similar but with weaker bounds overall.
The bounds on symmetron parameter space are plotted in Fig. 2. Unlike the chameleon, the symmetron has a mass parameter µ that fixes it to a specific length scale µ −1 . For an experiment at a length scale L, if L µ −1 then the fifth force would be exponentially suppressed, as is clear in Eq. (114). Likewise, in an enclosed experiment if L µ −1 then the energy considerations in the previous subsection imply that the field simply remains in the symmetric phase where φ = 0. The coupling to matter is quadratic, so in the symmetric phase where φ = 0 the coupling to matter switches off and the fifth force vanishes. Therefore, to test a symmetron with mass parameter µ one must test it with an experiment on a length scale L ≈ µ −1 .

F. Quantum constraints
Classical physics effects induced by light scalar field have been detailed so far. It turns out that laboratory experiments can also be sensitive to the quantum properties of the scalar field. This can typically be seen in two types of situations. In particle physics, the scalars are so light compared to accelerator scales that light scalars can be produced and have a phenomenology very similar to dark matter, i.e. they would appear as missing mass. They could also play a role in the precision tests of the standard model. As we already mentioned above, when the scalars are light compared to the inverse size of the laboratory scales, we can expect that they will induce quantum interactions . The tan curves derive from atom interferometry, while the green region is ruled out by torsion balances. Both the tan and green regions apply only to µ = 10 −1 meV. On the right, only the red curves have been conclusively ruled out by bouncing neutrons. Reproduced from [113] and [122].
which tells us that the light scalar must originate from a completion at energies much larger than the standard model scale.
Quantum effects are also important when the light scalars are strongly coupled to the walls of the Casimir or the Eötwash experiment and light enough in the vacuum between the plates. The mass of the scalar field is given by The density is piece-wise constant and labelled ρ 1,2,3 in the case of a Casimir experiment. Here ρ 2 is the density between the plates. Notice that as φ is continuous, the mass jumps across boundaries as ρ varies from the vacuum density to the plate one. The force between two objects can be calculated using a path integral formalism which takes into account both the classical effects already investigated in this review and the quantum effects akin to the Casimir interaction [124] where the integration is taken over all space and ∂ d ρ is the derivative in the direction defined by the parameter d which specifies the position of one of the bodies. Varying d is equivalent to changing the distance between the objects. For instance, in the case of a plate of density ρ 3 positioned along the x-axis between x = d and x = d + L, the vacuum of density ρ 2 between x = 0 and L, a plate of density ρ 1 for x ≤ 0 and finally again the vacuum for ) and the force is along the x-axis. The quantum average A(φ) is taken over all the quantum fluctuations of φ. When the field has a classical profile φ clas , this quantum calculation can be performed in perturbation theory The first contribution leads to the classical force that we have already considered. The second term is the leading quantum contribution. Notice that the linear coupling in A is absent as the quantum fluctuations involve the fluctuations around a background which satisfies the equations of motion of the system. The higher order terms in the expansion of A(φ) in a power series are associated to higher loop contribution to the force when the first term is given by a one-loop diagram. The Feynman propagator ∆(x, x) at coinciding points is fraught with divergences. Fortunately, they cancel in the force calculation as we will see. Let us focus on the one dimensional force as befitting Casimir experiments. The quantum pressure on a plate of surface area A is then given by where we have considered that the derivative A (φ clas ) A is nearly constant. This is exact for symmetron models and chameleon models with φ M . As the classical solution is continuous at the boundary between the plates, the quantum force is in fact given by where m 3 is the mass of the scalar close to the boundary and inside the plate whereas m 2 is the mass close to the boundary and in the vacuum. As the quantum divergence of ∆(x, x) are x-independent, we see immediately that they cancel in the force (127) which is finite. Moreover, the limit L → ∞ is finite and corresponds to the case of an infinitely wide plate. Notice that the contribution in −∆(d + L, d + L) is the usual renormalisation due to the quantum pressure exerted to the right of the very wide plate of width L.
In the case of a Casimir experiment between two plates, the Feynman propagator with three regions (plate-vacuumplate) must be calculated. In the case of the Eötwash experiment where a thin electrostatic shield lies between the plate, the Feynman propagator is obtained by calculating a Green's function involving five regions. In practice this can only be calculated analytically by assuming that the mass of the scalar field is nearly constant in each of the regions. This leads to the expression with γ 2 i = ρ 2 + m 2 i . When the density in the plates becomes extremely large compared to the one in the vacuum, the limit m 1,3 → ∞ gives the finite result For massless fields in the vacuum m 2 = 0, this gives the Casimir interaction (88) as expected.
When applying these results to screened models, care must be exerted as they assume that the mass of the field is constant between the plates. The quantum contributions to the pressure F x /A can be constrained by the Casimir experiments and the resulting torque between plates by the Eötwash results. These are summarised in figures 3 for symmetrons. In a nutshell, when the µ parameter of the symmetron model becomes lower than 1/d, the field typically vanishes everywhere. The linear coupling to matter vanishes but A = 1/M 2 is non-vanishing thus providing the quadratic coupling to the quantum fluctuations. As the density between the plate is small but non-zero, the mass of the scalar remains positive and the quantum calculation is not plagued with quantum instabilities. For chameleons, the coupling can be taken as A 1/M 2 too. The main difference is that when the density between the plates is low, the mass of the scalar cannot become much lower than 1/d, see (116), implying that the quantum constraints are less strong than in the symmetron case.
As the expansion of A(φ) involves higher order terms suppressed by the strong coupling scale M and contributing to higher loops, they can be neglected on distances between the plates d > ∼ 1/ m 1,3 M . As the density in the plates is very large, this is always a shorter distance scale than 1/M where the calculations of the effective field theory should not be trusted naively. In the limit m 1,3 → ∞ the one loop result becomes exact and coincide with (half) the usual Casimir force expression for electrodynamics as obtained when the coupling to the boundaries is also very strong and Dirichlet boundary conditions are imposed.
Finally, measurements of fermions' anomalous magnetic moments are sensitive to the effects of new scalar fields coupled to matter. The anomalous magnetic moment is where g f is the fermion's g-factor. There are two effects to consider. First is the well-known result that at 1-loop the scalar particle corrects the QED vertex, modifying the anomalous magnetic moment by an amount [117,125,126] where m f is the mass of the fermion. Second, the classical fifth force introduces systematic effects in the experiment, such as a modified cyclotron frequency, that must be accounted for in order to infer the correct measured value of a f [126,127].
In the case of the electron, the measurement of a e and the Standard Model prediction agree at the level of 1 part in 10 12 [128]. Setting δa e ≤ 10 −12 yields the constraint [126] β(φ) < ∼ 10 16 .
In the case of the chameleon where β = m Pl /M , this rules out M < 80 GeV.
Cosmologically, a chameleon with these parameters has an effective mass m cosmo > 10 −13 eV and Compton wavelength < 10 3 km, so this theory does not significantly influence our universe on large scales.

IV. ASTROPHYSICAL CONSTRAINTS AND PROSPECTS
In this section we discuss the ways in which screened fifth forces may be searched for using astrophysical objects beyond the solar system, specifically stars, galaxies, voids and galaxy clusters. We describe the tests that have already been conducted and the ways in which they may be strengthened in the future. Astrophysical constraints are most often phrased in terms of the n = 1 Hu-Sawicki model of f (R) (taken as a paradigmatic chameleon-screened theory; [69,151,152]) and nDGP or a more general galileon model (taken as paradigmatic Vainshtein-screened theories; [59,73]).
Testing screening in astrophysics requires identifying unscreened objects where the fifth force should be manifest. Ideally this would be determined by solving the scalar's equation of motion given the distribution of mass in the universe, although the uncertainties in this distribution and the model-dependence of the calculation make more approximate methods expedient. This may be done by identifying proxies for the degree of screening in certain theories which can be estimated from the observed galaxy field. In thin-shell screening mechanisms (chameleon, symmetron and the environmentally-dependent dilaton) it is the surface Newtonian potential of an object relative to the background scalar field value that determines whether it is screened (as discussed in Sec. II). This screening criterion may be derived analytically for an object in isolation or in the linear cosmological regime (e.g. [48,49] for the chameleon), while N-body simulations in modified gravity have shown that it is also approximately true in general when taking account of both environmental and self-screening [153][154][155] (see Fig. 4). The threshold value of potential for screening is given by Eq. 52: in n = 1 Hu-Sawicki f (R), χ 3 2 f R0 so that probing weaker modified gravity (lower f R0 ) requires testing objects in weaker-field environments [69]. Rigorous observational screening criteria are not so easy to derive in other screening mechanisms, although heuristically one would expect that in kinetic mechanisms governed by nonlinearities in the first derivative of the scalar field it is the first derivative of the Newtonian potential (i.e. acceleration) that is relevant while in Vainshtein theories governed by the second derivative of the field it is instead the spacetime curvature (Sec. II E).
Several methods have been developed to build "screening maps" of the local universe to identify screened and unscreened objects. Shao et al. [156] apply an f (R) scalar field solver to a constrained N-body simulation to estimate directly the scalar field strength as a function of position. Cabre et al. [153] use galaxy group and cluster catalogues to estimate the gravitational potential field and hence the scalar field by the equivalence described above. Desmond et al. [157] adopt a similar approach but include more contributions to the potential, also model acceleration and curvature, and build a Monte Carlo pipeline for propagating uncertainties in the inputs to uncertainties in the gravitational field. By identifying weak-field regions these algorithms open the door to tests of screening that depend on local environment, with existing tests using one of the final two.

A. Stellar tests
Gravitational physics affects stars through the hydrostatic equilibrium equation, which describes the pressure gradient necessary to prevent a star from collapsing under its own weight. In the Newtonian limit of GR, this is given . The x-axis is the total halo mass in M , the y-axis is the Newtonian potential sourced at the halo's position by mass within one Compton wavelength of the scalar field, and the points are colour-coded by the degree of screening with red fully unscreened and dark blue fully screened. The vertical and horizontal lines mark where the internal and external potentials equal 3 2 fR0, showing that these cuts can reliably separate screened from unscreened galaxies. Reproduced from [153].
In the presence of a thin-shell-screened fifth force this becomes with Θ(x) the Heaviside step function, β the coupling coefficient of the scalar field and r s the screening radius of the star beyond which it is unscreened. In the case of chameleon theories, the factor 1 − M (rs) M (r) corresponds to the screening factor and is associated to the mass ratio of the thin shell which couples to the scalar field. The stronger inward gravitational force due to modified gravity requires that the star burns fuel at a faster rate to support itself than it would in GR, making the star brighter and shorter-lived. The magnitude of this effect depends on the mass of the star: on the main sequence, low-mass stars have L ∝ G 4 N while high-mass stars have L ∝ G N [158]. Thus in the case that the star is fully unscreened (r s = 0), low-mass stars have L boosted by a factor (1 + 2β 2 ) 4 , and high-mass stars by (1 + 2β 2 ).
To explore the full effect of a fifth force on the behaviour of stars, Eq. 134 must be coupled with the equations describing stellar structure and energy generation. This has been done by modifying the stellar structure code MESA [158][159][160][161], enabling the heuristic expectations described above to be quantified (see Fig. 5). The expectation that stars are brighter in modified gravity-and low-mass stars more so than high-mass-also leads to the prediction that unscreened galaxies would be more luminous and redder than otherwise identical screened ones. No quantitative test has been designed around this though because no galaxy formation simulation including the effect of modified gravity on stars has yet been run.
Fifth forces also have important effects in astroseismology, the study of stellar oscillations. The equation of motion for small perturbations of mass elements in stars is¨ with a the force per unit mass, which is a = − ∇Φ in GR but a = − ∇Φ − β m Pl ∇φ in the presence of a scalar field. Combining this equation with the other stellar structure equations gives the frequency of linear radial adiabatic oscillations so that enhancing the effective value of G due to the addition of a fifth force causes the pulsation period Π to change according to where Q is the star's scalar charge. Stellar oscillations are useful observationally because they provide several methods of determining distances to galaxies [162]. These afford a test of gravity when multiple distance indicators with different screening properties are combined. In particular, if a distance indicator is sensitive to G N and calibrated assuming GR, it will fail to give the correct distance to an unscreened galaxy in a fifth-force theory. This will lead to a discrepancy with the distance estimated using an indicator that is not sensitive to G N , e.g. because it is based on the physics of a high-density, screened object. This test has been carried out by comparing Cepheid and TRGB (Tip of the Red Giant Branch) distance indicators. Cepheids are post-main-sequence stars that oscillate radially by the κ-mechanism [163] when crossing the instability strip in the Hertzsprung-Russell diagram. The period of this pulsation is tightly correlated with the luminosity of the star, allowing Cepheids to be used as standard candles. TRGB stars are red giants that have become sufficiently hot for helium fusion to occur, moving the star onto the horizontal branch of the Hertzsprung-Russell diagram and leaving an observable discontinuity in the I-band magnitude. This occurs at an almost fixed absolute magnitude, making the TRGB feature another standard candle. The TRGB luminosity is sourced by a thin hydrogen-burning shell surrounding the helium-burning core, so if the core is screened then TRGBs exhibit regular GR behaviour. This occurs for χ < ∼ 10 −6 , which is the case for thin-shell theories that pass the tests described below. With Cepheids unscreened down to much lower values of χ, this means that TRGB and Cepheid distances would be expected to disagree in unscreened galaxies. The fact that they seem not to-and that any discrepancy between them is uncorrelated with galaxy environment-has yielded the constraint f R0 < ∼ 10 −7 [164,165]. Notice that astrophysical constraints yield tighter bounds on f (R) models than solar system tests.
Variable stars are also useful for more general tests of gravity. [166] showed that the consistency between the mass estimates of Cepheids from stellar structure vs astroseismology allows a constraint to be placed on the effective gravitational constant within the stars. Using just 6 Cepheids in the Large Magellanic Cloud afforded a 5% constraint on G N , and application of this method to larger datasets spanning multiple galaxies will allow a test of the environment-dependence of G N predicted by screening. Screening may also provide a novel local resolution of the Hubble tension [165,167,168].
Finally, other types of star are useful probes of the phenomenon of "Vainshtein breaking" whereby the Vainshtein mechanism may be ineffective inside astrophysical objects. An unscreened fifth force inside red dwarf stars would impact the minimum mass for hydrogen burning, and a constraint can be set by requiring that this minimum mass is below the below the lowest mass of any red dwarf observed [169,170]. It would also affect the radii of brown dwarf stars and the mass-radius relation and Chandresekhar mass of white dwarfs [171].

B. Galaxy and void tests
Screened fifth forces have interesting observable effects on the dynamics and morphology of galaxies. The most obvious effect is a boost to the rotation velocity and velocity dispersion beyond the screening radius due to the enhanced gravity. This is strongly degenerate with the uncertain distribution of dark matter in galaxies, although the characteristic upturn in the velocity at the screening radius helps to break this. In the case of chameleon screening, Naik et al. [172] fitted the rotation curves of 85 late-type galaxies with an f (R) model, finding evidence for f R0 ≈ 10 −7 assuming the dark matter follows an NFW profile but no evidence for a fifth force if it instead follows a cored profile as predicted by some hydrodynamical simulations. This illustrates the fact that a fifth force in the galactic outskirts can make a cuspy matter distribution appear cored when reconstructed with Newtonian gravity, of potential relevance to the "cusp-core problem" [173] (see also [174]). Screening can also generate new correlations between dynamical variables; for example Burrage et al. [175] use a symmetron model to reproduce the Radial Acceleration Relation linking the observed and baryonic accelerations in galaxies [176]. Further progress here requires a better understanding of the role of baryonic effects in shaping the dark matter distributions in galaxies, e.g. from cosmological hydrodynamical simulations in ΛCDM.
One way to break the degeneracy between a fifth force and the dark matter distribution is to look at the relative kinematics of galactic components that respond differently to screening. Since main-sequence stars have surface Newtonian potentials of ∼ 10 −6 , they are screened for viable thin-shell theories. Diffuse gas on the other hand may be unscreened in low-mass galaxies in low-density environments, causing it to feel the fifth force and hence rotate faster [177,178]: where M (< r) is the enclosed mass and v g and v * are the gas and stellar velocities respectively. We see that comparing stellar and gas kinematics at fixed galactocentric radius factors out the impact of dark matter, which is common to both. Comparing the kinematics of stellar Mgb absorption lines with that of gaseous Hβ and [OIII] emission lines in 6 low-surface brightness galaxies, Vikram et al. [179] place the constraint f R0 < ∼ 10 −6 . This result can likely be significantly strengthened by increasing the sample size using data from IFU surveys such as MaNGA or CALIFApotentially combined with molecular gas kinematics, e.g. from ALMA-and by modelling the fifth force within the galaxies using a scalar field solver rather than an analytic approximation. A screened fifth force also generates asymmetries in galaxies' rotation curves when they fall nearly edge-on in the fifth-force field, although modelling this effect quantitatively is challenging so no concrete results have yet been achieved with it [180].
The strongest constraints to date on a thin-shell-screened fifth force with astrophysical range come from galaxy morphology. Consider an unscreened galaxy situated in a large-scale fifth-force field a φ sourced by surrounding structure. Since main-sequence stars self-screen, the galaxy's stellar component feels regular GR while the gas and dark matter also experience a φ . This causes them to move ahead of the stellar component in that direction until an equilibrium is reached in which the restoring force on the stellar disk due to its offset from the halo centre exactly compensates for its insensitivity to a φ so that all parts of the galaxy have the same total acceleration [177,178]: where r * is the displacement of the stellar and gas centroids. This effect is illustrated in Fig. 6 (a), and can be measured by comparing galaxies' optical emission (tracing stars) to their HI emission (tracing neutral hydrogen gas). A second observable follows from the stellar and halo centres becoming displaced: the potential gradient this sets up across the stellar disk causes it to warp into a characteristic cup shape in the direction of a φ . This is shown in Fig. 6 (b). The shape of the warp can be calculated as a function of the fifth-force strength and range, the environment of the galaxy and the halo parameters that determine the restoring force: which can be simplified on assuming a halo density profile. Desmond et al. [181][182][183] create Bayesian forward models for the warps and gas-star offsets for several thousand galaxies observed in SDSS and ALFALFA, including Monte Carlo propagation of uncertainties in the input quantities and marginalisation over an empirical noise model describing non-fifth-force contributions to the signals. This method yields the constraint f R0 < 1.4 × 10 −8 at 1σ confidence, as well as tight constraints on the coupling coefficient of a thin-shell-screened fifth force with any range within 0.3-8 Mpc [184] (see Fig. 7 (a)). Subsequent work has verified using hydrodynamical simulations that the baryonic noise model used in these analyses is accurate [185]. 10 −8 is around the lowest Newtonian potential probed by any astrophysical object, so it will be very hard to reach lower values of f R0 . Lower coupling coefficients may however be probed using increased sample sizes from upcoming surveys such as WFIRST, LSST and SKA, coupled with estimates of the environmental screening field out to higher redshift using deeper wide photometric surveys. The above tests target thin-shell-screened fifth forces. The Vainshtein mechanism is harder to probe due to the efficiency of its screening on small scales and the difficulty of developing robust observational proxies for objects' degrees of screening. While LLR is sensitive to cubic galileons with small crossover scale r c ∼ O(100) kpc [37], the larger values r c ∼ 6000 Mpc required for self-acceleration [186] must be probed on galactic or cosmological scales. The most promising method for doing this utilises the breaking of the Strong Equivalence Principle (SEP) that galileons imply [187] in the presence of black holes. Galileons couple to the trace of the stress-energy tensor, which is equivalent to density but excludes gravitational binding energy. This means that non-relativistic objects (e.g. stars, gas and dark matter in galaxies) have a scalar charge-to-mass ratio equal to the coupling coefficient β, while black holes are purely relativistic objects with Q = 0. Thus in the presence of an unscreened large-scale galileon field, the supermassive black holes at galaxies' centres will lag behind the rest of the galaxy, which is measurable by comparing the galaxies' light with radio or X-ray emission from the Active Galactic Nuclei (AGN) powered by the black hole. Two situations can lead to an unscreened galileon field. The first is in galaxy clusters environments: an extended distribution of mass does not Vainshtein-screen perfectly in its interior [188], so a residual fifth-force field is present in cluster outskirts. This leads to O(kpc) offsets between black holes and satellite galaxy centres for realistic cluster parameters. Sakstein et al. [189] solve the galileon equation of motion for a model of the Virgo cluster and use the fact that the black hole in the satellite galaxy M87 is within 0.03 arcsec of the galaxy centre to rule out O(1) coupling coefficients for r c < ∼ 1 Gpc. Second, the galileon symmetry implies that the linear contribution to the field on cosmological scales is unscreened [190,191], allowing black hole offsets to develop even for field galaxies. Assuming a constant density ρ 0 in the centre of the halo, the black hole offset in this case is given by [187] where ∇Φ ext N is the unscreened large-scale gravitational field, proportional to the galileon fifth-force field. Bartlett et al. [192] modelled this field using constrained N-body simulations of the local ∼200 Mpc and forward-modelled the offsets in 1916 galaxies with AGN, including a more sophisticated model for the halo density profiles, to set the bound β < 0.28 for r c > ∼ 1/H 0 (see Fig. 7 (b)). This probes the cosmologically-relevant region of the galileon parameter space, complementing cosmological probes such as the Integrated Sachs Wolfe (ISW) effect (see Sec. V). It could be improved to probe smaller r c values by modelling the full, non-linear dynamics of the galileon within the test galaxies. Another possible signature is "wandering" black holes seemingly unassociated with a galaxy [189].
While galaxies are the directly observable tracers of the cosmic web, much dynamical information can be found in voids, the underdense regions that comprise most of the universe's volume. These are particularly promising for   [184]. (b) Constraints on the galileon coupling coefficient (on the plot denotes by α) as a function of the theory's cross-over scale. The orange region is excluded by Lunar Laser Ranging, the green region from the position of the supermassive black hole in M87, and the blue region from a statistical analysis of the black hole positions in field galaxies. Reproduced from [192].
testing screening because they are the regions where it is least efficient. Their usefulness is however hampered by the ambiguity that exists in defining voids, and by the fact that voids must be identified observationally using biased tracers of the underlying density field (galaxies). Voids in modified gravity have been studied through both analytic [193,194] and simulational [195,196] methods. Typically, the enhanced growth of structure in the presence of a fifth force causes voids to become larger and emptier. In addition, when voids are identified through lensing the modified relation between mass and lensing potential can affect the lensing signal [78,197]. Voids can also be cross-correlated with galaxies to infer the growth rate of structure [198], used in the ISW effect [199], integrated along the line of sight to produce projected 2D voids [200], and used as a means of splitting samples of galaxies into high density (screened) and low density (unscreened) environments or in marked correlation functions [201,202]. Finally, the redshift-space anisotropy of voids is a powerful probe of the nature of gravity through redshift space distortions [203]. Future surveys will improve 3D spectroscopic voidfinding and the calibration of photometric void finders with robust photometric redshifts.

C. Galaxy cluster tests
A fifth force causes structure to grow more quickly, leading to more cluster-sized halos at late times. This is however counteracted by screening and the Yukawa cutoff due to the mass of the scalar field so that cluster abundance only deviates significantly from the ΛCDM expectation at lower masses and in sparser environments [204]. The excursion set formalism for halo abundance provides a good description under chameleon gravity as well [205], albeit with a modified collapse threshold δ c , and has been used to constrain f R0 < ∼ 10 −5 in the Hu-Sawicki model [206,207]. Similar constraints are achievable using the peaks in the weak lensing convergence field, which trace massive halos [208]. Other formalisms for calculating cluster abundance in the presence of a fifth force have also been developed [209][210][211][212]. Qualitatively similar results hold for Vainshtein-screened theories, where, although the centres of clusters are efficiently screened, massive halos grow at an increased rate because of enhanced accretion due to the fifth force in the surrounding matter [213]. This can be significantly altered for K-mouflage models where clusters are not screened so we expect massive halos to be more abundant than in ΛCDM. This is illustrated in Fig. 8; the "arctan" models are particularly interesting because they pass the solar system tests.
The internal structures of cluster halos are also altered by modified gravity, particularly through an increase in the concentration of the Navarro-Frenk-White profile [215][216][217], although this is hard to use to set constraints due to degeneracies with the impact of baryons. Another important effect is on the boundary of cluster halos, namely the splashback radius where accreting dark matter turns around after first infall [218]. This is marked by a sharp drop in the logarithmic density slope, and consequently in the lensing signal and subhalo profile. Adhikari et al. [219] studied the splashback feature in both chameleon and symmetron models, finding that for viable and interesting values of the fifth-force properties the splashback radius is increased relative to GR in Vainshtein models and reduced in chameleon. This results from competition between the enhanced acceleration of accreting matter and reduced dynamical friction within the halos. There is however controversy observationally about the location of the cluster splashback radius [220][221][222], so these predictions cannot be used to set concrete constraints. Further out, the halo-matter correlation function is enhanced by a fifth force [206,223]. A powerful and general method for probing modified gravity leverages the inequality between the two weak-field metric potentials, a violation of the weak equivalence principle. This leads to a difference between the dynamical and lensing masses of objects: while photons respond to the sum of the metric potentials, nonrelativistic tracers are affected solely by the time-time potential. Thin-shell screening alters the Newtonian potential but not the lensing one, which in the Parametrised Post-Newtonian framework is captured by the parameter γ. Although γ may be constrained on O(kpc) scales by comparing strong lensing and stellar dynamical mass estimates [224,225], it has found most use on the scale of clusters. An approximation for chameleon theories of the Jordan-Brans-Dicke type is [226] M dyn (r) where Θ is the Heaviside step function, ω BD is the JBD parameter (see Sec. II G 5) and the radius at which the scalar field transitions to its background value is given by Here φ env is the cosmological boundary condition for the field far from the cluster (e.g. 1 − A −2 (φ) f R0 in the f (R) case) and r s is the scale length of the cluster's assumed-NFW density profile. The difference between dynamical and "true" masses of clusters in f (R) gravity has also been calibrated from N-body simulations in [227]: where p 1 = 2.21 and p 2 = 1.503 log 10 |f R (z)−1| 1+z + 21.64. This works well for f R0 ∈ 10 −6.5 , 10 −4 and z ∈ [0, 1]. To test this effect, cluster strong lensing may be compared to X-ray masses or the velocity dispersions of the cluster galaxies [188,228], and stacked weak lensing can be compared to Sunyaev-Zel'dovich masses or infall motions at the cluster outskirts [229]. Dynamical masses can also be estimated from X-ray data of cluster temperature and pressure profiles. The combination of weak lensing measurements with measurements of the X-ray brightness, temperature and Sunyaev-Zel'dovich signal from the Coma cluster [230] (or from multiple clusters' weak lensing and X-ray signals [231]) implies f R0 < ∼ 6 × 10 −5 , and this test has also been applied to galileons [232]. The modification to clusters' dynamical masses under a fifth force can be probed without requiring lensing data by assuming that the gas fractions of clusters are constant in order to estimate the true total mass. This is capable of constraining f (R) to the f R0 ∼ 5 × 10 −5 level [233]. All of these tests will benefit from enlarged cluster samples in the future.

A. Screening and cosmic acceleration
Screened fifth forces coupled to matter also have interesting cosmological consequences. In the modified gravity models studied above, the screening mechanisms are necessary to make the models consistent with observations at small scales. As detailed in sections II C and II D we can classify the screening types into non-derivative and derivative screening mechanisms. From the former the chameleon is the most popular example, appearing in popular models such as Hu-Sawicki f (R). For the latter, the Vainshtein and K-mouflage mechanisms are the characteristic ones, appearing in subsets of Horndeski theory, such as models with a modified kinetic term (for K-mouflage) or models like Cubic Galileons, which feature the Vainshtein screening as a way to evade small scale gravitational constraints.
No-go theorems [35,87,234] were developed for chameleon-screened theories, and they state namely that i) the Compton wavelength of such scalars can be at most 1Mpc at the present cosmic density, which means that the effective range of these theories is restricted to nonlinear scales in large scale structure formation and they have no effect on the linear growth of structures; and (ii) that the conformal factor (64) relating the Einstein and Jordan frames of these theories is essentially constant in one Hubble time, therefore these scalar fields cannot be responsible for self-acceleration and one needs to invoke either a cosmological constant term or another form of dark energy to explain the acceleration of the expansion of the Universe. More precisely, in the context of chameleon-screened models one can show that the equation of state of dark energy at late times is of order [53] where m is the mass of the light scalar. The bound from solar systems on the mass ratio m/H > ∼ 10 3 coming from solar system tests, see (80), implies that the equation of state is indistinguishable from the one of a cosmological constant. On the other hand, these theories have effects on large scale structures and then irrespective of what would drive the acceleration one could test the screening effects at the cosmological perturbation level.
In the second class of models, the scalar field evolves significantly on cosmic timescales, as in the case of cubic Galileons, kinetic gravity braiding models and K-mouflage models. These models present either K-mouflage or Vainshtein screenings and therefore are not affected by the no-go theorems.
In the following sections we will present the different ways in which these screened modified gravity theories affect cosmological observables and the current and future bounds that can be placed on their parameters.

B. Screening and structure formation
The formation of large scale structure is affected by the presence of modified gravity. Screening could play a role too as we will see below as the growth of structure depends on the type of screening mechanisms. For derivative screening, the growth is affected at the linear level in a scale independent way. For non-derivative screenings, the linear growth is modified in a scale dependent way. The latter can be easily understood as there is a characteristic length scale, i.e. the Compton wavelength of the scalar field, beyond which modified gravity is Yukawa-suppressed. Non-linear effects are also important and tend to dampen the effects of modifies gravity on small scales.
As an example and on cosmological scales the f (R) modification of the Einstein-Hilbert action leads to a modified Poisson equation, which can be expressed as in comoving coordinates and the term δρ is the matter density fluctuation compared to the cosmological background and Φ the modified Newtonian potential. Furthermore, the fluctuation of the Ricci scalar, δR = R −R compared to the cosmological backgroundR and is expressed as The variation of the function f (R) is given by δf R = f R (R) − f R (R). In these equations we have assumed a quasistatic approximation. It can be shown [235] that despite the fact that these equations are nonlinear in δR, they are self-averaging. This means that on large scales one recovers δR → 0. Using these governing equations one can solve perturbatively the Vlasov-Poisson system of equations, which consists in first approximation (no vorticity and single-stream regime) of the continuity, Euler and Poisson equations, in powers of the linear growth factor. The results of these computations at 1-loop order and beyond can be seen in references [235][236][237][238][239][240].
In scalar-tensor theories with screening and a conformal factor A(φ) particles feel a total gravitational potential Φ which is the sum of the standard Newtonian term Φ N and an additional contribution Φ A , where the governing equations are given by where it is assumed that A(φ) 1, to satisfy constraints on the variation of fermionic masses. As a result ln A A − 1 and the dependence on ln A of the Newtonian potential Φ becomes linear in A. This additional gravitational potential implies that matter particles of mass m are sensitive to a "fifth force" given by This fifth force is the one which leads to a modification of the growth of structures.

C. Cosmological probes: CMB and large scale structure
Historically, the background expansion of the Universe has been the traditional way of testing cosmological models and this has been developed mostly through the study of standard candles, especially with the use of observations of supernovae SNIa [3,241]. However, recent constraints on the equation of state parameter of dark energy, are overall consistent with a cosmological constant w ≈ −1 [242]. This, plus the fact that self-acceleration is mostly ruled out in the most popular screened scalar field models, has led to the tendency in the literature to look for features of dark energy and modified gravity in the formation of structures and the modification of gravitational lensing. Moreover, other interesting tensions in the data, such as the H 0 tension [243], cannot be satisfactorily resolved with late-time dynamics of a dark energy field, according to the latest analysis [244,245] and therefore will not be covered in this section. Therefore, in the following we will concentrate mostly on the Integrated Sachs Wolfe effect in the CMB, lensing of the CMB and the formation of structures probed by the Galaxy power spectrum and its effect on Weak Lensing (cosmic shear).

ISW and CMB Lensing
The relic radiation from the early Universe that we observe in the GHz frequency range, called the Cosmic Microwave Background is one of the most powerful cosmological probes. It constrains not only the background of the Universe, but also its growth of structure. Its primary anisotropies, imprinted at the time of recombination, provide plenty information about the constituents of the Universe; while its secondary anisotropies, that happen later when the CMB photons are traversing the Universe, provide information about the intervening medium, the expansion of the Universe and the large scale structures. For studying late modified gravity and dark energy, these secondary anisotropies are the most important probes, namely the Integragted Sachs-Wolfe effect (ISW) ( [246][247][248] that affect the power spectrum at low multipoles (large scales) and lensing of the CMB [249,250] that affects the spectrum at small scales (high multipoles).
In the case of ISW, the effect is observed as a temperature fluctuation caused by time variations in the gravitational potentials that are felt by photons when they enter and leave potential wells (or potential hills) when entering dark matter halos (or voids). The effect on the CMB temperature T is given by where η * is the conformal time at the last scattering surface and η 0 at the observer. By changing the time evolution of the gravitational potentials, MG models affect the large scales of the CMB power spectrum through the ISW effect.
The ISW effect played a major role in ruling out cubic Galileon models, which are the only non-trivial part left from the Horndeski theory after GW170817. In [213] cubic Galileons were analyzed and it was found that in the presence of massive neutrinos (model dubbed νGalileon, in red in Fig. 9), the models were still a very good fit to CMB temperature, lensing and Baryon Acoustic Oscillation (BAO) data, using Planck-2013 temperature and lensing [251] and WMAP-9 polarization [252] data. For BAO they used 6dF, SDSS DR7 and BOSS DR9 data ( [253][254][255]). In the absence of massive neutrinos (model dubbed Galileon in Fig.9), however, ΛCDM was favoured by data. Nevertheless, they showed that the νGalileon model shows a negative ISW effect that is hard to reconcile with current observations. More recently, a paper by [256] performed a detailed study of self-accelerating Galileon models using CMB data from Planck-15 in temperature and polarization and CMB lensing [257]. They also included BAO data, H 0 data and ISW data. As in the older analysis, they showed that the cubic Galileon predicts a negative ISW effect and therefore it is in a 7.8σ tension with observations, effectively ruling this model out. Furthermore, in [258] the effect of different neutrino masses and hierarchies was analyzed and it was also found out that all cubic, quartic and quintic Galileons remain ruled out by CMB and ISW observations.
FIG. 9: CMB temperature power spectrum. In black dots, data from Planck-2013, in blue the cubic Galileon model without massive neutrinos, in red the same model in presence of massive neutrinos and in green baseline ΛCDM with standard neutrino mass. The different between solid and dashed lines corresponds to an analysis of Planck with and without BAO data, respectively. Reproduced from [213] thankfully provided by Alex Barreira.

Cosmological perturbations in large scale structure
As mentioned above in the corresponding sections for f (R) and scalar field models, the dynamics of the field at large scales is given by the Poisson equation and the corresponding Klein-Gordon equation. However, when including the full energy-momentum tensor, the first-order perturbed Einstein equations in Fourier space give two equations that describe the evolution of the two gravitational potentials Φ and Ψ. In the quasistatic approximation, these equations read Where ρ(a) is the average dark matter density and ∆(a, k) = δ + 3aHθ is the comoving density contrast with δ the fractional overdensity, and θ the peculiar velocity. We have denoted by the lensing potential. The ratio of the two gravitational potentials is denoted as η, gravitational anisotropic stress or gravitational slip η(a, k) ≡ Ψ(a, k)/Φ(a, k) .
The scale and time-dependent functions η(a, k), µ(a, k) and Σ(a, k) stand for all possible deviations of Einstein gravity in these equations, being equal to unity when standard GR is recovered, and can encompass any modification by a scalar-tensor theory at the linear level in perturbations. Given that there are only two scalar degrees of freedom, it means that of course there is a relationship between µ, Σ and η and they are related by The µ(a, k) function is usually probed best by Galaxy Clustering experiments that directly trace the evolution of the Φ potential, since this one affects non-relativistic particles. µ is directly related to the effective Newtonian constant defined above in (9) as in the linear regime and in Fourier space. On the other hand, relativistic particles, and therefore light, follows the equation for Φ(a, k) + Ψ(a, k), meaning that gravitational weak lensing is mostly sensitive to the function Σ (a, k). a. f (R) models and chameleon theories : For the f (R) theories described above these expressions reflect the presence of an extra fifth force. In particular, it is convenient to introduce the mass of the scalaron field, i.e. the scalar field associated to the f (R) models [53] where we have used that R ρ/m 2 Pl at late time in Universe. Neglecting the anisotropic stress the expressions for µ and η read [? ] Given the constraints on f R,0 mentioned above, the modifications of lensing are practically non-existent and Σ(a, k) 1 with great precision. It is convenient to rewrite the above expressions as µ(a, k) = A 2 (a) 1 + 2β(a) 2 where in the case of f (R) models we have β(a) = β = 1/ √ 6 and m f R (a) = m(a) where β(a) is the coupling at the minimum of the effective potential V eff (φ) = V (φ) + (A(φ) − 1)ρ as a function of a with ρ ∝ a −3 and similarly for the mass of the scalar field m(a). These expressions are valid for any chameleon theories.
In all chameleon theories, there is a one-to-one correspondence between the coupling and mass variations as a function of the scale factor and the potential V (φ) and coupling function A(φ) which is called the tomographic map. This allows to parameterise the chameleon models with the function m(a) and β(a). The mapping reads [87] where Ω m is the matter fraction of the Universe. In this expression, the matter fraction and the Hubble rate can be taken as the ones of the standard model as solar system tests imply that chameleon models essentially coincide with ΛCDM at the background level. The potential itself is given by This provides a parametric reconstruction of V (φ). For the Hu-Sawicki models of f (R), we have [53] m fR (a) = m 0 Ωm0 where Ω Λ is the dark energy fraction and Ω m0 the matter fraction now. The mass of the scalaron now is given by which is greater than H 0 for small |f R0 | 1. Finally, the µ parameterisation allows one to see how screening works on cosmological linear scales [50]. Defining the comoving Compton wavelength λ c (a) = 1 am(a) (168) we find that for scales outside the Compton wavelength, i.e. k < ∼ λ c we have µ(a, k) A 2 (a) 1 (169) and GR is retrieved. This corresponds to the Yukawa suppression of the fifth force induced by the light scalar. On the contrary when k > ∼ λ c , we have an enhancement of the gravitational interaction as µ(a, k) 1 + 2β 2 (a) (170) which is simply due to the exchange of the nearly-massless scalar field between over densities. As a result, we can have qualitative description of chameleon models such as f (R) on the growth of structures [259]. First of all on very large scales, GR is retrieved and no deviation from Λ-CDM is expected. On intermediate scales, deviations are present as (170) is relevant. Finally on much smaller scales the screening mechanism prevents any deviations and GR is again retrieved. The onset of the modified gravity regime is set by the mass of the scalar now which is constrained by the solar system tests to be in the sub-Mpc range. This falls at the onset of the non-linear regime of growth formation and therefore one expects the effects of modified gravity to be intertwined with non-linear effects in the growth process.
b. Jordan-Brans-Dicke models : For the JBD models with a mass term these functions are given by [260,261] µ(a, k) = 1 φ so that for cosmological purposes Σ(a) = 1 .
In this case lensing is not affected at all. c. Horndeski theory : For a generic Horndeski theory (of second order in the equations of motion) these two functions µ and η can be expressed as a combination of five free functions of time p 1,2,3,4,5 , which are related to the free functions G i in the Horndeski action [260,262] µ(a, k) = η(a, k) = 1 + p 3 (a)k 2 p 4 (a) + p 5 (a)k 2 .
There is another physically more meaningful parametrization of the linear Horndeski action, given by [263] which is related to the Effective Field Theory of dark energy [32,264,265], where small deviations to the background cosmology are parameterised linearly. This parametrization is of great help when discussing current cosmological constraints. It is defined using four functions of time α M , α K , α B and α T plus the effective Planck mass M 2 and a function of time for a given background specified by the time variation of the Hubble rate H(a) as a function of the scale factor a. The term α T measures the excess of speed of gravitational waves compared to light and therefore as we previously mentioned, after the event GW170817, this term is constrained to be effectively zero. The term α K quantifies the kineticity of the scalar field and therefore appears in models like K-mouflage, which require the K-mouflage screening.The coefficient α B quantifies the braiding or mixing of the kinetic terms of the scalar field and the metric and can cause dark energy clustering. It appears in all modified gravity models where a fifth force is present [266]. It receives contributions also from terms related to the cubic Galileons, which present the Vainshtein screening. Finally, α M quantifies the running rate of the effective Planck mass and it is generated by a non-minimal coupling. This parameter modifies the lensing terms, since it directly affects the lensing potential. It appears in f (R) models, where the chameleon screening is necessary, as we have seen. d. DGP models : Cosmological linear perturbations for DGP have been worked out in [267]. In the paper by [260], it is assumed that the small-scale (quasi-static) approximation is valid, i.e. k/a r 5 H and is obtained and where γ = 1 + 2 Hr 5 w eff . This corresponds to and for all practical purposes we can set Σ = 1 within the cosmological horizon (see [260]).

D. Large scale structure observations: Galaxy Clustering and Weak Lensing
The most important probes for large scale structure, especially in the upcoming decade with the advent of new observations by DESI [268] 12 , Euclid [269,270] 13 , Vera Rubin [271] 14 and WFIRST [272] 15 , will be galaxy clustering and weak lensing. Galaxy clustering measures the 2-point-correlation function of galaxy positions either in 3 dimensions, i.e. angular positions and redshift or in effectively 2 dimensions (angular galaxy clustering) when the redshift information is not particularly good. In Fourier space, this correlation function of galaxies, known as the observed galaxy power spectrum P obs gg is directly related to the power spectrum of matter density perturbations P δδ,zs in redshift space by P obs gg (z, k, µ θ ) = AP(z)P δδ,zs (k, z)E err (z, k) + P shot (z) , where AP (z) corresponds to the Alcock-Paczynski effect, E err (z, k) is a damping term given by redshift errors and P shot (z) is the shot noise from estimating a continuous distribution out of a discrete set of points. µ θ is the cosine of the angle between the line of sight and the wave vector k. Furthermore, the redshift space power spectrum, is given by where FoG(z, k, µ θ ) is the "Fingers of God" term that accounts for non-linear peculiar velocity dispersions of the galaxies and K is the redshift space distortion term that depends -in linear theory, where it is known as the Kaiser term [273] -on the growth rate f (z) and the bias b(z), but can be more complicated when takling into account nonlinear perturbation theory at mildly nonlinear scales. For a detailed explanation of these terms, we refer the reader to [274] and the many references therein. Relativistic effects in galaxy clustering may provide a particularly sensitive probe of fifth forces and screening. With relativistic effects included, the cross-correlation of two galaxy populations with different screening properties yields a dipole and octopole in the correlation function due to the effective violation of the weak equivalence principle -as encapsulated in Euler's equation -as the galaxies in the two groups respond differently to an external fifth-force field [275,276]. This may be observable in upcoming spectroscopic surveys such as DESI [277]. Ref. [278] showed that the octopole is a particularly clean probe of screening per se (as opposed to the background modification that screened theories also imply) because it is not degenerate with the difference in bias between the galaxy sub-populations.
The second probe, weak lensing, is the 2-point correlation function of cosmic shear, which emerges when galaxy shapes get distorted, their ellipticities increased and their magnitudes changed, due to light travelling though large scale structures in the Universe, from the source to the observer [279]. These ellipticities and magnitudes are correlated through the distribution of matter in the Universe and the expansion. Therefore they can provide very valuable information about the formation of structures from high redshifts until today. This angular correlation function can be expressed as where E(z) = H(z)/H 0 is the dimensionless Hubble function,Ŵ γ j (z) are window functions, or lensing kernels, that project the redshift distributions and the power spectrum into angular scales and finally P Φ+Ψ (k , z) is the Weyl power spectrum, which is related to the matter power spectrum P δδ by In this equation we can see clearly the observational signature of the Σ lensing function defined above in (157) and (154). We refer the reader again to [274] and the many references therein for details on the formulae of weak lensing. In figure 10 we show the non-linear matter power spectrum P (k, z) for ΛCDM (in light blue), K-mouflage (in green), JBD (in orange) and nDGP (in red) computed with privately modified versions of MGCAMB, hi class and EFTCAMB. The models and their fiducial values have been chosen to be close enough to ΛCDM , to be still allowed by observations, but far enough so that distinctive changes can be measured with next-generation surveys. The standard cosmological parameters set for this specific prediction are Ω m,0 = 0.315, Ω b,0 = 0.05, h = 0.674, ns = 0.966 and σ 8 = 0.8156. For JBD, the model parameter ω BD is set to ω BD = 800, while for nDGP the observational parameter is Ω rc = c 2 /(4r 2 5 H 2 0 ), where r 5 is the cross-over scale defined above in (61) and we set here Ω rc = 0.25. For K-mouflage the physical parameter is 2 = d ln A(a)/d ln a and it is related to the fifth force enacted by the scalar field, which comes from the conformal transformation of the metric (see [280] for more details). The prediction shown here is done for the case 2 = −0.04.
These distinctive features can be observed when taking the ratio to ΛCDM for the three cosmological models considered above. While at linear scales k < ∼ 0.07h/Mpc the models show only a slight change in amplitude compared to ΛCDM (with nDGP showing the largest amplitude increase of about 10%), it is clear that for small scales there are distinctive features at play that dampen the power spectrum. In the right panels of figure 10 we show the angular cosmic shear (weak lensing) power spectra in the 1,1 bin (lower redshifts) defined in (181) for all three screened models defined above. Also here, the ratio of the weak lensing C with respect to ΛCDM is shown in the lower panel. In this case, the very sharp features observed in the matter power spectrum get smoothed out by the projection along the line of sight and into angular multipoles.

E. Going beyond linear scales
At the linear level of perturbations in the matter density and the scalar field these equations above can be computed very efficiently using modified versions of Einstein-Boltzmann codes, in particular of CAMB 16 (Code for Anisotropies in the Microwave Background) (see [281]), which is written mainly in fortran, and CLASS 17 (see [282,283]), which is mainly written in the C programming language. Both of these codes come with user-friendly python wrappers. The most common modifications of these codes accounting for theories of modified gravity and dark energy are based on two types; the first one are codes in which generic parametrizations of the deviations of GR as in (153) to (156) are used. And the second one are codes in which specific modified gravity (MG) models or generic class of models are implemented and their full scalar field equations are solved, beyond the quasi-static approximation. From the first Lower left: Ratio to ΛCDM of the matter power spectra for the three cosmological models considered above. While at linear scales k < ∼ 0.07h/Mpc the models considered show only a slight change in amplitude compared to ΛCDM , at smaller scales there are some distinctive features, like shifts in the BAO peaks and damping of power at small scales. Upper right: Angular cosmic shear (weak lensing) power spectra for the 1,1 bin (lower redshifts) defined in (181) for the models mentioned above. Lower right: Ratio of the weak lensing C for screened modified gravity with respect to ΛCDM . The distinctive features observed in the matter power spectrum get smoothed out by the projection along the line of sight and into angular multipoles.
type, the two more common are ISitGR 18 (see [284,285], and MGCAMB 19 [286,287] and more recently a branch of CLASS, called QSA CLASS (see [288]). For the second type we will mention here the two most important ones, namely hi class 20 (see [289] and EFTCAMB 21 (see [290,291]).
Up to now, we have only developed the formalism to compute the perturbations of matter and the field at the linear level. However, in order to study correctly and accurately the power spectrum and compare it with observations of galaxy clustering and weak lensing, one must go beyond linear scales. For galaxy clustering, the region around k ≈ 0.1Mpc −1 , where Baryon Acoustic Oscillations (BAO) and redshift space distortions (RSD) are important, needs 18 https://labs.utdallas.edu/mishak/isitgr/ 19 https://github.com/sfu-cosmo/MGCAMB 20 http://miguelzuma.github.io/hi class public/ 21 http://eftcamb.org to be treated perturbatively, in order to make accurate predictions. This involves using either Eulerian or Lagrangian perturbation theory [292,293] and furthermore using resummation techniques to capture accurately the large scale contributions [294,295]. For smaller scales, formalisms such as the effective field theory of large scale structures are needed in order to take into account the UV divergences of the perturbative models [296]. For the models we are interested in here, there has been some recent work by [236,297], some new work on GridSPT by [298] and some more foundational work by [299,300].
To obtain meaningful constraints with weak lensing, the power spectrum needs to be calculated at even higher kvalues, for up to k ≈ 10Mpc −1 , which is only possible using N-body simulations, which capture the full evolution of the nonlinear gravitational dynamics. In scalar field models and especially in models that invoke screening mechanisms, these simulations are extremely computationally expensive and numerically complicated, since the nonlinear evolution of the field needs to be taken into account. Several interesting approaches have been taken in the literature such as COLA [301], Ramses [302], Ecosmog [303], MG-evolution [304], φ-enics (an interesting finite element method approach that can capture the nonlinear evolution of the scalar field and reproduce very accurately the Vainshtein screening) [305] and the simulation work on f (R) theories by several groups [306][307][308][309]. Since these simulations are time-consuming, faster approaches that allow for an efficient exploration of the parameter space would be extremely valuable and would be included into forecasts and Monte Carlo parameter estimation. Several approaches include fitting formulae based on simulations [310], emulators for f (R) theories [311,312] and hybrid approaches in which the halo model, perturbation theory and simulations are calibrated to create a model, such as REACT (see [313,314]). This code can compute predictions for nDGP and f (R) models which are roughly 5% accurate at scales k < ∼ 5h/Mpc.

F. Constraints on screened models with current data
In this section we will focus on the constraints on different screened scalar-field models with current observations from CMB, background expansion and large scale structure.

Constraints on f (R) models
From the CMB, constraints have been placed by the Planck collaboration on the f (R) model in terms of the Compton wavelength parameter, which is defined as and its value today B 0 is related to the fundamental parameter f R,0 . Indeed we have the relation where ω is the equation of state of the Universe and m f R is the mass of the scalaron (166). Notice that the denominator 1 + ω is very small and therefore B is less suppressed than the ratio H 2 /m 2 f R In the analysis of [315] the datasets used were Planck TT+lowP+BAO+SNIa + local H0 measurements (these last three observations are usually abbreviated as BSH), while CMB lensing was used to remove the degeneracy between B 0 and the optical depth τ . At the 95% confidence level, they found B 0 < 0.12 with Planck data alone and when BAO, weak lensing (WL) and RSD were added, a much more stringent bound of B 0 < 0.79 × 10 −4 was found, which forces the model to be very close to ΛCDM.
A very comprehensive, but by now relatively outdated, analysis by [316] using WMAP5 CMB data [317] and crosscorrelations of ISW with galaxy clustering data provided interesting bounds on the variations of the gravitational potentials on an interesting redshift range 0.1 < z < 1.5. For f (R) models that follow the same expansion of the universe as ΛCDM they obtained a bound of B 0 < 0.4 at the 95% confidence level (CL). In the analysis by [318] large scale structure data coming from WiggleZ, BAO (from 6dF, SDSS DR7 and BOSS DR9, see [253][254][255]) was combined with Planck-2013 CMB [319] data and WMAP polarization data [317] to find log 10 B 0 < −4.07 at the 95% CL. A more recent paper [320] uses the designer approach to f (R) and tests it with Planck and BAO data. In this designer approach one can fix the evolution of the background and then find the corresponding scalar field model that fits these constraints. With this the bound of B 0 < 0.006 (95%CL) for the designer models with w = −1 is obtained, and a bound of B 0 < 0.0045 for models with varying equation of state is reached, which was then constrained to be |w + 1| < 0.002 (95%CL). All these bounds imply that f (R) models cannot be self-accelerating and also if they are present, their background expansion will be very close to the one of ΛCDM according to observational bounds. This confirms the known results from gravitational tests in the solar system.

Constraints on nDGP models
The self-accelerating branch of DGP (sDGP) has been plagued with the presence of ghost fields, nevertheless it has been compared to observations, most recently in [321,322] where it was found, after using Planck temperature data, ISW and ISW-galaxy-cross-correlations, together with distance measurements that these models are much disfavoured compared to ΛCDM. The normal branch of DGP (nDPG) is non-self-accelerating, but it is still of interest since it shows clear deviations at scales important for structure formation. In [323] it was shown that the growth rate values estimated from the BOSS DR12 data [324] constrains the crossover scale r 5 of DGP gravity in the combination [r 5 H 0 ] −1 which has to be < 0.97 at the 2σ level, which amounts to r 5 > 3090Mpc/h, meaning that r 5 ∼ H −1 0 , therefore making this model very similar to GR within our Hubble horizon. Further tests of this model against simulations and large scale structure data have been performed in [325,326].

Constraints on Brans-Dicke theory
As mentioned previously, the most stringent constraint on JBD comes from solar system tests, where the Cassini mission put the bound of ω BD > 40, 000 (see [327,328]). However, under an efficient screening mechanism (invoking a specific potential), the theory could still depart considerably from GR at cosmological scales. In an analysis by [329] the authors used Planck [319], WMAP [317], SPT and ACT [330,331] data plus constraints on BBN to set bounds on the JBD parameter. They assumed the scalar field to have initial conditions such, that the gravitational constant would be the Newton constant today. With this they find ω BD > 692 at 99% C.L. When the scalar is free and varied as a parameter they find ω BD > 890, which amounts to 0.981 < G eff /G N < 1.285 at the 99% C.L. In a more recent analysis by [261], the authors used the combined data of the Planck CMB temperature, polarization, and lensing reconstruction, the Pantheon supernova distances, BOSS measurements of BAO, along with the joint 3 × 2pt dataset of cosmic shear, galaxy-galaxy lensing, and galaxy clustering from KiDS and 2dFLenS. They took into account perturbation theory and N-body calculations from COLA and RAMSES to compute the theoretical predictions for the power spectrum. They constrain the JBD coupling constant to be ω BD > 1540 at the 95% C.L. and the effective gravitational constant, G eff /G = 0.997 ± 0.029. They also found that the uncertainty in the gravitational theory alleviates the tension between KiDS, 2dFLenS and Planck to below 1σ and the tension [332] in the Hubble constant between Planck and the local measurements to 3σ. Despite these improvements, a careful model selection analysis, shows no substantial preference for JBD gravity relative to ΛCDM.

Constraints on Horndeski theories and beyond
For Horndeski models, there has been a great effort by the Planck collaboration to test the parametrized deviations of GR such as in (174) and (175) or in the α-formalism of [263]. However, in order to do so, certain conditions and restrictions on these parameters have to be met, given the relatively limited constraining power of current data. The code used in this case is the EFTCAMB code mentioned in section V E.
In the Planck 2015 modified gravity paper [315], the authors considered Horndeski models with α M = −α B , α T = α H = 0, and α K was fixed to a constant. This amounts to consider non-minimally coupled K-mouflage type models as in [263] with the only free function being α M . Additionally, the analysis used the ansatz, where α today M is a constant and p > 0 determines its backward time evolution. Furthermore, they relate the evolution of α M to a linear (p=1) and exponential (p > 1, varying free) parametrization [315]. Using the Planck TT+TE+EE+BSH data set combination (BSH standing again for BAO, SN and local Hubble constraints) they find α today M < 0.043 (95% confidence level) for the linear case and α today M < 0.062 and p = 0.92 0.53 0.24 (95% confidence level) for the exponential case. ΛCDM is recovered for p = 1 and α today M = 0, therefore placing relatively strong limits on possible deviations of Einstein's GR.
As we discussed above, the gravitational wave event GW170817 constrained the Horndeski theory to be effectively composed only of Brans-Dicke models and cubic Galileons, and the latter are effectively ruled out by ISW observations. This limits then the interest on an overall analysis of Horndeski models in general. However in [333] the authors analysed Horndeski models that still can have non-trivial modifications to GR, possible at the level of linear perturbations and they confirmed the conjecture by [266] that (Σ − 1)(µ − 1) ≥ 0 for surviving models.
As an extension beyond this review, DHOST models, as mentioned above, can also provide an interesting phenomenology and are able to evade certain constraints affecting the Horndeski theories. [334,335] studied DHOST models that present self-acceleration and [336], among others, have studied the astrophysical signatures of these models. However, their theoretical modelling has not been implemented yet in computational tools capable of analyzing the full Planck CMB dataset. Finally, [337] performed a cosmological constraint analysis, assuming the form α i = α i,0 a κ on these surviving Horndeski models and using Planck and BICEP2/Keck [338] CMB data and galaxy clustering data from SDSS and BOSS, they found that when setting the kineticity to the following value α K = 0.1a 3 , the α M,0 parameter has an upper limit of 0.38 when α B,0 = 0 and 0.41 when α B,0 = 0 at the 95% C.L. More importantly, they conclude that the effects of Horndeski theory on primordial B-modes (which at the time were expected to be measured accurately by BICEP/KECK2) are constrained by CMB and LSS data to be insignificant at the 95% C.L. However, they draw the attention to the fact that the assumptions on some parameters, for example the assumed form of the kineticity, have major and dramatic effects on these results. In conclusion, the theory space of Horndeski models has been mostly ruled out by measurements of the ISW effect and the combination of CMB and large scale structure, when considering the gravitational wave event GW170817 and its electromagnetic counterpart GRB170817A . On the other hand, beyond Horndeski theories, such as DHOST, seem promising, but computational tools required to do a proper cosmological analysis are not available yet, so the models can only be constrained by astrophysical observations so far.

VI. CONCLUSIONS AND PERSPECTIVES
Scalar-tensor theories are among the most generic and plausible extensions to ΛCDM, with potential relevance to much of astrophysics and cosmology. They must be screened to pass solar system tests of fifth forces. In this review we have presented the most common screened modified gravity mechanisms and introduced them using an effective field theory point of view. The effective point of view is taken by first selecting a background which could be cosmological, astrophysical or local in the solar system. The coefficients of the different operators depend on the environment. This is a feature of all the screening mechanisms-physics is dependent on the distribution of matter-and gives them relevance to various different types of environment on a range of scales.
The screening mechanisms can be divided in two categories. The non-derivative mechanisms consist of the chameleon and Damour-Polyakov cases. The derivative ones are the K-mouflage and Vainshtein scenarios. The latter lead to scale-independent modifications of gravity on large scales. For models with derivative screening and having effects on large cosmological scales, the effects on smaller scales are reduced due to the strong reduction of fifth force effects inside the K-mouflage and Vainshtein radii. Nonetheless, the force laws on short scales in these scenarios deviates from 1/r 2 and leads to effects such as the advance of the perisastron of planets and effective violation of the strong equivalence principle in galaxies, both of which afford tight constraints. There is still some capability for ground-based experiments to test Vainshtein-screened theories however [90]. The time dependence induced by the cosmological evolution is not screened in K-mouflage and Vainshtein screened models, which also leads to tight bounds coming from solar system tests of gravitation.
The chameleon and Damour-Polyakov mechanisms, on the other hand, have effects on scales all the way from the laboratory to the cosmos, and must be taken on a case-by-case basis for each experimental set up and astrophysical observation. This makes the comparison between the short and large scale physics richer, and leads to more complementarity between astrophysical and laboratory tests. For the symmetron, an experiment with a length scale between objects d typically best constrains theories with mass parameter µ ≈ d −1 . If the mass were larger, then the scalar force between objects would be exponentially suppressed (as in (114)), while if it were smaller the field would remain near φ = 0 where it is effectively decoupled from matter. It is therefore desirable to employ a range of tests across as many length scales as possible. There is a notable exception to this general rule: if the ambient matter density between objects is of order the symmetry-breaking value ρ amb ≈ µ 2 M 2 , then the symmetron is essentially massless. This enables even long-ranged experiments to test symmetron theories with µ d −1 at that particular value of M . The chameleon does not have a fixed mass parameter and hence there is more overlap between various experiments' capabilities to test the theory. Here the differentiating feature tends to be when objects of a particular size become screened. If a given experiment's source and/or test mass is screened, then the experiment's capability to test the theory is strongly suppressed. Small values of the chameleon parameters {M, Λ} correspond to even microscopic objects being screened, so only small-scale experiments are able to test that region of parameter space. One can observe this general trend in Fig. 1: the bottom-left corner is constrained by particle physics experiments, the middle region by atomic-scale experiments, and the upper-right region by experiments employing macroscopic test masses like a torsion balance. This trend continues with astrophysical tests constraining the region further above and to the right of the parameter space illustrated in the figure.
We have seen that, although screening mechanisms are easily classified, empirical testing is most often performed at the model level. Some of these models are archetypal, such as the f (R) models of the Hu-Sawicki type for chameleons, the symmetrons for Damour-Polyakov, and the nDGP model for Vainshtein. For K-mouflage, there is no such template although specific models such as the "arctan" are promising because they pass the solar system tests. On cosmological scales it is easier to test many theories at once, e.g. through the effective field theory of dark energy. Unfortunately, the link between the large scales and the small scales where screening must be taken into account is then lost. This is also a problem on cosmological scales where non-linear effects must be taken into account for weak lensing for instance and bridging the gap beyond perturbation theory and highly non-linear scales necessitates tools such as N-body simulations, which may be computationally expensive. A parameterisation of screening mechanisms valid from laboratory scales to the cosmological horizon would certainly be welcome. In the realm of non-derivative screenings, such a parameterisation exist and depends only on the mass and coupling dependence as a function of the scale factor the Universe. This allows to reconstruct the whole dynamics of the models, on all scales [53,87]. The same type of parameterisation exists for K-mouflage where the coupling function A(a) and the screening factor Z(a) are enough to reconstruct the whole dynamics too [339]. For Vainshtein and generalised cubic models defined by the function G 3 , this should also be the case although it has not yet been developed.
Fortunately, the space of theories which still need to be tested has drastically shrunk in the last few years. The models with the Vainshtein mechanisms and some influence on large scales are restricted to theories parameterised by one function G 3 which must be non-trivial as the simplest case, the cubic Galileon, has been excluded by observations of the Integrated Sachs Wolfe effect. Quartic and quintic galileons are powerfully constrained by GW170817, the observation of a gravitational wave event with near-simultaneous optical counterpart. Of course, theories with the Vainshtein property and no link with the cosmological physics of late-time acceleration of the Universe's expansion are fine although the parameter space is restricted by galaxy-scale tests. On the thin-shell-screening side, wide regions of chameleon and symmetron parameter space are ruled out by laboratory probes, and a largely complementary part by astrophysical tests involving stars and galaxies. The n = 1 Hu-Sawicki theory-the workhorse chameleon-screened model for over a decade-is now constrained by galaxy morphology to the level f R0 < 1.4 × 10 −8 [184], such that it can no longer have appreciable astrophysical or cosmological effects. The phenomenological description of DHOST models is less developed, and it would be interesting to see whether and how these models could answer some of the pressing questions in cosmology such as the origin of the acceleration.
Future observations on cosmological scales from upcoming surveys such as Euclid will certainly provide a host of new results on screened models such as K-mouflage or nDGP. Only recently has it been realised that galactic scales afford strong probes of screening, and many more tests will likely be developed in the future. In the solar system, future satellite tests [340], which will test the equivalence principle down to a level of 10 −17 in the Eötvos parameter, should also constrain screening mechanisms of the non-derivative type [341,342]. Finally, laboratory experiments ranging from the search for new interaction with Casimir configurations to atom interferometry should also provide new possibilities for the detection of screened modified gravity. While we have focused in this review on the relevance of screened scalar fields to the physics of dark energy, it may also be relevant to the other missing pillar of ΛCDM, dark matter. This is a key target for many upcoming astrophysical and cosmological surveys. Much less is known about screening is this regard, although fifth forces are clearly degenerate with dark matter in determining diverse objects' dynamics.
In conclusion, screening is a crucial ingredient in the physics of light scalar fields. Testing it with the next generation of experiments and observations may well lead to surprises and new discoveries.