Non-parametric reconstruction of dark energy and cosmic expansion from the Pantheon compilation of type Ia supernovae

The late time cosmic acceleration is one of the most important discoveries in modern cosmology and it revives Einstein's cosmological constant hypothesis. Since the first discovery of cosmic acceleration from type Ia supernovae (SNe Ia) in the late 1990s [1, 2], it has now been confirmed by various other independent observations, such as the large scale structure [3], growth function [4] and cosmic microwave background radiations [5, 6]. This leads to the final foundation of the standard model of cosmology, i.e. the cold dark matter plus a cosmological constant ($ \Lambda {\rm{CDM}}$) model. Here, $ \Lambda $ denotes the cosmological constant (or an alternative name dark energy), which provides a negative pressure and is responsible for the acceleration of the universe. According to the $ \Lambda {\rm{CDM}}$ model, the equation of state (EoS) of dark energy is a constant and does not evolve with cosmos, i.e., $ w = p/\rho\equiv -1 $. Although the $ \Lambda {\rm{CDM}}$ model has achieved great success, it still confronts some problems, among which the most important ones are the coincidence problem and the fine-turning problem [7, 8]. The EoS of dark energy plays an essential role in the evolution of the universe. In terms of different EoS of dark energy, several alternative models have been proposed, such as the model with a constant w but does not necessarily equate to −1, the evolving dark energy models, e.g., the Chevallier-Polarski-Linder parametrization [9, 10] and various other parameterizations [11, 12]. In some models the dark energy is replaced by a scalar field, such as the quintessence field [13], phantom field [14] and tachyon field [15]. The effective EoS of these scalar fields is also evolving with the cosmos.

The above models depend on the specific parametrization of dark energy or the scalar fields; thus, they are model dependent. Moreover, most parameterizations lack physical interpretation; hence, they are just phenomenological. Because we have no prior knowledge on the explicit form of dark energy, reconstructing it in a non-parametric manner is of great importance. Therefore, some model-independent methods have been proposed, among which the Gaussian processes (GP) is one of the most widely used methods. Unlike the best-fitting method that must have a concrete model to fit the data, the GP method can reconstruct a theoretical curve from the discrete data points without evolving any specific model. Since Ref. [16] first applied the GP method to investigate the dark energy, it has been widely used and has shown its applicability in cosmology [17–26]. The advantage of the GP method is that it does not need the concrete form of the model; the only assumption is that the observational data points are drawn from the multivariate Gaussian distribution.

In this study, we attempt to reconstruct the EoS of dark energy using the GP method from the latest dataset of SNe Ia, i.e., the Pantheon compilation [27], which consist of 1048 finely calibrated SNe Ia. The EoS of dark energy has strong influence on the Hubble expansion rate $ H(z) $ and the deceleration parameter $ q(z) $, which will be simultaneously obtained in the reconstruction. The reconstructed dark energy depends on the second-order derivative of the GP (see the next section for details); hence, it has a large uncertainty. To improve the significance, we use the direct measurement of Hubble parameters $ H(z) $ as an additional constraint in the GP reconstruction. Because tensin is greater than $ 3\sigma $ is the between the Hubble constant $ H_0 $ from the local distance ladders [28] and from the global CMB radiation [6], we will also investigate the impact of different values of $ H_0 $ on the reconstruction.

The remainder of this paper is organized as follows. In section 2, we introduce the methodology and the relevant data used in our analysis. The results along with some discussions are presented in section 3. Finally, a short summary is given in section 4.

The Hubble expansion rate $ H(z) $ strongly depends on the contents of the universe and the EoS of dark energy. In a spatially flat Friedmann-Robertson-Walker universe dominated by non-relativistic matter (including baryons and dark matter) and dark energy, the evolution of $ H(z) $ is governed by the Friedmann equation [29]

where $ H_0 $ is the Hubble constant, $ \Omega_M $ and $ \Omega_{\Lambda} $ are the normalized density of non-relativistic matter and dark energy at current epoch, respectively, and $ w(z) = p(z)/\rho(z) $ is the EoS of dark energy. The normalized co-moving distance is related to the Hubble expansion rate by [30]

where $ E(z)\equiv H(z)/H_0 $ is the normalized Hubble parameter. From Eq. (2) we have

where the prime denotes the derivative with respect to redshift z. Combining Eqs. (1) and (3), we can solve for $ w(z) $,

The acceleration of the universe is often represented by the so-called "deceleration parameter", which is defined by $ q(z) = -\ddot{a}a/\dot{a}^2 $, where $ a = 1/(1+z) $ is the scale factor of the universe, and the dot denotes the derivative with respect to cosmic time. Using $ H = \dot{a}/a $, the deceleration parameter can be rewritten as

A positive or negative q implies a decelerating or accelerating universe, respectively.

If we know $ d_c(z) $ as a function of z, $ E(z) $ can be obtained from Eq. (3). Similarly $ w(z) $ and $ q(z) $ can be obtained from Eqs. (4) and (5), respectively. In a spatially flat universe, the normalized co-movig distance $ d_c(z) $ is related to the luminosity distance $ D_{\rm L}(z) $ by

The luminosity distance can be measured from the brightness of SNe Ia. SNe Ia are often assumed to have an approximately constant absolute magnitude after the color and stretch corrections; therefore, they are widely regarded as the standard candles. The distance modulus of SNe Ia can be derived from the observation of light curves through the empirical relation [31–33]

where $ m_B^* $ is the B-band apparent magnitude, $ X_1 $ and $ {\cal C} $ are the stretch and color parameters respectively, and $ M_{\rm B} $ is the absolute magnitude. $ \alpha $ and $ \beta $ are two nuisance parameters. The luminosity distance of SNe Ia can be calculated from the distance modulus through the following relation

Several SNe Ia samples have been released, among which the most up-to-date is the Pantheon compilation [27]. The Pantheon sample is at present the largest sample that consists of different supernovae surveys, including SDSS, SNLS, various low-z samples, and some high-z samples from the HST. The total number of SNe in the Pantheon dataset is 1048, which is about twice that of the Union2.1 sample [34] and is approximately 40% more than that of the JLA sample [35]. Moreover, the furthest SNe reaches approximately $ z\sim 2.3 $, and the systematic uncertainty is further reduced compared with the previous samples. Usually, the nuisance parameters $ \alpha $ and $ \beta $ are optimized simultaneously with the cosmological parameters or are marginalized over. However, this method is model dependent; thus, the distance calibrated in one cosmological model could not be directly used to constrain the other models. The Pantheon sample applies a new method called the BEAMS with Bias Corrections (BBC) to calibrated the SNe. According to the BBC method, the nuisance parameters $ \alpha $ and $ \beta $ are determined by fitting to an randomly chosen reference cosmology. There is no special requirement on the reference cosmology; however, it should not deviate too far from the data. Once $ \alpha $ and $ \beta $ are determined, we can fix them in other cosmology fits. In the Pantheon sample, the corrected apparent magnitude $ m_{B, {\rm corr}}^* = m_B^*+\alpha X_1-\beta {\cal C} $ are reported. Therefore, we do not need to do the color and stretch corrections any more; hence, we fix $ \alpha = \beta = 0 $ in Eq. (7). The statistical uncertainty $ { D}_{\rm stat} $ and systematic uncertainty $ { C}_{\rm sys} $ are also given in Ref. [27]. The total uncertainty matrix of distance modulus is given by

We convert the distance modulus of SNe to the normalized co-moving distance through the relation

The uncertainty of $ d_c $ is propagated from the uncertainties of $ \mu $ and $ H_0 $ using the standard error propagation formula,

where $ \sigma_{H_0} $ is the uncertainty of Hubble constant, the superscript "T" denotes the transpose of a matrix, $ { D}_1 $ and $ { D}_2 $ are the Jacobian matrices,

where $ { d}_c $ is a vector whose components are the normalized co-moving distances of all the SNe Ia in Pantheon, and $ {\rm diag}({ v}) $ is the square diagonal matrix with the elements of vector v on the main diagonal.

To obtain the co-moving distance - redshift relation from the discrete data points, we use the GP method [19] to reconstruct the $ d_c(z) $ function, and the derivatives $ d_c'(z) $ and $ d_c''(z) $ can be obtained simultaneously in the reconstruction procedure. The GP can reconstruct a function $ y = f(x) $ from the discrete data points $ (x_i,y_i $) without assuming a particular parametrization of the function $ f(x) $. It assumes that the data points are drawn from the multivariate Gaussian distribution,

where $ { x} = \{x_i\} $, $ { y} = \{y_i\} $, $ { \mu} $ is the mean of the Gaussian distribution, C is the covariance matrix of the data, and $ [{ K}({ x},{ x})]_{ij} = k(x_i,x_j) $ is another covariance matrix that controls the behavior of the reconstructed function. All the freedoms of GP originate from the choice of the covariance function $ k(x_i,x_j) $. There are several covariance functions available; however, any covariance function should be symmetric, positive definite, and monotonously decreasing with $ |x_i-x_j| $. In this study, we use the simplest and most widely used squared-exponential covariance function defined by

The hyperparameters $ \sigma_f $ and l are optimized by maximizing the marginalized likelihood. For more details on the GP, please refer to Ref. [19].

It has been noted that although the GP method can reconstruct the function $ f(x) $ with a relatively high precision, the reconstructed derivatives, especially the higher order derivatives of $ f(x) $ have large uncertainty. If there are observational constraints on the derivatives of the function, then the uncertainty can be reduced. From Eq. (3) we can know that $ d'_c(z) = 1/E(z) = H_0/H(z) $; thus, the direct measurement of Hubble parameters can be used as an additional constraint on the first-order derivative of $ d_c(z) $. Here, we use the 51 $ H(z) $ data points (except for $ H_0 $) complied in Ref. [36], which is, to the best our knowledge, the largest data sample available at present. These $ H(z) $ data points are measured from two different methods, i.e., the differential age of galaxies (DAG) method and baryon acoustic oscillations (BAO) method. Because the BAO method relies on the cosmological model, to avoid the model dependence we only use the remaining 31 data points measured from the DAG method.

To normalize $ H(z) $, the precise measurement of $ H_0 $ is necessary. It is well known that the is greater than $ 3\sigma $ between the values measured from the local distance ladders and that from the global CMB radiation, where the former gives $ H_0 = 73.24\pm 1.74\; {\rm km\; s^{-1}\; Mpc^{-1}} $ [28] and the latter gives $ H_0 = 67.4\pm 0.5\; {\rm km\; s^{-1}\; Mpc^{-1}} $ [6]. To investigate the influence of Hubble constant on the reconstruction, we consider these two different values as the prior on $ H_0 $.

The publicly available python package $ {GaPP}$ [19] is used to do the GP reconstructions. First, we reconstruct from the Pantheon dataset. Then, we add the $ H(z) $ data to make a combined reconstruction. Because $ w(z) $ depends on the matter density parameter $ \Omega_M $, we fix it to the value of Planck 2018 results, i.e., $ \Omega_M = 0.315 $ [6]. The impact of different $ \Omega_M $ values on the reconstruction of $ w(z) $ will be discussed later. Note that the other quantities ($ d_c(z) $, $ E(z) $ and $ q(z) $) are independent of $ \Omega_M $. The absolute magnitude of SNe Ia is degenerated with the Hubble constant, and we fix it to $ M_{\rm B} = -19.35 $, the best-fitting value of $ \Lambda {\rm{CDM}}$ model.

The GP reconstructions of $ d_c(z) $, $ E(z) $, $ q(z) $, and $ w(z) $ from SNe data with prior $ H_0 = 67.4\pm 0.5\; {\rm km\; s^{-1}\; Mpc^{-1}} $ plotted in Fig. 1. The blue curves are the reconstructed central values, and the shaded regions are the $ 1\sigma $ and $ 2\sigma $ uncertainties. For comparison, we also plot the best-fitting curves of $ \Lambda {\rm{CDM}}$ model, with the best-fitting parameters $ \Omega_M = 0.3 $, $ H_0 = 68.9\; {\rm km\; s^{-1}\; Mpc^{-1}} $. Figure 1 shows that the reconstructions of all the four quantities are well matched to the $ \Lambda {\rm{CDM}}$ model in low redshift (z < 0.5) region. However, in the redshift region higher than 0.5, the reconstructed curves do not agree with the $ \Lambda {\rm{CDM}}$ model. This disagreement is especially obvious for $ w(z) $, which is more than $ 2\sigma $ deviation from −1 in the intermediate redshift region. Owing to the large uncertainties at z > 0.5, the remaining three quantities ($ d_c(z) $, $ E(z) $, and $ q(z) $) are still consistent with the $ \Lambda {\rm{CDM}}$ model within $ 2\sigma $ confidence level.

Figure 1.  (color online) GP reconstructions of $d_c(z)$, $E(z)$, $q(z)$, and $w(z)$ from SNe, with prior $H_0 = 67.4\pm 0.5\; {\rm km\; s^{-1}\; Mpc^{-1}}$. The black dashed curves are the predictions of flat $\Lambda$CDM model with parameters $\Omega_M = 0.3$, $H_0 = 68.9\; {\rm km\; s^{-1}\; Mpc^{-1}}$.

The GP reconstructions from the same data but with prior $ H_0 = 73.24\pm 1.74\; {\rm km\; s^{-1}\; Mpc^{-1}} $ are plotted in Fig. 2. Similar to Fig. 1, the reconstructions of $ d_c(z) $, $ E(z) $, and $ q(z) $ are consistent with the $ \Lambda {\rm{CDM}}$ model within $ 2\sigma $ confidence level, especially in the low redshift region; they are excellently in agreement with $ \Lambda {\rm{CDM}}$. However, the reconstruction of $ w(z) $ shows a discrepancy greater than $ 2\sigma $ from −1 in the low redshift region; however, it is consistent with −1 in the intermediate and high redshift regions. This is contrary to Fig. 1, which shows the discrepancy from $ \Lambda {\rm{CDM}}$ in the intermediate redshift region. This implies that the $ H_0 $ value has a significant impact on the reconstruction of $ w(z) $.

Figure 2.  (color online) Same as Fig.1 but with prior $H_0 = 73.24\pm $$ 1.74\; {\rm km\; s^{-1}\; Mpc^{-1}}$.

The tension is greater than $ 3\sigma $ tension between the two $ H_0 $ priors we used here. A wrong $ H_0 $ prior may lead to a wrong result on the reconstruction of $ w(z) $. It is interesting to see which of the two $ H_0 $ priors is more consistent with the SNe Ia data. To this end, instead of reconstructing the normalized comoving distance $ d_c(z) $, we first directly reconstruct the un-normalized co-moving distance $ D_c(z) $ (and its derivatives) as the function of redshift, where

Because $ D_{\rm L}(z) $ is directly measured from SNe Ia, $ D_c(z) $ is independent of $ H_0 $. Because $ d'_c(0) = 1/E(0)\equiv 1 $, we can infer the Hubble constant by

where $ D'_c(0) $ is the derivative of $ D_c(z) $ at z = 0. Using this method, we obtain $ H_0 = 70.6\pm 0.5\; {\rm km\; s^{-1}\; Mpc^{-1}} $, which is approximately the mean value of the local and global measurements of $ H_0 $. Then, we normalize $ D_c(z) $ (and its derivatives) with the inferred $ H_0 $ and calculate $ E(z) $, $ q(z) $, and $ w(z) $ in the same manner as the previous cases. The results are plotted in Fig. 3. Now, the reconstructed $ w(z) $ is consistent with −1 within $ 2\sigma $ confidence level in the entire redshift region, except for a small region near $ z\sim 1.2 $. In addition, from the $ q(z) $ subfigure we see a turn point at $ z = 0.59_{-0.06}^{+0.08} $, where the universe changes from accelerating to decelerating. The location of the turn point is in good agreement with the prediction of $ \Lambda {\rm{CDM}}$ model.

Figure 3.  (color online) Same as Fig.1 but with no prior on $ H_0$.

The reconstructed quantities, especially $ E(z) $, $ q(z) $, and $ w(z) $, which depend on the derivatives of $ d_c(z) $, have a large uncertainty in the high redshift region. To reduce the uncertainty, we combine the SNe Ia with the 31 DAG $ H(z) $ data in the reconstruction, where the inverse of the normalized $ H(z) $ data are treated as an additional constraint on the first-order derivative of $ d_c(z) $. The reconstruction from the SNe+$ H(z) $ data with prior $ H_0 = 67.4\pm$$ 0.5\; {\rm km\; s^{-1}\; Mpc^{-1}} $ and prior $ H_0 = 73.24\pm1.74\; {\rm km\; s^{-1}\; Mpc^{-1}} $ are plotted in Fig. 4 and Fig. 5, respectively. Compared with Fig. 1 and Fig. 2, we may see that adding the $ H(z) $ data can partially reduce the uncertainty. However, now the constructed $ d_c(z) $ is not consistent with $ \Lambda {\rm{CDM}}$ within $ 2\sigma $ confidence level in the intermediate redshift region. The $ H_0 $ prior directly affects $ E(z) $; thus, it affects the slope of $ d_c(z) $. With the small $ H_0 $ prior, the slope of $ d_c(z) $ is small; therefore, the reconstructed $ d_c(z) $ increases slower than the prediction of $ \Lambda {\rm{CDM}}$. On the contrary, with the large $ H_0 $ prior, the reconstructed $ d_c(z) $ increases faster than the prediction of $ \Lambda {\rm{CDM}}$. Similar to the SNe only case, with a small $ H_0 $ prior, $ w(z) $ deviates from −1 in the intermediate redshift region, while with large $ H_0 $ prior, $ w(z) $ deviates from −1 in the low redshift region.

Figure 4.  (color online) GP reconstruction of $d_c(z)$, $E(z)$, $q(z)$, and $w(z)$ from SNe+$H(z)$, with prior $H_0 = 67.4\pm 0.5\; {\rm km\; s^{-1}\; Mpc^{-1}}$. The black dashed curves are the predictions of flat $\Lambda$CDM model with parameters $\Omega_M = 0.3$, $H_0 = 68.9\; {\rm km\; s^{-1}\; Mpc^{-1}}$.

Figure 5.  (color online) Same to Fig.4 but with prior $ H_0 = 73.24\pm $$1.74\; {\rm km\; s^{-1}\; Mpc^{-1}}$.

Finally, we reconstruct $ d_c(z) $, $ E(z) $, $ q(z) $, and $ w(z) $ without $ H_0 $ prior using a method similar to the one mentioned above. We directly reconstruct the un-normalized $ D_c(z) $, treating the $ c/H(z) $ data as an additional constraint on the first derivative of $ D_c(z) $. Then, we infer $ H_0 $ and its uncertainty from the reconstructed $ D_c(z) $ curve using Eq. (17). The results are plotted in Fig. 6. The inferred Hubble constant is $ H_0 = 70.5\pm 0.5\; {\rm km\; s^{-1}\; Mpc^{-1}} $, which agrees with the value inferred from SNe Ia only. The $ d_c(z) $ is excellently in agreement with $ \Lambda {\rm{CDM}}$ model in the entire redshift region, and $ E(z) $ is consistent with the $ \Lambda {\rm{CDM}}$ model within $ 1\sigma $ confidence level. Compared with Fig. 3, an obvious difference can be seen in the $ q(z) $ subfigure. Except for an unambiguous accelerating-to-decelerating turn point at $ z = 0.59_{-0.05}^{+0.05} $, there is another possible, but with large uncertainty, turn point near $ z\sim 1.8 $. In addition, $ q(z) $ and $ w(z) $ are consistent with $ \Lambda {\rm{CDM}}$ within $ 2\sigma $ confidence level. Therefore, we conclude that the combined data of SNe+$ H(z) $ shows no evidence for the deviation from the standard cosmological model.

Figure 6.  (color online) Same as Fig.4 but with no prior on $ H_0$.

To investigate the influence of $ \Omega_M $ on $ w(z) $, we reconstruct $ w(z) $ from SNe+$ H(z) $ data with different $ \Omega_M $ values and without $ H_0 $ prior. Define the significance of deviation of $ w(z) $ from −1 as

where w and $ \sigma_w $ are the reconstructed central value and $ 1\sigma $ uncertainty, respectively. We plot $ \sigma(z) $ for different $ \Omega_M $ values ($ \Omega_M = [0.27,0.30,0.315,0.33] $) in Fig. 7. It is shown that the $ \Omega_M $ value has a significant effect on the reconstruction of $ w(z) $. For $ \Omega_M = 0.27 $, $ w(z) $ deviates from −1 at more than $ 3\sigma $ confidence level in the intermediate redshift region. For $ \Omega_M = 0.3 $ and $ \Omega_M = 0.33 $, $ w(z) $ deviates from −1 at more than $ 2\sigma $ confidence level near $ z\sim 1.0 $ and $ z\sim 0.3 $, respectively. However, for $ \Omega_M = 0.315 $,$ w(z) $ is consistent with −1 within $ 2\sigma $ confidence level in the entire redshift region. It is interesting that the deviation of $ w(z) $ from −1 reaches its peak value at redshift $ z \sim 1.0 $ for any $ \Omega_M $. Owing to the large uncertainty of $ w(z) $ in the high redshift region, $ w(z) $ is consistent with −1 within $ 1\sigma $ uncertainty for all $ \Omega_M $ values. Therefore, a precise measurement of the matter density parameter $ \Omega_M $ is necessary to tightly constrain the dark energy.

Figure 7.  (color online) The deviation of $ w(z)$ from −1 for different $ \Omega_M$ values.

Note that the uncertainty of $ E(z) $, $ q(z) $, and $ w(z) $ increases sharply beyond redshift $ z\sim 1.5 $ owing to the lack of data points. Adding the $ H(z) $ data can partially reduce the uncertainty; however, it is still unacceptably large. Especially, the $ w(z) $ reconstructed from SNe+$ H(z) $ has a sudden break and the uncertainty blows up near $ z\sim 2 $, which is unreasonable. This limiatation of the GP method has already been noted in Ref. [19]. To overcome this flaw, more data points in the high redshift region must be considered.

Recently, Zhang and Li [21] used the Union2.1 and JLA compilations of SNe Ia combined with the $ H(z) $ data to reconstruct the dark energy in the redshift region z < 1.5. They found that the Union2.1+$ H(z) $ and JLA+$ H(z) $ data give similar results, i.e., slight dynamical dark energy can be observed for both; however, the constant dark energy cannot be excluded. They also investigated the effect of $ H_0 $ and $ \Omega_M $ on the construction, and they showed that $ H_0 $ has a notable influence on the results; however, the influence of $ \Omega_M $ is slight. In our work, with the most recent SNe Ia data, we reconstructed the dark energy up to redshift $ z\sim 2.5 $ and obtained similar results to those in Ref. [21]. However, our results show that the effect of both $ H_0 $ and $ \Omega_M $ is non-negligible on the reconstruction of dark energy. This difference may be caused by the reduction of the uncertainty at z < 1.5. From Fig. 7, we see that $ \Omega_M $ value only affects the result bellow $ z\sim 1.5 $. Beyond $ z\sim 1.5 $, owing to the large uncertainty, the influence of $ \Omega_M $ is negligible.

In this study, we reconstructed the EoS of dark energy and the cosmic expansion from SNe Ia using the non-parametric method. To improve the significance, we also added the direct measurement of Hubble parameter, $ H(z) $ data, to construct a combined reconstruction. However, even with the $ H(z) $ data, the reconstruction still has a large uncertainty in the high redshift region. It was found that the $ H_0 $ value has a strong effect on the reconstruction. Without $ H_0 $ prior, the inferred Hubble constant from the combination of SNe+$ H(z) $ data is $H_0 = 70.5\pm $$ 0.5\; {\rm km\; s^{-1}\; Mpc^{-1}} $, thus alleviating the tension between the local and global measurements of $ H_0 $. We have also investigate if the matter density parameter, $ \Omega_M $, has some influence on the reconstruction. It was shown that the reconstruction of $ w(z) $ strongly depends on $ \Omega_M $. With the inferred Hubble constant and the Planck 2018 matter density parameter ($ \Omega_M = 0.315 $), the reconstructed $ w(z) $ is consistent with $ \Lambda {\rm{CDM}}$ model within $ 2\sigma $ confidence level. With the current observational accuracy, it is still premature to distinguish between the evolving and constant dark energy.