Do cosmological data rule out f(R) with w≠−1 ?
DDo cosmological data rule out f ( R ) with w (cid:54) = − ? Richard A. Battye, ∗ Boris Bolliet, † and Francesco Pace ‡ Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy,The University of Manchester, Manchester, M13 9PL, U.K. (Dated: June 18, 2018)We review the Equation of State (EoS) approach to dark sector perturbations and apply it to f ( R ) gravitymodels of dark energy. We show that the EoS approach is numerically stable and use it to set observationalconstraints on designer models. Within the EoS approach we build an analytical understanding of the dynamicsof cosmological perturbations for the designer class of f ( R ) gravity models, characterised by the parameter B and the background equation of state of dark energy w . When we use the Planck cosmic microwave backgroundtemperature anisotropy, polarisation and lensing data as well as the baryonic acoustic oscillation data from SDSSand WiggleZ, we find B < . (95% C.L.) for the designer models with w = − . Furthermore, we find B < . and | w + 1 | < . (95% C.L.) for the designer models with w (cid:54) = − . Previous analysesfound similar results for designer and Hu-Sawicki f ( R ) gravity models using the Effective Field Theory (EFT)approach [Raveri et al. , Phys. Rev. D , 043513 (2014); Hu et al. , Mon. Not. R. Astron. Soc. , 3880(2016)]; therefore this hints for the fact that generic f ( R ) models with w (cid:54) = − can be tightly constrained bycurrent cosmological data, complementary to solar system tests [Brax et al. , Phys. Rev. D , 104021 (2008);Faulkner et al. , Phys. Rev. D , 063505 (2007)]. When compared to a w CDM fluid with the same sound speed,we find that the equation of state for f ( R ) models is better constrained to be close to -1 by about an order ofmagnitude, due to the strong dependence of the perturbations on w . PACS numbers: 04.50.Kd, 95.36.+x, 98.80.-kKeywords: Cosmology; modified gravity; dark energy; f(R) gravity
I. INTRODUCTION
With the observational campaign of Supernovae type Ia [1–4], followed by observations of the Cosmic Microwave Back-ground (CMB) anisotropy [5, 6], the Baryon Acoustic Oscil-lations (BAO) [7, 8] and large scale structure [9–11], it hasbecome widely accepted that the expansion of the universeis accelerating. The current observational data is consistentwith the standard Λ cold dark matter (CDM) model, wherethe accelerated expansion is caused by the cosmological con-stant Λ , and indicates no statistically significant evidence fordark energy and modified gravity models (see, e.g., [12] andreferences therein).Nevertheless, the cosmological constant suffers from im-portant conceptual issues when it is interpreted in the contextof quantum field theory (see, e.g., [13] for a recent review).This has led part of the community to question the physicalorigin of the accelerated expansion and to investigate dark en-ergy and modified gravity models (see, e.g., [14]). Whetherthese models do not suffer the same type of issues as the cos-mological constant often remains under debate.Moreover, the forthcoming galaxy surveys and stage IVCMB experiments will measure the acceleration of the uni-verse and its consequences on structure formation at a levelof accuracy never achieved before. Hence, research on darkenergy and modified gravity is well motivated by the follow-ing question: In the light of this forthcoming data, will thecosmological constant still be the best answer to cosmic ac- ∗ [email protected] † [email protected] ‡ [email protected] celeration? In other words, is there a modified gravity or darkenergy model that will account for the observational data ina better way than the cosmological constant? Of course, thishas to be formulated in a precise statistical manner, see [15]for an example in the context of inflationary models.Recently, the Horndeski models [16–18] have received agrowing attention due to their generality. They include ascalar field coupled to gravity. The Horndeski Lagrangianis the most general one that leads to second order equationsof motion for the scalar field. It is fully represented by fourarbitrary time dependent functions of the scalar field and itskinetic term. Notable subclasses of the Horndeski models, ob-tained by specifying the unknown functions, are Quintessence[19–26], k -essence [27–32], Brans-Dicke theory [33, 34], Ki-netic Gravity Braiding (KGB) [35, 36] and f ( R ) models [37–40]. The latter can also be constructed by replacing the Ricciscalar R in the Einstein-Hilbert Lagrangian by an arbitraryfunction, f ( R ) , and are the main focus of this paper.Here, we are interested in f ( R ) models that mimic the Λ CDM (or the w CDM) cosmological expansion history butdiffer at the level of the dynamics of cosmological perturba-tions. Different approaches have been developed to study thephenomenology of cosmological perturbations in dark energyand modified gravity in a unified way, with the ultimate ob-jective of deriving observational constraints. These includethe Parameterized Post Friedmaniann (PPF) approach [41–44], the Equation of State for perturbations (EoS) approach[45–47] (see also [48] for an earlier and similar approach),the Effective Field Theory (EFT) approach [49–53] and [54]for an alternative method. They are in principle equivalent(see [55] for a numerical consistency analysis), although theydiffer with respect to the choices of the phenomenologicalparametrisation of dark energy and modified gravity. So far, a r X i v : . [ a s t r o - ph . C O ] J un the EFT approach has been applied to generic Horndeski mod-els [51, 52], while the EoS approach has been applied specifi-cally to quintessence, k -essence and KGB models [47], f ( R ) gravity [56] and Generalised Einstein-Aether theories [57]. Inthis paper we use the EoS approach, for which the dark en-ergy and modified gravity models are specified in terms of theanisotropic stress and pressure of the perturbed dark energyfluid.The paper is organised as follows. In section II we re-view the EoS approach and its numerical implementation in aBoltzmann code for arbitrary dark energy and modified grav-ity models. In Sec. III we recall the features of the designer f ( R ) models that are relevant to our analysis. In Sec. IV westudy the phenomenology of cosmological perturbations prop-agating in the dark energy fluid of the models with constant w de , numerically as well as analytically. In Sec. V we presentthe linear matter power spectrum, the CMB temperature angu-lar anisotropy power spectrum and the CMB power spectrumof the lensing potential, computed for several designer modelsand we derive observational constraints on the free parametersof the designer models, i.e., w de and B , from current CMBand BAO data. In Sec. VI we compare f ( R ) and w CDM grav-ity models and their observational constraints. We discuss ourresults and conclude in Sec. VII. In Appendix we present acomparison between the perturbed equations of state obtainedwithin the EoS [58] and EFT approaches [51, 52].Unless otherwise stated, we use π G = 1 , throughout thepaper. II. NUMERICAL IMPLEMENTATION OF THEEQUATION OF STATE APPROACH
In the EoS approach, modifications to general relativity arewritten in the right hand side of the field equations. Then,they can be interpreted as a stress energy tensor, mapping anymodified gravity theory to a corresponding dark energy fluid.More precisely, we have G µν = T µν + D µν , (1)where G µν is the Einstein tensor, T µν is the stress energytensor of the matter components, i.e., baryonic matter, radi-ation and dark matter, and D µν is the stress-energy tensorof the dark energy fluid. The background geometry is as-sumed to be isotropic and spatially flat, with a line element ds = − dt + a δ ij dx i dx j , where a is the scale factor and t is the cosmic time. Due to the Bianchi identities and the localconservation of energy for the matter components, the stressenergy tensor of the dark sector is covariantly conserved, ∇ µ D µν = 0 . (2)The linear perturbation of the conservation equations (2)yields the general relativistic version of the Euler and continu-ity equations for the velocity and density perturbation. Theycharacterise the dynamics of cosmological perturbations andcan be written in terms of a gauge invariant density perturba-tion, ∆ , and a rescaled velocity perturbation, Θ . These two quantities are defined as ∆ ≡ δ + 3(1 + w ) Hθ , Θ ≡ w ) Hθ , (3)where w ≡ P/ρ is the background equation of state, ρ and P are the homogeneous density and pressure, δρ is the densityperturbation, θ is the divergence of the velocity perturbation,and H ≡ ( d ln a/dt ) is the Hubble parameter.The rescaled velocity perturbation, Θ , is not a gauge invari-ant quantity, in the sense that its value depends on the choiceof the coordinate system, see, e.g. [59]. To see this, say that Θ is evaluated in the conformal Newtonian gauge (CNG), i.e., Θ c = Θ , where the superscript c indicates the CNG. Then thevalue of the rescaled velocity perturbation in the synchronousgauge (SG), Θ s , is given, in Fourier space, by Θ s = Θ c − w ) T , (4)with T ≡ (cid:40) ( h (cid:48) + 6 η (cid:48) ) / (2K ) in the SG , . (5)where K ≡ k/ ( aH ) and k is the wavenumber of the pertur-bation, h and η are the scalar metric perturbations in the SG,and where a prime denotes a derivative with respect to ln a .Since the SG is defined as the rest frame of the CDM fluid,we see that T is nothing else than the velocity perturbation ofthe CDM fluid evaluated in the CNG.To work in a gauge invariant way, with respect to the syn-chronous and conformal Newtonian gauges, we can define agauge invariant velocity perturbation as ˆΘ ≡ Θ + 3(1 + w ) T , (6)In the same line of thought, using the variable T , the evolu-tion equations for the gauge invariant density perturbation andrescaled velocity perturbation can be written in a way that isvalid for both gauges [56]. These are the so-called perturbedfluid equations and are given by ∆ (cid:48) − w ∆ −
2Π + g K (cid:15) H ˆΘ = 3(1 + w ) X , ˆΘ (cid:48) + 3 (cid:0) c a − w + (cid:15) H (cid:1) ˆΘ − c a ∆ − −
3Γ = 3(1 + w ) Y , (7)where c a ≡ dP/dρ is the adiabatic sound speed and g K ≡ / (3 (cid:15) H ) , with (cid:15) H ≡ − H (cid:48) /H and where X ≡ (cid:40) η (cid:48) + (cid:15) H T in the SG ,φ (cid:48) + ψ in the CNG , (8a) Y ≡ (cid:40) T (cid:48) + (cid:15) H T in the SG ,ψ in the CNG . (8b)Finally, Π is the perturbed scalar anisotropic stress and Γ is the gauge invariant entropy perturbation . The gauge in-variant entropy perturbation can be expressed in terms of the Note that our θ and Π differ from θ MB and σ MB (anisotropic stress) as de-fined in [59], by θ MB = k a θ and ( ρ + P ) σ MB = − ρ Π . perturbed pressure, density and rescaled velocity as Γ = δPρ − c a (∆ − Θ) . (9)The perturbed fluid equations (7) are valid for both matter(that we shall denote with a subscript ‘m’) and dark energy(that we shall denote with subscript ‘de’) fluid variables.The Einstein-Boltzmann code CLASS [60, 61] written in Cprovides the infrastructure required to solve the dynamics ofmatter perturbations. We have incorporated the EoS approachfor dark energy perturbations into
CLASS and dubbed the modified code
CLASS_EOS_FR . The code is publicly avail-able on the internet . We have implemented the perturbedfluid equations (7) for dark energy perturbations in this exactsame form. We now describe the remaining technical stepsnecessary to close the system of equation (7) and integrate itin the code.As prescribed by the EoS approach, we expand the per-turbed dark energy anisotropic stress and gauge invariant en-tropy perturbation in terms of the perturbed fluid variables.These are the so-called equations of state for dark energy per-turbations and are written as Π de = c Π∆de ∆ de + c ΠΘde ˆΘ de + c Π∆m ∆ m + c ΠΘm ˆΘ m + c ΠΠm Π m , Γ de = c Γ∆de ∆ de + c ΓΘde ˆΘ de + c Γ∆m ∆ m + c ΓΘm ˆΘ m + c ΓΓm Γ m , (10)where the coefficients c αβ are a priori scale and time depen-dent functions, but shall only depend on the homogeneousbackground quantities, such as the Hubble parameter, thebackground equation of state of dark energy, or the adiabaticsound speeds. These functions are specified for each dark en-ergy and modified gravity model, e.g., see [56] for f ( R ) grav-ity and [57] for Generalised Einstein-Aether. Note that theequations of state for perturbations for generic f ( R ) modelscan also be obtained starting from a general Horndeski modeland specifying the appropriate free functions to match with f ( R ) theories. In this case, the expressions for the coeffi-cients of c αβ are as reported in appendix .Initial conditions for dark sector perturbations are set atan early time, a ini , when dark energy is subdominant, i.e., Ω de ( a ini ) (cid:28) where Ω de is the dark energy density parameter.If not specified from the specific dark energy model, appro-priate initial conditions for the dark energy perturbations aregenerally: ∆ de ( a ini ) = Θ de ( a ini ) = 0 . Note that when there ex-ists an attractor for the dark energy perturbations during mat-ter domination, it is numerically more efficient to set initialconditions that match the attractor (see Sec. IV).In order to evaluate the equation of state (10) and integrateequations (7), we collect the perturbed matter fluid variablesat every time step. In our code, we do this in the followingway. First, we obtain the total matter gauge invariant densityperturbation via Ω m ∆ m = − K Z − Ω de ∆ de with Z ≡ (cid:40) η − T in the SG φ in the CNG (11)and the gauge invariant matter velocity perturbation via Ω m ˆΘ m = 2 X − Ω de ˆΘ de , see [56] where these equations arederived. Next, the matter pressure perturbation δP m and the website:https://github.com/borisbolliet/class eos fr public matter anisotropic stress σ classm are available in CLASS . We usethem to compute the matter anisotropic stress perturbation (inour convention) Π m and the matter gauge invariant entropyperturbation as ρ m Π m = − (cid:104) ( ρ m + P m ) σ classm (cid:105) , (12) ρ m Γ m = (cid:104) δP m (cid:105) − c a, m (∆ m − Θ m ) , (13)where the brackets mean a sum over all the matter fluid com-ponents, i.e., baryons, CDM, photons and neutrinos, and c a, m = w m Ω m + (cid:10) w m Ω m (cid:11) (1 + w m ) Ω m , (14)is the matter adiabatic sound speed, where Ω m ≡ − Ω de and w m ≡ (cid:104) w m Ω m (cid:105) / Ω m are the matter density parameter and back-ground equation of state respectively. Last, we update the totalstress energy tensor accordingly as δρ tot = (cid:104) δρ m (cid:105) + ρ de ∆ de − ρ de Θ de ( ρ tot + P tot ) θ classtot = (cid:104) ( ρ m + P m ) θ classm (cid:105) + K aHρ de Θ de ( ρ tot + P tot ) σ classtot = (cid:104) ( ρ m + P m ) σ classm (cid:105) − ρ de Π de δP tot = (cid:104) δP m (cid:105) + ρ de Γ de + c a, de ρ de (∆ de − Θ de ) . See footnote 1 for the
CLASS perturbed velocity, which fol-lows the conventions of [59].Although the numerical integration can be carried out ei-ther in the conformal Newtonian gauge or in the synchronousgauge in
CLASS_EOS_FR , we find that, in the super-Hubbleregime, i.e., K (cid:28) , the synchronous gauge performs betterthan the conformal Newtonian gauge. III. A BRIEF REMINDER ON THE DESIGNER f ( R ) GRAVITY MODELS In f ( R ) gravity, the f ( R ) functions are solutions to a sec-ond order differential equation given by the projection of the − − z − − − − B B = 10 B = 1 B = 0 . B = 0 . − − z − − − − | B | B = 0 . w de = − . w de = − . w de = − . w de = − . w de = − . FIG. 1. The redshift evolution of B = − ( f (cid:48)R / [ (cid:15) H (1 + f R )]) fordifferent designer f ( R ) models. Unless otherwise written, we chose w de = − and B = 1 . A grey line indicates negative values. Thebackground cosmology was set to h = 0 . , Ω de = 0 . and Ω b h =0 . , where h = H / is the reduced Hubble parameter. stress-energy tensor of f ( R ) on the time direction, which canbe written as [37, 62–65] f (cid:48)(cid:48) + (cid:18) (cid:15) H − − ¯ (cid:15) (cid:48) H ¯ (cid:15) H (cid:19) f (cid:48) − ¯ (cid:15) H f = 6 H ¯ (cid:15) H Ω de , (15)where the prime still denotes a derivative with respect to ln a and ¯ (cid:15) H = (cid:15) (cid:48) H + 4 (cid:15) H − (cid:15) H (see Eq. (2.6a) of [56] for the deriva-tion in our conventions). This equation holds for any f ( R ) gravity model and at any time during the expansion history.During the non-relativistic matter era, i.e., w m = 0 , this equa-tion simplifies because (cid:15) H = 3 / , ¯ (cid:15) (cid:48) H = 0 and ¯ (cid:15) H = (cid:15) H = 3 / (see Eq. (2.5) of [56]). In this regime, the solutions to (15) are f ( a ) = C (cid:110) b + a n + + b − a n − + e − (cid:82) w de )d ln a (cid:111) , (16)with n ± = ( − ± (cid:112) / and C = H w de +5 w de − . So-lutions with b − (cid:54) = 0 are not admissible because they breakthe condition lim a → f R = 0 [66–68], where a subscript ‘ R ’means a derivative with respect to the Ricci scalar. We con-clude that any viable f ( R ) gravity model can be parameter-ized, in the non-relativistic matter era, by the a priori timedependent equation of state w de ( a ) and a constant number b + . We then trade b + for the more commonly used parameter B ≡ − f (cid:48)R (cid:15) H (1 + f R ) , (17)evaluated today and dubbed B , since there is a one-to-onecorrespondence between b + and B . From here, there are twoways to proceed. The first possibility is to specify explicitlya f ( R ) function at all time, and then extract the time evolu-tion of Ω de and w de from the time derivatives of f . The secondpossibility is to specify a time evolution for Ω de and w de andthen integrate Eq. (15) to get f ( R ) at all time. This latter ap-proach is the so-called designer, or mimetic, f ( R ) approachand leads to the f ( R ) gravity models that we are interestedin. Designer models are particularly interesting because theirfunctional form is dictated by the chosen background evolu-tion of the dark fluid and therefore there is no arbitrariness inhow the f ( R ) Lagrangian looks like. In this way the wantedbackground evolution is achieved exactly and the model hasless degrees of freedom: the only value to be determined is B ,which ultimately will dictate the strength of the perturbations.In CLASS_EOS_FR , we have implemented the designermodels with constant equation of state w de . The user speci-fies a value for w de and B , then the code explores a range of b + solving (15), between a ini and today, until it finds the valueof b + that leads to the desired value of B . Note that the so-lution in (16) is singular for w de (cid:39) . and w de (cid:39) − . ,however as long as one avoids the two poles, the numericalintegration is efficient.In [37, 56], the designer models with w de = − were stud-ied at both the background and perturbation levels. Here, weconsider as well the designer models with w de (cid:54) = − (and w (cid:48) de = 0 ), i.e., the ones that mimic a w CDM expansion his-tory.In Fig. 1 we show the redshift evolution of a set of solutionsto (15) for different values of B and w de . We present B , ratherthan f ( R ) itself, because this is the main quantity entering theequations of state for perturbation Π de and Γ de [56]. On thebottom panel we fix B = 1 and vary w de . For models with w de < − , B starts being negative and eventually becomespositive at late time. This can be described analytically withEq. (16), see, e.g., [37]. On the top panel we fix w de = − andvary B . As can be seen, as soon as dark energy dominates,i.e., z (cid:46) . , B settles to its final value B . Changing thevalue of w de essentially amounts to a shift of the curves onthis plot because for a less negative w de dark energy dominatesearlier. The bottom panel shows that when we keep B fixed, B grows more slowly for less negative w de . More precisely,with Eq. (16) in the matter era, one finds B ∼ z w de . IV. EVOLUTION OF PERTURBATIONS IN THE DARKENERGY FLUID OF f ( R ) GRAVITY
In this section we investigate numerically and analyticallythe evolution of cosmological perturbations for the designer f ( R ) gravity models described in Sec. III. To this aim, weuse the formalism of the EoS approach described in Sec. II.To gain some understanding about the behaviour of the cos-mological perturbations, we consider the expressions of theequations of state for perturbations for a f ( R ) fluid with con-stant equation of state parameter, i.e., c a, de = w de , and whenthe matter sector is dominated by non-relativistic species, i.e., w m = Π m = Γ m = 0 , as is the case after radiation domination.Furthermore, we focus on modes that enter the Hubble hori-zon before dark energy dominates so that we have K (cid:29) at all time. This assumption holds for wavenumbers in theobservational range of interest to us (see top panel of Fig. 2).Finally we assume B (cid:28) , which is true at all times if B (cid:28) and is equivalent to M (cid:29) , with M ≡ (cid:15) H / ( (cid:15) H B ) . In thisregime, the equations of state for dark energy perturbationssimplify to Π de = ∆ de , (18a) Γ de = (cid:110) − w de + M K (cid:111) ∆ de +
13 Ω m Ω de ∆ m . (18b)Using the field equation (3.11a) and (3.11b) in [56], theperturbed fluid equations (7) can be rewritten as a system oftwo coupled second order differential equations for the gaugeinvariant density perturbations, ∆ (cid:48)(cid:48) m + (2 − (cid:15) H )∆ (cid:48) m − Ω m ∆ m = − Ω de ∆ de , (19a) ∆ (cid:48)(cid:48) de + (2 − (cid:15) H )∆ (cid:48) de + (K + M )∆ de = −
13 Ω m Ω de K ∆ m . (19b)For the modes of interest, this set of equations provides a faith-ful description of the dynamics of cosmological perturbationsas long as B (cid:28) . Again, this is always the case before darkenergy dominates (irrespective of B ). In addition if B (cid:28) ,then these equations are also valid during dark energy domi-nation, because B is always smaller than B (see Fig. 1). Letus assume B = O (1) , or equivalently M (cid:29) , from nowon. As we shall see in Sec. V, this is a reasonable assumptiongiven current observational constraints.The differential equation (19b) for the gauge invariant en-ergy density perturbation is similar to an harmonic oscillatorwith a time dependent frequency ω = K + M (cid:29) . Sincethe oscillatory time scale is much smaller than the dampingtime scale. i.e., the expansion rate, the homogeneous solu-tion to (19b) becomes rapidly subdominant compared to theparticular solution. This confirms that the specific values forthe initial dark energy perturbations are not important. Moreprecisely, the dark energy density perturbation relates to thematter density perturbation via Ω de ∆ de = −
13 K K + M Ω m ∆ m . (20)We refer to [69] for the same result formulated in a differentlanguage. In our code, we set the initial conditions for ∆ de and Θ de according to (20) at a time such that K / [3(K + M )] = | Ω de ∆ de / Ω m ∆ m | = 0 . . Note that given this criterion, theinitial starting time for dark energy perturbation depends onthe wavenumber.We deduce from (20) the two regimes for the behaviour ofsub-horizon modes: (i) the general relativistic (GR) regimewhen K (cid:28) M , i.e., at early time, and (ii) the scalar-tensor(ST) regime when K (cid:29) M , i.e., at late time. This im-plies Ω de ∆ de = − K M Ω m ∆ m in the GR regime, and Ω de ∆ de = − − z K , M k = 1Mpc − k = 0 . − k = 0 . − M − − z Ω d e ∆ d e k = 1Mpc − k = 0 . − k = 0 . − .
00 0 .
25 0 .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . z − . − . − . − . − . . . . . . γ ≡ l n ( ∆ m / ∆ m ) / l n Ω m γ ΛCDM = 6 / γ ST k = 0 .
01 Mpc − k = 0 .
10 Mpc − k = 1 .
00 Mpc − (1 + w de ) × . . . . . . σ log B = − . B = − . B = − . B = − . FIG. 2. The redshift evolution of K for three wavenumbers and M (dashed line in the top panel), Ω de ∆ de and γ (middle panels) and σ asa function of w de (bottom panel) for different designer f ( R ) models.The attractor solution (20) and the growth index γ ST (25) are the thickgrey lines. Unless otherwise written, we chose w de = − and B =0 . as well as the same cosmology as in Fig. 1 with A s = 2 . × − and n s = 0 . . − Ω m ∆ m in the ST regime. Moreover, in both regimes,the differential equation for the matter perturbation (19a) be-comes ∆ (cid:48)(cid:48) m + (2 − (cid:15) H )∆ (cid:48) m − ε Ω m ∆ m = 0 , (21)where ε ≡ (4K + 3M ) / (3K + 3M ) can be interpretedas a modification to the gravitational constant [70]. One has ε = 4 / in the ST regime and ε = 1 in the GR regime. Sinceone has K ∼ z − and M ∼ z − w de during the matter era,the ST regime starts earlier for less negative w de .Eq. (20) and (21) enable a clear discussion of the dynamicsof cosmological perturbations in f ( R ) gravity. Before doingso, we go one step further and obtain the growth index γ ≡ ln f / ln Ω m of the matter perturbation [71], where f ≡ ∆ (cid:48) m / ∆ m is the growth rate.Taking the time derivative of the growth rate and using (21)we find γ (cid:48) + 3 w de Ω de ln Ω m γ + Ω γ m ln Ω m − − γ m m ε = 3 w de Ω de −
12 ln Ω m , (22)for the growth index. To linearise this equation, we use theapproximations ln Ω m ≈ − Ω de and Ω γ m ≈ − γ Ω de which arevalid when Ω de = O (1) . We get γ (cid:48) + (cid:0) − w de + ε (cid:1) γ = (cid:16) − ε Ω de + ε − w de (cid:17) . (23)This can be solved analytically for a constant ε . We find γ = 3(1 − ε )2 + 3 ε Ω m,0 Ω de,0 (1 + z ) − w de + 3( ε − w de )2 + 3 ε − w de . (24)In the GR regime the first term on the right hand side vanishes,and the second term gives a constant γ w CDM = 3(1 − w de ) / (5 − w de ) , i.e., the w CDM growth index. If in addition w de = − ,the the growth index is γ ΛCDM = 6 / ≈ . , i.e., the well-known Λ CDM result. In the ST regime, the growth index isnot constant any more due to the first term on the right handside. We find γ ST = 12 + 16(1 − w de ) − Ω m,0 de,0 (1 + z ) − w de . (25)Since the first term on the right hand side of (24) is alwaysnegative, we have γ ST < γ w CDM as well as γ ST < γ ΛCDM .We now summarise the important consequences for the dy-namics of perturbations in f ( R ) gravity that are deduced fromthe above considerations.1. For B < (or M < ), the homogeneous solution to(19b) is unstable. Therefore, the gauge invariant densityperturbation for both matter and dark energy grows ex-ponentially with time. This is not compatible with thedynamics of matter perturbations in the matter domi-nated era, and consequently f ( R ) models with B < or w de < − are not viable, see Fig. 1.2. The gauge invariant density perturbation in the dark en-ergy component relates to that of the matter componentin a simple way given in (20). In the GR regime, thedark energy perturbation is negligible compared to thematter perturbation, while in the ST regime both are ofthe same magnitude, see Fig. 2. 3. In w CDM, for less negative w de structures are less grav-itationally bounded compared to Λ CDM because darkenergy starts dominating earlier. Hence there is an anti-correlation between w de and the amplitude of clustering,i.e. σ , in w CDM models (see, e.g. Fig. 16 of [72]). In f ( R ) gravity matter perturbations grow at a faster ratethan in w CDM and Λ CDM because γ ST < γ w CDM [seeEq. (25)]. This, combined with the fact that the STregime starts earlier for less negative w de , implies a cor-relation between w de and σ (see bottom panel of Fig. 2),and can be used to discriminate between f ( R ) gravityand w CDM models of dark energy [see also the nextsection for a comparison between w CDM and f ( R ) models].In the next section we compute relevant observables that weuse to set observational constraints on the designer f ( R ) grav-ity models. V. IMPACT OF f ( R ) GRAVITY ON OBSERVABLES ANDCONSTRAINTS
The CMB angular anisotropy power spectrum is a snapshotof the acoustic waves in the photon-baryon fluid at decou-pling, distorted by the Integrated Sachs-Wolfe effect (ISW)and the lensing due to the subsequent gravitational collapseof the matter. How and when can dark energy perturbationsin f ( R ) gravity affect the CMB anisotropy? Since in viable f ( R ) gravity models, dark energy perturbations are subdom-inant at early time (see point 2 on page 6), they can not haveany impact on the physical phenomena at play at the epoch ofdecoupling. However, they alter the growth of structure fromthe end of the matter dominated era (see point 3 on page 6).Therefore, they may have an impact on the late ISW effect(see, e.g. [73]) and lensing of the CMB anisotropy (see, e.g.[74]). The late ISW effect is contributing to the CMB tem-perature anisotropy on large angular scales ( (cid:96) (cid:46) ) and thelensing power spectrum of the CMB probes structure forma-tion on a wider range of scales ( (cid:96) (cid:46) ). So we expect theCMB angular anisotropy power spectrum to be affected bydark energy perturbations only at low multipoles, i.e., wherethe cosmic variance limits the constraining power of the CMBdata. Hence, the lensing power spectrum shall be a more com-pelling probe of dark energy perturbation than the CMB tem-perature anisotropy angular power spectrum.In the left panels of Fig. 3 we show the CMB temperatureangular anisotropy power spectrum computed for several de-signer models with different w de and B , against the Λ CDMprediction. We see that significant differences appear when B (cid:38) and that the late ISW effect can be strongly enhancedfor larger values of B . Moreover, at fixed B the late ISWcontribution is more significant for less negative w de , as can beunderstood with the results of Sec. IV (see point 3 on page 6).In the middle panels we show the CMB lensing power spec-trum computed in the same settings. Its amplitude is largerfor larger B and less negative w de , again in agreement with ‘ − − ‘ ( ‘ + ) C TT ‘ / π ΛCDM B = 5 B = 1 B = 0 . ‘ − − − [ ‘ ( ‘ + ) ] C φφ ‘ / π ΛCDM B = 5 B = 1 B = 0 . − − − k [ h/ Mpc] P ( k ) [ M p c / h ] ΛCDM B = 5 B = 1 B = 0 . ‘ − − ‘ ( ‘ + ) C TT ‘ / π ΛCDM w de = − . w de = − . w de = − . ‘ − − − [ ‘ ( ‘ + ) ] C φφ ‘ / π ΛCDM w de = − . w de = − . w de = − . − − − k [ h/ Mpc] P ( k ) [ M p c / h ] ΛCDM w de = − . w de = − . w de = − . FIG. 3. Effects of f R gravity on the CMB angular temperature power spectrum (left), lensing power spectrum (middle) and the linear matterpower spectrum (right) for different designer f ( R ) models against the Λ CDM predictions. Unless otherwise written, we chose w de = − and B = 0 . as well as the same cosmology as in Fig. 2.TABLE I. Posterior mean (68% C.L.) for log B , σ and w de for designer f ( R ) models that mimic a Λ CDM and a w CDM expansion. Theellipses indicate the absence of 68% C.L. constraints, in this case only the 95% C.L. upper limits are relevant (see Table II).CMB+BAO CMB+BAO+Lensing CMB+BAO CMB+BAO+Lensing( Λ CDM) ( Λ CDM) ( w CDM) ( w CDM) log B − . +1 . − . · · · · · · · · · σ . +0 . − . · · · . +0 . − . . +0 . − . (1 + w de ) × . +1 . − . . +0 . − . TABLE II. Posterior upper limits (95% C.L.) for log B , σ and w de for designer f ( R ) models that mimic a Λ CDM and a w CDM expansion.CMB+BAO CMB+BAO+Lensing CMB+BAO CMB+BAO+Lensing( Λ CDM) ( Λ CDM) ( w CDM) ( w CDM) log B < − . < − . < − . < − . σ < . < . < . < . w de ) × < < . the analysis of Sec. IV. Similar conclusions apply to the linearmatter power spectrum presented in the right panels of Fig. 3.In particular, for scales which are still in the GR regime to-day ( k ≈ − h Mpc − ), the amplitude of the matter powerspectrum is close to the Λ CDM prediction, while for scalesthat entered the ST regime during the matter dominated era( k (cid:38) − h Mpc − ), its amplitude is enhanced. For observational constraints, we consider the follow-ing combinations of data sets: CMB+BAO and CMB,BAO+Lensing. For CMB and Lensing we refer to the Planck2015 public likelihoods for low- (cid:96) and high- (cid:96) temperature aswell as polarisation and lensing data [5]. For BAO we referto the distance measurements provided by the WiggleZ DarkEnergy Survey [75] and SDSS [76]. We use Montepython [77] for the Monte Carlo Markov chain sampling of the pa-rameter space. We varied the six base cosmological parame-
TABLE III. Posterior mean (68% C.L.) for σ and w de for a w CDMmodel. CMB+BAO CMB+BAO+Lensing σ . +0 . − . . +0 . − . (1 + w de ) × − . +7 . − . − . +6 . − . TABLE IV. Posterior upper limits (95% C.L.) for σ and w de for a w CDM model. CMB+BAO CMB+BAO+Lensing σ < . < . w de ) × < . < ters as well as all the Planck nuisance parameters. For those,we used the same priors as the Planck Collaboration [5]. Inaddition we varied the background dark energy equation ofstate w de and log B that characterise the designer f ( R ) mod-els. For w de we used a uniform prior between − and . For log B we used a uniform prior between − and . In Tables Iand II, we show the 68% C.L. and 95% C.L. constraints fromour analyses.For designer models with w de = − , B and σ are deter-mined at 68% C.L. for CMB+BAO. We get B ≈ . and σ (cid:39) . ± . . If we add the information relative to clus-tering at late time, via the CMB lensing data, B and σ arenot determined, but constrained to B (cid:46) . and σ < . (95% C.L.).For designer models with w de (cid:54) = − , B is not determinedany more by CMB+BAO. Moreover, due to the correlation be-tween w de and the amplitude of clustering (see bottom panelof Fig. 2), σ takes substantially larger values than with the w de = − models. When we add CMB lensing data, the pos-terior mean value of σ is brought down by fifteen percentand more importantly the 68% C.L. region for the dark en-ergy background equation of state is reduced by a factor often. We get (1 + w de ) < . , in other words the expansionhistory has to be very close to Λ CDM.
VI. COMPARISON WITH w CDM MODELS
To quantify the relative importance of perturbations in (de-signer) f ( R ) models, we can compare their observationalconstraints with a w CDM model where the background equa-tion of state w is free to vary (but constant in time) and wekeep the sound speed (defined in the frame comoving with thefluid) c = δp/δρ = 1 fixed. To study the perturbations ofsuch a model, we use the CLASS implementation of the pa-rameterized post-Friedmaniann (PPF) framework as describedin [78]. When w de ≥ − , this framework recovers the be-haviour of canonical minimally coupled scalar field modelsand it is accurate also when w de ≈ − . A welcome aspectof the PPF formalism is that it allows to study the evolution of perturbations in the phantom regime ( w de < − ), which isusually preferred by Supernovae data [3, 79]. In addition, thecrossing of the “phantom barrier” ( w de = − ) is allowed, cov-ering therefore also the more general case of non-canonicalminimally coupled models, such as k -essence. The PPF for-malism allows also sound speeds c (cid:54) = 1 , as in k -essence mod-els, but here we limit ourselves to the standard case of luminalsound speed, as this is also the value in f ( R ) models.We note, in principle, that in w CDM models w de can takevalues smaller than -1, this is the so-called phantom regime,while in the designer f ( R ) models we consider in this workthe phantom crossing is not allowed due to instabilities, seeSec. IV.Moreover we saw that in f ( R ) gravity small variations of w de lead to large variations in σ (see bottom panel of figure 2),while in w CDM models small variations of w de lead to smallvariations in σ : in the range of w de presented in the bottompanel of figure 2, for the same cosmological parameters, σ would vary by less than 1%.Using the same data sets described before, in Tables III andIV we show the 68% and 95% C.L. constraints on σ and w de for the w CDM fluid, respectively. For w de , we use a uniformprior between − and .Our results agree with [5]. In particular, the preferred valuefor w de is in the phantom regime. It means that these datasets favour a higher value of σ with respect to the Λ CDMcosmology, as was the case for the f ( R ) models.Our last remark is that since σ depends weakly on w de in w CDM compared to f ( R ) , the constraints on w de in w CDMare weaker than in f ( R ) by one order of magnitude, see ta-bles III and IV. VII. DISCUSSION AND CONCLUSION
Intense observational and theoretical efforts are being de-ployed to unveil the nature of the cosmic acceleration of theuniverse. Going beyond the cosmological constant Λ , twomain hypotheses can be explored: dark energy and modifiedgravity. Many models belonging to these two broad groupscan be described in terms of the Horndeski Lagrangian. In thiswork we concentrated on a well studied sub-class of Horn-deski theories, the so-called f ( R ) gravity models. Such mod-ifications to GR may affect both the background expansionhistory and the evolution of cosmological perturbations. Inthis paper we considered the designer f ( R ) gravity modelsfor which the f ( R ) function is tuned to reproduce the w CDMexpansion history.We used the EoS approach to study analytically the dynam-ics of linear cosmological perturbations in this context, andwe implemented it numerically in our
CLASS_EOS_FR code.To prove the reliability of our numerical implementation, wecompared our results with several other f ( R ) codes publiclyavailable such as MGCAMB [80, 81],
FRCAMB [82],
EFTCAMB [83–85] and found agreement at the sub-percent level for all ofthem [55], except for
MGCAMB which disagreed by more thanfive percent relative error with the other codes for the compu-tation of the matter power spectrum for k > h Mpc − .Unlike for the simple w CDM dark energy model, we foundthat for designer f ( R ) gravity models a less negative w de leadsto a larger σ (see point 3 on page 6). To arrive at this conclu-sion we derived an analytical formula for the growth index γ (see Eq. (25)).Using CMB lensing data we found that designer f ( R ) models with (1 + w de ) > . and B > . are dis-favoured at 95% C.L. Note that similar constraints wereobtained for the designer f ( R ) models also by [84], usingcosmological data as we did here. The authors of [86]performed a similar analysis on the Hu-Sawicki f ( R ) modelsand found, as we did, a higher value of σ with respect to the Λ CDM value . Moreover, for the screening mechanism tohappen on solar system scales the authors of [88, 89] found | w de | (cid:46) − for generic f ( R ) models.The results we obtained are consistent with these previousanalyses and hint for the fact that generic f ( R ) models with w de (cid:54) = − can be ruled out based on current cosmologicaldata, complementary to solar system tests. VIII. ACKNOWLEDGEMENTS
Appendix: Comparison between the EoS and EFT approachesfor dark energy perturbations
In this section we compare the expressions for the entropyperturbations and the perturbed anisotropic stress of [56] withthe corresponding expressions from [52] in the conformalNewtonian gauge. In the following we will denote with thesuperscript “BBP” variables in [56] and with“GLV” variablesin [52]. In addition we use the subscript ‘m’ for all matterspecies and ζ i = g K (cid:15) H − ¯ (cid:15) H g K (cid:15) H − dP i dρ i . (A.1)In f ( R ) gravity the equations of states for scalar perturba-tions, in both formalisms are [52, 56] Π BBPde = K g K (cid:15) H (cid:26) ∆ de − f (cid:48)R f R ) Θ de + Ω m Ω de f R f R ∆ m − Ω m Ω de f (cid:48)R f R ) Θ m (cid:27) − f R f R Ω m Ω de Π m , (A.2a) Γ BBPde = (cid:26) ζ de − ¯ (cid:15) H g K (cid:15) H f R ) − f (cid:48)R f (cid:48)R (cid:27) ∆ de − ζ de Θ de + Ω m Ω de (cid:26) ζ m − ¯ (cid:15) H g K (cid:15) H f R − f (cid:48)R f (cid:48)R (cid:27) ∆ m − Ω m Ω de ζ m Θ m − Ω m Ω de Γ m , (A.2b) P GLVde Γ GLVde = γ γ + γ α B K γ + α B K ( δρ de − Hq de ) + γ γ + α B K γ + α B K H ( q de + q m ) + 13 ( δρ m − Hq m ) − dP de dρ de δρ de − δp m , (A.2c) P GLVde Π GLVde = γ α B K γ + α B K ) ( δρ de − Hq de ) − γ K γ + α B K ) H ( q de + q m ) , (A.2d)where the functions γ i are given by γ = 3 α B (cid:15) H , γ = , γ = − ¯ (cid:15) H (cid:15) H α B , γ = 1 − ¯ (cid:15) H (cid:15) H , γ = − , γ = − α B . Wefurther define α B = f (cid:48)R f R ) . Note that with respect to [52],we defined P de Π GLVde = − k a σ GLVde .Unlike [58], the authors of [52] use a non-standard con-tinuity equation for the effective dark energy fluid which implies ρ GLVde = ρ BBPde + 3 M pl H f R ,P GLVde = P BBPde − M pl H (3 − (cid:15) H ) f R , for the background and δρ GLVde = (1 + f R ) δρ BBPde + f R δρ BBPm ,δP
GLVde = (1 + f R ) δP BBPde + f R δP BBPm ,q GLVm + q GLVde = − f R H { ρ BBPde Θ BBPde + ρ BBPm Θ BBPm } ,q GLVde = − H { (1 + f R ) ρ BBPde Θ BBPde + f R ρ BBPm Θ BBPm } ,P GLVde Π GLVde = [(1 + f R ) P BBPde Π BBPde + f R P BBPm Π BBPm ] . for the perturbed fluid variables. From this, we conclude thatboth formalisms are equivalent.0 [1] A. G. Riess, A. V. Filippenko, P. Challis, and et al., AJ ,1009 (1998), arXiv:astro-ph/9805201.[2] S. Perlmutter, G. Aldering, G. Goldhaber, and et al., Astrophys.J. , 565 (1999), arXiv:astro-ph/9812133.[3] A. G. Riess, L.-G. Strolger, J. Tonry, S. Casertano, H. C. Fer-guson, and et al., Astrophys. J. , 665 (2004), arXiv:astro-ph/0402512.[4] A. G. Riess, L.-G. Strolger, S. Casertano, H. C. Ferguson,B. Mobasher, B. Gold, P. J. Challis, and et al., Astrophys. J. , 98 (2007), arXiv:astro-ph/0611572.[5] Planck Collaboration XIII, A&A , A13 (2016),arXiv:1502.01589.[6] G. Hinshaw, D. Larson, E. Komatsu, D. N. Spergel, C. L. Ben-nett, J. Dunkley, M. R. Nolta, M. Halpern, R. S. Hill, and et al.,ApJS , 19 (2013), arXiv:1212.5226 [astro-ph.CO].[7] W. J. Percival, B. A. Reid, D. J. Eisenstein, and et al., MNRAS , 2148 (2010), arXiv:0907.1660 [astro-ph.CO].[8] D. Parkinson, S. Riemer-Sørensen, C. Blake, G. B. Poole, T. M.Davis, S. Brough, M. Colless, C. Contreras, W. Couch, andet al., Phys. Rev. D , 103518 (2012), arXiv:1210.2130 [astro-ph.CO].[9] S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek,A. S. Bolton, J. R. Brownstein, A. Burden, and et al., MNRAS , 2617 (2017), arXiv:1607.03155.[10] S. Rota, B. R. Granett, J. Bel, L. Guzzo, J. A. Peacock, M. J.Wilson, A. Pezzotta, S. de la Torre, B. Garilli, and et al., A&A , A144 (2017), arXiv:1611.07044.[11] M. Ata, F. Baumgarten, J. Bautista, F. Beutler, D. Bizyaev,M. R. Blanton, J. A. Blazek, A. S. Bolton, J. Brinkmann, andet al., MNRAS , 4773 (2018), arXiv:1705.06373.[12] Planck Collaboration XIV, A&A , A14 (2016),arXiv:1502.01590.[13] D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata,A. G. Riess, and E. Rozo, Physics Reports , 87 (2013),arXiv:1201.2434.[14] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, PhysicsReports , 1 (2012), arXiv:1106.2476.[15] T. Giannantonio and E. Komatsu, Phys. Rev. D , 023506(2015), arXiv:1407.4291.[16] G. W. Horndeski, International Journal of Theoretical Physics , 363 (1974).[17] C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade, Phys. Rev.D , 064039 (2011), arXiv:1103.3260 [hep-th].[18] T. Kobayashi, M. Yamaguchi, and J. Yokoyama, Progress ofTheoretical Physics , 511 (2011), arXiv:1105.5723 [hep-th].[19] L. H. Ford, Phys. Rev. D , 2339 (1987).[20] P. J. E. Peebles and B. Ratra, ApJL , L17 (1988).[21] B. Ratra and P. J. E. Peebles, Phys. Rev. D , 3406 (1988).[22] C. Wetterich, Nuclear Physics B , 668 (1988).[23] R. R. Caldwell, R. Dave, and P. J. Steinhardt, Physical ReviewLetters , 1582 (1998), arXiv:astro-ph/9708069.[24] E. J. Copeland, A. R. Liddle, and D. Wands, Phys. Rev. D ,4686 (1998), gr-qc/9711068.[25] P. J. Steinhardt, L. Wang, and I. Zlatev, Phys. Rev. D ,123504 (1999), arXiv:astro-ph/9812313. Using CFHTLenS data [87], the normalisation of the matter power spec-trum is significantly closer to the Λ CDM value, implying a lower value of f R , ( B in our notation) and, using their Eq. 10 (see also their figure 2), w de ≈ − . [26] T. Barreiro, E. J. Copeland, and N. J. Nunes, Phys. Rev. D ,127301 (2000), arXiv:astro-ph/9910214.[27] C. Armend´ariz-Pic´on, T. Damour, and V. Mukhanov, PhysicsLetters B , 209 (1999), hep-th/9904075.[28] T. Chiba, T. Okabe, and M. Yamaguchi, Phys. Rev. D ,023511 (2000), astro-ph/9912463.[29] N. A. Hamed, H. S. Cheng, M. A. Luty, and S. Mukohyama,Journal of High Energy Physics , 074 (2004), hep-th/0312099.[30] F. Piazza and S. Tsujikawa, JCAP , 004 (2004), hep-th/0405054.[31] R. J. Scherrer, Physical Review Letters , 011301 (2004),astro-ph/0402316.[32] V. Mukhanov and A. Vikman, JCAP , 004 (2006), astro-ph/0512066.[33] C. Brans and R. H. Dicke, Physical Review , 925 (1961).[34] A. De Felice and S. Tsujikawa, JCAP , 024 (2010),arXiv:1005.0868 [astro-ph.CO].[35] C. Deffayet, O. Pujol`as, I. Sawicki, and A. Vikman, JCAP ,026 (2010), arXiv:1008.0048 [hep-th].[36] O. Pujol`as, I. Sawicki, and A. Vikman, Journal of High EnergyPhysics , 156 (2011), arXiv:1103.5360 [hep-th].[37] Y.-S. Song, W. Hu, and I. Sawicki, Phys. Rev. D , 044004(2007), astro-ph/0610532.[38] A. Silvestri and M. Trodden, Reports on Progress in Physics ,096901 (2009), arXiv:0904.0024 [astro-ph.CO].[39] T. P. Sotiriou and V. Faraoni, Reviews of Modern Physics ,451 (2010), arXiv:0805.1726 [gr-qc].[40] A. De Felice and S. Tsujikawa, Living Reviews in Relativity , 3 (2010), arXiv:1002.4928 [gr-qc].[41] C. Skordis, Phys. Rev. D , 123527 (2009), arXiv:0806.1238[gr-qc].[42] T. Baker, P. G. Ferreira, C. Skordis, and J. Zuntz, Phys. Rev. D , 124018 (2011), arXiv:1107.0491 [astro-ph.CO].[43] T. Baker, P. G. Ferreira, and C. Skordis, Phys. Rev. D ,024015 (2013), arXiv:1209.2117 [astro-ph.CO].[44] P. G. Ferreira, T. Baker, and C. Skordis, General Relativity andGravitation , 1788 (2014).[45] R. A. Battye and J. A. Pearson, JCAP , 019 (2012),arXiv:1203.0398 [hep-th].[46] R. A. Battye and J. A. Pearson, Phys. Rev. D , 061301 (2013),arXiv:1306.1175.[47] R. A. Battye and J. A. Pearson, JCAP , 051 (2014),arXiv:1311.6737.[48] M. Kunz and D. Sapone, Physical Review Letters , 121301(2007), astro-ph/0612452.[49] J. Bloomfield, ´E. ´E. Flanagan, M. Park, and S. Watson, JCAP , 010 (2013), arXiv:1211.7054 [astro-ph.CO].[50] J. Bloomfield, JCAP , 044 (2013), arXiv:1304.6712 [astro-ph.CO].[51] J. Gleyzes, D. Langlois, F. Piazza, and F. Vernizzi, JCAP ,025 (2013), arXiv:1304.4840 [hep-th].[52] J. Gleyzes, D. Langlois, and F. Vernizzi, International Jour-nal of Modern Physics D , 1443010 (2014), arXiv:1411.3712[hep-th].[53] J. Gleyzes, D. Langlois, M. Mancarella, and F. Vernizzi, JCAP , 054 (2015), arXiv:1504.05481.[54] D. Bertacca, N. Bartolo, and S. Matarrese, JCAP , 021 (2012),arXiv:1109.2082.[55] E. Bellini, A. Barreira, N. Frusciante, B. Hu, S. Peirone,M. Raveri, M. Zumalac´arregui, A. Avilez-Lopez, M. Ballardini,and et al., Phys. Rev. D , 023520 (2018), arXiv:1709.09135. [56] R. A. Battye, B. Bolliet, and J. A. Pearson, Phys. Rev. D ,044026 (2016), arXiv:1508.04569.[57] R. A. Battye, F. Pace, and D. Trinh, Phys. Rev. D , 064041(2017), arXiv:1707.06508.[58] R. A. Battye and F. Pace, Phys. Rev. D , 063513 (2016),arXiv:1607.01720.[59] C.-P. Ma and E. Bertschinger, Astrophys. J. , 7 (1995),astro-ph/9506072.[60] J. Lesgourgues, ArXiv e-prints (2011), arXiv:1104.2932 [astro-ph.IM].[61] D. Blas, J. Lesgourgues, and T. Tram, JCAP , 034 (2011),arXiv:1104.2933.[62] L. Pogosian and A. Silvestri, Phys. Rev. D , 023503 (2008),arXiv:0709.0296.[63] S. Nojiri, S. D. Odintsov, and D. S´aez-G´omez, Physics LettersB , 74 (2009), arXiv:0908.1269 [hep-th].[64] P. K. S. Dunsby, E. Elizalde, R. Goswami, S. Odintsov,and D. Saez-Gomez, Phys. Rev. D , 023519 (2010),arXiv:1005.2205 [gr-qc].[65] L. Lombriser, A. Slosar, U. Seljak, and W. Hu, Phys. Rev. D , 124038 (2012), arXiv:1003.3009.[66] J. Khoury and A. Weltman, Phys. Rev. D , 044026 (2004),arXiv:astro-ph/0309411.[67] J. Khoury and A. Weltman, Physical Review Letters , 171104(2004), astro-ph/0309300.[68] W. Hu and I. Sawicki, Phys. Rev. D , 064004 (2007),arXiv:0705.1158.[69] S. Tsujikawa, R. Gannouji, B. Moraes, and D. Polarski, Phys.Rev. D , 084044 (2009), arXiv:0908.2669 [astro-ph.CO].[70] A. A. Starobinsky, Soviet Journal of Experimental and Theoret-ical Physics Letters , 157 (2007), arXiv:0706.2041.[71] P. J. E. Peebles, The large-scale structure of the universe (Princeton University Press, Princeton, 1980).[72] E. Komatsu, J. Dunkley, M. R. Nolta, and et al., ApJS , 330(2009), arXiv:0803.0547.[73] J. Lesgourgues,
Proceedings, Theoretical Advanced Study Insti-tute in Elementary Particle Physics: Searching for New Physics at Small and Large Scales (TASI 2012): Boulder, Colorado,June 4-29, 2012 , ArXiv e-prints (2013), arXiv:1302.4640[astro-ph.CO].[74] A. Lewis and A. Challinor, Physics Reports , 1 (2006),astro-ph/0601594.[75] E. A. Kazin, J. Koda, C. Blake, N. Padmanabhan, S. Brough,M. Colless, C. Contreras, W. Couch, S. Croom, and et al., MN-RAS , 3524 (2014), arXiv:1401.0358.[76] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden,and M. Manera, MNRAS , 835 (2015), arXiv:1409.3242.[77] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP , 001 (2013), arXiv:1210.7183.[78] W. Fang, W. Hu, and A. Lewis, Phys. Rev. D , 087303(2008), arXiv:0808.3125.[79] J. L. Tonry, B. P. Schmidt, B. Barris, P. Candia, P. Challis,A. Clocchiatti, A. L. Coil, A. V. Filippenko, P. Garnavich,C. Hogan, and et al., Astrophys. J. , 1 (2003), arXiv:astro-ph/0305008.[80] G.-B. Zhao, L. Pogosian, A. Silvestri, and J. Zylberberg, Phys.Rev. D , 083513 (2009), arXiv:0809.3791.[81] A. Hojjati, L. Pogosian, and G.-B. Zhao, JCAP , 005 (2011),arXiv:1106.4543 [astro-ph.CO].[82] J.-h. He, Phys. Rev. D , 103505 (2012), arXiv:1207.4898.[83] B. Hu, M. Raveri, N. Frusciante, and A. Silvestri, Phys. Rev. D , 103530 (2014), arXiv:1312.5742.[84] M. Raveri, B. Hu, N. Frusciante, and A. Silvestri, Phys. Rev. D , 043513 (2014), arXiv:1405.1022.[85] B. Hu, M. Raveri, N. Frusciante, and A. Silvestri, ArXiv e-prints (2014), arXiv:1405.3590 [astro-ph.IM].[86] B. Hu, M. Raveri, M. Rizzato, and A. Silvestri, MNRAS ,3880 (2016), arXiv:1601.07536.[87] C. Heymans, E. Grocutt, A. Heavens, M. Kilbinger, T. D. Kitch-ing, F. Simpson, J. Benjamin, T. Erben, H. Hildebrandt, andet al., MNRAS , 2433 (2013), arXiv:1303.1808.[88] P. Brax, C. van de Bruck, A.-C. Davis, and D. J. Shaw, Phys.Rev. D , 104021 (2008), arXiv:0806.3415.[89] T. Faulkner, M. Tegmark, E. F. Bunn, and Y. Mao, Phys. Rev.D76