1 Introduction

It is a great mystery since the current accelerated expansion of our universe was discovered in 1998 [1, 2]. More than 20 years passed, and we still do not know the very nature of the cosmic acceleration by now. Usually, an unknown energy component with negative pressure (dark energy) is introduced to interpret this mysterious phenomenon in general relativity (GR). Alternatively, one can make a modification to GR (modified gravity). In fact, modified gravity can also successfully explain the cosmic acceleration without invoking dark energy. So far, these two scenarios are both competent [3,4,5,6].

In order to understand the nature of the cosmic accelerated expansion, one of the most important tasks is to distinguish between the scenarios of dark energy and modified gravity. If the observational data could help us to confirm or exclude one of these two scenarios as the real cause of this mysterious phenomenon, it will be a great step forward. However, many cosmological observations merely probe the cosmic expansion history. Unfortunately, as is well known (see e.g. [7]), one can always build models sharing a same cosmic expansion history, and hence these models cannot be distinguished by using the observational data of the expansion history only. So, some independent and complementary probes are required. Later, it is proposed that if the cosmological models share a same cosmic expansion history, they might have different growth histories, which are characterized by the matter density contrast \(\delta (z)\equiv \delta \rho _m/\rho _m\) as a function of redshift z. Therefore, they might be distinguished from each other by combining the observations of both the expansion and growth histories (see e.g. [8,9,10,11,12,13,14,15,16,17,18, 113,114,115,116] and references therein).

It is convenient to introduce the growth rate \(f\equiv d\ln \delta /d\ln a\), where \(a=(1+z)^{-1}\) is the scale factor. As is well known, a good parameterization for the growth rate is given by [19,20,21,22,23]

$$\begin{aligned} f\equiv \frac{d\ln \delta }{d\ln a}=\Omega _m^\gamma , \end{aligned}$$
(1)

where \(\gamma \) is the growth index, and \(\Omega _m\) is the fractional energy density of matter. Beginning in e.g. [8, 9], it was advocated that the growth index \(\gamma \) is useful to distinguish the scenarios of dark energy and modified gravity. For example, it is found that \(\gamma =6/11\simeq 0.545\) for \(\Lambda \)CDM model [8, 9], and \(\gamma \simeq 0.55\) for most of dark energy models in GR [8]. In fact, they are clearly distinct from the ones of modified gravity theories. For instance, it is found that \(\gamma \simeq 0.68\) for Dvali–Gabadadze–Porrati (DGP) braneworld model [9, 15], and \(\gamma \simeq 0.42\) for most of viable f(R) theories [24,25,26,27]. In general, the growth index \(\gamma \) is a function of redshift z. It is argued that \(\gamma (z)\) lies in a relatively narrow range around the above values respectively, and hence one might distinguish between them.

In the literature, most of the relevant works assumed a particular cosmological model to obtain the growth index \(\gamma \). Thus, the corresponding results are model-dependent in fact. However, robust results should be model-independent. So, it is of interest to obtain the growth index \(\gamma \) from the observational data by using a model-independent approach. In fact, recently we have made an effort in [18] to obtain a non-parametric reconstruction of the growth index \(\gamma \) via Gaussian processes by using the latest observational data. Although the approach of Gaussian processes is clearly model-independent, its reliability at high redshift might be questionable. So, it is of interest to test the growth index \(\gamma \) by using a different method, and cross-check the corresponding results with the ones from Gaussian processes.

As is well known, one of the powerful model-independent approaches is cosmography [28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49, 118]. In fact, the only necessary assumption of cosmography is the cosmological principle. With cosmography, one can analyze the evolution of the universe without assuming any particular cosmological model. Essentially, cosmography is the Taylor series expansion of the quantities related to the cosmic expansion history (especially the luminosity distance \(d_L\)), and hence it is model-independent indeed. In the present work, we will constrain the growth index \(\gamma \) by using the latest observational data via the cosmographic approach. However, there are several shortcomings in the usual cosmography (see e.g. [48]). For instance, it is plagued with the problem of divergence or an unacceptably large error, and it fails to predict the future evolution of the universe. Thus, some generalizations of cosmography inspired by the Padé approximant were proposed in [48] (see also e.g. [50,51,52,53,54,55,56, 117, 119, 120]), which can avoid or at least alleviate the problems of ordinary cosmography. So, we also consider the Padé cosmography in this work.

The rest of this paper is organized as follows. In Sect. 2, we describe the methodology to constrain the growth index \(\gamma \) by using the latest observational data. In Sects. 3 and 4, we obtain the corresponding constraints on \(\gamma \) with the z-cosmography, the y-cosmography, and the Padé cosmography, respectively. In Sect. 5, conclusion and discussion are given.

2 Methodology

In the literature, there are many approaches to deal with the growth history. For example, one can consider a Lagrangian derived from an effective field theory (EFT) expansion [57, 58] (see also e.g. [59]), and implement the full background and perturbation equations for this action in the Boltzmann code EFTCAMB/EFTCosmoMC [60,61,62]. The second approach is more phenomenological [63,64,65,66,67,68,69] (see also e.g. [59, 70]), by directly parameterizing the functions of the gravitational potentials \(\Phi \) and \(\Psi \), such as \(\mu =G_\mathrm{eff}/G\), \(\eta =\Phi /\Psi \), and/or \(\Sigma \), Q, in the modified relativistic Poisson equations. It can be implemented by using the code MGCAMB [63, 69] integrated in CosmoMC [71]. The third approach is the simplest one, by directly parameterizing the growth rate f as in Eq. (1), with no need for numerically solving the perturbation equations. For simplicity, we choose this approach in the present work.

By definition \(f\equiv d\ln \delta /d\ln a\), it is easy to get (see e.g. [18, 72, 73])

$$\begin{aligned} \frac{\delta }{\delta _0} =\exp \left( \int _1^a \frac{f\,d\tilde{a}}{\tilde{a}}\right) = \exp \left( -\int _0^z\frac{f\,d\tilde{z}}{1+\tilde{z}}\right) , \end{aligned}$$
(2)

where the subscript “0” indicates the present value of the corresponding quantity, namely \(\delta _0=\delta (z=0)\). On the other hand, the cosmic expansion history can be characterized by the luminosity distance \(d_L=\left( c/H_0\right) D_L\), where c is the speed of light, \(H_0\) is the Hubble constant, and (see e.g. the textbooks [28, 29])

$$\begin{aligned} D_L\equiv (1+z)\int _0^z\frac{d\tilde{z}}{E(\tilde{z})}, \end{aligned}$$
(3)

in which \(E\equiv H/H_0\), and the Hubble parameter \(H\equiv \dot{a}/a\) (where a dot denotes the derivative with respect to cosmic time t). Note that we consider a flat Friedmann–Robertson–Walker (FRW) universe in this work. As is well known, E(z) is free of \(H_0\) actually. Differentiating Eq. (3), we get [49]

$$\begin{aligned} \frac{1+z}{E(z)} = \frac{d D_L}{dz}-\frac{D_L}{1+z}. \end{aligned}$$
(4)

If the luminosity distance \(d_L\) (or equivalently \(D_L\)) is known (in fact it will be given by the cosmography as below), we can obtain the dimensionless Hubble parameter E(z) by using Eq. (4). Then, the fractional energy density of matter is given by

$$\begin{aligned} \Omega _m(z)\equiv \frac{8\pi G\rho _m}{3H^2}= \frac{\Omega _{m0}(1+z)^3}{E^2(z)}. \end{aligned}$$
(5)

So, the growth rate \(f=\Omega _m^\gamma \) is on hand, and hence \(\delta /\delta _0\) in Eq. (2) is ready.

Table 1 The observational data of \(S_8\equiv \sigma _{8,\,0}(\Omega _{m0}/0.3)^{0.5}\). See the Refs. for details
Full size table

The data of the growth rate f can be obtained from redshift space distortion (RSD) measurements. In fact, the observational \(f_{obs}\) data have been used in some relevant works (e.g. [15, 74, 75]). However, it is sensitive to the bias parameter b which can vary in the range \(b\in [1,\,3]\). This makes the observational \(f_{obs}\) data unreliable [76]. Instead, the combination \(f\sigma _8(z)\equiv f(z)\,\sigma _8(z)\) is independent of the bias, and hence is more reliable, where \(\sigma _8(z)=\sigma _8(z=0) \,\delta (z)/\delta _0=\sigma _{8,\,0}\,\delta (z)/\delta _0\) is the redshift-dependent rms fluctuations of the linear density field within spheres of radius \(8h^{-1}\)Mpc [76]. In fact, the observational \(f\sigma _{8,\,obs}\) data can be obtained from weak lensing or RSD measurements [76, 77]. In the present work, we use the sample consisting of 63 observational \(f\sigma _{8,\,obs}\) data published in [77], which is the largest \(f\sigma _8\) compilation in the literature by now. As mentioned above, once \(D_L\) is given, we can get the theoretical \(f\sigma _8\) by using Eqs. (4), (5), and (1), (2). Thus, the \(\chi ^2\) from the \(f\sigma _8\) data is given by

$$\begin{aligned} \chi ^2_{f\sigma _8}=\sum _i \frac{\left[ \,f\sigma _{8,\,obs}(z_i) -f\sigma _{8,\,\mathrm{mod}}(z_i)\,\right] ^2}{\sigma _{f\sigma _8}^2 (z_i)}. \end{aligned}$$
(6)

It is easy to see that only using the observational \(f\sigma _8\) data is not enough to constrain the model parameter \(\Omega _{m0}\), and the cosmographic parameters \(q_0\), \(j_0 \ldots \) in \(D_L\). Since they mainly affect the cosmic expansion history, we also use such kinds of observations. Obviously, the type Ia supernovae (SNIa) data is useful. Here, we consider the Pantheon sample [78,79,80] consisting of 1048 SNIa, which is the largest spectroscopically confirmed SNIa sample by now. The corresponding \(\chi ^2\) is given by

$$\begin{aligned} \chi ^2_\mathrm{Pan}=\Delta \varvec{m}^{\,T}\cdot \varvec{C}^{-1}\cdot \Delta \varvec{m}, \end{aligned}$$
(7)

where for the ith SNIa, \(\Delta m_i=m_i-m_{\mathrm{mod},\,i}\,\), and \(\varvec{C}\) is the total covariance matrix,

$$\begin{aligned} m_\mathrm{mod} = 5\log _{10} D_L + \mathcal{M}, \end{aligned}$$
(8)

in which \(\mathcal M\) is a nuisance parameter corresponding to some combination of the absolute magnitude M and \(H_0\). We refer to [78,79,80] for technical details (see also e.g. [81]). Since \(H_0\) is absorbed into \(\mathcal M\) in the analytic marginalization, the Pantheon SNIa sample is free of the Hubble constant \(H_0\).

We further consider the observational data from the baryon acoustic oscillation (BAO). Note that there exist many kinds of BAO data in the literature, such as \(D_V(z)\), \(d_z\equiv r_s(z_d)/D_V(z)\), \(D_A(z)/r_d\), \(D_M(z)/r_d\), \(H(z)\,r_s(z_d)\) and A. However, in the former ones, they will introduce one or more extra model parameters, for instance \(H_0\), and/or \(\Omega _b h^2\). Since the \(f\sigma _8\) data, the SNIa data, the cosmography for \(D_L\), and other data are all free of \(H_0\) and \(\Omega _b h^2\), we choose to avoid introducing extra model parameters here. Thus, in this work, we use the BAO data only in the form of (see e.g. [82, 83])

$$\begin{aligned} A\equiv \Omega _{m0}^{1/2}H_0 D_V/(cz) =\frac{\Omega _{m0}^{1/2}}{z}\left[ \frac{D_L^2}{(1+z)^2}\cdot \frac{z}{E(z)}\,\right] ^{1/3}, \end{aligned}$$
(9)

which does not introduce extra model parameters since the factor \(c/H_0\) in \(D_V\) is canceled. We consider the six data of the acoustic parameter A(z) compiled in the last column of Table 3 of [83]. The first data point from 6dFGS is uncorrelated with other five ones, and hence its \(\chi ^2_\mathrm{6dFGS}=(A_{obs}-A_\mathrm{mod})^2/\sigma ^2\) directly. The 2nd and 3rd data points from SDSS are correlated with coefficient 0.337, and hence the inverse covariance matrix of these two data points is given by

$$\begin{aligned} \varvec{C}^{-1}_\mathrm{SDSS}= \left( \begin{array}{cc} 4406.72 &{} -1485.06 \\ -1485.06 &{} 4406.72 \end{array} \right) . \end{aligned}$$
(10)

The inverse covariance matrix of the last three data points from WiggleZ is given in Table 2 of [83],

$$\begin{aligned} \varvec{C}^{-1}_\mathrm{WiggleZ}= \left( \begin{array}{ccc} 1040.3 &{} -807.5 &{} 336.8 \\ -807.5 &{} 3720.3 &{} -1551.9 \\ 336.8 &{} -1551.9 &{} 2914.9 \end{array} \right) . \end{aligned}$$
(11)

The \(\chi ^2\) from the data of SDSS and WiggleZ are both given in the form of \(\chi ^2=\Delta \varvec{A}^{\,T}\cdot \varvec{C}^{-1}\cdot \Delta \varvec{A}\). Thus, the total \(\chi ^2\) from the BAO data is \(\chi ^2_\mathrm{BAO}= \chi ^2_\mathrm{6dFGS}+\chi ^2_\mathrm{SDSS}+\chi ^2_\mathrm{WiggleZ}\).

Table 2 The mean with \(1\sigma \), \(2\sigma \), \(3\sigma \) marginalized uncertainties of the model parameters for the cases with the z-cosmography and \(\gamma =\gamma _0\) (labeled as “ z-0 ”), \(\gamma =\gamma _0+\gamma _1\,z\) (labeled as “ z-1 ”), \(\gamma =\gamma _0+\gamma _1\,z+\gamma _2\,z^2\) (labeled as “ z-2 ”). See the text for details
Full size table

On the other hand, the free parameter \(\sigma _{8,\,0}\) cannot be well constrained by using the \(f\sigma _8\) data and the observations of the expansion history. Fortunately, in the literature there are many observational data of the combination \(S_8\equiv \sigma _{8,\,0}(\Omega _{m0}/0.3)^{0.5}\) from the cosmic shear observations [84], which can be used to constrain both the free parameters \(\sigma _{8,\,0}\) and \(\Omega _{m0}\). Here, we consider the ten data points given in Table 1. The corresponding \(\chi ^2_{S_8}= \sum _i \, (S_{8,\,obs,\,i}-S_{8,\,\mathrm{mod},\,i})^2 / \sigma _{S_8,\,i}^2\,\). Note that if the upper and the lower uncertainties of the data are not equal, we choose the bigger one as \(\sigma _{S_8,\,i}\) conservatively.

In fact, there are other kinds of observational data in the literature. However, we do not use them here, to avoid introducing extra model parameters, as mentioned above. For instance, if we want to use the 51 observation H(z) data compiled in [95] (the largest sample by now to our best knowledge), an extra free parameter \(H_0\) is necessary. So, we give up. On the other hand, since the usual cosmography cannot work well at very high redshift, we also do not consider the observational data from cosmic microwave background (CMB) at redshift \(z\sim 1090\). Otherwise, the cosmographic parameters should be fine-tuned. However, the Padé cosmography works well at very high redshift, and hence we can use the CMB data in this case (see Sect. 4).

All the model parameters can be constrained by using the observational data to perform a \(\chi ^2\) statistics. Here, the total \(\chi ^2_\mathrm{tot}=\chi ^2_{f\sigma _8}+\chi ^2_\mathrm{Pan}+ \chi ^2_\mathrm{BAO}+\chi ^2_{S_8}\). In the following, we use the Markov Chain Monte Carlo (MCMC) code CosmoMC [71] to this end.

3 Observational constraints with the ordinary cosmography

3.1 The case of z-cosmography

At first, we consider the case of z-cosmography. Introducing the so-called cosmographic parameters, namely the Hubble constant \(H_0\), the deceleration \(q_0\), the jerk \(j_0\), the snap \(s_0 \ldots \),

$$\begin{aligned}&H_0\equiv \left. \frac{1}{a}\frac{da}{dt}\right| _{t=t_0},\quad q_0\equiv \left. -\frac{1}{aH^2}\frac{d^{2}a}{dt^{2}}\right| _{t=t_0},\nonumber \\&j_0\equiv \left. \frac{1}{aH^3}\frac{d^{3}a}{dt^{3}}\right| _{t=t_0},\quad s_0\equiv \left. \frac{1}{aH^4}\frac{d^{4}a}{dt^{4}}\right| _{t=t_0}, \ldots \end{aligned}$$
(12)

one can express the quantities related to the cosmic expansion history, e.g. the scale factor a(t), the Hubble parameter H(z), and the luminosity distance \(d_L(z)\), as a Taylor series expansion (see e.g. [28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49] and references therein). The most important one is the luminosity distance \(d_L(z)\), and its Taylor series expansion with respect to redshift z reads (see e.g. [28,29,30,31,32, 48, 49] for details)

$$\begin{aligned} d_L(z)= & {} \frac{cz}{H_0}\left[ 1+\frac{1}{2}\left( 1-q_0\right) z -\frac{1}{6}\left( 1-q_0-3q_0^2+j_0\right) z^2\right. \nonumber \\&\left. +\frac{1}{24}\left( 2-2q_0-15q_0^2-15q_0^3\right. \right. \nonumber \\&\left. \left. +5j_0 +10q_0 j_0+s_0\right) z^3+\mathcal{O}\left( z^4\right) \right] . \end{aligned}$$
(13)

Since the constraints become loose if the number of free parameters increases, we only consider the cosmography up to third order. Thus, the dimensionless luminosity distance \(D_L=H_0 d_L/c\) is given by

$$\begin{aligned} D_L(z)= & {} z+\frac{1}{2}\left( 1-q_0\right) z^2 -\frac{1}{6}\left( 1-q_0-3q_0^2+j_0\right) z^3\nonumber \\&+\mathcal{O}(z^4), \end{aligned}$$
(14)

in which only two free cosmographic parameters \(q_0\) and \(j_0\) are involved. Note that the Hubble constant \(H_0\) does not appear, since the factor \(c/H_0\) in \(d_L\) is canceled.

In the literature, the growth index \(\gamma \) is often assumed to be constant (see e.g. [19, 20, 70, 96, 97]). However, in general it is varying as a function of redshift z. To be model-independent, we can also expand \(\gamma (z)\) as a Taylor series with respect to redshift z, namely \(\gamma (z)=\gamma _0+\gamma _1\,z+\gamma _2\,z^2+\cdots \), where the coefficients \(\gamma _0\), \(\gamma _1\), \(\gamma _2 \ldots \) are constants. Here, we consider three cases, labeled as “ z-0 ”, “ z-1 ”, “ z-2 ”, in which \(\gamma (z)\) is Taylor expanded up to zeroth, first, second orders, respectively.

Substituting Eq. (14) into Eq. (4), we can get the dimensionless Hubble parameter E(z). Using Eqs. (5), (1), (2), and \(\gamma \), we obtain \(f=\Omega _m^\gamma \) and then \(f\sigma _8\). Substituting \(D_L(z)\) and E(z) into Eqs. (8) and (9), we find \(m_\mathrm{mod}\) and \(A_\mathrm{mod}\). Finally, the total \(\chi ^2_\mathrm{tot}\) is ready.

By using the latest observational data, we obtain the constraints on all the model parameters involved, and present them in Table 2, for the z-0, z-1, z-2 cases. Since we mainly concern the parameters related to the growth index \(\gamma \), namely \(\gamma _0\), \(\gamma _1\) and \(\gamma _2\), we also present their 1D marginalized probability distributions in Fig. 1. Obviously, in all cases, \(q_0<0\) and \(j_0>0\) far beyond \(3\sigma \) confidence level (C.L.), and these mean that today the universe is accelerating, while the acceleration is still increasing. From Tabel 2 and Fig. 1, it is easy to see that for all cases, \(\gamma _0\simeq 0.42\) is inconsistent with the latest observational data far beyond \(3\sigma \) C.L. Note that \(\gamma _0\simeq 0.55\) is consistent with the latest observational data within the \(1\sigma \) region for the z-2 case, but is inconsistent beyond \(2\sigma \) C.L. for both the z-0 and z-1 cases (while it is still consistent within the \(3\sigma \) region). On the other hand, a varying \(\gamma \) with non-zero \(\gamma _1\) and/or \(\gamma _2\) is favored. In the linear case with \(\gamma =\gamma _0+\gamma _1\,z\) (namely the z-1 case), \(\gamma _1<0\) in the \(1\sigma \) region, and hence the growth index \(\gamma \) decreases as redshift z increases. In the quadratic case with \(\gamma =\gamma _0+\gamma _1\,z+ \gamma _2\,z^2\) (namely the z-2 case), \(\gamma _2<0\) beyond \(3\sigma \) C.L., and hence the function \(\gamma (z)\) is a parabola opening down, namely \(\gamma \) increases and then decreases as redshift z increases. There exists an arched structure in the moderate redshift range. This result is quite similar to the one of [18].

Fig. 1
figure 1

The 1D marginalized probability distributions of the parameters related to \(\gamma \). The \(1\sigma \), \(2\sigma \), \(3\sigma \) uncertainties are shown by the green, blue, red regions, respectively. The top, middle, bottom panels correspond to the z-0, z-1, z-2 cases, respectively. See the text and Table 2 for details

Full size image

3.2 The case of y-cosmography

Let us turn to the case of y-cosmography. As is well known, a Taylor series with respect to redshift z converges only at low redshift z around 0, and it might diverge at high redshift \(z>1\). In the literature (see e.g. [33,34,35,36,37, 39, 48, 49]), a popular alternative to the z-cosmography is replacing z with the so-called y-shift, \(y=1-a=z/(1+z)\). Obviously, \(y<1\) holds in the whole cosmic past \(0\le z<\infty \), and hence the Taylor series with respect to y-shift converges. In this case, we can expand the dimensionless luminosity distance \(D_L=H_0 d_L/c\) as a Taylor series with respect to y (see e.g. [32, 34, 48, 49] for details),

$$\begin{aligned} D_L(y)= & {} y+\frac{1}{2}\left( 3-q_0\right) y^2+ \frac{1}{6}(11-5q_0+3q_0^2-j_0)y^3\nonumber \\&+\mathcal{O}(y^4), \end{aligned}$$
(15)

in which only two free cosmographic parameters \(q_0\) and \(j_0\) are involved, since we only consider the cosmography up to third order in this work as mentioned above. Accordingly, here we also expand the growth index \(\gamma \) as a Taylor series with respect to y, namely \(\gamma (y)=\gamma _0+\gamma _1\,y+\gamma _2\,y^2+\cdots \). Similarly, we consider three cases, labeled as “ y-0 ”, “ y-1 ”, “ y-2 ”, in which \(\gamma (y)\) is Taylor expanded up to zeroth, first, second orders, respectively. Noting \(y=z/(1+z)\) and \(dF/dz=(1+z)^{-2}dF/dy\) for any function F, the formalism in Sect. 2 is still valid in the case of y-cosmography.

By using the latest observational data, we obtain the constraints on all the model parameters involved, and present them in Table 3, for the y-0, y-1, y-2 cases. In Fig. 2, we also present the 1D marginalized probability distributions of the parameters related to the growth index \(\gamma \), namely \(\gamma _0\), \(\gamma _1\) and \(\gamma _2\). Obviously, the y-0 case is fairly different from the y-1, y-2 cases. In fact, beyond \(3\sigma \) C.L., and \(q_0>0\) also beyond \(3\sigma \) C.L. in the y-0 case. The unusual result that the universe is decelerating (\(q_0>0\)) suggests that the y-0 case with a constant \(\gamma =\gamma _0\) is not competent to describe the real universe, and consequently \(\gamma \) should be varying instead. This conclusion is also supported by the abnormal \(\chi ^2_{min}=1281.5972\) of the y-0 case, which is significantly larger than the ones of the y-1, y-2 cases (see Tabel 5). In both the y-1, y-2 cases, \(q_0<0\) beyond \(3\sigma \) C.L., and this means that the universe is undergoing an acceleration. On the other hand, \(\gamma _0\simeq 0.42\) is inconsistent with the latest observational data far beyond \(3\sigma \) C.L. in both the y-1, y-2 cases. \(\gamma _0\simeq 0.55\) is well consistent with the latest observational data within \(1\sigma \) region in the y-1 case, but it is inconsistent with the latest observational data beyond \(3\sigma \) C.L. in the y-2 case. A varying \(\gamma \) with non-zero \(\gamma _1\) and/or \(\gamma _2\) is favored. It is easy to see that \(\gamma _1<0\) far beyond \(3\sigma \) C.L. in both the y-1, y-2 cases, and \(\gamma _2>0\) far beyond \(3\sigma \) C.L. in the y-2 case. However, \(\gamma =\gamma (y)=\gamma (z/(1+z))\), and hence one should be careful to treat \(\gamma \) as a function of redshift z.

Table 3 The mean with \(1\sigma \), \(2\sigma \), \(3\sigma \) marginalized uncertainties of the model parameters for the cases with the y-cosmography and \(\gamma =\gamma _0\) (labeled as “ y-0 ”), \(\gamma =\gamma _0+\gamma _1\,y\) (labeled as “ y-1 ”), \(\gamma =\gamma _0+\gamma _1\,y+\gamma _2\,y^2\) (labeled as “ y-2 ”). See the text for details
Full size table
Fig. 2
figure 2

The same as in Fig. 1, except for the y-0, y-1, y-2 cases. See the text and Table 3 for details

Full size image

4 Observational constraints with the Padé cosmography

In the previous section, two types of ordinary cosmography are considered. As mentioned above, the z-cosmography might diverge at high redshift z. So, the y-cosmography was proposed as an alternative in the literature, which converges in the whole cosmic past \(0\le z<\infty \). However, there still exist several problems in the y-cosmography. In practice, the Taylor series should be truncated by throwing away the higher order terms, since it is difficult to deal with infinite series. So, the error of a Taylor approximation with lower order terms will become unacceptably large when \(y=z/(1+z)\) is close to 1 (say, when \(z>9\)). On the other hand, the y-cosmography cannot work well in the cosmic future \(-1<z<0\). The Taylor series with respect to \(y=z/(1+z)\) does not converge when \(y<-1\) (namely \(z<-1/2\)), and it drastically diverges when \(z\rightarrow -1\) (it is easy to see that \(y\rightarrow -\infty \) in this case). Therefore, in [48], we proposed some generalizations of cosmography inspired by the Padé approximant, which can avoid or at least alleviate the problems of ordinary cosmography.

The so-called Padé approximant can be regarded as a generalization of the Taylor series. For any function F(x), its Padé approximant of order \((m,\,n)\) is given by the rational function [98,99,100,101,102] (see also e.g. [50,51,52,53,54,55,56])

$$\begin{aligned} F(x)=\frac{\alpha _0+\alpha _1 x+\cdots +\alpha _m x^m}{1+ \beta _1 x+\cdots +\beta _n x^n}, \end{aligned}$$
(16)

where m and n are both non-negative integers, and \(\alpha _i\), \(\beta _i\) are all constants. Obviously, it reduces to the Taylor series when all \(\beta _i=0\). Actually in mathematics, a Padé approximant is the best approximation of a function by a rational function of given order [101]. In fact, the Padé approximant often gives a better approximation of the function than truncating its Taylor series, and it may still work where the Taylor series does not converge [101].

One can directly parameterize the dimensionless luminosity distance based on the Padé approximant with respect to redshift z [48],

$$\begin{aligned} D_L=\frac{H_0 d_L}{c}=\frac{\alpha _0+\alpha _1 z +\cdots +\alpha _m z^m}{1+\beta _1 z+\cdots +\beta _n z^n}. \end{aligned}$$
(17)

Following [48], we consider a moderate order \((2,\,2)\) in this work, and then

$$\begin{aligned} D_L\equiv \frac{H_0 d_L}{c}=\frac{\alpha _0+\alpha _1 z+\alpha _2 z^2}{1+\beta _1 z+\beta _2 z^2}. \end{aligned}$$
(18)

Obviously, it can work well in the whole redshift range \(-1<z<\infty \), including not only the past but also the future of the universe. In particular, it is still finite even when \(z\gg 1\). In fact, this \(D_L\) was confronted with Union2.1 SNIa data and Planck 2015 CMB data in [48], and the parameters \(\alpha _0\) and \(\beta _2\) were found to be very close to 0 even in the \(3\sigma \) region. So, in the present work, it is safe to directly set

$$\begin{aligned} \alpha _0=\beta _2=0, \end{aligned}$$
(19)

and then the free parameters are now \(\alpha _1\), \(\alpha _2\) and \(\beta _1\). Note that in fact \(\alpha _0=0\) is required by \(d_L(z=0)=0\) theoretically. On the other hand, we can also expand \(\gamma (z)\) as a Taylor series with respect to redshift z, namely \(\gamma (z)=\gamma _0+\gamma _1\,z+\gamma _2\,z^2+...\). Again, we consider three cases, labeled as “ P-0 ”, “ P-1 ”, “ P-2 ”, in which \(\gamma (z)\) is Taylor expanded up to zeroth, first, second orders, respectively.

Table 4 The mean with \(1\sigma \), \(2\sigma \), \(3\sigma \) marginalized uncertainties of the model parameters for the cases with the Padé cosmography and \(\gamma =\gamma _0\) (labeled as “ P-0 ”), \(\gamma =\gamma _0+\gamma _1\,z\) (labeled as “ P-1 ”), \(\gamma =\gamma _0+\gamma _1\,z+\gamma _2\,z^2\) (labeled as “ P-2 ”). See the text for details
Full size table

Since the Padé cosmography still works well at very high redshift \(z\gg 1\) in contrast to the ordinary cosmography as mentioned above, in this section, we further use the latest CMB data in addition to the observational data mentioned in Sect. 2. However, using the full data of CMB to perform a global fitting consumes a large amount of computation time and power. As a good alternative, one can instead use the shift parameter R [103, 104] from CMB, which has been used extensively in the literature (including the works of the Planck and the WMAP Collaborations themselves). It is argued in e.g. [105,106,107] that the shift parameter R is model-independent and contains the main information of the observation of CMB. As is well known, the shift parameter R is defined by [103,104,105,106,107]

$$\begin{aligned} R\equiv \sqrt{\Omega _{m0}H^2_0}\,\left( 1+z_*\right) d_A(z_*)/c =\frac{\sqrt{\Omega _{m0}}\, D_L(z_*)}{1+z_*}, \end{aligned}$$
(20)
Fig. 3
figure 3

The same as in Fig. 1, except for the P-0, P-1, P-2 cases. See the text and Table 4 for details

Full size image

where the redshift of the recombination \(z_*=1089.92\) from the Planck 2018 data [94], and the angular diameter distance \(d_A\) is related to the luminosity distance \(d_L\) through \(d_A=d_L(1+z)^{-2}\) (see e.g. the textbooks [28, 29]). Here we adopt the value \(R_{obs}=1.7502\pm 0.0046\) [108] derived from the Planck 2018 data. Thus, the corresponding \(\chi ^2\) from the latest CMB data is given by \(\chi ^2_R=(R_\mathrm{mod}-R_{obs})^2/\sigma _R^2\). Although the number of data points \(\mathcal N\) and the number of free parameters \(\kappa \) both increase by 1, the degree of freedom \(dof=\mathcal{N}-\kappa \) is unchanged in this case. It is worth noting that the acoustic scale \(l_A\), and \(\Omega _b h^2\), the scalar spectral index \(n_s\) are commonly used with the shift parameter R in the literature, but they will introduce extra model parameters as mentioned above, and hence we do not use them here.

By using the latest observational data, we obtain the constraints on all the model parameters involved, and present them in Table 4, for the P-0, P-1, P-2 cases. In Fig. 3, we also present the 1D marginalized probability distributions of the parameters related to the growth index \(\gamma \), namely \(\gamma _0\), \(\gamma _1\) and \(\gamma _2\). In all cases, \(\gamma _0\simeq 0.42\) is inconsistent with the latest observational data beyond \(2\sigma \) C.L. (but it could be consistent in the \(3\sigma \) region). On the other hand, in all cases, \(\gamma _0\simeq 0.55\) is well consistent with the latest observational data within the \(1\sigma \) region. Note that in all cases, a constant \(\gamma =\gamma _0\) (namely \(\gamma _1=0\) and \(\gamma _2=0\)) is well consistent with the latest observational data within the \(1\sigma \) region (but see below).

5 Conclusion and discussion

In this work, we consider the constraints on the growth index \(\gamma \) by using the latest observational data. To be model-independent, we use cosmography to describe the cosmic expansion history, and also expand the general \(\gamma (z)\) as a Taylor series with respect to redshift z or y-shift, \(y=1-a=z/(1+z)\). We find that the present value \(\gamma _0=\gamma (z=0)\simeq 0.42\) (for most of viable f(R) theories) is inconsistent with the latest observational data beyond \(3\sigma \) C.L. in the six cases with the usual cosmography, or beyond \(2\sigma \) C.L. in the three cases with the Padé cosmography. This result supports our previous work [18]. On the other hand, \(\gamma _0\simeq 0.55\) (for dark energy models in GR) is consistent with the latest observational data at \(1\sigma \) C.L. in five of the nine cases under consideration, but is inconsistent beyond \(2\sigma \) C.L. in the other four cases (while it is still consistent within the \(3\sigma \) region). Therefore, we can say nothing firmly about \(\gamma _0\simeq 0.55\). This result is still consistent with the reconstructed \(\gamma (z)\) at \(z=0\) obtained in our previous work [18]. A varying \(\gamma \) with non-zero \(\gamma _1\) and/or \(\gamma _2\) is favored in the cases with the usual cosmography, while in the cases with the Padé cosmography, a constant \(\gamma =\gamma _0\) (namely \(\gamma _1=0\) and \(\gamma _2=0\)) can still be consistent with the latest observational data (but this might be artificial, see below).

Table 5 Comparing the nine cases considered in the present work. See the text for details and caution
Full size table

It is of interest to compare the 9 cases considered here. We adopt several goodness-of-fit criteria used extensively in the literature to this end, such as \(\chi ^2_{min}/dof\), \(P(\chi ^2>\chi ^2_{min})\) (see e.g. [109, 110]), Bayesian Information Criterion (BIC) [111] and Akaike Information Criterion (AIC) [112], where the degree of freedom \(dof=\mathcal{N}-\kappa \), while \(\mathcal N\) and \(\kappa \) are the number of data points and the number of free model parameters, respectively. The BIC is defined by [111]

$$\begin{aligned} \mathrm{BIC}=-2\ln \mathcal{L}_{max}+\kappa \ln \mathcal{N}, \end{aligned}$$
(21)

and the AIC is defined by [112]

$$\begin{aligned} \mathrm{AIC}=-2\ln \mathcal{L}_{max}+2\kappa , \end{aligned}$$
(22)

where \(\mathcal{L}_{max}\) is the maximum likelihood. In the Gaussian cases, \(\chi ^2_{min}=-2\ln \mathcal{L}_{max}\). The difference in BIC or AIC between two models makes sense. We choose the P-0 case to be the fiducial model when we calculate \(\Delta \)BIC and \(\Delta \)AIC. In Table 5, we present \(\chi ^2_{min}/dof\), \(P(\chi ^2>\chi ^2_{min})\), \(\Delta \)BIC and \(\Delta \)AIC for the nine cases considered in this work. Clearly, the cases with y-cosmography are the worst, while the cases P-0 and z-2 are the best. In fact, the goodness-of-fit criteria for the cases P-0 and z-2 are fairly close. A caution should be mentioned here. All the criteria given in Table 5 are based on \(\chi ^2_{min}\), which are read from the output .likestats files of the CosmoMC program GetDist. However, as the CosmoMC [71] readme file puts it, “ file\(_{-}\)root.likestats gives the best fit sample model, its likelihood, and\(\ldots \) Note that MCMC does not generally provide accurate values for the best-fit.” Keeping this in mind, we could say that the cases P-0 and z-2 are equally good, since their not so accurate \(\chi ^2_{min}\) are very close actually. In the P-0 case, the growth index \(\gamma =\gamma _0\) is constant. However, in the z-2 case, \(\gamma _2<0\) beyond \(3\sigma \) C.L. (see Table 2 and Fig. 1), and hence the function \(\gamma (z)\) is a parabola opening down, namely \(\gamma \) increases and then decreases as redshift z increases. There exists an arched structure in the moderate redshift range. This result is quite similar to the one of [18]. In Fig. 4, we show a demonstration of \(\gamma =\gamma _0+\gamma _1\,z+\gamma _2\,z^2\) with \(\gamma _0=0.6\), \(\gamma _1=0.45\), \(\gamma _2=-0.24\), which are all well within the \(1\sigma \) regions of their observational constraints for the z-2 case (see Tabel 2 and Fig. 1).

It is worth noting that throughout this work, we always consider the growth index \(\gamma \) as a Taylor series with respect to z or y, namely \(\gamma (z)=\gamma _0+\gamma _1\,z+\gamma _2\,z^2+\cdots \), or \(\gamma (y)=\gamma _0+\gamma _1\,y+\gamma _2\,y^2+\cdots \). However, in Sect. 4, we parameterize the dimensionless luminosity distance \(D_L\) by using the Padé approximant, and hence it can still work well at very high redshift \(z\sim 1090\). Obviously, it is better to also parameterize the growth index \(\gamma (z)\) by using the Padé approximant (we thank the referee for pointing out this issue). But the cost is expensive to do this. If we want to catch the arched structure in \(\gamma (z)\), at least a Padé approximant of order \((2,\,2)\) is needed, which has 5 free parameters (n.b. Eq. (16)), and almost double the number of free parameters in a 2nd order Taylor series. So, in the P-2 case, the total number of free model parameters will be ten. It will consume significantly more computation power and time, but the corresponding constraints will be very loose. Therefore, we choose not to do this at a great cost. But one should be aware of the possible artificial results from this choice. For example, \(\gamma (z)=\gamma _0+\gamma _1\,z+\gamma _2\,z^2\) will diverge at \(z\sim 1090\), and hence the values of \(\gamma _1\) and \(\gamma _2\) tend to be zero to fit the high-z CMB data in the P-1, P-2 cases (we thank the referee for pointing out this issue).

Fig. 4
figure 4

A demonstration of \(\gamma =\gamma _0+\gamma _1\,z+\gamma _2\,z^2\) with \(\gamma _0=0.6\), \(\gamma _1=0.45\), \(\gamma _2=-0.24\), which are all well within the \(1\sigma \) regions of their observational constraints for the z-2 case (see Tabel 2 and Fig. 1)

Full size image

Some remarks are in order. First, the growth rate f and then the growth index \(\gamma \) for modified gravity scenarios (especially f(R) theories) in principle are not only time-dependent but also scale-dependent (see e.g. [24, 25]). However, as is shown in e.g. [24, 25], the behavior of \(\gamma \) is nearly scale-independent at low redshift , and \(\gamma _0=\gamma (z=0)\) is also nearly independent of scale. So, this issue does not change the main conclusions of the present work, although it may be studied carefully in the future work. Second, as is mentioned in the beginning of Sect. 2, there exist other two approaches dealing with the growth history, which numerically solve the perturbation equations by using the code CAMB integrated in CosmoMC. We will also consider these approaches in the future work. Third, in the present work, we do not use some types of observational data (for example, the observational H(z) data, and other kinds of BAO data) to avoid introducing extra model parameters. However, in principle, it is not terrible to do so, although the constraints might be loose and the calculations might be complicated. Finally, in this work, we only consider the Taylor series expansion of the growth index \(\gamma \) up to 2nd order, and the usual cosmography up to 3rd order. In fact, one can also further consider higher orders in these cases. We anticipate that our main conclusions will not change significantly.