Designing Horndeski and the effective fluid approach

TL;DR

This paper develops a simplified effective fluid approach for Horndeski theories, introduces a designer model matching $ ext{Λ}$CDM background, and uses it to analyze cosmological data, showing potential to alleviate tensions.

Contribution

It extends the effective fluid approach to Horndeski theories, implements it in a modified Boltzmann code, and introduces a designer model with well-behaved perturbations for cosmological analysis.

Findings

The effective fluid approach achieves 0.1% accuracy compared to hi_CLASS.

The kinetic braiding model's quasistatic approximation can break down.

The designer HDES model can alleviate the $ ext{2}\sigma$ tension in growth data.

Abstract

We present a family of designer Horndeski models, i.e. models that have a background exactly equal to that of the $\Lambda$CDM model but perturbations given by the Horndeski theory. Then, we extend the effective fluid approach to Horndeski theories, providing simple analytic formulae for the equivalent dark energy effective fluid pressure, density and velocity. We implement the dark energy effective fluid formulae in our code EFCLASS, a modified version of the widely used Boltzmann solver CLASS, and compare the solution of the perturbation equations with those of the code hi_CLASS which already includes Horndeski models. We find that our simple modifications to the vanilla code are accurate to the level of $\sim 0.1\%$ with respect to the more complicated hi_CLASS code. Furthermore, we study the kinetic braiding model both on and off the attractor and we find that even though the full case has a proper $\Lambda$CDM model limit for large $n$, it is not appropriately smooth, thus causing the quasistatic approximation to break down. Finally, we focus on our designer model (HDES), which has both a smooth $\Lambda$CDM limit and well-behaved perturbations, and we use it to perform Markov Chain Monte Carlo analyses to constrain its parameters with the latest cosmological data. We find that our HDES model can also alleviate the soft $2\sigma$ tension between the growth data and Planck 18 due to a degeneracy between $\sigma_8$ and one of its model parameters that indicates the deviation from the $\Lambda$CDM model.

Figures

Tables

Table 1. Table 1: The H ​ ( z ) 𝐻 𝑧 H(z) data used in the current analysis (in units of km ​ s − 1 ​ Mpc − 1 km superscript s 1 superscript Mpc 1 \textrm{km}~{}\textrm{s}^{-1}\textrm{Mpc}^{-1} ). This compilation is partly based on those of Refs.  Moresco et al. ( 2016 ) and Guo and Zhang ( 2016 ) .

$z$	$H\left(z\right)$	${\sigma }_H$	Ref.
$0.07$	$69.0$	$19.6$	Zhang et al. (2014)
$0.09$	$69.0$	$12.0$	Stern et al. (2010)
$0.12$	$68.6$	$26.2$	Zhang et al. (2014)
$0.17$	$83.0$	$8.0$	Stern et al. (2010)
$0.179$	$75.0$	$4.0$	Moresco et al. (2012)
$0.199$	$75.0$	$5.0$	Moresco et al. (2012)
$0.2$	$72.9$	$29.6$	Zhang et al. (2014)
$0.27$	$77.0$	$14.0$	Stern et al. (2010)
$0.28$	$88.8$	$36.6$	Zhang et al. (2014)
$0.35$	$82.7$	$8.4$	Chuang and Wang (2013)
$0.352$	$83.0$	$14.0$	Moresco et al. (2012)
$0.3802$	$83.0$	$13.5$	Moresco et al. (2016)
$0.4$	$95.0$	$17.0$	Stern et al. (2010)
$0.4004$	$77.0$	$10.2$	Moresco et al. (2016)
$0.4247$	$87.1$	$11.2$	Moresco et al. (2016)
$0.44$	$82.6$	$7.8$	Blake et al. (2012)
$0.44497$	$92.8$	$12.9$	Moresco et al. (2016)
$0.4783$	$80.9$	$9.0$	Moresco et al. (2016)

Table 2. Table 2: Compilation of the f ​ σ 8 ​ ( z ) 𝑓 subscript 𝜎 8 𝑧 f\sigma_{8}(z) measurements used in this analysis along with the reference matter density parameter Ω m 0 subscript Ω subscript 𝑚 0 \Omega_{m_{0}} (needed for the growth correction) and related references.

$z$	$f{\sigma }_8\left(z\right)$	${\sigma }_{f{\sigma }_8}\left(z\right)$	${\Omega }_{m,0}^{\text{ref}}$	Ref.
0.02	0.428	0.0465	0.3	Huterer et al. (2016)
0.02	0.398	0.065	0.3	Turnbull et al. (2012),Hudson and Turnbull (2013)
0.02	0.314	0.048	0.266	Davis et al. (2011),Hudson and Turnbull (2013)
0.10	0.370	0.130	0.3	Feix et al. (2015)
0.15	0.490	0.145	0.31	Howlett et al. (2015)
0.17	0.510	0.060	0.3	Song and Percival (2009)
0.18	0.360	0.090	0.27	Blake et al. (2013)
0.38	0.440	0.060	0.27	Blake et al. (2013)
0.25	0.3512	0.0583	0.25	Samushia et al. (2012)
0.37	0.4602	0.0378	0.25	Samushia et al. (2012)
0.32	0.384	0.095	0.274	Sanchez et al. (2014)
0.59	0.488	0.060	0.307115	Chuang et al. (2016)
0.44	0.413	0.080	0.27	Blake et al. (2012)
0.60	0.390	0.063	0.27	Blake et al. (2012)
0.73	0.437	0.072	0.27	Blake et al. (2012)
0.60	0.550	0.120	0.3	Pezzotta et al. (2016)
0.86	0.400	0.110	0.3	Pezzotta et al. (2016)
1.40	0.482	0.116	0.27	Okumura et al. (2016)
0.978	0.379	0.176	0.31	Zhao et al. (2018)
1.23	0.385	0.099	0.31	Zhao et al. (2018)
1.526	0.342	0.070	0.31	Zhao et al. (2018)
1.944	0.364	0.106	0.31	Zhao et al. (2018)

Table 3. Table 3: Λ Λ \Lambda CDM parameters with 68 % percent 68 68\% limits based on TT,TE,EE+lowP and a flat Λ Λ \Lambda CDM model (middle column) or a w 𝑤 w CDM model (right column); see Ref.  Aghanim et al. ( 2018 ) and the Planck chains archive.

Parameter	Value ($\Lambda$CDM)	Value ($w$CDM)
${\Omega }_b{h}^2$	$0.02225\pm 0.00016$	$0.02229\pm 0.00016$
${\Omega }_c{h}^2$	$0.1198\pm 0.0015$	$0.1196\pm 0.0015$
$n_s$	$0.9645\pm 0.0049$	$0.9649\pm 0.0048$
$H_0$	$67.27\pm 0.66$	$>81.3$
${\Omega }_m$	$0.3156\pm 0.0091$	${0.203}_{-0.065}^{+0.022}$
$w$	$-1$	$-{1.55}_{-0.38}^{+0.19}$
${\sigma }_8$	$0.831\pm 0.013$	${0.983}_{-0.055}^{+0.100}$

Table 4. Table 4: The best-fit parameters for the Λ Λ \Lambda CDM and the HDES ( n = 2 ) 𝑛 2 (n=2) models respectively.

Model	${\Omega }_{m,0}$	$100{\Omega }_b{h}^2$	$\widetilde{J_c}$	$h$	${\sigma }_8$
Best-fit values
$\Lambda$CDM	$0.311\pm 0.006$	$2.243\pm 0.014$	$0$	$0.680\pm 0.004$	$0.758\pm 0.025$
HDES	$0.313\pm 0.006$	$2.240\pm 0.014$	$-0.309\pm 0.244$	$0.678\pm 0.004$	$0.911\pm 0.068$

Table 5. Table 5: The χ 2 superscript 𝜒 2 \chi^{2} and AIC parameters for the Λ Λ \Lambda CDM and the HDES models respectively.

Model	${\chi }^2$	AIC	$\Delta$AIC
$\Lambda$CDM	$1087.64$	$1095.68$	$0$
HDES	$1086.30$	$1096.35$	$0.678$

Equations

$$S=\int d^{4}x\sqrt{-g}\left[\frac{1}{2\kappa }R+{\mathcal{L}}_m\right],$$

$$G_{\mu \nu }=\kappa T_{\mu \nu }^{(m)},$$

$$d{s}^2=-\left(1+2\Psi (\vec{x},t)\right)d{t}^2+a(t{)}^2(1+2\Phi (\vec{x},t))d{\vec{x}}^2,$$

$$T_{\nu }^{\mu }=P{\delta }_{\nu }^{\mu }+(\rho +P)U^{\mu }{U}_{\nu },$$

$$T_0^0$$

$$T_i^0$$

$$T_j^i$$

$$H^2$$

$$H^2+\dot{H}$$

$$-\frac{k^2}{a^2}\Phi +3\frac{\dot{a}}{a}\left(\frac{\dot{a}}{a}\Psi -\dot{\Phi }\right)=\frac{\kappa }{2}\delta T_0^0,$$

$$k^2\left(\frac{\dot{a}}{a}\Psi -\dot{\Phi }\right)=\frac{\kappa }{2}a(\bar{\rho }+\bar{P})\theta ,$$

$$-\frac{k^2}{3a^2}(\Phi +\Psi )+\left(2\frac{\ddot{a}}{a}+\frac{{\dot{a}}^2}{a^2}\right)\Psi +\frac{\dot{a}}{a}\left(\dot{\Psi }-3\dot{\Phi }\right)-\ddot{\Phi }=\frac{\kappa }{6}\delta T_i^i,$$

$$-k^2(\Phi +\Psi )=\frac{3\kappa }{2}{a}^2(\bar{\rho }+\bar{P})\sigma ,$$

$$\dot{\delta }=-(1+w)(\frac{\theta }{a}+3\dot{\Phi })-3\frac{\dot{a}}{a}\left(c_s^2-w\right)\delta ,$$

$$\dot{\theta }=-\frac{\dot{a}}{a}(1-3w)\theta -\frac{\dot{w}}{1+w}\theta +\frac{c_s^2}{1+w}\frac{k^2}{a}\delta -\frac{k^2}{a}\sigma +\frac{k^2}{a}\Psi ,$$

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

$$V^{\prime }$$

$$-3\cdot 10^{-15}\le c_g/c-1\le 7\cdot 10^{-16},$$

$$G_{4X}\approx 0,G_5\approx \text{const.},$$

$$c_T^2=\frac{G_4-X{G}_{5\varphi }-X{G}_{5X}\ddot{\varphi }}{G_4-2X{G}_{4X}-X\left(G_{5X}\dot{\varphi }H-G_{5\varphi }\right)}.$$

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

$${\mathcal{L}}_2$$

$${\mathcal{L}}_3$$

$${\mathcal{L}}_4$$

$$K$$

$$G_4$$

$$K$$

$$G_4$$

$$K$$

$$G_3$$

$$G_4$$

$$K$$

$$G_4$$

$$G_3$$

$$K$$

$$G_3$$

Full text

††preprint: IFT-UAM/CSIC-19-051

I Introduction

Not long ago, measurements of the distance-redshift relation from distant Supernovae type Ia (SNIa) revealed that the Universe is not only expanding as time goes by, it is actually accelerating Riess et al. (1998); Perlmutter et al. (1999). The consequences of these observations are far-reaching and several analyses of the data sets have been carefully carried out ever since (see Ref. Nielsen et al. (2016) and references therein). Although there were some concerns about possible systematic errors, analyses of new and improved data sets have shown that results in previous works are robust Riess et al. (2007); Betoule et al. (2014). Moreover, most recent astrophysical measurements of the Cosmic Microwave Background (CMB) anisotropies and the distribution of galaxies in the Universe, when interpreted in the context of the cosmological constant cold dark matter model ($\Lambda$CDM), are in very good agreement with a late-time accelerating phase Aghanim et al. (2018); Abbott et al. (2018).

Current Bayesian analyses of astrophysical measurements indicate that $\Lambda$CDM beats alternative models Heavens et al. (2017). In spite of being successful at fitting most data sets, $\Lambda$CDM is just a very good phenomenological model as its main constituents are either unknown or misunderstood. First, Cold Dark Matter (CDM) has not been directly detected thus far despite the huge effort this research field has attracted over the past years Bertone and Hooper (2018). Second, there exists an important disagreement between both predicted and inferred values of the cosmological constant $\Lambda$ whose solution will possibly lead to new physics Weinberg (1989); Carroll (2001).

Even though reconciling the quantum field theory prediction with the observed value of the cosmological constant seems unlikely, it has become clear that a Dark Energy (DE) component resembling a cosmological constant not only can alleviate several problems present in a CDM model but also can drive the current accelerating expansion of the Universe Kofman and Starobinsky (1985). Although several mechanisms have been proposed in the literature which could be responsible for speeding up the Universe, nowadays there are two main approaches. On the one hand, one finds Modified Gravity (MG) models Clifton et al. (2012). Einstein’s Theory of General Relativity (GR), the theory of gravity that is assumed in the $\Lambda$CDM model, seems to break down on tiny scales and possibly will require modifications on large scales to account for current observations Bertschinger (2011). However, modifying GR can be laborious as several tests carried out up to cosmological scales agree very well with GR Bertotti et al. (2003); Reyes et al. (2010); Collett et al. (2018); Abbott et al. (2016); He et al. (2018); Delva et al. (2018); Herrmann et al. (2018); Ishak (2019); Luna et al. (2018); Basilakos et al. (2018); Perez-Romero and Nesseris (2018); Basilakos and Nesseris (2017); Nesseris et al. (2017); Basilakos and Nesseris (2016). On the other hand, there are DE models Copeland et al. (2006) which rely on yet unobserved scalar fields that would dominate the energy content of the Universe at late times and also avoid fine-tuning issues Ratra and Peebles (1988); Armendariz-Picon et al. (2000).

Although DE and MG models are clearly motivated by different underlying physics, it is possible to study both kinds of models on the same footing. In an effective fluid approach departures from GR can be interpreted as an effective fluid contribution in such a way that comparison with DE models might become relatively simple Kunz and Sapone (2007); Pogosian et al. (2010); Capozziello et al. (2006a, b); Capozziello et al. (2019). When interpreted as fluids, MG models can be described by an equation of state $w(a)$, a sound speed $c_{s}^{2}(a,k)$, and an anisotropic stress $\pi(a,k)$: background is affected by the behavior of $w(a)$ while perturbations are mainly governed by $c_{s}^{2}(a,k)$ and $\pi(a,k)$. Since both DE and MG models predict different behavior for these three functions, in an effective fluid approach different models can be, to a certain degree, distinguished.

It is well known that both DE and MG models can accommodate background astrophysical observations as well as the standard cosmological model $\Lambda$CDM (e.g., the so-called designer $f(R)$ models Multamaki and Vilja (2006); de la Cruz-Dombriz and Dobado (2006); Pogosian and Silvestri (2008); Nesseris (2013)). As a consequence, these models are degenerated at the background level even though there have been various attempts to disentangle them by using model independent approaches Nesseris and Shafieloo (2010); Nesseris and Garcia-Bellido (2012). Fortunately, the study of linear order perturbations might break this degeneracy because DE and MG models predict different growths of structures and could in principle be distinguishable from $\Lambda$CDM Tsujikawa (2007); Pogosian and Silvestri (2008).

Given the wide range of both DE and MG models it is useful to have a unified framework which encompasses several of them. It turns out that such a theory exists since 1974 when Horndeski found the most general Lorentz-invariant extension of GR in four dimensions Horndeski (1974). This theory can be obtained by using a single scalar field and restricting the equations of motion to being second order in time derivatives. The Horndeski Lagrangian comprehends theories such as Kinetic Gravity Braiding, Brans-Dicke and scalar tensor gravity, single field quintessence and K-essence theories, as well as $f(R)$ theories in their scalar-tensor formulation Baker et al. (2013). Although the range of models encompassed by the Horndeski Lagrangian was severely reduced (see, for instance, Creminelli and Vernizzi (2017); Sakstein and Jain (2017); Ezquiaga and Zumalacarregui (2017); Baker et al. (2017); Amendola et al. (2018); Crisostomi and Koyama (2018); Frusciante et al. (2018); Kase and Tsujikawa (2018); McManus et al. (2016); Lombriser and Taylor (2016); Copeland et al. (2019); Noller and Nicola (2018a); de Rham and Melville (2018)) with the recent discovery of gravitational waves by the LIGO Collaboration Abbott et al. (2017a), an interesting remaining subclass of models (including $f(R)$ theories Sotiriou and Faraoni (2010); De Felice and Tsujikawa (2010); Nojiri et al. (2017); Nojiri and Odintsov (2011) and Kinetic Gravity Braiding Deffayet et al. (2010)) is well worth an investigation.

Recently we employed an effective fluid approach to study $f(R)$ theories Arjona et al. (2019). Even though it is not easy to obtain expressions for quantities describing perturbations (e.g., pressure perturbation $\delta P$) in MG models Kunz and Sapone (2006), by using the quasistatic and subhorizon approximations we found analytical expressions for the effective DE perturbations as well as the quantities describing the effective DE fluid, namely, $w(a)$, $c_{s}^{2}(a,k)$, and $\pi(a,k)$. We implemented our approach in the code CLASS111http://class-code.net/ Blas et al. (2011) and found excellent agreement with the so-called Equation of State (EOS) approach Battye et al. (2016, 2018), which does not use any approximation. In this paper we extend our work Arjona et al. (2019) to the remaining part of the Horndeski Lagrangian which contains $f(R)$ theories as a special case. Horndeski theories have been implemented in the code hi_CLASS Zumalacarregui et al. (2017) which solves the full set of dynamical equations without using the quasistatic approximation. In our approach we find analytical expressions for the effective DE perturbations that give us a better understanding of the underlying physics and also allow us to compare with our numerical implementation. Moreover, we show that it is possible to find ‘designer Horndeski theories’ matching a given background evolution. We implement one such a model in the hi_CLASS code and show there is good agreement with our approach, namely, our effective fluid approach assuming both quasistatic and subhorizon approximations performs quite well.

The paper is organized as follows. In Sec. II we discuss the equations for perturbations in a Friedmann-Lemaitre-Robertson-Walker (FLRW) metric and set our notation. Then, we introduce the Horndeski Lagrangian and discuss both background and perturbation equations in Sec. III. In Sec. IV we study the remaining subclass of Horndeski theories by utilizing the effective fluid approach, we discuss the subhorizon and quasistatic approximations and present analytical results for two classes of models, those in which we have dark energy anisotropic stress and those in which we do not. In Sec. V we show analytical results for a family of models named ‘designer Horndeski’ which mimic the $\Lambda$CDM background and in Sec. VI we compare our analytical solutions for DE perturbations with a fully numerical solution of the system of differential equations and show they are in very good agreement. We then constrain the parameter space for a viable designer Horndeski model in Sec. VII and in Sec. VIII we present our conclusions. In Appendices A and B we give details about our analytical computations.

II Theoretical framework

In the standard cosmological model one assumes the Einstein-Hilbert action

[TABLE]

where $g$ is the determinant of the metric $g_{\mu\nu}$, $R$ is the Ricci scalar, $\kappa\equiv\frac{8\pi G_{N}}{c^{4}}$ and $\mathcal{L}_{m}$ is the Lagrangian for matter fields.222Throughout this paper we set the speed of light $c=1$ and $\kappa=8\pi G_{N}$ with $G_{N}$ being the bare Newton’s constant. Our conventions are: $(-+++)$ for the metric signature, the Riemann and Ricci tensors are given respectively by $V_{b;cd}-V_{b;dc}=V_{a}R^{a}_{bcd}$ and $R_{ab}=R^{s}_{asb}$. Applying the principle of least action to Eq. (1) one obtains the field equations

[TABLE]

where $G_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R$ is the Einstein tensor and $T_{\mu\nu}^{(m)}$ is the energy-momentum tensor for matter fields. At this point one needs to make more assumptions about the geometrical properties and the matter content in the Universe. First, since observations indicate the Universe on large scales is statistically homogeneous and isotropic Hogg et al. (2005); Ade et al. (2016) (also having tiny inhomogeneities which can be treated within linear perturbation theory), one further assumes a perturbed FLRW metric

[TABLE]

where $a$ is the scale factor, $\vec{x}$ represents spatial coordinates, $t$ is the cosmic time and $\Psi$ and $\Phi$ are the gravitational potentials in the Newtonian gauge. Second, one can suppose the matter fields are ideal fluids (with small perturbations) having an energy-momentum tensor given by

[TABLE]

where $P$ is the pressure, $\rho$ is the energy density, and $U^{\mu}=\left(1-\Psi,\frac{\vec{u}}{a(t)}\right)$ is the velocity four-vector. As a result, the elements of the energy-momentum tensor up to first order are given by :

[TABLE]

where $\bar{\rho}$ is the background energy density, $\bar{P}$ is the background pressure, $u_{i}=a(t)\dot{x_{i}}$, $\Sigma^{i}_{j}(\vec{x},\tau)\equiv T^{i}_{j}-\delta^{i}_{j}T^{k}_{k}/3$ is an anisotropic stress tensor, and $\delta\rho(\vec{x},\tau)$ and $\delta P(\vec{x},\tau)$ are the density and pressure perturbations, respectively.333In our notation, a dot over a function $f$ denotes the derivative with respect to the cosmic time : $\dot{f}\equiv\frac{df}{dt}$. In addition, Greek indices run from [math] to $3$ whereas Latin indices take on values from $1$ to $3$.

II.1 Background

If one only considers zero order quantities in the Einstein field equations (2), then there are two independent Friedmann equations describing the background evolution of the Universe:

[TABLE]

where $H\equiv\frac{\dot{a}}{a}$ is the cosmic Hubble parameter.444The conformal Hubble parameter $\mathcal{H}$ and the Hubble parameter $H$ are related via $\mathcal{H}=aH$.

II.2 Linear perturbations

Considering just the first order perturbations in the Einstein field equations (2) we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we defined the velocity $\theta\equiv ik^{j}u_{j}$ and wrote the anisotropic stress as $(\bar{\rho}+\bar{P})\sigma\equiv-(\hat{k}_{i}\hat{k}_{j}-\frac{1}{3}\delta_{ij})\Sigma^{ij}$.

From the conservation of the energy-momentum tensor $T^{\mu\nu}_{;\nu}=0$ one obtains the equations for the evolution of perturbations. Defining the equation of state parameter as $w\equiv\frac{\bar{P}}{\bar{\rho}}$ and the sound speed $c_{s}^{2}\equiv\frac{\delta P}{\delta\rho}$ we find the equations governing the evolution of density and pressure perturbations are given by

[TABLE]

[TABLE]

The system of differential equations (14)-(15) presents problems when the equation of state crosses $-1$ because there is a singularity. However, a simple change of variable turns out to be helpful in solving this inconvenience. We will use the scalar velocity perturbation $V\equiv ik_{j}T^{j}_{0}/\rho=(1+w)\theta$ instead of the velocity $\theta$. In terms of this new variable the evolution equations (14)-(15) become

[TABLE]

[TABLE]

where a prime ′ denotes a derivative with respect to the scale factor and we defined the anisotropic stress parameter $\pi\equiv\frac{3}{2}(1+w)\sigma$.

III Horndeski

Horndeski theory constitutes the most general Lorentz-invariant extension of GR in four dimensions and encompasses several DE and MG models. Although in its most general form the Horndeski Lagrangian has several free functions, the recent discovery of gravitational waves by the LIGO Collaboration significantly constrained the allowed models. In particular, it has been shown that the constraint on the speed of Gravitational Waves (GWs) must satisfy Ezquiaga and Zumalacarregui (2017)

[TABLE]

which for Horndeski theories implies that

[TABLE]

as can be seen from the sound speed formula for tensor perturbations Kobayashi et al. (2011)

[TABLE]

In this section we will derive evolution equations for the remaining parts of the Horndeski Lagrangian, namely,

[TABLE]

where

[TABLE]

Here $\phi$ is a scalar field, $X\equiv-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi$ is a kinetic term, and $\Box\phi\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\phi$; $K$, $G_{3}$ and $G_{4}$ are free functions of $\phi$ and $X$.555From now on we define $G_{i}\equiv G_{i}\left(\phi,X\right)$, $G_{i,X}\equiv G_{iX}\equiv\frac{\partial G_{i}}{\partial X}$ and $G_{i,\phi}\equiv G_{i\phi}\equiv\frac{\partial G_{i}}{\partial\phi}$ where $i=2,3,4$. Since we are mainly interested in the late-time dynamics of the Universe, hereafter we will further assume $\mathcal{L}_{m}$ is the Lagrangian of a CDM component. As has been mentioned in Frusciante et al. (2018), although the functions $K$, $G_{3}$ and $G_{4}$ are able to modify the background with a general dependence on $X$ and $\phi$, this does not hold at the perturbations level. For instance, $K(\phi,X)$ encloses the k-essence and quintessence theory and is partly responsible for the background and the perturbations, however $K(\phi)$ does not contribute to the perturbations.

The term $G_{3}(\phi,X)$ includes the kinetic gravity braiding with $G_{3X}\neq 0$ being in charge of combining the kinetic term of the scalar and the metric, but the term $G_{3}(\phi)$ only modifies the background as a dynamical dark energy. Finally, $G_{4}$ is the only function that is able to modify the non-minimal coupling of the scalar to the Ricci curvature.

Among the theories embedded in the action (21) one finds, for example:

• •

f(R) theories. When interpreted as a non-minimal coupled scalar field, these theories can be written using Chiba (2003)

[TABLE]

where $\phi\equiv\frac{f_{,R}}{\sqrt{\kappa}}$ has units of mass and $f_{,R}\equiv\dfrac{df}{dR}$.

• •

Brans-Dicke theories. In our notation we have

[TABLE]

where $V(\phi)$ is the field potential and $\omega_{BD}$ is the Brans-Dicke parameter Brans and Dicke (1961).

• •

Kinetic gravity braiding. This kind of scalar-tensor models exhibit mixing of scalar and tensor kinetic terms Deffayet et al. (2010) and can be written as

[TABLE]

• •

Non-minimal coupling (NMC) model Quiros (2019). In our notation and for a coupling constant $\zeta$

[TABLE]

In the context of inflation, a Higgs-like inflation model corresponds to $\omega(\phi)=1$, $V(\phi)=\lambda\left(\phi^{2}-\nu^{2}\right)^{2}/4$.

• •

Cubic Galileon Quiros (2019). The simplest case is when

[TABLE]

• •

4-dimensional static and spherical symmetric solution of Black Hole with scalar hair Fang et al. (2018).

[TABLE]

As previously done for the Einstein-Hilbert action (1), here we apply the principle of least action to (21) in order to find evolution equations for both the gravitational field and the scalar field. Varying Eq. (21) with respect to the metric and the scalar field one finds666See Appendix A for a derivation of the field equations. Kobayashi et al. (2011)

[TABLE]

which allows us to find the field equations. First, the gravitational field equation is given by

[TABLE]

where we have defined

[TABLE]

and $T_{\mu\nu}^{(m)}$ is the energy-momentum tensor of a CDM component. Note that from Eq. (42) we retrieve the GR field equations (2) if we set $K=G_{3}=0$ and $G_{4}=\frac{1}{2\kappa}$. Second, the scalar field equation reads

[TABLE]

where

[TABLE]

As it is mentioned in Ref. Kobayashi et al. (2011), one could think $\nabla^{\mu}J^{i}_{\mu}$ leads to higher than second-order derivatives. However, this is not the case since commutations of higher derivatives can be substituted by the curvature tensor and are hence canceled. In particular, one can prove that

[TABLE]

which will be of paramount importance when we will discuss perturbation equations.

It is possible to find a relatively simple expression for the scalar field equation (46) if we consider the case $i=3$, namely,

[TABLE]

The terms on top of the brace in Eq. (III) can be expanded as

[TABLE]

and the terms on top of the brace in Eq. (55) can in turn be written as

[TABLE]

Using Eq. (III) in Eq. (III) we find

[TABLE]

and the scalar field equation (46) can be written as

[TABLE]

In what follows, in order to simplify the notation we will denote the kinetic term of the scalar field evaluated at the background simply by $X$ and its linear order perturbation as $\delta X$.

III.1 Background

Thus far the discussion of the field equations has been quite general. Now, as previously done in Sec. II, we assume a perturbed FLRW as given in Eq. (3). If we consider only zero order quantities in the gravitational field equation (42), we obtain

[TABLE]

where

[TABLE]

Eqs. (59)-(60) are the modified Friedmann equations describing the background evolution of the Universe. Collecting terms they respectively read

[TABLE]

[TABLE]

Note that from Eqs. (67)-(68) we respectively retrieve the Friedmann equations (8)-(9) if we set $K=G_{3}=0$ and $G_{4}=\frac{1}{2\kappa}$. Rearranging terms in Eqs. (67)-(68) we can define an effective DE density

[TABLE]

and an effective DE pressure

[TABLE]

in such a way that we can write the modified Friedmann equations Eqs. (67)-(68) as

[TABLE]

where we are assuming that matter is pressureless $\bar{P}_{m}=0$ as indicated by current constraints Kopp et al. (2018). The effective DE density and pressure in Eqs. (III.1)-(III.1) allow us to define an effective DE equation of state as

[TABLE]

Let us now consider the scalar field equation (58) and only keep zero order quantities, that is to say,

[TABLE]

which fully agrees with Kimura and Yamamoto (2011). Note that defining

[TABLE]

we can write the scalar field equation (46) as

[TABLE]

and it becomes clear that there exists a Noether current for Lagrangians invariant under constant shifts of the field $\phi\rightarrow\phi+c$ Deffayet et al. (2010), namely,

[TABLE]

Taking into consideration that $X=\frac{1}{2}\dot{\phi}^{2}$, the charge density of the Noether current can be written as

[TABLE]

so that the scalar field equation is given by the simple expression

[TABLE]

When $P_{\phi}=0$ then it is easy to see that the solution to the previous equation is

[TABLE]

where $J_{c}$ is a constant. When $J_{c}=0$, then the system is on the attractor solution, but when $J_{c}\neq 0$ then the system is not on the attractor and as we will see in Sec. IV.2.3 interesting dynamics may arise.

III.2 Linear perturbations

Considering only first order quantities in the gravitational field equations (42) one obtains De Felice et al. (2011); Matsumoto (2019)

[TABLE]

Note that when $K=G_{3}=0$ and $G_{4}=\frac{1}{2\kappa}$, Eqs. (82)-(85) respectively correspond to the GR limit given by Eqs. (10)-(13) with $\sigma_{m}=0$.

If we now consider the scalar field equation (58) and take into account only first order quantities we find

[TABLE]

Expressions for the coefficients $A_{i}$, $\mu$, $\nu$, $B_{i}$, $C_{i}$ and $D_{i}$ can be found in Appendix B and are in agreement with those found in De Felice et al. (2011); Matsumoto (2019), except for $D_{8}$ which is actually equal to zero as can be seen by using the expression found in Matsumoto (2019) and using the background equations of motion for the scalar field.

IV The effective fluid approach

We have seen in the previous section that the gravitational field equations for the Horndeski Lagrangian can be written in such a way that they resemble those found in Sec. II where we assumed GR and a perfect fluid. Indeed, defining an effective DE density and pressure given by Eqs. (III.1)-(III.1) makes it possible to obtain an effective DE equation of state (see Eq. (73)). As mentioned in Sec. I, a fluid can be described by its equation of state, sound speed, and anisotropic stress, so in what follows we will explicitly derive those quantities.

In this section we will present relatively simple expressions for the effective DE sound speed and anisotropic stress under the subhorizon and quasistatic approximations. Actually, by defining an effective DE fluid we are considering a DE effective energy-momentum tensor $T_{\mu\nu}^{DE}$ obtained via the gravitational field equations (42) and defined explicitly as follows:

[TABLE]

Since we are taking into consideration expressions up to linear order, $T_{\mu\nu}^{DE}$ also contains small perturbations which allow us to define quantities such as DE effective perturbations in the pressure, density, and velocity. These can be extracted from the DE effective energy-momentum tensor $T_{\mu\nu}^{DE}$ by considering the decomposition of the tensor into its components, given by Eqs. (5)-(7). Qualitatively, these expressions have the following structure:

[TABLE]

where $(...)$ indicates expressions which might be cumbersome. It is therefore very helpful to work out these expressions under the subhorizon and quasistatic approximations in order to gain a better understanding.

We have explained in great detail the way we carry out the subhorizon and quasistatic approximations in our previous paper (see Sec. II.A.1 in Ref. Arjona et al. (2019)), but in a nutshell, the former refers to only considering modes deep in the Hubble radius, i.e. those for which $k^{2}\gg a^{2}H^{2}$, while the latter refers to neglecting derivatives of the potentials during matter domination as they are roughly constant but also terms of similar order as $\partial_{\eta}\sim 1/\eta\sim aH(a)$. For example, the perturbation in the Ricci scalar is

[TABLE]

Following the same procedure and applying the subhorizon approximation to the linearized gravitational field equations (82),(84), and to the linearized scalar field equation (III.2), one finds, respectively,

[TABLE]

Note that since $B_{7}=4G_{4\phi}$ and $B_{6}=B_{8}$ (see Appendix B), Eq. (92) leads to no anisotropic stress $\Phi=-\Psi$ when $G_{4}$ is a constant.

Solving Eqs. (91)-(93) for $\Phi$, $\Psi$ and $\delta\phi$ one finds

[TABLE]

where $G_{\textrm{eff}}$ and $Q_{\textrm{eff}}$ are Newton’s effective constant

[TABLE]

and we make use of the following correspondence $A_{3}=B_{6}=B_{8}$, $D_{7}=B_{7}$ and $D_{10}=A_{6}$ (see Appendix B). One can also define the following anisotropic stress parameters

[TABLE]

The aforementioned expressions for Newton’s effective constant and the anisotropic stress parameters are in agreement with the ones in Ref. De Felice et al. (2011).

The subhorizon approximation is also useful as the evolution equations for the growth of matter perturbations $\delta_{m}$ given by Eqs. (16)-(17) can be reduced to a single differential equation, where the variable $G_{\textrm{eff}}$ plays a primary role:

[TABLE]

with $G_{\textrm{eff}}$ given by Eq. (97) and initial conditions $\delta_{m}(a_{i})=a_{i}$ and $\delta_{m}^{\prime}(a_{i})=1$ for an initial value for the scale factor $a_{i}$ deep in the matter era.

In what follows, we will present the effective DE perturbations under the subhorizon and quasistatic approximations for two classes of models: those in which there is DE anisotropic stress and those where DE anisotropic stress vanishes.

IV.1 Horndeski models with DE anisotropic stress

We now apply the subhorizon and quasistatic approximations in Eqs. (88)-(90) using the same prescription as in Ref. Arjona et al. (2019). We also found, in agreement with Ref. Kimura and Yamamoto (2011), that the quasistatic approximation breaks down for this model due to the rapid oscillations of the scalar field, so if we eliminate the scalar field, then this can slightly increase the accuracy of the numerical solutions of the effective fluid equations. To eliminate $\delta\phi$ and its derivatives, we use Eq. (85) and insert the resulting equations in Eqs. (88)-(90).

Then, by keeping the dominant $k^{2}$ terms (the subhorizon approximation) and dropping time derivatives of the potentials (the quasistatic approximation) in Eqs. (88)-(90) we find

[TABLE]

for the effective DE pressure perturbation, effective DE density perturbation, and effective DE velocity perturbation, respectively (the interested reader can find the expressions for $\mathcal{F}_{i}$ in Appendix B). It is now also possible to obtain an expression for the effective DE anisotropic stress under the subhorizon approximation

[TABLE]

Having found expressions for the effective DE equation of state (see Eq. (73)) and the effective DE anisotropic stress (Eq. (105)), the only missing ingredient for an effective fluid description of the Horndeski Lagrangian is the sound speed. This quantity can easily be found using our equations for the effective DE pressure perturbation (102) and the effective DE density perturbation (103). The DE sound speed reads

[TABLE]

Due to the presence of anisotropic stress, perturbations on subhorizon scales in the effective DE fluid are not driven by the sound speed (106), but by an effective DE sound speed defined as Cardona et al. (2014); Arjona et al. (2019)

[TABLE]

Finally, it is clear that for the cosmological constant model, i.e. $\mathcal{L}_{2}=-\frac{\Lambda}{\kappa}$, $\mathcal{L}_{3}=0$, $\mathcal{L}_{4}=\frac{1}{2\kappa}R$, $\mathcal{L}_{5}=0$, we have $K=-\frac{\Lambda}{\kappa}$, $G_{3}=0$, $G_{4}=\frac{1}{2\kappa}$ and $G_{5}=0$, which implies that $w_{DE}=-1$ and $\left(\delta P_{DE},\delta\rho_{DE},\pi_{DE}\right)=\left(0,0,0\right)$ as expected.

IV.1.1 f(R) models

Thus far we have kept the discussion quite general, that is to say, we did not specify any function in the Horndeski Lagrangian (21). To mention an example, we will present the results for $f(R)$ models. With the definitions in Eqs. (25)-(26) and using units where $\kappa=1$, one obtains

[TABLE]

where

[TABLE]

and $F=f_{,R}$, $F_{,R}=f_{,RR}$. Then, the effective DE fluid quantities read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

These results are in perfect agreement with our previous work Arjona et al. (2019).

IV.2 Horndeski models with no dark energy anisotropic stress

With the same approach that we followed in (IV.1) we compute the DE perturbations for models where the is no DE anisotropic stress, i.e $\Phi=-\Psi$. With this restriction it is easy to see from Eq. (85) that $G_{4\phi}=0$. Then applying this condition under the subhorizon approximation in Eqs. (88)-(90) leads to

[TABLE]

and since $\Phi=-\Psi$ the anisotropic parameters read

[TABLE]

as expected, while the DE anisotropic stress parameter is zero $\pi_{DE}=0$. Our general expression for the DE sound speed (106) reduces in this case to

[TABLE]

which is equal to the DE effective sound speed since $\pi_{DE}=0$. Here we will show results for a few specific models embedded in the Horndeski Lagrangian.

IV.2.1 Quintessence

We can recover the Lagrangian of Quintessence by choosing the following functions

[TABLE]

where $\phi$ is the scalar field, $X$ is the kinetic term defined as $X=-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$ and $V(\phi)$ is the potential. Using a variational approach one finds that the effective pressure, density and velocity perturbations for Quintessence theories are given by

[TABLE]

and these expressions are in agreement with Amendola and Tsujikawa (2015). Also, the DE anisotropic stress parameter $\pi_{DE}$ is zero since for Quintessence $\Psi=-\Phi$. We find that under the subhorizon approximation

[TABLE]

so that the effective pressure, density and velocity perturbations for Quintessence theories are given by

[TABLE]

It is thus straightforward, using Eqs. (127) and (128), to see that the DE sound speed is given by

[TABLE]

Moreover, we also find that in the subhorizon approximation

[TABLE]

IV.2.2 K-essence

In our notation the Lagrangian of K-essence theories is specified by the functions Scherrer (2004); de Putter and Linder (2007)

[TABLE]

and as usual through the variation of the action it is possible to find expressions for the pressure, density, and velocity perturbations

[TABLE]

Since for K-essence $\Psi=-\Phi$ the DE anisotropic stress parameter $\pi_{DE}$ vanishes. We find that under the subhorizon approximation

[TABLE]

and therefore the DE perturbations for K-essence theories are given by

[TABLE]

and the DE sound speed reads

[TABLE]

in agreement with Refs. Scherrer (2004); de Putter and Linder (2007). The perturbations of the scalar field and the gravitational potential are respectively given by

[TABLE]

IV.2.3 Kinetic gravity braiding

An interesting DE model is the kinetic gravity braiding (KGB) which is characterized by the following Lagrangian

[TABLE]

Since $G_{4}$ is constant it is easily shown from Eq. (92) that the KGB model has no DE anisotropic stress and therefore the anisotropic parameters

[TABLE]

Furthermore, it follows that the effective Newton’s constant $G_{\textrm{eff}}/G_{N}$ is given by

[TABLE]

The effective DE density and pressure $\bar{\rho}_{DE}$ and $\bar{P}_{DE}$ read, respectively,

[TABLE]

and therefore the DE equation of state is given by

[TABLE]

We also find that the scalar field equation at the background level is

[TABLE]

As a specific example we now discuss the KGB model of Ref. Kimura and Yamamoto (2011) defined by

[TABLE]

where $n$ and $\alpha$ are parameters in the model. A number of reasons make the KGB an attractive model. First, it passes the recent observational constraints from gravitational waves. Second, it is known that this model connects the original Galileon model Deffayet et al. (2010) and the $\Lambda\text{CDM}$ model by the parameter $n$, at least for the background and first order perturbations: linear perturbations of the KGB model reduce to those of $\Lambda$CDM (original Galileon) for $n=\infty$ ($n=1$) Kimura and Yamamoto (2011).

The charge density of the Noether current Eq. (79) is in this case

[TABLE]

and satisfies the differential equation

[TABLE]

whose solution reads

[TABLE]

with $J_{c}$ a constant. It is therefore clear that $J_{0}$ approaches zero as the Universe expands. The simplest attractor solution is located at $J_{0}=0$ and has two branches, namely,

[TABLE]

and

[TABLE]

Because the first case has ghostly perturbations, as it is shown in Kimura and Yamamoto (2011), we will focus on the attractor solution Eq. (156). Using Eqs. (71) and (152) we find that the modified Friedmann equation is given by

[TABLE]

where we have neglected radiation. The background equation of the KGB model reduces to that of $\Lambda\text{CDM}$ for $n=\infty$ as can be seen from Eq. (157). Also, one can easily find an expression for the parameter $\alpha$ by using Eq. (157) at the present epoch

[TABLE]

The DE equation of state becomes

[TABLE]

and through Eq. (156) it is also possible to find an analytical expression for the kinetic term

[TABLE]

where the prime stands for the derivative with respect to the scale factor.

To derive the $\Lambda$CDM limit for the perturbations in this model we rewrite Eqs. (84) and (III.2) in terms of the kinetic term perturbation $\delta X=\dot{\phi}\dot{\delta\phi}-\dot{\phi}^{2}\Psi$. Then, for $n\rightarrow\infty$ the former equation reduces to $\delta X\left(-\frac{2}{a}-\frac{H^{\prime}(a)}{H(a)}+O(1/n)\right)$, while the latter equation gives

[TABLE]

which implies that the kinetic term perturbation decays as $\delta X\sim 1/a^{3}$ and thus can be ignored at late time. Since DE perturbations in the KGB model are proportional to $\delta X$ for large $n$, then they reduce to zero as expected for the $\Lambda$CDM model.

Finally, it should be noted that a standard hydrodynamical description of the KGB in terms of an effective fluid, has been studied in Ref. Pujolas et al. (2011). There, it was shown that the KGB model can also be described in terms of an imperfect fluid with a chemical potential, in which the equations of motion reduce to the standard diffusion equation. However, in our current analysis we will only focus on the ideal fluid approach, which is totally equivalent, as we are interested in finding simple analytic solutions and with comparing with our previous work.

V Designer Horndeski

In this section we will address the shortcomings found in the KGB model defined by Eqs. (150)-(151). We will show that it is possible, starting from the Lagrangian (142), to find a model corresponding to a given background but yet having different perturbations. Using the modified Friedmann equation and the scalar field conservation equation, we can find specific designer models such that the background is always that of the $\Lambda$CDM model, namely, having $w_{DE}=-1$. This is particularly useful in detecting deviations from $\Lambda$CDM at the perturbations level and is a natural expansion of our earlier work Nesseris (2013); Arjona et al. (2019). We start with the modified Friedmann equation, which can be written as

[TABLE]

while the scalar field conservation equation can be written as

[TABLE]

where $J_{c}$ is a constant which quantifies our deviation from the attractor, as in the case of the KGB model Kimura and Yamamoto (2011). We now have two equations given by (V) and (163), but three unknown functions $(G_{3X}(X),K(X),H(a))$ thus the system is undetermined. Therefore, we need to specify one of the three unknown functions $(G_{3X}(X),K(X),H(a))$ and determine the other two using Eqs. (V) and (163). To facilitate this, we express the Hubble parameter as a function of the kinetic term $X$, ie $H=H(X)$ and then solve the previous equations to find $(G_{3X}(X),K(X))$. Doing so yields:

[TABLE]

With Eqs. (V) we can make a whole family of designer models that behave as $\Lambda$CDM at the background level but have different perturbations. We now proceed to specify some examples using our formalism.

V.1 Example 1

Choosing $K(X)=0$ and solving Eqs. (V) we find

[TABLE]

and the derivative of the scalar field $\phi^{\prime}(a)$ is

[TABLE]

where the prime is the derivative with respect to the scale factor. However, this model has the problem that it does not have a smooth limit to $\Lambda$CDM when $J_{c}=0$.

V.2 Example 2

On the other hand, specifying $G_{3}(X)$ leads to another interesting designer model, defined as

[TABLE]

where the kinetic term is defined as

[TABLE]

and $X_{0}$ is an integration constant. However, this model has the problem that at early times the perturbations do not go to zero and we do not recover GR, since the kinetic term goes to infinity as it grows as $X\sim H(a)$.

V.3 Example 3 (HDES)

To solve the previous shortcomings we follow a different approach. First, we demand that the kinetic term behaves as $X=\frac{c_{0}}{H(a)^{n}}$, where $c_{0}>0$ and $n>0$. Then, from Eqs. (163) and (V) we find:

[TABLE]

This specific model solves both previous problems, i.e., it has a smooth limit to $\Lambda$CDM and it also recovers GR when $J_{c}\sim 0$, thus we will designate this model as HDES and focus on it in what follows.

V.4 Comparison with the $\alpha$ parameters

To facilitate comparisons with the literature we also provide the expressions for our designer HDES model in terms of the $\alpha_{i}$ functions, where $i=M,K,B,T$. The functions $G_{i}(\phi,X)$ and $\alpha_{i}$ are connected in the following manner Zumalacarregui et al. (2017):

[TABLE]

where the dot is the derivative with respect to the cosmic time, $M_{\ast}^{2}(\tau)$ is the cosmological strength of gravity, $\alpha_{T}$ is the tensor speed excess, $\alpha_{B}$ is called the braiding and $\alpha_{K}$ is referred to as the kineticity. For more information on these $\alpha_{i}$ functions see Ezquiaga and Zumalacarregui (2018). At all times we require $D=\alpha_{K}+\frac{3}{2}\alpha_{B}^{2}>0$ so that there are no ghostly instabilities and that $\alpha_{M,K,B,T}\simeq 0$ at early times, so as to recover GR.

For our HDES designer model given by Eqs. (169), we have that the $\alpha_{i}$ functions of Eq. (170) are given by

[TABLE]

Since in Eqs. (173)-(174) we have a degeneracy with the coefficients $c_{0}$ and $J_{c}$, they appear together as $\sqrt{c_{0}}J_{c}$, we can choose to absorb $c_{0}$ in the definition of $J_{c}$. Finally, it is straightforward to see that our $\alpha_{i}$ functions are dimensionless since through dimensional analysis we found that $[c_{0}]=H^{n+2}_{0}$, $[J_{c}]=H_{0}$, the kinetic term $[X]=H^{2}_{0}$, $[K]=H^{2}_{0}$ and $[G_{3X}]=H^{-2}_{0}$.

Notice that not all designer models satisfy the above conditions, so in what follows we consider only HDES, given by Eq. (169). Then, the stability condition $D=\alpha_{K}+\frac{3}{2}\alpha_{B}^{2}>0$ for our model Eq. (169) gives

[TABLE]

where we have set $\tilde{J}_{c}=J_{c}/H_{0}$ and $\tilde{c}_{0}=c_{0}/H^{n+2}_{0}=1$.

Then, inequality (176) implies that in order for the system to be stable we must have either $\tilde{J}_{c}>0$ for $0<n\leq 2$ or a complicated set of expressions that can however be easily derived from Eq. (176) with algebraic manipulations. For $n=2$ the inequality is automatically satisfied for any value of $\tilde{J}_{c}$ as $\alpha_{K}=0$ as can be seen from Eq. (173). We show the complicated parameter space that is allowed for $n=1$ and $n=3$ as a function of scale factor $a$ but also as a function of $n$ for $a=1$, in Fig. 1.

V.5 Analytic solutions for the growth

Furthermore, in this case we can also find approximate solutions to the growth equation Eq. (101) in matter domination for $n=2$. To do this, we first do a series expansion around $a=0$ to the $G_{\textrm{eff}}$ of Eq. (97), which gives:

[TABLE]

which we can use to solve Eq. (101) in matter domination, where $H(a)/H_{0}\simeq\sqrt{\Omega_{m,0}a^{-3}}$. Then, we get

[TABLE]

where $I_{n}(z)$ is the modified Bessel function of the first kind and $\Gamma(n)$ is the usual Gamma function. Using Eq. (178) and the definition of the growth rate $f\sigma_{8}(a)\equiv f(a)\cdot\sigma(a)=\sigma_{8}a\delta_{m}^{\prime}(a)/\delta_{m}(a=1)$, we can calculate the latter exactly. However, it is instructive to perform a series expansion around $a=1$, which gives:

[TABLE]

where we have defined the parameters

[TABLE]

where $\,{}_{0}F_{1}(c_{1},z)$ is a hypergeometric function.

As can be seen from Eq. (179) there is a strong degeneracy between $\tilde{J}_{c}$ and $\sigma_{8}$, which can also be demonstrated by doing a series expansion of $f\sigma_{8}(a=1)$ for small $\tilde{J}_{c}$, which gives

[TABLE]

which implies that if we keep the growth today given constant, i.e., $f\sigma_{8}(a=1)=C_{0}=\textrm{const}.$ then $\sigma_{8}$ will scale roughly as

[TABLE]

Since $\Omega_{m,0}$ is strongly constrained from Planck, we expect that the low redshift $f\sigma_{8}$ data will exhibit a degeneracy between $\tilde{J}_{c}$ and $\sigma_{8}$. More specifically, by inspecting Eq. (183) we expect a strong negative correlation between the two parameters and this is exactly what we see from the actual Markov Chain Monte Carlo (MCMC) that we present in later sections. This degeneracy is interesting as it can potentially alleviate the soft $2\sigma$ tension between the growth rate data ($\sigma_{8}=0.88$) and Planck ($\sigma_{8}=0.831)$, which has been extensively discussed in the literature, see Ref. Nesseris et al. (2017); Sagredo et al. (2018) and references therein.

VI Numerical solutions

Here we present the numerical solutions of the two models, the KGB and HDES, that we described in the previous section.

VI.1 The KGB model

VI.1.1 The attractor

To explore the possibility of working outside the attractor we only need to use Eqs. (71) and (152), as these constrain $J_{c}$ and $\alpha$ with $H(a=1)=H_{0}$. To parameterize the deviation from the attractor we will use the parameter $J_{c}$. An illustrative example is found in Fig. 2 where we plot the dark energy density $\Omega_{DE}$ with respect to the scale factor for several values of $n$ (left) and $J_{c}$ (right). The values of values for $J_{c}$ were chosen so as to highlight the differences of these models with respect to GR.

In the KGB model the DE density can be written via Eq. (III.1) as

[TABLE]

From Fig. 2 we can see that working outside the attractor for the KGB model $(n=1)$ we might find new parts of the parameter space and new phenomenology. In the right panel of Fig. 2, we see that the orange line can be ruled out because it predicts a very high value for the DE density at early times. The red and green lines, although outside the attractor solution, are plausible solutions that are interesting to analyze in more depth.

VI.1.2 Numerical solution

In this section we present the results of the numerical solution of the evolution equations. In all cases we will assume $\Omega_{m,0}=0.3$, $k=300H_{0}$ and $\sigma_{8,0}=0.8$, unless otherwise specified. The reason we choose the specific value of $k=300H_{0}\sim 0.1\;h/\textrm{Mpc}$ for the wave-number is that it corresponds to the largest value of $k$ we can choose without entering the non-linear regime. Finally, we set the initial conditions for the DE variables to zero at $a_{i}=10^{-3}$, when we are well inside the matter dominated regime.

Next we will also present our results for the growth rate of matter perturbations parameter $f\sigma_{8}(a)\equiv f(a)\cdot\sigma(a)$, where $f(a)=\frac{dln\delta}{dlna}$ is the growth rate and $\sigma(a)=\sigma_{8,0}\frac{\delta(a)}{\delta(1)}$ is the redshift-dependent rms fluctuations of the linear density field within spheres of radius $R=8h^{-1}\textrm{{Mpc}}$, while the parameter $\sigma_{8,0}$ is its value today. The $f\sigma_{8}(a)$ parameter is important as it can be shown to be not only independent of the bias $b_{1}$, but also a good discriminator of DE models. The reason for this is that in linear theory the quadrupole contribution to the galaxy power spectrum in redshift space is sensitive only to the combination $f\sigma_{8}(a)$.

Specifically, here we will compare the numerical solutions for the following cases:

• •

The numerical solution of the full system of equations given by Eqs. (82)-(85), which however we rewrite in terms of $\delta X=\dot{\phi}\dot{\delta\phi}-\dot{\phi}^{2}\Psi$ as the system is more stable this way. We call this case “Full KGB”.

• •

The numerical solution of the effective fluid approach given by Eqs. (16)-(17). We call this case “Eff. Fluid”.

• •

The numerical solution of the growth factor equation (101). We call this case “ODE-Geff”.

• •

The $\Lambda$CDM model.

In the left panel of Fig. 3 we show the evolution of the matter and effective DE perturbation variables $(\delta_{m},V_{m},\delta_{DE},V_{DE})$ for the KGB for $n=2$. In the right panel we show the evolution of the $f\sigma_{8}(z)$ parameter for the KGB model for $n=2$ and $\sigma_{8,0}=0.8$ versus the $f\sigma_{8}$ data compilation from Ref. Sagredo et al. (2018). We show the theoretical curves for the “Full KGB” brute-force solution, the effective fluid approach, the $\Lambda$CDM model and the numerical solution of the $G_{\textrm{eff}}$ equation. As can be seen, the agreement with all approaches is excellent.

An interesting thing to note in Fig. 3 is that $V_{DE}>\delta_{DE}$ and $V_{DE}\sim V_{m}$ at intermediate redshift. The reason for this is that in the effective fluid approach the DE velocity perturbations are not always subdominant, as it would be expected in a general DE fluid. This can be seen by remembering that the velocity perturbations are actually a component of the effective energy momentum tensor, namely the $T^{0}_{i}$ part, thus they contain some of the main contributions of the Modified Gravity (MoG) theory and can be in some cases rather large. See, for example, Eqs. (6) and (17) for the definition of $V_{DE}$ and Eqs. (83) and (90) for all of the extra terms that are rewritten as $V_{DE}$.

As an example, also consider the case of quintessence and k-essense, where $V_{DE}$ is proportional to the scalar field perturbations, see Eqs. (125) and (135) respectively. In the case of $f(R)$, $V_{DE}$ is given by (113) and is proportional to $\dot{F}/F$, which parameterizes the deviations from GR, so it is a proxy for the $f(R)$ modified gravity perturbations.

However, in the case of the KGB model the subhorizon approximation fails when the parameter $n$ is large. This can easily be seen by calculating the large $n$ limit of the $G_{\textrm{eff}}$ parameter via Eq. (97):

[TABLE]

which at $a=1$ tends to $G_{\textrm{eff}}/G_{N}\simeq\frac{3}{5}+\frac{2}{5\Omega_{m,0}}$, which is different from unity as expected at this limit. However, in general deviations of $G_{\textrm{eff}}/G_{N}$ from unity on such scales are not problematic as screening mechanisms play an important role. In any case, our finding is in agreement with what was previously found in Ref. Kimura and Yamamoto (2011), namely: the quasistatic approximation breaks down for the model due to the rapid oscillations of the scalar field. As a result, in what follows we will only focus on our new designer model, which does not suffer from this issue.

VI.2 Designer Model

We now focus on our designer model HDES, given by Eq. (169). Again, we will consider the numerical solutions for the following cases:

• •

The numerical solution of the full system of equations given by Eqs. (82)-(85), which however we rewrite in terms of $\delta X=\dot{\phi}\dot{\delta\phi}-\dot{\phi}^{2}\Psi$ as the system is more stable this way. We call this case “Full-DES”.

• •

The numerical solution of the effective fluid approach given by Eqs. (16)-(17). We call this case “Eff. Fluid”.

• •

The numerical solution of the growth factor equation (101). We call this case “ODE-Geff”.

• •

The $\Lambda$CDM model.

As mentioned in the previous sections, we can absorb the constant $c_{0}$ in that of $J_{c}$, so we will only vary the latter, i.e., we set $\tilde{c}_{0}=1$. Furthermore, since the model is stable for all values of $J_{c}$ when $n=2$, we will consider this case when studying cosmological constraints. Again, we use $\Omega_{m,0}=0.3$, $k=300H_{0}$ and $\sigma_{8,0}=0.8$, unless otherwise specified.

In the left panel of Fig. 4 we show the evolution of the $f\sigma_{8}(z)$ parameter for the HDES model with $n=2$, $\tilde{J}_{c}=5\cdot 10^{-2}$ and $\sigma_{8,0}=0.8$. The values of values for $\tilde{J}_{c}$ were chosen so as to highlight the differences of these models with respect to GR. We show the theoretical curves for the HDES model for the “Full-DES” brute-force numerical solution, the effective fluid approach, the $\Lambda$CDM model and the numerical solution of the $G_{\textrm{eff}}$ equation. As can be seen, the agreement with all approaches is excellent. In the right panel of the same figure we show the percent difference between the “Full-DES” brute-force numerical solution and the effective fluid approach (magenta dot dashed line) and the numerical solution of the growth factor equation (101) (green dotted line).

VI.3 Modifications to CLASS and the ISW effect.

Here we will present our modifications to the CLASS Boltzmann code, which we call EFCLASS. We will compare the outcome with the hi_CLASS code, which solves the full set of dynamical equations but at the cost of significantly more complicated modifications. At the same time, we will also compare with a brute force calculation of the ISW effect as in our previous paper Arjona et al. (2019).

In order to modify the CLASS code in our effective fluid approach we only need two functions, the DE velocity and the anisotropic stress Arjona et al. (2019). In the case of the HDES model, the anisotropic stress $\pi_{DE}$ is zero, as can be seen from Eq. (85), since $G_{4\phi}=0$. Therefore, we only need the DE velocity which we can easily be obtained from Eq. (119), however we found that this approach is not very stable numerically. Hence, in order to have a consistent solution, we solve Eq. (17) for $V_{DE}$ and since $w_{DE}=-1$, the only variable we need is the effective pressure $\delta P_{DE}$ given by Eq. (117). The expressions are rather cumbersome, but for $n=1$ we have

[TABLE]

In the left panel of Fig. 5 we show the low-$\ell$ multipoles of the TT CMB spectrum for a flat universe with $\Omega_{m,0}=0.3$, $n_{s}=1$, $A_{s}=2.3\cdot 10^{-9}$, $h=0.7$ and $(\tilde{c_{0}},\tilde{J_{c}},n)=(1,2\cdot 10^{-3},1)$. Our EFCLASS code is denoted by the green line, hi_CLASS by the orange line and for reference the $\Lambda$CDM with a blue line. On the right panel of Fig. 5 we show the percent difference of our code with hi_CLASS as a reference777In this case we did not use $n=2$ as we found that in this case hi_CLASS crashes and we cannot compare with that code.. As can be seen, our simple modification achieves roughly $\sim 0.1\%$ accuracy across all multipoles.

We also compare our results with a brute force calculation of the Integrated Sachs-Wolfe (ISW) effect. In this case the power spectrum is given by Song et al. (2007):

[TABLE]

where $I_{\ell}^{\textrm{ISW}}(k)$ is a kernel that depends on the line of sight integral of the growth and a bessel function and $P_{\zeta}$ is the power spectrum (see Ref. Song et al. (2007) and Appendix A of Ref. Arjona et al. (2019)), and is given by the primordial power spectrum times a transfer function

[TABLE]

where $A_{s}$ is the primordial amplitude, $k_{0}$ is the pivot scale and $T(k)$ is the usual matter-radiation Bardeen, Bond, Kaiser and Szalay (BBKS) transfer function (see Eq. (7.71) in Ref. Dodelson (2003)).

In Fig. 6 we present the results for the calculation of the ISW effect and a comparison with CLASS/hi_CLASS for the $\Lambda$CDM model (left) and the HDES model (right), for the same parameters as in Fig. 5. We see that there is excellent agreement for all multipoles, except $\ell=2$. The reason for this is that we have used the BBKS formula for the transfer function $T(k)$ which is very accurate at small scales, but only at the level of $10\%$ on large scales, i.e., small multipoles.

VII Cosmological constraints

Here we present the cosmological constraints for the $n=2$ HDES and $\Lambda$CDM models discussed in previous sections. We use the latest cosmological observations including the supernovae type Ia (SnIa), Baryon Acoustic Oscillations (BAO), CMB and the Hubble expansion H(z) data. Specifically, we use the Pantheon SnIa compilation of Ref. Scolnic et al. (2018), the BAO measurements from 6dFGS Beutler et al. (2011), SDDS Anderson et al. (2014), BOSS CMASS Xu et al. (2012), WiggleZ Blake et al. (2012), MGS Ross et al. (2015), BOSS DR12 Gil-Marin et al. (2016) and DES Y1 Abbott et al. (2017b). For the CMB we use the shift parameters $(R,l_{a})$ based on the Planck 2018 release Aghanim et al. (2018) and as derived by Ref. Zhai and Wang (2018). We assume the existence of three families of neutrinos with $N_{\textrm{eff}}=3.046$.

Furthermore, we also incorporate the direct measurements of the Hubble expansion $H(z)$ data. These can be derived in two ways: by the clustering of galaxies or quasars and by the differential age method. The former provides direct measurements of the Hubble parameter by measuring the BAO peak in the radial direction from the clustering of galaxies or quasars Gaztanaga et al. (2009). The latter method obtains the Hubble parameter via the redshift drift of distant objects over significant time periods, usually a decade or longer. This is possible as in GR the Hubble parameter can be expressed via the rate of change of the redshift $H(z)=-\frac{1}{1+z}\frac{dz}{dt}$ Jimenez and Loeb (2002). These methods result in a compilation of 36 Hubble parameter $H(z)$ data points, which for clarity we show in Table 1 along with their corresponding references.

The growth-rate data used here are obtained via the redshift-space distortions (RSD). These are sensitive probes of the Large Scale Structure (LSS) and can measure the quantity $f\sigma_{8}(a)\equiv f(a)\cdot\sigma(a)$, which is a product of the growth rate $f(a)=\frac{dln\delta}{dlna}$ and the redshift-dependent rms fluctuations $\sigma(a)=\sigma_{8,0}\frac{\delta(a)}{\delta(1)}$ of the linear density field within spheres of radius $R=8h^{-1}\textrm{{Mpc}}$. In this notation the parameter $\sigma_{8,0}$ is the value of the rms fluctuations today and is a direct measure of the amplitude of fluctuations in linear scales.

We should mention that $f\sigma_{8}(a)$ can be estimated via the ratio of the monopole to the quadrupole of the redshift-space power spectrum $P(k)$. The latter is sensitive on the quantity $\beta=f/b_{1}$, where $f$ is the growth-rate as defined earlier and $b_{1}$ is the galaxy bias Percival and White (2009); Song and Percival (2009); Nesseris and Perivolaropoulos (2007). In all cases we assume linear theory. The combination $f\sigma_{8}(a)$ not only is independent of bias, as the latter completely cancels out, but it has also been demonstrated to be an excellent discriminator of DE models as it probes the dynamics of a given gravitational theory and not only the geometric of space-time Song and Percival (2009). The covariances of the data and how to make the necessary corrections for the Alcock-Paczynski effect are given in Refs. Sagredo et al. (2018); Nesseris et al. (2017); Kazantzidis and Perivolaropoulos (2018), while other related analyses with these data can be found in Refs. Basilakos et al. (2018); Basilakos and Nesseris (2017, 2016); Arjona et al. (2019).

In this paper we use the growth-rate data compilation of Ref. Sagredo et al. (2018), which we show in Table 2 for completeness, along with the corresponding references for each point. This dataset was analyzed in Ref. Sagredo et al. (2018) with the “Internal Robustness method” of Ref. Amendola et al. (2013), by examining combinations of subsets and it was shown that this specific dataset is indeed internally robust.

With these in mind, our total likelihood function $L_{\rm tot}$ can be given as the product of the separate likelihoods of the data (we assume they are statistically independent) as follows:

[TABLE]

which is related to the total $\chi^{2}$ via $\chi^{2}_{\rm tot}=-2\log{L_{\rm tot}}$ or

[TABLE]

Calculating the best-fit is not enough, but we also need to study the statistical significance of our constraints. To achieve this we make use of the well known Akaike Information Criterion (AIC) Akaike (1974). The AIC estimator is given (assuming Gaussian errors) by

[TABLE]

where $k_{p}$ and $N_{\rm dat}$ stand for the number of free parameters and the total number of data points respectively. For other similar statistical tools see also Ref. Liddle (2007). In this analysis we have 1048 data points from the Pantheon set, 3 from the CMB shift parameters, 10 from the BAO measurements, 22 from the growth measurements and finally 36 $H(z)$ points, for a total of $N_{\rm dat}=1118$.

The AIC can be interpreted similarly to the $\chi^{2}$, i.e. a smaller relative value signifies a better fit to the data. To apply this statistic to model selection we take the pair difference between models $\Delta{\rm AIC}={\rm AIC}_{\rm model}-{\rm AIC}_{\rm min}$. This can in principle be interpreted with the Jeffreys’ scale in the following manner: when $4<\Delta{\rm AIC}<7$ this indicates positive evidence against the model with higher value of ${\rm AIC}_{\rm model}$, while in the case when $\Delta{\rm AIC}\geq 10$ it can be interpreted as strong evidence. On the other hand, when $\Delta{\rm AIC}\leq 2$, then this means that the two models are statistically equivalent. However, in Ref. Nesseris and Garcia-Bellido (2013) it has been shown that in general the Jeffreys’ scale can sometimes lead to misleading conclusions, and thus it should be interpreted with care.

Finally, our total $\chi^{2}$ is given by Eq. (190) while the parameter vectors (assuming a spatially flat Universe) are given by: $p_{\Lambda\textrm{CDM}}=\left(\Omega_{m,0},100\Omega_{b}h^{2},h,\sigma_{8}\right)$ for the $\Lambda$CDM and $p_{\textrm{HDES}}=\left(\Omega_{m,0},100\Omega_{b}h^{2},h,\tilde{J_{c}},\sigma_{8}\right)$ for the HDES model. Using the aforementioned cosmological data and methodology, we can obtain the best-fit parameters and their uncertainties via the MCMC method based on a Metropolis-Hastings algorithm. The codes used in the analysis were written by one of the authors.888The MCMC code for Mathematica used in the analysis is freely available at http://members.ift.uam-csic.es/savvas.nesseris/. The priors we assumed for the parameters are given by $\Omega_{m,0}\in[0.1,0.5]$, $\Omega_{b}h^{2}\in[0.001,0.08]$, $\tilde{J_{c}}\in[-1,12]$, $h\in[0.4,1]$, $\sigma_{8}\in[0,2]$ and we sample $\sim 10^{5}$ MCMC points for each of the two models.

VII.1 Results

In Figs. 7 and 8 we show the 68.3$\%$, 95.4$\%$ and 99.7$\%$ confidence contours for the $\Lambda$CDM and the HDES models, respectively, along with the one-dimensional (1D) marginalized likelihoods for all parameter combinations in the familiar triangle plot. We also highlight with a black point the mean MCMC values and with a red point the Planck 2018 concordance cosmology. The latter is based on the TT,TE,EE+lowP spectra, a flat $\Lambda$CDM model and the values are shown in Table 3.

In Tables 4 and 5 we show the best-fit values of the model parameters and the values for the $\chi^{2}$ and AIC parameters for the $\Lambda$CDM and the HDES model respectively. As can be seen from Tables 4 and 5, we find that as the difference in the AIC parameters is roughly $\sim 0.68$, then both models seem to be statistically equivalent with each other. Furthermore, as seen in Fig. 8, there is a clear negative correlation between $\tilde{J_{c}}$ and $\sigma_{8}$ as we saw in Sec. V.5 and Eq. (183) due to the strong degeneracy between the parameters. This degeneracy is useful as it can potentially alleviate and relax the tension that has been recently observed, see Refs. Nesseris et al. (2017); Sagredo et al. (2018). In particular, in Fig 9 we show the 68.3$\%$, 95.4$\%$ and 99.7$\%$ confidence contours for the $\Lambda$CDM (left) and the HDES $(n=2)$ (right) models respectively in the $(\Omega_{m,0},\sigma_{8})$ plane. As can be seen, for the HDES model, the best-fit in the $(\Omega_{m,0},\sigma_{8})$ plane moves toward higher values of $\sigma_{8}$, closer to those of Planck.

VIII Conclusions

The recent discovery of gravitational waves emission from a binary neutron star merger with an optical counterpart, signified a major breakthrough in astrophysics and cosmology as it provided a direct measurement of the speed of propagation of gravitational waves. This observation not only represented an important advance for astronomy, but it also served to greatly reduce the number of alternative models aiming at explaining the current accelerating phase of the Universe. In particular, since the constraint on the speed of propagation of gravitational waves is extremely close to the speed of light, the Horndeski Lagrangian simplified to only three functions. Although this means a notable progress in constraining cosmological models, degeneracies with the $\Lambda$CDM model remain and must be further investigated.

In this paper we used an effective fluid approach to study the remaining Horndeski Lagrangian. This formalism makes it possible to compare models with different underlying physics (e.g., DE and MG models) in a relatively easy way: each model is mapped to three functions describing the effective fluid, namely, the equation of state $w$, the sound speed $c_{s}^{2}$, and the anisotropic stress $\pi$. Even though the remaining Horndeski Lagrangian is now simpler than its original version, finding exact analytical solutions can be quite laborious. Nevertheless, the subhorizon and quasistatic approximations are pretty helpful at overcoming this difficulty.

One of our main results is the set of Eqs. (102)-(IV.1). These equations along with the equation of state Eq. (73) describe the remaining Horndeski Lagrangian in an effective fluid approach under the subhorizon and quasistatic approximations. In this paper, we provide explicit expressions for the effective fluid description of several DE and MG models.

In order to exemplify our results and since we focused on explanations to the late-time accelerating universe, we carried out an analysis where only DM and an effective DE fluid are taken into consideration. A particularly interesting model also included in our formalism is the KGB model. In Sec. VI we show our analytical solutions agree pretty well with a full numerical solution of the system of differential equations describing the DM and effective DE perturbations. We also confirm that the subhorizon approximation breaks down for the KGB model due to the rapid oscillations of the scalar field in the large $n$ limit, in agreement with Ref. Kimura and Yamamoto (2011). Also, for the KGB model the background equation for the expansion history $H(a)$ can only be found numerically for $n>1$, thus slowing down the codes significantly.

Due to these problems, we propose a completely new class of Horndeski models based on the designing principle, i.e., fixing the background to a specific model, usually that of the $\Lambda$CDM and then determining the Lagrangian. Given the freedom in specifying the remaining functions of the Horndeski Lagrangian, we propose a way to find families of models which match a particular background expansion, i.e., the equation of state $w_{DE}$. Since current observations are in good agreement with the standard $\Lambda$CDM at the background level, we provide equations specifying a $w_{DE}=-1$ designer Horndeski model (see Eqs. (169)), which we call HDES. Furthermore, for this model we are able to find exact solutions for the growth $\delta_{m}(a)$ in the matter domination epoch by solving Eq. (101). The solutions we found are given by Eq. (179) and they imply a degeneracy between $\sigma_{8}$ and the parameter of the HDES model $\tilde{J_{c}}$, which can approximately be described via Eq. (183).

Although fixing the background to $\Lambda$CDM is a common practice, the treatment of the perturbations might not be rigorous enough in current studies. Public codes solving the perturbation equations for the Horndeski Lagrangian (e.g., hi$\_$CLASS) use ad hoc parametrizations for the $\alpha_{i}$ functions which differ significantly from our findings that approximate a realistic model (see Eqs. (171)-(175)), see for example Refs. Perenon et al. (2019); Noller and Nicola (2018b); Spurio Mancini et al. (2019).

We implemented the parametrized version for the DE effective fluid of our $w_{DE}=-1$ designer Horndeski HDES model in the public code CLASS, which we call EFCLASS, by following the straightforward implementation explained in our previous paper Arjona et al. (2019). For the sake of comparison and in order to check the validity of our effective fluid approach, we compared results from our code EFCLASS with the public code hi_CLASS, which solves numerically the full perturbation equations.

In Fig. 5 we show the CMB angular power spectrum computed with both codes and as can be seen in the right panel of Fig. 5, the agreement is remarkable and on average on the order of $\sim 0.1\%$. Since the hi$\_$CLASS code does not utilize either the subhorizon or the quasistatic approximation, but our EFCLASS does it, we conclude our effective fluid approach is quite accurate and powerful. Furthermore, the main advantage of our method is that while hi_CLASS requires significant and non-trivial modifications, our EFCLASS code practically only requires the implementation of Eq. (187), which is trivial.

We further investigated our $w_{DE}=-1$ designer Horndeski HDES model by computing cosmological constraints with recent data sets using an MCMC analysis. The results of our MCMC analysis are shown in Tables 4 and 5, where we present the best-fit values of the model parameters and the values for the $\chi^{2}$ and AIC parameters for the $\Lambda$CDM and the HDES model respectively. We find that as the difference in the AIC parameters is roughly $\sim 0.68$, then both models seem to be statistically equivalent with each other. Furthermore, as seen in Fig. 8, there is a clear negative correlation between $\tilde{J_{c}}$ and $\sigma_{8}$. This can be understood, as we saw in Sec. V.5, due to the strong degeneracy between the parameters described by Eq. (183). This degeneracy is useful as it can potentially alleviate the $\sigma_{8}$ tension that has been recently observed, see Ref. Nesseris et al. (2017); Sagredo et al. (2018).

Numerical Analysis Files: The numerical codes used by the authors in the analysis of the paper and our modifications to the CLASS code, which we call EFCLASS, will be released upon publication of the paper on the websites of the EFCLASS here and here.

Acknowledgements

The authors would like to thank Hector Gil Marín and Tomás Ortín for useful discussions. They also acknowledge support from the Research Projects FPA2015-68048-03-3P [MINECO-FEDER], PGC2018-094773-B-C32 and the Centro de Excelencia Severo Ochoa Program SEV-2016-0597. S.N. also acknowledges support from the Ramón y Cajal program through Grant No. RYC-2014-15843.

Appendix A Scalar and Gravitational field equations

For completeness, in this Appendix we show how to compute both the gravitational and the scalar-field equations derived from the Horndeski action (21).

A.1 Scalar field equation

For a function of a single variable with higher derivatives, the stationary values of the functional Courant and Hilbert (1953)

[TABLE]

can be obtained from the Euler-Lagrange equation

[TABLE]

Since our Lagrangian $\mathcal{L}_{i}$ functions defined in the Horndeski action (21) depend on the scalar field $\phi$ and its first and second derivatives, we can use the Euler-Lagrange equation (193) to compute the scalar field equation for $\mathcal{L}_{2}$, $\mathcal{L}_{3}$ and $\mathcal{L}_{4}$. For $\mathcal{L}_{2}$ we have

[TABLE]

[TABLE]

Since $\mathcal{L}_{2}=K\left(\phi,X\right)$, applying Eq. (195) leads to

[TABLE]

[TABLE]

where we have replaced the partial derivatives by covariant derivatives and we are using the fact that $X=-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi$. Hence, for $\mathcal{L}_{2}$ the scalar field equation reads

[TABLE]

For the term $\mathcal{L}_{3}$ we follow the same approach

[TABLE]

[TABLE]

Knowing that $\mathcal{L}_{3}=-G_{3}\left(\phi,X\right)\left[\Box\phi=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\phi\right]$, applying Eq. (200) gives

[TABLE]

[TABLE]

[TABLE]

where we have replaced again the partial derivatives by covariant derivatives. We can then conclude that, for $\mathcal{L}_{3}$ the scalar field equation reads

[TABLE]

and we make the following assignment

[TABLE]

For $\mathcal{L}_{4}$ we have

[TABLE]

[TABLE]

Since $\mathcal{L}_{4}=G_{4}\left(\phi\right)R$, applying Eq. (208) leads to

[TABLE]

Our result for the scalar field equation considering $G_{4X}=0$ and $G_{5}=0$ is in full agreement with Ref. Kobayashi et al. (2011). Hence, the scalar-field equation can be written as

[TABLE]

A.2 Gravitational field equations

Defining the arbitrary functions $\mathcal{L}_{i}$ from the action (21) as

[TABLE]

we can then vary the action with respect to the metric tensor; using the principle of least action, this leads to

[TABLE]

For $\delta S_{2}$ we have

[TABLE]

and using the fact that

[TABLE]

and that the variation of $K$ with respect to the metric can be written as

[TABLE]

we get

[TABLE]

For $\delta S_{3}$ we have

[TABLE]

The variations of $G_{3}$ with respect to the metric can be written as

[TABLE]

hence

[TABLE]

The last term of the above equation can be expanded in the following way

[TABLE]

since

[TABLE]

and

[TABLE]

Also we have that $g^{ab}\Gamma^{\gamma}_{ab}=\ldots=-\nabla_{a}\delta g^{\gamma a}+\frac{1}{2}g_{ab}g^{\gamma\lambda}\nabla_{\lambda}\delta g^{ab}$, so we get for the last term in Eq. (221):

[TABLE]

Combining all terms we have

[TABLE]

For $\delta S_{4}$ we have

[TABLE]

where

[TABLE]

Since $\delta\Gamma^{\lambda}_{\mu\nu}$ is the difference of two connections, it should transform as a tensor. Therefore, it can be written as

[TABLE]

Then, substituting Eq. (229) into (228), we get

[TABLE]

hence

[TABLE]

where

[TABLE]

[TABLE]

Since the energy-momentum tensor is defined as

[TABLE]

the gravitational field equation can be written

[TABLE]

Appendix B Coefficients

Here we show the coefficients for the perturbations in the Horndeski theory in Eq. (21). They are given by:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and using Eq.(68) to eliminate $G_{4}$ in favor of $K$ we can express $B_{9}$ as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and using Eqs. (68) and (III.1) we find that

[TABLE]

[TABLE]

For the DE effective perturbation equations we found the following coefficients

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The coefficients for the KGB DE effective perturbation equations are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]