Measuring gravity at cosmological scales

TL;DR

This paper provides a pedagogical overview of how to test modifications of gravity, especially Horndeski theories, using cosmological observations of large-scale structure and the cosmic microwave background.

Contribution

It reviews model-independent methods for constraining Horndeski scalar-tensor theories through cosmological data, emphasizing minimal assumptions.

Findings

Identifies key measurable quantities to test gravity at cosmological scales.

Highlights the importance of late and early universe observations for testing gravity.

Provides a framework for model-independent constraints on modified gravity theories.

Abstract

This paper is a pedagogical introduction to models of gravity and how to constrain them through cosmological observations. We focus on the Horndeski scalar-tensor theory and on the quantities that can be measured with a minimum of assumptions. Alternatives or extensions of General Relativity have been proposed ever since its early years. Because of Lovelock theorem, modifying gravity in four dimensions typically means adding new degrees of freedom. The simplest way is to include a scalar field coupled to the curvature tensor terms. The most general way of doing so without incurring in the Ostrogradski instability is the Horndeski Lagrangian and its extensions. Testing gravity means therefore, in its simplest term, testing the Horndeski Lagrangian. Since local gravity experiments can always be evaded by assuming some screening mechanism or that baryons are decoupled, or even that the effects of modified gravity are visible only at early times, we need to test gravity with cosmological observations in the late universe (large-scale structure) and in the early universe (cosmic microwave background). In this work we review the basic tools to test gravity at cosmological scales, focusing on model-independent measurements.

Tables

Table 1. Table 1: Comparison table for the different definition of energy density and pressure of the scalar field

	This work	Ref. Zumalacarregui et al. (2017)	Ref. Bellini and Sawicki (2014)
Energy density	${\rho }_{HL}$	${\rho }_{DE}=\frac{{\rho }_{HL}}{3}-H^2\left(M_{\ast }^2-1\right)$	$\widetilde{ℰ}=\frac{{\rho }_{HL}}{M_{\ast }^2}\equiv {\tilde{\rho }}_{HL}$
Pressure	$p_{HL}$	$p_{DE}=\frac{p_{HL}}{3}+\left(3H^2+2\dot{H}\right)\left(M_{\ast }^2-1\right)$	$\widetilde{𝒫}=\frac{p_{HL}}{M_{\ast }^2}\equiv {\tilde{p}}_{HL}$

Table 2. Table 2: The reconstructed or measured model-independent variables E , E ′ , P 2 , P 3 , η ​ ( z ) 𝐸 superscript 𝐸 ′ subscript 𝑃 2 subscript 𝑃 3 𝜂 𝑧 E,E^{\prime},P_{2},P_{3},\eta(z) at three different redshifts z = ( 0.294 , 0.58 , 0.86 ) 𝑧 0.294 0.58 0.86 z=(0.294,0.58,0.86) , together with their 1 ​ σ 1 𝜎 1\sigma errors, for each of the reconstruction methods. The polynomial regression method is compabitle with the Λ Λ \Lambda CDM scenario while the other two methods show some tension at lower redshift.

Method	Parameter		Redshift bins		Weighted mean
		$z_1=0.294$	$z_2=0.58$	$z_3=0.86$
	$E\left(z\right)$	$1.12\pm 0.01$	$1.27\pm 0.02$	$1.51\pm 0.02$
	$E^{\prime }\left(z\right)$	$-0.56\pm 0.07$	$-0.60\pm 0.36$	$-1.75\pm 0.66$
Binning	$P_2\left(z\right)$	$0.75\pm 0.10$	$0.54\pm 0.07$	$0.18\pm 0.14$
	$P_3\left(z\right)$	$-0.17\pm 0.35$	$0.53\pm 0.61$	$-1.27\pm 1.52$
	${\eta }_{obs}\left(z\right)$	$0.48\pm 0.45$	$-0.03\pm 0.34$	$-2.78\pm 6.84$	$0.15\pm 0.27$
	$E\left(z\right)$	$1.10\pm 0.01$	$1.30\pm 0.02$	$1.55\pm 0.03$
	$E^{\prime }\left(z\right)$	$-0.73\pm 0.05$	$-1.30\pm 0.10$	$-1.89\pm 0.16$
Gaussian Process	$P_2\left(z\right)$	$0.74\pm 0.09$	$0.53\pm 0.06$	$0.23\pm 0.11$
	$P_3\left(z\right)$	$-0.10\pm 0.20$	$-0.03\pm 0.21$	$-0.21\pm 0.30$
	${\eta }_{obs}\left(z\right)$	$0.38\pm 0.23$	$0.91\pm 0.36$	$0.58\pm 0.93$	$0.53\pm 0.19$
	$E\left(z\right)$	$1.12\pm 0.01$	$1.29\pm 0.02$	$1.50\pm 0.02$
	$E^{\prime }\left(z\right)$	$-0.73\pm 0.04$	$-1.06\pm 0.04$	$-1.45\pm 0.04$
Polynomial Regression	$P_2\left(z\right)$	$0.76\pm 0.15$	$0.55\pm 0.15$	$0.18\pm 0.14$
	$P_3\left(z\right)$	$-0.09\pm 0.80$	$0.14\pm 0.78$	$-0.17\pm 3.02$
	${\eta }_{obs}\left(z\right)$	$0.57\pm 1.05$	$0.48\pm 0.96$	$-0.11\pm 3.21$	$0.49\pm 0.69$

Equations

$$S=\int d^{4}x\sqrt{-g}\sum _{i=2}^5{\mathcal{L}}_i+S_m$$

$$\displaystyle\begin{aligned} \mathcal{L}_{2}=&K(\phi,X)\,,\\ \mathcal{L}_{3}=&-G_{3}(\phi,X)\Box\phi\,,\\ \mathcal{L}_{4}=&G_{4}(\phi,X)R+G_{4,X}\left[\left(\Box\phi\right)^{2}-\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}\right]\,,\\ \mathcal{L}_{5}=&G_{5}(\phi,X)G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi-\frac{G_{5,X}}{6}\Bigl{[}\left(\Box\phi\right)^{3}-3\left(\Box\phi\right)\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{2}+2\left(\nabla_{\mu}\nabla_{\nu}\phi\right)^{3}\Bigr{]}\,.\end{aligned}$$

$$g_{\mu \nu }={\eta }_{\mu \nu }+h_{\mu \nu }$$

$$d{s}^2$$

$$S_g=\frac{1}{2}\int d^{3}xdt[8B_{,i}^{,i}\dot{\Phi }+4{\Phi }_{,i}{\Psi }_{,i}-4\dot{\Phi }{\dot{E}}_{,i}^{,i}+2{\Phi }_{,i}^2-6{\dot{\Phi }}^2]$$

$$d{s}^2=-d{t}^2+({\delta }_{ij}+h_{ij})d{x}^{i}d{x}^j$$

$$\bar{h}_{ij}=h_{ij}-\frac{1}{2}\eta_{ij}h=\left(\begin{array}[]{ccc}h_{+}&h_{\times}&0\\ h_{\times}&-h_{+}&0\\ 0&0&0\end{array}\right)$$

$$S=\int d^{3}xdt\{Q_S[{\dot{\phi }}^2-\frac{c_s^2}{a^2}({\partial }_i\phi {)}^2]+\sum _{\alpha =1}^2{Q}_T[{\dot{h_{\alpha }}}^2-\frac{c_T^2}{a^2}({\partial }_i{h}_{\alpha }{)}^2]\}$$

$$d{s}^2$$

$${\delta }^{\prime }$$

$${\theta }^{\prime }$$

$${\hat{k}}^2\Psi$$

$$\Psi$$

$${\rho }_m=3H^2{\Omega }_m$$

$${\delta }^{\prime }$$

$${\theta }^{\prime }$$

$${\delta }^{\prime \prime }+\left(2+\frac{H^{\prime }}{H}\right){\delta }^{\prime }-\frac{3}{2}{\Omega }_m\delta =0$$

$${\varphi }^{\prime \prime }+F{\varphi }^{\prime }+c_S^2{k}^2\varphi =S$$

$$k_{\mathrm{phys}}=\frac{k_{\mathrm{com}}}{a}$$

$$k^2\Phi$$

$$k^2\Psi$$

$$\eta =-\frac{\Phi }{\Psi }$$

$$\eta$$

$$Y\eta$$

$${\alpha }_t$$

$${\alpha }_s$$

$$m^2$$

$$\eta$$

$$\eta =1-\frac{2{\alpha }_{t}F{k}^2}{(1+{\alpha }_t)F{k}^2+M^2},Y=\frac{1}{F}\left(1+\frac{{\alpha }_t{k}^2}{k^2+M^2}\right).$$

$$Y=h_1\left(1+\frac{{\alpha }_t{k}^2}{m^2+k^2}\right),$$

$${\rho }_m(k)=\int \rho (r)e^{i\mathbf{kr}}{d}^{3}r$$

$$k^2\Psi =-\frac{1}{2}Y(k){\rho }_m(k)$$

$$r^{-2}\frac{\partial }{\partial r}\left(r^2\frac{\partial }{\partial r}\right)\Psi =\frac{1}{2}\sigma (r)$$

$$\sigma (r)=\frac{V}{(2\pi {)}^3}\int e^{i\mathbf{kr}}Y(k){\rho }_m(k)d^{3}k$$

Full text

Alternatives or extensions of General Relativity have been proposed ever since its early years. Because of Lovelock theorem, modifying gravity in four dimensions typically means adding new degrees of freedom. The simplest way is to include a scalar field coupled to the curvature tensor terms. The most general way of doing so without incurring in the Ostrogradski instability is the Horndeski Lagrangian and its extensions. Testing gravity means therefore, in its simplest term, testing the Horndeski Lagrangian. Since local gravity experiments can always be evaded by assuming some screening mechanism or that baryons are decoupled, or even that the effects of modified gravity are visible only at early times, we need to test gravity with cosmological observations in the late universe (large-scale structure) and in the early universe (cosmic microwave background). In this work we review the basic tools to test gravity at cosmological scales, focusing on model-independent measurements.

I Introduction

Gravity is the force that shapes the overall temporal and spatial structure of the Universe. There is no much need then to explain why it is important to test its validity at all scales and regimes. The huge progress in collecting cosmological data achieved in the last couple of decades makes possible, for the first time, to test gravity and measure its properties at astrophysical and cosmological scales. In order to test a theory one either has to build a set of alternatives against which to compare the standard model, or to parametrize the deviations from it in some meaningful and general way: both approaches are referred to as “modified gravity”.

Lovelock’s theorem Lovelock (1972) states that Einstein’s gravity is the unique local diffeomorphism invariant theory of a tensor field in 4D with second-order equations of motion. It is clear then that modifying gravity often implies adding new degrees of freedom, either scalars, vectors, or tensors. Adding a mass to the graviton, for instance, requires an additional tensor field; including more derivatives is also equivalent to adding more propagating degrees of freedom. Other options based on torsion, non-metricity, or non-locality can also be contemplated, see for instance the review Clifton et al. (2012).

In this paper we review the main properties of an important class of modified gravity based on a single scalar field, the so-called Horndeski Lagrangian (HL). This model is general enough to display most of the phenomenology of non-Einsteinian gravity: generalized Poisson equation, Yukawa corrections to Newton’s potential, presence of anisotropic stress, change in the gravitational wave speed, instabilities, ghosts. Still, the HL is relatively simple in that contains a single propagating degree of freedom in addition to General Relativity. Using the HL as a paradigm of modified gravity, we focus on its observability at various scales, from the local environment to galaxy clusters, with emphasis on cosmological observations. Although the recent measurement of the gravitational wave speed Abbott et al. (2017a) severely constrain one of the phenomenological time-dependent parameters of the Horndeski model, as we will see, the other parameters are still mostly unconstrained and open to theoretical and observational investigation. A major topic of this review is the question of which properties of gravity can be measured as model-independently as possible.

We will not try to cover exhaustively the field of research in modified gravity; good reviews are already available Clifton et al. (2012); Papantonopoulos (2015). Rather, we wish to discuss pedagogically some aspects or issues that are generic to the quest for traces of modified gravity.

We assume units such that $c=8\pi G_{N}=M_{\textrm{Planck}}^{-2}=1$ and metric signature $-+++$. An overdot denotes derivation with respect to cosmic time $t$, a prime with respect to $\log a$. A comma will refer to partial derivative, i.e. $\partial_{\mu}\phi\equiv\phi_{,\mu}$. Also $\Box\phi=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\phi$ where $\nabla_{\mu}$ is the covariant derivative. Greek indexes run over space and time coordinates, Latin indexes over space coordinates only.

II Beyond Einstein

The re-discovery of the most general scalar-tensor theory that gives second order equations of motion, Horndeski action Horndeski (1974) or Covariant Galileons Deffayet et al. (2009), and their extensions Zumalacarregui and García-Bellido (2014); Gleyzes et al. (2015a, b); Crisostomi et al. (2016a, b); Ben Achour et al. (2016); Langlois and Noui (2016) provides a very general framework for such theories (see Kobayashi (2019) for a recent review). The HL is defined as the sum of four terms $\mathcal{L}_{2}$ to $\mathcal{L}_{5}$. Defining with $X=-g_{\mu\nu}\phi^{,\mu}\phi^{,\nu}/2$ the canonical kinetic term, the four terms are specified by two non-canonical kinetic functions $K(\phi,X)$ and $G_{3}(\phi,X)$ and by two coupling functions $G_{4,5}(\phi,X)$, all of them in principle arbitrary:

[TABLE]

where $S_{m}$ is the action for matter fields – dark matter, baryons and radiation – and

[TABLE]

Note that $G_{3}$ and $G_{5}$ must have an $X$ dependence, otherwise they are total derivatives and could be rewritten – after integration by parts – as $K$ and $G_{4}$ respectively.111Notice that the number of these functions cannot be reduced by fields redefinitions without going beyond Horndeski action Bettoni and Liberati (2013); Zumalacarregui and García-Bellido (2014). As usual, each term in the HL has dimension $mass^{4}$. Often one chooses the scalar field to have dimensions of mass, but this is not necessary. As already mentioned, the Horndeski Lagrangian is the most general Lagrangian for a single scalar which gives second-order equations of motion for both the scalar and the metric on an arbitrary background. This is a necessary, but not sufficient, condition for the absence of instabilities, as we will see later on. The terms $\mathcal{L}_{4},\mathcal{L}_{5}$ couple the field $\phi$ to the Ricci scalar $R$ and the Einstein tensor $G_{\mu\nu}=R_{\mu\nu}-Rg_{\mu\nu}/2$. As a consequence, $G_{4,5}$ are the gravity-modifying coupling function. The background equations of motion of the HL are given for completeness in the Appendix, although we do not need them in the following. It is enough to realize that the large freedom offered by the HL allows one to find a background evolution that satisfies all observational constraints.

Let us now briefly discuss some useful limits of the HL.

• •

If $G_{4}=1/2$ and $G_{5}=0$ (it is actually sufficient $G_{5}=const$) the HL reduces to the Einstein–Hilbert Lagrangian with a scalar field having a non-canonical kinetic sector given by $\mathcal{L}_{2},\mathcal{L}_{3}$. The canonical form is obtained for $K=X-V(\phi)$ and $G_{3}=0$ ($G_{3}=const$ is sufficient). $\Lambda$CDM is recovered for $K=-2\Lambda$.

• •

The “minimal” form of modified gravity within the HL is provided by $G_{4}=G_{4}(\phi)$ and $G_{5}=const$: this is then equivalent to a Brans–Dicke scalar-tensor model, again with a non-canonical kinetic sector.

• •

The original Brans-Dicke model is recovered assuming a kinetic sector, $K=(\omega_{BD}/\phi)X,G_{3}=0$, and $G_{4}(\phi)=\phi/2$.

• •

If the kinetic sector vanishes, $K_{,X}=G_{3}=0$, then we reduce ourselves to a $f(R)$ model Sotiriou and Faraoni (2010), whose Lagrangian is $\mathcal{L}_{R}=(R+f(R))/2$. In fact, this model is equivalent to a scalar-tensor theory with $G_{4}(\phi)=e^{2\phi/\sqrt{6}}/2$ and a potential $K(\phi)=-(Rf_{,R}-f)/2$ where $\phi=\sqrt{6}/2\log(1+f_{,R})$. This relation should then be inverted to get $R=R(\phi)$ and used to replace $R$ with $\phi$ in $K(\phi)$.

• •

If one sets $G_{i}(\phi,X)=G_{i}(X)$ then the Lagrangian is invariant under the shift $\phi\rightarrow\phi+c$ with $c=const$. This shift-symmetric version of the HL is connected to the Covariant Galileon when the functional dependence of the $G_{i}$ is fixed Deffayet et al. (2009) and is able to produce the accelerated expansion without a potential that makes the field slow roll.

In general, the equations of motion for the scalar will couple it to the matter energy density. The full set of equations of motion has been studied in several papers, for instance in De Felice et al. (2011); De Felice and Tsujikawa (2012). Any modification of the HL, or addition of terms (except the so-called Beyond Horndeski terms), based on the same scalar field, will introduce higher order equations of motion and associated instabilities, as a consequence of the Ostrogradsky theorem Ostrogradski (1850); Woodard (2015).222See Chen et al. (2013) for a discussion on how to exorcise Ostrogradski ghosts in non-degenerate theories. Of course one can in principle add several scalar fields, but on grounds of simplicity this is rather unnatural. Notice that we do not demand that the $\phi$ drives the present-day accelerated expansion. It could be, after all, that the modification of gravity and the accelerated expansion are independent phenomena. It would be very interesting, though, to explain the latter in terms of the former.

III Decomposition in modes and stability

Einstein’s gravity is carried by a massless spin-2 field, the metric. Being represented by a symmetric matrix, a metric in four dimensions has 10 degrees of freedom (DOF). These DOF can be collected according to how they behave under spatial rotations, i.e. as scalars, vectors and tensors. There are then four scalars (4 DOF), two divergence-free vectors (4 DOF) and one traceless, divergence-free tensor (2 DOF). However, only the tensor DOF propagate, that is, they are subject to linearized equations of motion second-order in the time derivatives. The other DOF obey constraint equations, fully determined by the matter content. This should have been expected, since a massless tensor field, as the gravitational field, has only two independent degrees of freedom.

The two propagating degrees of freedom are associated to the two polarizations $+,\times$ of the gravitational waves. To see that there are no other propagating DOF, one can proceed by linearly expanding the metric around Minkowski

[TABLE]

and keeping only scalar terms, i.e. functions that can be obtained from scalar or from derivatives of scalars. The most general such metric is then

[TABLE]

Inserting this metric into the Einstein-Hilbert Lagrangian without matter and developing to second order, one finds the second-order action in Minkowski space

[TABLE]

The linearly perturbed equations of motion can be obtained then by the Euler-Lagrange equations with respect to $\Phi,\Psi,E,B$, but here we need only identify the degrees of freedom. When one varies the action with respect to $B$, one gets the constraint $\dot{\Phi}=0$ which then shows that $\Phi$ is not a propagating DOF. The same is true for $\Psi$, since there are no time derivatives for it. As we know, in fact, the potentials $\Phi,\Psi$ are determined by the matter distribution through two constraints, the Poisson equations, which do not involve time derivatives. So there are no scalar propagating DOF in Einstein gravity without matter.

The same holds for the vector degrees of freedom. If one considers instead the tensor DOF in $h_{\mu\nu}$

[TABLE]

where, after imposing the traceless, divergence-less conditions, and considering a wave propagating in direction $x_{3}$,

[TABLE]

one finds that the two modes $h_{\alpha}=\{h_{+},h_{\times}\}$ obey in vacuum the same gravitational wave equation, $\Box h_{\alpha}=0$, analogous to electromagnetic waves. GWs propagate therefore with speed $c_{T}$ equal to unity.

The same procedure can be applied to the HL. One finds then in absence of matter fields De Felice and Tsujikawa (2012); Bellini and Sawicki (2014)

[TABLE]

where $\varphi$ is the scalar mode perturbation and $h_{\alpha}$ the two tensor modes, and $c_{S},c_{T}$ their speed of propagation, respectively. As expected, HL has now three propagating DOF, plus those belonging to the matter sector.

The four coefficients $Q_{S},c_{S},Q_{T},c_{T}$ depend on the HL functions. Their expression will be given in Sec. VI. From the classical point of view, stability is guaranteed when $Q_{x},c_{x}^{2}$ (with $x=S,T$) have the same sign. In this case in fact the equations of motion are well-behaved wave equations with speed $c_{x}$, whose amplitude is constant (or decaying in an expanding space), rather than growing exponentially as it would happen for $c_{x}^{2}<0$ (gradient instability). For the quantum stability, however, one must also require $Q_{x}>0$ (or more exactly, the same sign of the kinetic energy of matter particles, assumed by convention to be positive), since otherwise the Hamiltonian is unbounded from below, which means particles can decay into lower and lower energy states, without limit, generating so-called ghosts. Therefore, for the overall stability of the theory one requires $Q_{x},c_{x}^{2}>0$.

IV The quasi-static approximation

In what follows, we put ourselves in Fourier space. That is, we replace every perturbation variable $X(\vec{x},t)$ with a plane wave parametrized by the comoving wavevector $\vec{k}$, $X(\vec{x},t)=X_{k}(t)e^{i\vec{k}\cdot\vec{x}}$. Since we deal only with linearized equations, this simply means replacing every perturbation variable or their time-derivative with its corresponding Fourier coefficient $X_{k}$ or its time derivative $\dot{X}_{k}$, and every space derivative $\partial_{i}^{(n)}$of order $n$ with $(ik_{i})^{n}X_{k}$. We drop from now on the $k$ subscripts. We then assume that the so-called *quasi-static approximation *(QSA) is valid for the evolution of perturbations. This implies that we are observing scales well inside the cosmological horizon, $\hat{k}\equiv k_{\text{}}/(aH)\gg 1$, where $k$ is the comoving wavenumber, and also inside the Jeans length of the scalar, $c_{\text{S}}\hat{k}\gg 1$, such that the terms containing $k$ (i.e., the spatial derivatives) dominate over the time-derivative terms. For the scalar field, this means we neglect its wavelike nature, and convert its Klein-Gordon differential equation into a Poisson-like constraint equation. If $c_{S}\approx 1$, the scales at which the QSA is valid correspond to all sub-horizon scales, which are also the observed scales in the recent Universe. For models with $c_{S}\to 0$, the QSA might be valid only in a narrow range of scales, or even be completely lost in the non-linear regime.

Let us explain in more detail the QSA procedure by using standard gravity as an example. Let us write down the perturbation equations for a single pressureless matter fluid in $\Lambda$CDM. From now on, we adopt the FLRW perturbed metric in the so-called longitudinal gauge, namely

[TABLE]

If we use $N=\log a$ as time variable, so that $\dot{x}=Hx^{\prime},$ the coefficients of the perturbation variables become dimensionless, and we are left with Amendola and Tsujikawa (2010)

[TABLE]

where instead of the matter density $\rho_{m}$, we use $\Omega_{m}$,

[TABLE]

and where $\delta\equiv\delta\rho_{m}/\rho_{m}$, $\theta=ik_{i}v^{i}/aH$ if $v^{i}=adx^{i}/dt$ is the peculiar velocity, so that $\theta=ik_{i}(x^{i})^{\prime}$. A glance at these equations tells us that, as an order of magnitude, $\delta\sim\theta\sim\hat{k}^{2}\Psi\sim\hat{k}^{2}\Phi$. Moreover, we assume $X\sim X^{\prime},X^{\prime\prime}$ for every perturbation variable $X=\{\delta,\theta,\Psi,\Phi\}$ (unless there is an instability, see below) and, consequently, $\hat{k}^{2}X\gg X^{\prime},X^{\prime\prime}$. Therefore for $\hat{k}\gg 1$ the equations become

[TABLE]

and one can derive the well-known second-order growth equation with dimensionless coefficients

[TABLE]

The same QSA procedure can be followed for more complicate systems. When a coupled scalar field is present, its perturbation is of the same order as the gravitational potentials, $\delta\phi\sim\Psi\sim\Phi$.

The QSA says nothing about the background behavior. Additional conditions might be imposed, for instance that the background scalar field slow rolls so that the kinetic terms, proportional to the derivatives $\phi^{\prime},\phi^{\prime\prime}$, are negligible with respect to the potential ones. This is indeed expected in order to produce an accelerated regime not too dissimilar from $\Lambda$CDM but, first, one can have acceleration driven by purely kinetic terms, and second, acceleration can be produced even with a significant fraction of energy in the kinetic terms. So slow-roll approximation and QSA should be kept well distinguished. However, in some formula below we will explicitly make use of the slow-roll approximation on top of the QSA.

Let us emphasize that the QSA applies only for classically stable systems. Imagine a scalar field obeying a second-order equation in Fourier space

[TABLE]

where from now on we use the physical wavenumber

[TABLE]

instead of the comoving one, and where $F,S$ are the friction and the source, respectively, depending in general on the background solution and on other coupled fields. If $c_{S}^{2}<0$, the solution $\phi$ will increase asymptotically as $e^{|c_{S}|k\log a}=a^{|c_{S}|k}$ and in this case $\phi^{\prime\prime}\sim|c_{S}^{2}|k^{2}\phi$ will not be negligible with respect to $c_{S}^{2}k^{2}\phi$ as we assumed in the QSA.

For simplicity, from now on we assume that the space curvature has been found to be vanishing, so $|\Omega_{k0}|\ll 1$. Using Einstein’s field equations and a pressureless perfect fluid for matter, we can derive from the HL two generalized Poisson equations in Fourier space, one for $\Phi$ and one for $\Psi$:

[TABLE]

(we remind that in our units $8\pi G_{N}=1$) where $z$ is the redshift, $k$ the physical wavenumber, $\delta_{m}$ the matter density contrast and $\eta$ and $Y$ are two functions of scale and time that parametrize deviations from standard gravity. In some papers the function $Y$ is also called $\mu$. Comparing with eqs. (12,13), we see that in Einstein’s General Relativity they reduce to $\eta=Y=1$. Clearly, the anisotropic stress $\eta$, or gravitational slip, is defined as

[TABLE]

(From now on, all the perturbation quantities are meant to be root-mean-squares of the corresponding random variables, and therefore positive definite; we can therefore define ratios like $\eta$). A value of $\eta\not=1$ can be generated in standard General Relativity only by off-diagonal spatial elements of the energy-momentum tensor. For a perturbed fluid, these elements are quadratic in the velocity, $T_{ij}\sim\rho v_{i}v_{j}$, and therefore vanish at first order for non-relativistic particles. Free-streaming relativistic particles can instead induce a deviation from $\eta=1$: this is the case of neutrinos. However, they play a substantial role only during the radiation era and are negligible today Weinberg (2004). Therefore, $\eta\not=1$ in the late universe means that gravity is modified, unless there is some hitherto unknown abundant form of hot dark matter.

In the QSA one can show that for HL De Felice et al. (2011); Amendola et al. (2013a)

[TABLE]

for suitably defined functions $h_{1-5}$ of time alone that depend only on $K,G_{3,4,5}$. Their full form will be given in Sec. VI along with another popular parametrization of the HL equations proposed in Bellini and Sawicki (2014). In general, the functions $h_{3,4,5}$ are proportional to $\mu^{-2}$, where $\mu$ is a mass scale. In the simplest cases $\mu$ corresponds to the standard mass $m$, i.e. the second derivative of the scalar field potential, plus other terms proportional to $\phi^{\prime}$ or $\phi^{\prime\prime}$. These kinetic terms are expected to be subdominant if $\phi$ drives acceleration today or, more in general, during an evolution that is not strongly oscillating, so often we can assume that $h_{3,4,5}$ scale simply as $m^{-2}$. This approximation will be adopted in the explicit expression for $f(R)$ and conformal coupling that are given below. If the scalar field drives acceleration one expects $m$ to be very small, of order $H_{0}\approx 10^{-33}$eV. In this case, at the observable sub-horizon scales, $\eta\to h_{2}h_{4}/h_{5}$ and $Y\to h_{1}h_{5}/h_{3}$. If instead this scale is of the order of the linear scales that can be directly observed (e.g. 100 Mpc), then one could observationally detect the $k$-dependence of $Y,\eta$ and find that at sufficiently large scale, such that $k\ll m$, $\eta\to h_{2},Y\to h_{1}$.

The same form of $Y,\eta$ can be obtained also in other theories not based on scalars that produce second-order equations of motion, namely, in bimetric models Könnig and Amendola (2014) and in vector models De Felice et al. (2016).

It is worth stressing the fact that the time dependence of $Y,\eta$, expressed by the functions $h_{1-5}$, is essentially arbitrary. Given observations at several epochs, one can always design a HL that exactly fits the data, no matter how precise they are. In contrast, the space dependence, which in Fourier space becomes the $k$ dependence, is very simple and fixed. The reason is that the HL equations are by definition second-order, and therefore contain at most factors of $k^{2}$. The $k$ dependence is therefore potentially a more robust test for the validity of the HL than the time one. A model with two coupled scalar fields would instead generate for $Y,\eta$ a ratio of polynomials of order $k^{4}$ (see e.g. Vardanyan and Amendola (2015)). Clearly, one has to remember that all this is valid at linear scales: if the $k$ dependence is important only at non-linear scales, e.g. for $k>1$ Mpc*-1*, then it might be completely lost.

Another equivalent form that we will employ often is

[TABLE]

where

[TABLE]

This form has a simple physical interpretation. $Y\eta$ is the modifier of the $\Phi$-Poisson equation, just as $Y$ is the modifier of the $\Psi$-Poisson equation. The parameters $\alpha_{t},\alpha_{s}$ are the strengths of the fifth-force mediated by the scalar field for $\Psi$ (the metric time-time perturbed component) and for $\Phi$ (the metric space-space perturbed component), respectively. Finally, $m$ is the effective mass of the scalar field, and $\lambda\equiv 1/m$ its spatial range. This interpretation will be discussed in the next section.

A particularly simple case is realized with the $f(R)$ models, where $f(R)$ is the function of the curvature $R$ that is to be added to the Einstein-Hilbert Lagrangian. In this case in fact

[TABLE]

The derivative $f_{,R}$ is often negligible at the present epoch, in order to reproduce a viable cosmology. In this case, $m^{2}_{R}=(3f_{,RR})^{-1}$ and, for large $k$, $\eta\to 1/2$ and $Y\to 4/3$, regardless of the specific $f(R)$ model.

Another simple case is conformal scalar-tensor theory, with $G_{3}=G_{5}=0$, $G_{4}=F(\phi)/2$, and $K=(1-3\alpha_{t})F(\phi)X-V(\phi)$, where $\alpha_{t}=(F_{,\phi}/F)^{2}/2$. In this form, the strength of the fifth force is $\alpha_{t}$. In this case we have

[TABLE]

where $M^{2}=V_{,\phi\phi}$. When $M$ is vanishingly small, $\eta\to 1-2\alpha_{t}/(1+\alpha_{t})$ and $Y=(1+\alpha_{t})/F$. Comparing with eq. (26), we see that for $f(R)$, $\alpha_{t}=-\alpha_{s}=1/3$.

V Potentials in real space

In real space, one can derive the modified Newtonian potential for a radial mass density distribution $\rho(r)$ by inverse Fourier transformation. Let us start with eq. (24)

[TABLE]

For a non-linear static structure (e.g. the Earth or a galaxy) the local density is much higher than the background average density, so $\delta_{m}(k)=[\rho_{m}(k)-\rho_{m}]/\rho_{m}\approx\rho_{m}(k)/\rho_{m}$ where $\rho_{m}=\langle\rho_{m}(k)\rangle$ is the background density and

[TABLE]

is the Fourier transform of $\rho(r)$. The Poisson equation (21) becomes then

[TABLE]

In real space and for a radial configuration, this reads

[TABLE]

(since we use the physical $k$, now $r$ refers to the physical distance) where

[TABLE]

is the inverse Fourier transform of $Y(k)\rho_{m}(k)$ and $V$ is an arbitrary large volume that encompasses the structure. Assuming that $\Psi$ vanishes at infinity, eq. (31) has the general solution

[TABLE]

where $z=\cos\theta$ and where $\Psi_{N}$ is the standard Newtonian potential, while $\Psi_{Y}$ is the Yukawa correction proportional to $\alpha_{t}$. This can be solved for any given radial density distribution $\rho(r)$. For $m\to\infty$ (or $\alpha_{t}\to 0)$ we are back to the Newtonian case.

Let us focus now on the modified gravity part. This can be analytically integrated in some simple cases. We write

[TABLE]

where $F$ has two parts

[TABLE]

For a mass point at the origin, for instance, one has $\rho(r)=M\delta_{D}^{(3)}(r)$ where $\delta_{D}^{(3)}$ is the Dirac delta function in 3D, defined for any regular function $f(r)$ as $\int d^{3}rf(r)\delta_{D}^{(3)}(r)=4\pi\int f(r)r^{2}dr\delta_{D}^{(3)}=f(0)$, and therefore

[TABLE]

i.e. the so-called Yukawa correction. The total potential is then

[TABLE]

where we reintroduced for a moment Newton’s constant $G_{N}$. As anticipated, $\alpha_{t}$ gives the strength of the Yukawa interaction and $\lambda\equiv 1/m$ its spatial range. The prefactor $h_{1}$ renormalizes the product $G_{N}M$, so that only the product $h_{1}G_{N}M$ is then observable (beside $\alpha_{t},m$). Sometimes $h_{1}G_{N}$ is denoted $G_{\mathrm{eff}}$ because it can be seen as a renormalization of Newton’s constant.

A typical dark matter halo can be approximated by a Navarro-Frenk-White profile Navarro et al. (1996) with scale $r_{s}$ and density parameter $\rho_{0}$,

[TABLE]

In this case we have Pizzuti et al. (2017)

[TABLE]

where $\mathrm{Ei}(x)$ is the ExpIntegral function,

[TABLE]

Exactly the same procedure can be applied to the second potential $\Phi$, which obeys another Poisson equation

[TABLE]

One has now

[TABLE]

Notice that the mass $m$ is the same for $\Psi,\Phi$: there is just one boson, not two. The real-space expression for $\Phi$ for a point-mass $M$ is then identical to the one for $\Psi$ with $\alpha_{s}$ in place of $\alpha_{t}$ and $h_{1}h_{2}$ in place of $-h_{1}$,

[TABLE]

Finally, the so-called lensing potential $\psi(r)=\Psi(r)-\Phi(r)$ is responsible for the gravitational lensing of source images in the linear regime. In this regime, given an elliptical source at distance $r_{s}$ characterized by semiaxes of angular extent $\theta^{s}_{i}={(\theta^{s}_{x},\theta^{s}_{y})}$, the image we see is distorted by intervening matter into a new set of semiaxes $\theta_{j}=(\delta_{ij}+D_{ij})^{-1}\theta^{s}_{i}$ where the distortion matrix is proportional to $\psi(r)$

[TABLE]

All observations of gravitational lensing lead therefore ultimately to an estimation of $\psi(r)$. What is observed in practice is the power spectrum of ellipticities, i.e. the correlation of ellipticities of galaxies in the sky due to a non-zero $\psi(r)$ along the line of sight (see e.g. Dodelson (2003), chap. 10).

From eqs. (20,21) we see then that

[TABLE]

In our formalism, the lensing potential in real space amounts then to

[TABLE]

where

[TABLE]

Since $h_{2}$ is in general different from unity, the mass $M^{(\Psi)}\equiv h_{1}M$ one infers at infinity from the $\Psi$ potential (often called dynamical mass) is different from the mass $M^{(\Phi)}\equiv h_{1}h_{2}M$ one infers from the $\Phi$ potential or the one from the lensing combination $\Psi-\Phi$, i.e. $M^{(\Psi-\Phi)}=h_{1}(1+h_{2})M/2$ (lensing mass). These masses of course coincide in standard gravity. As we will see below, one can indeed compare observationally the estimations and extract $\eta$ by taking suitable ratios.

VI The parameters of the Yukawa correction

In Bellini and Sawicki (2014) it has been shown that the HL perturbation equations can be entirely written in terms of four functions of time only, $\alpha_{K,B,M,T}$, given as

[TABLE]

This parametrization (collectively called $\alpha_{i}$) is linked to the physical properties of the HL. Briefly, $\alpha_{T}$ expresses the deviation of the GW speed from $c$, $c_{T}^{2}=1+\alpha_{T}$; $\alpha_{K}$ is connected to the field kinetic sector; $\alpha_{B}$ to the mixing ("braiding") of the scalar and gravitational kinetic terms; $M_{\star}$ is the time-dependent effective reduced Planck mass and $\alpha_{M}$ its running. They are designed so that $\alpha_{i}=0$ for $\Lambda$CDM. They do not vanish, in general, for standard gravity with a non-$\Lambda$CDM background expansion, nor for non-standard gravity with a $\Lambda$CDM expansion. Several observational limits on these parameters in specific models have already been obtained (see e.g. Kreisch and Komatsu (2018)).

It is clear that cancellations can occur among terms belonging to different $G_{i}$ sectors. However, one should distinguish between dynamical cancellations, i.e. involving a particular background solution for $\phi(t),H(t)$, and algebraic cancellations, which only depend on some special choice for the functions $G_{i}(\phi,X)$. The former ones, if they exist at all and are not unstable, can be guaranteed only for some particular set of initial conditions, and might occur only for some period, unless the solution happens to be an attractor. The algebraic cancellations, however, are independent of the background evolution and therefore valid at all times. Therefore, usually only the second class is regarded as an interesting one.

We can now express the four coefficients introduced in eq. (8) that determine the stability of the HL as Bellini and Sawicki (2014)

[TABLE]

where

[TABLE]

and where $\xi=H^{\prime}/H$ and $\tilde{\Omega}_{m}=\frac{\rho_{m}}{3M_{\star}^{2}H^{2}}=1-\frac{\rho_{HL}}{3M_{\star}^{2}H^{2}}$ (with this last relation one can get rid of $\rho_{m}$ everywhere). Here "matter" represents all the components beside the scalar field, i.e. baryons, dark matter, neutrinos, radiation. The matter equation of state $w_{m}=\sum_{i}w_{i}\Omega_{i}/\Omega_{m}$ is then an effective value for all the matter components. Note that $c_{S}^{2}=1$ in the standard minimally coupled scalar field case $K=X-V(\phi),G_{3}=G_{5}=0$ and $G_{4}=1/2$.

The relation between the "observable" parameters $h_{1-5}$ that enter the Yukawa correction and the "physical" parameters $\alpha_{K,B,M,T}$ is

[TABLE]

where333With respect to the mass defined in Bellini and Sawicki (2014), we have $\mu^{2}=M^{2}/H^{2}$.

[TABLE]

Two remarks are in order. First, the quantity $\mu^{2}$ acts as an effective squared mass in the perturbation equation of motion for $\phi$; we need to assume therefore that it is non-negative to avoid instability below some finite value of $k$. Second, the expressions for $\alpha_{i}$ and $h_{i}$ are completely general and do not assume the QSA. The QSA is needed only when we connect the theory to observations through $Y,\eta$.

Considering now only pressureless matter, from the background equations in Appendix A we see that,

[TABLE]

where $w_{HL}=p_{HL}/\rho_{HL}$. In a $\Lambda$CDM background, $2\xi+3\Omega_{m}=0$, and $\mu^{2}$ simplifies to $\mu^{2}=-3\xi(\alpha_{M}\alpha_{B}+\alpha_{2})=-3\xi(\alpha_{M}\alpha_{B}+\alpha_{B}\xi+\alpha_{B}^{\prime}+3(\Omega_{m}-\tilde{\Omega}_{m}))$. Notice that $\alpha_{K}$ does not appear in the $h_{i}$-$\alpha_{i}$ relation: this means that the kinetic parameter $\alpha_{K}$ is not an observable in the QSA linear regime. In Sec. XIII we will discuss which combinations of $\alpha_{i}$ are really model independent (MI) observables in cosmology.

Assuming Einstein–Hilbert action for the gravitational sector and a canonical kinetic term for the scalar field, we have $M_{\star}^{2}=1$ and $\alpha_{B,M,T}=0$, so $\alpha_{1}=0$ and

[TABLE]

Therefore $h_{1}=h_{2}=1$ and

[TABLE]

so that, as per construction, $Y,\eta\to 1$.

It is worth noticing that the stability conditions $Q_{T},Q_{S},c_{S}^{2}>0$ imply $(2-\alpha_{B})\alpha_{1}+2\alpha_{2}>0$ and therefore $h_{3}>0$ if one also requires $\mu^{2}>0$. As we have seen, $\lambda=\sqrt{h_{3}}$ is the range of the fifth-force interaction, so it makes sense that it is positive definite for stable systems. In the standard Brans-Dicke model with a potential $V(\phi)$, for instance, and neglecting several subdominant kinetic terms, we have

[TABLE]

where $m_{\phi}^{2}=V_{,\phi\phi}\,$, and therefore finally

[TABLE]

where $\phi=M^{2}_{*}$ (notice that in Brans-Dicke $\phi$ has dimensions mass2 and therefore $m_{\phi}$ is dimensionless), so the fifth-force range is

[TABLE]

Assuming a $\Lambda$CDM expansion and $\alpha_{T}=0$, the conditions for stability during the matter era simplify to $\alpha_{K}>-3\alpha_{B}^{2}/2$ and

[TABLE]

Generalizing, we have that for a background parametrized by a (possibly time-dependent) EOS $w_{HL}$ and for matter with an effective $w_{m}$, one has

[TABLE]

To these stability conditions, arising from Eqs. (60) and (61), one should add the requirement that the friction term in the perturbation equations for $\delta\phi$, or equivalently, for the gravitational potentials $\Phi,\Psi$, is positive. This condition is quite milder than those from Eqs. (60) and (61). While a negative $c_{s}^{2}$, for instance, even for a short period, induces a unbounded growth for $k\to\infty$, a negative friction term typically leads to a power-law growth $a^{p}$, which might be a problem only if it lasts for too long. However, in order to obtain the friction instability condition one should carefully investigate the existence of growing modes also when the various coefficient are time-dependent and no simple criteria have been identified so far. Therefore we just quote the condition for negative friction (i.e. stability) for the gravitational waves, best obtained by writing down the equation in conformal time, since in this case the $k^{2}$ term is time-independent (provided $\alpha_{T}=const$). The condition is simply $\alpha_{M}>-2$.

From the $h_{i}-\alpha_{i}$ relations (64) we can derive the Yukawa strengths

[TABLE]

The Yukawa strength $\alpha_{t}$ is always positive, and therefore the fifth force is attractive, if $Q_{T},Q_{S},c_{S}^{2}>0$. We also notice that if $\alpha_{M}=\alpha_{T}=0$, then $\alpha_{1}=\alpha_{B}$ and the two strengths become equal, and $h_{2}=1$. Therefore $\Psi=-\Phi$ and, finally, $\eta=1$, even if both potentials do actually have a non-vanishing Yukawa correction, so that $Y\not=1$. In order for the parameters $\alpha_{M},\alpha_{T}$ to vanish, the gravity sector of the HL must be standard, $G_{4}=const,G_{5}=0$, barring the case of accidental dynamical cancellation for some particular background evolution. Therefore, we conclude that $\eta\not=1$ implies, and is implied by, modified gravity, at least when matter is represented by a perfect fluid Saltas et al. (2014). One cannot make a similar statement for $Y$. This is a crucial statement for what follows. Notice however that, as we show below, although modified gravity implies $\eta\not=1$, a value $\eta=1$ does not necessarily implies standard gravity, but only scale-free gravity, at least at the quasi-static level. In ref. Sawicki et al. (2017) it has been shown that $\eta=1$ at all scales implies indeed standard gravity.

We can draw more conclusions from eqs. (79).

• •

The two strengths $\alpha_{t},\alpha_{s}$ are equal also if $\alpha_{M}=\alpha_{T}$. In this case, $\eta=(1+\alpha_{T})^{-1}$ and has no scale dependence.

• •

The $k\to\infty$ limit of the modified gravity parameters (provided we are still in the linear regime) is

[TABLE]

This coincides with eqs. (4.9) of Bellini and Sawicki (2014). If $\alpha_{T}=0$ then

[TABLE]

It turns out that if one imposes stability, $c_{s}^{2}>0$, then $Y$ is always larger than, or equal to, $1/M_{\star}^{2}$, so that matter perturbations in Horndeski with $\alpha_{T}=0$ always grow faster, in the quasi-static regime, than any standard gravity model with the same $M_{\star}$ and the same background. It also follows that the lensing combination that appears in Eq. (52) amounts to

[TABLE]

Since the denominator has to be positive for stability, the sign of the effect on the gravitational lensing depends only on $\alpha_{M},\alpha_{B}$.

• •

The Yukawa corrections disappear completely if $\alpha_{1}=0$, i.e. for

[TABLE]

This is therefore the general condition to have a scale-free gravity, corresponding to $h_{3}=h_{4}=h_{5}$444We recently noticed that this relation was first provided in an unpublished draft by Mariele Motta in early 2016. If we also assume $\alpha_{T}=0$ and consequently $G_{4X}=G_{5}=0$ (conformal coupling) in the HL, as required by the GW speed constraints we discuss in Sec. VIII, it follows $\alpha_{B}=-2\alpha_{M}$ Lombriser and Lima (2017); Linder (2018)555Note that in Ref. Lombriser and Lima (2017) $\alpha_{B}$ is defined as our $-\alpha_{B}/2$ and

[TABLE]

which gives an algebraic cancellation for $G_{3}=-f^{\prime}(\phi)\log X$ and $G_{4}=f(\phi)$. In this particular model, the local gravity experiments would not detect a Yukawa correction even if gravity actually couples to the scalar field. Gravity becomes then scale free. The Planck mass would still vary with time, though. So in this model $\eta\to 1$ even if gravity is actually modified. Assuming a $\Lambda$CDM background, for this model to be stable, $c_{s}^{2}>0$ implies the condition $(\alpha_{B}H)^{\prime}>0$. For $\alpha_{B}$ constant or slowly-varying, the stability condition amounts to $\alpha_{B}<0$, so $\alpha_{M}>0$ and therefore $Y$, or the effective Newton’s constant, will decrease with time. A larger $Y$ in the past means faster perturbation growth for the same $\Omega_{m}$. Once again, however, since $\Omega_{m}$ is not a MI observable quantity, whether this means that perturbations grow faster than in $\Lambda$CDM or not is a model-dependent statement.

• •

From eq. (52) we find also that the lensing potential lacks a Yukawa term whenever $A_{L}=0$, defined in (54), i.e. $h_{2}\alpha_{s}+\alpha_{t}=0$, which amounts to

[TABLE]

Then we see that $A_{L}=0$ not only when $\alpha_{1}=0$, but also for $\alpha_{1}=-\alpha_{B}$. Again imposing the GW speed constraint, this becomes $\alpha_{B}=-\alpha_{M}$. On the HL functions, this implies

[TABLE]

which actually means that the $G_{3}$ sector, after an integration by parts, can be absorbed in $K(\phi,X)$. So for the conformal coupling and when the $G_{3}$ term is absent or does not depend on $X$, the lensing potential becomes simply twice the standard Newtonian potential

[TABLE]

This means radiation, being conformally-invariant,666The electromagnetic Lagrangian $\sqrt{-g}F_{\alpha\beta}g^{\beta\mu}g^{\alpha\nu}F_{\mu\nu}$ does not change for $g_{\mu\nu}\to f(\phi)g_{\mu\nu}$. does not feel the modification of gravity, except for the overall factor $h_{1}$ which, if time dependent, induces a time-dependent mass or Newton’s constant.

• •

In the same case as above, $\alpha_{B}=-\alpha_{M}$ and $\alpha_{T}=0$, one has

[TABLE]

which becomes

[TABLE]

when the kinetic component $\alpha_{K}$ is small. Similarly, $\alpha_{s}=-1/(3c_{s}^{2})$. For $c_{s}=1$ one obtains a Yukawa strength of $1/3\,(-1/3)$ for the $\Psi\,(\Phi)$ potential. This case is exactly realized for the $f(R)$ models.

• •

Finally, in the uncoupled case $\alpha_{M}=\alpha_{T}=0$, in which only the kinetic sector of the scalar field is modified, one has that $\alpha_{s}=\alpha_{t}>0$, so that there is a Yukawa correction, but $\eta=1$ at all quasi-static scales.

VII Local tests of gravity

Gravity has been tested since a long time in the laboratory and within the solar system (see e.g. Will (1993)). The generic outcome of these experiments is that Einsteinian gravity works well at all the scales that have been probed so far. In many experiments one assumes the existence of the same type of "fifth-force" Yukawa correction to the static Newtonian potential predicted by the HL model,

[TABLE]

(here we drop the subscript from $\alpha_{t}$ since we need consider only $\Psi$; moreover, any overall parameter can be absorbed in $G_{N}M$). Current limits on $\alpha$ and $\lambda=1/m$ have been obtained in a range of scales from micrometers to astronomical units. The constraints on the strength $\alpha$ obviously weakens for very small $\lambda$. To give an idea, at the smallest scales probed in laboratory, one has Adelberger et al. (2003) $|\alpha|\leq 10^{6}$ at $\lambda\sim 10^{-5}$m and $|\alpha|\leq 10^{-2}$ at $\lambda\sim 10^{-3}$m (Casimir-force experiments probe even shorter scales, but the constraints on $|\alpha|$ get correspondingly weaker). At planetary scales, one has $|\alpha|\leq 10^{-6}$ for $\lambda\sim 10^{6}$m (Earth-Moon distance), and $|\alpha|\leq 10^{-8}$ at $\lambda\sim 10^{11}$m (planetary orbits). Beyond this distance, the constraints from direct tests of the Newtonian $1/r$ potential weaken again.

However, the scalar field responsible for the Yukawa term induces also two post-Newtonian corrections to the Minkowski metric. For a mass distribution with velocity field $v_{i}(\mathbf{x},t)$ and density distribution $\rho(\mathbf{x},t)$, we define $U$ as the potential that solves the standard Poisson equation for non-relativistic particles, i.e. Will (1993)

[TABLE]

and $V_{i}$ as a velocity-weighted potential

[TABLE]

Then we can write down the parametrized post-Newtonian metric as follows

[TABLE]

(the full post-Newtonian metric includes several other terms which however are not excited by a conformally coupled scalar field, see e.g. Clifton et al. (2012)). Clearly, $\gamma=\beta=0$ produces the standard weak-field metric. Taking the extreme case of $\lambda\to\infty$, one has

[TABLE]

where $\beta_{0}=d\sqrt{\alpha(\phi)}/d\phi$. The parameter $1+\gamma$ can be seen as the local-gravity analogue of the anisotropic stress $\eta$, both being the ratio of $(g_{ii}-1)/(g_{00}+1)$ at linear level.

Local tests of gravity can therefore measure the Yukawa correction for both $\Phi,\Psi$, i.e. $\alpha_{t},\alpha_{s}$ and $\lambda$, and the ratio $\Phi/\Psi$, in a model-independent way. The parameter $|\gamma|$, for instance, is constrained to be less than $10^{-5}$ Patrignani et al. (2016), inducing a similar constraint on $\alpha$ at large scales. A similar constraint applies also to $\beta$. With such a small strength, there would hardly be any interesting effect in cosmology.

However, all these tests are performed within a limited range of scales, both spatial and temporal. Moreover, the tests are performed with (some of the) standard matter particles and not with, say, dark matter. Therefore, they are completely escaped if standard model particles do not feel modified gravity, for instance because the scalar field that carries the modification of gravity does not couple to them or because of screening effects, as we discuss next.

So far we considered only linear scales. At strongly non-linear scales, e.g. in the galaxy or in the solar system, the effects of modified gravity depend on the actual configuration of the scalar field. If such a configuration is static and homogeneous within a scale $r_{s}$, then the effects of modified gravity can be screened within $r_{s}$, since they are proportional to the variation of $\phi$. This is the so-called chameleon effect Khoury and Weltman (2004a, b). On the other hand, screening can occur also because of non-linearities in the kinetic part of the Klein-Gordon equation: this is the Vainshtein effect Vainshtein (1972); Babichev and Deffayet (2013). Finally, a third mechanism appears if the coupling $\alpha$ sets on a vanishing value in structures (high density regions), via a symmetry restoration, while being different from zero at the background (low density) Pietroni (2005); Olive and Pospelov (2008); Hinterbichler and Khoury (2010); Hinterbichler et al. (2011). In all cases, the strong deviation from standard gravity that we might see in cosmology are no longer visible by local experiments. In this sense, one can always build models that escape the local gravity constraints. This can be achieved also by assuming the baryons are completely decoupled from the scalar field.

In the light of these arguments, let us consider for instance in more detail the constraint on $G_{N}$ associated with the big bang nucleosynthesis (BBN), sometimes quoted as one of the most stringent cosmological bound. The yields of light elements during the primordial expansion depends on the baryon-to-photon constant ratio $\eta_{b}$ and on the cosmic expansion rate during nucleosynthesis, which in turn depend on $G_{N}$ at that time and on various standard model parameters. Fixing the standard model parameters and estimating $\eta_{b}$ by CMB measurements, one can find constraints on $G_{N}(t_{\textrm{BBN}})/G_{N}(t_{0})\approx 1\pm 0.2$ Copi et al. (2004) by comparing the predicted abundances with the observed ones, for instance deuterium in quasar absorption systems. This means that $G_{N}$ at nucleosynthesis was close to $G_{N}$ on Earth today. The easiest explanation, that $G_{N}$ did not vary at all or anyway less than 0.2 throughout the expansion, implies $|\alpha_{M}(t_{0})|<0.2(H_{0}T_{0})^{-1}$ (equal to $\approx 0.2$ in $\Lambda$CDM), where $T_{0}$ is the cosmic age. However, $G_{N}$ in the solar system might be screened, as we have mentioned, and therefore equal to the "bare" $G_{N}$ of standard gravity. Therefore, any model which is standard general relativity in the early Universe, like essentially all models built to explain present day’s acceleration, will automatically pass the BBN constraint. Moreover, one should notice that this constraint depends on a estimate of $\Omega_{b}h^{2}$ from CMB that assumed $\Lambda$CDM. Also, $G_{N}$ is in fact degenerate with the number $g_{*}$ of relativistic degrees of freedom at nucleosynthesis, so that the bound applies to $G_{N}g_{*}$ rather than to $G_{N}$ alone. Finally, a simultaneous change of the other standard-model parameters might weaken considerably the constraint, see Dent et al. (2007); Iocco et al. (2009).

VIII The impact of gravitational waves

The Horndeski model predict an anomalous propagation speed $c_{T}$ for gravitational waves (or rather, $c_{T}/c$), since the scalar field is coupled in a non-conformal way. As already mentioned, one has Amendola et al. (2013a); Bettoni et al. (2017),

[TABLE]

The almost simultaneous detection of GWs and the electromagnetic counterparts tells us that within $40$ Mpc (at $z\sim 0.008$) from us, GWs propagate essentially at the speed of light Abbott et al. (2017a). Since the signals arrived within $1$s difference and light took $10^{15}$ s to reach us, we have that $|c_{T}^{2}/c^{2}-1|<10^{-15}$. Such tight constraint immediately ruled out most of the scalar-tensor theories containing derivative couplings to gravity or at least those models which show this effect in the nearby universe (in cosmological scales) Lombriser and Taylor (2016); Lombriser and Lima (2017); Ezquiaga and Zumalacarregui (2017); Creminelli and Vernizzi (2017); Baker et al. (2017); Sakstein and Jain (2017). That is, we need to have $G_{5}=0$ and $G_{4,X}=0$. In other words, the surviving Lagrangian has arbitrary $K,G_{3}$ but vanishing $G_{5}$ and $X$-independent $G_{4}$. This kind of Lagrangian is just a form of Brans-Dicke gravity (plus a scalar field potential and a non-canonical kinetic term). It is also equivalent to standard gravity with matter conformally coupled to a scalar field, i.e. coupled to a metric $\hat{g}_{\mu\nu}=f(\phi)g_{\mu\nu}$. A dynamical cancellation among the terms depending on $G_{5}$ and $G_{4,X}$ appears extremely fine-tuned. A possible way out is to design a model with an attractor on which the conformal coupling holds, as proposed in Amendola et al. (2018a). In this case after the attractor is reached we measure $c_{T}=1$, but this does not have to be true in the past. Deviations from the speed of light in the past could be detected in B-mode CMB polarization Amendola et al. (2014).

The constraints on $c_{T}$ also affect directly $\eta$. From eq. (65) one has in fact

[TABLE]

Hence, the GW constraint $h_{2}=1$ implies that $\eta$ should also be equal to unity for sufficiently large scales (small $k$) Amendola et al. (2018b), i.e. it should recover its General Relativity value. The obvious exception are theories without a mass scale beside the Planck mass Nersisyan et al. (2018), in which case $\eta=h_{4}/h_{5}$ at all scales. On the other hand, no obvious GW constraint affects $Y$.

Gravitational waves might in principle measure another HL parameter, the running of the Planck mass, $\alpha_{M}$. In fact, as it has been shown for instance in Saltas et al. (2014), the GW amplitude $h$ obeys the equation

[TABLE]

Assuming $c_{T}=1$, this equation in the sub-horizon limit is solved by Amendola et al. (2017)

[TABLE]

where the prefactor is the ratio of the Planck mass values at emission and at observation, and $h_{s}$ is the standard amplitude expression that, for merging binaries, can be approximated as (see e.g. Maggiore (2008), eq. 4.189)

[TABLE]

Here, $d_{L}$ is the luminosity distance, $\mathcal{M}_{c}$ the so-called chirp mass and $f_{\rm GW}$ the GW frequency measured by the observer.

GWs in standard gravity can measure the luminosity distance $d_{L}$ because the chirp mass and the frequency can be independently measured by the interferometric signal. In modified gravity, what is really measured is therefore a GW distance Nishizawa (2018); Amendola et al. (2017); Belgacem et al. (2017)

[TABLE]

Comparing this with an optical determination of $d_{L}$ leads to a direct measurement of $M_{\star}$ at various epochs, and therefore of $\alpha_{M}$.

It is however likely that both the emission and the observation occur in heavily screened environments. In this case, $M_{\star}$ is the same at both ends, and no deviation from $d_{L}$ would be observed. If emission occurs in a partially unscreened environment, then one should see instead some deviation, although not necessarily connected to the cosmological, unscreened, value of $\alpha_{M}$.

IX Model dependence

The standard model of cosmology, $\Lambda$CDM, is amazingly simple. It consists of a flat, homogeneous and isotropic background space with perturbations that, at scales above some Megaparsec, have been evolving linearly until recently. The initial conditions for perturbations are set by the inflationary mechanism, and provide an initially linear and scale-invariant spectrum of scalar, vector, and tensor perturbations, that is, power-law spectra $k^{n_{x}}$, where $x$ stands for the three types of perturbations that can be excited in General Relativity. These are encoded in a spin-2 massless field, that mediates gravitational interactions via Einstein’s equations. The energy content is shared among relativistic particles, photons, and quasi-relativistic ones, neutrinos, pressureless “cold” dark matter particles, standard model particles (“baryons”), and a cosmological constant. The density of photons can be directly measured via the CMB temperature: it amounts to $0.005\%$; the density of neutrinos depends on their mass and is known to be less than 1% of the total content today. Therefore, today, only the last three components are important. The density of baryons can be fixed by the primordial Big Bang Nucleosynthesis Cyburt et al. (2016). Since the space curvature has been measured (although so far only in a model-dependent way) to be negligible, only a single parameter is left free, the present fraction of the total energy density in pressureless matter, $\Omega_{m0}$. The fraction in the cosmological constant is then $\Omega_{\Lambda 0}=1-\Omega_{m0}$.

With this one free parameter, the fraction of energy in the cosmological constant, $\Omega_{\Lambda}\approx 0.7$, $\Lambda$CDM fits all the current cosmological data: the cosmic microwave background (CMB), the weak lensing data, the redshift distortion data (RSD), the distance indicators (supernovae Ia, SNIa; baryon acoustic oscillations, BAO; cosmic chronometers, CCH; gravitational waves, GW).

There are actually a few discrepancies. Two in particular seem to be more robust. The first is the value of $H_{0}$ obtained through local measurements, in particular through Cepheids, $H_{0}=(73.45\pm 1.66)$ km/s/Mpc Riess et al. (2018), independent of cosmology, that deviates from the Planck Planck Collaboration et al. (2016) value obtained through an extrapolation from the last scattering epoch performed assuming $\Lambda$CDM, $H_{0}=(67.51\pm 0.64)$ km/s/Mpc Planck Collaboration et al. (2016). The second one is the level of linear matter clustering embodied in the normalization parameter $\sigma_{8}$: here again, the value from CMB ($\sigma_{8}=0.82\pm 0.014)$ Planck Collaboration et al. (2016) differs from the late-universe value delivered by weak lensing, $\sigma_{8}=0.745\pm 0.039$ Hildebrandt et al. (2017), and by RSD data Barros et al. (2018) $\sigma_{8}=0.75\pm 0.024$.

Another source of discrepancies is related to the dark matter clustering Weinberg et al. (2014). Dark matter-only simulations fail at reproducing some of the observed properties of the DM distribution. Although the inclusion of baryon physics may solve this, so far there is no conclusive statement and some of these issues may in fact be due to modification of gravity.

These conflicting results already display a basic problem of cosmological parameter estimation, namely the fact that it is very often model-dependent. The Planck satellite estimates of the cosmological parameters, from $\Omega_{m}$ to $h$, from the equation of state of dark energy $w_{0}$ to the clustering amplitude $\sigma_{8}$, can be obtained only by assuming, among others, a particular model of initial conditions (inflation) and of later evolution ($\Lambda$CDM). For instance, if we assume $w_{0}=-0.9$ instead of the cosmological constant value $w_{0}=-1$, one obtains $H_{0}\approx 65.5$ km/sec/Mpc(Planck Collaboration et al. (2016), fig. 27), outside the error range given above.

Another example of model-dependency comes from distance indicators and the dark energy equation of state. Cosmological distance indicators, whether based on SNIa, BAO or other, measure basically the comoving distance

[TABLE]

where $\Omega_{k0}=-k/H_{0}^{2}$ is the present amount of spatial curvature $k$ expressed as a fraction of total energy density. We see that $r(z)$ depends only on $H_{0},\Omega_{k0}$ and $E(z)=H(z)/H_{0}$. However, since distance indicators depend on the assumption of a standard candle or ruler or clock, whose absolute value we do not know, the absolute scale of $r(z)$, that is $H_{0}$, cannot be measured (except for "standard sirens", i.e. gravitational waves, Abbott et al. (2017b)). Assuming for simplicity that $\Omega_{k0}=0$, the only direct observable is $E(z)=H(z)/H_{0}$. If we neglect also radiation (a very good approximation for observations at redshift less than a few) and assume that beside pressureless matter we have a dark energy component with EOS $w(z)$, we have

[TABLE]

where

[TABLE]

We can then invert the relation (107) and obtain

[TABLE]

(here $E_{,z}$ means differentiation with respect to redshift). It appears then that in order to reconstruct $w(z)$ one needs to know $\Omega_{m0}$, beside $E(z)$. For instance, if the true cosmology is $\Lambda$CDM with $\Omega_{m0}=0.3,$ and we assume erroneously that $\Omega_{m0}=0.31$, we would infer $w(z=0)=-0.986$ and $w(z=1)=-0.897$, way different from the true value $-1$.

The problem is that $\Omega_{m0}$ is not a model-independent observable. Whenever an estimate of $\Omega_{m0}$ is given, e.g. from CMB or lensing or SNIa, it always depends on assuming a model. The reason is that there is no way, with phenomena based on gravity alone (clustering and velocity of galaxies, lensing, integrated Sachs-Wolfe, etc) to distinguish between various components of matter, since matter responds to gravity in a universal manner, unless one breaks the equivalence principle (see the “dark degeneracy” of Ref. Kunz (2009)). So in order to measure $w(z)$ one has to assume a model, that is, a parametrization, even with extremely precise measurements. For instance, if $w(z)=w_{0}+w_{a}(1-a)$, then we reduce the complexity to just two parameters, and a measurement of $E(z)$ at at least three different redshifts can fix simultaneously $w_{0},w_{a},\Omega_{m0}$. Without a parametrization, $w(z)$ cannot be reconstructed. With a parametrization, the result depends on the parametrization itself.

On the other hand, it is clear that we can always perform null tests on $w(z)$, as for most other cosmological parameters. That is, we can assume a specific $w(z)$, e.g. $w=-1$, and test whether it is consistent with the data. In this case, in flat space, one needs just three distance measurements at three different redshifts, since there are only two parameters, $H_{0}$ and $\Omega_{m0}$. If the system of three equations in two parameters has no solution, the $\Lambda$CDM model is falsified. While it is relatively easy to test, i.e. falsify, a model of gravity, it is much more complicate to measure the properties of gravity in a way that does not demand too many assumptions. This explains why the title of this paper mentions "measuring", and not "testing," gravity.

The rest of the paper will discuss what kind of model-independent measurements we can perform in cosmology, with emphasis on parameters of modified gravity. As it is obvious, one cannot claim absolute model independence. The point is rather to isolate clearly which are the assumptions, and see how far can one reach with a minimum amount of them. In the following, we will assume this set of assumptions:

a) the Universe is well-described by a homogeneous and isotropic geometry with small (linear) perturbations;

b) gravity is universal;

c) standard model particles behave from inflation onwards in the same way as we test them in our laboratories;

d) dark matter is “cold”.

One can replace the last statement with the assumption that we know the equation of state and sound speed of dark matter, provided it is not relativistic, and that the fluid remains barotropic, i.e. $p=p(\rho)$, as we will show later on. Unless otherwise specified, however, for the rest of this work we assume pressureless, cold dark matter.

Notice that we are not assuming any particular form of gravity, standard or otherwise: in fact, we refer to “gravity” as to one or more forces that act universally and without screening, at least beyond a certain scale. So we include in our treatment also gravity plus one or more scalar, vector and tensor fields. Later on we will use the Horndeski generalized scalar-tensor model for a specific example but the methods discussed here are not restricted to this case.

X Model-independent determination of the homogeneous and isotropic

geometry

What we observe in cosmology are redshifts and angular positions of sources. What we need to build and test models, however, are distances. Can we convert redshifts and angles into distances in a model-independent (MI) way? If this turns out not to be possible then there is no reason to continue our investigation to the perturbation level. Fortunately it appears we can.

The FLRW metric of a homogeneous and isotropic Universe in spherical coordinates is

[TABLE]

where $a(t)$ is the scale factor normalized at present time to $a(t_{0})=1$. If we measure $s,t,r$ in units of the natural scale length $H_{0}^{-1}$, the metric can be rewritten as

[TABLE]

The value of $\Omega_{k0}$ has been estimated by Planck to be extremely small, $|\Omega_{k0}|<0.004$ (Planck Collaboration et al. (2016), tab. 5, last column) but, again, this is a model-dependent estimate, and for now we consider it as a free parameter. We see than up to an overall scale, the FLRW metric depends only on $\Omega_{k0}$ and on $E(z)=\dot{a}/a$, from which $a(t)$ is obtained by inverting

[TABLE]

(where again $t$ is in units of $H_{0}^{-1}$).

Baryon-acoustic oscillations are the remnant of the primordial pressure waves propagating through the plasma of baryons and photons before their decoupling. By assumption $c)$, we assume their interaction at all times is the same as in our laboratories. Therefore, we can predict that the comoving scale of the BAO today is a constant $R$ independent of the redshift at which it is observed. For instance, in $\Lambda$CDM, $R$ (in units of $H_{0}^{-1}$) is equal to

[TABLE]

where the indexes $\gamma,b,m$ refer to radiation, baryons, and dark matter, respectively; moreover, $R_{s}(z)=(3\Omega_{b0}/4\Omega_{\gamma 0})/(1+z)$, $z_{drag}\approx 1000$ is the redshifts of the drag epoch (see the numerical formula given in Eisenstein and Hu (1998)), and $z_{eq}=2.396\times 10^{4}\Omega_{m0}h^{2}\approx 3200$ is the redshift at equivalence. The value $R$ can be used as a standard ruler: as for SNIa, we do not need to know $R$, but just to assume that it is constant. Therefore, we can search in the clustering of galaxies for such a scale, in particular by identifying a peak in the correlation function. The angle under which we observe $R$ gives us the ’transverse BAO’. In turn, this angle gives us the dimensionless angular diameter distance

[TABLE]

The correlation function however depends both on the angle between sources and on their redshift difference. That is, one can observe also a ’longitudinal BAO’ scale which, for a small redshift separation $dz$, amounts to

[TABLE]

This means that BAO can estimate at every redshift two combinations involving $E(z)$ and $\Omega_{k0}$, and therefore determine both in a MI way. Therefore the FLRW metric can in principle be reconstructed, within the range covered by BAO observations, without assumptions beside a). Clearly, SNIa and other distance indicators can contribute to the statistics, but do not offer information on alternative combinations of cosmological parameters. Once we have the FLRW metric, the redshifts and angles can be converted to distance by solving $ds=0$ . Given two sources at redshifts $z_{1},z_{2}$ separated by an angle $\theta$, their relative distance $r_{12}$ is Liske (2000)

[TABLE]

where the comoving distance $r(z)$ is defined in (106). The background geometry is then recoverable in a MI way. But this is not a test of gravity.

We move then to the next layer, perturbations.

XI Measuring gravity: the anisotropic stress

We have seen that the gravitational slip $\eta$ is defined as the ratio between the two gravitational potentials

[TABLE]

The lensing potential $\Psi-\Phi$ is the combination that exerts a force on the relativistic particles (i.e., for our purposes, light), while $\Psi$ exerts a force on non-relativistic particles (i.e., for our purposes, galaxies). The explicit form of the equation of motion for a generic particle moving with velocity $\boldsymbol{v}$ and relativistic factor $\gamma^{2}=(1-v^{2})^{-1}$ in a weak-field Minkowski metric is in fact Baldi et al. (2011)

[TABLE]

For small velocities, only the last term on the rhs survives; for relativistic velocities, only the square bracket term. This means that in order to test gravity at cosmological scales we need to combine observations of lensing and of clustering and velocity of galaxies.

The linear gravitational perturbation theory gives the growth of the matter density contrast $\delta_{m}(k,z)$ at any redshift $z$ and any wavenumber $k$, given a background cosmology and a gravity model. It is convenient to define also the growth function

[TABLE]

and the growth rate

[TABLE]

where, as usual, the prime stands for derivative with respect to $N=\log_{e}a$.

However, what we observe is the galaxy number density contrast in redshift space, usually expressed in terms of the galaxy power spectrum as a function of wavenumber $k$ and redshift, $P_{gal}(k,z)$. Since galaxies are expected to be a biased tracer of mass, we need to introduce a bias function

[TABLE]

that in general depends on time and space (that is, on $k,z$). If $b=1$ then the number density of galaxies in a given place is proportional to the amount of underlying total matter, $\rho_{gal}=const\times\rho_{m}$. If $b>1$ ($<1$) then galaxies are more (less) clustered than matter. Moreover, since we observe in redshift space, which means we observe sum of cosmic expansion and the radial component of the local peculiar velocity, to convert to real space we need the Kaiser transformation Kaiser (1987), which induces a correction factor $(1+f\mu^{2}/b)^{2}$ that depends on the cosine $\mu$ of the angle between the line of sight and the wavevector $\vec{k}$.

This means that the relation between what we observe, namely the galaxy power spectrum in redshift space, and what we need to test gravity, namely $\delta_{m}$, can be written as Seo and Eisenstein (2003)

[TABLE]

where

[TABLE]

($A,R$ are mnemonics for amplitude and redshift, respectively) where $\delta_{\text{m0}}(k)=\delta_{m}(k,0)$ is the root-mean-square matter density contrast today. With this definition, $\delta_{m0}$ is normalized as

[TABLE]

where $W_{8}(kR_{8})$ is the window function for a 8 $h^{-1}$Mpc sphere, $W(x)=3(\sin{x}-x\cos{x})/x^{3}$. Sometimes one defines $\hat{\delta}_{m0}=\delta_{m0}/\sigma_{8}$, which is then normalized to unity. $\hat{\delta}_{m0}$ can be referred to as the shape of the present power spectrum. Eq. (122) shows that $A,R$ are the only two observables one can derive from linear galaxy clustering. This dataset is often collectively called redshift distortion, RSD.

There is then a third observable that one can obtain from weak lensing. From eq. (52) we see that by estimating the shear distortion one can measure the quantity

[TABLE]

Since $E(z)$ can be estimated independently, we define another observable, to be denoted $L$ Amendola et al. (2013a), as follows

[TABLE]

Together with $E=H(z)/H_{0}$, the quantities $A,R,L$ are the only cosmological information one can directly gather at the linear level777As we have seen, also $\Omega_{k0}$ is a direct observable, but for simplicity we have assumed that is negligible at all relevant epochs.. Other observations, like the integrated Sachs-Wolfe or velocity fields, only give combinations of $A,R,L,E$, rather than new information. A direct measurement of the peculiar velocity field and its time derivative, for instance, would produce through the Euler equation (11) an estimation of the combination $V=\Omega_{m0}YG\delta_{m0}$, which however is equivalent to $2RE^{2}(2+(\log ER)^{\prime})/3(1+z)^{3}$. That is, at least at the linear level, one can add more statistics, but will always end up with these four quantities rather than, say, a direct estimate of $\Omega_{m0}$ or $Y$. A preliminary non-linear analysis Rampf et al. (2017) shows that employing higher-order statistics we can obtain more MI information, but we will not consider this here.

We can now write down the lensing equation in Fourier space in the following way (see Motta et al. (2013))

[TABLE]

where $\hat{k}=k/aH$. The linearized matter conservation equations, i.e. the continuity equation and the Euler equation, can be combined in a single second-order equation

[TABLE]

that depends only on the pressureless assumption $d$) and not on the gravitational model. In terms of our observational variables and for slowly varying potentials, this becomes

[TABLE]

These equations show clearly that lensing and matter growth can measure some combination of $R,L,E$ and their derivatives, as will be seen explicitly below. For now, let us just rewrite eq. (129), employing also eq. (21) as

[TABLE]

We see then that $Y$ is not, unfortunately, a MI quantity. Even if we have precise information on $R,E$, we would still need at any $k,z$ the combination $\Omega_{m0}/f$, which is not an observable. Only a null test of standard gravity plus a specific cosmological model, say $\Lambda$CDM, is possible: in this case in fact $Y=1$, and $f\approx\Omega_{m}^{0.55}$ are known, and we have that $\Omega_{m0}$ is uniquely measured by a combination of $R,R^{\prime},E,E^{\prime}$. Any two measures at different $k$ or $z$ must then give the same $\Omega_{m0}$. We show below that $\eta$, in contrast to $Y$, is a MI quantity.

Although $A,R,L$ might be interesting statistics on their own, our goal here is to test gravity. Now, the bias function depends on complicate, possibly non-linear and hydrodynamical processes, so that, even if $b$ depends on gravity, we do not know how. Also the shape $\delta_{m0}$ of the power spectrum depends on initial conditions (inflation) and, possibly, on processes that distorted the initial spectrum during the cosmic evolution. In fact, even if we could exactly measure the power spectrum shape from CMB without a parametrization like $n_{s}$ or its “running”, nothing prevents that an unknown process, for instance the presence of early dark energy or early modified gravity, distorted the spectrum at some intermediate redshift between last scattering and today. Therefore, in order to obtain model-independent measures, we should get rid of both $b$ and $\delta_{m0}$. It was shown in Amendola et al. (2012); Motta et al. (2013) that one can obtain only three statistics where the effects of the shape of the primordial power spectrum is canceled out, namely

[TABLE]

In the last equation we introduced the often-employed quantity

[TABLE]

Notice that we are not defining $\sigma_{8}(k,z)$ as an integral over the power spectrum at $z$, as in eq. (124), because we are interested in the $k$-dependence. These quantities depend in general on $k,z$ in an arbitrary way. Every other ratio of $A,R,L$ or their derivatives can be obtained through $P_{1-3}$ or their derivatives.

Let us discuss the three statistics $P_{1-3}$ in turn. The first quantity, $P_{1}$, often called $\beta$ in the literature, contains the bias function. Since we do not know how to extract gravitational information, if any, from the bias, we do not consider it any longer.

Concerning $P_{3}$, we notice that, although related, what is observed is $R$ and not $f\sigma_{8}(k,z)$. In order to determine the latter from the observable $R$, one has to assume a value of $\delta_{m0}/\sigma_{8}=\hat{\delta}_{m0}$ (typically chosen to be given by $\Lambda$CDM), so that it is not a model-independent observable. One could imagine that $P_{3}$ alone is instead a direct test of gravity, since it depends only on $f$. However, in order to predict the theoretical value of $R$ (or $f$) as a function of the gravity parameter $Y$ from eq. (130) one needs to choose a value of $\Omega_{m0}$ and the initial condition $f(k,z_{i})$ at some epoch $z_{i}$ for every $k$. In almost all the papers on this topic since Lahav et al. (1991), this initial condition is assumed to be given by a purely matter-dominated universe at some high redshift (this is, for instance, how the well-known approximated formula $f\approx\Omega_{m}^{\gamma}(z)$ is obtained). However, in models of early dark energy or early modified gravity, this assumption is broken. Therefore, once again, $P_{3}$ alone cannot provide a MI measurement of gravity. Clearly, exactly as we have seen for the dark energy EOS $w(z)$, if one parametrizes $Y(k,z)$ with a sufficiently small number of free parameters, then the RSD data alone, which provide $R(k,z)$, can fix both $\Omega_{m0}$ and $Y(k,z)$.

We can also see that $P_{2}$ is trivially related to the $E_{G}$ statistics, whose expected value at a scale $k$ is (see Leonard et al. (2015) and references therein) as

[TABLE]

In $\Lambda$CDM and with Planck 2015 parameters, its present value is $E_{g0}\approx\Omega_{m0}/f_{0}\approx 0.58$. With our definitions, the relation with $P_{2}$ is given by

[TABLE]

The $E_{g}$ statistics has been used several times as a test of modified gravity Zhang et al. (2007); Reyes et al. (2010); de la Torre et al. (2017); Leonard et al. (2015). However, it is not per se a model-independent test. In fact, the theoretical value of $E_{g}$ depends on $\Omega_{m0}$ and on $f$. As already stressed, $\Omega_{m0}$ is not an observable quantity. Moreover, the growth rate $f$ is estimated by solving the differential equation of the perturbation growth and this requires initial conditions and $Y$. As a consequence of this, when we compare $E_{g}$ to the predicted value (132), we can never know whether any discrepancy is due to a different value of $\Omega_{m0}$ or different initial conditions, or non-standard modified gravity parameters $Y,\eta$. As previously, one can employ $E_{g}$ only to perform a null test of standard gravity plus $\Lambda$CDM, or other specific models. This is of course a task of primary importance, but is different from measuring the properties of gravity in a model-independent way.

In contrast, we can define a MI statistics to measure gravity, in particular the parameter $\eta,$ by combining the equation for the growth of structure formation eq. (129) with the lensing equation (127), and with eqs. (21,14). We see then that the gravitational slip as a function of model-independent observables is given by

[TABLE]

In order to distinguish the observables from the theoretical expectations, we denoted the combination on the left-hand-side of this equation as $\eta_{\textrm{obs}}$. The statistics $\eta_{\textrm{obs}}$ is model-independent because it estimates directly $\eta$ without any need to assume a model for the bias, nor to guess $\sigma_{8}$ or $\Omega_{m0}$, nor to assume initial conditions for $f$. So if observationally one finds $\eta_{\textrm{obs}}\not=$1, then $\Lambda$CDM and all the models in standard gravity and in which dark energy is a perfect fluid are ruled out. As a consequence, cautionary remarks like those in Amon et al. (2017), namely that their results about $E_{g}$ cannot be employed until the tension between $\Omega_{m0}$ in different observational dataset is resolved, do not apply to $\eta_{\textrm{obs}}$. The price to pay is that eq. (137) depends on derivatives of $E$ and, through $P_{3}$, of $f\sigma_{8}(z)$. Derivatives of random variables are notoriously very noisy. In the next sections we will compare several methods to extract the signal.

If we abandon the linear regime then of course new observables can be devised, see e.g. Rampf et al. (2017). One interesting case is provided by relaxed galaxy clusters, for which we can reasonable expect that the virial theorem is at least approximately respected. In this case, we can directly measure the potential $\Psi$ by the Jeans equation, i.e. the equilibrium equation between the motion of the member galaxies and the gravitational force (note that the potential remains linear for galaxies and clusters, even for a non-linear distribution of matter). The lensing potential can instead be mapped through weak and strong lensing of background galaxies. In this case, one can gather much more information on the modified gravity parameters than in the linear regime Pizzuti et al. (2017). However, the validity of this approach relies entirely on two important assumptions. First, we must assume the validity of the virial theorem, which can be more or less reasonable, but cannot be proved independently. Second, since we have access only to the radial component of the member galaxy velocities, we must assume a model for the velocity anisotropy, i.e. how the other components are distributed within the cluster.

Concluding this section, we recap and emphasize the main points. $A,R,L,E$ are the only independent linear observables in cosmology. The ratios $P_{1-3}$ are independent of the initial conditions (i.e., of the power spectrum shape). $P_{2},P_{3}$ are also independent of the galaxy bias. The combination $\eta_{\textrm{obs}}(P_{2,}P_{3,}E)$ is therefore a model-independent test of gravity: it does not depend on bias, on initial conditions, nor on other unobservable quantities like $\Omega_{m0}$ or $\sigma_{8}$. If $\eta_{\textrm{obs}}\not=1$, gravity is not Einsteinian; if $\eta_{\textrm{obs}}$ does not have the same $k^{2}$ dependence as the Horndeski theory, the entire Horndeski model is rejected. All this, of course, provided our conditions $a)-d)$ are verified.

XII General perfect fluid

What happens if we remove condition $d$), namely, that matter is pressureless? If matter is a perfect fluid and we know or hypothesize a different equation of state and sound speed, then eq. (128) is modified since the continuity and Euler equations, which come directly from the conservation of the energy-momentum tensor, now read

[TABLE]

where the sound speed is $c^{2}_{s}\equiv\delta p/\delta\rho$ and $\sigma$ is the matter anisotropic stress. Here we are assuming that $\delta$ represents the density contrast of matter, both baryons and dark matter, whose microphysical properties are described by with some effective parameters $\sigma,c_{s},w$. Assuming a zero anisotropic stress, since we are dealing with non-relativistic matter, and for small and constant $w$ and $c_{s}^{2}$, we obtain the following second order differential equation

[TABLE]

which reduces to eq. (17) for cold dark matter, where $\sigma=w=c_{s}^{2}=0$. In the case of a constant $w$, the matter would not follow an $a^{-3}$ behavior as a function of time, but it would scale with $(1+z)^{3(1+w)}$, so that the lensing equation (127) would now read

[TABLE]

Taking the appropriate ratios of the two equations above, we can obtain $\eta$ as we did for eq. 137, but this time some extra term appears

[TABLE]

where $\mathcal{W}_{1}=3(c_{s}^{2}-2w)$ and $\mathcal{W}_{2}=6(c_{s}^{2}-w)$. Both $\mathcal{W}_{1}$ and $\mathcal{W}_{2}$ reduce to zero for standard cold dark matter, such that we recover eq. (137) exactly in that case. For a barotropic fluid such that $p=p(\rho)$, $c_{s}^{2}=w$ and $\mathcal{W}_{2}=0$. In this case, we have again a MI estimator for $\eta$, provided we know $c_{s},w$. On the other hand, if $\mathcal{W}_{2}\not=0$, we see that this estimation of $\eta$ contains the growth rate $f$, which we argued not to be a model-independent observable in the linear regime. However, an extension of this formalism to the quasilinear scales Rampf et al. (2017) has shown that $f$ can indeed be recovered in a model-independent way, using observations of the bispectrum.

XIII The linear, scalar, quasi-static, model-independent Horndeski observables

For the previous sections we can draw a remarkable conclusion. Since $\eta$ is the only linear, quasi-static, MI cosmological observable, we see that, among the HL parameters, only the time-dependent functions $h_{2},h_{4},h_{5}$ (see eq. 23) share the same property. The GW speed constraint has already measured $h_{2}(t_{0})=1$. Assuming this can be extended at all times, so that $\alpha_{T}=0$, and assuming $H$ is also measured in a MI way, we see that what can still be measured at the linear perturbation level are the combinations

[TABLE]

that corresponds to the two scales one can measure in $\eta$. If $O_{2}$ vanishes, $h_{4}=h_{5}$ and $\eta=h_{2}=1$ as in the standard case. As we have already seen, this happens only in two cases, for $\alpha_{M}=0$ and for $\alpha_{M}=-\alpha_{B}/2$.

XIV Data

In the next sections we obtain an estimate of $\eta_{\textrm{obs}}$ using all the currently available data888This section and the next two are a summary of the following published paper, A. M. Pinho, S. Casas, and L. Amendola, Model-independent reconstruction of the linear anisotropic stress $\eta$, arXiv:1805.00025, JCAP11(2018)027. The first step is to reconstruct $E(z)$ (and therefore $E^{\prime}(z)$), $P_{2}(z)$ and $P_{3}(z)$ using the data all the currently relevant available data, shown in Fig. 1, where we also plot the $\Lambda$CDM curves of the different functions using the cosmological parameters from the TT+lowP+lensing Planck 2015 best-fits Planck Collaboration et al. (2015). A similar analysis, with the much smaller dataset then available, was carried out also in Ref. Trilleras (2015).

For the Hubble parameter measurements, we have used the most recent compilation of $H(z)$ data from Yu et al. (2017), including the measurements from Simon et al. (2004); Stern et al. (2009); Moresco et al. (2012); Moresco (2015), Baryon Oscillation Spectroscopic Survey (BOSS) Delubac et al. (2015); Font-Ribera et al. (2014); Moresco et al. (2016) and the Sloan Digital Sky Survey (SDSS) Zhang et al. (2012); Alam et al. (2016). In this compilation, the majority of the measurements was obtained using the cosmic chronometric technique. This method infers the expansion rate $dz/dt$ by taking the difference in redshift of a pair passively-evolving galaxies. The remaining measurements were obtained through the position of the BAO peaks in the power spectrum of a galaxy distribution for a given redshift. For this case, the measurements from Delubac et al. (2015) and Font-Ribera et al. (2014) are obtained using the BAO signal in the Lyman-$\alpha$ forest distribution alone or cross correlated with Quasi-Stellar Objects (QSO) (for the details of the method, we refer the reader to the original papers). Ref. (Alam et al., 2016) provides the covariance matrix of three $H(z)$ measurements from the radial BAO galaxy distribution. To this compilation we add the results from WiggleZ Blake et al. (2012). In addition to these, we use the recent results from Riess et al. (2017) where a compilation of Type Ia Supernovae from CANDELS and CLASH Multi-cycle Treasury programs were analyzed providing a few tight measurements of the expansion rate $E(z)$.

The $E_{g}$ data include the results from KiDS+2dFLenS+GAMA Amon et al. (2017), i.e, a joint analysis of weak gravitational lensing, galaxy clustering and redshift space distortions. We also include image and spectroscopic measurements of the Red Cluster Sequence Lensing Survey (RCSLenS) Blake et al. (2016) where the analysis combines the the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS), the WiggleZ Dark Energy Survey and the Baryon Oscillation Spectroscopic Survey (BOSS). Finally the work of VIMOS Public Extragalactic Redshift Survey (VIPERS) de la Torre et al. (2016) is also accounted for in our data. The latter reference uses redshift-space distortions and galaxy-galaxy lensing.

These sources provide measurements in real space within the scales $3<R_{p}<60h^{-1}$Mpc and in the linear regime, which is the one we are interested in. They have been obtained over a relatively narrow range of scales $\lambda$ meaning that we can consider them relative to the $k=2\pi/\lambda$-th Fourier component, as a first approximation. In any case, the discussion about the $k$-dependence of $\eta$ is beyond the scope of this work, so the final result can be seen as an average over the range of scales effectively employed in the observations. Moreover, in the estimation of $E_{g}$, based on Leonard et al. (2015), one assumes that the redshift of the lens galaxies can be approximated by a single value. With these approximations, indeed $E_{g}$ is equivalent to $P_{2}/2$, otherwise $E_{g}$ represents some sort of average value along the line of sight. We caution that these approximations can have a systematic effect both on the measurement of $E_{g}$ and on our derivation of $\eta$. In a future work we will quantify the level of bias possibly introduced by these approximations in our estimate.

Finally, the quantity $f\sigma_{8}(z)$ is connected to the $P_{3}$ parameter. Our data include measurements from the 6dF Galaxy Survey Beutler et al. (2012), the Subaru FMOS galaxy redshift survey (FastSound) Okumura et al. (2015), WiggleZ Blake et al. (2012), VIMOS-VLT Deep Survey (VVDS) Song and Percival (2008), VIMOS Public Extragalactic Redshift Survey (VIPERS) de la Torre et al. (2016); Hawken et al. (2016); de la Torre et al. (2013); Mohammad et al. (2017) and the Sloan Digital Sky Survey (SDSS) Howlett et al. (2015); Samushia et al. (2012); Tojeiro et al. (2012); Chuang and Wang (2013); Alam et al. (2016); Gil-Marín et al. (2016, 2016); Chuang et al. (2016). The values from Cabré and Gaztañaga (2009) and Guzzo et al. (2008) will not be considered since the $f\sigma_{8}(z)$ value is not directly reported.

XV Reconstruction of functions from data

The only difficulty in obtaining $\eta_{\textrm{obs}}$ is that we need to take the ratios $P_{2},P_{3}$ at the same redshift, while we have datapoints at different redshifts, and that we need to take derivatives of $E(z)$ and $f\sigma_{8}(z)$. This essentially means we need to have a reliable way to interpolate the data to reconstruct the underlying behavior.

There is no universally accepted method to interpolate data. Depending on how many assumptions one makes regarding the theoretical model, e.g. whether the reconstructed functions need just to be continuous, or smooth, depending on few or many parameters, etc., one gets unavoidably different results, especially in the final errors. Here, we consider and compare three methods to obtain the value of $\eta_{\textrm{obs}}$: binning, Gaussian Process (GP), and generalized linear regression.

The first, and simplest, method assembles the data into bins. This consists in dividing the data into particular redshift interval (bin) and for each of these intervals one calculates the average value of the subset of the data contained in that bin. The corresponding redshift and error of each bin are computed as weighted averages.

Another way to reconstruct a continuous function from a dataset is using a Gaussian Process algorithm as explained in Rasmussen and Williams (2006). This process can be regarded as the generalization of Gaussian distributions to function space since it provides a distribution over functions instead of a distribution of a random variable. Considering a dataset $\mathcal{D}=\{(x_{i},y_{i})|i=1,...n\}$, where $x_{i}$ are deterministic variables and $y_{i}$ random variables, the goal is to obtain a continuous function $f(x)$ that best describes the dataset. A function $f$ evaluated at a point $x$ is a Gaussian random variable with mean $\mu(x)$ and variance $\mathrm{Var}(x)$. The $f(x)$ values depend on the function value evaluated at other $\bar{x}$ point (particularly if they are close points). The relation between these can be given by a covariance function $\mathrm{cov}(f(x),f(\bar{x}))=k(x,\bar{x})$. The covariance function $k(x,\bar{x})$ is in principle arbitrary. Since we are interested in reconstruct the derivative of data, a Gaussian covariance function as

[TABLE]

is the chosen function since it is the most common having the least number of additional parameters. This function depends on the hyperparameters $\sigma_{f}$ and $\ell$ that allow to set the strength of the covariance function. These hyperparameters can be regarded as the typical scale and change in the $x$ and $y$ direction. The full covariance function takes the data covariance matrix $C$ into account by $M(x,\bar{x})=k(x,\bar{x})+C$. The log likelihood is then

[TABLE]

where $|M|$ is the determinant of $M(x_{i},x_{j})$. The distribution eq. (146) is usually sharply peaked and so we maximize the distribution to optimize the hyperparameters, although this is an approximation to the marginalization process and it may not be the best approach for all datasets. We employ the Python publicly available GaPP code from Seikel et al. (2012) Seikel et al. (2012).

As a third method, we use a generalized linear regression. Let us assume we have $N$ data $y_{i}$, one for each value of the independent variable $x_{i}$ and that

[TABLE]

where $e_{i}$ are errors (random variables) which are assumed to be distributed as Gaussian variables. Here $f_{i}$ are theoretical functions that depend linearly on a number of parameters $A_{\alpha}$

[TABLE]

where $g_{i\alpha}(x_{i})$ are functions of the variable $x_{i}$, chosen to be simple powers, $g_{i\alpha}=x_{i}^{\alpha}$, so that $f_{i}$ are polynomials of order $n$.

The order of the polynomial is in principle arbitrary, up to the number $N$ of datapoints. However, it is clear that with too many free parameters the resulting $\chi^{2}$ will be very close to zero, that is, statistically unlikely. At the same time, too many parameters also render the numerical Fisher matrix computationally unstable (producing, e.g., a non-positive definite matrix) and the polynomial wildly oscillating. On the other hand, too few parameters restrict the allowed family of functions. Therefore, we select the order of the polynomial function by choosing the degree for which the reduced chi-squared $\chi_{red}^{2}=\frac{\chi_{min}^{2}}{N-n-1}$, is closest to unity and such that the Fisher matrix is positive definite.

XVI Results

Let us now discuss the results of the final observable $\eta_{{\rm obs}}$ for each of these methods. The binning method contains the least number of assumptions compared to the polynomial regression or the Gaussian Process method. It is essentially a weighted average over the data points and its error bars at each redshift bin. Since we need to take derivatives in order to calculate $P_{3}$ and $E^{\prime}$, and we have few data points, we opt to compute finite difference derivatives. This has the caveat that it introduces correlations among the errors of the function and its derivatives, that we cannot take into account with this simple method. Moreover, for the binning method, we do not take into account possible non-diagonal covariance matrices for the data, which we do for polynomial regression and the Gaussian Process reconstruction.

Figure 2 shows the reconstructed functions obtained by the binning method, the Gaussian Process and with polynomial regression, alongside with the theoretical prediction of the standard $\Lambda$CDM model. In all cases the error bars or the bands represent the $1\sigma$ uncertainty.

With the binning method, the number of bins is limited by the maximum number of existing data redshifts from the smallest data set corresponding to one of our model-independent observables. In this case, this is the quantity $E_{g}$, for which we have effectively only three redshift bins. There are nine $E_{g}$ data points, but most of them are very close to each other in redshift, due to being measured by different collaborations or at different scales in real space for the same $z$. As explained in the data section above, we just regard this data as an average over different scales, assuming that non-linear corrections have been correctly taken into account by the respective experimental collaboration. Since we do not have to take derivatives of $E_{g}$, or equivalently $P_{2}$, this leaves us with three possible redshift bins, centered at $z_{1}=0.294$, $z_{2}=0.580$ and $z_{3}=0.860$, all of them with an approximate bin width of $\Delta z\approx 0.29$. At these redshifts we obtain $\eta_{{\rm obs}}(z_{1})=0.48\pm 0.45$, $\eta_{{\rm obs}}(z_{2})=-0.03\pm 0.34$ and $\eta_{{\rm obs}}(z_{3})=-2.78\pm 6.84$. These values and the estimation of the intermediate model-independent quantities can be seen in Table 2.

Regarding the Gaussian Process method, we have computed the normalized Hubble function and its derivative, $E(z)$ and $E^{\prime}(z)$ with the dgp module of the GaPP code. We reconstructed the $E(z)$ and $E^{\prime}(z)$ for the redshift interval of the data using the Gaussian function as the covariance function and initial values of the hyperparameters $\theta=[\sigma_{f}=0.5,\ell_{f}=0.5]$ that later are estimated by the code. The same procedure was done for the $P_{2}(z)$ data. We obtain for $E(z)$ and $E^{\prime}(z)$ functions the hyperparameters $\sigma_{f}=2.12$ and $\ell_{f}=2.06$ and for the $P_{2}$ function, $\sigma_{f}=0.58$ and $\ell_{f}=0.67$.

For the $P_{3}(z)$ observable, the hyperparameters obtained by the GaPP code led to a very flat and unrealistic reconstruction, that suggested us to take another approach for obtaining the optimal hyperparameters. We sampled the logarithm of the marginal likelihood on a grid of hyperparameters $\sigma_{f}$, $\ell_{f}$ from 0.01 to 2, setting this way a prior with the redshift range of the dataset, and 300 points equally separated in log-space for each dimension. Remember that the hyperparameter $\ell_{f}$ constrains the typical scale on the independent variable $z$. Thus, as an additional prior, we impose that $\ell_{f}$ needs to be smaller than the redshift range of the data, which was not guaranteed by the default GaPP code. Then we chose the pair of hyperparameters corresponding to the maximum of the log-marginal likehood. Therefore, for the $\ln(f\sigma_{8}(z))$ data, we obtain $\sigma_{f}=0.549$ and $\ell_{f}=1.361$. Its reconstructed derivative $P_{3}$ can be seen in the lower right panel of Figure 2. The function remains relatively flat, compared to the one given by other methods, but this approach has improved the determination of this observable.

Regarding the choice of the kernel function, several functions were compared, each of them with a different number of parameters to see the impact on the output. We tested the Gaussian kernel with two parameters, ($\sigma_{f},\ell_{f}$); the rational quadratic kernel with three parameters and the double Gaussian kernel with four parameters (see the original reference for the explicit implemented formula Seikel et al. (2012)). We performed tests using the $H(z)$ data obtained with the cosmic chronometer technique and the $f\sigma_{8}(z)$ data. Our tests show that the different choices shift the reconstructed function up to $6\%$ on its central value compared to the Gaussian kernel function. This happens for $H(z)$ while the effect is negligible for $f\sigma_{8}(z)$. Taking into account the above choices and procedure, we report that with the Gaussian Process method we obtain $\eta_{{\rm obs}}(z_{1})=0.38\pm 0.23$, $\eta_{{\rm obs}}(z_{2})=0.91\pm 0.36$ and $\eta_{{\rm obs}}(z_{3})=0.58\pm 0.93$.

For the polynomial regression method, we find $\eta_{{\rm obs}}(z_{1})=0.57\pm 1.05$, $\eta_{{\rm obs}}(z_{2})=0.48\pm 0.96$ and $\eta_{{\rm obs}}(z_{3})=-0.11\pm 3.21$. Note that we applied the criteria of a $\chi_{red}^{2}$ closest to one and a positive definite Fisher matrix to chose the order of the polynomial for each of the datasets. These criteria led to a choice of a polynomial of order 3 for the $E(z)$ and $E_{g}(z)$ data and order 6 for the $\ln(f\sigma_{8}(z))$ data. These polynomials can be seen in Figure 2 as solid yellow lines, together with their $1\sigma$ uncertainty bands. The higher order of the polynomial of $\ln(f\sigma_{8}(z))$ explains the "bumpiness" of the reconstruction of $P_{3}$, leading to larger errors on this observable in comparison to the GP method.

In Fig. 3 we show the reconstructed $\eta_{\rm obs}$ as a function of redshift with the three different methods, again with GP in a green dashed line, polynomial regression in a yellow solid line and the binning method in blue squares with error bars. It is possible to conclude that the methods are consistent with each other, within their $1\sigma$ uncertainties and that in most bins the results are consistent with the standard gravity scenario. We find that the error bars of the Gaussian Process reconstruction are generally smaller than the other methods, such that at the lowest redshift, GP is not compatible with $\eta_{\rm obs}=1$ at nearly $2\sigma$, while in the case of the binning method at the intermediate redshift, $z=0.58$, the tension is nearly $3\sigma$.

Finally, we can combine the estimates at three redshifts of Table 2 into a single value. Assuming a constant $\eta_{\rm obs}$ in this entire observed range and performing a simple weighted average, we find finally $\eta_{\rm obs}=0.15\pm 0.27$ (binning), $\eta_{\rm obs}=0.53\pm 0.19$ (Gaussian Process) and $\eta_{\rm obs}=0.49\pm 0.69$ (polynomial regression). The Gaussian Process method yields the smallest error and would exclude standard gravity. However, despite being sometimes advertised as “model-independent”, we believe that this method actually makes a strong assumption, since it compresses the ignorance about the reconstruction into a kernel function that depends on two or a small number of parameters, which are often not even fully marginalized over, as we did in our case. Also the binning method taken at face value would rule out standard gravity. However, as already mentioned, we did not take into account the correlation induced by the finite differences, and this might have decreased the overall error. Overall, we think the polynomial regression method is the most satisfactory one, providing the best compromise between the least number of assumptions and the best estimation of the data derivative. Therefore, we consider it as our “fiducial” result.

XVII Conclusion

Measuring the properties of gravity at large scales is one of the main tasks of cosmology for the next years. Several large observational campaigns that are underway, or will soon be Laureijs et al. (2011); Amendola et al. (2013b); Rawlings (2011); LSST Science Collaboration et al. (2009); Levi et al. (2013), will collect enough data on galaxy clustering and lensing to render this task possible to a high level of accuracy.

In order to test gravity one has to provide an alternative, either a full model or at least some parametrization that goes beyond Einstein’s gravity. Here we chose to consider the Horndeski Lagrangian because, although based on a single scalar field, displays most of the properties that make the field of modified gravity models such a rich area of research. We connected a more theoretical-oriented parameterization, the $\alpha_{i}$ parameters of Ref. Bellini and Sawicki (2014) with more phenomenologically-oriented ones, the $h_{i}$ parameters. To gain a deeper physical understanding, we discuss also some interesting limiting cases, how the Newtonian potentials look in real space, and the impact of the constraints from gravitational wave speed.

The first practical goal in cosmology is to test and possibly rule out specific models of gravity and background expansion. For instance, one can rule out $\Lambda$CDM in a number of way, the simplest of which being measuring a deviation from the predicted $H(z)$ behavior (which is not equivalent, as we have seen, to simply finding deviations from a $w=-1$ equation of state). Models in which gravity is modified can often be designed to have a perfect $\Lambda$CDM background, so it is necessary to test them at perturbation level.

Here, however, a problem arises, namely that many more assumptions need generally to be made. Some of them, listed in sec. IX, are in some sense of fundamental character, and we follow them in this work. However, most cosmological analyses that test gravity assume in addition one or more of the following assumptions: 1) that the initial conditions are given by a simple inflationary spectrum described by one or two parameters; 2) that the cosmological evolution at $z$ larger than a few is given by a pure CDM dominated Universe living in standard gravity; 3) that the linear bias depends only on time and not on scale; and 4) that the value of some parameters, like $\Omega_{m0}$ obtained from CMB analyses assuming $\Lambda$CDM, can be applied also to different models.

We have shown that a statistics called $\eta_{\textrm{obs}}$ can be measured without any of the assumptions 1-4. This statistics is an estimator of the anisotropic stress parameter $\eta$, one of the two phenomenological functions of linearized, scalar, sub-horizon modified gravity. In this sense, we say that $\eta_{\textrm{obs}}$ is (relatively) model independent. If $\eta_{\textrm{obs}}$ differs from unity, either gravity is modified, or at least one of the four "fundamental" assumptions of sec. IX are false.

We provided a preliminary estimate of $\eta_{\textrm{obs}}$ based on currently available data, $\eta_{\textrm{obs}}=0.49\pm 0.69$ in the redshift range $z=(0.2-0.8)$. The full $k$- and $z$-dependence is still unaccessible with current data. According to Amendola et al. (2013a), the Euclid mission will be able to measure $\eta$ to a few percent, so almost two orders of magnitude better than current values, and begin to put interesting limits on the $k$- and $z$-dependence. As the philosopher of science Alexandre Koyré said concerning the emergence of modern science999Koyré, A. (1948). Du monde de l’à peu près à l’univers de la précision. In A. Koyré (Ed.), Etudes d’histoire de la pensée philosophique (pp. 341-362). Paris: Gallimard. , also in measuring gravity at cosmological scales, we will finally move "from the world of approximation to the Universe of precision".

Acknowledgments

We acknowledge DFG for support through the TR33 project "The Dark Universe". We thank Alejandro Guarnizo Trilleras, Jiaming Zhao, Adrià Gómez, Guillem Domènech, Martin Kunz, Ippocratis Saltas, Mariele Motta, and Ignacy Sawicki for discussions and collaboration on these topic. We also wish to thank Jiaming Zhao for his valuable comments. DB acknowledges support from the Attraccíon del Talento Científico en Salamanca programme during the last stages of this work. A.M.P. gratefully acknowledges the support by the Landesgraduiertenförderung (LGF) grant of the Graduiertenakademie Universität Heidelberg. S.C. acknowledges support from CNES and CNRS grants.

Appendix A Background equations of motion of the Horndeski Lagrangian

The equations of motion for a flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric are De Felice et al. (2011); Bellini and Sawicki (2014)

[TABLE]

Here $\rho_{m},p_{m}$ are the conserved background energy density and the pressure of matter, respectively, and analogously ${\rho}_{HL}$ and ${p}_{HL}$ are the background energy density and pressure of the Horndeski field, defined as

[TABLE]

where $M_{\star}$ is defined in eq. (55) and $\alpha_{B}$ in eq. (58). Since in the literature there appear various definitions of the energy density associated to the scalar field, we report in table 1 the relation between the one adopted in this work and others.101010We thank Jiaming Zhao for pointing out potential issues related to this. Following Bellini and Sawicki (2014), we write the equation of motion of the scalar field $\phi$ as a (non-)conservation equation of the "shift-charge density" $n$,

[TABLE]

with

[TABLE]

and the non-conservation term, driven by a violation of the shift symmetry $\phi\to\phi+\mathrm{const}$, given by the $\phi-$derivative of $p_{HL}$

[TABLE]