Testing Screened Modified Gravity

Abstract

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.

1. Introduction: Why Light Scalars?

Light scalar fields are mainly motivated by two unexplained phenomena, the existence of astrophysical effects associated with dark matter [1,2] and the apparent acceleration of the expansion of the Universe [3,4]. In both cases, traditional explanations exist. Dark matter could be a Beyond the Standard Model (BSM) particle (or particles) with weak interactions with ordinary matter (WIMPs) [5]. The acceleration of the Universe could be the result of the pervading presence of a constant energy density, often understood as a pure cosmological constant term [6] in the Lagrangian governing the dynamics of the Universe on large scales, whose origin remains mysterious [7,8]. Lately, this standard scenario, at least on the dark matter side, has been challenged due to the lack of direct evidence in favour of WIMPS at accelerators or in large experiments dedicated to their search (Xenon1T and similar experiments) [9,10,11,12,13,14,15,16,17,18,19,20]. In this context, the axion or its related cousins, the ALP’s (Axion-Like Particles) have come back to the fore [21]. More generally, (pseudo-)scalars could play the role of dark matter thanks to the misalignment mechanism, i.e., they behave as oscillating fields, as long as their mass m is low, typically $m\lesssim 1$ eV [22,23,24].

On the late acceleration side, the cosmological constant is certainly a strong contender, albeit a very frustrating one. The complete absence of dynamics required by a constant vacuum energy is at odds with what we know about another phase of acceleration, this time in the very early Universe, i.e., inflation [25,26,27,28]. This is the leading contender to unravel a host of conundrums, from the apparent isotropy of the Cosmic Microwave Background (CMB) to the generation of primordial fluctuations. The satellite experiment Planck [29] has taught us that the measured non-flatness of the primordial power spectrum of fluctuations is most likely due to a scalar rolling down its effective potential. This and earlier results have prompted decades of research on the possible origin of the late acceleration of the Universe.

In most of these models, scalar fields play a leading role and appear to be very light on cosmological scales, with masses sometimes as low as the Hubble rate now, ${10}^{-33}$ eV [30,31]. Quantum mechanical considerations and in particular the presence of gravitational interactions always generate interactions between these scalars and matter. The existence of such couplings is even de rigueur from an effective field theory point of view (in the absence of any symmetry guaranteeing their absence) [32]. Immediately this leads to a theoretical dead end, however, as natural couplings to matter would inevitably imply strong violations of the known bounds on the existence of fifth forces in the solar system, e.g., from the Cassini probe [33]. As a typical example, $f\left(R\right)$ models [34] with a normalised coupling to matter of $\beta =1/\sqrt{6}$ belong to this category of models, which would be excluded if non-linearities did not come to the rescue [35].

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 $\Lambda$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 (Section 2) before discussing tests of screening on laboratory (Section 3), astrophysical (Section 4) and cosmological (Section 5) scales. Section 6 summarises and discusses the complementarity (and differences) between these classes of tests.

2. Screening Light Scalars

2.1. 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 overdensity. 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 $\varphi$, although the analysis is straightforwardly generalised to multiple scalars. The scalar’s background configuration is denoted by ${\varphi }_0$ and the induced perturbation of the scalar field due to the perturbation by an overdensity will be denoted by $\phi \equiv \varphi -{\varphi }_0$. At the 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]

$${\mathcal{L}}_2=\frac{Z\left({\varphi }_0\right)}{2}{\left({\partial }_{\mu }\phi \right)}^2+\frac{m_{\varphi }^2\left({\varphi }_0\right)}{2}{\phi }^2-\delta g_{\mu \nu }\delta T^{\mu \nu }+\cdots$$

(1)

at the second-order in the scalar perturbation $\phi$ and the matter perturbation $\delta T_{\mu \nu }$. The latter is the perturbed energy-momentum tensor of matter compared to the background. In this Lagrangian, matter is minimally coupled to the perturbed Jordan metric $\delta g_{\mu \nu }$ and the Jordan frame energy-momentum tensor is therefore conserved ${\partial }_{\mu }\delta T^{\mu \nu }=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\left({\varphi }_0\right)$, i.e., the mass of the scalar field. The second is the wave function normalisation $Z\left({\varphi }_0\right)$. The third is the composite metric $\delta g_{\mu \nu }$, which is not the local metric of space-time but the leading 2-tensor mediating the interactions between the scalar field and matter. This Jordan metric can be expanded as

$$\delta g_{\mu \nu }=\frac{\beta \left({\varphi }_0\right)}{m_{\mathrm{Pl}}}\phi {\eta }_{\mu \nu }-\gamma \left({\varphi }_0\right){\partial }_{\mu }{\partial }_{\nu }\phi +\delta \left({\varphi }_0\right){\partial }_{\mu }\phi {\partial }_{\nu }\phi +\cdots .$$

(2)

At the 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 $\beta \left({\varphi }_0\right)$1. One can also introduce a term in second derivatives of $\phi$, 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 $\phi$ squared and a coupling constant of dimension minus four. These two terms can be seen as disformal interactions [47].

The equations of motion for $\phi$ are given by

$${\partial }_{\mu }\left(Z\left({\varphi }_0\right){\partial }^{\mu }\phi \right)-m^2\left({\varphi }_0\right)\phi -2{\partial }_{\mu }\left(\delta \left({\varphi }_0\right){\partial }_{\nu }\phi \right)\delta T^{\mu \nu }=-\frac{\beta \left({\varphi }_0\right)}{m_{\mathrm{Pl}}}\delta T+{\partial }_{\mu }{\partial }_{\nu }\left(\gamma \left({\varphi }_0\right)\right)\delta T^{\mu \nu }$$

(3)

where $\delta T\equiv \delta T_{\nu }^{\mu }$, and we have used the conservation of matter ${\partial }_{\mu }\delta T^{\mu \nu }=0$. This equation will allow us to describe the different screening mechanisms.

2.2. Modified Gravity

Let us now specialise the Klein–Gordon equation to experimental or observational cases where $\delta T^{00}=\rho$ is a static matter overdensity locally and the background is static too. This corresponds to a typical experimental situation where overdensities are the test masses of a gravitational experiment. In this case, we can focus on the case where ${\varphi }_0$ can be considered to be locally constant. As a result, we have

$$K^{\mu \nu }\left({\varphi }_0\right){\partial }_{\mu }{\partial }_{\nu }\phi -m^2\left({\varphi }_0\right)\phi =-\frac{\beta \left({\varphi }_0\right)}{m_{\mathrm{Pl}}}\delta T.$$

(4)

The kinetic terms are modified by the tensor

$$K^{\mu \nu }\left({\varphi }_0\right)=Z\left({\varphi }_0\right){\eta }^{\mu \nu }-2\delta \left({\varphi }_0\right)\delta T^{\mu \nu }.$$

(5)

When the overdensities are static, the disformal term in $K^{\mu \nu }$, which depends on the matter energy momentum tensor, does not contribute and we have $K^{ij}\left({\varphi }_0\right)\simeq Z\left({\varphi }_0\right){\delta }^{ij}$ leading to the modified Yukawa equation

$$\Delta \phi -\frac{m^2\left({\varphi }_0\right)}{Z\left({\varphi }_0\right)}\phi =\frac{\hat{\beta }\left({\varphi }_0\right)}{m_{\mathrm{Pl}}}\rho .$$

(6)

where $\hat{\beta }\left({\varphi }_0\right)=\frac{\beta \left({\varphi }_0\right)}{Z\left({\varphi }_0\right)}$. For nearly massless fields we can neglect the mass term within the Compton wavelength of size $m^{-1}\left({\varphi }_0\right)$, which is assumed to be much larger than the experimental setting. In this case, the Yukawa equation becomes a Poisson equation

$$\Delta \phi =\frac{\hat{\beta }\left({\varphi }_0\right)}{m_{\mathrm{Pl}}}\rho .$$

(7)

As a result, the scalar field behaves similarly to the Newtonian potential, and the matter interacts with the effective Newtonian potential2

$$\Phi ={\Phi }_N+\frac{\beta \left({\varphi }_0\right)}{m_{\mathrm{Pl}}}\phi =\left(1+\frac{2{\beta }^2\left({\varphi }_0\right)}{Z\left({\varphi }_0\right)}\right){\Phi }_N$$

(8)

i.e., gravity is modified with an increase in Newton’s constant by

$$G_{\mathrm{eff}}=\left(1+\frac{2{\beta }^2\left({\varphi }_0\right)}{Z\left({\varphi }_0\right)}\right)G_N.$$

(9)

Notice that the scalar field does not couple to photons as $\delta T=0$; hence, matter particles deviated with a larger Newtonian interaction than photons, $G_{\mathrm{eff}}\ge G_N$. As a result, the modification of $G_N$ into $G_{\mathrm{eff}}$ is not just a global rescaling and gravitational physics is 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:

$$G_{\mathrm{eff}}=\left(1+\frac{2{\beta }^2\left({\varphi }_0\right)}{Z\left({\varphi }_0\right)}{e}^{-m\left({\varphi }_0\right)r}\right)G_N,$$

(10)

where r is the distance to the object sourcing the field. This equation allows us to classify the screening mechanisms.

2.3. The Non-Derivative Screening Mechanisms: Chameleon and Damour–Polyakov

The first mechanism corresponds to an environment-dependent mass $m\left({\varphi }_0\right)$. 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_{\varphi }\left({\varphi }_0\right)r}$, where r is the distance from the mass element. This implies that only a $thinshell$ of mass $\Delta 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

$$\beta \left({\varphi }_0\right)\to \frac{\Delta M}{M}\beta \left({\varphi }_0\right)$$

(11)

where M is the mass of the object. As long as $\Delta M/M\ll 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 $\beta \left({\varphi }_0\right)=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 Equation (7), we obtain

$$\frac{\phi }{m_{\mathrm{Pl}}}=\frac{2\beta \left({\varphi }_0\right)}{Z\left({\varphi }_0\right)}{\Phi }_N.$$

(12)

Let us first consider the case of a normalised scalar field with $Z\left({\varphi }_0\right)=1$. The scalar field is screened when its response to the presence of an overdensity is suppressed compared to the Newtonian case. This requires that

$$\frac{\left|\phi \right|}{2m_{\mathrm{Pl}}\left|{\Phi }_N\right|}\le \beta \left({\varphi }_{\mathrm{out}}\right)$$

(13)

where $\phi ={\varphi }_{\mathrm{in}}-{\varphi }_0$ is the variation of the scalar field inside the dense object. Here ${\Phi }_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 that are sufficiently dense, the field ${\varphi }_{\mathrm{in}}$ nearly vanishes and $\phi \simeq -{\varphi }_0$ only depends on the environment. As a result, for such dense objects, screening occurs when $|{\Phi }_N|\ge \frac{{\varphi }_0}{2m_{\mathrm{Pl}}\beta \left({\varphi }_{\mathrm{out}}\right)}$, 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

$${\beta }_A=\frac{|{\phi }_A|}{2m_{\mathrm{Pl}}\left|{\Phi }_N^A\right|}$$

(14)

for a screened object labelled by A is the scalar charge of this object3, 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 $\beta \left({\varphi }_0\right)$.

2.4. The Derivative Screening Mechanisms: K-Mouflage and Vainshtein

The third case, in fact, covers two mechanisms. If locally in a region of space, the normalisation factor

$$Z\left({\varphi }_0\right)\gg 1$$

(15)

then obviously the effective coupling $\hat{\beta }\left({\varphi }_0\right)\ll \beta \left({\varphi }_0\right)$ and gravitational tests can be evaded. Notice that we define screening as reducing the effective coupling. This case covers the K-mouflage4 and Vainshtein mechanisms.

The normalisation factor is a constant at the leading order. Going beyond leading order, i.e., including higher order operators in the effective field theory, $Z\left({\varphi }_0\right)$ can be expanded in a power series

$$Z\left({\varphi }_0\right)=1+a\left({\varphi }_0\right)r_c^2\frac{\square \phi }{m_{\mathrm{Pl}}}+b\left({\varphi }_0\right)\frac{{(\partial \phi )}^2}{{\Lambda }^4}+c\left(\varphi \right)\frac{{\square }^2\phi }{{\Lambda }^5}+\cdots$$

(16)

where $r_c$ is a cross over scale and has the dimension of length and $m_{\mathrm{Pl}}$ is the Planck scale. The scale $\Lambda$ 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.

2.4.1. K-Mouflage

The K-mouflage screening mechanism [46,54,55] is at play when $Z\left({\varphi }_0\right)\ge 1$ and the term in ${(\partial \phi )}^2/{\Lambda }^4$ dominates in (16), i.e.,

$$|\overrightarrow{\nabla }\phi |\ge {\Lambda }^2$$

(17)

and therefore, the Newtonian potential must satisfy

$$|\overrightarrow{\nabla }{\Phi }_N|\ge \frac{{\Lambda }^2}{2\beta \left({\varphi }_0\right)m_{\mathrm{Pl}}}.$$

(18)

Hence, K-mouflage screening occurs where the gravitational acceleration ${\overrightarrow{a}}_N=-\overrightarrow{\nabla }{\Phi }_N$ is large enough. Let us consider two typical situations. First, the Newtonian potential of a pointlike source of mass M has a gradient satisfying (18) inside a radius $R_K$

$$R_K={\left(\frac{\beta \left({\varphi }_0\right)M}{4\pi m_{\mathrm{Pl}}{\Lambda }^2}\right)}^{1/2}.$$

(19)

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 sourced by overdensities $\delta \rho$ compared to the background energy density $\overline{\rho }$. In this case, screening takes place for wave-numbers k such that

$$\beta \left({\varphi }_0\right)\delta \le \frac{{\Lambda }^2{m}_{\mathrm{Pl}}k}{\overline{\rho }}$$

(20)

where $\delta \equiv \frac{\delta \rho }{\overline{\rho }}$. In particular for models motivated by dark energy ${\Lambda }^4\simeq m_{\mathrm{Pl}}^2{H}_0^2$ screening occurs on scales such that $k/H_0\gtrsim \beta \left({\varphi }_0\right)\delta$, i.e., large scale structures such as galaxy clusters are not screened as they satisfy $k/H_0\lesssim 1$ [56,57].

2.4.2. Vainshtein

The Vainshtein mechanism [58,59] follows the same pattern as K-mouflage. The main difference is that now the dominant term in $Z\left({\varphi }_0\right)$, i.e., (16), is given by the $\square \phi$ term. This implies that

$$\Delta {\Phi }_N\ge \frac{1}{2\beta \left({\varphi }_0\right)r_c^2},$$

(21)

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

$$R_V={\left(\frac{3\beta \left({\varphi }_0\right)r_c^2}{4\pi m_{\mathrm{Pl}}^{2}M}\right)}^{1/3}.$$

(22)

Notice the power $1/3$ compared to the $1/2$ in the K-mouflage case. Similarly, on large scales where the density contrast is $\delta$, the scalar field is screened for wave numbers such that

$$\delta \ge \frac{1}{3{\Omega }_m\beta \left({\varphi }_0\right)}$$

(23)

where ${\Omega }_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 $\delta \gtrsim 1$ as long as $\beta \left({\varphi }_0\right)\gtrsim 1$.

Finally, let us consider the case where the term in ${\square }^2\phi$ dominates in (16). This corresponds to6

$${\Delta }^2{\Phi }_N\ge \frac{{\Lambda }^5}{2\beta \left({\varphi }_0\right)m_{\mathrm{Pl}}}$$

(24)

For a point source, the transition happens at the radius

$$R_{\mathrm{MG}}={\left(\frac{5M}{4\pi m_{\mathrm{Pl}}{\Lambda }^5}\right)}^{1/5}.$$

(25)

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]

$${\Lambda }^5=m_{\mathrm{Pl}}{m}_G^4$$

(26)

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 than the chameleon and Damour–Polyakov ones under radiative effects.

2.5. Screening Criteria: The Newtonian Potential and Its Derivatives

Finally, let us notice that the screening mechanisms can be classified by inequalities of the type

$$\left({\nabla }^k\right){\Phi }_N\gtrsim C$$

(27)

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 space-time 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\left({\varphi }_0\right)$ 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\le 2$ are forbidden due to the diffeomorphism invariance of the theory, see the discussion in Section 2.7.1.

2.6. 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 Section 4.2 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,

$$\Delta \phi =\ddot{\gamma }\left({\varphi }_0\right)\rho .$$

(28)

As a result, we see that a scalar charge is induced by the cosmological evolution of the background [64,65]

$${\beta }_{\mathrm{ind}}=m_{\mathrm{Pl}}({\gamma }_2{\dot{\varphi }}_0^2+{\gamma }_1{\ddot{\varphi }}_0)$$

(29)

where ${\gamma }_1=\frac{d\gamma }{d\varphi }$ and ${\gamma }_2=\frac{d^2\gamma }{d{\varphi }^2}$. This is particularly relevant to black holes solutions with a linear dependence in time ${\dot{\varphi }}_0=q$. In this case, the induced charge is strictly constant

$${\beta }_{\mathrm{ind}}={\gamma }_2{m}_{\mathrm{Pl}}{q}^2$$

(30)

which could lead to interesting phenomena in binary systems.

2.7. Examples of Screened Models

2.7.1. 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^{\mu \nu }$, the response of the weak gravitational field $h_{\mu \nu }=g_{\mu \nu }-{\eta }_{\mu \nu }$ is given in momentum space by 7

$$h_{\mu \nu }\left(p^{\lambda }\right)=\frac{16\pi G_N}{p^2}\left(T^{\mu \nu }-{\eta }_{\mu \nu }\frac{T}{2}\right)$$

(31)

where two features are important. The first is that $1/p^2=1/p^{\lambda }{p}_{\lambda }$ 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

$${\mathcal{L}}_{\mathrm{FP}}=-\frac{m_G^2{m}_{\mathrm{Pl}}^2}{8}(h_{\mu \nu }{h}^{\mu \nu }-h^2).$$

(32)

The tensorial structure of the Fierz–Pauli mass term guarantees the absence of ghosts in a flat background. The response to matter becomes

$$h_{\mu \nu }\left(p^{\lambda }\right)=\frac{16\pi G_N}{p^2+m_G^2}\left(T^{\mu \nu }-{\eta }_{\mu \nu }\frac{T}{3}\right).$$

(33)

The factor in $1/(p^2+m_G^2)$ is the propagator of a massive field of mass $m_G$. More surprising is the change $1/2\to 1/3$ in the tensorial expression. In particular, in the limit $m_G\to 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 unravelled as follows. Writing

$$h_{\mu \nu }={\overline{h}}_{\mu \nu }+\frac{\beta }{m_{\mathrm{Pl}}}\phi {\eta }_{\mu \nu }$$

(34)

where

$${\overline{h}}_{\mu \nu }\left(p^{\lambda }\right)=\frac{16\pi G_N}{p^2}(T^{\mu \nu }-{\eta }_{\mu \nu }\frac{T}{2})$$

(35)

and

$$\phi =\frac{\beta }{m_{\mathrm{Pl}}}\frac{T}{p^2+m_g^2}$$

(36)

corresponding to a scalar satisfying

$$\square \phi -m_G^2\phi =-\frac{\beta }{m_{\mathrm{Pl}}}T,$$

(37)

we find that (33) is satisfied provided that

$$\beta =\frac{1}{\sqrt{6}}.$$

(38)

Hence, we have decomposed the massive graviton into a helicity two part ${\overline{h}}_{\mu \nu }$ and a scalar part $\phi$ coupled to matter with a scalar charge $\beta =1/\sqrt{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\le {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]

$${\mathcal{L}}_3\sim \frac{{(\square \phi )}^3}{{\Lambda }^5}$$

(39)

where ${\Lambda }^5=m_{\mathrm{Pl}}{m}_G^4$. The structure in $\square \phi$ follows from the symmetry $\phi \to \phi +{\lambda }_{\mu }{x}^{\mu }$, which can be absorbed into a diffeomorphism ${\overline{h}}_{\mu \nu }\to {\overline{h}}_{\mu \nu }+{\partial }_{\mu }{\xi }_{\nu }+{\partial }_{\nu }{\xi }_{\mu }$ where ${\xi }_{\mu }=\frac{\beta }{4m_{\mathrm{Pl}}}(-2(\lambda .x)x_{\mu }+{\lambda }_{\mu }{x}^2)$. The Klein–Gordon equation is modified by terms in $\square \left({(\square \phi )}^2\right)$. As a result, the normalisation factor is dominated by $Z\left({\varphi }_0\right)\sim {\square }^2\phi /{\Lambda }^5$ as mentioned in the previous section. This leads to the Vainshtein mechanism inside $R_{\mathrm{MG}}$, which allows massive gravity to evade gravitational tests in the solar system, for instance.

2.7.2. Cubic Galileon Models

The cubic Galileon models [59] provide an example of a Vainshtein mechanism with the $1/5$ power instead of the $1/3$. They are defined by the Lagrangian

$${\mathcal{L}}_{\mathrm{Gal}3}=\frac{1}{2}{(\partial \varphi )}^2+\frac{{(\partial \varphi )}^2\square \varphi }{{\Lambda }^3}.$$

(40)

The normalisation factor for the kinetic terms involves $\square \varphi$ 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 2.8 for more details.

2.7.3. Quartic K-Mouflage

The simplest example of K-mouflage model is provided by the Lagrangian [64]

$${\mathcal{L}}_{\mathrm{KM}4}=\frac{1}{2}{(\partial \varphi )}^2-\frac{{(\partial \varphi )}^4}{{\Lambda }^4}$$

(41)

which is associated with a normalisation factor containing a term in ${(\partial \varphi )}^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 2.9 for more details.

2.7.4. Ratra–Peebles and f(R) Chameleons

Chameleons belong to a type of scalar-tensor theory [68] specified entirely by two functions of the field. The first one is the interaction potential $V\left(\varphi \right)$ and the second one is the coupling function $A\left(\varphi \right)$. The dynamics are driven by the effective potential [48,49]

$$V_{\mathrm{eff}}\left(\varphi \right)=V\left(\varphi \right)+(A\left(\varphi \right)-1)\rho$$

(42)

where $\rho$ is the conserved matter density. When the effective potential has a minimum ${\varphi }_0\equiv \varphi \left(\rho \right)$, its second derivative defines the mass of the chameleon

$$m^2\left(\rho \right)=\frac{d^2{V}_{\mathrm{eff}}}{d{\varphi }^2}{|}_{\varphi \left(\rho \right)}$$

(43)

Cosmologically, the chameleon minimum of the effective potential is an attractor when $m\left(\rho \right)\gg 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 2.9. A typical example of chameleon theory is provided by [48,49]

$$V\left(\varphi \right)=\frac{{\Lambda }^{n+4}}{{\varphi }^n},A\left(\varphi \right)=e^{\beta \varphi /m_{\mathrm{Pl}}}$$

(44)

associated with a constant coupling constant $\beta$. More generally, the coupling becomes density dependent as

$$\beta \left(\rho \right)=m_{\mathrm{Pl}}\frac{dlnA}{d\varphi }{|}_{\varphi \left(\rho \right)}.$$

(45)

Chameleons with $n=1$ are extremely well constrained by laboratory experiments, see Section 3.5.

Surprisingly, models of modified gravity defined by the Lagrangian [34]

$${\mathcal{L}}_{f\left(R\right)}=\frac{\sqrt{-g}f\left(R\right)}{16\pi G_N}$$

(46)

which is a function of the Ricci scalar R, can be transformed into a scalar-tensor setting. First of all, the field equations of $f\left(R\right)$ gravity can be obtained after a variation of the Lagrangian (46) with respect to the metric $g_{\mu \nu }$ and they read

$$f_R{R}_{\mu \nu }-\frac{1}{2}f{g}_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }{f}_R+g_{\mu \nu }\square f_R=8\pi G{T}_{\mu \nu },$$

(47)

where $f_R\equiv df\left(R\right)/dR$ and □ is the d’Alembertian operator. These equations naturally reduce to Einstein’s field equations when $f\left(R\right)=R$. This theory can be mapped to a scalar field theory via

$$\frac{df}{dR}=e^{-2\beta \varphi {/}_{{\mathrm{m}}_{\mathrm{Pl}}}}$$

(48)

where $\beta =\frac{1}{\sqrt{6}}$. The coupling function is given by the exponential

$$A\left(\varphi \right)=e^{\beta \varphi /m_{\mathrm{Pl}}}$$

(49)

leading to the same coupling to matter as massive gravity. Contrary to the massive gravity case, $f\left(R\right)$ models evade solar system tests of gravity thanks to the chameleon mechanism when the potential

$$V\left(\varphi \right)=\frac{m_{\mathrm{Pl}}^2}{2}\frac{R\frac{df}{dR}-R}{{\left(\frac{df}{dR}\right)}^2}$$

(50)

is appropriately chosen.

A popular model has been proposed by Hu–Sawicki [69] and reads

$$f\left(R\right)=-2\Lambda -\frac{f_{R_0}{c}^2}{n}\frac{R_0^{n+1}}{R^n},$$

(51)

which involves two parameters, the exponent (positive definite) $n>0$ and the normalisation $f_{R_0}$, which is constrained to be $|f_{R_0}|\lesssim {10}^{-6}$ by the requirement that the solar system is screened [69] (see Section 2.9). On large scales, structures are screened for which

$${\Phi }_N\gtrsim \chi \equiv \frac{3}{2}\left|f_{R_0}\right|$$

(52)

for the $n=1$ Hu–Sawicki model, where $\chi$ is the “self-screening parameter”. This follows directly from the fact that

$$f_{R_0}=-\frac{2\beta \delta \varphi }{m_{\mathrm{Pl}}}$$

(53)

where $\delta \varphi$ 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=\frac{\left|\delta \varphi \right|}{2\beta m_{\mathrm{Pl}}{\Phi }_N}$ must be less than $\beta =1/\sqrt{6}$—gives a useful criterion for identifying unscreened objects (see Section 4).

2.7.5. $f\left(R\right)$ and Brans-Dicke

The $f\left(R\right)$ models can be written as a scalar-tensor theory of the Brans–Dicke type. The first step is to replace the $f\left(R\right)$ Lagrangian density by

$$\mathcal{L}=\sqrt{-g}[f\left(\lambda \right)+(R-\lambda )\frac{\mathrm{d}f\left(\lambda \right)}{\mathrm{d}\lambda }]$$

(54)

which reduces to the original model by solving for $\lambda$. Then an auxiliary field

$$\psi \equiv \frac{\mathrm{d}f\left(\lambda \right)}{\mathrm{d}\lambda }$$

(55)

can be introduced, together with the potential $V\left(\psi \right)=m_{\mathrm{Pl}}^2(f\left(\lambda \left(\psi \right)\right)-\lambda \left(\psi \right)\psi )/2$, which corresponds to the Legendre transform of the function $f\left(\lambda \right)$. After replacing back into the original action, one recovers a scalar field action for $\psi$ in the Jordan frame that reads

$$S=\int {\mathrm{d}}^{4}x\sqrt{-g}\left[\psi R-\frac{{\omega }_{BD}\left(\psi \right)}{\psi }{\nabla }_{\mu }\psi {\nabla }^{\mu }\psi -V\left(\psi \right)\right]+16\pi G\int {\mathrm{d}}^{4}x\sqrt{-\tilde{g}}{\mathcal{L}}_m({\psi }_m^{\left(i\right)},g_{\mu \nu }).$$

(56)

This theory corresponds to the well-known generalised Jordan–Fierz–Brans–Dicke [70] theories with ${\omega }_{BD}=0$. When the ${\omega }_{BD}$ parameter is non-vanishing and a constant, this reduces to the popular Jordan–Brans–Dicke theory. Exact solutions of these theories have been tested against observations of the Solar System [33,71], and the Cassini mission sets the constraint ${\omega }_{BD}>40,000$, so that JBD has to be very close to GR. This bound is a reformulation of (88), see Section 2.9 for more details. After going to the Einstein frame, the theory must be a scalar-tensor with the Chameleon or Damour–Polyakov mechanisms in order to evade the gravitational tests in the solar system.

2.7.6. The Symmetron

The symmetron [72] is another scalar–tensor theory with a Higgs-like potential

$$V\left(\varphi \right)=-\frac{{\mu }^2}{2}{\varphi }^2+\frac{\lambda }{4}{\varphi }^4$$

(57)

and a non-linear coupling function

$$A\left(\varphi \right)=1+\frac{{\varphi }^2}{2M^2}+\cdots .$$

(58)

where the quadratic term is meant to be small compared to unity. The coupling is given by

$$\beta \left(\varphi \right)=\frac{m_{\mathrm{Pl}}\varphi }{M^2}$$

(59)

which vanishes at the minimum of the effective potential when $\rho >{\mu }^2{M}^2$. This realises the Damour–Polyakov mechanism.

2.7.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 non-linear 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

$$S=M_5^3\int d^{5}x\sqrt{-g}R+\int d^{4}x\sqrt{-g}\left[-2M_5^{3}K+\frac{M_4^2}{2}R-\sigma +{\mathcal{L}}_{\mathrm{matter}}\right],$$

(60)

where ${\mathcal{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 the bulk, respectively. The metric $g_{\mu \nu }$ is on the brane, R its Ricci scalar, and $K=g^{\mu \nu }{K}_{\mu \nu }$ is the trace of extrinsic curvature, $K_{\mu \nu }$. Finally, $\sigma$ 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

$$r_5\approx \frac{M_4^2}{M_5^3}.$$

(61)

For scales larger than $r_5$, the five-dimensional physics contributes to the dynamics, while for scales smaller than $r_5$, gravity is four-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

$$w_{\mathrm{eff}}=\frac{P}{\rho }=-\frac{\frac{dH}{dt}+3H^2+\frac{2\kappa }{a^2}}{3H^2+\frac{3\kappa }{a^2}},$$

(62)

where $\kappa$ is the three-dimensional spatial curvature. During the self-accelerating phase $w_{\mathrm{eff}}\to -1$ in (62), therefore emulating a cosmological constant.

2.8. 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

$$S=\int {\mathrm{d}}^{4}x\sqrt{-g}\left[\frac{m_{Pl}^2}{2}R-\frac{1}{2}{(\nabla \varphi )}^2-V\left(\varphi \right)\right]+\int {\mathrm{d}}^{4}x\sqrt{-\tilde{g}}{\mathcal{L}}_m({\psi }_m^{\left(i\right)},{\tilde{g}}_{\mu \nu })$$

(63)

where $\varphi$ is the scalar field, $V\left(\varphi \right)$, its potential and it couples to the matter fields ${\psi }_m^{\left(i\right)}$ through the Jordan frame metric ${\tilde{g}}_{\mu \nu }$, which is related to the metric $g_{\mu \nu }$ as

$${\tilde{g}}_{\mu \nu }=A^2(\varphi ,X)g_{\mu \nu }+B^2(\varphi ,X){\partial }_{\mu }\varphi {\partial }_{\mu }\varphi .$$

(64)

The disformal factor term in $B^2(\varphi ,X)$ leads to the derivative interactions in (2). In the previous discussions, see Section 2.7.4, we focused on the conformal parameter $A\left(\varphi \right)$ chosen to be X-independent where $X=-{(\partial \varphi )}^2/2{\Lambda }^2$ and $\Lambda$ 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(\varphi ,X)$ is crucial [45].

As can be expected, (63) can be generalised 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

$$S=\int {\mathrm{d}}^{4}x\sqrt{-g}\left[\sum _{i=2}^5{\mathcal{L}}_i+{\mathcal{L}}_m({\psi }_m^{\left(i\right)},g_{\mu \nu })\right]$$

(65)

where the four Lagrangian terms correspond to different combinations of four functions $G_{2,3,4,5}$ of the scalar field and its kinetic energy $\chi =-{\partial }^{\mu }{\partial }_{\mu }\varphi /2$, the Ricci scalar and the Einstein tensor $G_{\mu \nu }$ and are given by

$$\begin{array}{ccc}\hfill {\mathcal{L}}_2& =& K(\varphi ,\chi ),\hfill \\ \hfill {\mathcal{L}}_3& =& -G_3(\varphi ,\chi )\square \varphi ,\hfill \\ \hfill {\mathcal{L}}_4& =& G_4(\varphi ,\chi )R+G_{4,\chi }[{(\square \varphi )}^2-\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)\left({\nabla }^{\mu }{\nabla }^{\nu }\varphi \right)],\hfill \\ \hfill {\mathcal{L}}_5& =& G_5(\varphi ,\chi )G_{\mu \nu }\left({\nabla }^{\mu }{\nabla }^{\nu }\varphi \right)\hfill \\ \hfill & -& \frac{1}{6}{G}_{5,\chi }{[(\square \varphi )}^3-3(\square \varphi )\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi \right)\left({\nabla }^{\mu }{\nabla }^{\nu }\varphi \right)\hfill \\ & +& 2\left({\nabla }^{\mu }{\nabla }_{\alpha }\varphi \right)\left({\nabla }^{\alpha }{\nabla }_{\beta }\varphi \right)\left({\nabla }^{\beta }{\nabla }_{\mu }\varphi \right)],\hfill \end{array}$$

(66)

After the gravitational wave event GW170817 ([76,77], and as already anticipated in [78], 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 ${\mathcal{L}}_2$ and Cubic Galileons (Horndeski with Lagrangians up to ${\mathcal{L}}_3$) as the surviving class of models [79,80,81]. 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 [82], the Vainshtein mechanism can break within massive sources [83]. This phenomenology was studied further in [84], and may be used to constrain such theories, as described in Section 4.1.

2.9. 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

$${\beta }_{\mathrm{eff}}^2\le 4·{10}^{-5}$$

(67)

as long as the range of the scalar force exceeds several astronomical units and ${\beta }_{\mathrm{eff}}^2$ corresponds to the strength of the fifth force acting on the satellite. As we have mentioned, this translates into the effective bound

$$\frac{\beta \left({\varphi }_0\right){\beta }_{\odot }}{Z\left({\varphi }_0\right)}\le {10}^{-5}$$

(68)

where ${\varphi }_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 $\beta \left({\varphi }_0\right)$ for the satellite and ${\beta }_{\odot }$ 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\left({\varphi }_0\right)$ factor. As an example, for cubic Galileon models, the ratio between the fifth force and the Newtonian force behaves such as

$$\frac{F_N}{F_{\varphi }}\simeq 2{\beta }^2{\left(\frac{r}{R_V}\right)}^{3/2}$$

(69)

where $\beta \left({\varphi }_0\right)=\beta$ and $R_V$ is the Vainshtein radius. For cosmological models where $L\sim H_0^{-1}$, the Vainshtein radius of the Sun is around 0.1 kpc. As a result, for the planets of the solar system, $r/r_V\ll 1$ and the fifth force is negligible. K-mouflage models of cosmological interest with $\Lambda \simeq {\Lambda }_{\mathrm{DE}}\sim {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

$${\beta }_{\odot }\lesssim {10}^{-5}$$

(70)

where we have assumed that $Z\left({\varphi }_0\right)=1$ and $\beta \left({\varphi }_0\right)\simeq 1$. This is a stringent bound, which translates into

$${\varphi }_0\lesssim {10}^{-11}{m}_{\mathrm{Pl}}$$

(71)

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 ${\Phi }_{N\odot }\simeq {10}^{-6}$.

In fact, chameleon-screened theories are constrained even more strongly by the Lunar Ranging experiment (LLR) [85,86]. This experiment constrains the Eötvos parameter

$${\eta }_{AB}=\frac{|{\overrightarrow{a}}_A-{\overrightarrow{a}}_B|}{|{\overrightarrow{a}}_A+{\overrightarrow{a}}_B|}$$

(72)

for two bodies falling in the presence of a third one C. The accelerations ${\overrightarrow{a}}_{A,B}$ are towards C and due to C. For bodies such as the Earth $A=\oplus$, the moon $B=\mathrm{moon}$ and the Sun $C=\odot$, 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 ${\beta }_A$, which is dependent on the Newtonian potential of the body ${\beta }_A\propto {\Phi }_A^{-1}$ implying a strong dependence on the nature of the objects. As the strength of the gravitational interaction between two screened bodies is given by

$$G_{AB}=G_N(1+2{\beta }_A{\beta }_B)$$

(73)

as long as the two objects are closer than the background Compton wavelength $m^{-1}\left({\varphi }_0\right)$, the Eötvos parameter becomes

$${\eta }_{AB}\simeq {\beta }_C|{\beta }_A-{\beta }_B|$$

(74)

In the case of the LLR experiment, we have ${\beta }_A\simeq {\varphi }_0/2\beta \left({\varphi }_0\right){\Phi }_A{m}_{\mathrm{Pl}}$ and therefore ${\beta }_{\mathrm{moon}}\gg {\beta }_{\oplus }$. Using ${\Phi }_{N\mathrm{moon}}\simeq {10}^{-11}$ and ${\Phi }_{\oplus }\simeq {10}^{-9}$ we find that the LLR constraint

$${\eta }_{\mathrm{LLR}}\le {10}^{-13}.$$

(75)

This becomes for the scalar charge of the Earth [49]

$${\beta }_{\oplus }\lesssim {10}^{-6}$$

(76)

which is stronger than the Cassini bound, i.e., we must impose that

$${\varphi }_0\lesssim {10}^{-15}{m}_{\mathrm{Pl}}.$$

(77)

This corresponds to the energy scale of particle accelerators such as the Large Hadron Collider (LHC). This bound leads to a relevant constraint on the parameter space of popular models. Let us first consider the $n=1$ inverse power-law chameleon model with

$$V\left(\varphi \right)={\Lambda }_{\mathrm{DE}}^4{e}^{\frac{{\Lambda }_{\mathrm{DE}}}{\varphi }}\approx {\Lambda }_{\mathrm{DE}}^4+\frac{{\Lambda }_{\mathrm{DE}}^5}{\varphi }+\cdots .$$

(78)

This model combines the screening property of inverse power-law chameleon and the cosmological constant term ${\Lambda }_{\mathrm{DE}}^4$ leading to the acceleration of the expansion of the Universe. The mass of the scalar is given by

$$m^2\left(\varphi \right)\approx \frac{2{\Lambda }_{\mathrm{DE}}^5}{{\varphi }^3}$$

(79)

implying that in the solar system $m_0\gtrsim {10}^6{H}_0$. Now, as long as the chameleon sits at the minimum of its effective potential, we have $m_{\mathrm{cosmo}}\simeq {\left(\frac{{\rho }_{\mathrm{cosmo}}}{{\rho }_G}\right)}^{1/2}{m}_0$ where ${\rho }_{\mathrm{cosmo}}$ is the cosmological matter density and ${\rho }_G\simeq {10}^6{\rho }_{\mathrm{cosmo}}$ is the one in the Milky Way. As a result, we have the constraints on the cosmological mass of the chameleon [87,88]

$$m_{\mathrm{cosmo}}\gtrsim {10}^3{H}_0.$$

(80)

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 a large scale. In the context of the Hu–Sawicki model, and as $m_{\mathrm{cosmo}}/H_0\propto f_{R_0}^{-1/2}$, the solar system tests imply typically that $f_{R_0}\lesssim {10}^{-6}$[69]. For models with the Damour–Polyakov mechanism such as the symmetron, and if ${\rho }_G\le {\mu }^2{M}^2$, the field value in the solar system is close to ${\varphi }_0\simeq \frac{\mu }{\sqrt{\lambda }}$. The mass of the scalar is also of order $\mu$, implying that the range of the symmetron is very short unless $\mu \lesssim {10}^{-18}$ eV. In this case, the LLR bound applies and leads to

$$M\lesssim {10}^{-3}{m}_{\mathrm{Pl}}$$

(81)

which implies that the symmetron models must be an 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 [89]. Indeed, the interaction potential induced by a screening object does not vary as $1/r$ anymore. As a result, Bertrand’s theorem8 is violated and the planetary trajectories are not closed anymore. For K-mouflage models defined by a Lagrangian $\mathcal{L}={\Lambda }^{4}K\left(X\right)$ where $X=-{(\partial \varphi )}^2/2{\Lambda }^4$ and $\Lambda \simeq {\Lambda }_{\mathrm{DE}}$, the periastron is given by [90]

$$\delta \theta =2\pi {\beta }^2\frac{K^{\prime 2}}{K^{″}}x\frac{dX}{dx}$$

(82)

where $x=r/R_K$ is the reduced radius ($R_K$ is the K-mouflage radius). For the Moon, the LLR experiment implies that $\delta \theta \le 2·{10}^{-11}$ which constrains the function $K\left(X\right)$ and its derivatives $K^{\prime }\left(X\right)$ and $K^{″}\left(X\right)$. A typical example of models passing the solar system tests is given by $K\left(X\right)=-1+X+K_{\star }(X-X_{\star }arctan\left(\frac{X}{X_{\star }}\right))$ with $X_{\star }\simeq 1$ and $K_{\star }\gtrsim {10}^3$. In these models, the screening effect is obtained as $K^{\prime }\sim K_{\star }\gg 1$ as long as $\left|X\right|\gtrsim |X_{\star }|$. For cubic Galileons, the constraint from the periastron of the Moon reduces to a bound on the suppression scale [89,91]

$${\Lambda }^3\gtrsim m_{\mathrm{Pl}}{H}_0^2.$$

(83)

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\left(\varphi \right)$, i.e.

$$G_{\mathrm{eff}}=A^2\left(\varphi \right)G_N.$$

(84)

For chameleon-screened objects, the difference between the Jordan and Einstein values of the Newton constant is irrelevant as deep inside screened objects $\varphi$ are constant and $A\left(\varphi \right)$ 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

$$\varphi (\overrightarrow{x},t)\simeq \varphi \left(\overrightarrow{x}\right)+{\dot{\varphi }}_{\mathrm{cosmo}}(t-t_0)+{\varphi }_{\mathrm{cosmo}}\left(t_0\right)$$

(85)

where $t_0$ is the present time. Here $\varphi \left(\overrightarrow{x}\right)$ 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 [92]

$$\frac{dln{G}_{\mathrm{eff}}}{dt}=\frac{\beta }{m_{\mathrm{Pl}}}\dot{\varphi }$$

(86)

where we have taken the scalar to be coupled conformally with a constant strength $\beta$. The LLR experiment has given a bound in the solar system [85]

$$|\frac{dln{G}_{\mathrm{eff}}}{dt}|\le 0.02H_0,$$

(87)

i.e., Newton’s constant must vary on timescales larger than the age of the Universe. This can be satisfied by K-mouflage or Vainshtein models provided $\beta \le 0.1$ as long as $\dot{\varphi }\sim m_{\mathrm{Pl}}{H}_0$, i.e., the scalar field varies by an order of magnitude around the Planck scale in one Hubble time [90].

3. 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 that 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 be considered to be massless. Quantum effects [93] of the Casimir type imply that two metallic plates separated by a distance d will then interact and attract according to

$$F\left(d\right)=-\frac{{\pi }^2}{480d^4}$$

(88)

as long as the coupling between the scalar and matter is large enough9. In the Casimir or Eötwash context, this would mean that the usual quantum effects due to electromagnetism would 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 setup. 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 gradient10. In the following, we will mostly focus on the classical case and the resulting constraints.

3.1. Casimir Interaction and Eötwash Experiment

We now turn to the Casimir effect [94], 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 [95]

$$\left|\frac{F_{\varphi }}{A}\right|=V_{\mathrm{eff}}\left(\varphi \left(0\right)\right)-V_{\mathrm{eff}}\left({\varphi }_0\right),$$

(89)

where A is the surface area of the plates and $V_{\mathrm{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 ${\varphi }_0$ and in the vacuum chamber halfway between the plates. In general, the field acquires a bubble-like profile between the plates and $\varphi \left(0\right)$ 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 $\beta$, the resulting pressure between two plates separated by distance d is given by

$$\left|\frac{F_{\varphi }}{A}\right|=\frac{{\beta }^2{\rho }_{\mathrm{plate}}^2}{2m_{\mathrm{Pl}}^2{m}^2}{e}^{-md},$$

(90)

which makes the Yukawa suppression of the interaction between the two plates explicit. In the screened case, the situation can be very different.

Let us first focus on the symmetron case. As long as $\mu \gtrsim d^{-1}$, the value $\varphi \left(0\right)$ is very close to the vacuum value $\varphi \left(0\right)\simeq \frac{\mu }{\sqrt{\lambda }}$ implying that $F_{\varphi }/A\simeq 0$, i.e., the Casimir effect does not efficiently probe symmetrons with large masses compared to the inverse distance between the plates. On the other hand, when $\mu \lesssim d^{-1}$, the field essentially vanishes in the plates and between the plates [96]. As a result, the classical pressure due to the scalar becomes

$$\left|\frac{F_{\varphi }}{A}\right|=\frac{{\mu }^4}{4\lambda }.$$

(91)

Notice that in this regime, the symmetron decouples from matter inside the experiment as $\beta \left(\varphi \right(0\left)\right)=0$. We will see how this compares to the quantum effects in Section 3.6. 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 [97]

$$\left|\frac{F_{\varphi }}{A}\right|={\Lambda }^4{\left(\frac{\sqrt{2}B(\frac{1}{2},\frac{2+n}{2n})}{n\Lambda d}\right)}^{2n/(2+n)}.$$

(92)

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 the photon Casimir pressure at large distances. Quantum effects can also be taken into account when the chameleon’s mass is small enough, see Section 3.6.

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 $|\frac{\Delta F_{\varphi }}{A}|\le 0.35$ mPa [98]. The plate density is of the order of ${\rho }_{\mathrm{plate}}=10\mathrm{g}{\mathrm{cm}}^{-3}$. The constraints deduced from the Casimir experiment can be seen in Section 3.5. It should be noted that realistic experiments sometimes employ a plate-and-sphere configuration, which can have an $O\left(1\right)$ modification to  (92) [99].

The Eöt-Wash experiment [100] 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,101]. Defining by $A\left(\theta \right)$ the surface area of the two plates that 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 $\theta$ and is given by

$$T\sim a_{\theta }{\int }_d^{\infty }dx\frac{F_{\varphi }}{A},$$

(93)

where $a_{\theta }=\frac{dA}{d\theta }$ 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

$$T=a_{\theta }\frac{{\beta }^2{\rho }_{\mathrm{plate}}^2}{2m_{\mathrm{Pl}}^2{m}^3}{e}^{-md},$$

(94)

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\lesssim {\mu }^{-1}$, implying that

$$\begin{array}{ccc}& & d\lesssim {\mu }^{-1},T\simeq a_{\theta }\frac{{\mu }^{4}d}{4\lambda }\hfill \\ & & d\gtrsim {\mu }^{-1},T\simeq a_{\theta }\frac{{\mu }^3}{4\lambda }.\hfill \end{array}$$

(95)

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

$$T=a_{\theta }\left(\frac{2+n}{n-2}\right){\Lambda }^{4}d{\left(\frac{\sqrt{2}B(\frac{1}{2},\frac{2+n}{2n})}{n\Lambda d}\right)}^{2n/(2+n)}$$

(96)

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_{\star }\sim m_0^{-1}$, where $m_0$ is the mass in the vacuum between the plates. This becomes

$$T\sim a_{\theta }\left(\frac{2+n}{n-2}\right){\Lambda }^4\left[d_{\star }{\left(\frac{\sqrt{2}B(\frac{1}{2},\frac{2+n}{2n})}{n\Lambda d_{\star }}\right)}^{2n/(2+n)}-d{\left(\frac{\sqrt{2}B(\frac{1}{2},\frac{2+n}{2n})}{n\Lambda d}\right)}^{2n/(2+n)}\right]$$

(97)

for $d\lesssim d_{\star }$ and essentially vanishes for larger distances. In the case $n=2$, a logarithmic behaviour appears.

The 2006 Eöt-Wash experiment [102] provided the bound for a separation between the plates of $d=55 \mathsf{\mu }\mathrm{m}$, which is

$$|T|\le a_{\theta }{\Lambda }_T^3,$$

(98)

where ${\Lambda }_T=0.35{\Lambda }_{\mathrm{DE}}$[35] and ${\Lambda }_{\mathrm{DE}}=2.4$ meV. We must also modify the torque calculated previously in order to take into account the effects of a thin electrostatic shielding sheet of width $d_s=10 \mathsf{\mu }\mathrm{m}$ between the plates in the Eöt-Wash experiment. This reduces the observed torque, which becomes $T_{obs}=e^{-m_c{d}_s}T$. As a result, we have that

$$e^{-m_c{d}_s}{\int }_{55\mu \mathrm{m}}^{\infty }\left|\frac{F_{\varphi }}{A}\right|dx\le {\Lambda }_T^3.$$

(99)

Surprisingly, the Eötwash experiment tests the dark energy scale in the laboratory as ${\Lambda }_T\approx {\Lambda }_{\mathrm{DE}}$.

3.2. 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 [103]. 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

$$\left(-\frac{{\hslash }^2}{2m_N}\frac{d^2}{d{z}^2}+V\left(z\right)\right){\psi }_n=E_n{\psi }_n$$

(100)

where $m_N$ is the neutron’s mass and the potential over a horizontal plate is

$$V\left(z\right)=m_{N}gz+m_N(A\left(\varphi \left(z\right)\right)-1)$$

(101)

where $\varphi \left(z\right)$ is the vertical profile of the scalar field. We put the mirror at $z\le 0$. The contribution due to the scalar field is

$$\delta V\left(z\right)=m_N(A\left(\varphi \left(z\right)\right)-1)$$

(102)

which depends on the model. In the absence of any scalar field, the wavefunctions are Airy functions

$${\psi }_k\left(z\right)=c_k\mathrm{Ai}\left(\frac{z}{z_0}-{ϵ}_k\right)$$

(103)

where $c_k$ is a normalisation constant, $z_0={({\overline{h}}^2/2m_N^{2}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\cdots \}$ for the first levels $k=1\dots$. At the first-order of perturbation theory, the energy levels are

$$E_k=E_k^{\left(0\right)}+\delta E_k$$

(104)

where $E_k^{\left(0\right)}={ϵ}_k{m}_{N}g{z}_0$, and the perturbed energy level

$$\delta E_k=m_N\langle {\psi }_k|(A\left(\varphi \left(z\right)\right)-1)|{\psi }_k\rangle$$

(105)

is the averaged value of the perturbed potential in the excited states.

Let us see what this entails for chameleon models [104,105]. In this case, the perturbation depends on

$$A\left(\varphi \right)-1\simeq \frac{\beta }{m_{\mathrm{Pl}}}\varphi +\cdots$$

(106)

where the profile of the chameleon over the plate is given by

$$\varphi \left(z\right)=\Lambda {\left(\frac{2+n}{\sqrt{2}}\Lambda z\right)}^{2/(n+2)}.$$

(107)

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 $\beta$, one can deduce useful constraints on the parameters of the model. Thus far, we have assumed that the neutrons are not screened. When they are screened, the correction to the energy levels is easily obtained by replacing $\beta \to \lambda \beta$ where $\lambda$ is the corresponding screening factor.

In the case of symmetrons, the correction to the potential energy depends on

$$A\left(\varphi \right)-1=\frac{{\varphi }^2}{2M^2}$$

(108)

whilst the symmetron profile is given by [106]

$$\varphi \left(z\right)=\frac{\mu }{\sqrt{\lambda }}tanh\frac{\mu z}{\sqrt{2}}$$

(109)

where we assume that the plate is completely screened. The averaged values of $\delta V$ are constrained by

$$|\delta E_3-\delta E_1|\lesssim 2·{10}^{-15}\mathrm{eV}$$

(110)

which leads to strong constraints on symmetron models. See Section 3.5.

3.3. Atomic Interferometry

Atomic interferometry experiments are capable of very precisely measuring the acceleration of an atom in free fall [107,108]. 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_{\mathrm{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:

$${\overrightarrow{F}}_{\mathrm{source}}={\overrightarrow{F}}_{\mathrm{N}}+{\overrightarrow{F}}_{\varphi }.$$

(111)

As such, atom interferometry is a sensitive probe of new theories that predict a classical fifth force ${\overrightarrow{F}}_{\varphi }$. In experiments such as [109,110] the source is a ball of matter and the extra acceleration $a_{\varphi }=F_{\varphi }/m_{\mathrm{atom}}$ is determined at the level $a_{\varphi }\lesssim 5.5 \mathsf{\mu }\mathrm{m}/{\mathrm{s}}^2$ at a distance $d=R_B+d_B$ where $R_B=0.95 \mathrm{cm}$ is the radius of the ball and $d_B=0.88 \mathrm{cm}$ is the distance to the interferometer. The whole setup is embedded inside a cavity of radius $R_c=6.1 \mathrm{cm}$.

Scalar fields, of the type considered in this review, generically predict such a force. The fifth force is of the form

$${\overrightarrow{F}}_{\varphi }=-m_{\mathrm{atom}}\overrightarrow{\nabla }lnA\left(\varphi \right),$$

(112)

where $A\left(\varphi \right)$ 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 interacts with the scalar field. In screened theories, it is often written as the product $q_A={\beta }_A{m}_A$ of the mass $m_A$ of the objects and the reduced scalar charge ${\beta }_A$. The reduced scalar charge can be factorised as ${\beta }_A={\lambda }_A\beta \left({\varphi }_0\right)$ where $\beta \left({\varphi }_0\right)$ is a the coupling of a point-particle to the scalar field in the background environment characterised by the scalar field value ${\varphi }_0$. The screening factor ${\lambda }_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\left(\varphi \right)$, the size, mass, and geometry of the object, as well as the ambient scalar field value ${\varphi }_0$. For a spherical object, the screening factor of object A is given by

$${\lambda }_A=\frac{|{\varphi }_A-{\varphi }_0|}{2m_{\mathrm{Pl}}{\Phi }_A}$$

(113)

when the object is screened, otherwise ${\lambda }_A=\beta \left({\varphi }_0\right)$. Here ${\varphi }_A$ is the value of the scalar field deep inside the body A, ${\Phi }_A$ the Newtonian potential at its surface and ${\varphi }_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

$$|F_{\varphi }|={\left(\frac{\beta \left({\varphi }_0\right)}{m_{\mathrm{Pl}}}\right)}^2\frac{\left({\lambda }_A{m}_a\right)\left({\lambda }_B{m}_B\right)}{4\pi r^2}{e}^{-m_{c}r},$$

(114)

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 [111,112], as small objects such as 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\left({\varphi }_0\right)$ of the scalar particles is small, and distance scales of order $\sim \mathrm{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 Equation (114), one only needs to determine the ambient field value ${\varphi }_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 $\rho$, the field will roll to minimise its effective potential $V_{\mathrm{eff}}\left(\varphi \right)$ given by (42) for a value $\varphi \left(\rho \right)$. In a dense region such as the vacuum chamber walls, $\varphi \left(\rho \right)$ is small, while in the rarefied region inside the vacuum chamber $\varphi \left(\rho \right)$ is large. The field thus starts at a small value ${\varphi }_{\min ,\mathrm{wall}}$ near the walls and rolls towards a large value ${\varphi }_{\min ,\mathrm{vac}}$ near the centre. However, the field will only reach ${\varphi }_{\min ,\mathrm{vac}}$ if the vacuum chamber is sufficiently large. The energy of the scalar field depends upon both potential energy $V\left(\varphi \right)$ and gradient energy ${\left(\overrightarrow{\nabla }\varphi \right)}^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]

$$m_0\left({\varphi }_{\mathrm{vac}}\right)\equiv \frac{1}{R}.$$

(115)

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

$$m_0\left({\varphi }_{\min ,\mathrm{vac}}\right)\gg \frac{1}{R}.$$

(116)

As such, the field inside the vacuum chamber is

$${\varphi }_0=min({\varphi }_{\min ,\mathrm{vac}},{\varphi }_{\mathrm{vac}}).$$

(117)

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\%$ [110]. More exotic shapes have been shown to boost the sensitivity even further, by up to a factor $\sim 3$ [113].

Atom interferometry experiments of this type, with an in-vacuum source mass, have now been performed by two separate groups [109,114,115]. In these experiments, the acceleration between an atom and a marble-sized source mass has been constrained to $a\lesssim 50\mathrm{nm}/{\mathrm{s}}^2$ at a distance of $r\lesssim 1\mathrm{cm}$. These experiments have placed strong bounds on the parameters of chameleon and symmetron [116] modified gravity, as will be detailed in Section 3.5.

3.4. Atomic Spectroscopy

In the previous section, we saw that the scalar field mediates a new force, Equation (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 Equation (114) perturbs the electron’s Hamiltonian [117,118]

$$\delta H={\left(\frac{\beta \left({\varphi }_0\right)}{m_{\mathrm{Pl}}}\right)}^2\frac{{\lambda }_p{m}_p{m}_e}{r},$$

(118)

where ${\lambda }_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

$$\delta E_n=\langle {\psi }_n|\widehat{\delta H}|{\psi }_n\rangle ,$$

(119)

where ${\psi }_n$ are the unperturbed electron’s eigenstates.

This was first computed for a generic scalar field coupled to matter with a strength $\beta \left(\varphi \right)=m_{\mathrm{Pl}}/M$ [119], using measurements of the hydrogen 1s-2s transition [120,121,122] to rule out

$$M\lesssim \mathrm{TeV}.$$

(120)

However, that study did not account for the screening behaviour exhibited by chameleon and symmetron theories. That analysis was recently extended to include screened theories and the work will be published shortly, resulting in the bound that is illustrated in Figure 1.

3.5. Combined Laboratory Constraints

Combined bounds on theory parameters derived from the experimental techniques detailed in this section are plotted in Figure 1 and Figure 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 [91] 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_{\mathrm{Pl}}/\beta$. 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\gtrsim m_{\mathrm{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\ne 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 Figure 2. Unlike the chameleon, the symmetron has a mass parameter $\mu$ that fixes it to a specific length scale ${\mu }^{-1}$. For an experiment at a length scale L, if $L\gg {\mu }^{-1}$ then the fifth force would be exponentially suppressed, as is clear in Equation (114). Likewise, in an enclosed experiment, if $L\ll {\mu }^{-1}$ then the energy considerations in the previous subsection imply that the field simply remains in the symmetric phase where $\varphi =0$. The coupling to matter is quadratic,

$${\mathcal{L}}_{\mathrm{symm}}\supset -\frac{\beta \left(\varphi \right)}{m_{\mathrm{Pl}}}\rho \equiv -\frac{{\varphi }^2}{M^2}\rho ,$$

(121)

so in the symmetric phase, where $\varphi =0$, the coupling to matter switches off and the fifth force vanishes. Therefore, to test a symmetron with mass parameter $\mu$, one must test it with an experiment on a length scale $L\approx {\mu }^{-1}$.

3.6. Quantum Constraints

Classical physics effects induced by the 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 due to their vacuum fluctuations. This typically occurs in Casimir experiments where two plates attract each other or the Eötwash setting where two plates face each other.

Particle physics experiments test the nature of the interactions of new states to the standard model at very high energy. In particular, the interactions of the light scalars to matter and the gauge bosons of the standard are via the Higgs portal, i.e., the Higgs field couples both to the standard model particles and the light scalar and as such mediates the interactions of the light scalar to the standard model. This mechanism is tightly constrained by the precision tests of the standard model. For instance, the light scalars will have an effect on the fine structure constant, the mass of the Z boson or the Fermi interaction constant $G_F$. The resulting bound on $\beta =\frac{m_{\mathrm{Pl}}}{M}$ is [124]

$$M\gtrsim {10}^3\mathrm{GeV}$$

(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

$$m^2=V^{″}\left(\varphi \right)+A^{″}\left(\varphi \right)\rho .$$

(123)

The density is piece-wise constant and labelled ${\rho }_{1,2,3}$ in the case of a Casimir experiment. Here ${\rho }_2$ is the density between the plates. Notice that as $\varphi$ is continuous, the mass jumps across boundaries as $\rho$ 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 [125]

$${\overrightarrow{F}}_{\varphi }=-\int d^{3}x{\overrightarrow{\partial }}_d\rho \langle A\left(\varphi \right)\rangle$$

(124)

where the integration is taken over all space and ${\overrightarrow{\partial }}_d\rho$ 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 ${\rho }_3$ positioned along the x-axis between $x=d$ and $x=d+L$, the vacuum of density ${\rho }_2$ between $x=0$ and L, a plate of density ${\rho }_1$ for $x\le 0$ and finally again the vacuum for $x>d+L$ we have ${\partial }_d\rho =({\rho }_3-{\rho }_2)(\delta (x-d-L)-\delta (x-d))$ and the force is along the x-axis. The quantum average $\langle A\left(\varphi \right)\rangle$ is taken over all the quantum fluctuations of $\varphi$. When the field has a classical profile ${\varphi }_{\mathrm{clas}}$, this quantum calculation can be performed in perturbation theory

$$\langle A\left(\varphi \right)\rangle =A\left({\varphi }_{\mathrm{clas}}\right)+\frac{A^{″}\left({\varphi }_{\mathrm{clas}}\right)}{2}\Delta (x,x)+\cdots$$

(125)

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^{\prime }$ 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\left(\varphi \right)$ in a power series are associated with higher loop contributions to the force when the first term is given by a one-loop diagram. The Feynman propagator $\Delta (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

$$\frac{F_x}{A}=-\frac{A^{″}}{2}({\rho }_3-{\rho }_2)(\Delta (d,d)-\Delta (d+L,d+L))$$

(126)

where we have considered that the derivative $A^{″}\left({\varphi }_{\mathrm{clas}}\right)\simeq A^{″}$ is nearly constant. This is exact for symmetron models and chameleon models with ${\varphi }_{}\ll M$. As the classical solution is continuous at the boundary between the plates, the quantum force is in fact given by

$$\frac{F_x}{A}=\frac{m_2^2-m_3^2}{2}(\Delta (d,d)-\Delta (d+L,d+L))$$

(127)

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 $\Delta (x,x)$ are x-independent, we see immediately that they cancel in the force (127), which is finite. Moreover, the limit $L\to \infty$ is finite and corresponds to the case of an infinitely wide plate. Notice that the contribution in $-\Delta (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-vacuum-plate) 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

$$\frac{F_x}{A}=-\frac{1}{2{\pi }^2}{\int }_0^{\infty }d\rho {\rho }^2\frac{{\gamma }_2({\gamma }_2-{\gamma }_1)({\gamma }_2-{\gamma }_3)}{e^{2d{\gamma }_2}({\gamma }_1+{\gamma }_2)({\gamma }_2+{\gamma }_3)-({\gamma }_2-{\gamma }_1)({\gamma }_2--{\gamma }_3)}$$

(128)

with ${\gamma }_i^2={\rho }^2+m_i^2$. When the density in the plates becomes extremely large compared to the one in the vacuum, the limit $m_{1,3}\to \infty$ gives the finite result

$$\frac{F_x}{A}=-\frac{1}{2{\pi }^2}{\int }_0^{\infty }d\rho {\rho }^2\frac{{\gamma }_2}{e^{2d{\gamma }_2}-1}.$$

(129)

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 Figure 3 for symmetrons. In a nutshell, when the $\mu$ 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 nonzero, 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^{″}\simeq 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\left(\varphi \right)$ 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\gtrsim 1/\sqrt{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}\to \infty$, 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

$$a_f=\frac{g_f-2}{2},$$

(130)

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 [119,126,127]

$$\delta a_f\approx 2{\left(\frac{\beta \left(\varphi \right)m_f}{4\pi m_{\mathrm{Pl}}}\right)}^2,$$

(131)

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$ [127,128].

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}$ [129]. Setting $\delta a_e\le {10}^{-12}$ yields the constraint [127]

$$\beta \left(\varphi \right)\lesssim {10}^{16}.$$

(132)

In the case of the chameleon where $\beta =m_{\mathrm{Pl}}/M$, this rules out $M<80\mathrm{GeV}$.

In the case of the muon, the experimental measurement of the magnetic moment [130,131] and the standard model prediction [132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151] differ by 1 part in ${10}^9$ at $4.2\sigma$. A generic scalar field without a screening mechanism cannot account for this discrepancy without also being in tension with Solar System tests of gravity. However, it has recently been shown that both the chameleon and symmetron are able to resolve this anomaly while also satisfying all other experimental bounds [128]. The chameleon parameters that accomplish this are

$$M=500\mathrm{GeV}\mathrm{and}\Lambda <\mathrm{meV}.$$

(133)

Cosmologically, a chameleon with these parameters has an effective mass $m_{\mathrm{cosmo}}>{10}^{-13}\mathrm{eV}$ and Compton wavelength $<{10}^3\mathrm{km}$, so this theory does not significantly influence our universe on large scales.

4. 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\left(R\right)$ (taken as a paradigmatic chameleon-screened theory; [69,152,153]) 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 achieved 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 Section 2). 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 [154,155,156] (see Figure 4). The threshold value of potential for screening is given by Equation 52: in $n=1$ Hu–Sawicki $f\left(R\right)$, $\chi \simeq \frac{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 non-linearities 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 space-time curvature (Section 2.5).

Several methods have been developed to build “screening maps” of the local universe to identify screened and unscreened objects. Shao et al. [157] apply an $f\left(R\right)$ scalar field solver to a constrained N-body simulation to estimate directly the scalar field strength as a function of position. Cabre et al. [154] 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. [158] adopt a similar approach but include more contributions to the potential, 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 the local environment, with existing tests using one of the final two.

4.1. 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 by

$$\frac{\mathrm{d}P}{\mathrm{d}r}=-\frac{G_{N}M\left(r\right)\rho \left(r\right)}{r^2}.$$

(134)

In the presence of a thin-shell-screened fifth force, this becomes

$$\frac{\mathrm{d}P}{\mathrm{d}r}=-\frac{G_{N}M\left(r\right)\rho \left(r\right)}{r^2}\left[1+2{\beta }^2\left(1-\frac{M\left(r_{\mathrm{s}}\right)}{M\left(r\right)}\right)\Theta (r-r_{\mathrm{s}})\right],$$

(135)

with $\Theta \left(x\right)$ the Heaviside step function, $\beta$ the coupling coefficient of the scalar field and $r_{\mathrm{s}}$ the screening radius of the star beyond which it is unscreened. In the case of chameleon theories, the factor $1-\frac{M\left(r_{\mathrm{s}}\right)}{M\left(r\right)}$ corresponds to the screening factor and is associated with 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\propto G_N^4$ while high-mass stars have $L\propto G_N$[159]. Thus in the case that the star is fully unscreened ($r_{\mathrm{s}}=0$), low-mass stars have L boosted by a factor ${(1+2{\beta }^2)}^4$, and high-mass stars by $(1+2{\beta }^2)$.

To explore the full effect of a fifth force on the behaviour of stars, Equation (134) must be coupled with the equations describing stellar structure and energy generation. This has been achieved by modifying the stellar structure code MESA [159,160,161,162], enabling the heuristic expectations described above to be quantified (see Figure 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

$$\ddot{\overrightarrow{\delta r}}=-\frac{1}{\rho }\frac{\mathrm{d}P}{\mathrm{d}r}+\overrightarrow{a},$$

(136)

with $\overrightarrow{a}$ the force per unit mass, which is $\overrightarrow{a}=-\overrightarrow{\nabla }\Phi$ in GR but $\overrightarrow{a}=-\overrightarrow{\nabla }\Phi -\frac{\beta }{m_{\mathrm{Pl}}}\overrightarrow{\nabla }\varphi$ in the presence of a scalar field. Combining this equation with the other stellar structure equations gives the frequency of linear radial adiabatic oscillations

$${\omega }^2\sim \frac{G_{N}M}{R^3},$$

(137)

so that enhancing the effective value of G due to the addition of a fifth force causes the pulsation period $\Pi$ to change according to

$$\frac{\Delta \Pi }{\Pi }=-\beta Q,$$

(138)

where Q is the star’s scalar charge.

Stellar oscillations are useful observationally because they provide several methods of determining distances to galaxies [163]. 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 $\kappa$-mechanism [164] 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 $\chi \lesssim {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 $\chi$, 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 the galaxy environment—has yielded the constraint $f_{R0}\lesssim {10}^{-7}$[165,166]. Notice that astrophysical constraints yield tighter bounds on $f\left(R\right)$ models than solar system tests.

Variable stars are also useful for more general tests of gravity. [167] 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 six 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 [166,168,169].

Finally, other types of stars 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 lowest mass of any red dwarf observed [170,171]. It would also affect the radii of brown dwarf stars and the mass–radius relation and Chandresekhar mass of white dwarfs [172].

4.2. 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. [173] fitted the rotation curves of 85 late-type galaxies with an $f\left(R\right)$ model, finding evidence for $f_{R0}\approx {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” [174] (see also [175]). Screening can also generate new correlations between dynamical variables; for example, Burrage et al. [176] use a symmetron model to reproduce the Radial Acceleration Relation linking the observed and baryonic accelerations in galaxies [177]. 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 $\Lambda$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 $\sim {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 [178,179]:

$$\frac{v_g^2}{r}=\frac{G_N(1+2{\beta }^2)M(<r)}{r^2},\frac{v_{*}^2}{r}=\frac{G_{N}M(<r)}{r^2}\Rightarrow \frac{v_g}{v_{*}}=\sqrt{1+2{\beta }^2},$$

(139)

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 $\beta$ and [OIII] emission lines in six low-surface brightness galaxies, Vikram et al. [180] place the constraint $f_{R0}\lesssim {10}^{-6}$. This result can likely be significantly strengthened by increasing the sample size using data from IFU surveys such as MaNGA or CALIFA—potentially 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 [181].

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 ${\overrightarrow{a}}_{\varphi }$ 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 ${\overrightarrow{a}}_{\varphi }$. 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 ${\overrightarrow{a}}_{\varphi }$ so that all parts of the galaxy have the same total acceleration [178,179]:

$$\frac{G_{N}M(<r_{*})}{r_{*}^2}{\hat{r}}_{*}=2\beta {\overrightarrow{a}}_{\varphi },$$

(140)

where ${\overrightarrow{r}}_{*}$ is the displacement of the stellar and gas centroids. This effect is illustrated in Figure 6a 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 ${\overrightarrow{a}}_{\varphi }$. This is shown in Figure 6b. 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:

$$z=\frac{2\beta a_{\varphi }{r}^3}{G_{N}M(<r)},$$

(141)

which can be simplified by assuming a halo density profile. Desmond et al. [182,183,184] 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\times {10}^{-8}$ at $1\sigma$ 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 [185] (see Figure 7a). Subsequent work has verified using hydrodynamical simulations that the baryonic noise model used in these analyses is accurate [186]. The value of ${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\sim \mathcal{O}\left(100\right)$ kpc [37], the larger values $r_c\sim 6000$ Mpc required for self-acceleration [187] must be probed on galactic or cosmological scales. The most promising method for this utilises the breaking of the Strong Equivalence Principle (SEP) that Galileons imply [188] 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 $\beta$, 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 cluster environments: An extended distribution of mass does not Vainshtein-screen perfectly in its interior [189], so a residual fifth-force field is present in cluster outskirts. This leads to $\mathcal{O}$(kpc) offsets between black holes and satellite galaxy centres for realistic cluster parameters. Sakstein et al. [190] 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 $\mathcal{O}$(1) coupling coefficients for $r_c\lesssim 1$ Gpc. Second, the Galileon symmetry implies that the linear contribution to the field on cosmological scales is unscreened [191,192], allowing black hole offsets to develop even for field galaxies. Assuming a constant density ${\rho }_0$ in the centre of the halo, the black hole offset in this case is given by [188]

$$R=0.1\mathrm{kpc}\left(2{\beta }^2\right)\left(\frac{|\nabla {\Phi }_N^{\mathrm{ext}}|}{20{(\mathrm{km}/\mathrm{s})}^2/\mathrm{kpc}}\right)\left(\frac{0.01M_{\odot }{/\mathrm{pc}}^3}{{\rho }_0}\right),$$

(142)

where $\nabla {\Phi }_N^{\mathrm{ext}}$ is the unscreened large-scale gravitational field, proportional to the Galileon fifth-force field. Bartlett et al. [193] 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 $\beta <0.28$ for $r_c\gtrsim 1/H_0$ (see Figure 7b). This probes the cosmologically-relevant region of the Galileon parameter space, complementing cosmological probes such as the Integrated Sachs Wolfe (ISW) effect (see Section 5). 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 [190].

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 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 [194,195] and simulation [196,197] 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 [79,198]. Voids can also be cross-correlated with galaxies to infer the growth rate of structure [199], used in the ISW effect [200], integrated along the line of sight to produce projected 2D voids [201], and used as a means of splitting samples of galaxies into high-density (screened) and low-density (unscreened) environments or in marked correlation functions [202,203]. Finally, the redshift-space anisotropy of voids is a powerful probe of the nature of gravity through redshift space distortions [204]. Future surveys will improve 3D spectroscopic void finding and the calibration of photometric void finders with robust photometric redshifts.

4.3. Galaxy Cluster Tests

A fifth force causes a structure to grow more quickly, leading to more cluster-sized halos at late times. This is, however, counteracted by screening and the Yukawa cut off due to the mass of the scalar field so that cluster abundance only deviates significantly from the $\Lambda$CDM expectation at lower masses and in sparser environments [205]. The excursion set formalism for halo abundance provides a good description under chameleon gravity as well [206], albeit with a modified collapse threshold ${\delta }_c$, and has been used to constrain $f_{R0}\lesssim {10}^{-5}$ in the Hu–Sawicki model [207,208]. Similar constraints are achievable using the peaks in the weak lensing convergence field, which trace massive halos [209]. Other formalisms for calculating cluster abundance in the presence of a fifth force have also been developed [210,211,212,213]. 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 [214]. 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 $\Lambda$CDM. This is illustrated in Figure 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 the 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 [207,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, non-relativistic 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 $\gamma$. Although $\gamma$ may be constrained on $\mathcal{O}$(kpc) scales by comparing strong lensing and stellar dynamical mass estimates [224,225], it has found the most use on the scale of clusters. An approximation for chameleon theories of the Jordan–Brans–Dicke type is [226]

$$M_{\mathrm{dyn}}\left(r\right)\simeq \left(1+\frac{\Theta (r-r_{\mathrm{c}})}{3+2{\omega }_{\mathrm{BD}}}\left[1-\frac{M\left(r_{\mathrm{c}}\right)}{M\left(r\right)}\right]\right)M{\left(r\right)}_{\mathrm{lens}},$$

(143)

where $\Theta$ is the Heaviside step function, ${\omega }_{\mathrm{BD}}$ is the JBD parameter (see Section 2.7.5) and the radius at which the scalar field transitions to its background value is given by

$$r_{\mathrm{c}}\simeq \frac{32\pi G\rho \left(r_{\mathrm{s}}\right)r_{\mathrm{s}}^3}{3+2{\omega }_{\mathrm{BD}}}\frac{1}{1-A^{-2}\left({\varphi }_{\mathrm{env}}\right)}-r_{\mathrm{s}}.$$

(144)

Here ${\varphi }_{\mathrm{env}}$ is the cosmological boundary condition for the field far from the cluster (e.g., $1-A^{-2}\left(\varphi \right)\simeq f_{R0}$ in the $f\left(R\right)$ 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\left(R\right)$ gravity has also been calibrated from N-body simulations in [227]:

$$\frac{M_{\mathrm{dyn}}}{M_{\mathrm{lens}}}=\frac{7}{6}-\frac{1}{6}tanh\left[p_1\left({log}_{10}(M_{\mathrm{lens}}/M_{\odot })-p_2\right)\right],$$

(145)

where $p_1=2.21$ and $p_2=1.503{log}_{10}\left[\frac{|f_R\left(z\right)-1|}{1+z}\right]+21.64$. This works well for $f_{R0}\in \left[{10}^{-6.5},{10}^{-4}\right]$ and $z\in [0,1]$. To test this effect, strong cluster lensing may be compared to X-ray masses or the velocity dispersions of the cluster galaxies [189,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}\lesssim 6\times {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\left(R\right)$ to the $f_{R0}\sim 5\times {10}^{-5}$ level [233]. All of these tests will benefit from enlarged cluster samples in the future.

5. Cosmological Consequences

5.1. 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 Section 2.3 and Section 2.4, 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\left(R\right)$. 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 such as Cubic Galileons, which feature the Vainshtein screening as a way to evade small scale gravitational constraints.

No-go theorems [35,87,88] were developed for chameleon-screened theories, and they state namely that (i) the Compton wavelength of such scalars can be at most $\simeq 1\mathrm{Mpc}$ at the present cosmic density, which means that the effective range of these theories is restricted to non-linear 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]

$${\omega }_{\varphi }+1\simeq \mathcal{O}\left(\frac{H^2}{m^2}\right)$$

(146)

where m is the mass of the light scalar. The bound from solar systems on the mass ratio $m/H\gtrsim {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.

5.2. Screening and Structure Formation

The formation of a 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 modifying gravity on small scales.

As an example and on cosmological scales, the $f\left(R\right)$ modification of the Einstein–Hilbert action leads to a modified Poisson equation, which can be expressed as

$${\nabla }^2\Phi =\frac{16\pi G}{3}{a}^2\delta \rho -\frac{a^2}{6}\delta R,$$

(147)

in comoving coordinates and the term $\delta \rho$ is the matter density fluctuation compared to the cosmological background and $\Phi$ the modified Newtonian potential. Furthermore, the fluctuation of the Ricci scalar, $\delta R=R-\overline{R}$ compared to the cosmological background $\overline{R}$ and is expressed as

$${\nabla }^2\delta f_R=\frac{a^2}{3}[\delta R-8\pi G\delta \rho ].$$

(148)

The variation of the function $f\left(R\right)$ is given by $\delta f_R=f_R\left(R\right)-f_R\left(\overline{R}\right)$. In these equations, we have assumed a quasi-static approximation. It can be shown [234] that despite the fact that these equations are non-linear in $\delta R$, they are self-averaging. This means that on large scales one recovers $\delta R\to 0$. Using these governing equations, one can solve perturbatively the Vlasov–Poisson system of equations, which consists in the 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 [234,235,236,237,238,239].

In scalar-tensor theories with screening and a conformal factor $A\left(\varphi \right)$, particles feel a total gravitational potential $\Phi$, which is the sum of the standard Newtonian term ${\Phi }_N$ and an additional contribution ${\Phi }_A$,

$$\Phi ={\Phi }_N+{\Phi }_A,$$

(149)

where the governing equations are given by

$$\begin{array}{cc}\hfill \frac{1}{a^2}\Delta {\Phi }_N=4\pi G\delta \rho ,& {\Phi }_A=ln\frac{A}{\overline{A}}\simeq (A-\overline{A})\end{array}$$

(150)

where it is assumed that $A\left(\varphi \right)\simeq 1$ to satisfy constraints on the variation of fermionic masses. As a result, $lnA\simeq A-1$ and the dependence on $lnA$ of the Newtonian potential $\Phi$ becomes linear in A. This additional gravitational potential implies that matter particles of mass m are sensitive to a “fifth force” given by

$${\overrightarrow{F}}_{\varphi }=-m\overrightarrow{\nabla }lnA.$$

(151)

This fifth force is the one that leads to a modification of the growth of structures.

5.3. 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,240]. However, recent constraints on the equation-of-state parameter of dark energy are overall consistent with a cosmological constant $w\approx -1$ [241]. 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 [242], cannot be satisfactorily resolved with late-time dynamics of a dark energy field, according to the latest analysis [243,244] and therefore will not be covered in this section. Therefore, in the following section, 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).

5.3.1. 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 of information about the constituents of the Universe; while its secondary anisotropies, which 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) ([245,246,247] that affect the power spectrum at low multipoles (large scales) and lensing of the CMB [248,249] 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

$$\frac{\delta T}{T}\left(\hat{n}\right)=-{\int }_{{\eta }_0}^{{\eta }_{*}}\mathrm{d}\eta \frac{\partial (\Psi +\Phi )}{\partial \eta },$$

(152)

where ${\eta }_{*}$ is the conformal time at the last scattering surface and ${\eta }_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 parts left from the Horndeski theory after GW170817. In [214], cubic Galileons were analysed, and it was found that in the presence of massive neutrinos (model dubbed $\nu$ Galileon, in red in Figure 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 [250] and WMAP-9 polarisation [251] data. For BAO they used 6dF, SDSS DR7 and BOSS DR9 data ([252,253,254]). In the absence of massive neutrinos (model dubbed Galileon in Figure 9), however, $\Lambda$CDM  was favoured by the data. Nevertheless, they showed that the $\nu$ Galileon model shows a negative ISW effect that is hard to reconcile with current observations. More recently, a paper by [255] performed a detailed study of self-accelerating Galileon models using CMB data from Planck-15 in temperature and polarisation and CMB lensing [256]. 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$\sigma$ tension with observations, effectively ruling this model out. Furthermore, in [257], the effect of different neutrino masses and hierarchies was analysed, and it was also found out that all cubic, quartic and quintic Galileons remain ruled out by CMB and ISW observations.

5.3.2. Cosmological Perturbations in Large Scale Structure

As mentioned above in the corresponding sections for $f\left(R\right)$ 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 $\Phi$ and $\Psi$. In the quasistatic approximation, these equations read

$$\begin{array}{ccc}\hfill -k^2\Phi (a,k)& =& 4\pi G{a}^2\mu (a,k)\rho \left(a\right)\Delta (a,k);\hfill \end{array}$$

(153)

$$\begin{array}{ccc}\hfill -k^2{\Phi }_K(a,k)& =& 4\pi G{a}^2\Sigma (a,k)\rho \left(a\right)\Delta (a,k).\hfill \end{array}$$

(154)

where $\rho \left(a\right)$ is the average dark matter density and $\Delta (a,k)=\delta +3aH\theta$ is the comoving density contrast with $\delta$, the fractional overdensity, and $\theta$ the peculiar velocity. We have denoted by

$${\Phi }_L(a,k)=\frac{\Phi (a,k)+\Psi (a,k)}{2}$$

(155)

the lensing potential. The ratio of the two gravitational potentials is denoted as $\eta$, gravitational anisotropic stress or gravitational slip

$$\eta (a,k)\equiv \Psi (a,k)/\Phi (a,k).$$

(156)

The scale and time-dependent functions $\eta (a,k)$, $\mu (a,k)$ and $\Sigma (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 $\mu$, $\Sigma$ and $\eta$ and they are related by

$$\Sigma (a,k)=\frac{\mu (a,k)}{2}\left[1+\eta (a,k)\right].$$

(157)

The $\mu (a,k)$ function is usually probed best by galaxy clustering experiments that directly trace the evolution of the $\Phi$ potential, since this one affects non-relativistic particles. $\mu$ is directly related to the effective Newtonian constant defined above in (9) as

$$G_{eff}/G_N=\mu$$

(158)

in the linear regime and in Fourier space. On the other hand, relativistic particles, and therefore light, follows the equation for $\Phi (a,k)+\Psi (a,k)$, meaning that gravitational weak lensing is mostly sensitive to the function $\Sigma (a,k)$.

$f\left(R\right)$ Models and Chameleon Theories

For the $f\left(R\right)$ 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 with the $f\left(R\right)$ models [53]

$$m_{f_R}^2=\frac{f_R}{3f_{RR}}\sim \frac{1}{3f_{RR}}$$

(159)

where we have used that $R\simeq \rho /m_{\mathrm{Pl}}^2$ at a late time in the Universe. Neglecting the anisotropic stress, the expressions for $\mu$ and $\eta$ read [258]

$$\begin{array}{cc}\hfill \mu (a,k)& =\frac{1}{f_R\left(a\right)}\frac{4+3k^2{a}^{-2}{m}_{f_R}^{-2}\left(a\right)}{3(1+k^2{a}^{-2}{m}_{f_R}^{-2}\left(a\right))},\hfill \end{array}$$

(160)

$$\begin{array}{cc}\hfill \Sigma \left(a\right)& =\frac{1}{f_R\left(a\right)},\hfill \end{array}$$

(161)

Given the constraints on $f_{R,0}$ mentioned above, the modifications of lensing are practically non-existent and $\Sigma (a,k)\simeq 1$ with great precision.

It is convenient to rewrite the above expressions as

$$\begin{array}{ccc}\hfill \mu (a,k)& & =A^2\left(a\right)\left(1+\frac{2\beta {\left(a\right)}^2}{1+\frac{m^2\left(a\right)a^2}{k^2}}\right),\hfill \end{array}$$

(162)

$$\begin{array}{ccc}\hfill \Sigma (a,k)& & =A^2\left(a\right),\hfill \end{array}$$

(163)

where in the case of $f\left(R\right)$ models, we have $\beta \left(a\right)=\beta =1/\sqrt{6}$ and $m_{f_R}\left(a\right)=m\left(a\right)$, where $\beta \left(a\right)$ is the coupling at the minimum of the effective potential $V_{\mathrm{eff}}\left(\varphi \right)=V\left(\varphi \right)+(A\left(\varphi \right)-1)\rho$ as a function of a with $\rho \propto a^{-3}$ and similarly for the mass of the scalar field $m\left(a\right)$. 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\left(\varphi \right)$ and coupling function $A\left(\varphi \right)$, which is called the tomographic map. This allows to parameterise the chameleon models with the function $m\left(a\right)$ and $\beta \left(a\right)$. The mapping reads [88]

$$\frac{\varphi \left(a\right)}{m_{\mathrm{Pl}}}=\frac{{\varphi }_{\mathrm{ini}}}{m_{\mathrm{Pl}}}+9{\int }_{a_{\mathrm{ini}}}^{a}dx\frac{\beta \left(x\right){\Omega }_m\left(x\right)H^2\left(x\right)}{x{m}^2\left(x\right)}$$

(164)

where ${\Omega }_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 $\Lambda$CDM  at the background level. The potential itself is given by

$$V\left(a\right)=V_{\mathrm{ini}}-\frac{3}{m_{\mathrm{Pl}}^2}{\int }_{a_{\mathrm{ini}}}^{a}dx\frac{{\beta }^2\left(x\right){\rho }^2\left(x\right)}{x^2{m}^2\left(x\right)}.$$

(165)

This provides a parametric reconstruction of $V\left(\varphi \right)$. For the Hu–Sawicki models of $f\left(R\right)$, we have [53]

$$m_{{\mathrm{f}}_{\mathrm{R}}}\left(a\right)=m_0{\left(\frac{\frac{{\Omega }_{m0}}{a^3}+4{\Omega }_{\Lambda 0}}{{\Omega }_{m0}+4{\Omega }_{\Lambda 0}}\right)}^{(n+2)/2}$$

(166)

where ${\Omega }_{\Lambda }$ is the dark energy fraction and ${\Omega }_{m0}$ the matter fraction now. The mass of the scalaron now is given by

$$m_0=\frac{H_0}{\sqrt{(n+1)|f_{R0}|}}{\left(4{\Omega }_{\Lambda 0}+{\Omega }_{m0}\right)}^{1/2}.$$

(167)

which is greater than $H_0$ for small $|f_{R0}|\ll 1$.

Finally, the $\mu$ parameterisation allows one to see how screening works on cosmological linear scales [50]. Defining the comoving Compton wavelength

$${\lambda }_c\left(a\right)=\frac{1}{am\left(a\right)}$$

(168)

we find that for scales outside the Compton wavelength, i.e., $k\lesssim {\lambda }_c$ we have

$$\mu (a,k)\simeq A^2\left(a\right)\simeq 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\gtrsim {\lambda }_c$, we have an enhancement of the gravitational interaction as

$$\mu (a,k)\simeq 1+2{\beta }^2\left(a\right)$$

(170)

which is simply due to the exchange of the nearly-massless scalar field between overdensities.

As a result, we can have a qualitative description of chameleon models such as $f\left(R\right)$ on the growth of structures [259]. First of all, on very large scales, GR is retrieved and no deviation from $\Lambda$-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.

Jordan–Brans–Dicke Models

For the JBD models with a mass term, these functions are given by [260,261]

$$\begin{array}{cc}\hfill \mu (a,k)& =\frac{1}{\overline{\varphi }}\frac{2(2+{\omega }_{BD})}{3+2{\omega }_{BD}},\hfill \end{array}$$

(171)

$$\begin{array}{cc}\hfill \eta (a,k)& =\frac{2+{\omega }_{BD}}{1+{\omega }_{BD}},\hfill \end{array}$$

(172)

so that for cosmological purposes

$$\Sigma \left(a\right)=1.$$

(173)

In this case, lensing is not affected at all.

Horndeski Theory

For a generic Horndeski theory (of second-order in the equations of motion), these two functions $\mu$ and $\eta$ 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 [258,260]

$$\begin{array}{cc}\hfill \mu (a,k)& =\frac{p_1\left(a\right)+p_2\left(a\right)k^2}{1+p_3\left(a\right)k^2},\hfill \end{array}$$

(174)

$$\begin{array}{cc}\hfill \eta (a,k)& =\frac{1+p_3\left(a\right)k^2}{p_4\left(a\right)+p_5\left(a\right)k^2}.\hfill \end{array}$$

(175)

There is another physically more meaningful parametrisation of the linear Horndeski action, given by [262], which is related to the effective field theory of dark energy [32,263,264], where small deviations to the background cosmology are parameterised linearly. This parametrisation is of great help when discussing current cosmological constraints. It is defined using four functions of time ${\alpha }_M$, ${\alpha }_K$, ${\alpha }_B$ and ${\alpha }_T$ plus the effective Planck mass $M_{\star }^2$ and a function of time for a given background specified by the time variation of the Hubble rate $H\left(a\right)$ as a function of the scale factor a. The term ${\alpha }_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 ${\alpha }_K$ quantifies the kineticity of the scalar field and therefore appears in models such as K-mouflage, which require the K-mouflage screening. The coefficient ${\alpha }_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 [265]. It receives contributions also from terms related to the cubic Galileons, which present the Vainshtein screening. Finally, ${\alpha }_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\left(R\right)$ models, where the chameleon screening is necessary, as we have seen.

DGP Models

Cosmological linear perturbations for DGP have been worked out in [266]. In the paper by [260], it is assumed that the small-scale (quasi-static) approximation is valid, i.e., $k/a\gg r_5\mathcal{H}$ and obtains

$$-k^2\Psi =4\pi G_N\left(1-\frac{1}{3\gamma }\right)\overline{\rho }{a}^2\delta ,$$

(176)

and

$$-k^2\Phi =4\pi G_N\left(1+\frac{1}{3\gamma }\right)\overline{\rho }{a}^2\delta ,$$

(177)

where $\gamma =1+2ϵH{r}_5{w}_{\mathrm{eff}}$. This corresponds to

$$\mu \left(a\right)=1+\frac{1}{3\gamma },$$

(178)

and for all practical purposes, we can set $\Sigma =1$ within the cosmological horizon (see [260]).

5.4. Large Scale Structure Observations: Galaxy Clustering and Weak Lensing

The most important probes for a large scale structure, especially in the upcoming decade with the advent of new observations by DESI [267]12, Euclid [268,269]13, Vera Rubin [270]14 and WFIRST [271]15, will be galaxy clustering and weak lensing. Galaxy clustering measures the 2-point-correlation function of galaxy positions either in three dimensions, i.e., angular positions and redshift, or in effectively two 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_{gg}^{obs}$ is directly related to the power spectrum of matter density perturbations $P_{\delta \delta ,zs}$ in redshift space by

$$P_{\mathrm{gg}}^{\mathrm{obs}}(z,k,{\mu }_{\theta })=\mathrm{AP}\left(z\right)P_{\delta \delta ,\mathrm{zs}}(k,z)E_{\mathrm{err}}(z,k)+P_{\mathrm{shot}}\left(z\right),$$

(179)

where $AP\left(z\right)$ corresponds to the Alcock–Paczynski effect, $E_{err}(z,k)$ is a damping term given by redshift errors and $P_{\mathrm{shot}}\left(z\right)$ is the shot noise from estimating a continuous distribution out of a discrete set of points. ${\mu }_{\theta }$ is the cosine of the angle between the line of sight and the wave vector $\mathbf{k}$. Furthermore, the redshift space power spectrum is given by

$$P_{\delta \delta ,\mathrm{zs}}(z,k,{\mu }_{\theta })=\mathrm{FoG}(z,k,{\mu }_{\theta })K^2(z,{\mu }_{\theta };b\left(z\right);f\left(z\right))P_{\delta \delta }(k,{\mu }_{\theta },z),$$

(180)

where $\mathrm{FoG}(z,k,{\mu }_{\theta })$ 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 [272]—on the growth rate $f\left(z\right)$ and the bias $b\left(z\right)$, but can be more complicated when taking into account non-linear perturbation theory at mildly non-linear scales. For a detailed explanation of these terms, we refer the reader to [273] 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 [274,275]. This may be observable in upcoming spectroscopic surveys such as DESI [276]. Reference [277] 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 degenerated 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 become distorted, their ellipticities increased and their magnitudes changed, due to light travelling through large scale structures in the Universe, from the source to the observer [278]. 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

$$C_{ij}^{\gamma \gamma }(\ell )=\frac{c}{H_0}\int \frac{{\hat{W}}_i^{\gamma \left(z\right)}{\hat{W}}_j^{\gamma }\left(z\right)}{E\left(z\right)r^2\left(z\right)}{P}_{\Phi +\Psi }\left(k_{\ell },z\right)dz,$$

(181)

where $E\left(z\right)=H\left(z\right)/H_0$ is the dimensionless Hubble function, ${\hat{W}}_j^{\gamma \left(z\right)}$ are window functions, or lensing kernels, that project the redshift distributions and the power spectrum into angular scales, and finally, $P_{\Phi +\Psi }\left(k_{\ell },z\right)$ is the Weyl power spectrum, which is related to the matter power spectrum $P_{\delta \delta }$ by

$$P_{\Phi +\Psi }={\Sigma }^2(k,z){\left[3{\left(\frac{H_0}{c}\right)}^2{\Omega }_{\mathrm{m}}^0(1+z)\right]}^2{P}_{\delta \delta }.$$

(182)

In this equation, we can see clearly the observational signature of the $\Sigma$ lensing function defined above in (154) and (157). We refer the reader again to [273] 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 $\Lambda$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 $\Lambda$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 ${\Omega }_{m,0}=0.315$, ${\Omega }_{b,0}=0.05$, $h=0.674$, $ns=0.966$ and ${\sigma }_8=0.8156$. For JBD, the model parameter ${\omega }_{BD}$ is set to ${\omega }_{BD}=800$, while for nDGP, the observational parameter is ${\Omega }_{rc}=c^2/\left(4r_5^2{H}_0^2\right)$, where $r_5$ is the crossover scale defined above in (61) and we set here ${\Omega }_{rc}=0.25$. For K-mouflage, the physical parameter is ${ϵ}_2=dlnA\left(a\right)/dlna$, and it is related to the fifth force enacted by the scalar field, which comes from the conformal transformation of the metric (see [279] for more details). The prediction shown here is made for the case ${ϵ}_2=-0.04$.

These distinctive features can be observed when taking the ratio to $\Lambda$CDM  for the three cosmological models considered above. While at linear scales $k\lesssim 0.07\mathrm{h}/\mathrm{Mpc}$, the models show only a slight change in amplitude compared to $\Lambda$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. Furthermore, the ratio of the weak lensing $C_{\ell }$ with respect to $\Lambda$CDM  is shown in the lower panel. In this case, the very sharp features observed in the matter power spectrum are smoothed out by the projection along the line of sight and into angular multipoles.

5.5. 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 CAMB16 (Code for Anisotropies in the Microwave Background) (see [280]), which is written mainly in fortran, and CLASS17 (see [281,282]), 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 is codes in which generic parametrisations of the deviations of GR as in (153) to (156) are used. The second one is codes in which specific modified gravity (MG) models or a generic class of models are implemented and their full scalar field equations are solved, beyond the quasi-static approximation. From the first type, the two more common are ISitGR18 (see [283,284], and MGCAMB19 [285,286] and more recently a branch of CLASS, called QSA_CLASS (see [287]). For the second type, we will mention here the two most important ones, namely hi_class20 (see [288] and EFTCAMB21 (see [289,290]).

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\approx 0.1{\mathrm{Mpc}}^{-1}$, where Baryon Acoustic Oscillations (BAO) and redshift space distortions (RSD) are important, needs to be treated perturbatively in order to make accurate predictions. This involves using either Eulerian or Lagrangian perturbation theory [291,292] and furthermore using resummation techniques to capture accurately the large scale contributions [293,294]. 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 [295]. For the models we are interested in here, there has been some recent work by [235,296], some new work on GridSPT by [297] and some more foundational work by [298,299].

To obtain meaningful constraints with weak lensing, the power spectrum needs to be calculated at even higher k-values, for up to $k\approx 10{\mathrm{Mpc}}^{-1}$, which is only possible using N-body simulations, which capture the full evolution of the non-linear 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 non-linear evolution of the field needs to be taken into account. Several interesting approaches have been taken in the literature such as COLA [237], Ramses [300], Ecosmog [301], MG-evolution [302], $\varphi$-enics (an interesting finite element method approach that can capture the non-linear evolution of the scalar field and reproduce very accurately the Vainshtein screening) [303] and the simulation work on $f\left(R\right)$ theories by several groups [304,305,306,307]. 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 in forecasts and Monte Carlo parameter estimation. Several approaches include fitting formulae based on simulations [308], emulators for $f\left(R\right)$ theories [309,310] and hybrid approaches in which the halo model, perturbation theory and simulations are calibrated to create a model, such as REACT (see [311,312]). This code can compute predictions for nDGP and $f\left(R\right)$ models, which are roughly 5% accurate at scales $k\lesssim 5\mathrm{h}/\mathrm{Mpc}$.

5.6. 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 structures.

5.6.1. Constraints on $f\left(R\right)$ Models

From the CMB, constraints have been placed by the Planck collaboration on the $f\left(R\right)$ model in terms of the Compton wavelength parameter, which is defined as

$$\begin{array}{ccc}\hfill B& \equiv & \frac{f_{RR}}{f_R}{R}^{\prime }\frac{H}{H^{\prime }},\hfill \end{array}$$

(183)

and its value today $B_0$ is related to the fundamental parameter $f_{R,0}$. Indeed, we have the relation

$$B=\frac{{\Omega }_m}{1+\omega }\frac{H^2}{m_{f_R}^2}$$

(184)

where $\omega$ is the equation-of-state of the Universe and $m_{f_R}$ is the mass of the scalaron (166). Notice that the denominator $1+\omega$ is very small, and therefore, B is less suppressed than the ratio $H^2/m_{f_R}^2$. In the analysis of [313], 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 $\tau$. 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\times {10}^{-4}$ was found, which forces the model to be very close to $\Lambda$CDM.

A very comprehensive, but by now relatively outdated, analysis by [314] using WMAP5 CMB data [315] and cross-correlations 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\left(R\right)$ models that follow the same expansion of the universe as $\Lambda$CDM  they obtained a bound of $B_0<0.4$ at the 95% confidence level (CL). In the analysis by [316], large scale structure data coming from WiggleZ, BAO (from 6dF, SDSS DR7 and BOSS DR9, see [252,253,254]) were combined with Planck-2013 CMB [250] data and WMAP polarisation data [315] to find ${log}_{10}{B}_0<-4.07$ at the 95% CL. A more recent paper [317] uses the designer approach to $f\left(R\right)$ 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 equations-of-state is reached, which was then constrained to be $|w+1|<0.002$ (95%CL). All these bounds imply that $f\left(R\right)$ models cannot be self-accelerating, and also, if they are present, their background expansion will be very close to the one of $\Lambda$CDM  according to observational bounds. This confirms the known results from gravitational tests in the solar system.

5.6.2. 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 [318,319] 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 $\Lambda$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 [320], it was shown that the growth rate values estimated from the BOSS DR12 data [321] constrains the crossover scale $r_5$ of DGP gravity in the combination ${\left[r_5{H}_0\right]}^{-1}$, which has to be $<0.97$ at the $2\sigma$ level, which amounts to $r_5>3090\mathrm{Mpc}/\mathrm{h}$, meaning that $r_5\sim H_0^{-1}$, 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 [322,323].

5.6.3. 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 ${\omega }_{BD}>40,000$ (see [33,324]). 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 [325], the authors used Planck [250], WMAP [315], SPT and ACT [326,327] 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 found ${\omega }_{BD}>692$ at 99% C.L. When the scalar was free and varied as a parameter, they found ${\omega }_{BD}>890$, which amounts to $0.981<G_{\mathrm{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, polarisation, and lensing reconstruction, the Pantheon supernova distances, BOSS measurements of BAO, along with the joint $3\times 2$pt 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 ${\omega }_{BD}>1540$ at the 95% C.L. and the effective gravitational constant, $G_{\mathrm{eff}}/G=0.997\pm 0.029$. They also found that the uncertainty in the gravitational theory alleviates the tension between KiDS, 2dFLenS and Planck to below $1\sigma$ and the tension [242] in the Hubble constant between Planck and the local measurements to $3\sigma$. Despite these improvements, a careful model selection analysis shows no substantial preference for JBD gravity relative to $\Lambda$CDM.

5.6.4. Constraints on Horndeski Theories and Beyond

For Horndeski models, there has been a great effort by the Planck collaboration to test the parametrised deviations of GR such as in (174) and (175) or in the $\alpha$-formalism of [262]. 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 5.5.

In the Planck 2015 modified gravity paper [313], the authors considered Horndeski models with ${\alpha }_M=-{\alpha }_B$, ${\alpha }_T={\alpha }_H=0$, and ${\alpha }_K$ was fixed to a constant. This amounts to consider non-minimally coupled K-mouflage type models as in [262], with the only free function being ${\alpha }_M$. Additionally, the analysis used the ansatz,

$${\alpha }_M={\alpha }_M^{\mathrm{today}}{a}^p$$

(185)

where ${\alpha }_M^{\mathrm{today}}$ is a constant and $p>0$ determines its backward time evolution. Furthermore, they relate the evolution of ${\alpha }_M$ to a linear (p=1) and exponential ($p>1$, varying free) parametrisation [313]. Using the Planck TT+TE+EE+BSH data set combination (BSH standing again for BAO, SN and local Hubble constraints) they find ${\alpha }_M^{\mathrm{today}}<0.043$ (95% confidence level) for the linear case and ${\alpha }_M^{\mathrm{today}}<0.062$ and $p=0.{92}_{0.24}^{0.53}$ (95% confidence level) for the exponential case. $\Lambda$CDM is recovered for $p=1$ and ${\alpha }_M^{\mathrm{today}}=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 then limits the interest on an overall analysis of Horndeski models in general. However, in [328], 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 [265] that $(\Sigma -1)(\mu -1)\ge 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. References [329,330] studied DHOST models that present self-acceleration, and Reference [331], among others, have studied the astrophysical signatures of these models. However, their theoretical modelling has not been implemented yet in computational tools capable of analysing the full Planck CMB dataset. Finally, the authors of [332] performed a cosmological constraint analysis, assuming the form ${\alpha }_i={\alpha }_{i,0}{a}^{\kappa }$ on these surviving Horndeski models, and using Planck and BICEP2/Keck [333] CMB data and galaxy clustering data from SDSS and BOSS, they found that when setting the kineticity to the following value ${\alpha }_K=0.1a^3$, the ${\alpha }_{M,0}$ parameter has an upper limit of 0.38 when ${\alpha }_{B,0}\ne 0$ and 0.41 when ${\alpha }_{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.

6. Conclusions and Perspectives

Scalar-tensor theories are among the most generic and plausible extensions to $\Lambda$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 commonly 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 environments on a range of scales.

The screening mechanisms can be divided into 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 in 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 periastron of planets and effective violation of the strong equivalence principle in galaxies, both of which afford tight constraints. However, there is still some capability for ground-based experiments to test Vainshtein-screened theories [91]. 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 setup 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 $\mu \approx 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 $\varphi =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 the order of the symmetry-breaking value ${\rho }_{\mathrm{amb}}\approx {\mu }^2{M}^2$, then the symmetron is essentially massless. This enables even long-ranged experiments to test symmetron theories with $\mu \gg 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,\Lambda \}$ 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 Figure 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 such as 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\left(R\right)$ 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, a parameterisation that exists and depends only on the mass and coupling dependence as a function of the scale factor the Universe allows to reconstruct the whole dynamics of the models, on all scales [53,88]. The same type of parameterisation exists for K-mouflage where the coupling function $A\left(a\right)$ and the screening factor $Z\left(a\right)$ are enough to reconstruct the whole dynamics too [279]. 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 that 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 because 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 a 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\times {10}^{-8}$ [185], 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 [334], 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 [335,336]. 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 $\Lambda$CDM, dark matter. This is a key target for many upcoming astrophysical and cosmological surveys. Much less is known about screening in 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.