# Measuring gravity at cosmological scales

Luca Amendola, Dario Bettoni, Ana Marta Pinho, Santiago Casas

► To cite this version:

Luca Amendola, Dario Bettoni, Ana Marta Pinho, Santiago Casas. Measuring gravity at cosmological scales. Universe, 2020, 6 (2), pp.20. 10.3390/universe6020020. hal-02058559

HAL Id: hal-02058559

https://hal.science/hal-02058559

Submitted on 27 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Measuring gravity at cosmological scales

Luca Amendola, Dario Bettoni, Ana Marta Pinho, Santiago Casas

► To cite this version:

Luca Amendola, Dario Bettoni, Ana Marta Pinho, Santiago Casas. Measuring gravity at cosmological scales. Universe, 2020, 6 (2), pp.20. 10.3390/universe6020020. hal-02058559

HAL Id: hal-02058559

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

Submitted on 27 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Review

# Measuring Gravity at Cosmological Scales

Luca Amendola<fnref n="1" />, Dario Bettoni <fnref n="1" />,<fnref n="2" />,* ID, Ana Marta Pinho <fnref n="1" /> and Santiago Casas <fnref n="3" /> ID

Received: 18 December 2019; Accepted: 20 January 2020; Published: 22 January 2020

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

check for
updates

**Abstract:** This review is a pedagogical introduction to models of gravity and how they are constrained 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 the 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.

**Keywords:** dark energy; modified gravity; large-scale structure

## 1. Introduction

Gravity is the force that shapes the overall temporal and spatial structure of the universe. There is not much need then to explain why it is important to test its validity at all scales and regimes. The substantial progress in collecting cosmological data achieved in the last couple of decades has made it possible, for the first time, to test gravity and measure its properties at astrophysical and cosmological scales. 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 [1] 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 [2]).

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: the generalized Poisson equation, Yukawa corrections to Newton's potential, the presence of anisotropic stress, changes in the gravitational wave speed, instabilities, and ghosts. Still, the HL is relatively simple in that it contains a

<footnote n="1">¹ ITP, Ruprecht-Karls-Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany; l.amendola@thphys.uni-heidelberg.de (L.A.); pinho@thphys.uni-heidelberg.de (A.M.P.)</footnote>

<footnote n="2">² Departamento de Física Fundamental, Universidad de Salamanca, E 37008 Salamanca, Spain</footnote>

<footnote n="3">³ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France; santiago.casas@cea.fr</footnote>

<footnote>* Correspondence: bettoni@usal.es</footnote>
Review

# Measuring Gravity at Cosmological Scales

Luca Amendola<fnref n="1" />, Dario Bettoni<fnref n="1" />,<fnref n="2" />,* , Ana Marta Pinho<fnref n="1" /> and Santiago Casas<fnref n="3" />

Received: 18 December 2019; Accepted: 20 January 2020; Published: 22 January 2020

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

check for
updates

**Abstract:** This review is a pedagogical introduction to models of gravity and how they are constrained 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 the 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.

**Keywords:** dark energy; modified gravity; large-scale structure

## 1. Introduction

Gravity is the force that shapes the overall temporal and spatial structure of the universe. There is not much need then to explain why it is important to test its validity at all scales and regimes. The substantial progress in collecting cosmological data achieved in the last couple of decades has made it possible, for the first time, to test gravity and measure its properties at astrophysical and cosmological scales. 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 [1] 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 [2]).

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: the generalized Poisson equation, Yukawa corrections to Newton's potential, the presence of anisotropic stress, changes in the gravitational wave speed, instabilities, and ghosts. Still, the HL is relatively simple in that it contains a

<footnote n="1">¹ ITP, Ruprecht-Karls-Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany; l.amendola@thphys.uni-heidelberg.de (L.A.); pinho@thphys.uni-heidelberg.de (A.M.P.)</footnote>

<footnote n="2">² Departamento de Física Fundamental, Universidad de Salamanca, E 37008 Salamanca, Spain</footnote>

<footnote n="3">³ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France; santiago.casas@cea.fr</footnote>

<footnote>* Correspondence: bettoni@usal.es</footnote>
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 [3] severely constrains 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 [2,4]. 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_{\text{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 a partial derivative, i.e., $\partial_{\mu}\phi \equiv \phi_{,\mu}$. Additionally, $\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, while Latin indexes run over space coordinates only.

## 2. Beyond Einstein

The re-discovery of the most general scalar-tensor theory that gives second-order equations of motion, Horndeski action [5], or covariant Galileons [6], and their extensions [7–13] provides a very general framework for such theories (see [14] 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:

$$S = \int d^4 x \sqrt{-g} \sum_{i=2}^{5} \mathcal{L}_i + S_m \quad (1)$$

where $S_m$ is the action for matter fields—dark matter, baryons, and radiation—and

$$\begin{align}
\mathcal{L}_2 &= K(\phi, X), \nonumber \\
\mathcal{L}_3 &= -G_3(\phi, X)\Box\phi, \nonumber \\
\mathcal{L}_4 &= G_4(\phi, X)R + G_{4,X}[(\Box\phi)^2 - (\nabla_\mu\nabla_\nu\phi)^2], \tag{2} \\
\mathcal{L}_5 &= G_5(\phi, X)G_{\mu\nu}\nabla^\mu\nabla^\nu\phi - \frac{G_{5,X}}{6}[(\Box\phi)^3 - 3(\Box\phi)(\nabla_\mu\nabla_\nu\phi)^2 + 2(\nabla_\mu\nabla_\nu\phi)^3]. \nonumber
\end{align}$$

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. (Notice that the number of these functions cannot be reduced by field redefinitions without going beyond Horndeski action [7,15]). 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 a 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.
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 [3] severely constrains 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 [2,4]. 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_{\text{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 a partial derivative, i.e., $\partial_{\mu}\phi \equiv \phi_{,\mu}$. Additionally, $\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, while Latin indexes run over space coordinates only.

## 2. Beyond Einstein

The re-discovery of the most general scalar-tensor theory that gives second-order equations of motion, Horndeski action [5], or covariant Galileons [6], and their extensions [7–13] provides a very general framework for such theories (see [14] 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:

$$S = \int d^4 x \sqrt{-g} \sum_{i=2}^{5} \mathcal{L}_i + S_m \quad (1)$$

where $S_m$ is the action for matter fields—dark matter, baryons, and radiation—and

$$ \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}[(\Box\phi)^2 - (\nabla_\mu\nabla_\nu\phi)^2], \\ \mathcal{L}_5 &= G_5(\phi, X)G_{\mu\nu}\nabla^\mu\nabla^\nu\phi - \frac{G_{5,X}}{6}[(\Box\phi)^3 - 3(\Box\phi)(\nabla_\mu\nabla_\nu\phi)^2 + 2(\nabla_\mu\nabla_\nu\phi)^3]. \end{aligned} \quad (2) $$

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. (Notice that the number of these functions cannot be reduced by field redefinitions without going beyond Horndeski action [7,15]). 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 a 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 = \text{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 = \text{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 = \text{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 an $f(R)$ model [16], 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 obtain $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 \to \phi + c$ with $c = \text{const}$. This shift-symmetric version of the HL is connected to the covariant Galileon when the functional dependence of the $G_i$ is fixed [6] 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 [17,18]. 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 [19,20]. (See [21] for a discussion on how to exorcise Ostrogradsky ghosts in non-degenerate theories.) Of course one can in principle add several scalar fields; however, 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.

### 3. 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 (DOFs). These DOFs can be collected according to how they behave under spatial rotations, i.e., as scalars, vectors and tensors. There are then four scalars (4 DOFs), two divergence-free vectors (4 DOFs), and one traceless, divergence-free tensor (2 DOFs). However, only the tensor DOFs propagate; that is, they are subject to linearized equations of motion that are second-order in the time derivatives. The other DOFs obey constraint equations, fully determined by the matter content. This should have been expected, since a massless tensor field, like the gravitational field, has only two independent degrees of freedom.

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

$$g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu} \qquad (3)$$

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

$$ds^2 = -(1+2\Psi)dt^2 + 2B_{,i}dx^i dt + ((1+2\Phi)\delta_{ij} + 2E_{,ij})dx^i dx^j. \qquad (4)$$

Inserting this metric into the Einstein-Hilbert Lagrangian without matter and developing it to a second order, one finds the second-order action in Minkowski space:
* If $G_4 = 1/2$ and $G_5 = 0$ (it is actually sufficient $G_5 = \text{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 = \text{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 = \text{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 an $f(R)$ model [16], 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 obtain $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 \to \phi + c$ with $c = \text{const}$. This shift-symmetric version of the HL is connected to the covariant Galileon when the functional dependence of the $G_i$ is fixed [6] 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 [17,18]. 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 [19,20]. (See [21] for a discussion on how to exorcise Ostrogradsky ghosts in non-degenerate theories.) Of course one can in principle add several scalar fields; however, 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.

### 3. 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 (DOFs). These DOFs can be collected according to how they behave under spatial rotations, i.e., as scalars, vectors and tensors. There are then four scalars (4 DOFs), two divergence-free vectors (4 DOFs), and one traceless, divergence-free tensor (2 DOFs). However, only the tensor DOFs propagate; that is, they are subject to linearized equations of motion that are second-order in the time derivatives. The other DOFs obey constraint equations, fully determined by the matter content. This should have been expected, since a massless tensor field, like the gravitational field, has only two independent degrees of freedom.

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

$$g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu} \qquad (3)$$

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

$$ds^2 = -(1+2\Psi)dt^2 + 2B_{,i}dx^i dt + ((1+2\Phi)\delta_{ij} + 2E_{,ij})dx^i dx^j. \qquad (4)$$

Inserting this metric into the Einstein–Hilbert Lagrangian without matter and developing it to a second order, one finds the second-order action in Minkowski space:
$$S_g = \frac{1}{2} \int d^3 x dt [8 B_{,i}^i \dot{\Phi} + 4 \Phi_{,i} \Psi_{,i} - 4 \dot{\Phi} \dot{E}_{,i}^i + 2 \Phi_{,i}^2 - 6 \dot{\Phi}^2]. \quad (5)$$

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 $\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. Therefore, there are no scalar propagating DOFs in Einstein gravity without matter.

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

$$ds^2 = -dt^2 + (\delta_{ij} + h_{ij})dx^i dx^j \quad (6)$$

where, after imposing the traceless, divergence-less conditions, and considering a wave propagating in direction $x_3$,

$$\bar{h}_{ij} = h_{ij} - \frac{1}{2} \eta_{ij} h = \begin{pmatrix} h_+ & h_\times & 0 \\ h_\times & -h_+ & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad (7)$$

one finds that the two modes $h_\alpha = \{h_+, h_\times\}$ obey in vacuum the same gravitational wave equation, $\square 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 the absence of matter fields [18,22]

$$S = \int d^3 x dt \{Q_S[\dot{\varphi}^2 - \frac{c_S^2}{a^2}(\partial_i \varphi)^2] + \sum_{\alpha=1}^{2} Q_T[\dot{h}_\alpha^2 - \frac{c_T^2}{a^2}(\partial_i h_\alpha)^2]\} \quad (8)$$

where $\varphi$ is the scalar mode perturbation, $h_\alpha$ represents the two tensor modes, and $c_S, c_T$ represents their speed of propagation, respectively. As expected, HL has now three propagating DOFs, 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 Section 6. From the classical point of view, stability is guaranteed when $Q_x, c_x^2$ (with $x = S, T$) have the same sign. In fact, in this case, 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$.

## 4. 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/(aH) \gg 1$, where $k$ is the comoving wavenumber, and $c_S \hat{k} \gg 1$ inside the Jeans length of the scalar, 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
$$S_g = \frac{1}{2} \int d^3 x dt [8 B_{,i}^i \dot{\Phi} + 4 \Phi_{,i} \Psi_{,i} - 4 \dot{\Phi} \dot{E}_{,i}^i + 2 \Phi_{,i}^2 - 6 \dot{\Phi}^2]. \quad (5)$$

The linearly perturbed equations of motion can be obtained then by the Euler–Lagrange equations with respect to $Φ, Ψ, 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 $Φ = 0$, which then shows that $Φ$ is not a propagating DOF. The same is true for $Ψ$, since there are no time derivatives for it. As we know, in fact, the potentials $Φ, Ψ$ are determined by the matter distribution through two constraints, the Poisson equations, which do not involve time derivatives. Therefore, there are no scalar propagating DOFs in Einstein gravity without matter.

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

$$ds^2 = -dt^2 + (\delta_{ij} + h_{ij})dx^i dx^j \quad (6)$$

where, after imposing the traceless, divergence-less conditions, and considering a wave propagating in direction $x_3$,

$$\bar{h}_{ij} = h_{ij} - \frac{1}{2} \eta_{ij} h = \begin{pmatrix} h_+ & h_\times & 0 \\ h_\times & -h_+ & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad (7)$$

one finds that the two modes $h_\alpha = \{h_+, h_\times\}$ obey in vacuum the same gravitational wave equation, $\square 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 the absence of matter fields [18,22]

$$S = \int d^3 x dt \{Q_S[\dot{\varphi}^2 - \frac{c_S^2}{a^2}(\partial_i \varphi)^2] + \sum_{\alpha=1}^{2} Q_T[\dot{h}_\alpha^2 - \frac{c_T^2}{a^2}(\partial_i h_\alpha)^2]\} \quad (8)$$

where $\varphi$ is the scalar mode perturbation, $h_\alpha$ represents the two tensor modes, and $c_S, c_T$ represents their speed of propagation, respectively. As expected, HL has now three propagating DOFs, 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 Section 6. From the classical point of view, stability is guaranteed when $Q_x, c_x^2$ (with $x = S, T$) have the same sign. In fact, in this case, 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$.

## 4. 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/(aH) \gg 1$, where $k$ is the comoving wavenumber, and $c_S \hat{k} \gg 1$ inside the Jeans length of the scalar, 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 ΛCDM. From now on, we adopt the Friedmann–Lemaître–Robertson–Walker (FLRW) perturbed metric in the so-called longitudinal gauge, namely

$$ds^2 = -(1+2\Psi)dt^2 + a^2(t)(1+2\Phi)\delta_{ij}dx^i dx^j. \qquad (9)$$

If we use $N = \log a$ as a time variable, such that $\dot{x} = Hx'$, the coefficients of the perturbation variables become dimensionless, and we are left with [23]

$$\delta' = -\theta - 3\Phi' \qquad (10)$$

$$\theta' = -\left(2 + \frac{H'}{H}\right)\theta + \hat{k}^2\Psi \qquad (11)$$

$$\hat{k}^2\Psi = -\frac{3}{2}\Omega_m(\delta + 3\hat{k}^{-2}\theta) \qquad (12)$$

$$\Psi = -\Phi \qquad (13)$$

where, instead of the matter density $\rho_m$, we use $\Omega_m$,

$$\rho_m = 3H^2\Omega_m \qquad (14)$$

and where $\delta \equiv \delta\rho_m/\rho_m$, and $\theta = ik_i v^i / aH$ if $v^i = adx^i/dt$ is the peculiar velocity, such that $\theta = ik_i(x^i)'$. 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', X''$ for every perturbation variable $X = \{\delta, \theta, \Psi, \Phi\}$ (unless there is an instability, see below) and, consequently, $\hat{k}^2X \gg X', X''$. Therefore, for $\hat{k} \gg 1$, the equations become

$$\delta' = -\theta \qquad (15)$$

$$\theta' = -\left(2 + \frac{H'}{H}\right)\theta - \frac{3}{2}\Omega_m\delta, \qquad (16)$$

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

$$\delta'' + \left(2 + \frac{H'}{H}\right)\delta' - \frac{3}{2}\Omega_m\delta = 0. \qquad (17)$$

The same QSA procedure can be followed for more complicated 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 are such that the kinetic terms, proportional to the derivatives $\phi', \phi''$, are negligible with respect to the potential ones. This is indeed expected in order to produce an accelerated regime not too dissimilar from ΛCDM; however, 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. Therefore, slow-roll approximation and QSA should be kept well distinguished. However, in some formulas 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

$$\phi'' + F\phi' + c_S^2 k^2 \phi = S \qquad (18)$$
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 ΛCDM. From now on, we adopt the Friedmann–Lemaître–Robertson–Walker (FLRW) perturbed metric in the so-called longitudinal gauge, namely

$$ds^2 = -(1+2\Psi)dt^2 + a^2(t)(1+2\Phi)\delta_{ij}dx^i dx^j. \qquad (9)$$

If we use $N = \log a$ as a time variable, such that $\dot{x} = Hx'$, the coefficients of the perturbation variables become dimensionless, and we are left with [23]

$$\delta' = -\theta - 3\Phi' \qquad (10)$$

$$\theta' = -\left(2 + \frac{H'}{H}\right) \theta + \hat{k}^2 \Psi \qquad (11)$$

$$\hat{k}^2 \Psi = -\frac{3}{2} \Omega_m (\delta + 3\hat{k}^{-2}\theta) \qquad (12)$$

$$\Psi = -\Phi \qquad (13)$$

where, instead of the matter density $\rho_m$, we use $\Omega_m$,

$$\rho_m = 3H^2\Omega_m \qquad (14)$$

and where $\delta \equiv \delta\rho_m/\rho_m$, and $\theta = ik_i v^i / aH$ if $v^i = adx^i/dt$ is the peculiar velocity, such that $\theta = ik_i(x^i)'$. 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', X''$ for every perturbation variable $X = \{\delta, \theta, \Psi, \Phi\}$ (unless there is an instability, see below) and, consequently, $\hat{k}^2X \gg X', X''$. Therefore, for $\hat{k} \gg 1$, the equations become

$$\delta' = -\theta \qquad (15)$$

$$\theta' = -\left(2 + \frac{H'}{H}\right) \theta - \frac{3}{2} \Omega_m \delta, \qquad (16)$$

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

$$\delta'' + \left(2 + \frac{H'}{H}\right) \delta' - \frac{3}{2} \Omega_m \delta = 0. \qquad (17)$$

The same QSA procedure can be followed for more complicated 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 are such that the kinetic terms, proportional to the derivatives $\phi', \phi''$, are negligible with respect to the potential ones. This is indeed expected in order to produce an accelerated regime not too dissimilar from ΛCDM; however, 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. Therefore, slow-roll approximation and QSA should be kept well distinguished. However, in some formulas 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

$$\phi'' + F\phi' + c_S^2 k^2 \phi = S \qquad (18)$$
where from now on we use the physical wavenumber

$$k_{\text{phys}} = \frac{k_{\text{com}}}{a} \qquad (19)$$

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'' \sim |c_S^2|k^2\phi$ will not be negligible with respect to $c_S^2k^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$:

$$k^2\Phi = \frac{1}{2}Y(k,z)\eta(z,k)\rho_m(z)\delta_m(z,k) \qquad (20)$$

$$k^2\Psi = -\frac{1}{2}Y(k,z)\rho_m(z)\delta_m(z,k) \qquad (21)$$

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 (we remind the reader that, in our units, $8\pi G_N = 1$). In some papers, the function $Y$ is also called $\mu$. Comparing with Equations (12) and (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

$$\eta = -\frac{\Phi}{\Psi}. \qquad (22)$$

(From now on, all the perturbation quantities are meant to be root-mean-squares of the corresponding random variables and are therefore positive definite; we can therefore define ratios such as $\eta$.) A value of $\eta \neq 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 [24]. Therefore, $\eta \neq 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 [17,25]

$$\eta = h_2 \left( \frac{1 + k^2 h_4}{1 + k^2 h_5} \right), \quad Y = h_1 \left( \frac{1 + k^2 h_5}{1 + k^2 h_3} \right), \qquad (23)$$

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 Section 6 along with another popular parametrization of the HL equations proposed in [22]. 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'$ or $\phi''$. 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 the 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_2h_4/h_5$ and $Y \to h_1h_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$.
where from now on we use the physical wavenumber

$$k_{\text{phys}} = \frac{k_{\text{com}}}{a} \qquad (19)$$

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'' \sim |c_S^2|k^2\phi$ will not be negligible with respect to $c_S^2k^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$:

$$k^2\Phi = \frac{1}{2}Y(k,z)\eta(z,k)\rho_m(z)\delta_m(z,k) \qquad (20)$$

$$k^2\Psi = -\frac{1}{2}Y(k,z)\rho_m(z)\delta_m(z,k) \qquad (21)$$

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 (we remind the reader that, in our units, $8\pi G_N = 1$). In some papers, the function $Y$ is also called $\mu$. Comparing with Equations (12) and (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

$$\eta = -\frac{\Phi}{\Psi} \qquad (22)$$

(From now on, all the perturbation quantities are meant to be root-mean-squares of the corresponding random variables and are therefore positive definite; we can therefore define ratios such as $\eta$.) A value of $\eta \neq 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 [24]. Therefore, $\eta \neq 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 [17,25]

$$\eta = h_2 \left( \frac{1 + k^2 h_4}{1 + k^2 h_5} \right), \quad Y = h_1 \left( \frac{1 + k^2 h_5}{1 + k^2 h_3} \right), \qquad (23)$$

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 Section 6 along with another popular parametrization of the HL equations proposed in [22]. 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'$ or $\phi''$. 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 the 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_2h_4/h_5$ and $Y \to h_1h_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 [26] and in vector models [27].

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 an 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., [28]). 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 \text{ Mpc}^{-1}$, then it might be completely lost.

Another equivalent form that we will employ often is

$$Y\eta = h_1h_2 \left(1 + \frac{\alpha_s k^2}{m^2 + k^2}\right), \quad Y = h_1 \left(1 + \frac{\alpha_t k^2}{m^2 + k^2}\right) \qquad (24)$$

where

$$\alpha_t \equiv (h_5 - h_3)/h_3 \\ \alpha_s \equiv (h_4 - h_3)/h_3 \qquad (25)$$

$$m^2 \equiv 1/h_3.$$

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,

$$\eta = 1 - \frac{1}{2} \frac{k^2}{(3/4)m_R^2 + k^2}, \quad Y = \frac{1}{1 + f_{,R}} \left( 1 + \frac{1}{3} \frac{k^2}{m_R^2 + k^2} \right). \qquad (26)$$

The derivative $f_{,R}$ is often negligible at the present epoch, in order to reproduce a viable cosmology. In this case, $m_R^2 = (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

$$\eta = 1 - \frac{2\alpha_t F k^2}{(1 + \alpha_t) F k^2 + M^2}, \quad Y = \frac{1}{F} \left( 1 + \frac{\alpha_t k^2}{k^2 + M^2} \right) \qquad (27)$$

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 Equation (26), we see that, for $f(R), \alpha_t = -\alpha_s = 1/3$.

## 5. 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 Equation (24):

$$Y = h_1 \left( 1 + \frac{\alpha_t k^2}{m^2 + k^2} \right). \qquad (28)$$
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 [26] and in vector models [27].

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 an 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., [28]). 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 \text{ Mpc}^{-1}$, then it might be completely lost.

Another equivalent form that we will employ often is

$$Y\eta = h_1 h_2 \left( 1 + \frac{\alpha_s k^2}{m^2 + k^2} \right), \quad Y = h_1 \left( 1 + \frac{\alpha_t k^2}{m^2 + k^2} \right) \qquad (24)$$

where

$$\alpha_t \equiv (h_5 - h_3)/h_3 \\ \alpha_s \equiv (h_4 - h_3)/h_3 \qquad (25)$$

$$m^2 \equiv 1/h_3.$$

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,

$$\eta = 1 - \frac{1}{2} \frac{k^2}{(3/4)m_R^2 + k^2}, \quad Y = \frac{1}{1 + f_{,R}} \left( 1 + \frac{1}{3} \frac{k^2}{m_R^2 + k^2} \right). \qquad (26)$$

The derivative $f_{,R}$ is often negligible at the present epoch, in order to reproduce a viable cosmology. In this case, $m_R^2 = (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

$$\eta = 1 - \frac{2\alpha_t F k^2}{(1 + \alpha_t) F k^2 + M^2}, \quad Y = \frac{1}{F} \left( 1 + \frac{\alpha_t k^2}{k^2 + M^2} \right) \qquad (27)$$

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 Equation (26), we see that, for $f(R), \alpha_t = -\alpha_s = 1/3$.

## 5. 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 Equation (24):

$$Y = h_1 \left( 1 + \frac{\alpha_t k^2}{m^2 + k^2} \right). \qquad (28)$$
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

$$ \rho_m(k) = \int \rho(r) e^{i\mathbf{k}\mathbf{r}} d^3\mathbf{r} \qquad (29) $$

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

$$ k^2 \Psi = -\frac{1}{2} Y(k) \rho_m(k). \qquad (30) $$

In real space and for a radial configuration, this reads

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

where

$$ \sigma(r) = \frac{V}{(2\pi)^3} \int e^{i\mathbf{k}\mathbf{r}} Y(k) \rho_m(k) d^3\mathbf{k} \qquad (32) $$

is the inverse Fourier transform of $Y(k)\rho_m(k)$, and $V$ is an arbitrary large volume that encompasses the structure (since we use the physical $k, r$ now refers to the physical distance). Assuming that $\Psi$ vanishes at infinity, Equation (31) has the general solution

$$ \Psi(y) = -\frac{h_1}{4} \int_0^\infty dr \int_{-1}^1 dz \left( \frac{1}{|\mathbf{r}-\mathbf{y}|} + \frac{\alpha_t e^{-m|\mathbf{r}-\mathbf{y}|}}{|\mathbf{r}-\mathbf{y}|} \right) \rho(r) r^2 \qquad (33) $$

$$ = \Psi_N + \Psi_Y \qquad (34) $$

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

$$ \Psi_Y(y) = -\frac{h_1}{4} \alpha_t \int_0^\infty dr \int_{-1}^1 dz e^{-m|\mathbf{r}-\mathbf{y}|} |\mathbf{r}-\mathbf{y}| \rho(r) r^2 \qquad (35) $$

$$ = -\frac{h_1}{4} \alpha_t \int_0^\infty \rho(r) r^2 dr \int_{-1}^1 \frac{e^{-m\sqrt{r^2+y^2-2ryz}}}{\sqrt{r^2+y^2-2ryz}} dz \qquad (36) $$

$$ = -\frac{h_1}{4} \alpha_t \int_0^\infty \rho(r) r^2 dr F(y, r) \qquad (37) $$

where $F$ has two parts

$$ F_1(y, r) = \frac{e^{-m(r-y)} - e^{-m(r+y)}}{mry}, \quad r > y \qquad (38) $$

$$ F_2(y, r) = \frac{e^{m(r-y)} - e^{-m(r+y)}}{mry}, \quad r < y. \qquad (39) $$

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^3r f(r)\delta_D^{(3)}(r) = 4\pi \int f(r)r^2dr\delta_D^{(3)}(r) = f(0)$; therefore,
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

$$ \rho_m(k) = \int \rho(r) e^{i\mathbf{k}\mathbf{r}} d^3\mathbf{r} \qquad (29) $$

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

$$ k^2 \Psi = -\frac{1}{2} Y(k) \rho_m(k). \qquad (30) $$

In real space and for a radial configuration, this reads

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

where

$$ \sigma(r) = \frac{V}{(2\pi)^3} \int e^{i\mathbf{k}\mathbf{r}} Y(k) \rho_m(k) d^3\mathbf{k} \qquad (32) $$

is the inverse Fourier transform of $Y(k)\rho_m(k)$, and $V$ is an arbitrary large volume that encompasses the structure (since we use the physical $k, r$ now refers to the physical distance). Assuming that $\Psi$ vanishes at infinity, Equation (31) has the general solution

$$ \Psi(y) = -\frac{h_1}{4} \int_0^\infty dr \int_{-1}^1 dz \left( \frac{1}{|\mathbf{r}-\mathbf{y}|} + \frac{\alpha_t e^{-m|\mathbf{r}-\mathbf{y}|}}{|\mathbf{r}-\mathbf{y}|} \right) \rho(r) r^2 \qquad (33) $$

$$ = \Psi_N + \Psi_Y \qquad (34) $$

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

$$ \Psi_Y(y) = -\frac{h_1}{4} \alpha_t \int_0^\infty dr \int_{-1}^1 dz e^{-m|\mathbf{r}-\mathbf{y}|} |\mathbf{r}-\mathbf{y}| \rho(r) r^2 \qquad (35) $$

$$ = -\frac{h_1}{4} \alpha_t \int_0^\infty \rho(r) r^2 dr \int_{-1}^1 \frac{e^{-m\sqrt{r^2+y^2-2ryz}}}{\sqrt{r^2+y^2-2ryz}} dz \qquad (36) $$

$$ = -\frac{h_1}{4} \alpha_t \int_0^\infty \rho(r) r^2 dr F(y, r) \qquad (37) $$

where $F$ has two parts

$$ F_1(y, r) = \frac{e^{-m(r-y)} - e^{-m(r+y)}}{mry}, \quad r > y \qquad (38) $$

$$ F_2(y, r) = \frac{e^{m(r-y)} - e^{-m(r+y)}}{mry}, \quad r < y. \qquad (39) $$

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^3r f(r)\delta_D^{(3)}(r) = 4\pi \int f(r)r^2dr\delta_D^{(3)}(r) = f(0)$; therefore,
$$ \Psi_Y(y) = -\frac{h_1}{4} \alpha_t \int_0^\infty dr \int_{-1}^1 dz \frac{e^{-m|r-y|}}{|r-y|} \rho(r) r^2 \quad (40) $$

$$ = -\frac{h_1}{8\pi} \alpha_t M \frac{e^{-m|y|}}{|y|} \quad (41) $$

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

$$ \Psi(r) = -h_1 \frac{G_N M}{r} (1 + \alpha_t e^{-mr}) \quad (42) $$

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$, such that only the product $h_1 G_N M$ is then observable (besides $\alpha_t, m$). Sometimes $h_1 G_N$ is denoted $G_{\text{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 [29] with scale $r_s$ and density parameter $\rho_0$,

$$ \rho(r) = \frac{\rho_0}{\frac{r}{r_s}(1 + \frac{r}{r_s})^2}. \quad (43) $$

In this case, we have [30]

$$ \Psi_Y(y) = -2\pi h_1 \alpha_t \int_0^\infty \rho(r) r^2 dr F(y, r) \quad (44) $$

$$ = \frac{2\pi h_1 \alpha_t \rho_0}{y} r_s^3 [e^{-m(r_s+y)} (\text{Ei}(mr_s) - \text{Ei}[m(r_s+y)])] \quad (45) $$

$$ -e^{m(r_s+y)} \text{Ei}[-m(r_s+y)] + e^{m(r_s-y)} \text{Ei}(-mr_s)] \quad (46) $$

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

$$ \text{Ei}(x) = - \int_{-x}^{\infty} \frac{e^{-t}}{t} dt. \quad (47) $$

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

$$ k^2 \Phi = \frac{1}{2} Y \eta \rho_m \delta_m. \quad (48) $$

One has now

$$ Y\eta = h_1 h_2 \left( 1 + \frac{\alpha_s k^2}{m^2 + k^2} \right). \quad (49) $$

Notice that the mass $m$ is the same for $\Psi, \Phi$: there is only 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$,

$$ \Phi(r) = h_1 h_2 \frac{G_N M}{r} (1 + \alpha_s e^{-mr}). \quad (50) $$

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_i^s = (\theta_x^s, \theta_y^s)$, the image we see is distorted by intervening matter into a new set of semiaxes $\theta_j = (\delta_{ij} + D_{ij})^{-1} \theta_i^s$ where the distortion matrix is proportional to $\psi(r)$

$$ D_{ij} = \int_{0}^{r_s} dr' (1 - \frac{r'}{r_s}) r' \psi_{,ij}. \quad (51) $$
$$ \Psi_Y(y) = -\frac{h_1}{4} \alpha_t \int_0^\infty dr \int_{-1}^1 dz \frac{e^{-m|r-y|}}{|r-y|} \rho(r) r^2 \quad (40) $$

$$ = -\frac{h_1}{8\pi} \alpha_t M \frac{e^{-m|y|}}{|y|} \quad (41) $$

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

$$ \Psi(r) = -h_1 \frac{G_N M}{r} (1 + \alpha_t e^{-mr}) \quad (42) $$

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$, such that only the product $h_1 G_N M$ is then observable (besides $\alpha_t, m$). Sometimes $h_1 G_N$ is denoted $G_{\text{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 [29] with scale $r_s$ and density parameter $\rho_0$,

$$ \rho(r) = \frac{\rho_0}{\frac{r}{r_s}(1 + \frac{r}{r_s})^2}. \quad (43) $$

In this case, we have [30]

$$ \Psi_Y(y) = -2\pi h_1 \alpha_t \int_0^\infty \rho(r) r^2 dr F(y, r) \quad (44) $$

$$ = \frac{2\pi h_1 \alpha_t \rho_0}{y} r_s^3 [e^{-m(r_s+y)} (\text{Ei}(mr_s) - \text{Ei}[m(r_s+y)])] \quad (45) $$

$$ -e^{m(r_s+y)} \text{Ei}[-m(r_s+y)] + e^{m(r_s-y)} \text{Ei}(-mr_s)] \quad (46) $$

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

$$ \text{Ei}(x) = - \int_{-x}^{\infty} \frac{e^{-t}}{t} dt. \quad (47) $$

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

$$ k^2 \Phi = \frac{1}{2} Y \eta \rho_m \delta_m . \quad (48) $$

One has now

$$ Y\eta = h_1 h_2 \left( 1 + \frac{\alpha_s k^2}{m^2 + k^2} \right) . \quad (49) $$

Notice that the mass $m$ is the same for $\Psi, \Phi$: there is only 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$,

$$ \Phi(r) = h_1 h_2 \frac{G_N M}{r} (1 + \alpha_s e^{-mr}). \quad (50) $$

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_i^s = (\theta_x^s, \theta_y^s)$, the image we see is distorted by intervening matter into a new set of semiaxes $\theta_j = (\delta_{ij} + D_{ij})^{-1} \theta_i^s$ where the distortion matrix is proportional to $\psi(r)$

$$ D_{ij} = \int_{0}^{r_s} dr' (1 - \frac{r'}{r_s}) r' \psi_{,ij}. \quad (51) $$
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., [31], chap. 10).

From Equations (20) and (21), we see then that

$$k^2\psi = -\frac{1}{2}\gamma(k, z)(1 + \eta(z, k))\rho_m(z)\delta_m(z, k). \qquad (52)$$

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

$$\Psi(r) - \Phi(r) = -h_1(1 + h_2)\frac{G_N M}{r}(1 + A_L e^{-mr}) \qquad (53)$$

where

$$A_L = \frac{h_2\alpha_s + \alpha_t}{1 + h_2}. \qquad (54)$$

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$ that one infers from the $\Phi$ potential and 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.

## 6. The Parameters of the Yukawa Correction

In [22], 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

$$M_*^2 \equiv 2(G_4 - 2XG_{4X} + XG_{5\phi} - \dot{\phi}HXG_{5X}) \qquad (55)$$

$$HM_*^2\alpha_M \equiv (\dot{M}_*)^2 \qquad (56)$$

$$H^2 M_*^2 \alpha_K \equiv 2X (K_X + 2XK_{XX} - 2G_{3\phi} - 2XG_{3\phi X}) + \qquad (57) \\ + 12\phi X H (G_{3X} + XG_{3XX} - 3G_{4\phi X} - 2XG_{4\phi XX}) + \\ + 12X H^2 (G_{4X} + 8XG_{4XX} + 4X^2 G_{4XXX}) - \\ - 12X H^2 (G_{5\phi} + 5XG_{5\phi X} + 2X^2 G_{5\phi XX}) + \\ + 4\phi X H^3 (3G_{5X} + 7XG_{5XX} + 2X^2 G_{5XXX})$$

$$HM_*^2\alpha_B \equiv 2\dot{\phi}(XG_{3X} - G_{4\phi} - 2XG_{4\phi X}) + \qquad (58) \\ + 8XH(G_{4X} + 2XG_{4XX} - G_{5\phi} - XG_{5\phi X}) + \\ + 2\dot{\phi}XH^2(3G_{5X} + 2XG_{5XX})$$

$$M_*^2\alpha_T \equiv 2X(2G_{4X} - 2G_{5\phi} - (\ddot{\phi} - \dot{\phi}H)G_{5X}). \qquad (59)$$

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, and $\alpha_B$ to the mixing (“braiding”) of the scalar and gravitational kinetic terms; $M_*$ is the time-dependent effective reduced Planck mass, and $\alpha_M$ its running. They are designed such 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., [32]).

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 a special choice for the
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., [31], chap. 10).

From Equations (20) and (21), we see then that

$$k^2\psi = -\frac{1}{2}\gamma(k, z)(1 + \eta(z, k))\rho_m(z)\delta_m(z, k). \qquad (52)$$

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

$$\Psi(r) - \Phi(r) = -h_1(1 + h_2)\frac{G_N M}{r}(1 + A_L e^{-mr}) \qquad (53)$$

where

$$A_L = \frac{h_2\alpha_s + \alpha_t}{1 + h_2}. \qquad (54)$$

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$ that one infers from the $\Phi$ potential and 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.

## 6. The Parameters of the Yukawa Correction

In [22], 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

$$M_*^2 \equiv 2(G_4 - 2XG_{4X} + XG_{5\phi} - \dot{\phi}HXG_{5X}) \qquad (55)$$

$$HM_*^2\alpha_M \equiv (\dot{M}_*)^2 \qquad (56)$$

$$H^2 M_*^2 \alpha_K \equiv 2X (K_X + 2XK_{XX} - 2G_{3\phi} - 2XG_{3\phi X}) + \qquad (57) \\ + 12\phi XH (G_{3X} + XG_{3XX} - 3G_{4\phi X} - 2XG_{4\phi XX}) + \\ + 12XH^2 (G_{4X} + 8XG_{4XX} + 4X^2G_{4XXX}) - \\ - 12XH^2 (G_{5\phi} + 5XG_{5\phi X} + 2X^2G_{5\phi XX}) + \\ + 4\phi XH^3 (3G_{5X} + 7XG_{5XX} + 2X^2G_{5XXX})$$

$$HM_*^2\alpha_B \equiv 2\dot{\phi}(XG_{3X} - G_{4\phi} - 2XG_{4\phi X}) + \qquad (58) \\ + 8XH(G_{4X} + 2XG_{4XX} - G_{5\phi} - XG_{5\phi X}) + \\ + 2\dot{\phi}XH^2(3G_{5X} + 2XG_{5XX})$$

$$M_*^2\alpha_T \equiv 2X(2G_{4X} - 2G_{5\phi} - (\ddot{\phi} - \dot{\phi}H)G_{5X}). \qquad (59)$$

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, and $\alpha_B$ to the mixing (“braiding”) of the scalar and gravitational kinetic terms; $M_*$ is the time-dependent effective reduced Planck mass, and $\alpha_M$ its running. They are designed such 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., [32]).

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 a 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 Equation (8) that determine the stability of the HL as [22]

$$Q_S = \frac{M_*^2(2\alpha_K + 3\alpha_B^2)}{(2 - \alpha_B)^2}, \qquad (60)$$

$$c_S^2 = \frac{(2 - \alpha_B) \alpha_1 + 2\alpha_2}{2\alpha_K + 3\alpha_B^2},$$

$$Q_T = \frac{M_*^2}{8}, \qquad (61)$$

$$c_T^2 = 1 + \alpha_T$$

where

$$\alpha_1 \equiv \alpha_B + (\alpha_B - 2) \alpha_T + 2\alpha_M \qquad (62)$$

$$\alpha_2 \equiv \alpha_B \zeta + \alpha'_B - 2\zeta - 3(1 + w_m)\tilde{\Omega}_m \qquad (63)$$

and where $\zeta = H'/H$ and $\tilde{\Omega}_m = \frac{\rho_m}{3M_*^2H^2} = 1 - \frac{\rho_{HL}}{3M_*^2H^2}$ (with this last relation one can get rid of $\rho_m$ everywhere). Here, "matter" represents all the components in addition to the scalar field, i.e., baryons, dark matter, neutrinos, and 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

$$h_1 = \frac{\alpha_T + 1}{M_*^2}, \qquad (64)$$

$$h_2 = \frac{1}{\alpha_T + 1}, \qquad (65)$$

$$h_3 = \frac{1}{2H^2\mu^2}((2 - \alpha_B)\alpha_1 + 2\alpha_2) \qquad (66)$$

$$h_4 = \frac{1}{H^2\mu^2}(\alpha_1 + \alpha_2) \qquad (67)$$

$$h_5 = \frac{1}{H^2\mu^2}\left(\frac{\alpha_M + 1}{\alpha_T + 1}\alpha_1 + \alpha_2\right) \qquad (68)$$

where

$$\mu^2 \equiv -3[2\zeta^2 + \zeta' + \zeta(3 + \alpha_M)]\alpha_B - 3\zeta\alpha_2. \qquad (69)$$

(With respect to the mass defined in [22], we have $\mu^2 = M^2/H^2$.)

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$.
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 Equation (8) that determine the stability of the HL as [22]

$$Q_S = \frac{M_*^2(2\alpha_K + 3\alpha_B^2)}{(2 - \alpha_B)^2}, \qquad (60)$$

$$c_S^2 = \frac{(2 - \alpha_B) \alpha_1 + 2\alpha_2}{2\alpha_K + 3\alpha_B^2},$$

$$Q_T = \frac{M_*^2}{8}, \qquad (61)$$

$$c_T^2 = 1 + \alpha_T$$

where

$$\alpha_1 \equiv \alpha_B + (\alpha_B - 2) \alpha_T + 2\alpha_M \qquad (62)$$

$$\alpha_2 \equiv \alpha_B \zeta + \alpha'_B - 2\zeta - 3(1 + w_m)\tilde{\Omega}_m \qquad (63)$$

and where $\zeta = H'/H$ and $\tilde{\Omega}_m = \frac{\rho_m}{3M_*^2H^2} = 1 - \frac{\rho_{HL}}{3M_*^2H^2}$ (with this last relation one can get rid of $\rho_m$ everywhere). Here, "matter" represents all the components in addition to the scalar field, i.e., baryons, dark matter, neutrinos, and 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

$$h_1 = \frac{\alpha_T + 1}{M_*^2}, \qquad (64)$$

$$h_2 = \frac{1}{\alpha_T + 1}, \qquad (65)$$

$$h_3 = \frac{1}{2H^2\mu^2}((2 - \alpha_B)\alpha_1 + 2\alpha_2) \qquad (66)$$

$$h_4 = \frac{1}{H^2\mu^2}(\alpha_1 + \alpha_2) \qquad (67)$$

$$h_5 = \frac{1}{H^2\mu^2}\left(\frac{\alpha_M + 1}{\alpha_T + 1}\alpha_1 + \alpha_2\right) \qquad (68)$$

where

$$\mu^2 \equiv -3[2\zeta^2 + \zeta' + \zeta(3 + \alpha_M)]\alpha_B - 3\zeta\alpha_2. \qquad (69)$$

(With respect to the mass defined in [22], we have $\mu^2 = M^2/H^2$.)

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 the Appendix A, we see that,

$$\xi = -\frac{3}{2} - \frac{p_{HL}}{2H^2 M_*^2} = -\frac{3}{2}(1 + w_{HL}\tilde{\Omega}_{HL}) \quad (70)$$

where $w_{HL} = p_{HL}/\rho_{HL}$. In a Λ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 + 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 Section 13, we will discuss which combinations of $\alpha_i$ are 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_*^2 = 1$ and $\alpha_{B,M,T} = 0$, so $\alpha_1 = 0$ and

$$\mu^2 = -9\xi\tilde{\Omega}_{HL}(1 + w_{HL}) \quad (71)$$

$$\alpha_2 = 3\tilde{\Omega}_{HL}(1 + w_{HL}). \quad (72)$$

Therefore, $h_1 = h_2 = 1$ and

$$h_{3,4,5} = -\frac{1}{3\xi H^2} = \frac{2}{9H^2(1 + w_{HL}\tilde{\Omega}_{HL})} \quad (73)$$

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

$$\mu^2 = \frac{3\alpha_M m_\phi^2 \phi'}{3H^2} \quad (74)$$

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

$$h_3 = \frac{3 + 2\omega}{2\phi m_\phi^2} \quad (75)$$

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

$$\lambda = m^{-1} = (m_\phi M_*)^{-1} \sqrt{\frac{3 + 2\omega}{2}}. \quad (76)$$

Assuming a ΛCDM expansion and $\alpha_T = 0$, the conditions for stability during the matter era simplify to $\alpha_K > -3\alpha_B^2/2$ and

$$(2 - \alpha_B)(\alpha_B + 2\alpha_M + 3\Omega_m) + 2\alpha'_B - 6\tilde{\Omega}_m > 0. \quad (77)$$

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

$$(2 - \alpha_B)(\alpha_B + 2\alpha_M) - 3[1 + w_{HL} + (w_m - w_{HL})\tilde{\Omega}_m]\alpha_B + 2\alpha'_B + 6(w_{HL} + 1)(1 - \tilde{\Omega}_m) > 0. \quad (78)$$

To these stability conditions, arising from Equations (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 Equations (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
Considering now only pressureless matter, from the background equations in the Appendix A, we see that,

$$ \xi = -\frac{3}{2} - \frac{p_{HL}}{2H^2 M_*^2} = -\frac{3}{2}(1 + w_{HL}\tilde{\Omega}_{HL}) \quad (70) $$

where $w_{HL} = p_{HL}/\rho_{HL}$. In a Λ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 + 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 Section 13, we will discuss which combinations of $\alpha_i$ are 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_*^2 = 1$ and $\alpha_{B,M,T} = 0$, so $\alpha_1 = 0$ and

$$ \mu^2 = -9\xi\tilde{\Omega}_{HL}(1 + w_{HL}) \quad (71) $$

$$ \alpha_2 = 3\tilde{\Omega}_{HL}(1 + w_{HL}). \quad (72) $$

Therefore, $h_1 = h_2 = 1$ and

$$ h_{3,4,5} = -\frac{1}{3\xi H^2} = \frac{2}{9H^2(1 + w_{HL}\tilde{\Omega}_{HL})} \quad (73) $$

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

$$ \mu^2 = \frac{3\alpha_M m_\phi^2 \phi'}{3H^2} \quad (74) $$

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

$$ h_3 = \frac{3 + 2\omega}{2\phi m_\phi^2} \quad (75) $$

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

$$ \lambda = m^{-1} = (m_\phi M_*)^{-1} \sqrt{\frac{3 + 2\omega}{2}}. \quad (76) $$

Assuming a ΛCDM expansion and $\alpha_T = 0$, the conditions for stability during the matter era simplify to $\alpha_K > -3\alpha_B^2/2$ and

$$ (2 - \alpha_B)(\alpha_B + 2\alpha_M + 3\Omega_m) + 2\alpha'_B - 6\tilde{\Omega}_m > 0. \quad (77) $$

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

$$ (2 - \alpha_B)(\alpha_B + 2\alpha_M) - 3[1 + w_{HL} + (w_m - w_{HL})\tilde{\Omega}_m]\alpha_B + 2\alpha'_B + 6(w_{HL} + 1)(1 - \tilde{\Omega}_m) > 0. \quad (78) $$

To these stability conditions, arising from Equations (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 Equations (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 = \text{const}$). The condition is simply $\alpha_M > -2$.

From the $h_i - \alpha_i$ relations (64), we can derive the Yukawa strengths

$$ 
\begin{aligned}
\alpha_s &= \frac{h_4 - h_3}{h_3} = \frac{\alpha_1 \alpha_B}{(2 - \alpha_B) \alpha_1 + 2 \alpha_2} \\
\alpha_t &= \frac{h_5 - h_3}{h_3} = \frac{\alpha_1^2}{((2 - \alpha_B) \alpha_1 + 2 \alpha_2) (\alpha_T + 1)}.
\end{aligned} 
\qquad (79) $$

The Yukawa strength $\alpha_t$ is always positive, so 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, such that $Y \neq 1$. In order for the parameters $\alpha_M, \alpha_T$ to vanish, the gravity sector of the HL must be standard, $G_4 = \text{const}$, $G_5 = 0$, barring the case of accidental dynamical cancellation for some particular background evolution. Therefore, we conclude that $\eta \neq 1$ implies, and is implied by, modified gravity, at least when matter is represented by a perfect fluid [33]. 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 \neq 1$, a value $\eta = 1$ does not necessarily imply standard gravity, but only scale-free gravity, at least at the quasi-static level. In [34], it has been shown that $\eta = 1$ at all scales implies indeed standard gravity.

We can draw more conclusions from Equations (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

$$ Y_{\infty} = \frac{h_1 h_5}{h_3} = \frac{2 \alpha_1 (1 + \alpha_M) + \alpha_2 (1 + \alpha_T)}{M_{\star}^2 (2\alpha_2 + \alpha_1 (2 - \alpha_B))} \qquad (80) $$

$$ \eta_{\infty} = \frac{h_2 h_4}{h_5} = \frac{\alpha_1 + \alpha_2}{\alpha_1 (1 + \alpha_M) + \alpha_2 (1 + \alpha_T)}. \qquad (81) $$

This coincides with Equation (4.9) of [22]. If $\alpha_T = 0$, then

$$ Y_{\infty} = \frac{1}{M_{\star}^2} \left[ 1 + \frac{(2\alpha_M + \alpha_B)^2}{2\alpha_2 + (2\alpha_M + \alpha_B)(2 - \alpha_B)} \right] \qquad (82) $$

$$ \eta_{\infty} = 1 - \frac{(2\alpha_M + \alpha_B)\alpha_M}{(2\alpha_M + \alpha_B)(1 + \alpha_M) + \alpha_2}. \qquad (83) $$

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$, such 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 Equation (52) amounts to

$$ \Sigma \equiv Y(1 + \eta) = 2 \left[ 1 + \frac{(2\alpha_M + \alpha_B)(\alpha_B + \alpha_M)}{2\alpha_2 + (2 - \alpha_B)(2\alpha_M + \alpha_B)} \right]. \qquad (84) $$

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$.
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 = \text{const}$). The condition is simply $\alpha_M > -2$.

From the $h_i - \alpha_i$ relations (64), we can derive the Yukawa strengths

$$ 
\begin{aligned}
\alpha_s &= \frac{h_4 - h_3}{h_3} = \frac{\alpha_1 \alpha_B}{(2 - \alpha_B) \alpha_1 + 2 \alpha_2} \\
\alpha_t &= \frac{h_5 - h_3}{h_3} = \frac{\alpha_1^2}{((2 - \alpha_B) \alpha_1 + 2 \alpha_2) (\alpha_T + 1)}.
\end{aligned} 
\qquad (79) $$

The Yukawa strength $\alpha_t$ is always positive, so 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, such that $Y \neq 1$. In order for the parameters $\alpha_M, \alpha_T$ to vanish, the gravity sector of the HL must be standard, $G_4 = \text{const}$, $G_5 = 0$, barring the case of accidental dynamical cancellation for some particular background evolution. Therefore, we conclude that $\eta \neq 1$ implies, and is implied by, modified gravity, at least when matter is represented by a perfect fluid [33]. 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 \neq 1$, a value $\eta = 1$ does not necessarily imply standard gravity, but only scale-free gravity, at least at the quasi-static level. In [34], it has been shown that $\eta = 1$ at all scales implies indeed standard gravity.

We can draw more conclusions from Equations (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

$$ Y_{\infty} = \frac{h_1 h_5}{h_3} = \frac{2 \alpha_1 (1 + \alpha_M) + \alpha_2 (1 + \alpha_T)}{M_*^2 (2\alpha_2 + \alpha_1 (2 - \alpha_B))} \qquad (80) $$

$$ \eta_{\infty} = \frac{h_2 h_4}{h_5} = \frac{\alpha_1 + \alpha_2}{\alpha_1 (1 + \alpha_M) + \alpha_2 (1 + \alpha_T)}. \qquad (81) $$

This coincides with Equation (4.9) of [22]. If $\alpha_T = 0$, then

$$ Y_{\infty} = \frac{1}{M_*^2} \left[ 1 + \frac{(2\alpha_M + \alpha_B)^2}{2\alpha_2 + (2\alpha_M + \alpha_B)(2 - \alpha_B)} \right] \qquad (82) $$

$$ \eta_{\infty} = 1 - \frac{(2\alpha_M + \alpha_B)\alpha_M}{(2\alpha_M + \alpha_B)(1 + \alpha_M) + \alpha_2}. \qquad (83) $$

It turns out that, if one imposes stability, $c_S^2 > 0$, then $Y$ is always larger than, or equal to, $1/M_*^2$, such 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_*$ and the same background. It also follows that the lensing combination that appears in Equation (52) amounts to

$$ \Sigma \equiv Y(1 + \eta) = 2 \left[ 1 + \frac{(2\alpha_M + \alpha_B)(\alpha_B + \alpha_M)}{2\alpha_2 + (2 - \alpha_B)(2\alpha_M + \alpha_B)} \right]. \qquad (84) $$

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

$$ \alpha_B = 2 \frac{\alpha_T - \alpha_M}{1 + \alpha_T}. \qquad (85) $$

This is therefore the general condition to have a scale-free gravity, corresponding to $h_3 = h_4 = h_5$ (we 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 Section 8, it follows that $\alpha_B = -2\alpha_M$ [35,36] (note that in [35] $\alpha_B$ is defined as our $-\alpha_B/2$), and that

$$ G_{4\phi} = -XG_{3X}, \qquad (86) $$

which gives an algebraic cancellation for $G_3 = -f'(\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 then becomes scale free. The Planck mass would still vary with time, however. Therefore, 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)' > 0$. For $\alpha_B$ constant or slowly varying, the stability condition amounts to $\alpha_B < 0$, so $\alpha_M > 0$; 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 an MI observable quantity, whether this means that perturbations grow faster than in $\Lambda$CDM or not is a model-dependent statement.

* From Equation (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

$$ \alpha_1\alpha_B + \alpha_1^2 = 0. \qquad (87) $$

We then 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

$$ G_{3X} = 0, \qquad (88) $$

which actually means that the $G_3$ sector, after an integration by parts, can be absorbed in $K(\phi, X)$. Therefore, 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

$$ \Psi(r) - \Phi(r) = 2h_1 \frac{G_N M}{r}. \qquad (89) $$

This means that radiation, being conformally invariant (the 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

$$ \alpha_t = \frac{\alpha_M^2}{c_s^2(2\alpha_K + 3\alpha_M^2)}, \qquad (90) $$

which becomes

$$ \alpha_t = \frac{1}{3c_s^2} \qquad (91) $$
* The Yukawa corrections disappear completely if $\alpha_1 = 0$, i.e., for

$$ \alpha_B = 2 \frac{\alpha_T - \alpha_M}{1 + \alpha_T}. \qquad (85) $$

This is therefore the general condition to have a scale-free gravity, corresponding to $h_3 = h_4 = h_5$ (we 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 Section 8, it follows that $\alpha_B = -2\alpha_M$ [35,36] (note that in [35] $\alpha_B$ is defined as our $-\alpha_B/2$), and that

$$ G_{4\phi} = -XG_{3X}, \qquad (86) $$

which gives an algebraic cancellation for $G_3 = -f'(\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 then becomes scale free. The Planck mass would still vary with time, however. Therefore, 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)' > 0$. For $\alpha_B$ constant or slowly varying, the stability condition amounts to $\alpha_B < 0$, so $\alpha_M > 0$; 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 an MI observable quantity, whether this means that perturbations grow faster than in $\Lambda$CDM or not is a model-dependent statement.

* From Equation (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

$$ \alpha_1\alpha_B + \alpha_1^2 = 0. \qquad (87) $$

We then 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

$$ G_{3X} = 0, \qquad (88) $$

which actually means that the $G_3$ sector, after an integration by parts, can be absorbed in $K(\phi, X)$. Therefore, 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

$$ \Psi(r) - \Phi(r) = 2h_1 \frac{G_N M}{r}. \qquad (89) $$

This means that radiation, being conformally invariant (the 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

$$ \alpha_t = \frac{\alpha_M^2}{c_s^2(2\alpha_K + 3\alpha_M^2)}, \qquad (90) $$

which becomes

$$ \alpha_t = \frac{1}{3c_s^2} \qquad (91) $$
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$, such that there is a Yukawa correction, but $\eta = 1$ at all quasi-static scales.

## 7. Local Tests of Gravity

Gravity has been tested for a long time in the laboratory and within the solar system (see, e.g., [37]). The generic outcome of these experiments is that Einsteinian gravity works well at all of 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,

$$ \Psi(r) = -\frac{G_N M}{r} (1 + \alpha e^{-mr}). \qquad (92) $$

(Here we drop the subscript from $\alpha_t$ since we need to 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 the laboratory, one has [38] $|\alpha| \le 10^6$ at $\lambda \sim 10^{-5}$ m and $|\alpha| \le 10^{-2}$ at $\lambda \sim 10^{-3}$ m (Casimir-force experiments probe even shorter scales, but the constraints on $|\alpha|$ become correspondingly weaker). At planetary scales, one has $|\alpha| \le 10^{-6}$ for $\lambda \sim 10^6$ m (Earth–Moon distance), and $|\alpha| \le 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., [37]

$$ U = \int \frac{\rho(\mathbf{x}, t)}{|\mathbf{x} - \mathbf{x}'|} d^3 \mathbf{x}', \qquad (93) $$

and $V_i$ as a velocity-weighted potential

$$ V_i = \int \frac{\rho(\mathbf{x}, t) v_i}{|\mathbf{x} - \mathbf{x}'|} d^3 \mathbf{x}'. \qquad (94) $$

We can then write the parametrized post-Newtonian metric as follows:

$$ g_{00} = -1 + 2U - 2(1 + \beta)U^2 \qquad (95) $$

$$ g_{0i} = -\frac{1}{2}(3 + 4(1 + \gamma))V_i \qquad (96) $$

$$ g_{ij} = (1 + 2(1 + \gamma)U)\delta_{ij}. \qquad (97) $$

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

$$ \beta = \frac{\beta_0 \alpha}{2 (1+\alpha)^2} \qquad (98) $$

$$ \gamma = -2 \frac{\alpha}{(1+\alpha)} \qquad (99) $$

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 a linear level.
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$, such that there is a Yukawa correction, but $\eta = 1$ at all quasi-static scales.

## 7. Local Tests of Gravity

Gravity has been tested for a long time in the laboratory and within the solar system (see, e.g., [37]). The generic outcome of these experiments is that Einsteinian gravity works well at all of 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,

$$ \Psi(r) = -\frac{G_N M}{r} (1 + \alpha e^{-mr}). \qquad (92) $$

(Here we drop the subscript from $\alpha_t$ since we need to 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 the laboratory, one has [38] $|\alpha| \le 10^6$ at $\lambda \sim 10^{-5}$ m and $|\alpha| \le 10^{-2}$ at $\lambda \sim 10^{-3}$ m (Casimir-force experiments probe even shorter scales, but the constraints on $|\alpha|$ become correspondingly weaker). At planetary scales, one has $|\alpha| \le 10^{-6}$ for $\lambda \sim 10^6$ m (Earth–Moon distance), and $|\alpha| \le 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., [37]

$$ U = \int \frac{\rho(\mathbf{x}, t)}{|\mathbf{x} - \mathbf{x}'|} d^3 \mathbf{x}', \qquad (93) $$

and $V_i$ as a velocity-weighted potential

$$ V_i = \int \frac{\rho(\mathbf{x}, t) v_i}{|\mathbf{x} - \mathbf{x}'|} d^3 \mathbf{x}'. \qquad (94) $$

We can then write the parametrized post-Newtonian metric as follows:

$$ g_{00} = -1 + 2U - 2(1 + \beta)U^2 \qquad (95) $$

$$ g_{0i} = -\frac{1}{2}(3 + 4(1 + \gamma))V_i \qquad (96) $$

$$ g_{ij} = (1 + 2(1 + \gamma)U)\delta_{ij}. \qquad (97) $$

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

$$ \beta = \frac{\beta_0 \alpha}{2 (1+\alpha)^2} \qquad (98) $$

$$ \gamma = -2 \frac{\alpha}{(1+\alpha)} \qquad (99) $$

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 a linear level.
Local tests of gravity can, therefore, measure the Yukawa correction for both $Φ, Ψ$, i.e., $α_t, α_s$ and $λ$, and the ratio $Φ/Ψ$, in a model-independent way. The parameter $|γ|$, for instance, is constrained to be less than $10^{-5}$ [39], inducing a similar constraint on $α$ at large scales. A similar constraint applies also to $β$. 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 have 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 $φ$. This is the so-called chameleon effect [40,41]. 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 [42,43]. Finally, a third mechanism appears if the coupling $α$ sets on a vanishing value in structures (high density regions), via a symmetry restoration, while being different from zero at the background (low density) [44–47]. 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 that the baryons are completely decoupled from the scalar field.

In 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 both the baryon-to-photon constant ratio $η_b$ and 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 $η_b$ by CMB measurements, one can find constraints on $G_N(t_{BBN})/G_N(t_0) \approx 1 \pm 0.2$ [48] 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, which is that $G_N$ did not vary at all or at least in any way less than 0.2 throughout the expansion, implies $|α_M(t_0)| < 0.2(H_0T_0)^{-1}$ (equal to $\approx 0.2$ in Λ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 an estimate of $Ω_b h^2$ from CMB that assumes ΛCDM. Additionally, $G_N$ is, in fact, degenerate with the number $g_*$ of relativistic degrees of freedom at nucleosynthesis, such that the bound applies to $G_N g_*$ rather than to $G_N$ alone. Finally, a simultaneous change in the other standard-model parameters might considerably weaken the constraint (see [49,50]).

## 8. 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 [25,51],

$$c_T^2 = 1 + \alpha_T. \qquad (100)$$

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 [3]. Since the signals arrived within a 1s difference and since the light took $10^{15}$ s to reach us, we have that $|c_T^2/c^2 - 1| < 10^{-15}$. Such a 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
Local tests of gravity can, therefore, measure the Yukawa correction for both $Φ, Ψ$, i.e., $α_t, α_s$ and $λ$, and the ratio $Φ/Ψ$, in a model-independent way. The parameter $|γ|$, for instance, is constrained to be less than $10^{-5}$ [39], inducing a similar constraint on $α$ at large scales. A similar constraint applies also to $β$. 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 have 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 $φ$. This is the so-called chameleon effect [40,41]. 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 [42,43]. Finally, a third mechanism appears if the coupling $α$ sets on a vanishing value in structures (high density regions), via a symmetry restoration, while being different from zero at the background (low density) [44–47]. 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 that the baryons are completely decoupled from the scalar field.

In 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 both the baryon-to-photon constant ratio $η_b$ and 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 $η_b$ by CMB measurements, one can find constraints on $G_N(t_{BBN})/G_N(t_0) \approx 1 \pm 0.2$ [48] 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, which is that $G_N$ did not vary at all or at least in any way less than 0.2 throughout the expansion, implies $|α_M(t_0)| < 0.2(H_0T_0)^{-1}$ (equal to $\approx 0.2$ in Λ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 an estimate of $Ω_b h^2$ from CMB that assumes ΛCDM. Additionally, $G_N$ is, in fact, degenerate with the number $g_*$ of relativistic degrees of freedom at nucleosynthesis, such that the bound applies to $G_N g_*$ rather than to $G_N$ alone. Finally, a simultaneous change in the other standard-model parameters might considerably weaken the constraint (see [49,50]).

## 8. 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 [25,51],

$$c_T^2 = 1 + \alpha_T. \qquad (100)$$

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 [3]. Since the signals arrived within a 1s difference and since the light took $10^{15}$ s to reach us, we have that $|c_T^2/c^2 - 1| < 10^{-15}$. Such a 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) [35,52–56]. That is, we need to have $G_5 = 0$ and $G_{4,X} = 0$. In other words, the surviving Lagrangian has an arbitrary $K$, $G_3$ but a vanishing $G_5$ and an 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 [57]. 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 [58].

The constraints on $c_T$ also affect directly $\eta$. From Equation (65), one has in fact

$$h_2 = \frac{1}{c_T^2}. \qquad (101)$$

Hence, the GW constraint $h_2 = 1$ implies that $\eta$ should also be equal to unity for sufficiently large scales (small $k$) [59]; i.e., it should recover its general relativity value. The obvious exception are theories without a mass scale in addition to the Planck mass [60], in which case $\eta = h_4/h_5$ at all scales. On the other hand, no obvious GW constraint affects $\eta$.

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 [33], the GW amplitude $h$ obeys the equation

$$\ddot{h} + (3 + \alpha_M)H\dot{h} + c_T^2 \frac{k^2}{a^2} h = 0. \qquad (102)$$

Assuming $c_T = 1$, this equation in the sub-horizon limit is solved by [61]

$$h_a = \left( \frac{M_{*,\text{em}}}{M_{*,\text{obs}}} \right) \times h_s, \qquad (103)$$

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., [62], Equation (4.189))

$$h_s = \frac{4}{d_L} \left( \frac{GM_c}{c^2} \right)^{5/3} \left( \frac{\pi f_{\text{GW}}}{c} \right)^{2/3}. \qquad (104)$$

Here, $d_L$ is the luminosity distance, $M_c$ the so-called chirp mass, and $f_{\text{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 [61,63,64]

$$d_{\text{GW}} = \left( \frac{M_{*,\text{obs}}}{M_{*,\text{em}}} \right) d_L. \qquad (105)$$

Comparing this with an optical determination of $d_L$ leads to a direct measurement of $M_*$ 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_*$ 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$.
(in cosmological scales) [35,52–56]. That is, we need to have $G_5 = 0$ and $G_{4,X} = 0$. In other words, the surviving Lagrangian has an arbitrary $K$, $G_3$ but a vanishing $G_5$ and an 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 [57]. 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 [58].

The constraints on $c_T$ also affect directly $\eta$. From Equation (65), one has in fact

$$h_2 = \frac{1}{c_T^2}. \qquad (101)$$

Hence, the GW constraint $h_2 = 1$ implies that $\eta$ should also be equal to unity for sufficiently large scales (small $k$) [59]; i.e., it should recover its general relativity value. The obvious exception are theories without a mass scale in addition to the Planck mass [60], 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 [33], the GW amplitude $h$ obeys the equation

$$\ddot{h} + (3 + \alpha_M)H\dot{h} + c_T^2 \frac{k^2}{a^2}h = 0. \qquad (102)$$

Assuming $c_T = 1$, this equation in the sub-horizon limit is solved by [61]

$$h_a = \left( \frac{M_{*,\text{em}}}{M_{*,\text{obs}}} \right) \times h_s, \qquad (103)$$

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., [62], Equation (4.189))

$$h_s = \frac{4}{d_L} \left( \frac{GM_c}{c^2} \right)^{5/3} \left( \frac{\pi f_{\text{GW}}}{c} \right)^{2/3}. \qquad (104)$$

Here, $d_L$ is the luminosity distance, $M_c$ the so-called chirp mass, and $f_{\text{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 [61,63,64]

$$d_{\text{GW}} = \left( \frac{M_{*,\text{obs}}}{M_{*,\text{em}}} \right) d_L. \qquad (105)$$

Comparing this with an optical determination of $d_L$ leads to a direct measurement of $M_*$ 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_*$ 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$.
## 9. Model Dependence

The standard model of cosmology, Λ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, i.e., power-law spectra $k^{h_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) as well as quasi-relativistic particles, 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 BBN [65]. 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$, ΛCDM fits all the current cosmological data: the cosmic microwave background (CMB), the weak lensing data, the redshift distortion data (RSD), and the distance indicators (supernovae Ia, SNIa; baryon acoustic oscillations, BAOs; cosmic chronometers, CCH; gravitational waves, GW).

There are indeed a few discrepancies. In particular, two 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 [66], independent of cosmology, which deviates from the Planck [67] value obtained through extrapolation from the last scattering epoch performed assuming ΛCDM, $H_0 = (67.51 \pm 0.64)$ km/s/Mpc [67]. 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$) [67] differs from the late-universe value delivered by weak lensing, $\sigma_8 = 0.745 \pm 0.039$ [68] and by RSD data [69], $\sigma_8 = 0.75 \pm 0.024$.

Another source of discrepancies is related to the dark matter clustering [70]. 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 a 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 (Λ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 ([67], Figure 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, BAOs, or otherwise, basically measure the comoving distance

$$r(z) = \frac{1}{H_0 \sqrt{\Omega_{k0}}} \sinh \left( \sqrt{\Omega_{k0}} \int_0^z \frac{d\hat{z}}{E(\hat{z})} \right) \qquad (106)$$

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)$, i.e., $H_0$, cannot be measured (except for “standard sirens”,
## 9. Model Dependence

The standard model of cosmology, Λ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, i.e., power-law spectra $k^{h_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) as well as quasi-relativistic particles, 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 BBN [65]. 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$, ΛCDM fits all the current cosmological data: the cosmic microwave background (CMB), the weak lensing data, the redshift distortion data (RSD), and the distance indicators (supernovae Ia, SNIa; baryon acoustic oscillations, BAOs; cosmic chronometers, CCH; gravitational waves, GW).

There are indeed a few discrepancies. In particular, two 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 [66], independent of cosmology, which deviates from the Planck [67] value obtained through extrapolation from the last scattering epoch performed assuming ΛCDM, $H_0 = (67.51 \pm 0.64)$ km/s/Mpc [67]. 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$) [67] differs from the late-universe value delivered by weak lensing, $\sigma_8 = 0.745 \pm 0.039$ [68] and by RSD data [69], $\sigma_8 = 0.75 \pm 0.024$.

Another source of discrepancies is related to the dark matter clustering [70]. 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 a 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 (Λ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 ([67], Figure 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, BAOs, or otherwise, basically measure the comoving distance

$$r(z) = \frac{1}{H_0 \sqrt{\Omega_{k0}}} \sinh \left( \sqrt{\Omega_{k0}} \int_0^z \frac{d\hat{z}}{E(\hat{z})} \right) \quad (106)$$

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)$, i.e., $H_0$, cannot be measured (except for “standard sirens”,
i.e., gravitational waves [71]). Assuming for simplicity that $\Omega_{k0} = 0$, the only direct observable is $E(z) = H(z)/H_0$. If we also neglect radiation (a very good approximation for observations at redshift less than a few) and assume that, besides pressureless matter, we have a dark energy component with EOS $w(z)$, we have

$$E^2(z) = \Omega_{m0}(1+z)^3 + (1-\Omega_{m0})(1+z)^{3(1+\bar{w})} \quad (107)$$

where

$$\bar{w}(z) = \frac{1}{\ln(1+z)} \int_{0}^{z} \frac{w(\hat{z})}{1+\hat{z}} d\hat{z}. \quad (108)$$

We can then invert the relation (107) and obtain

$$w(z) = \frac{2(1+z)E_{E,z} - 3E^2}{3E^2 - 3\Omega_{m0}(1+z)^3} \quad (109)$$

where $E_{E,z}$ means differentiation with respect to redshift. It appears then that, in order to reconstruct $w(z)$, one needs to know $\Omega_{m0}$, in addition to $E(z)$. For instance, if the true cosmology were $\Lambda$CDM with $\Omega_{m0} = 0.3$, and we assumed erroneously that $\Omega_{m0} = 0.31$, we would infer $w(z = 0) = -0.986$ and $w(z = 1) = -0.897$, which is much 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 universally to gravity, unless one breaks the equivalence principle (see the “dark degeneracy” of [72]). Therefore, to measure $w(z)$, one has to assume a model, i.e., a parametrization, 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 simultaneously fix $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 complicated 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 is obvious, one cannot claim absolute model independence. The point is rather to clearly isolate the assumptions, and see how far can one reach with a minimum amount of them. In the following, we will assume the following:

(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.
i.e., gravitational waves [71]). Assuming for simplicity that $\Omega_{k0} = 0$, the only direct observable is $E(z) = H(z)/H_0$. If we also neglect radiation (a very good approximation for observations at redshift less than a few) and assume that, besides pressureless matter, we have a dark energy component with EOS $w(z)$, we have

$$E^2(z) = \Omega_{m0}(1+z)^3 + (1-\Omega_{m0})(1+z)^{3(1+\bar{w})} \quad (107)$$

where

$$\bar{w}(z) = \frac{1}{\ln(1+z)} \int_{0}^{z} \frac{w(\hat{z})}{1+\hat{z}} d\hat{z}. \quad (108)$$

We can then invert the relation (107) and obtain

$$w(z) = \frac{2(1+z)E_{E,z} - 3E^2}{3E^2 - 3\Omega_{m0}(1+z)^3} \quad (109)$$

where $E_{E,z}$ means differentiation with respect to redshift. It appears then that, in order to reconstruct $w(z)$, one needs to know $\Omega_{m0}$, in addition to $E(z)$. For instance, if the true cosmology were $\Lambda$CDM with $\Omega_{m0} = 0.3$, and we assumed erroneously that $\Omega_{m0} = 0.31$, we would infer $w(z = 0) = -0.986$ and $w(z = 1) = -0.897$, which is much 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 universally to gravity, unless one breaks the equivalence principle (see the “dark degeneracy” of [72]). Therefore, to measure $w(z)$, one has to assume a model, i.e., a parametrization, 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 simultaneously fix $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 complicated 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 is obvious, one cannot claim absolute model independence. The point is rather to clearly isolate the assumptions, and see how far can one reach with a minimum amount of them. In the following, we will assume the following:

(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. Therefore, we include in our treatment gravity plus at least one scalar, vector, and tensor field. 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.

## 10. Model-Independent Determination of the Homogeneous and Isotropic Geometry

What we observe in cosmology are redshifts and angular positions of sources. However, we need to build models for and test distances. Can we convert redshifts and angles into distances in an 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

$$ds^2 = -dt^2 + a^2(t) \left[ \frac{dr^2}{1 - kr^2} + r^2 d\theta^2 + r^2 \sin^2 \theta d\phi^2 \right] \quad (110)$$

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

$$ds^2 = -dt^2 + a^2(t) \left[ \frac{dr^2}{1 + \Omega_{k0}r^2} + r^2 d\theta^2 + r^2 \sin^2 \theta d\phi^2 \right]. \quad (111)$$

The value of $\Omega_{k0}$ has been estimated by Planck to be extremely small, $|\Omega_{k0}| < 0.004$ ([67], Table 5, last column), but, again, this is a model-dependent estimate, and for now we consider it as a free parameter. We see that, 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

$$t - t_0 = \int_1^a \frac{d\bar{a}}{\bar{a}E(\bar{a})} \quad (112)$$

where again $t$ is in units of $H_0^{-1}$.

BAOs are the remnant of 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

$$R = \frac{4}{3} \sqrt{\frac{\Omega_{\gamma 0}}{\Omega_{m 0} \Omega_{b 0}}} \ln \left( \frac{\sqrt{R_s(z_{drag}) + R_s(z_{eq})} + \sqrt{1 + R_s(z_{drag})}}{1 + \sqrt{R_s(z_{eq})}} \right) \quad (113)$$

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 [73]), and $z_{eq} = 2.396 \times 10^4 \Omega_{m 0} 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$; we need only 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

$$H_0 d_A(z) \equiv \frac{R}{\theta} = \frac{1}{(1+z)\sqrt{\Omega_{k0}}} \sinh(\sqrt{\Omega_{k0}} \int_0^z \frac{d\hat{z}}{E(\hat{z})}). \quad (114)$$
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. Therefore, we include in our treatment gravity plus at least one scalar, vector, and tensor field. 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.

## 10. Model-Independent Determination of the Homogeneous and Isotropic Geometry

What we observe in cosmology are redshifts and angular positions of sources. However, we need to build models for and test distances. Can we convert redshifts and angles into distances in an 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

$$ds^2 = -dt^2 + a^2(t) \left[ \frac{dr^2}{1 - kr^2} + r^2 d\theta^2 + r^2 \sin^2 \theta d\phi^2 \right] \quad (110)$$

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

$$ds^2 = -dt^2 + a^2(t) \left[ \frac{dr^2}{1 + \Omega_{k0} r^2} + r^2 d\theta^2 + r^2 \sin^2 \theta d\phi^2 \right]. \quad (111)$$

The value of $\Omega_{k0}$ has been estimated by Planck to be extremely small, $|\Omega_{k0}| < 0.004$ ([67], Table 5, last column), but, again, this is a model-dependent estimate, and for now we consider it as a free parameter. We see that, 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

$$t - t_0 = \int_1^a \frac{d\bar{a}}{\bar{a}E(\bar{a})} \quad (112)$$

where again $t$ is in units of $H_0^{-1}$.

BAOs are the remnant of 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 ΛCDM, $R$ (in units of $H_0^{-1}$) is equal to

$$R = \frac{4}{3} \sqrt{\frac{\Omega_{\gamma 0}}{\Omega_{m 0} \Omega_{b 0}}} \ln \left( \frac{\sqrt{R_s(z_{drag}) + R_s(z_{eq})} + \sqrt{1 + R_s(z_{drag})}}{1 + \sqrt{R_s(z_{eq})}} \right) \quad (113)$$

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 [73]), and $z_{eq} = 2.396 \times 10^4 \Omega_{m 0} 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$; we need only 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

$$H_0 d_A(z) \equiv \frac{R}{\theta} = \frac{1}{(1+z)\sqrt{\Omega_{k0}}} \sinh(\sqrt{\Omega_{k0}} \int_0^z \frac{d\hat{z}}{E(\hat{z})}). \quad (114)$$
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

$$E(z) = \frac{dz}{R}. \qquad (115)$$

This means that BAOs can estimate at every redshift two combinations involving $E(z)$ and $\Omega_{k0}$ and therefore determine both in an MI way. Therefore, the FLRW metric can in principle be reconstructed within the range covered by BAO observations, without assumptions in addition to (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 [74]

$$r_{12}^2(z_1, z_2, \theta) = r^2 \left( \frac{z_1 + z_2}{2} \right) \sin^2 \frac{\theta}{2} + r^2 \left( \frac{z_2 - z_1}{2} \right) \cos^2 \frac{\theta}{2} \qquad (116)$$

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

We move then to the next layer: perturbations.

## 11. Measuring Gravity: The Anisotropic Stress

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

$$\eta = -\frac{\Phi}{\Psi}. \qquad (117)$$

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 $v$ and relativistic factor $\gamma^2 = (1 - v^2)^{-1}$ in a weak-field Minkowski metric is in fact [75]

$$\gamma^2 \dot{v} = \gamma^2 [2v(v \cdot \nabla)(\Psi - \Phi) - v^2 \nabla(\Psi - \Phi)] - \nabla \Psi. \qquad (118)$$

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 the clustering and velocity of galaxies.

The linear gravitational perturbation theory provides 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 also define the growth function

$$G(k, z) = \frac{\delta_m(k, z)}{\delta_m(k, 0)} \qquad (119)$$

and the growth rate

$$f(k, z) = \frac{\delta'_{m}(k, z)}{\delta_{m}(k, z)} \qquad (120)$$

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_{\text{gal}}(k, z)$. Since galaxies are expected to be a biased tracer of mass, we need to introduce a bias function

$$b(k, z) = \frac{\delta_{\text{gal}}(k, z)}{\delta_m(k, z)} \qquad (121)$$
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

$$E(z) = \frac{dz}{R}. \qquad (115)$$

This means that BAOs can estimate at every redshift two combinations involving $E(z)$ and $\Omega_{k0}$ and therefore determine both in an MI way. Therefore, the FLRW metric can in principle be reconstructed within the range covered by BAO observations, without assumptions in addition to (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 [74]

$$r_{12}^2(z_1, z_2, \theta) = r^2 \left( \frac{z_1 + z_2}{2} \right) \sin^2 \frac{\theta}{2} + r^2 \left( \frac{z_2 - z_1}{2} \right) \cos^2 \frac{\theta}{2} \qquad (116)$$

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

We move then to the next layer: perturbations.

## 11. Measuring Gravity: The Anisotropic Stress

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

$$\eta = -\frac{\Phi}{\Psi}. \qquad (117)$$

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 $v$ and relativistic factor $\gamma^2 = (1 - v^2)^{-1}$ in a weak-field Minkowski metric is in fact [75]

$$\gamma^2 \dot{v} = \gamma^2 [2v(v \cdot \nabla)(\Psi - \Phi) - v^2 \nabla(\Psi - \Phi)] - \nabla \Psi. \qquad (118)$$

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 the clustering and velocity of galaxies.

The linear gravitational perturbation theory provides 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 also define the growth function

$$G(k, z) = \frac{\delta_m(k, z)}{\delta_m(k, 0)} \qquad (119)$$

and the growth rate

$$f(k, z) = \frac{\delta'_{m}(k, z)}{\delta_{m}(k, z)} \qquad (120)$$

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_{\text{gal}}(k, z)$. Since galaxies are expected to be a biased tracer of mass, we need to introduce a bias function

$$b(k, z) = \frac{\delta_{\text{gal}}(k, z)}{\delta_{m}(k, z)} \qquad (121)$$
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} = \text{const} \times \rho_{m}$. If b > 1 (< 1), then galaxies are more (or less) clustered than matter. Moreover, since we observe in redshift space, which means we observe a sum of cosmic expansion and the radial component of the local peculiar velocity, to convert to real space, we need the Kaiser transformation [76], 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 [77]

$$P_{gal}(k, z, \mu) = (A + R\mu^{2})^{2} \qquad (122)$$

where

$$A(k, z) = Gb\delta_{m0}, \quad R(k, z) = Gf\delta_{m0} \qquad (123)$$

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

$$\frac{1}{2\pi^2} \int \delta_{m0}^2 W_8^2(kR_8) d^3k = \sigma_8^2 \qquad (124)$$

where $W_8(kR_8)$ is the window function for an 8 h$^{-1}$Mpc sphere, and $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. Equation (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 Equation (52), we see that by estimating the shear distortion, one can measure the quantity

$$Y(1 + \eta)\rho_{m}\delta_{m} = Y(1 + \eta)\frac{\Omega_{m0}(1 + z)^{3}}{E^{2}}G\delta_{m0}. \qquad (125)$$

Since $E(z)$ can be estimated independently, we define another observable, to be denoted $L$ [25], as follows:

$$L(k, z) = \Omega_{m0}Y(1 + \eta)G\delta_{m0}. \qquad (126)$$

Together with $E = H(z)/H_{0}$, the quantities $A, R, L$ are the only cosmological information one can directly gather at the linear level (as we have seen, $\Omega_{k0}$ is also 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}$; however, this is equivalent to $2RE^{2}(2 + (\log ER)')/3(1 + z)^{3}$. At least at the linear level, one could 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 [78] shows that, employing higher-order statistics, we can obtain more MI information, but we will not consider this here.

We can now write the lensing equation in Fourier space in the following way (see [79]):

$$-\hat{k}^2(\Psi - \Phi) = \frac{3(1+z)^3L}{2E^2} \qquad (127)$$
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} = \text{const} \times \rho_m$. If b > 1 (< 1), then galaxies are more (or less) clustered than matter. Moreover, since we observe in redshift space, which means we observe a sum of cosmic expansion and the radial component of the local peculiar velocity, to convert to real space, we need the Kaiser transformation [76], 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 [77]

$$P_{gal}(k, z, \mu) = (A + R\mu^2)^2 \qquad (122)$$

where

$$A(k, z) = Gb\delta_{m0}, \quad R(k, z) = Gf\delta_{m0} \qquad (123)$$

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

$$\frac{1}{2\pi^2} \int \delta_{m0}^2 W_8^2(kR_8) d^3k = \sigma_8^2 \qquad (124)$$

where $W_8(kR_8)$ is the window function for an 8 h⁻¹Mpc sphere, and $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. Equation (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 Equation (52), we see that by estimating the shear distortion, one can measure the quantity

$$Y(1 + \eta)\rho_m\delta_m = Y(1 + \eta)\frac{\Omega_{m0}(1 + z)^3}{E^2}G\delta_{m0}. \qquad (125)$$

Since $E(z)$ can be estimated independently, we define another observable, to be denoted $L$ [25], as follows:

$$L(k, z) = \Omega_{m0}Y(1 + \eta)G\delta_{m0}. \qquad (126)$$

Together with $E = H(z)/H_0$, the quantities $A, R, L$ are the only cosmological information one can directly gather at the linear level (as we have seen, $\Omega_{k0}$ is also 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}$; however, this is equivalent to $2RE^2(2 + (\log ER)')/3(1 + z)^3$. At least at the linear level, one could 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 [78] shows that, employing higher-order statistics, we can obtain more MI information, but we will not consider this here.

We can now write the lensing equation in Fourier space in the following way (see [79]):

$$-\hat{k}^2(\Psi - \Phi) = \frac{3(1 + z)^3 L}{2E^2} \qquad (127)$$
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,

$$ \delta_m'' + \delta_m'(2 + \frac{E'}{E}) = -\hat{k}^2 \Psi, \qquad (128) $$

which 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

$$ -\hat{k}^2 \Psi = R' + R(2 + \frac{E'}{E}) . \qquad (129) $$

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 Equation (129), employing also Equation (21) as

$$ \frac{3}{2} \frac{\Omega_{m0} (1+z)^3 Y}{f E^2} = \frac{R'}{R} + (2 + \frac{E'}{E}). \qquad (130) $$

We see then that $Y$ is not, unfortunately, an 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_{m0}^{0.55}$ are known, and we have that $\Omega_{m0}$ is uniquely measured by a combination of $R, R', E, E'$. Any two measures at different $k$ or $z$ values must then give the same $\Omega_{m0}$. We show below that $\eta$, in contrast to $Y$, is an 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 complicated, possibly non-linear and hydrodynamical processes; thus, even if $b$ depends on gravity, we do not know how. Additionally, 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 such as $n_s$ or its “running”, nothing precludes the possibility that an unknown process, for instance the presence of early dark energy or early modified gravity, has distorted the spectrum at some intermediate redshift between the last scattering and today. Therefore, in order to obtain model-independent measures, we should eliminate both $b$ and $\delta_{m0}$. It was shown in [25,79] that one can obtain only three statistics where the effects of the shape of the primordial power spectrum is canceled out, namely

$$ P_1 \equiv \frac{R}{A} = \frac{f}{b'}, \qquad (131) $$

$$ P_2 \equiv \frac{L}{R} = \frac{\Omega_{m0} Y (1 + \eta)}{f}, \qquad (132) $$

$$ P_3 \equiv \frac{R'}{R} = f + \frac{f'}{f} = \frac{(f \sigma_8(k, z))'}{f \sigma_8(k, z)}. \qquad (133) $$

In the last equation, we introduced the often-employed quantity

$$ f \sigma_8(k, z) = \sigma_8 G(k, z) f(k, z) = R \frac{\sigma_8}{\delta_{m0}}. \qquad (134) $$

Notice that we are not defining $\sigma_8(k, z)$ as an integral over the power spectrum at $z$, as in Equation (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.
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,

$$ \delta_m'' + \delta_m'(2 + \frac{E'}{E}) = -\hat{k}^2 \Psi, \qquad (128) $$

which 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

$$ -\hat{k}^2 \Psi = R' + R(2 + \frac{E'}{E}) . \qquad (129) $$

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 Equation (129), employing also Equation (21) as

$$ \frac{3}{2} \frac{\Omega_{m0}(1+z)^3 Y}{f E^2} = \frac{R'}{R} + (2 + \frac{E'}{E}). \qquad (130) $$

We see then that $Y$ is not, unfortunately, an 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_{m0}^{0.55}$ are known, and we have that $\Omega_{m0}$ is uniquely measured by a combination of $R, R', E, E'$. Any two measures at different $k$ or $z$ values must then give the same $\Omega_{m0}$. We show below that $\eta$, in contrast to $Y$, is an 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 complicated, possibly non-linear and hydrodynamical processes; thus, even if $b$ depends on gravity, we do not know how. Additionally, 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 such as $n_s$ or its “running”, nothing precludes the possibility that an unknown process, for instance the presence of early dark energy or early modified gravity, has distorted the spectrum at some intermediate redshift between the last scattering and today. Therefore, in order to obtain model-independent measures, we should eliminate both $b$ and $\delta_{m0}$. It was shown in [25,79] that one can obtain only three statistics where the effects of the shape of the primordial power spectrum is canceled out, namely

$$ P_1 \equiv \frac{R}{A} = \frac{f}{b'}, \qquad (131) $$

$$ P_2 \equiv \frac{L}{R} = \frac{\Omega_{m0} Y(1 + \eta)}{f}, \qquad (132) $$

$$ P_3 \equiv \frac{R'}{R} = f + \frac{f'}{f} = \frac{(f\sigma_8(k, z))'}{f\sigma_8(k, z)}. \qquad (133) $$

In the last equation, we introduced the often-employed quantity

$$ f\sigma_8(k, z) = \sigma_8 G(k, z) f(k, z) = R \frac{\sigma_8}{\delta_{m0}}. \qquad (134) $$

Notice that we are not defining $\sigma_8(k, z)$ as an integral over the power spectrum at $z$, as in Equation (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), such 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, to predict the theoretical value of $R$ (or $f$) as a function of the gravity parameter $\Upsilon$ from Equation (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 [80], 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 an MI measurement of gravity. Exactly as we have seen for the dark energy EOS $w(z)$, if one parametrizes $\Upsilon(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 $\Upsilon(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 [81] and references therein)

$$E_g = \left\langle \frac{a\nabla^2(\Psi - \Phi)}{3H_0^2 f \delta_m} \right\rangle_k . \qquad (135)$$

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

$$P_2 = 2E_g. \qquad (136)$$

The $E_g$ statistics has been used several times as a test of modified gravity [81–84]. 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 $\Upsilon$. 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}$, to different initial conditions, or to non-standard modified gravity parameters $\Upsilon, \eta$. As previously shown, 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 MI statistics to measure gravity, in particular the parameter $\eta$, by combining the equation for the growth of structure formation, Equation (129), with the lensing Equation (127) and with Equations (14) and (21). We see then that the gravitational slip as a function of model-independent observables is given by

$$\eta_{\text{obs}} \equiv \frac{3P_2(1+z)^3}{2E^2\left(P_3+2+\frac{E'}{E}\right)} - 1 = \eta. \qquad (137)$$

In order to distinguish the observables from the theoretical expectations, we denoted the combination on the rhs of this equation as $\eta_{\text{obs}}$. The statistics $\eta_{\text{obs}}$ is model-independent because it estimates $\eta$ directly 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$. Therefore, if observationally one finds $\eta_{\text{obs}} \neq 1$, then $\Lambda$CDM and all of the models in standard gravity and in which dark energy is a perfect fluid are ruled out. As a consequence, cautionary remarks such as those in [85], 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_{\text{obs}}$. The price to pay is that Equation (137) depends on derivatives of $E$ and, through $P_3$, of $f\sigma_8(z)$.
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), such 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, to predict the theoretical value of $R$ (or $f$) as a function of the gravity parameter $\Upsilon$ from Equation (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 [80], 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 an MI measurement of gravity. Exactly as we have seen for the dark energy EOS $w(z)$, if one parametrizes $\Upsilon(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 $\Upsilon(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 [81] and references therein)

$$E_g = \left\langle \frac{a\nabla^2(\Psi - \Phi)}{3H_0^2 f \delta_m} \right\rangle_k . \qquad (135)$$

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

$$P_2 = 2E_g. \qquad (136)$$

The $E_g$ statistics has been used several times as a test of modified gravity [81–84]. 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 $\Upsilon$. 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}$, to different initial conditions, or to non-standard modified gravity parameters $\Upsilon, \eta$. As previously shown, 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 MI statistics to measure gravity, in particular the parameter $\eta$, by combining the equation for the growth of structure formation, Equation (129), with the lensing Equation (127) and with Equations (14) and (21). We see then that the gravitational slip as a function of model-independent observables is given by

$$\eta_{\text{obs}} \equiv \frac{3P_2(1+z)^3}{2E^2 \left(P_3 + 2 + \frac{E'}{E}\right)} - 1 = \eta. \qquad (137)$$

In order to distinguish the observables from the theoretical expectations, we denoted the combination on the rhs of this equation as $\eta_{\text{obs}}$. The statistics $\eta_{\text{obs}}$ is model-independent because it estimates $\eta$ directly 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$. Therefore, if observationally one finds $\eta_{\text{obs}} \neq 1$, then $\Lambda$CDM and all of the models in standard gravity and in which dark energy is a perfect fluid are ruled out. As a consequence, cautionary remarks such as those in [85], 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_{\text{obs}}$. The price to pay is that Equation (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., [78]. One interesting case is provided by relaxed galaxy clusters, for which we can reasonably expect that the virial theorem is at least approximately respected. In this case, we can directly measure the potential Ψ 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 [30]. 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_{\text{obs}}(P_2, P_3, E)$ is, therefore, a model-independent test of gravity: it does not depend on the bias, on initial conditions, nor on other unobservable quantities such as $\Omega_{m0}$ and $\sigma_8$. If $\eta_{\text{obs}} \neq 1$, gravity is not Einsteinian; if $\eta_{\text{obs}}$ does not have the same $k^2$ dependence as the Horndeski's theory, the entire Horndeski model is rejected. All this, of course, assumes that our conditions (a)–(d) are verified.

## 12. 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 Equation (128) is modified since the continuity and Euler equations, which come directly from the conservation of the energy-momentum tensor, now read

$$ \delta' = -\frac{1+w}{aH}(\theta - 3aH\Phi') - 3aH(c_s^2 - w)\delta \quad (138) $$

$$ \theta' = -(1-3w)\theta - \frac{w'}{1+w}\theta + \frac{c_s^2}{1+w}aH\hat{k}^2\delta - aH\hat{k}^2\sigma + aH\hat{k}^2\Psi \quad (139) $$

where the sound speed is $c_s^2 \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 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:

$$ \frac{\delta''}{1+w} + \left(2+\frac{H'}{H}+3(c_s^2-2w)\right) \frac{\delta'}{1+w} + 6(c_s^2-w) \left(1+\frac{H'}{H}\right) \frac{\delta}{1+w} = -\hat{k}^2\Psi, \quad (140) $$

which reduces to Equation (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)}$, such that the lensing Equation (127) would now read

$$ -\hat{k}^2(\Psi - \Phi) = \frac{3(1+z)^{3(1+w)}L}{2E^2}. \quad (141) $$

Taking the appropriate ratios of the two equations above, we can obtain $\eta$ as we did for Equation (137), but this time, some extra term appear:
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., [78]. One interesting case is provided by relaxed galaxy clusters, for which we can reasonably 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 [30]. 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_{\text{obs}}(P_2, P_3, E)$ is, therefore, a model-independent test of gravity: it does not depend on the bias, on initial conditions, nor on other unobservable quantities such as $\Omega_{m0}$ and $\sigma_8$. If $\eta_{\text{obs}} \neq 1$, gravity is not Einsteinian; if $\eta_{\text{obs}}$ does not have the same $k^2$ dependence as the Horndeski’s theory, the entire Horndeski model is rejected. All this, of course, assumes that our conditions (a)–(d) are verified.

## 12. 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 Equation (128) is modified since the continuity and Euler equations, which come directly from the conservation of the energy-momentum tensor, now read

$$ \delta' = -\frac{1+w}{aH}(\theta - 3aH\Phi') - 3aH(c_s^2 - w)\delta \quad (138) $$

$$ \theta' = -(1-3w)\theta - \frac{w'}{1+w}\theta + \frac{c_s^2}{1+w}aH\hat{k}^2\delta - aH\hat{k}^2\sigma + aH\hat{k}^2\Psi \quad (139) $$

where the sound speed is $c_s^2 \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 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:

$$ \frac{\delta''}{1+w} + \left(2 + \frac{H'}{H} + 3(c_s^2 - 2w)\right) \frac{\delta'}{1+w} + 6(c_s^2 - w) \left(1 + \frac{H'}{H}\right) \frac{\delta}{1+w} = -\hat{k}^2\Psi, \quad (140) $$

which reduces to Equation (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)}$, such that the lensing Equation (127) would now read

$$ -\hat{k}^2(\Psi - \Phi) = \frac{3(1+z)^{3(1+w)}L}{2E^2}. \quad (141) $$

Taking the appropriate ratios of the two equations above, we can obtain $\eta$ as we did for Equation (137), but this time, some extra term appear:
$$ \frac{3(1+w)P_2(1+z)^{3(1+w)}}{2E^2 \left( P_3 + 2 + \frac{E'}{E} + W_1 + \frac{W_2}{f} \left( 1 + \frac{E'}{E} \right) \right)} - 1 = \eta_{\text{obs}} \qquad (142) $$

where $W_1 = 3(c_s^2 - 2w)$ and $W_2 = 6(c_s^2 - w)$. Both $W_1$ and $W_2$ reduce to zero for standard cold dark matter, such that we recover Equation (137) exactly in that case. For a barotropic fluid such that $p = p(\rho)$, $c_s^2 = w$ and $W_2 = 0$. In this case, we have again an MI estimator for $\eta$, provided we know $c_s, w$. On the other hand, if $W_2 \neq 0$, we see that this estimation of $\eta$ contains the growth rate $f$, which we argued is not a model-independent observable in the linear regime. However, an extension of this formalism to the quasilinear scales [78] has shown that $f$ can indeed be recovered in a model-independent way, using observations of the bispectrum.

## 13. 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 Equation (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, such that $\alpha_T = 0$, and assuming $H$ is also measured in an MI way, we see that what can still be measured at the linear perturbation level are the combinations

$$ O_1 = \frac{\alpha_1 + \alpha_2}{\mu^2} \qquad (143) $$

$$ O_2 = \frac{\alpha_1 \alpha_M}{\mu^2}, \qquad (144) $$

which correspond 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$.

## 14. Data

In the next sections, we obtain an estimate of $\eta_{\text{obs}}$ using all the currently available data (This section and the next two are a summary of a published paper by A. M. Pinho, S. Casas, and L. Amendola, entitled 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'(z)$), $P_2(z)$, and $P_3(z)$ using all the currently available relevant data, shown in Figure 1, where we also plot the $\Lambda$CDM the curves of the different functions using the cosmological parameters from the TT+lowP+lensing Planck 2015 best-fits [67]. A similar analysis, with a much smaller dataset than is available, was carried out also in [86].

For the Hubble parameter measurements, we have used the most recent compilation of $H(z)$ data from [87], including the measurements from [88–91], the Baryon Oscillation Spectroscopic Survey (BOSS) [92–94], and the Sloan Digital Sky Survey (SDSS) [95,96]. In this compilation, the majority of the measurements were 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 [92] and [93] are obtained using the BAO signal in the Lyman-$\alpha$ forest distribution alone or cross-correlated with quasi-stellar objects (QSOs) (for the details of the method, we refer the reader to the original papers). The reference [96] provides the covariance matrix of three $H(z)$ measurements from the radial BAO galaxy distribution. To this compilation, we add the results from WiggleZ [97]. In addition to these, we use the recent results from [98] where a compilation of type Ia supernovae from CANDELS and CLASH Multi-cycle Treasury programs were analyzed, yielding a few tight measurements of the expansion rate $E(z)$.
$$ \frac{3(1+w)P_2(1+z)^{3(1+w)}}{2E^2 \left( P_3 + 2 + \frac{E'}{E} + W_1 + \frac{W_2}{f} \left( 1 + \frac{E'}{E} \right) \right)} - 1 = \eta_{\text{obs}} \quad (142) $$

where $W_1 = 3(c_s^2 - 2w)$ and $W_2 = 6(c_s^2 - w)$. Both $W_1$ and $W_2$ reduce to zero for standard cold dark matter, such that we recover Equation (137) exactly in that case. For a barotropic fluid such that $p = p(\rho)$, $c_s^2 = w$ and $W_2 = 0$. In this case, we have again an MI estimator for $\eta$, provided we know $c_s, w$. On the other hand, if $W_2 \neq 0$, we see that this estimation of $\eta$ contains the growth rate $f$, which we argued is not a model-independent observable in the linear regime. However, an extension of this formalism to the quasilinear scales [78] has shown that $f$ can indeed be recovered in a model-independent way, using observations of the bispectrum.

## 13. 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 Equation (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, such that $\alpha_T = 0$, and assuming $H$ is also measured in an MI way, we see that what can still be measured at the linear perturbation level are the combinations

$$ O_1 = \frac{\alpha_1 + \alpha_2}{\mu^2} \quad (143) $$

$$ O_2 = \frac{\alpha_1 \alpha_M}{\mu^2}, \quad (144) $$

which correspond 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$.

## 14. Data

In the next sections, we obtain an estimate of $\eta_{\text{obs}}$ using all the currently available data (This section and the next two are a summary of a published paper by A. M. Pinho, S. Casas, and L. Amendola, entitled 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'(z)$), $P_2(z)$, and $P_3(z)$ using all the currently available relevant data, shown in Figure 1, where we also plot the $\Lambda$CDM the curves of the different functions using the cosmological parameters from the TT+lowP+lensing Planck 2015 best-fits [67]. A similar analysis, with a much smaller dataset than is available, was carried out also in [86].

For the Hubble parameter measurements, we have used the most recent compilation of $H(z)$ data from [87], including the measurements from [88–91], the Baryon Oscillation Spectroscopic Survey (BOSS) [92–94], and the Sloan Digital Sky Survey (SDSS) [95,96]. In this compilation, the majority of the measurements were 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 [92] and [93] are obtained using the BAO signal in the Lyman-$\alpha$ forest distribution alone or cross-correlated with quasi-stellar objects (QSOs) (for the details of the method, we refer the reader to the original papers). The reference [96] provides the covariance matrix of three $H(z)$ measurements from the radial BAO galaxy distribution. To this compilation, we add the results from WiggleZ [97]. In addition to these, we use the recent results from [98] where a compilation of type Ia supernovae from CANDELS and CLASH Multi-cycle Treasury programs were analyzed, yielding a few tight measurements of the expansion rate $E(z)$.
The $E_g$ data include the results from KiDS+2dFLenS+GAMA [85], 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) [99], where the analysis combines the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS), the WiggleZ Dark Energy Survey, and the Baryon Oscillation Spectroscopic Survey (BOSS). The work of VIMOS Public Extragalactic Redshift Survey (VIPERS) [84] 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 < 60 \text{ h}^{-1}\text{Mpc}$ and in the linear regime, which is the one in which we are interested. 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 [81], 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 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 [100], the Subaru FMOS galaxy redshift survey (FastSound) [101], WiggleZ [97], the VIMOS-VLT Deep Survey (VVDS) [102], the VIMOS Public Extragalactic Redshift Survey (VIPERS) [84,103–105], and the Sloan Digital Sky Survey (SDSS) [96,106–112]. The values from [113,114] will not be considered since the $f\sigma_8(z)$ value is not directly reported.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Data sets used in this work (black dots with error bars), plotted with the corresponding theoretical ΛCDM prediction as a function of redshift (solid red line), using a Planck 2015 cosmology. **Left panel:** $E(z)$ data. We used the Planck 2015 value of $H_0$ to convert some of the data points from $H(z)$ to $E(z)$ (see main text). **Central panel:** Plot of the logarithm of the $f\sigma_8$ data points. **Right panel:** Data set for $P_2$, obtained using $E_g$ data and the relation (136) that converts between different notations. For $z > 0.5$, we see a larger discrepancy between ΛCDM and the data points, which was also noted in [85] and references therein.</figcaption></figure>

## 15. Reconstruction of Functions from Data

The only difficulty in obtaining $\eta_{\text{obs}}$ is that we need to take the ratios $P_2, P_3$ at the same redshift, while we have data points 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 obtains unavoidably different results, especially in the final errors. Here, we consider and compare three methods to obtain the value of $\eta_{\text{obs}}$: binning, the Gaussian process (GP), and generalized linear regression.
The $E_g$ data include the results from KiDS+2dFLenS+GAMA [85], 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) [99], where the analysis combines the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS), the WiggleZ Dark Energy Survey, and the Baryon Oscillation Spectroscopic Survey (BOSS). The work of VIMOS Public Extragalactic Redshift Survey (VIPERS) [84] 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 < 60 \text{ h}^{-1}\text{Mpc}$ and in the linear regime, which is the one in which we are interested. 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 [81], 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 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 [100], the Subaru FMOS galaxy redshift survey (FastSound) [101], WiggleZ [97], the VIMOS-VLT Deep Survey (VVDS) [102], the VIMOS Public Extragalactic Redshift Survey (VIPERS) [84,103–105], and the Sloan Digital Sky Survey (SDSS) [96,106–112]. The values from [113,114] will not be considered since the $f\sigma_8(z)$ value is not directly reported.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Data sets used in this work (black dots with error bars), plotted with the corresponding theoretical ΛCDM prediction as a function of redshift (solid red line), using a Planck 2015 cosmology. **Left panel:** $E(z)$ data. We used the Planck 2015 value of $H_0$ to convert some of the data points from $H(z)$ to $E(z)$ (see main text). **Central panel:** Plot of the logarithm of the $f\sigma_8$ data points. **Right panel:** Data set for $P_2$, obtained using $E_g$ data and the relation (136) that converts between different notations. For $z > 0.5$, we see a larger discrepancy between ΛCDM and the data points, which was also noted in [85] and references therein.</figcaption></figure>

## 15. Reconstruction of Functions from Data

The only difficulty in obtaining $\eta_{\text{obs}}$ is that we need to take the ratios $P_2, P_3$ at the same redshift, while we have data points 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 obtains unavoidably different results, especially in the final errors. Here, we consider and compare three methods to obtain the value of $\eta_{\text{obs}}$: binning, the Gaussian process (GP), and generalized linear regression.