Reconstruction of the Dark Energy equation of state from latest data: the impact of theoretical priors

TL;DR

This paper reconstructs the Dark Energy equation of state using non-parametric methods and theoretical priors, analyzing data from CMB, Supernovae, and BAO to explore its time dependence and model compatibility.

Contribution

It introduces a non-parametric reconstruction of the EoS with theoretical priors for viability and stability, comparing two methods and applying to broad classes of models.

Findings

Reconstruction methods yield consistent results.

EoS shows time dependence with model-dependent trends.

ΛCDM's constant EoS is within 1σ bounds.

Abstract

We reconstruct the Equation of State of Dark Energy (EoS) from current data using a non-parametric approach where, rather than assuming a specific time evolution of this function, we bin it in time. We treat the transition between the bins with two different methods, i.e. a smoothed step function and a Gaussian Process reconstruction, investigating whether or not the two approaches lead to compatible results. Additionally, we include in the reconstruction procedure a correlation between the values of the EoS at different times in the form of a theoretical prior that takes into account a set of viability and stability requirements that one can impose on models alternative to $\Lambda$CDM. In such case, we necessarily specialize to broad, but specific classes of alternative models, i.e. Quintessence and Horndeski gravity. We use data coming from CMB, Supernovae and BAO surveys. We find an overall agreement between the different reconstruction methods used; with both approaches, we find a time dependence of the mean of the reconstruction, with different trends depending on the class of model studied. The constant EoS predicted by the $\Lambda$CDM model falls anyway within the $1\sigma$ bounds of our analysis.

Tables

Table 1. Table 1 : Summary of the autocorrelation and correlation analytical fits obtained by [ 41 ] . The first parameterization of the correlation matrix is the exponential, while the second is the CPZ [ 43 ] . Here x 0 subscript 𝑥 0 x_{0} denotes the x 𝑥 x variable evaluated at present time, i.e. x 0 = a 0 = 1 subscript 𝑥 0 subscript 𝑎 0 1 x_{0}=a_{0}=1 (or x 0 = ln ⁡ ( a 0 ) = 0 subscript 𝑥 0 subscript 𝑎 0 0 x_{0}=\ln(a_{0})=0 ).

Autocorrelation function		$\alpha$	$\beta$	$\gamma$	$x$
$C\left(x\right)=\alpha +\beta \exp \left[\gamma \left(x-x_0\right)\right]$	Quintessence	$0.03$	$0.3$	$6.5$	$a$
	Horndeski	$0.05$	$0.8$	$2$	$\ln \left(a\right)$
Correlation matrix		$\xi$	$n$	$x$	$y$
$\tilde{C}(x,y)=\exp \left[{\left(-\frac{\left|x-y\right|}{\xi }\right)}^n\right]$	Quintessence	$0.7$	$1.8$	$a$	$a^{\prime }$
	Horndeski	$0.3$	$1.2$	$\ln \left(a\right)$	$\ln \left(a^{\prime }\right)$
$\tilde{C}(x,y)=\frac{1}{1+{\left(\frac{\left|x-y\right|}{\xi }\right)}^n}$	Quintessence	$0.6$	$2$	$a$	$a^{\prime }$
	Horndeski	$0.2$	$2$	$\ln \left(a\right)$	$\ln \left(a^{\prime }\right)$

Table 2. Table 2 : Mean values and 1 ​ σ 1 𝜎 1\sigma confidence levels or upper limits of the Quintessence case inferred parameters, using the datasets without and with the priors, reconstructing the Equation of State via smoothed step function. The input scale factors associated to the w → → 𝑤 \vec{w} are a → = ( 0.85 , 0.7 , 0.55 , 0.4 , 0.25 , 0.1 ) → 𝑎 0.85 0.7 0.55 0.4 0.25 0.1 \vec{a}=(0.85,~{}0.7,~{}0.55,~{}0.4,~{}0.25,~{}0.1) .

SMOOTHED STEP FUNCTION RECONSTRUCTION
Parameter	Planck+BAO+SN	Planck+BAO+SN (exp prior)	Planck+BAO+SN (CPZ prior)
${\Omega }_b{h}^2$………	$0.02244\pm 0.00021$	$0.02244\pm 0.00020$	$0.02242\pm 0.00021$
${\Omega }_c{h}^2$………	$0.1164\pm 0.0016$	$0.1166\pm 0.0015$	$0.1167\pm 0.0015$
$H_0$ ………..	${64.9}_{-1.1}^{+1.6}$	${65.2}_{-0.96}^{+1.6}$	${65.3}_{-0.95}^{+1.4}$
${\Omega }_{\Lambda }$…………	${0.669}_{-0.011}^{+0.017}$	${0.671}_{-0.010}^{+0.017}$	${0.672}_{-0.010}^{+0.015}$
${\Omega }_m$…………	${0.331}_{-0.017}^{+0.011}$	${0.329}_{-0.017}^{+0.010}$	${0.328}_{-0.015}^{+0.010}$
$w_1$……..
$w_2$……..
$w_3$……..
$w_4$……..
$w_5$……..
$w_6$……..

Table 3. Table 3 : Mean values and 1 ​ σ 1 𝜎 1\sigma confidence levels of the Quintessence case inferred parameters, using the datasets without and with the CPZ prior, reconstructing the Equation of State via Gaussian Process method. The input scale factors associated to the w → → 𝑤 \vec{w} are a → = ( 0.85 , 0.7 , 0.55 , 0.4 , 0.25 , 0.1 ) → 𝑎 0.85 0.7 0.55 0.4 0.25 0.1 \vec{a}=(0.85,~{}0.7,~{}0.55,~{}0.4,~{}0.25,~{}0.1) .

GAUSSIAN PROCESS RECONSTRUCTION
Parameter	Planck+BAO+SN	Planck+BAO+SN (CPZ prior)
${\Omega }_b{h}^2$………	${0.02242}_{-0.00020}^{+0.00022}$	$0.02241\pm 0.00021$
${\Omega }_c{h}^2$………	$0.1166\pm 0.0015$	$0.1167\pm 0.0014$
$H_0$ ………..	${62.8}_{-2.0}^{+1.6}$	${63.4}_{-1.6}^{+2.2}$
${\Omega }_{\Lambda }$…………	$0.645\pm 0.020$	${0.651}_{-0.017}^{+0.025}$
${\Omega }_m$…………	$0.355\pm 0.020$	${0.349}_{-0.025}^{+0.017}$
$w_1$……..
$w_2$……..
$w_3$……..
$w_4$……..
$w_5$……..
$w_6$……..

Table 4. Table 4 : Mean values and 1 ​ σ 1 𝜎 1\sigma confidence levels of the Horndeski case, reconstructing the Equation of State via smoothed step function and via Gaussian Process. The input scale factors associated to the w → → 𝑤 \vec{w} are a → = ( 0.9 , 0.8 , 0.7 , 0.6 , 0.5 , 0.4 ) → 𝑎 0.9 0.8 0.7 0.6 0.5 0.4 \vec{a}=(0.9,~{}0.8,~{}0.7,~{}0.6,~{}0.5,~{}0.4) . Here w 6 subscript 𝑤 6 w_{6} is unconstrained, with the 1 ​ σ 1 𝜎 1\sigma region extending over the full prior range.

RECONSTRUCTIONS: JLA+BAO (exp prior)
Parameter	Smoothed step function	Gaussian Process
${\Omega }_b{h}^2$………	${0.041}_{-0.031}^{+0.012}$	${0.049}_{-0.031}^{+0.017}$
${\Omega }_c{h}^2$………	${0.181}_{-0.12}^{+0.084}$	${0.157}_{-0.10}^{+0.067}$
$H_0$ ………..	$>76.6$	$>79.7$
${\Omega }_{\Lambda }$…………	$0.684\pm 0.086$	$0.713\pm 0.072$
${\Omega }_m$…………	$0.316\pm 0.086$	$0.287\pm 0.072$
$w_1$……..	$-{1.10}_{-0.22}^{+0.30}$	$-{1.10}_{-0.23}^{+0.28}$
$w_2$……..	$-{0.97}_{-0.22}^{+0.36}$	$-{0.86}_{-0.18}^{+0.26}$
$w_3$……..	$-{1.09}_{-0.30}^{+0.60}$	$-{0.98}_{-0.25}^{+0.44}$
$w_4$……..	$-{1.21}_{-0.39}^{+0.73}$	$-{1.05}_{-0.34}^{+0.58}$
$w_5$……..	$-{1.44}_{-0.63}^{+0.75}$	$-{1.36}_{-0.62}^{+0.74}$
$w_6$……..	unconstrained	unconstrained

Equations

$$d_L(a)=\frac{c}{a{H}_0}{\int }_a^1\frac{d{a}^{\prime }}{E(a^{\prime })}$$

$$E(a)\approx \sqrt{{\Omega }_m{a}^{-3}+{\Omega }_{\mathrm{DE}}\exp \left[{\int }_a^1\frac{3[1+w_{\mathrm{DE}}(a)]}{a}d{a}^{\prime }\right]}.$$

$$w_{\mathrm{DE}}(a)=w_1+\sum _{i=1}^{N-1}\left(w_{i+1}-w_1\right)\left[{\theta }_H(a-a_i)-{\theta }_H(a-a_{i+1})\right],$$

$$w_{\rm DE}(a)=w_{1}+\sum_{i=1}^{N-1}\dfrac{w_{i+1}-w_{i}}{2}\bigg{\{}1+\tanh\left[s\left(\dfrac{a-a_{i+1}}{a_{i+1}-a_{i}}\right)\right]\bigg{\}}$$

$$\vec{w}_{\rm DE}\\ =\begin{bmatrix}w_{1}\\ w_{2}\\ w_{3}\end{bmatrix}\sim\mathcal{N}(\vec{\mu},C)\\ =\mathcal{N}\bigg{(}\vec{\mu},\begin{bmatrix}C_{11}&C_{12}&C_{13}\\ C_{21}&C_{22}&C_{23}\\ C_{31}&C_{33}&C_{33}\end{bmatrix}\bigg{)}$$

$$C(a,a^{\prime })=e^{-\frac{|a-a^{\prime }{|}^2}{2{\xi }^2}}$$

$$C(a,a^{\prime})=\bigg{\{}\begin{array}[]{ll}0&~{}~{}|a-a^{\prime}|\gg\xi\\ 1&~{}~{}a=a^{\prime}\\ \end{array}$$

$$\begin{bmatrix}\vec{w}_{\rm DE}\\ w_{*}\\ \end{bmatrix}\sim\mathcal{N}\bigg{(}\begin{bmatrix}\vec{\mu}\\ \mu_{*}\\ \end{bmatrix},\begin{bmatrix}C(\vec{a},\vec{{a}}^{\prime})&C(\vec{a},a_{*})\\ C(a_{*},\vec{a})&C(a_{*},a_{*})\\ \end{bmatrix}\bigg{)}$$

$$w_{\mathrm{DE}}(\vec{a})\sim GP(\mu (\vec{a}),C(\vec{a},{\vec{a}}^{\prime }))$$

$$-\ln \mathcal{L}={\chi }^2={\chi }_{\mathrm{data}}^2+{\chi }_{\mathrm{prior}}^2,$$

$${\chi }_{\mathrm{prior}}^2=(\mathbf{w}-\bar{\mathbf{w}}{)}^T{\mathcal{C}}^{-1}(\mathbf{w}-\bar{\mathbf{w}}),$$

$${\bar{w}}_i=\{\begin{array}{ll}(w_i+w_{i+1})/2& i=1\\ (w_{i-1}+w_i+w_{i+1})/3& i=2,..,N-1\\ (w_{i-1}+2w_i)/3& i=N\end{array}$$

$$\mathcal{C}(a,a^{\prime })=\sqrt{C(a)C(a^{\prime })}\tilde{C}(a,a^{\prime })$$

Full text

1 Introduction

Since the discovery of the late time cosmic acceleration [1, 2], the physical mechanism underlying this accelerated phase has been an elusive one and still poses an open question for Cosmology. The simplest candidate that could drive this acceleration is a Cosmological Constant $\Lambda$, as in the $\Lambda$CDM model. Despite its success in describing cosmological observations [3], the $\Lambda$CDM model poses theoretical questions that still prompt the investigation of the nature of $\Lambda$, i.e. the value of this constant and the special moment of the Universe lifetime in which this starts to dominate over the other components [4, 5, 6] .

Furthermore, recent observations have reached a precision which started to highlight tensions between the measurements of $\Lambda$CDM parameters as inferred from different probes. The most striking one, appears in the values for the Hubble constant $H_{0}$ found with local measurements [7] and with Cosmic Microwave Background (CMB) observations [3]. While the former directly measure the current rate of expansion, the latter rely on the assumption of a specific cosmological model ($\Lambda$CDM) to extrapolate high redshift measurements to present time. A similar tension, although less statistically significant, is also found when measuring the growth of cosmological perturbations from a high or low redshift point of view, relying again on CMB [3] versus galaxy surveys [8, 9, 10]. Even though it is not excluded that these tensions could be driven by systematic errors, it could also be that they are due to the assumption of the $\Lambda$CDM model in the analysis of high redshift data.

Many alternatives to the $\Lambda$CDM have been put to test against observations (see e.g. [11, 12, 13, 14, 15, 16]) sometimes allowing to rule out specific theories [17, 18, 19, 20] , but more often leading only to constraints of their parameter space around the $\Lambda$CDM limit. Given the significant amount of available models and the difficulty in testing all of them against the data, a model independent approach is desirable and the first efforts in developing such an approach date back to almost twenty years ago, right after the discovery of the late acceleration phenomenon [21, 22, 23, 24]. A typical example is the so-called Chevallier-Polarski-Linder (CPL) parameterization [25, 26], a simple two parameters description of the time-dependence of the equation of state for the component that dominates at late times and sources the cosmic acceleration. We will generically refer to this component as Dark Energy (DE), and refer to its equation of state as $w_{\rm DE}$ (EoS). This could be an actual additional fluid contributing to the energy-momentum tensor, or, alternatively, result from modifications of gravity (MG). In either cases, generally $w_{\rm DE}$ is expected to be non-zero. The CPL parameterization has been used extensively in the literature to constrain the time evolution of DE (see e.g. [27, 3, 28, 29]), and also to quantify the ability of future experiments to shed light on the nature of DE [30]. Even though this approach allows to investigate the DE problem without restricting to a specific theory, it still relies on assumptions on the time evolution of $w_{\rm DE}$. An alternative approach is to reconstruct $w_{DE}$ constraining its value at different times. This so-called “non-parametric” approach has been often used in literature, and used to constrain the expansion history with data such as Type Ia SNe [31, 32, 33, 34] along with BAO [35, 36, 37, 38, 39, 40], but the amount of free parameters needed to reconstruct the EoS has often hindered its usefulness.

In this paper, we take one step further and combine theoretical priors into the non-parametric reconstruction of $w_{\rm DE}$, based on general requirements of physical viability [41]. These are derived for broad, but specific classes of theories. Hence, even though they are mild, including them in the reconstruction slightly limits the model independence. More importantly, the priors allow to introduce a, theoretically informed correlation between the values of $w_{\rm DE}$ at different times, aiding the reconstruction. Such an approach also guarantees that the reconstructed function will eventually correspond to a theoretically viable model, while it can be complicated to map a completely model independent reconstruction of the EoS to the one produced by any well behaved theory [42].

The paper is organized as follows. In Section \jnl@stylesec:npr we outline the reconstruction techniques we exploited to obtain a model independent $w_{\rm DE}$ from binned values constrained through observational data. In Section \jnl@stylesec:data_analysis_method we describe the general data analysis methodology employed, we describe the theoretical information coming from viability requirements for two classes of alternatives to $\Lambda$CDM (Quintessence and Horndeski), and we detail how this information is included in our analysis. We then present the results of this reconstruction in Section \jnl@stylesec:results, before drawing our conclusions in Section \jnl@stylesec:concl.

2 Non parametric reconstruction

We start discretizing the Dark Energy EoS, $w_{\rm DE}$, into several binned values, with the value of $w_{\rm DE}$ in each bin being a free parameter of our analysis. The evolution of $w_{\rm DE}$ in time can be equivalently expressed in terms of the redshift $z$ or of the scale factor $a$, which are connected to time through the Friedmann equations. In this paper we choose to discretize $w_{\rm DE}$ into bins equally spaced in the scale factor $a$ (with $w_{\rm DE}(a_{i})=w_{i}$ at the center of the $i$-th bin), in the interval $[a_{\rm min},1]$. The choice of binning the EoS in scale factor rather than in redshift, as well as the specific choice of the binning strategy, is due to the analysis method that we will detail further in this paper, in particular to the use of theoretical viability priors (see Section \jnl@stylesubsec:corr_prior). Indeed, accounting for theoretical considerations on the general behavior of $w_{\rm DE}(a)$ would set conditions on the reconstruction and on the correlation between the binned values of the function, e.g. motivating a correlation length $\xi$ between the $w_{i}$ and, consequently, the setting of the binning strategy [43]. We will describe the binning strategy more explicitly case by case in Section \jnl@stylesubsec:corr_prior .

Cosmological observables generally do not depend directly on the values of $w_{\rm DE}$ at each redshift, but rather on the evolution in time of this function between the observer and the measured source. Therefore, we need to join the binned values into a function of time that can be, e.g., integrated. Let us use the luminosity distance $d_{L}(a)$ inferred from observations of standard candles to illustrate this point. It can be expressed as

[TABLE]

with

[TABLE]

Clearly, if we are to fit a binned EoS to SNe data, we need to ensure a well defined evolution for $w_{\rm DE}$ by specifying how the function will behave within a certain bin and in moving from one to the other. Several reconstruction strategies have been used in past literature ranging from simple or smoothed step functions to more complex statistical tools (see e.g. [33, 37, 38, 34]); in our analysis we interpolate the binned $w_{i}$ values via two alternative techniques: a smoothed step function and Gaussian Process (GP), described in detail below. Indeed, once the binning properties for a specific case are given, we develop two different reconstructions, corresponding to the two different methods, and compare the results obtained.

2.1 Smoothed step function

The simplest choice to interpolate the binned values of $w_{\rm DE}$ is to use a step function. Defining $w_{\rm DE}(a_{i})=w_{i}$ as the values of the EoS at the center of each bin,

[TABLE]

where $\theta_{H}$ is the Heaviside function and $N$ is the number of bins. With this choice the EoS is constant within each bin and it has a sharp transition in its value when moving from one bin to the next (see blue line in Figure \jnl@stylefig:camb_binning).

Such a reconstruction can in principle cause numerical issues due to the fact that the function and its derivative are not well defined at the boundaries between bins. For this reason, we adopt a smoothed step function (see e.g. [44])

[TABLE]

where $N$ is the number of bins, $s$ is a smoothing factor used to control the slope of the transition from one bin to the adjacent ones. In this work we set $s=10$222The choice of this value is chosen by trial and error in order to find a smoothing parameter producing a smooth transition between neighbouring bins, but still preserving a constant value within the bins themselves. The value used here is the same used in previous reconstructions [44]; the resulting reconstruction can be seen in Figure \jnl@stylefig:camb_binning (green line).

2.2 Gaussian Process

In the smoothed step approach discussed above, the reconstructed function is assumed to be constant within each redshift bin; this can in principle bias the results obtained if the choice of the bin intervals is not done properly. To overcome such a problem, one can rely on more sophisticated techniques, e.g. Gaussian Process. A GP is defined as a collection of random variables , any finite number of which have a joint Gaussian distribution [45]. Considering the training points $\vec{a}=(a_{1}~{}a_{2}~{}a_{3})$, we can always think of the function $w_{\rm DE}(a)$ evaluated at these points as a vector and, at each point $a_{i}$, $w_{i}$ is a Gaussian random variable with mean $\mu(a_{i})$ and variance $\sigma_{a,i}^{2}$. The entire vector will then be modeled with a multivariate Gaussian distribution

[TABLE]

where $\mathcal{N}$ stands for Normal distribution and $C$ is the covariance matrix, with $C_{ij}=C(w_{i},w_{j})$ denoting the correlation between two points $a_{i}$ and $a_{j}$. We expect that the closer $a_{i}$ and $a_{j}$ are, the more the corresponding values of the EoS, $w_{i}$ and $w_{j}$, will be correlated. In order to interpolate the training points we use the python package sklearn, which provides several choices for the covariance $C$, and we work here with the Radial Basis Function (RBF) :

[TABLE]

where $\xi$ is the correlation length, such that

[TABLE]

Such a choice provides a Kernel for the GP which is both stationary, because it is a function of $a-a^{\prime}$, and isotropic, since it is a function of its module $|a-a^{\prime}|$. The reason why we chose this Kernel is that it is fully characterized by only one parameter, the correlation lenght $\xi$. Of course this results in a lower flexibility of our reconstruction, but it does not introduce new degeneracies associated to a high number of hyperparameters. Moreover, this kernel is infinitely differentiable, leading to the process being infinitely mean-square differentiable [45]. In this paper, we do not explore the dependence of our results on the choice of the GP Kernel, but rather leave this investigation for our ongoing work which extends to the reconstruction of the functions relevant for large scale structure.

Given a set of training points $(a_{i},w_{i})$, we want to obtain the value of the function at a point $a_{*}$, defined as $w_{*}$. Since we expect the function to be smooth, for a small variation of the $a$ variable we do not expect the function at that point to differ much from its values in the adjacent points.

Assuming that $w_{*}$ will be Gaussian distributed as well, i.e. $w_{*}\sim\mathcal{N}(\mu_{*},C(a_{*},a_{*}))$, where $C(a_{*},a_{*})$ is the self-covariance, the joint distribution will assume the form of

[TABLE]

After specifying the mean and correlation functions, considering now all the fitting points, the GP is defined as

[TABLE]

This reconstruction method does not require only the binned points $a_{i}$ and $w_{i}(a)$ as an input, but also specific choices for the mean and the correlation length of the Gaussian distribution. Usually, GPs are used on observational data where the training points are fixed (provided by the data), and therefore one can obtain $\mu$ and $\xi$ as the values providing the best fit of the reconstructed function to the observations (see e.g. [33, 46]). Here instead, while the $a_{i}$ are fixed, the $w_{i}$ values are free parameters of our analysis. Therefore we need to obtain different GP reconstructions for each sampled set of parameters $w_{i}$. In order to achieve this, we choose not to vary the hyperparameters of the Gaussian process, $\mu$ and $\xi$. Indeed we fix them to set values, where $\xi$ is directly obtained from the theoretical priors, imposed by the requirement that the reconstructed $w_{\rm DE}(a)$ can be eventually linked to a physically viable theoretical model. In Section \jnl@stylesubsec:corr_prior we will discuss in detail how specifying to classes of alternative models provides theoretical information on the correlation length, and we will provide the values used in the different cases.

This theoretical prior we will make use of will not provide however any information on the mean of the GP $\mu$. This means that a correct analysis would need to use $\mu$ as a free parameter. We point out however, that given our choice of datasets, with SN and BAO significantly helping CMB to constrain the redshift trend of $w_{\rm DE}$, we do not expect this function to be extremely different from the $\Lambda$CDM limit. Moreover, with a sensible choice of the bins used to reconstruct the function, we expect the choice of $\mu$ to not impact significantly the results: the GP will move the function towards the mean whenever the reconstruction is done at a scale factor that is more than a correlation length away from the training points. This can be an issue when using real data as training points, but in our case we decide the bin positions a priori, and our binning strategy is able to limit the impact of fixing the $\mu$ parameter.

As we will obtain our results using GP, but also the smoothed step function as an alternative reconstruction method, we will be able to verify these assumptions; the second method is not affected by the choice of the GP hyperparameters, therefore, when comparing the results, a significant discrepancy might imply that the choice of $\mu$ is affecting the reconstruction significantly. As such a discrepancy is not found (see Section \jnl@stylesec:results), we don’t investigate further the dependence of the GP reconstruction on these parameters, leaving it for a future dedicated work.

3 Data and analysis method

Given the binned $w_{\rm DE}$, a reconstruction method and the standard set of cosmological parameters, we can make predictions for the desired cosmological observables and constrain our parameter space against data. To this extent, we use the public code CAMB [47, 48], modifying it in such a way that we can provide as an input the binned values of $w_{\rm DE}$ and the chosen reconstruction method. We then sample the entire parameter space with the public Monte Carlo Markov Chain (MCMC) code CosmoMC [49]. For the cosmological parameters, we choose those of the minimal flat $\Lambda$CDM cosmology: the baryon and cold dark matter densities at present day, $\Omega_{b}h^{2}$ and $\Omega_{c}h^{2}$; the optical depth, $\tau$; the primordial power spectrum amplitude and tilt, $A_{s}$ and $n_{s}$, and the Hubble parameter $H_{0}$. Alongside the standard cosmological parameters we also sample the binned values of the EoS $w_{i}$. The prior range and bins positions in scale factors will depend on the specific cases investigated and will be specified in the following section; in general however, we assume that the EoS stays constant outside the binned interval, i.e. $w(a>a_{1})=w_{1}$ and $w(a<a_{N})=w_{N}$. We use flat priors on all the parameters.

We use the following datasets: the JLA (’Joint Light-curve Analysis’) dataset [50], that unifies Type Ia SNe observations of SDSS-II (Sloan Digital Sky Survey) [51] and SNLS (Supernova Legacy Survey) [52] collaborations, for a total of 740 Type Ia SNe up to redshift $z\sim 1$; the 6dFGS (6dF Galaxy Survey) [53] and SDSS Data Release 7 [54] for BAO and Planck 2015 data [55] for CMB. We will specify in the following subsection the prior ranges for $w_{i}$, as well as the different combinations of the data listed above, depending on the different classes of theories under consideration.

3.1 Theoretical correlation prior

The methods presented in Section \jnl@stylesec:npr allow to reconstruct the Dark Energy EoS as a function of redshift starting from its binned values; however the $w_{i}$ in each bin can in principle assume any value, with no relation to the values of the other bins. This does not take into account that viable theoretical models do not generally allow for extreme oscillations of this function [56], and therefore such a reconstruction can lead to a $w_{\rm DE}(a)$ that cannot be linked to any physically viable theoretical model. As already anticipated, we overcome this problem including a theoretically informed prior which account for the fact that a given $w_{i}$ will be correlated with the values assumed by the function in the adjacent bins. Following the approach of [43] we impose a correlation prior, which rescales the Likelihood function ($\mathcal{L}$) as

[TABLE]

where we assumed a Gaussian distribution $\mathcal{L}\propto\exp{[-\chi^{2}/2]}$. We define this prior as

[TABLE]

where $\mathbf{w}=\{w_{1},...,w_{N}\}$ is the vector containing the binned values of $w_{\rm DE}(a)$, $\mathcal{C}$ is a covariance matrix and $\mathbf{\bar{w}}$ is the vector composed by the expected values of $w_{\rm DE}$ in each bin. To avoid any specific choice for $\mathbf{\bar{w}}$, which could affect the final results, we define this as the mean value of the EoS in the bin under consideration and in the adjacent ones, i.e.

[TABLE]

where in the case of the last bin, we account for the fact that $w(a<a_{N})=w(a_{N})$. Including this prior distribution, means that we are assuming that our $w_{i}$ parameters are Gaussian distributed and their fluctuations around $\mathbf{\bar{w}}$ are described by the covariance matrix $\mathcal{C}$. Since we want this matrix to encode the theoretical correlation of the values of $w_{\rm DE}(a)$ at different redshifts, $\mathcal{C}$ needs to be obtained imposing conditions arising from the physical viability of theoretical models. We use therefore the results of [41], where such a correlation was obtained in different Dark Energy and Modified Gravity models, imposing a set of viability and stability conditions, e.g. requiring the avoidance of ghosts and gradient instabilities.

Following [41] we write the covariance matrix as

[TABLE]

with $C(a)$ an autocorrelation function and $\tilde{C}(a,a^{\prime})$ a correlation matrix, functions only of the scale factor.

We focus our analysis on the class of canonical single field Quintessence (from now on simply Quintessence), and on the broad class of Horndeski gravity. Despite somewhat limiting the model independence of our approach, we can still analyze a broad class of models, ensuring at the same time that the results we obtain will correspond to theoretical models which satisfy physical viability conditions.

We use $C(a)$ and $\tilde{C}(a,a^{\prime})$ obtained by [41], where the numerical correlations found for different classes of DE and MG models were encoded in the fitting formulas shown in Table \jnl@styletab:matrix_tab, where for $\tilde{C}$ two possible choices are taken into account, i.e. an exponential and a CPZ parameterization. Both these choices are built in such a way that if $\delta a=a_{i}-a_{j}$ (or $\delta\ln(a)=\ln(a_{i})-\ln(a_{j})$) is much higher than the correlation length $\xi$ then the correlation tends to zero. The prior that best minimizes the residuals in both the Quintessence and Horndeski cases is the exponential one [41].

We want to stress here that having a correlation length set by theoretical requirements significantly helps in the choice of the binning strategy, since the correlation length defines a number of effective degrees of freedom $N_{\rm{eff}}=(a_{max}-a_{min})/\xi$; as long as the number of bins $N$ satisfies $N>N_{\rm{eff}}$, the dependence of the reconstruction on the number of bins is negligible [43].

The specific analysis strategy will depend on the class of models under consideration. Quintessence models modify the background expansion of the Universe with respect to $\Lambda$CDM without modifying the equations for the evolution of cosmological perturbations; therefore we can use both background (SN and BAO) and perturbation (CMB) data. On the contrary, for Horndeski models the use of CMB data would require us to consider also their impact on the growth of cosmological structures, eventually jointly binning the phenomenological functions $\mu$ and $\Sigma$ [27, 57]. Since the reconstruction of these two functions is beyond the scope of this paper, we limit our analysis in this case to background data only.

The two classes of models also have different requirements for the prior range:

• •

Quintessence: the EoS lies in the region above the Phantom divide ($w=-1$), hence we impose $w_{\rm DE}(a)\geq-1$. To this purpose we sample the $w_{i}$ parameters in the range $[-1,0]$. The reconstruction of $w_{DE}(a)$ for Quintessence is done using both the CPZ and exponential correlation priors (Table \jnl@styletab:matrix_tab), to highlight how they differently affect the inferred $w_{i}$ parameters. Given the theoretical correlation length for this case ($\xi_{exp}=0.6\sim\xi_{CPZ}=0.7$) the $w_{i}$ parameters are defined in the bins $\vec{a}=(0.85,~{}0.7,~{}0.55,~{}0.4,~{}0.25,~{}0.1)$ in order to satisfy the requirement $N>N_{\rm{eff}}$.

• •

Horndeski: the EoS has no particular restriction in the values that it can assume, hence $w_{i}$ is sampled in the range $[-3,0]$. The reconstruction for the Horndeski case is obtained using the exponential prior (Table \jnl@styletab:matrix_tab). In this case the Planck CMB data are not included since we are not modelling the impact of this class of models on the evolution of perturbations, hence we only rely on background data, as Type Ia SNe and BAO. Concerning the latter, it is known that BAO depend on the assumption of a fiducial cosmology in order to compare data and theoretical predictions. We assume here that a change of the perturbations evolution as produced by Horndeski models does not affect significantly the use of this probe. This assumption is however deeply under scrutiny at present, and we refer the reader to recent works that investigate this problem [58, 59]. Also in this case, the choice of the scale factor range in which we define our bins is made in such a way to satisfy the requirement $N>N_{\rm{eff}}$; as the correlation length in this case is $\xi_{exp}=0.3$ we use $\vec{a}=(0.9,~{}0.8,~{}0.7,~{}0.6,~{}0.5,~{}0.4)$. In this case, since we are not using Planck data, the redshift range in which data are available is smaller than in the previous case, therefore we limit our reconstruction to $a_{6}=0.4\ (z_{6}=1.5)$.

We apply this prior information to our analysis both when using the smoothed step and GP reconstruction. It is worth stressing that in the latter case, the GP and prior both exploit a correlation between the values of $w_{DE}(a)$ at different scale factors. However, while GP returns the value of $w_{\rm DE}$ at a given $a_{*}$ given the $w_{i}$, without acting on the $w_{i}$ itself, the correlation prior is related to the relative values of the different $w_{i}$, penalizing those configurations of binned values which do not satisfy the viability conditions. In this sense, the GP and correlation priors are complementary and their combined use is possible.

4 Results

In this Section we present the results obtained following the analysis method described in the previous two Sections, discussing separately the reconstruction within the Quintessence realm and that for the more general Horndeski class.

4.1 Quintessence

We report the results for Quintessence in Table \jnl@styletab:quint_res. The results are obtained both with and without the inclusion of the theoretical prior (exponential and CPZ), for comparison. In all three cases, we find only upper limits on the $w_{i}$ parameters, with the lower bound $w=-1$ included within the $1-\sigma$ bound, thus indicating that these results are compatible with $\Lambda$CDM.

Figure \jnl@stylefig:quint_binned shows the results obtained reconstructing the DE EoS via the smoothed step function. We point out that the mean values of the reconstruction deviate from the $\Lambda$CDM limit given the requirement of a hard bound for the posterior at $w_{i}=-1$, which however does not prevent $\Lambda$CDM from being a good fit for this reconstruction.

The comparison of the three cases, bin by bin, can be quantified as a percentage difference of the mean values inferred adding the exponential (exp) and CPZ priors with respect to the mean values obtained with no theoretical prior. The main differences associated to the choice of a particular prior are related to the behavior of $w_{4}$ (exp: $-1.32\%$, CPZ: $-0.44\%$), $w_{5}$ (exp: $-3.3\%$, CPZ :$-2.01\%$) and $w_{6}$ (exp: $-4.3\%$, CPZ: $-2.47\%$) mean values and $1\sigma$ confidence levels, while $w_{1}$ (exp: $-0.249\%$, CPZ: $0\%$), $w_{2}$ (exp: $+0.214\%$, CPZ: $+0.107\%$) and $w_{3}$ (exp: $0\%$, CPZ: $+0.108\%$) do not differ significantly case by case. The larger impact of the correlation priors in the last three bins can be explained with the fact that these redshifts contain much less data; therefore the constraining power of the prior is comparable to that of the data.

The prior however, does not impact significantly the behavior of the function, showing how the requirement $w_{i}\geq-1$ already satisfies the physical viability conditions imposed through the correlation prior.

In order to test the dependence of the results on the reconstruction method, we perform the analysis for the Quintessence class, both with Planck+BAO+SN only and with the inclusion of the CPZ prior, using the Gaussian Process reconstruction method rather then the smoothed bins one. In this case, we fix the correlation length $\xi$ of the GP to that of the theoretical prior. Furthermore, we choose to force $w_{\rm DE}(a)$ to remain constant at the value of the first redshift bin $w_{1}$ for $a>a_{1}$.

Table \jnl@styletab:results_GP_Quint shows the results obtained on the cosmological and reconstruction parameters in the Planck+BAO+SN and Planck+BAO+SN+prior cases, while in Figure \jnl@stylefig:Quint_GP_single we show the reconstruction of the $w_{\rm DE}(a)$ function given the mean values of $w_{i}$. We notice how the results are compatible with the smoothed bins reconstruction, but with the prior affecting more significantly than before the earliest redshifts of the reconstruction: the mean values of the reconstruction obtained adding the prior respectively differ from the $w_{1}$, $w_{2}$ and $w_{3}$ means obtained with only the datasets by $+4.2\%$, $-1.5\%$ and $+0.86\%$. These results are emphasized by the comparison shown in Figure \jnl@stylefig:Quint_cfr_GP_smoothed.

In addition to these results obtained in the single field minimally coupled Quintessence case, a model implying a limitation $w_{i}\geq-1$, we show in Figure \jnl@stylefig:genquint the reconstruction obtained with the Planck+BAO+SN dataset when this requirement is not imposed, labelling it Generalized Quintessence. EoS crossing the phantom divide can be obtained by models with non minimal coupling or with kinetical braiding [60, 61, 62], as well as by scalar-tensor gravity [63]. Notice however that in general these models might introduce also modifications of perturbation evolution with respect to $\Lambda$CDM, which are not included in our reconstruction. We stress therefore that the results presented here assume that even with such a generalized expansion history, no modification of perturbations is produced. Moreover, as we have no available theoretical information on this reconstruction case, we do not make use of theoretical priors. We find that also in this case the two reconstruction methods are compatible with each other, obtaining the same trend in redshift for the EoS within the $1\sigma$ bound. The $\Lambda$CDM limit $w=-1$ is generally compatible with the results found within $1\sigma$ except for the bin at $a=0.25$, where the EoS is constrained to be below $w=-1$ at $1\sigma$. Nevertheless, no evidence for deviation from standard cosmological constant is found, as $\Lambda$CDM is always compatible within the $2\sigma$ limit.

4.2 Horndeski class of MG models

In Table \jnl@styletab:results_Horn we report the results obtained reconstructing the EoS via smoothed step function and GP in the Horndeski case, for which the $w_{i}$ parameters are free to vary within the range $[-3,0]$ and the exponential correlation prior is imposed. As already mentioned, in this case we restrict to background data, and do not use CMB measurements. In both cases the last bin considered $w_{6}$ is unconstrained due to the fact that very few data are available at this redshift ($z\sim[1,1.5]$). Because of the correlation that we have introduced between the binned values of $w_{\rm DE}$, this lack of constraining power propagates also to the constraints at higher scale factors (lower redshifts). For the other bins, we obtain compatible constraints from the two reconstruction methods, in agreement within $1\sigma$ with each other and with a $\Lambda$CDM expansion history. Notice that in this case, with $w_{\rm DE}(a)$ able to take values below $-1$ and with no CMB data, we do not have constraining power on $H_{0}$.

In Figure \jnl@stylefig:Horn_cfr_gpsmoo we show the reconstruction of $w_{DE}(a)$ given the inferred $w_{i}$ using both the smoothed and Gaussian process reconstruction methods. The function tends to decrease over redshift, independently of the reconstruction method, with the mean reconstructed function slightly closer to the $\Lambda$CDM value in the GP case for intermediate redshifts. We stress however that both reconstruction are compatible between each other and with the $\Lambda$CDM limit at approximately $1\sigma$.

The inclusion of CMB data also in the Horndeski case would increase the amount of information and potentially improve the constraints achievable with this method; however, as already pointed out, this would require the reconstruction of the $\mu(a)$ and $\Sigma(a)$ encoding the departures from the standard $\Lambda$CDM growth of cosmological perturbations. The possible degeneracies between the effects of these two functions and those of a varying $w_{\rm DE}(a)$ can therefore limit the improvement in the achievable constraints when CMB data are included. A full analysis of this scenario would therefore require the use of cosmological data able to disentangle the background and perturbation modifications, e.g. Cosmic Shear or Galaxy Clustering data from Large Scale Structure surveys [64, 65, 9].

5 Conclusions

We have reconstructed the Equation of State (EoS) of Dark Energy, $w_{\rm DE}(a)$, from the latest cosmological data with two different techniques: one in which the EoS is assumed constant within the range of each specified bin in scale factor, with smooth transitions between binned values; and one in which the function is reconstructed at each value of $a$ from its binned values using a Gaussian Process reconstruction. We modified the public code CAMB to produce predictions on cosmological observables given the reconstructed $w_{\rm DE}(a)$ and we compared these predictions with currently available data both for background (SN and BAO) and for perturbations (CMB) observables. Alongside the observational data, we took into account the contribution of theoretical conditions on the physical viability of the reconstructed EoS; this was possible with the inclusion of a correlation prior between the binned values of $w_{\rm DE}(a)$, obtained in two classes of models: single field Quintessence and Horndeski gravity [41].

The analysis of these two classes differs both in the binning strategy for $w_{\rm DE}(a)$, which is defined by the different correlation lengths $\xi$, and in the observables used to constrain the parameters: in the Quintessence case, the reconstruction of $w_{\rm DE}(a)$ was sufficient to fully characterize the modification on background and perturbations evolutions with respect to $\Lambda$CDM and we can therefore use all the data from BAO, SN and CMB. In the Horndeski case instead, the description of perturbations would have required also to reconstruct the $\mu(a)$ and $\Sigma(a)$, and therefore we limited our analysis to the background data coming from BAO and SN observations. Moreover, while for Horndeski the binned values $w_{i}$ were allowed to vary in the full prior range $[-3,0]$, this range was limited in Quintessence to $[-1,0]$ as these models cannot cross the phantom divide $w=-1$ without developing ghost instabilities.

In our results for the Quintessence case, we found that $\Lambda$CDM $w=-1$ is compatible with the reconstructed $w_{\rm DE}(a)$ within $1\sigma$. Moreover we found that the theoretical prior is not affecting the results significantly, due to the limitation of our analysis to the non-phantom part of the parameter space ($w_{i}\geq-1$) which already satisfies the most restrictive of the physical viability conditions [42]. The results of the two reconstruction methods are in agreement with each other, highlighting how for a high enough number of bins, the smoothed bins reconstruction is able to reproduce sufficient variations of $w$ in redshift to fit cosmological data.

In the Horndeski case, we also found that $w=-1$ is compatible with the reconstructed EoS within $1\sigma$. The smoothed bins and GP reconstruction methods are found to be compatible also in this case, although less than in the Quintessence analysis; this is possibly due to the fact that the oscillatory behavior hinted by the reconstruction could require a higher number of bins to be reproduced in the smoothed bins case, or that a different GP Kernel would be more appropriate for such behavior. Tighter constraints on the EoS, that would allow to quantify the impact of the theoretical prior also for this class of models, could be obtained including CMB data in the analysis, which requires to simultaneously reconstruct the phenomenological functions describing deviations from $\Lambda$CDM in clustering and lensing, respectively $\mu$ and $\Sigma$. As we discussed, given the degeneracies between the effects of $w$, $\mu$ and $\Sigma$, additional observables would be needed for this kind of analysis, e.g. coming from observations of Large Scale Structures.

Acknowledgments

We want to thank Simone Peirone, Marco Raveri and Robert Crittenden for useful discussions. MM acknowledges support from the D-ITP consortium, a program of the NWO that is funded by the OCW. AS acknowledges support from the NWO and the Dutch Ministry of Education, Culture and Science (OCW), and from the D-ITP consortium, a program of the NWO that is funded by the OCW.