Sweeping Horndeski canvas: new growth-rate parameterization for modified-gravity theories

We propose and numerically validate a new fitting formula that is sufficiently accurate to model the growth of structure in Horndeski theories of modified gravity for upcoming Stage IV and V large-scale structure surveys. Based on an analysis of more than 18,000 Horndeski models and adopting the popular parameterization of the growth rate f(z) = Ω M (z) γ , we generalize the constant growth index γ to a two-parameter redshift-dependent quantity, γ(z), that more accurately fits these models. We demonstrate that the functional form γ(z) = γ 0 + γ 0 z 2/(1 + z) improves the median χ 2 of the fit to viable Horndeski models by a factor of ∼ 40 relative to that of a constant γ, and is sufficient to obtain unbiased results even for precise measurements expected in Stage IV and V surveys. Finally, we constrain the parameters of the new fitting formula using current cosmological data.


Introduction
Over ∼13.7 billions years of cosmic evolution, the tiny primordial fluctuations seeded during inflation -under gravitational interaction -evolve into the large-scale structure observed and measured by galaxy surveys today.The temporal growth of cosmic structure has a rich and well-understood behavior in different epochs: it is robust in the matter-dominated era, but suppressed at late times, especially following the onset of dark energy.The clustering of galaxies, the weak gravitational lensing of distant background galaxies, and arguably the abundance of galaxy clusters as a function of their redshifts and mass proxies have all established themselves as powerful probes of structure growth.These measurements then translate into constraints on models of dark matter and dark energy or modified gravity (see, e.g.[1] for a recent, general review).
In the linear regime (corresponding to scales k ≲ 0.1h Mpc −1 today), the growth of density fluctuations is described by the linear growth function D(a) ≡ δ(a)/δ(a = 1).From it, we can define the growth rate as where a is the cosmic scale factor.For the standard, smooth dark-energy with an equation of state w, on sub-horizon-scales and in the absence of massive neutrinos, D and f are scale-independent 1 .The growth rate can be formally obtained by solving the second-order differential equation which, in standard gravity and on sub-horizon scales, reads D + 2H Ḋ = 4πGρ M D.Here ρ M is the background matter density, G is the Newton constant, and dots are derivatives with respect to time.
In a wide class of cosmological models, the growth rate is well-approximated by a fitting function 2f (z) = Ω M (z) γ(z) . (1.2) Here, Ω M (z) is the time-dependent matter density relative to critical, and the free function γ(z) is the so-called growth index.The latter has a long history in cosmology, dating back to [2][3][4][5], as it describes the growth rate in standard matter-dominated cosmologies.Ref. [6] proposed and verified that the growth-index parameterization γ = 0.55 + 0.02(1 + w(z = 1)) fits the true growth for all wCDM models to better than 0.2%, all the way from the matter-dominated era to the present time, for a wide range of Ω M values.Moreover, this parameterization also provides a good fit to the growth in some modified-gravity models with, for example, γ ≃ 0.67 for the DGP models [7,8].
A fitting formula for the growth rate such as Eq.(1.2) is useful for at least two reasons.First, it is easy to implement in cosmological analyses.Second, it is straightforward to test whether the growth agrees with the prediction of a cosmological model (say γ ≃ 0.55 for ΛCDM [9] (see also, e.g.[10,11]).Therefore, there has been considerable interest in developing phenomenological formulae for the growth rate and, in particular, investigating their robustness with respect to the choice of the cosmological model.Analytic works, e.g.[7,[12][13][14], have examined the best-fit γ(z) = γ 0 -the redshift-independent values of γ that minimized deviations of Eq. (1.2) from the exact growths defined in Eq. (1.1) -in various modified gravity and dark energy models.Further, [15][16][17][18][19][20][21][22][23][24] focused on the redshift evolution of the growth index in f (R) gravity, f (Q) gravity, and interacting dark energy models.The most common redshift-dependent growth-index descriptions are the linear parameterization γ(z) = γ 0 + γ 1 z and the Taylor expansion γ(z) = γ 0 + γ 1 z/(1 + z), the latter of which is motivated by the parameterization w(a) = w 0 + w a (1 − a) for the time-dependent equation of state of dark energy w(a) [25,26].
Our principal goal in this paper is to explore the accuracy of different γ(z) parameterizations within the viable space of Horndeski theory, in the context of future constraints on f (z)σ 8 (z) ≡ f σ 8 (z) by Stage-IV and Stage-V large-scale structure surveys.These upcoming surveys will yield f σ 8 (z) measurements with small error bars and extend to high redshifts.
To that effect, we adopt exact numerical calculations of Eq. (1.1) as the baseline, and then fit the growth rate f σ 8 (z) in ∼18,000 Horndeski models using a broad set of functional forms.We then compare their goodness of fit to the ground truth in Horndeski models in the context of errors predicted for future surveys, and propose the best new parameterization of γ(z).We further demonstrate the utility of the proposed parameterization by constraining its parameters using current observational data of type Ia supernovae, large-scale structure, and the cosmic microwave background (CMB), and briefly comment on the implications for stress-testing the standard cosmological model.
The rest of this paper is structured as follows.Sec. 2 reviews the effective field theory framework we exploit to evaluate the growth rate f σ 8 (z) in a given Horndeski model.Sec. 3 details our procedure to sample Horndeski models and to obtain each model's theoretical prediction on f σ 8 (z).Sec. 4 discusses various parameterizations of γ(z) focusing on their performance in fitting the theory models, and presents current constraints on parameters of the best fitting formula.Finally, Sec. 5 summarizes our analysis.

Horndeski models: theory background and growth of structure
We wish to consider the most general class of ΛCDM extensions for the accelerating universe that is not strongly disfavored by current data.We therefore must find an effective way to sample the model space.To this end, we adopt the Effective Field Theory (EFT) formalism for dark energy and modified gravity (henceforth EFTDE) [27,28].Within EFTDE, models with similar properties are established through a grouping of terms in the fundamental Lagrangian such that one can consider a class of models together (for more details, see [29] and references therein).

Effective Field Theory approach to Dark Energy
In general, the EFTDE action in unitary gauge can be written as, e.g.[30,31] where δg 00 = g 00 + 1 is the perturbation to the time component of the metric, R (3) is the perturbation to the spatial component, and δK µν is the perturbation of the extrinsic curvature.The background evolution depends on three EFTDE functions, c(t), Λ(t), and Ω(t).For any given expansion history, the first two functions, c(t) and Λ(t) can be constrained by the Friedmann equations and correspond to energy density and pressure.The effect of modified gravity is parameterized by the third function Ω(t).The other EFTDE functions in Eq. (2.1) represent perturbations around the background and correspond to observables that can be compared with observations.Tab. 1 in [32] gives a summary of all models that can be represented by the EFT formalism.
Within EFTDE, we focus on the Horndeski class of models (see, e.g.[33] for an in-depth review).The class of Horndeski theories is the most general scalar-tensor extension of general relativity, including but not limited to quintessence (see [34] for a review) and generalized Brans-Dicke (Jordan Brans-Dicke [35], f(R) [36], chameleons [37]) models.Moreover, within these scalar-tensor theories, a coupling between the derivative of a scalar field and the Einstein tensor (or the Ricci tensor alone) leads to an accelerated expansion of the cosmic background without demanding a scalar potential [38,39].
The Horndeski class is specified by imposing additional constraints on the EFTDE functions that describe perturbations around the background, as follows: To evaluate the growth rate in a given Horndeski model and cosmology, we employ the EFTCAMB framework3 [30,31,40].EFTCAMB characterizes a given Horndeski model by seven functions which we parameterize as follows.One aforementioned function, Ω(t), controls the background evolution.We henceforth relabel it Ω MG , for it not to be confused with an energy density parameter.
Inspired by f (R) gravity and (again) following the convention in [30,31,40], we further assume Ω MG (t) evolves in time as Ω MG (a) = Ω MG,0 a s 0 . (2.3) For ΛCDM, Ω MG (a) = 0. Further, there are six dimensionless, second-order EFTDE functions {γ MG 1 , . . ., γ MG 6 } that jointly define the perturbative properties of the model.These functions are related to the perturbation functions in the EFTDE action in Eq. (2.1) through We assume that the time evolution of these quantities follows a similar functional form to that of Ω MG (t) in Eq. (2.3) above, that is Note that Eq. (2.3) and Eq.(2.5) implicitly limit the Horndeski theory space accessible in our analysis.The constraint in Eq. = 0. Therefore, the Horndeski models are fully specified with six EFTDE parameters that control perturbations, γ MG,0 1,2,3 and s 1,2,3 , plus two EFTDE parameters that control the background, Ω MG,0 and s 0 .
In this work, we follow the "designer approach" (see e.g.[41]) to construct a Horndeski model, where we specify background cosmological parameters and construct the full theoretical model by specifying these standard cosmological parameters that control the expansion rate H(z) and the density perturbations.For the background expansion, we consider a flat ΛCDM cosmological model, specified by the physical baryon and cold-darkmatter densities (Ω b h 2 and Ω c h 2 respectively), and the constant dark-energy equation of state w = −1.Our choice of flat geometry implicitly implies that the dark energy density is given as Ω DE = 1 − Ω M .Our full model parameter space is therefore (2.6) Certain analyses, e.g.[42][43][44], fix γ MG 3 (a) = 0 to enforce that the speed of the propagation of gravitational waves (GW) be equal to the speed of light.This requirement is motivated by the constraints derived from the binary neutron-star merger events "observed" by both GW and optical instruments (see [45] for a review), e.g.GW170817 and GRB170817A [46].To better understand this, consider a simple model-independent parameterization of the speed of propagation of GW where c 2 T and c 2 are the squared speeds of GW and of light, respectively.Here α T (a) quantifies the GW speed's deviation from the speed of light, and can be mapped into the γ MG 3 function as [42] α . (2.8) From Eq. (2.8), it is clear that the GW constraint of −6 × 10 −15 < α T,0 < 1.4 × 10 −15 [45,46] translates into which is often simply taken to be γ MG 3 (a) = 0 within EFTDE (or α T = 0 in general).In this paper, we do not follow the above approach but rather, for full generality, allow for γ MG,0 3 ̸ = 0.This choice certainly merits a justification: [47] pointed out that current LIGO multi-messenger GW events have only been detected at the energy scale close to either the strong coupling scale or the EFT cut-off.They further explicitly showed that, within the EFTDE approach to Horndeski theories, the GW speed is generally a function of energy scale c T (k) (see their Eq.( 13)), and therefore can still potentially deviate from the speed of light when measured at lower frequencies (see their Fig. 1).Those Horndeski models hence do not necessarily obey the derived constraint in Eq. (2.9).Future observations of either GW events at a lower frequency, e.g. with LISA [48], or CMB B-mode polarization [49,50], will be able to place stringent constraints on the speed of GWs in these Horndeski models.Finally, we note that [41] did not find a qualitative difference between reconstructed Horndeski models with zero and non-zero γ MG 3 (a) when confronting models with current cosmological data4 .

Stability conditions in EFTDE and EFTCAMB
EFTCAMB further allows user to impose a set of consistency checks on the EFT functions in order to ensure that the EFTDE models being considered and evaluated meet the theoretical stability conditions [30].These so-called viability conditions [40], or rather viability priors in the context of cosmological inference [31], include 1. Physical stability: the EFTDE theory must have a background stable to perturbations.
In other words, the background must be free from ghost and gradient instabilities.The former corresponds to the situation where the model has a negative kinetic energy; the latter refers to the scenario where the squared sound speed is negative in some background regions.[51] (see Eqs. ( 42)-( 51) of [40] for details).
2. Mathematical stability: the EFTDE theory must have a well-defined π-field equation with no fast exponential growing modes of perturbations, as well as well-defined equations for tensor perturbations (see Eq. ( 52) of [40] for details).
3. Additional, model-specific stability: For Horndeski models, this enforces w(a) ≤ −1/3 at all time.Specifically to this work, this condition is automatically guaranteed as we consider only the case of a constant w = −1.
Generally speaking, the set of physical stability conditions is more restrictive than the mathematical ones.Further, the mathematical stability conditions implicitly assume that a) the π-field equation decouple from other field equations and b) its time-dependent coefficients evolve slowly (with time); these conditions are approximate and model-dependent.Therefore, in this work we only impose the physical conditions.We have explicitly verified through a number of pilot runs that running EFTCAMB with only physical conditions versus both physical and mathematical conditions does not qualitatively affect the range of Horndeski models successfully evaluated by EFTCAMB, hence the principal results of our work.

Growth prediction in Horndeski models
In order to draw connections between Horndeski models and observational data, we will focus on the prediction of each Horndeski model for the parameter combination f (z)σ 8 (z) ≡ f σ 8 (z).
This quantity plays a central role in describing galaxy peculiar velocities and redshift-space distortions; it is thus an excellent meeting place between observations and theories of modified gravity.For each theoretical model under consideration, we compute the exact f σ 8 (z) using EFTCAMB in bins of redshift.In this paper, we follow the convention in [52,53] and define f σ 8 (z) through where σ is the amplitude of (total) matter fluctuations obtained from the matter velocitydensity (cross-)correlation function, while σ (dd) 8 ≡ σ 8 is that same quantity obtained from the matter density-density (auto-)correlation function.Specifically,  Dividing this f σ 8 (z) by the value of σ 8 (z) (also calculated by EFTCAMB), gives the theoretically predicted growth rate f (z) of each Horndeski model, which will then be fit with the formula in Eq. (1.2) with a specified functional form of γ(z).For all Horndeski models considered in this work, Eqs.(2.10)-(2.11)or Eq.(1.1) yields the same numerical result and quantitative conclusion within the scales of interest, k ≃ 0.01 − 0.1 h Mpc −1 .This conclusion naturally follows under the assumption that growth rate is scale-independent.Even though this assumption may not hold in more generic Horndeski and modified gravity models (see e.g.[54,55]), it holds up rather well for Horndeski models we consider here, in particular within the scales probed by Stage IV and V surveys, i.e. k ≃ 0.01 − 0.1 h Mpc −1 .Within that range of k and each of the Horndeski models considered in this work, f σ 8 (z) only varies within sub-percent level at any given z.We further illustrate and discuss this point in App.B.

Testing growth parameterizations in Horndeski models
Our aim is to statistically chart a broad range of functional forms of γ(z), but only for Horndeski models that are compatible with current observational constraints.To do so, we first identify the sub-space of Horndeski theories in which models are both stable and compatible with current constraints on f σ 8 (z).Detailed description about how we carry this out can be found in App. A.
After we have determined a sub-space of Horndeski theories compatible with current data, we then follow the procedure outlined here: 1. We sample this theory sub-space, i.e. randomly draw Horndeski models from the subspace and calculate the theoretical prediction for f σ 8 (z) by each model; this is described in Sec.3.1.
2. For each model, we quantify the goodness-of-fit between the f σ 8 (z) computed using various proposed fitting formulae for γ(z) and the actual theory prediction.To compute the fit, we use the forecast constraints on f σ 8 (z) from Stage IV and V surveys; this is discussed in Sec.3.2.
Finally, we compare the goodness-of-fits and identify the best functional form for γ(z).Fig. 1 depicts the entire procedure.

Sampling and evaluating Horndeski models
Following the preliminary runs and cuts described in App.A, we identify a particular subspace of Horndeski theories where models that are both stable and consistent with current data constraints on f σ 8 : Eq. (3.1) specifies our priors on EFTDE parameters from which we sample Horndeski models in our main analysis.These priors are in broad agreement with the corresponding posteriors reported in [44].
From the Horndeski priors in Eq. (3.1) and cosmological priors in Tab. 1, we first draw and evaluate a sample of ∼20,000 Horndeski models using EFTCAMB.Of these, 19,908 models pass the stability conditions imposed within EFTCAMB (see Sec. 2).For each model successfully evaluated by EFTCAMB, we use its prediction of f σ 8 (z) to recompute the quantity χ 2 current data given in Eq. (A.2) and reject those with χ 2 current data > 5σ (same cut as in App.A).We thereby end up with a final set of 18,543 viable Horndeski models and we only use these models to test our fitting formulae and identify the most accurate one.

Testing the functional forms for γ(z)
In this section, we find the most accurate two-parameter description of the growth index that fits Horndeski models, assuming future f σ 8 data.
The final fitting formula must be sufficiently accurate even for future measurements of structure growth in the coming years and decades.We therefore assume optimistic constraints given by future data considered in this work.In this way, we impose a high burden of proof for any proposed fitting formula.Specifically, we choose the measurements of, and constraints on, f σ 8 (z) from several future surveys that together cover a redshift span up to z max = 5, as shown in Figure 2 and listed in Tab. 3. In the low-redshift region, we adopt forecasted error bars based on the Taipan Galaxy Survey [64], assuming a 5% error at z = 0.05 and a 2.7% error at z = 0.2.In the intermediate redshift range of 0.65 < z < 1.85, we adopt the DESI forecasts [65] where errors are in bins of size ∆z = 0.1 as given in Tab. 3. In the high redshift region, we use forecasts from MegaMapper5 where errors are in four bins of size ∆z = 0.75.
For every Horndeski model (generated following the protocols described in App.A), we assume future data centered on the predictions of that model, with errors representative of Stage IV and V surveys shown in Tab. 3. We then fit this simulated data with a number of distinct two-parameter fitting formulae for γ(z).To perform each fit, we employ the iminuit optimization package to find the best-fit values of the two fitting-formula parameters, γ 0 and Table 3: Constraints on f σ 8 from future surveys that cover a redshift range of up to z max = 5.This is also visualized in Fig. 2.

Redshift
Here (f σ 8 ) model (z) is the value obtained directly from EFTCAMB following the definition in Eq. (2.10); in (f σ 8 ) fit (z), f (z) was calculated by each of the two-parameter parameterizations listed in Tab. 4 and σ 8 (z) obtained from EFTCAMB following the definition in Eq. (2.11); σ future data i is the error on the i-th future measurement at redshift z i , both of which are given in Tab. 3. Finally, we select the fitting formula of γ(z) that gives the best goodness-of-fit across all sampled Horndeski models.

Results
We now present the two principal results of this paper.In Sec.4.1, we show the performance of different parameterizations of γ(z) in their fit to theoretical predictions of f σ 8 in Horndeski models and identify the best fitting formula.In Sec.4.2 we constrain the parameters of the best fitting function using current cosmological data.

And the winner is...
To first get a sense of what kind of γ(z) is typically predicted by Horndeski models, we numerically compute the redshift-dependent growth index for a limited number of models.Specifically, we evaluate the true growth index given f (z) and Ω M (z) in that model.Fig. 3 shows the exact γ(z) computed from Eq. (4.1) in 50 Horndeski models randomly selected from our prior.To guide the eye we also plot the ΛCDM growth index which, as expected, is very well approximated by γ(z) ≃ 0.55.
The general behavior of γ(z) in Horndeski models at z ≳ 1 can be easily understood from Eq. (4.1): as z increases, Ω M (z) → 1 in cosmological models without early dark energy (which is true for all models considered in this paper).Therefore, for any given (Horndeski) model, departures of the growth rate f (z ≳ 1) from the ΛCDM prediction are generally associated with relatively large fluctuations in the growth index simply because the latter is the exponent of Ω M (z), which is close to unity.Overall, the results in Fig. 3 not only show a clear redshift evolution of the growth index expected in these models but also demonstrate that this redshift dependence is fairly featureless.
The results in Fig. 3 motivate our selection of specific functional forms of the growth index and Tab. 4 enumerates these functions.In addition to the constant growth index and (for pure simplicity) the one going linearly with z, we also try several other forms that contain simple polynomials in redshift, as well as simple logarithmic or exponential terms.Since we would like to parameterize f σ 8 (z) up to z ≃ 5 where Stage IV and V data will constrain growth, trends in Fig. 3 emphasize that we need a nonlinear redshift-dependent parameterization, such as those suggested in Tab. 4.
The success (or failure) of each parameterization in fitting theoretical predictions of Horndeski models is further reported in Tab. 4. For each fitting function, we summarize the statistics of the quantity χ 2 fit , defined in Eq. (3.2), measured for our set of ∼18,000 Horndeski

ΛCDM Model Horndeski Models
Figure 3: The exact redshift-dependent growth index, γ(z) ≡ ln f (z)/ ln Ω m (z), numerically evaluated for 50 randomly selected Horndeski models from our prior out to z max ≃ 5.The red nearly horizontal line shows the exact γ(z) for the ΛCDM model.These results demonstrate that one needs a nonlinear multi-parameter parameterization to capture the features of growth index at high redshift in modified gravity.They also motivate functional forms for our trial fitting functions in Tab. 4. models.The summary is provided by two statistical measures: the median of the χ 2 fit values, and the 95-th percentile (i.e. the upper bound of 95% of values of χ 2 fit ).Since the theoretical data vector calculated by EFTCAMB is noiseless, a perfect fit of a fitting formula to true f σ 8 (z) of a theory model will have χ 2 fit = 0; this explains the generally small chi-squared values in Tab. 4. In general, we find that the distribution of the χ 2 fit values has a heavy tail in each instance; this explains why the 95% upper bounds are typically much larger than the corresponding medians in Tab. 4. The presence of the heavy tails reflects the improvement of f σ 8 constraints going from Stage III to Stage IV and V surveys: our model selection cut, Eq. (A.2), is only concerned with current constraints on f σ 8 , while our comparison by Eq. (3.2) is concerned with forecast constraints for Stage IV and V surveys.
As shown by the highlighted cell in Tab. 4, the best fitting formula for the growth index in the redshift range of 0 < z < 5 is This two-parameter fitting formula fits the future data with a median χ 2 fit of 0.028, which is about 40 times smaller than the median χ 2 fit with the constant growth index γ(z) = γ 0 .Furthermore, we find that with the new fitting formula from Eq. (4.2), the maximum deviation in f σ 8 (z) at any redshift between the fitting formula's approximation and Horndeski theory's true value for f σ 8 (z) has a median of 0.4% when averaged over all models.When using the traditional one-parameter growth index, γ(z) = γ 0 , the median of maximum differences per model is 2.5%.Our two-parameter, redshift-dependent fitting formula therefore approximates the theoretical predictions of f σ 8 (z) in Horndeski theories about six times better than the one-parameter, constant-γ case, leading to an improvement of ∼ 40 times in χ 2 .Table 4: Proposed fitting functions and the statistics of their fit to our sample of ∼18,000 Horndeski models.

Fitting function
Best-fit χ 2 for γ(z) Median 95% Percentile γ 0 1.16 36.6 γ 0 + γ 1 z 0.046 4.00 Several other fitting functions in Tab. 4, also do a good job, in particular γ(z) = γ 0 +γ 1 z comes close in the median, but falls short in fitting models near the tail; however, none are as good as the form in Eq. (4.2).
The result that the median χ 2 fit ≪ 1 with the best fitting function is very encouraging, as it implies that the contribution of the inaccurate fitting function to the bias in cosmological parameters will be subdominant.Specifically, our finding that χ 2 fit ≪ 1 implies that even the best-determined direction in parameter space will be biased by ≪ 1σ in our optimistic-data case.
In Fig. 4, we further showcase the performance of a few proposed fitting formulae on one randomly selected Horndeski model.It is evident that the prize-winning fitting form in Eq. (4.2) is the best of the fitting functions shown.We also observe that the best fitting function does a good job both at z ≪ 1 and at z > 1; both of these ranges are required to be accurately fit for the (γ 0 , γ 1 ) description to be a useful tool for the Stage IV and V surveys.
In the previous section, we have proposed and validated Eq. (4.2) as a new fitting function of the growth index γ(z) for future surveys.Using current cosmological data, we now demonstrate the applicability of this formula in consistency tests of general relativity and flat ΛCDM.
Building upon the work in [9], here we constrain the growth index γ(z), specifically (γ 0 , γ 1 ) in Eq. (4.2), from a combination of large-scale structure and CMB data sets: measurements of f σ 8 through peculiar velocities and redshift-space distortions6 [58-60, 66-73], measurements of baryon acoustic oscillation (BAO) from the Six-degree Field Galaxy Survey (6dFGS; [74]) and the Sloan Digital Sky Survey (SDSS; [73,75,76]), 3x2pt correlation functions from the Year-1 analysis of Dark Energy Survey (DES-Y1; [77]), and CMB measurements from Planck 2018 [53].In this work, we additionally include the type Ia supernovae data sets and likelihoods from Pantheon [78] which however make very little difference in our final constraints on (γ 0 , γ 1 ).To obtain constraints on the growth-index and cosmological pa- We also show the best-fitting f σ 8 (z) for three fitting formulae for the growth index γ(z), as well as the f σ 8 fit from a constant growth index, γ(z) = γ 0 .Bottom panel : Relative difference between the true f σ 8 (z) and the best-fit results of each fitting formula shown in the top panel.The shaded grey area shows the forecast statistical errors associated with each survey.rameters (after numerically marginalizing over nuisance parameters), we use cobaya 7 , which provides out-of-the-box access to most likelihoods for the aforementioned data sets, validated against their official analyses.The only exception is the f σ 8 likelihood for [59, 60, 66? -72], which we implement in [9], assuming a Gaussian likelihood and a diagonal covariance.
Our implementation of the growth index largely follows that of [9].Specifically, at any given redshift z, we re-scale the linear matter power spectrum as where D(γ, z) is numerically integrated from Eqs. (1.1)-(1.2),and P (k, z = 0) is the fiducial linear matter power spectrum evaluated at z = 0 which is specified by the standard set of cosmological parameters: where A s and n s are the amplitude and the spectral index of the primordial power spectrum, τ is the reionization optical depth, and θ MC is (an approximation to) the angular size of the sound horizon at recombination.We emphasize that this set of cosmological parameters is jointly constrained with (γ 0 , γ 1 ).When sampling with cobaya, we compute P (k, z) using the cosmological Boltzmann solver CAMB [79,80].We validate our implementation by reproducing, up to a high precision, the constraints on the standard cosmological parameters in the baseline analyses of Planck 2018 [53] and DES-Y1 [77].
Motivated by the fact that Eq. (1.2) has, so far, been validated only for sub-horizon perturbations, we exempt the primary CMB anisotropies from the rescaling in Eq. (4.3).In other words, the growth index never directly affects the unlensed CMB power spectra, but rather only the CMB lensing potential.Consequently, only the CMB lensing amplitude is sensitive to any change in the growth index8 γ.
For the cosmological parameters, we adopt the same priors as specified in the Planck 2018 baseline analysis [53], which considered flat ΛCDM at fixed neutrino mass m ν = 0.06 eV.Priors on all nuisance parameters also follow those in the official analyses of the corresponding data sets.We choose uniform priors on the two growth-index parameters: γ 0 ∈ U(0.0, 2.0), γ 1 ∈ U(−1.0, 1.0).
In Fig. 5 we present the constraints in the γ 0 − γ 1 plane, marginalized over all other cosmological and nuisance parameters.Allowing for redshift evolution, which is effectively controlled by the parameter γ 1 in γ(z) = γ 0 +γ 1 z 2 /(1+z), we observe the expected degeneracy between γ 1 and γ 0 in that parameterization.Specifically, we infer γ 0 = 0.621 ± 0.03 and γ 1 = 0.149 ± 0.235.We find evidence for a disagreement with the standard cosmological model -which predicts (γ 0 = 0.55, γ 1 = 0) -at approximately 99.8% level (corresponding to about "3.1-sigma" in a two-tailed test of statistical significance).Our finding is therefore in good statistical agreement with the conclusion from the analysis in [9] which however assumed γ = const., i.e. no redshift evolution.
Our γ 0 constraint suggests that the growth rate of large-scale structure is suppressed recently -with the onset of dark energy -relative to the prediction by flat ΛCDM and general relativity, while the γ 1 constraint implies no (strong) evidence of redshift evolution in the growth index.

Summary and Conclusions
In order to squeeze out stringent constraints on the growth of structure from data, it will be crucial to have precise parameterizations of the evolution of the growth of structure, specifically with the goal to cleanly separate it from the background evolution.One such parameterization is f (z) = Ω M (z) γ with a constant growth index γ which, while being highly accurate for dark-energy models close to ΛCDM, is no longer such for modified gravity.In this work, we have promoted γ to a function of redshift, i.e. γ(z).We have further identified and validated the best two-parameter fitting formula for γ(z) that accurately describes the growth of structure across the landscape of Horndeski theories of modified gravity: We have explicitly shown that Eq. (5.1) fits the theoretical predictions of f σ 8 (z) by Horndeski models with typical errors at the sub-percent level, well within the precision that will be reached by Stage IV and Stage V surveys.Further, as a demonstration, we have constrained the parameters of Eq. (5.1) using modern data from galaxy clustering, weak lensing, CMB, and type Ia supernovae.The result we obtained is in tension with the concordance cosmological model of (γ 0 , γ 1 ) = (0.55, 0)which was expected given such indications in our recent analysis which essentially assumed γ(z) = γ 0 [9].Specifically, we have found evidence that γ 0 > 0.55, while the posterior of γ 1 peaks at positive values but is statistically consistent with zero.
We conclude that forthcoming data from ongoing and upcoming large-scale structure surveys [64,[82][83][84][85][86] will dramatically expand the redshift coverage and increase the precision of the growth-of-structure sector.This, in turn, will enable new opportunities to test the selfconsistency of the standard cosmological model, and possibly detect deviations from general relativity from such measurements of the growth rate.Our new fitting formula should provide one reliable meeting point between data and theory.

Data and software products
The modified version of CAMB that implements our parameterization, γ(z) = γ 0 +γ 1 z 2 /(1+z), is available at: github.com/MinhMPA/CAMB_GammaPrime_Growth.git.To obtain cosmological constraints from data using this parameterization, one can use the fork of cobaya at: github.com/MinhMPA/cobaya.git.To request the f σ 8 (z) data evaluated for the Horndeski models considered in this paper, kindly send a message to Y.W. and N.-M.N.
2. In the second step, we then choose to exclude cosmological models -specified by the above cosmological parameters and EFT parameters -that are disfavored by current data at ≥ 5σ.The current f σ 8 data that we use are shown in Tab. 2. We define the goodness of fit to the theoretical model value as where (f σ 8 ) model (z) is the value obtained from EFTCAMB, and (f σ 8 ) current data , σ current data i , and z i are respectively the measurement, error, and redshift of current data, all of which are given in Tab. 2. Similar to [87], we typically find that -although model stability can strongly depend on γ MG,0 1 and s 1 -model prediction (in our case, for f σ 8 ) does not.With 20 f σ 8 measurements and six EFT parameters to be fit, we have N dof = 14 degrees of freedom.Assuming a normal distribution of individual f σ 8 measurements, keeping the models that are within 5σ from current measurements then requires χ 2 ≤ 60.
During this process, we also prune regions of the Horndeski parameter space where models largely fail the stability conditions specified in Sec.2.2.We thereby end up with the prior ranges specified in Eq. (3.1).

B Scale dependence of growth
Here we take a closer look at the scale dependence of the growth rate f (z, k) specifically in the context of Horndeski models.As mentioned near the end of Sec. 2, this scale dependence is not guaranteed to be negligible for growth in models beyond smooth dark energy.By using EFTCAMB we can straightforwardly investigate the effect by directly computing the linear growth as the ratio of the matter transfer functions at two different redshifts D(z, k) = T (z, k) T (z = 0, k) (B.1) and then numerically evaluating f (z, k) from Eq. (1.1).
We do find significant k dependence of f (z, k) near the horizon scale (k ≃ 0.0001 − 0.001h Mpc −1 ).However, note that most cosmological observations of large-scale structure come from smaller scales, roughly k ≃ 0.01 − 0.1h Mpc −1 (see, for example, Fig. 19 in the Planck legacy paper [88]).Therefore, we do not need to take into account the scale dependence of f (z) if we focus on this range of scales.This is clearly demonstrated in Fig. 6 where we show the k dependence of f σ 8 (z) on two scales, k = 0.01 and 0.1 h Mpc −1 , for a selection of a few Horndeski models from our priors (as well as for ΛCDM).When selecting the Horndeski models to showcase in Figure 6, we sequentially increase each EFT parameter to its largest value allowed by the priors specified in Eq. (3.1).
The maximum scale dependence that we observe in Fig. 6 is about 0.5% -much lower than the (statistical) errors in f σ 8 of all surveys considered in this work.We thus demonstrate that, for Horndeski models that deviate most from general relativity (and potentially have the most scale dependence) but are still within the our selected priors, the scale-dependent differences in f σ 8 (z) are well below 1% and much smaller than stringent constraints from future surveys.
We also quantitatively demonstrate the scale-independence of Horndeski models within our specified priors.Among ∼18,000 Horndeski models, we find that the difference between the values of f (z, k) evaluated at k = 0.01 and 0.1 h Mpc −1 (across all models and all redshifts) has a median of 0.3% and a 95% percentile of 0.5%, both at a sub-percent level.
Nevertheless, there are reasons why studying the scale dependence of the growth of structure is very interesting and should be pursued.First, modified-gravity models that are physically different from models in our Horndeski prior may lead to a much more significant scale dependence.Second, observations that probe larger spatial scales (say k ∼ 0.001 h Mpc −1 ) might also be able to observe this scale-dependent behavior.

( 2 . 2 ) corresponds to 2γ MG 5 = γ MG 3 = −γ MG 4 and γ MG 6
) where x denotes either the v or d component, W TH (k, R = 8h −1 Mpc) is the Fourier transform of the spherical top-hat window function of radius R = 8h −1 Mpc, T x (k) is the transfer

Figure 1 :
Figure 1: Flowchart describing how we sample Horndeski models from EFTDE theory space, evaluate the models and test the fitting functions for γ(z) against their predictions for f σ 8 (z).

Figure 2 :
Figure 2: Future f σ 8 (z) error forecasts that we use to assess the accuracy of the fitting formulae for the growth rate.The width of each line segment indicates the size of the redshift bin while the height shows the magnitude of the forecast error.

Figure 6 :
Figure 6: Top panel: f σ 8 (z) of a ΛCDM and three Horndeski models calculated by EFTCAMB at scales k = 0.1 h Mpc −1 and k = 0.01 h Mpc −1 .Bottom panel: The percent difference between f σ 8 (z) of each Horndeski model calculated at these two scales.The shaded area illustrates the constraints on f σ 8 (z) by future surveys in terms of percent error in each redshift bin, which we use in determining the goodness-of-fit of each proposed fitting formula.For visualization purposes, we scale the errors down to 1/5th of their true size in the bottom panel.

Table 1 :
[53]cial values of cosmological parameters and priors on them used in our sampling of Horndeski models.They closely follow the ΛCDM best-fit values and 68% in the Planck 2018 analysis[53].

Table 2 :
Current measurements of f σ 8 and errors at different redshifts.The data include the 6dF Galaxy Survey (6dFGS), peculiar velocities of type Ia supernovae (SNIa), Galaxy and Mass Assembly (GAMA), the WiggleZ Dark Energy Survey, the Baryon Oscillation Spectroscopic Survey (BOSS), the extended Baryon Oscillation Spectroscopic Survey (eBOSS) and the VIMOS Public Extragalactic Redshift Survey (VIPERS).