Dynamical dark energy in light of the latest observations

TL;DR

This paper explores whether the tensions in cosmological data can be explained by a dynamical dark energy model, which fits current observations better than the standard cosmological constant, and predicts it could be decisively detected by future surveys.

Contribution

The study demonstrates that a dynamical dark energy model can alleviate existing tensions in cosmological data and could be confirmed by upcoming observational surveys.

Findings

Dynamical dark energy reduces dataset tensions significantly.

Current data slightly favor dynamical dark energy over ΛCDM.

Future surveys like DESI could decisively detect dynamical dark energy.

Abstract

A flat Friedman-Roberson-Walker universe dominated by a cosmological constant ($\Lambda$) and cold dark matter (CDM) has been the working model preferred by cosmologists since the discovery of cosmic acceleration. However, tensions of various degrees of significance are known to be present among existing datasets within the $\Lambda$CDM framework. In particular, the Lyman-$\alpha$ forest measurement of the Baryon Acoustic Oscillations (BAO) by the Baryon Oscillation Spectroscopic Survey (BOSS) prefers a smaller value of the matter density fraction $\Omega_{\rm M}$ compared to the value preferred by cosmic microwave background (CMB). Also, the recently measured value of the Hubble constant, $H_0=73.24\pm1.74 \ {\rm km}\ {\rm s}^{-1} \ {\rm Mpc}^{-1}$, is $3.4\sigma$ higher than $66.93\pm0.62 \ {\rm km}\ {\rm s}^{-1} \ {\rm Mpc}^{-1}$ inferred from the Planck CMB data. In this work, we investigate if these tensions can be interpreted as evidence for a non-constant dynamical dark energy (DE). Using the Kullback-Leibler (KL) divergence to quantify the tension between datasets, we find that the tensions are relieved by an evolving DE, with the dynamical DE model preferred at a $3.5\sigma$ significance level based on the improvement in the fit alone. While, at present, the Bayesian evidence for the dynamical DE is insufficient to favour it over $\Lambda$CDM, we show that, if the current best fit DE happened to be the true model, it would be decisively detected by the upcoming DESI survey.

Figures

a. They indicate that the $H_{0}$, Ly$\alpha$FB and SNe measurements are in tension with the combined dataset. Introducing Tension $T$ as the number of standard deviations by which Surprise is greater than zero, we find values of $T=4.4$, $3.5$, and $1.7$ for the $H_{0}$, Ly$\alpha$FB and SNe measurements, respectively (shown in Fig. 1b), with the first two values signalling significant tension.

a. They indicate that the $H_{0}$, Ly$\alpha$FB and SNe measurements are in tension with the combined dataset. Introducing Tension $T$ as the number of standard deviations by which Surprise is greater than zero, we find values of $T=4.4$, $3.5$, and $1.7$ for the $H_{0}$, Ly$\alpha$FB and SNe measurements, respectively (shown in Fig. 1b), with the first two values signalling significant tension.

Tables

Table 1a	P15	JLA	gBAO	$P\left(k\right)$	WL	$H_0$	Ly$\alpha$FB	OHD
${\chi }_{w\left(z\right)\mathrm{CDM}}^2-{\chi }_{\Lambda \mathrm{CDM}}^2$	$-0.7$	$-1.6$	$-2.8$	$+1.1$	$-0.1$	$-2.9$	$-3.7$	$-2.3$
${\chi }_{w\mathrm{CDM}}^2-{\chi }_{\Lambda \mathrm{CDM}}^2$	$0.0$	$+0.5$	$+0.7$	$+0.4$	$+0.2$	$-2.9$	$-0.2$	$0.0$
${\chi }_{w_0{w}_a\mathrm{CDM}}^2-{\chi }_{\Lambda \mathrm{CDM}}^2$	$-0.7$	$+0.4$	$+0.9$	$+0.5$	$+0.4$	$-2.7$	$-0.3$	$0.0$
Table 1b	ALL12			ALL16			DESI++
S/N	$2.5\sigma$			$3.5\sigma$			$6.4\sigma$
$\ln \left[E_{w\left(z\right)\mathrm{CDM}}/E_{\Lambda \mathrm{CDM}}\right]$	$-6.7\pm 0.3$			$-3.3\pm 0.3$			$11.3\pm 0.3$
Table 1c	$M_1$		$M_2$	$M_3$		$M_4$		$M_5$
Definition	$M_{\nu }=0.06$ eV		$M_{\nu },N_{\mathrm{eff}}$	$A_{\mathrm{Lens}}$		$w$CDM		$w_0{w}_a$CDM
S/N	$3.6\sigma$		$3.4\sigma$	$3.4\sigma$		$3.5\sigma$		$3.4\sigma$
$\Delta \ln E$	$-2.8\pm 0.3$		$-3.9\pm 0.3$	$-3.6\pm 0.3$		$-0.9\pm 0.3$		$0.7\pm 0.3$

Acronym	Meaning	References
P15	The Planck 2015 CMB power spectra	[1]
JLA	The JLA supernovae	[2]
6dF	The 6dFRS (6dF) BAO	[3]
MGS	The SDSS main galaxy sample BAO	[4]
$P\left(k\right)$	The WiggleZ galaxy power spectra	[5]
WL	The CFHTLenS weak lensing	[6]
$H_0$	The Hubble constant measurement	[7]
OHD	$H\left(z\right)$ from galaxy age measurements	[8]
gBAO-$3z$	3-bin BAO from BOSS DR12 galaxies	[9]
gBAO-$9z$	9-bin BAO from BOSS DR12 galaxies	[10]
Ly$\alpha$FB	The Ly$\alpha$ forest BAO measurements	[11]
B	P15 + JLA + 6dF + MGS
ALL12	The combined dataset used in Zhao et al., (2012)[12]
ALL16-$3z$	B+$P\left(k\right)$+WL+$H_0$+OHD+gBAO-$3z$+Ly$\alpha$FB
ALL16	B+$P\left(k\right)$+WL+$H_0$+OHD+gBAO-$9z$+Ly$\alpha$FB
DESI++	P15 + mock DESI BAO[13] + mock SNe[14]

	ALL12	ALL16
SNe	SNLS3 (472) [15]	JLA (740) [2]
gBAO	BOSS DR9; $z_{\mathrm{eff}}=0.57$ [16]	BOSS DR12; nine $z_{\mathrm{eff}}\in [0.2,0.75]$ [10]
Ly$\alpha$FB	none	BOSS DR11; $z_{\mathrm{eff}}=2.34$ [11]
$H_0$ [km s-1 Mpc-1]	$74.2\pm 3.6$ [17]	$73.24\pm 1.74$ [7]
CMB	WMAP7 [18]	Planck 2015 [1]

		$\Lambda \mathrm{CDM}$			$w\left(z\right)\mathrm{CDM}$
${10}^3{\Omega }_b{h}^2$	$22.4\pm 0.13$	$22.5\pm 0.16$	$22.4\pm 0.14$	$22.3\pm 0.15$	$22.4\pm 0.21$	$22.5\pm 0.17$
${10}^3{\Omega }_c{h}^2$	$118.4\pm 0.9$	$121.9\pm 2.4$	$118.4\pm 1.0$	$119.2\pm 1.4$	$120.4\pm 2.6$	$117.7\pm 1.5$
${10}^3{\Theta }_s$	$1041\pm 0.3$	$1041\pm 0.4$	$1041\pm 0.3$	$1041\pm 0.3$	$1041\pm 0.4$	$1041\pm 0.3$
${10}^2\tau$	$7.9\pm 1.6$	$8.2\pm 1.7$	$4.4\pm 1.8$	$6.9\pm 1.7$	$7.3\pm 1.8$	$4.5\pm 1.9$
$\ln \left({10}^{10}{A}_s\right)$	$3.1\pm 0.03$	$3.1\pm 0.04$	$3.0\pm 0.04$	$3.1\pm 0.03$	$3.1\pm 0.04$	$3.0\pm 0.04$
${10}^2{n}_s$	$96.8\pm 0.4$	$97.5\pm 0.6$	$96.8\pm 0.4$	$96.6\pm 0.5$	$96.9\pm 0.8$	$97.0\pm 0.5$
$M_{\nu }\left(\mathrm{eV}\right)$	$0.06$		$0$	$0.06$		$0$
$N_{\mathrm{eff}}$	$3.04$	$3.27\pm 0.14$	$3.04$	$3.04$	$3.14\pm 0.16$	$3.04$
$A_{\mathrm{Lens}}$	$1$	$1$	$1.17\pm 0.07$	$1$	$1$	$1.17\pm 0.07$
${10}^3{\Omega }_m$	$307\pm 5.3$	$304\pm 5.8$	$301\pm 5.6$	$289\pm 11.2$	$290\pm 11.3$	$285\pm 10.9$
$H_0$	$67.9\pm 0.40$	$69.1\pm 0.85$	$68.5\pm 0.44$	$70.2\pm 1.3$	$70.5\pm 1.4$	$70.2\pm 1.3$

Equations

$$S=\frac{1}{2\ln 2}\left[({\theta }_1-{\theta }_2{)}^T{\mathcal{C}}_1^{-1}({\theta }_1-{\theta }_2)-\mathbf{Tr}\left({\mathcal{C}}_2{\mathcal{C}}_1^{-1}+\mathbb{I}\right)\right]$$

$$\Sigma =\frac{1}{\sqrt{2}\ln 2}\sqrt{\mathbf{Tr}{\left({\mathcal{C}}_2{\mathcal{C}}_1^{-1}+\mathbb{I}\right)}^2}.$$

Full text

The observational datasets considered in this work include the latest CMB temperature and polarisation anisotropy spectra, the supernovae (SNe) luminosity distance data, the BAO angular diameter distance data from the clustering of galaxies (gBAO) and from the Lyman-$\alpha$ forest (Ly$\alpha$FB), the measurement of $H_{0}$, $H(z)$ measurements using the relative age of old and passively evolving galaxies (OHD), the three-dimensional galaxy power spectra, and the two-dimensional weak lensing shear angular power spectra. Further details about the datasets and associated systematic effects can be found in Methods.

The KL divergence, also known as relative entropy, quantifies the proximity of two probability density functions (PDFs). Rather than focusing on particular model parameters, it is designed to compare the overall concordance of datasets within a given model. We use the difference between the actual and the expected KL divergence, called “Surprise”[15], as a measure of tension between datasets. Rather than comparing the PDFs for the $\Lambda$CDM parameters for every pair of datasets, we take the combined dataset, ALL16 (see Supplementary Table 1), and find the derived PDFs for the angular diameter distance $D_{A}(z)$ and the Hubble parameter $H(z)$ at redshifts corresponding to the available data. We then compute the KL divergence between the derived PDFs and the directly observed $D_{A}(z)$ and $H(z)$ from $H_{0}$, SNe, OHD, gBAO and Ly$\alpha$FB, and evaluate the corresponding Surprise and the standard deviation (see Methods for details). Results are shown with cyan bars in Fig. 1a. They indicate that the $H_{0}$, Ly$\alpha$FB and SNe measurements are in tension with the combined dataset. Introducing Tension $T$ as the number of standard deviations by which Surprise is greater than zero, we find values of $T=4.4$, $3.5$, and $1.7$ for the $H_{0}$, Ly$\alpha$FB and SNe measurements, respectively (shown in Fig. 1b), with the first two values signalling significant tension.

Next, we check if the tension within the $\Lambda$CDM model can be interpreted as evidence for a dynamical DE. The dynamics of DE can be probed in terms of its equation of state $w$, which is equal to $-1$ for $\Lambda$, but is different in dynamical DE models where it will generally be a function of redshift $z$. Commonly considered alternatives to $\Lambda$ are a model with a constant $w$ ($w$CDM), and one in which $w$ is linear function of the scale factor ($w_{0}w_{a}$CDM)[16]. We allow for a general evolution of the DE equation of state and use the correlated prior method[17] to perform a Bayesian non-parametric reconstruction of $w(z)$ using the Monte Carlo Markov Chain method with other cosmological parameters marginalised over (see Methods for details). Fig. 2 presents the reconstructed $w(z)$, along with the 68% confidence level (CL) uncertainty, shown with a light blue band, derived from the combined dataset ALL16. Table 1a shows the change in $\chi^{2}$ relative to $\Lambda$CDM for each individual dataset for the best fit $w(z)$CDM model derived from ALL16. Overall, the $\chi^{2}$ is improved by $-12.3$, which can be interpreted as the reconstructed dynamical DE model being preferred at $3.5\sigma$. The reconstructed DE equation of state evolves with time and crosses the $-1$ boundary, which is prohibited in single field minimally coupled quintessence models[18], but can be realised in models with multiple scalar fields, such as Quintom[19], or if the DE field mediates a new force between matter particles[20]. In the latter case, which is commonly classified as Modified Gravity, it is quite generic for the effective DE equation of state to be close to $-1$ around $z=0$, but evolve towards more negative values at intermediate redshifts, before eventually approaching [math] during matter domination. Such dynamics would be consistent with our reconstruction and could be tested in the future when BAO measurements at higher redshifts become available. In addition to the reconstruction from the combined ALL16 dataset presented in Fig. 2, we present reconstructions derived from ten different data combinations in Supplementary Fig. 1.

The results for tension between datasets, re-evaluated for the ALL 16 best fit $w(z)$CDM model, are shown with dark blue bars in Fig. 1. We find $T=0.7$, $1.1$ and $0.7$ for $H_{0}$, Ly$\alpha$FB and SNe, respectively, indicating that tensions that existed in the $\Lambda$CDM model are significantly released within $w(z)$CDM. A plot of the relevant data points along with the best fit predictions from the $\Lambda$CDM and the $w(z)$CDM model are provided in the Supplementary Fig. 2.

With a large number of additional $w$-bin parameters, one may be concerned that the improvement in the fit is achieved by $w(z)$CDM at the cost of a huge increase of the parameter space. However, correlations between the $w$-bins induced by the prior constrain most of that freedom. One way to estimate the effective number of additional degrees of freedom is to perform a principal component analysis (PCA) of the posterior covariance matrix of the $w$-bin parameters and compare it to that of the prior. Using this method, explained in detail in Methods, we find that our $w(z)$CDM model effectively has only four additional degrees of freedom compared to $\Lambda$CDM. We note that the demonstration that ALL16 is capable of constraining four principal components of $w(z)$ is one of the interesting results of this work.

It is interesting to compare $w(z)$ reconstructed from ALL16 to that obtained in Zhao et al. (2012)[21] using the same prior but a different dataset which we call ALL12 (a comparison of the ALL16 and the ALL12 datasets is provided in Supplementary Table 2). ALL16 contains about 40% new supernovae compared to ALL12, primarily provided by the SDSS-II survey. Moreover, in ALL12, the BAO measurement derived from the BOSS DR9 sample[22] was at a single effective redshift, while in ALL16 it is tomographic at nine redshifts from BOSS DR12[23], which contains four times more galaxies than DR9. In addition, ALL16 includes a high-redshift BAO measurement from Lyman-$\alpha$ forest, which was not available in 2012. This helps to constrain $w(z)$ at redshifts where the supernovae constraints are weak. The new 2016 $H_{0}$ measurement[12] is consistent with that in 2009[24], with the error bar halved. Comparing measurements of the expansion rate and the cosmic distances in ALL12 and ALL16, we find that those in ALL16 offer information at more redshift values, and with a greater signal-to-noise (S/N) ratio (see Supplementary Figure 4 for a visual comparison). Quantitatively, the S/N in measurements of $H(z)$, $D_{A}(z)$ and $d_{L}(z)$ in ALL16 is larger by 80%, 260% and 90%, respectively, compared to the ALL12 dataset. The Planck 2015 CMB data is also much more informative than the WMAP 7-year release[25], thanks to a higher angular resolution of the temperature and polarisation maps, and lower levels of statistical uncertainties.

Overall, ALL16 is more constraining due to a significant level of new and independent information in ALL16 compared to ALL12: the effective number of $w(z)$ degrees of freedom constrained by ALL12 was three, compared to four constrained by ALL16. Fig. 2 compares the two results and shows that they are highly consistent. We quantify the agreement by evaluating the dot-product of the $\hat{\bf w}$ vectors from the two reconstructions (the vectors are normalised so that a dot-product is unity if they are identical) and find that $\hat{\bf w}_{\rm ALL12}\cdot\hat{\bf w}_{\rm ALL16}=0.94\pm 0.02.$ We also evaluate the tension $T$ between the two reconstruction results and find that $T=-1.1$. This indicates an excellent alignment of the two results. This agreement, and the raised significance of an evolving $w(z)$ from $2.5\sigma$ to $3.5\sigma$ CL with more advanced observations, suggests the possibility of revealing the dynamics of DE at a much more statistically significant level in the near future, as we will present later.

To check whether the improvement in the fit warrants introducing additional effective degrees of freedom, we evaluate and compare the Bayesian evidence, $E\equiv\int d\theta\,{\cal L}({\bf D}|\theta)\,P(\theta)$, for $\Lambda$CDM and the $w(z)$CDM model. The Bayes’ factors, which are the differences in $\ln E$ between the two models, are shown in Table 1. The Bayes’ factors for both the ALL12 and the ALL16 DE models are negative, indicating that $\Lambda$CDM is favoured by this criterion. However, our forecast for a future dataset, DESI++, comprised of BAO measurements from DESI[14], around 4000 supernovae luminosity distances from future surveys[26] and CMB (assuming the Planck sensitivity), predicts a Bayes’ factor of $11.3\pm 0.3$ if the ALL16 $w(z)$CDM reconstruction happened to be the true model, which would be decisive according to the Jeffreys scale.

One may ask how much the evidence for DE depends on the particular choice of the prior parameters. In principle, the choice of the smoothing scale should be guided by theory. The value used in Zhao et al. (2012)[21] and this work, $a_{c}=0.06$, is a time-scale conservatively chosen to be sufficiently small not to bias reconstructions of $w(z)$ expected in quintessence DE models[17]. For the inference to be conclusive, the evidence for a dynamical DE should be strong over a wide range of the prior parameters. Therefore, we vary the strength of our prior by adjusting $\sigma_{\rm D}$, a parameter added to the diagonal of the inverse of the prior covariance matrix, and examine how the significance of the dynamical DE detection, as well as the Bayesian evidence, change with the variation of $\sigma_{\rm D}$. As shown in Fig. 3, and with additional details given in Methods, we find that neither ALL12 nor ALL16 provide evidence for a dynamical DE over the considered wide range of prior strengths. However, the Bayes factor for ALL16 is generally much less negative than that of ALL12 for all prior strengths, e.g., it increased from $-6.7\pm 0.3$ to $-3.3\pm 0.3$ for $\sigma_{\rm D}=3$, which is the prior used in this work. In fact, for ALL16, the Bayes factor remains close to zero for $\sigma_{\rm D}\lesssim 0.4$. We plot the model with $\sigma_{\rm D}=0.4$ as a light green band in Fig. 2 to demonstrate the impact of changing the prior strength, and also because it is a model that has the same Bayesian evidence as $\Lambda$CDM while deviating from $\Lambda$ at a $2.7\sigma$ CL. On the other hand, our forecast for DESI++ shows that, if the $w(z)$CDM model was true, it would be decisively supported by Bayesian evidence over a wide range of prior strengths, as shown in Fig. 3.

Various ways to relieve the tension between datasets have been proposed, including allowing for additional relativistic degrees of freedom[27], massive neutrinos[27], and interacting vacuum[29, 28]. In addition, to relieve the tension between the $\Lambda$CDM parameters required to fit the CMB temperature anisotropy spectrum at large and small scales, the Planck team introduced[7] a parameter $A_{\rm Lens}$ that rescales the amplitude of the weak lensing contribution to the temperature power spectrum. In the $w(z)$ reconstruction discussed earlier, we fixed $A_{\rm Lens}=1$, assumed that neutrinos are massless and set the effective number of relativistic species at the standard $\Lambda$CDM value of $N_{\rm eff}=3.04$. We have checked the effect of these parameters on the reconstructed $w(z)$ by considering model $M_{1}$, with $M_{\nu}$ fixed to $0.06$ eV, model $M_{2}$ , with $M_{\nu}$ and $N_{\rm eff}$ added as free parameters, and model $M_{3}$, with varied $A_{\rm Lens}$. In all these cases, we find that the shape of the reconstructed $w(z)$ and the significance of its deviation from $-1$ are practically the same. The inferred values of the cosmological parameters in these models are given in Supplementary Table 3, and the corresponding reconstructed $w(z)$ are shown in Supplementary Fig. 5. The Bayes factors for $M_{1}$, $M_{2}$ and $M_{3}$ relative to the corresponding $\Lambda$CDM models (with the same added parameters), are shown in Table 1c. We also checked that neither the constant $w$ model ($M_{4}$), nor the linear ($w_{0},w_{a}$) parametrisation of $w$ ($M_{5}$), are capable of releasing the tensions between all datasets simultaneously (see Table 1a). Interestingly, we find that our DE model with a non-parametrically reconstructed $w(z)$ has a larger Bayes factor compared to $w_{0}w_{a}$CDM despite the latter having only two parameters (see results for $M_{4}$ and $M_{5}$ in Table 1c).

There is always a possibility that the tensions between datasets, quantified in terms of the KL divergence in this work, are due to yet unknown systematic effects. However, it is intriguing that they persist with improvements in the quantity and the quality of the data[10, 11, 30]. If interpreted as a manifestation of DE dynamics, they suggest a $w(z)$ that crosses $-1$ and has a shape that is representative of modified gravity models. The commonly used ($w_{0},w_{a}$) parametrisation would have missed this behaviour and has a lower Bayesian evidence than the reconstructed $w(z)$ model, despite the latter having more degrees of freedom. Thus, our results demonstrate that the current data can provide non-trivial constraints on the DE dynamics. It is also intriguing that the evidence for $w(z)\neq-1$, while below that of $\Lambda$CDM, has become stronger with the new independent data added since 2012, and that the ALL16 reconstruction is consistent with the ALL12 $w(z)$. We emphasise that we have not optimised the prior to maximise either the Bayes ratio or the statistical significance of the departure from $-1$, as it would be contrary to the principles of Bayesian inference. Future data has the ability to conclusively confirm the DE evolution inferred in this work if it happened to be the one chosen by Nature.

0.1 Tension calculation

The Kullback-Leibler (KL) divergence[13], also known as relative entropy, has been extensively utilised as a way of quantifying the degree of tension between different datasets within the $\Lambda$CDM model[35, 34, 36, 15, 37, 31, 32, 33]. Rather than focusing on particular model parameters, it is designed to compare the overall concordance of datasets within a given model. Alternative methods of quantification of the tension have been discussed in the literature[38, 8].

The KL divergence quantifies the proximity of two probability density functions (PDFs), $P_{1}$ and $P_{2}$, of a multi-dimensional random variable $\theta$. If both $P_{1}$ and $P_{2}$ are assumed to be Gaussian[15], and data are assumed to be more informative than the priors, we can write the difference between the actual and the expected KL divergence, called the “Surprise”[15], as

[TABLE]

where $\theta_{1}$ and $\theta_{2}$ are the best-fit parameter vectors, ${\bf{\cal{C}}}_{1}$ and ${\bf{\cal{C}}}_{2}$ are the covariance matrices for $P_{1}$ and $P_{2}$, and $\mathbb{I}$ is the unity matrix. The standard deviation of the expected KL divergence is

[TABLE]

We can quantify the tension between $P_{1}$ and $P_{2}$ in terms of the signal-to-noise ratio $T=S/\Sigma$. If $T\lesssim 1$, then $P_{1}$ and $P_{2}$ are consistent with each other[32].

0.2 Datasets used

The datasets we consider include the Planck 2015 (P15) CMB temperature and polarization auto- and cross-angular power spectra[7], the JLA supernovae[39] (JLA); the 6dFRS (6dF)[40] and SDSS main galaxy sample (MGS)[41] BAO measurements, the WiggleZ galaxy power spectra[42] in four redshift slices, containing information about the Baryon Acoustic Oscillations (BAO) and Redshift Space Distortions (RSD) ($P(k)$), the weak lensing shear angular power spectra in six redshift slices from CFHTLenS[43] (WL), the recent estimate of the Hubble constant $H_{0}$ obtained from local measurements of Cepheids[12] ($H_{0}$), the $H(z)$ measurement using the relative age of old and passively evolving galaxies following a cosmic chronometer approach[44] (OHD), the BOSS DR12 “Consensus” BAO measurement (BAO-$3z$)[45], the BAO and RSD measurement using the complete BOSS DR12 sample covering the redshift range of $0.2<z<0.75$ at three effective redshifts, the BAO measurement using the same galaxy sample but at nine effective redshifts[23] (BAO-$9z$), and the Ly$\alpha$ BAO (Ly$\alpha$FB) measurements[3]. A summary of datasets and data combinations used in this work is shown in Supplementary Table 1.

We account for the systematic effects in our analysis as implemented in the public likelihood codes. However, we note that there may be additional systematic effects. For example, the relative velocity between baryons and dark matter may affect the BAO distance measurements[46]. This effect is estimated to be at sub-percent level for the galaxy BAO measurements of BOSS[47, 48], and is currently unknown for Ly$\alpha$FB. For SNe, we use the conventional $\chi^{2}$ statistics for the analysis, although alternative statistics may extract more information and reduce the systematic effects for the JLA sample to some extent[49, 50].

0.3 Non-parametric $w(z)$ reconstruction

To start, $w(z)$ is parameterised in terms of its values at discrete steps in $z$, or the scale factor $a$. Fitting a large number of uncorrelated bins would result in extremely large uncertainties and, in fact, would prevent the Monte Carlo Markov Chains (MCMC) from converging because of the many degenerate directions in parameter space. On the other hand, fitting only a few bins could significantly bias the result. Our approach is to introduce a prior covariance between the bins based on a specified two-point function that correlates values of $w$ at different $a$, $\xi_{w}(|a-a^{\prime}|)\equiv\left\langle[w(a)-w^{\rm fid}(a)][w(a^{\prime})-w^{\rm fid}(a^{\prime})]\right\rangle$, which can be taken to be of the form proposed in[51], $\xi_{\rm CPZ}(\delta a)=\xi_{w}(0)/[1+(\delta a/a_{c})^{2}]$, where $a_{c}$ describes the typical smoothing scale, and $\xi_{w}(0)$ is a normalisation factor set by the expected variance in the mean $w$, $\sigma^{2}_{\bar{w}}$. As shown in[17], results are largely independent of the choice of the correlation function. The prior covariance matrix ${\bf C}$ is obtained by projecting $\xi_{w}(|a-a^{\prime}|)$ onto the discrete $w$ bins[51, 17], and the prior PDF is taken to be of Gaussian form: ${\cal P}_{\rm prior}({\bf w})\propto\exp[-({\bf w}-{\bf w}^{\rm fid})^{T}{\bf C}^{-1}({\bf w}-{\bf w}^{\rm fid})/2]$, where ${\bf w}^{\rm fid}$ is the fiducial model. The reconstructed model is that which maximises the posterior probability, which by Bayes’ theorem is proportional to the likelihood of the data times the prior probability, ${\cal P}({\bf w}|{\bf D})\propto{\cal P}({\bf D}|{\bf w})\times{\cal P}_{\rm prior}({\bf w})$. Effectively, the prior results in a new contribution to the total $\chi^{2}$ of a model, which penalises models that are less smooth.

In our reconstruction of $w(z)$, we set $a_{c}=0.06$ and $\sigma_{\bar{w}}=0.04$, which is the “weak prior” used in Zhao et al. (2012)[21]. To calculate the observables, we use a version of ${\tt CAMB}$[52] modified to include DE perturbations for an arbitrary ${\bf w}$[53]. We use ${\tt PolyChord}$[54], a nested sampling plug-in for ${\tt CosmoMC}$[55], to sample the parameter space ${\bf P}\equiv(\omega_{b},\omega_{c},\Theta_{s},\tau,n_{s},A_{s},w_{1},...,w_{30},\mathcal{N})$ where $\omega_{b}$ and $\omega_{c}$ are the baryon and CDM densities, $\Theta_{s}$ is the angular size of the sound horizon at decoupling, $\tau$ is the optical depth, $n_{s}$ and $A_{s}$ are the spectral index and the amplitude of the primordial power spectrum, and $w_{1},...,w_{30}$ denote the 30 $w$-bin parameters. The first 29 $w$ bins are uniform in $a\in[0.286,1]$, corresponding to $z\in[0,2.5]$, and the last wide bin covers $z\in[2.5,1100]$. We marginalise over nuisance parameters such as the intrinsic SN luminosity.

0.4 Principal component analysis of $w(z)$

First, we diagonalise the posterior covariance of $w$-bins to find their uncorrelated linear combinations (eigenmodes), along with the eigenvalues, which quantify how well a given eigenmode is constrained[56]. We plot the inverse eigenvalues of the posterior covariance, ordered according to the number of nodes in the eignemodes, in Panel a of Supplementary Fig. 3. The number of nodes is representative of the smoothness in the evolution of eigenmodes, with the first four posterior eigenmodes shown in Panel b of Supplementary Fig. 3. Next, we perform a PCA of the prior covariance matrix and plot its inverse eigenvalues alongside those of the posterior. We see that the fifth and higher number eigenvalues of the two matrices coincide, which means that they are fully determined by the prior. However, the first four inverse eigenvalues of the posterior are significantly larger than that of the prior, indicating that they are constrained primarily by the data. This is precisely the intent of the correlated prior method: the smooth features in $w(z)$ are constrained by the data, with no bias induced by the prior, while the high frequency features are constrained by the prior. Thus, our $w(z)$CDM model effectively has only four additional degrees of freedom compared to $\Lambda$CDM.

0.5 Dependence of the result on the correlated prior

To investigate the dependence of our result on the strength of the correlated prior, in Fig. 3, we plot the Bayes factor and the statistical significance of $w\neq-1$ as a function of $\sigma_{\rm D}$, which is a parameter that is added to the inverse covariance matrix of the correlated prior to effectively strengthen it. We find that neither ALL12 nor ALL16 dataset provides evidence for a dynamical DE at all prior strengths, although the Bayes factor for ALL16 is generally much less negative than that of ALL12 for all priors, e.g., it grows from $-6.7\pm 0.3$ to $-3.3\pm 0.3$ for the prior used in this work, which corresponds to $\sigma_{\rm D}=3$. On the other hand, the plot shows that, if the $w(z)$CDM model was true, DESI++ would be able to provide a decisive Bayesian evidence over a wide range of prior strengths.

{addendum}

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.