Defying the laws of Gravity I: model-independent reconstruction of the Universe expansion from growth data

TL;DR

This paper employs model-independent methods to reconstruct the Universe's expansion history from growth data, confirming consistency with the standard cosmological model and General Relativity.

Contribution

It introduces two novel, model-independent techniques—crossing statistics and Gaussian processes—for reconstructing cosmic growth and expansion histories from observational data.

Findings

Results align with the flat-$\Lambda$CDM model

Constraints on matter density and fluctuation amplitude are consistent with current cosmology

Current data supports General Relativity as the gravity theory

Abstract

Using redshift space distortion data, we perform model-independent reconstructions of the growth history of matter inhomogeneity in the expanding Universe using two methods: crossing statistics and Gaussian processes. We then reconstruct the corresponding history of the Universe background expansion and fit it to type Ia supernovae data, putting constraints on $(\Omega_\mathrm{m,0},\sigma_{8,0})$. The results obtained are consistent with the concordance flat-$\Lambda$CDM model and General Relativity as the gravity theory given the current quality of the inhomogeneity growth data.

Figures

Equations

$$h^2(z)$$

$$\ddot{\delta }+2H\dot{\delta }$$

$$\dot{x}$$

$$\frac{\mathrm{d}{}^2\delta }{\mathrm{d}\ln a^2}$$

$$f$$

$$f{\sigma }_8$$

$${\sigma }_R^2(z)$$

$$h^2(z)$$

$$\frac{{\delta }^{\prime }(z)}{{\delta }_0}$$

$$\delta (z)$$

$$\begin{gathered} {\hat{h}}^2(z|\widehat{f{\sigma }_8},{\Omega }_{\text{m,0}},{\sigma }_{8,0})=\frac{(1+z{)}^4}{{\widehat{f{\sigma }_8}}^2(z)}({\widehat{f{\sigma }_8}}_0^2 \\ -3{\Omega }_{\text{m,0}}{\int }_0^z\left({\sigma }_{8,0}-{\int }_0^u\widehat{f{\sigma }_8}(v)\frac{\mathrm{d}v}{1+v}\right)\widehat{f{\sigma }_8}(u)\frac{\mathrm{d}u}{(1+u{)}^2}). \end{gathered}$$

$$d_{\text{L}}(z)$$

$$\mu (z)$$

$${\Omega }_{\text{de}}(z)$$

$$\delta (z)\propto (1+z{)}^{\frac{1}{4}\left[1-\sqrt{25-24{\Omega }_{\text{de}}}\right]}\approx (1+z{)}^{-1+\frac{3}{5}{\Omega }_{\text{de}}}.$$

$${\chi }_{\text{tot}}^2$$

$$\widehat{f{\sigma }_8}(x)$$

$$T_N(x|C_0,\dots ,C_N)$$

$$x$$

$$\left[\begin{array}{c}\text{y}\\ {\text{f}}_{\ast }\end{array}\right]$$

$$k_{{\sigma }_f,\ell }(\text{x},\text{y})$$

$$\forall z\in [0,2]{\Omega }_{\text{de}}(z)>0,$$

$$\{\begin{array}{l}\forall z\in [0,1]{\Omega }_{\text{de}}(z)>0\\ \exists z\in [1,2]{\Omega }_{\text{de}}(z)<0\end{array}$$

$$\{\begin{array}{l}\forall z\in [0,0.7]{\Omega }_{\text{de}}(z)>0\\ \exists z\in [0.7,1]{\Omega }_{\text{de}}(z)<0\end{array}$$

$$w(z)$$

$$\gamma (z)$$

Full text

††pubyear: 2019††pagerange: Defying the laws of Gravity I: Model-independent reconstruction of the Universe expansion from growth data–A

1 Introduction

The discovery at the end of last century of the late-time accelerated expansion rate of the Universe raised the question of its physical cause. There are two great avenues towards solving this problem: either that it is due to an unknown new physical component dubbed (physical) dark energy, or its origin lies in a modification of the laws of gravity (e.g., Sahni & Starobinsky, 2000; Yoo & Watanabe, 2012; Clifton et al., 2012). However, they represent particular cases of a more general situation like in scalar-tensor gravity, when both a new physical field is introduced for dark energy description and gravity is modified, too (see e.g. Boisseau et al. (2000); Copeland et al. (2006); Sahni & Starobinsky (2006)).

In the concordance model the role of dark energy is played by a cosmological constant $\Lambda$ while gravity is described by Einstein’s theory of General Relativity (GR). While GR has been remarkably successful to explain observations in the Solar system (see e.g. Gannouji et al. (2006)), its successful extrapolation to much larger cosmic scales remains unclear. Modified gravity models with modifications of gravity on cosmic scales are not excluded and could well be the solution to the recent acceleration of the Universe. The nature of dark energy and therefore also the correct model of gravity are burning issues of cosmology and theoretical physics in general.

The large-scale structures of the Universe are an ideal laboratory to test gravity, and to distinguish between physical dark energy and modified gravity (which may be also called geometrical dark energy as in Sahni & Starobinsky (2006)). In particular, redshift-space distortion (RSD) due to galaxy peculiar velocities can be used to estimate the growth factor $f$, which is a key to understanding gravity. However, as a first step, we restrict ourselves to GR. For a flat-FLRW Universe with dark energy as a perfect fluid with equation of state $w(z)$, the expansion history $h(z)=H(z)/H_{0}$ is described by

[TABLE]

In GR, the evolution of the matter overdensity $\delta(\textbf{{x}},z)=(\rho(\textbf{{x}},z)-\bar{\rho}(z))/\bar{\rho}(z)$ are governed in the Newtonian approximation by

[TABLE]

where dot stands for a derivative with respect to cosmic time $t$, and $\Omega_{\text{m}}(z)=\Omega_{\text{m,0}}(1+z)^{3}/h^{2}(z)$ is the matter density normalized by the critical density. From eqs. (1) and  (2), it is clear that changing the expansion will also affect the growth evolution. In fact, Starobinsky (1998) showed that the Universe expansion history $H(z)$ can be also reconstructed from $\delta(z)$ unambiguously in this case (the situation becomes more complicated in scalar-tensor gravity (Boisseau et al., 2000)).

In this paper, we aim first to reconstruct the growth history from data, and then to use it in order to deduce the expansion history (assuming GR) and to compare it with the supernovae data. We should note here that there are two different reconstructions involved here: a statistical reconstruction of the growth factor $f(z)$ from observational data on one hand, and on the other hand, the theoretical reconstruction of the background expansion $H(z)$ from $f(z)$. Lee (2014) applied the Starobinsky (1998) formalism by integrating from $z$ to a maximum redshift $z_{\mathrm{max}}$ and henceforth obtain $\Omega_{\text{m,0}}$, and studied the effect of $z_{\text{max}}$ on the error budget, showing it dominates over the data uncertainties. However, this approach differs from ours in that we integrate from [math] to $z$ and use $\Omega_{\text{m,0}}$ as a free parameter. A more similar approach to ours was recently applied in Yin & Wei (2018). The theoretical reconstruction is described in § 2, and the statistical reconstructions together with the results are presented in § 3, and our conclusions are drawn in § 4. We validate the method on mock data in § A.

2 Method

2.1 Theoretical Framework

From eq. (2), using that for any function $x(t)$,

[TABLE]

where the last equality defines the growth index $\gamma$. In general, $f=f(\textbf{{k}},z)$ and therefore $\gamma=\gamma(\textbf{{k}},z)$ (e.g., Gannouji et al., 2009). However, in GR and for dust-like matter, $\gamma$ is k-independent and has weak dependence on $z$, so that $\gamma(z)\simeq 0.55$, with a slight dependence on the equation of state of dark energy $w$ and the matter density parameter $\Omega_{\text{m,0}}$. Note, however, that $\gamma$ may not be exactly $z$-independent for quintessence (a scalar field with a potential minimally coupled to gravity) models of dark energy as was shown in Polarski et al. (2016).

Observational data on RSD provide us with the product $f\sigma_{8}$, where

[TABLE]

is the rms of the density fluctuations smoothed over a radius $R$, usually taken to be $8\,h^{-1}\mathrm{Mpc}$.

In Shafieloo et al. (2018) (see also L’Huillier et al., 2018), we used Pantheon and a compilation of growth data to put model-independent constraints on $(\Omega_{\text{m,0}},\sigma_{8,0},\gamma)$, where subscript 0 stands for the current value. However, in that paper, we treated $\gamma$ as a constant, effectively performing a consistency test of GR. In fact, one can solve the problem without assuming $\gamma=$ constant. Assuming $\delta$ and $\delta^{\prime}$ are known, the expansion history $h$ can be uniquely determined via (Starobinsky, 1998)

[TABLE]

where ′ denotes a derivative with respect to $z$ (not $\ln a$).

Thus, one can obtain from RSD measurements

[TABLE]

Therefore, for a given reconstruction $\widehat{f\sigma_{8}}(z)$ together with a given $(\Omega_{\text{m,0}},\sigma_{8,0})$, the expansion history $\hat{h}(z)$ is uniquely determined by

[TABLE]

Therefore, reconstructing $\widehat{f\sigma_{8}}(z)$, and exploring the $(\Omega_{\text{m,0}},\sigma_{8,0})$ parameter space, we can reconstruct $h(z)$ and the luminosity distance. In the following, we assume a flat-FLRW universe, which is consistent with current model-dependent (Farooq & Ratra, 2013; Planck Collaboration et al., 2018) and model-independent(Räsänen, 2014; Räsänen et al., 2015; L’Huillier & Shafieloo, 2017; Denissenya et al., 2018; Shafieloo et al., 2018) constraints. The luminosity distance is thus

[TABLE]

we can then fit the reconstructed $\mu$ to SNIa data, obtain a $\chi^{2}$, and hence put model-independent constraints on $(\Omega_{\text{m,0}},\sigma_{8,0})$.

2.2 Fit to the data

After obtaining $f\sigma_{8}$ from GP, we then obtain the expansion history $h$ via eq. (8), for a given choice of $(\Omega_{\text{m,0}},\sigma_{8,0})$. We note that, a priori, for a given $(f\sigma_{8},\Omega_{\text{m,0}},\sigma_{8,0})$, there is no guarantee for $h^{2}(z)$ to be larger than the matter term $\Omega_{\text{m,0}}(1+z)^{3}$. At low redshifts, local measurements yield $\Omega_{\text{m,0}}<1$, and thus, in the absence of spatial curvature, the energy density of dark energy is positive. However, the energy density of dark energy could be negative in the matter-dominated era. An effective negative $\rho_{\text{de}}$ was considered in e.g. Boisseau et al. (2015). Following these considerations, hereafter, we impose the positive dark energy condition that

[TABLE]

In addition, even when $\Omega_{\text{de}}<0$, we impose the conservative choice $\Omega_{\text{de}}>-0.05$. If we assume tracking dark energy during the matter era, CMB constraints which fix the perturbations amplitude at $z\approx 1100$ give a slightly stronger lower bound $\Omega_{\text{de}}>-0.025$ if we restrict the relative additional growth to no more than $10\%$. Indeed, we get for constant $\left|\Omega_{\text{de}}\right|\ll 1$ during the matter dominated stage

[TABLE]

This anomalous growth is reminiscent of the growth in the presence of massive neutrinos (where there is a weaker growth with expansion) with the essential difference that here the growth is boosted due to a negative $\Omega_{\text{de}}$.

In order to study the effect of $z_{\text{max}}$, we choose three values for ${z_{\text{max}}}=$ 0.7, 1, and 2. We note that $w$ will then diverge when $\Omega_{\text{de}}$ crosses zero (except for the degenerate case where the DE pressure also vanishes at the same moment).

Finally, calculate the $\chi^{2}$ for all the data:

[TABLE]

3 Results

We used the same data sets as in Shafieloo et al. (2018): the Pantheon compilation (Scolnic et al., 2018), and the compilation of RSD data including: 2dFGRS (Song & Percival, 2009), WiggleZ (Blake et al., 2011), 6dFGRS (Beutler et al., 2012), VIPERS (de la Torre & Peacock, 2013), the SDSS Main galaxy sample (Howlett et al., 2015), 2MTF (Howlett et al., 2017), BOSS DR12 (Gil-Marín et al., 2017), FastSound (Okumura et al., 2016), and eBOSS DR14Q (Zhao et al., 2019).

3.1 Crossing statistics

In this section, we study the effects of distorting the mean function on the final fit. In practice, this is equivalent to applying the Bayesian interpretation of the crossing statistics formalism to the reconstructed growth history (Shafieloo et al., 2011; Shafieloo, 2012a, b). In this formalism, the prediction from the theory to be tested, $\Lambda\text{CDM}$+GR in the present case, is multiplied by some hyperfunction $T_{N}(x|C_{0},\dots,C_{N})$:

[TABLE]

$P_{i}(x)$ is the $i$th order Chebyshev polynomial of the first kind, and $C_{i}$ are free parameters. The Chebyshev polynomials constitute an orthogonal basis for $x\in[-1,1]$, and as such, can represent any function. The zeroth order controls the absolute scaling, the first order the tilt, and higher order introduce a curvature and inflexion points. In the Bayesian interpretation of the crossing statistics, we are interested in the confidence intervals around the hyperparameters $C_{i}$. If the $C_{i}$ are consistent with $C_{0}=1,C_{i}=0(i\geq 1)$, the data have no preference for any departure from the mean function, meaning the model is consistent with the data. In case of significant deviation from $C_{0}=1,C_{i}=0(i\geq 1)$, the data suggest a preferred deformation of the mean function.

Distorting the starting $f\sigma_{8}$, we obtain $\hat{\mu}$ and fit both $\hat{\mu}$ and $\widehat{f\sigma_{8}}$ to the data. As noted by Hazra & Shafieloo (2014), going towards too high orders, one might miss the effects of the lower orders. Therefore, we start by limiting to the first order, i.e., tilting the mean function. We used the emcee Monte-Carlo Markov Chain (MCMC) package (Foreman-Mackey et al., 2013) to explore the parameter space $(\Omega_{\text{m,0}},\sigma_{8,0},C_{0},C_{1})$, and show the posteriors in Fig. 1. They are consistent with $C_{0}=1,C_{1}=0$, i.e., the data suggest no modifications, and are perfectly consistent with the best-fit $\Lambda\text{CDM}$+GR model. It is interesting to notice that the preferred $\Omega_{\text{m,0}}$ is rather high with respect to the Planck value, while the preferred $\sigma_{8,0}$ is low. We checked that when going to higher orders in the crossing functions, i.e., including $C_{2}$ and $C_{3}$, the contours and do not suggest further modifications (i.e., is still consistent with $C_{0}=1,C_{i>0}=0$).

Fig. 2 shows a random selection of 16 crossing functions and their effects on the reconstructed $\mu$ and on $(\Omega_{\text{m,0}},\sigma_{8,0})$. Their corresponding crossing parameters $C_{i}$ are shown as coloured points in Fig. 1. This gives the reader some intuition on how distorting the growth affects the reconstructed expansion. Small distortions from the best-fit $\Lambda\text{CDM}$+GR case can lead to significantly different $\hat{\mu}$.

3.2 Gaussian Process regression

We used Gaussian Process regression (GP, Rasmussen & Williams, 2006) to reconstruct $f\sigma_{8}(z)$ from the RSD measurements. GP have been increasingly used in cosmology (Holsclaw et al., 2010a, b, 2011; Shafieloo et al., 2012, 2013; Joudaki et al., 2018; L’Huillier et al., 2019) and other fields of astronomy (e.g. Iyer et al., 2019). A Gaussian process is effectively a random sampling on a function space, generalizing random numbers. GP can be used to reconstruct a smooth function $\textbf{{f}}_{*}$ at the test points $\textbf{{x}}_{*}$ given a discrete set of observations $(x_{i},y_{i})$ and a data covariance matrix C. For a given kernel, the covariance between pairs of random variables u and v is thus given by K(f(u),f(v)) = k(u,v), where k(u,v) is the covariance kernel. The joint-distribution of the training (observed) outputs y and the test (reconstructed) output $\textbf{{f}}_{*}$ is a Gaussian joint distribution given by

[TABLE]

where C is the covariance of the data.

We use the squared exponential kernel defined as

[TABLE]

where $(\sigma_{f},\ell)$ are two hyperparameters controlling the amplitude and the correlation scale, and thus the deviation from the mean function. For a given $(\sigma_{f}^{2},\ell)$, we can thus generate a number of samples of $f\sigma_{8}$ at any redshift $z$.

In practice, we start from the best-fit $\Lambda\text{CDM}$ as a mean function, use GP as a sampling of possible growth histories, and then apply the formalism from § 2.

In order to prevent fitting the noise, we impose a hard prior on $\ell\in[0.2,1]$.

A notable difference with the work of Yin & Wei (2018) is in the GP regression itself. While they obtain the mean and one-sigma contours, we follow each individual random sampling of the function space, therefore obtaining a set of plausible (and self-consistent) couples of expansion and growth histories. These two approaches are mathematically equivalent. However, in this case, since most of the constraining power comes from the SNIa, we do not train the GP on the $f\sigma_{8}$ data. This approach allow more flexibility in $f\sigma_{8}$, which in turn allow more flexibility in the distance moduli.

The left-hand panel of Fig. 3 shows the reconstructed $\widehat{f\sigma_{8}}$ as well as the best fit in dashed lines. The shadowed area shows the envelope of the reconstructed $\widehat{f\sigma_{8}}$. A vast majority of these reconstructions do not fit the data ($\chi^{2}>\chi^{2}_{\text{ref}}$), therefore, we show in solid lines those reconstructions of $\widehat{f\sigma_{8}}$ that, together with some appropriate $(\Omega_{\text{m,0}},\sigma_{8,0})$, yields a $\chi^{2}_{\text{tot}}<\chi^{2}_{\text{ref}}$. In order to study the effect of $z_{\text{max}}$ on the positive DE condition (14), we then separate the samples into three disjoint sets, according to the three values of $z_{\text{max}}$:

[TABLE]

and show these in different colours. The middle and right-hand panels show the reconstructed $h(z)$ and $\mu-\mu_{\text{ref}}$, with the same convention.

Due to the oscillations in the reconstructed $\widehat{f\sigma_{8}}$ from GP, the reconstructed shapes of $h$ and $\mu$ have more flexibility than $\Lambda\text{CDM}$, yielding possible better fit to the data, therefore representing a non-exhaustive set of plausible expansion and growth histories. It is worth mentioning here that, since the process of reconstructing $h^{2}$ from $\widehat{f\sigma_{8}}$ via eq. (8) involves two integrals, it is very sensitive to variations in the growth history and in the choice of $(\Omega_{\text{m,0}},\sigma_{8,0})$. In practice, most reconstructed $\widehat{f\sigma_{8}}$ cannot yield any $\hat{h}$ that fits the SNIa data. This can be seen by the large grey envelope in the three panels of Fig. 3 compared to the thinner band of allowed reconstructions. Therefore, it is important to explore the $(\sigma_{f}^{2},\ell)$ parameter space and generate a large number of random realizations.

Fig. 4 shows the area of the $(\Omega_{\text{m,0}},\sigma_{8,0})$ parameter space in cases A, B, and C that, for at least one reconstructed $\widehat{f\sigma_{8}}$, yields $\chi^{2}<\chi^{2}_{\mathrm{ref}}$, i.e., which is within the $1\sigma$ region of the best-fit $\Lambda\text{CDM}$ model. In case A, the extent is maximal: for a large range of the parameters $(\Omega_{\text{m,0}},\sigma_{8,0})$, one may find reconstructions with $\chi^{2}<\chi^{2}_{\text{ref}}$ (and $\Omega_{\text{de}}>0$). As we start allowing $\Omega_{\text{de}}$ to become negative, the contours shrink. It is apparent in case C, where $\Omega_{\text{de}}$ crosses 0 between $z=0.7$ and 1, which is only possible for large values of $\Omega_{\text{m,0}}$ (for a fixed $h(z)$, increasing $\Omega_{\text{m,0}}$ allows $\Omega_{\text{de}}$ to cross 0). It is worth noting that in cases B and C, $\Omega_{\text{de}}$ crosses 0 at least once, but can become positive again.

In Appendix A, we demonstrate the validity of the algorithm on a simulated data set, successfully recovering the input cosmology.

It is worth noting that, since we need to assume a value for $\Omega_{\text{m,0}}$ in order to reconstruct $h$, we can then obtain the (uniquely defined) equation of state $w(z)$ and growth rate $\gamma(z)$, as

[TABLE]

We note that, since two consecutive integrations are needed to obtain $h$, then in theory only one is needed to obtain $h^{\prime}$, making this method potentially less sensitive to numerical noise, provided that the quality of the $f\sigma_{8}$ data improve substantially. In addition, since the positive dark energy condition (14) ensures that $1-\Omega_{\text{m}}(z)\geq 0$ only up to $z_{\text{max}}$, $w$ might diverge. For these reasons, instead of $w$, we reconstruct $\Omega_{\text{de}}$, which does not involve derivatives of $h$. Fig. 5 shows our obtained reconstructions of $\Omega_{\text{de}}$ and $\gamma$. All these reconstructions have $\chi^{2}<\chi^{2}_{\Lambda\text{CDM}}$ and obey the positive dark energy condition for $z<z_{\text{max}}$ for $z_{\text{max}}=0.7,1,2$ as showed in different colours. The quality of current data does not constrain these functions beyond $z\gtrsim 0.1$. However, it is interesting that our reconstruction permits the DE energy density $\rho_{\mathrm{de}}$ to become negative (thus, $\Omega_{\text{m}}>1$) for $z\gtrsim 0.7$. This is not impossible and, in principle, may be realized by DE being a scalar field with a negative potential tracking dust-like matter for energy densities much exceeding the present critical one. On the other hand, it is clearly seen from Fig. 5 that the change of $\rho_{\mathrm{de}}$ from the positive present value to a negative one, if it occurs at all, can happen at sufficiently low redshifts $z<2$ only (the cases B and C).

It is interesting to compare these obtained model-independent constraints on $(\Omega_{\text{m,0}},\sigma_{8,0})$ (Fig. 4 with those obtained in Shafieloo et al. (2018), where the approach was different. In that paper, the authors reconstructed $h(z)$ from the Pantheon SNIa compilation via iterative smoothing, and then obtained $\widehat{f\sigma_{8}}(z)$ by assuming a constant $\gamma=0.55$ (or varying it as a free, but constant, parameter), and obtain $f=\Omega_{\text{m}}^{\gamma}(z)$.

4 Summary & Conclusion

Using the latest compilation of RSD measurements, we reconstruct the growth history using two model-independent approaches, namely, crossing statistics and Gaussian processes, only assuming a flat-FLRW Universe and general relativity. We then used the method introduced by Starobinsky (1998) to reconstruct the expansion history, and fit the corresponding distance moduli to the Pantheon SNIa compilation, and finally obtained model-independent constraints on $(\Omega_{\text{m,0}},\sigma_{8,0})$. In addition, it is possible to reconstruct the dark energy equation of state $w(z)$ and the growth rate $\gamma(z)$. Applying the crossing statistics formalism, i.e., multiplying the best-fit $\Lambda\text{CDM}$+GR growth by some hyperfunction, and obtained constraints on $\Omega_{\text{m,0}},\sigma_{8,0},C_{i}$, where the $C_{i}$ are hyperparameters of the crossing hyperfunctions, we find consistency with $C_{0}=1,C_{i}=0,i\geq 1$, i.e., the data do not call for any modification to the best-fit. However, the preferred values for $\Omega_{\text{m,0}}=0.381^{+0.049}_{-0.113}$ and $\sigma_{8,0}=0.68^{+0.15}_{-0.12}$ are respectively higher and lower than the Planck best-fit, although consistent with them. Using Gaussian processes gives similar results. Both approaches suggest no statistically significant departure from $\Lambda\text{CDM}$+GR. On the other hand, it is interesting that they do not prohibit DE energy density from reaching zero at some recent redshift $0.7\lesssim z<2$ and becoming negative for larger $z$.

Future surveys such as the Dark Energy Spectroscopic Instrument are expected to provide more accurate measurements of the growth, and thus, to further constrain the gravity model and its parameters. This approach can be thought of as the reciprocal approach of Shafieloo et al. (2018), which uses direct reconstructions of $h$ from the supernovae data to fit the RSD data. Both approaches can be thought of as a mutual consistency test of the data and theory: in the former paper, the reconstructed expansion histories are tested against the growth data, while here, the reconstructed growth is compared to the SNIa data.

Our analyses point towards the consistency of the reconstructed $\Lambda\text{CDM}$ background evolution (via the mean function) with the growth history inside GR. While in GR, $\gamma=0.55$ is a very good approximation as long as $\Omega_{\text{m,0}}$ is not too small (Polarski et al. (2016)), it is not valid anymore beyond GR. On the contrary, the present approach, can be applied to non-GR models provided that $G_{\text{eff}}$ is known, or can even be used to reconstruct the effective Newton constant $G_{\text{eff}}$ (L’Huillier et al. in prep), and therefore to constraining modified gravity.

Acknowledgements

We thank the anonymous referee for suggesting to include negative DE energy densities in our study. BL thanks Ryan E. Keeley for valuable comments about Gaussian Processes. This work benefited from the Supercomputing Center/ Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2017-C2-0021) and the high performance computing clusters Polaris and Seondeok at the Korea Astronomy and Space Science Institute. BL would like to acknowledge the support of the National Research Foundation of Korea (NRF-2019R1I1A1A01063740). A.A.S. was partly supported by the program KP19-270 “Questions of the origin and evolution of the Universe” of the Presidium of the Russian Academy of Sciences.

Appendix A Validation on mock data

In order to validate the method, we applied it to a controlled simulated realization of the data, where the input cosmology is known. We generated mock data, following the redshift distribution and errors following the data in Shafieloo et al. (2018), assuming a known cosmology of $(\Omega_{\text{m,0}},\sigma_{8,0})=(0.3,0.8)$. We fit a flat-$\Lambda\text{CDM}$ universe to the total (RSD+SNIa) data, and use this best-fit $\Lambda\text{CDM}$ model as a mean function for the GP, and use its $\chi^{2}$ as a reference. Hereafter, we use subscript ${}_{\text{ref}}$ to denote the best-fit model. We applied our pipeline, reconstruct $f\sigma_{8}$ and $h$, and calculate the $\chi^{2}$ to the data (mock growth and Pantheon-like SNIa). Fig. 6 is the same as Fig. 3 but with our simulated data, and without separating cases A, B, and C. In addition to the best-fit in dashed lines, the true cosmology is shown in dotted lines. Fig. 7 shows the allowed contours for $(\Omega_{\text{m,0}},\sigma_{8,0})$, and the true cosmology, denoted by the green point, is inside the contours, showing the validity of the method.