Priors on the effective Dark Energy equation of state in scalar-tensor theories

TL;DR

This paper derives theoretical prior covariances for the dark energy equation of state in scalar-tensor theories, aiding non-parametric reconstructions and reducing bias in cosmological data analysis.

Contribution

It introduces a method to compute prior covariances for w(z) based on scalar-tensor theories, including viability conditions, enhancing cosmological parameter estimation.

Findings

Prior favors tracking behaviors in scalar-tensor theories

Covariance matrices enable bias-reducing priors for w(z)

Method applicable to various parametrizations of dark energy

Abstract

Constraining the Dark Energy (DE) equation of state, w, is one of the primary science goals of ongoing and future cosmological surveys. In practice, with imperfect data and incomplete redshift coverage, this requires making assumptions about the evolution of w with redshift z. These assumptions can be manifested in a choice of a specific parametric form, which can potentially bias the outcome, or else one can reconstruct w(z) non-parametrically, by specifying a prior covariance matrix that correlates values of w at different redshifts. In this work, we derive the theoretical prior covariance for the effective DE equation of state predicted by general scalar-tensor theories with second order equations of motion (Horndeski theories). This is achieved by generating a large ensemble of possible scalar-tensor theories using a Monte Carlo methodology, including the application of physical viability conditions. We also separately consider the special sub-case of the minimally coupled scalar field, or quintessence. The prior shows a preference for tracking behaviors in the most general case. Given the covariance matrix, theoretical priors on parameters of any specific parametrization of w(z) can also be readily derived by projection.

Tables

Table 1. Table 1: Summary of the auto-correlation and correlation function fits.

	Best-fit auto-corr	residuals
Quintessence	$0.03+0.3\exp \left[6.5\left(a-1\right)\right]$	$0.01$
GBD	$0.05+0.8\exp \left[1.8\ln a\right]$	$0.007$
Horndeski	$0.05+0.8\exp \left[2\ln a\right]$	$0.007$
	Best-fit corr	residuals
Quintessence	$\exp \left[-{\left(\left|\delta a\right|/0.7\right)}^{1.8}\right]$	$9$
GBD	$\exp \left[-{\left(\left|\delta \ln a\right|/0.3\right)}^{1.3}\right]$	$6$
Horndeski	$\exp \left[-{\left(\left|\delta \ln a\right|/0.3\right)}^{1.2}\right]$	$6$
	Best-fit corr (fixed CPZ)	residuals
Quintessence	${\left(1+{\left(\left|\delta a\right|/0.6\right)}^2\right)}^{-1}$	$11$
GBD	${\left(1+{\left(\left|\delta \ln a\right|/0.2\right)}^2\right)}^{-1}$	$12$
Horndeski	${\left(1+{\left(\left|\delta \ln a\right|/0.2\right)}^2\right)}^{-1}$	$13$

Equations

$$\mathcal{S}=$$

$$+$$

$$\left(1+\Omega +\frac{1}{2}a{\Omega }^{\prime }\right)\frac{dy}{d\ln a}$$

$$w_{\mathrm{DE}}\equiv \frac{P_{DE}}{{\rho }_{DE}}=\frac{-2\dot{\mathcal{H}}-{\mathcal{H}}^2-P_m{a}^2/m_0^2}{3{\mathcal{H}}^2-{\rho }_m{a}^2/m_0^2},$$

$$w_{\mathrm{DE}}(a)=\frac{\Lambda }{2c-\Lambda }=\frac{K-P}{{\rho }_{\mathrm{DE}}(a)}$$

$$w_{\mathrm{DE}}(a)=\frac{\Lambda -\Omega P_m}{2c-\Lambda -\Omega {\rho }_m},$$

$$f(a)=\sum _{n=0}^N\frac{{\alpha }_n}{n!}{\left(a-a_0\right)}^n,$$

$$f(a)=\sum _{n=0}^N{\alpha }_n{\left(a-a_0\right)}^n.$$

$$f(a)=\frac{{\sum }_{n=0}^N{\alpha }_n{\left(a-a_0\right)}^n}{1+{\sum }_{m=1}^M{\beta }_m{\left(a-a_0\right)}^m},$$

$$P(\vec{w})=\sum _{\mathcal{M}}P\left(\vec{w}|\mathcal{M}\right)P\left(\mathcal{M}\right),$$

$${\bar{w}}_i$$

$$C_{ij}$$

$${\mathcal{C}}_{ij}=\frac{C_{ij}}{\sqrt{C_{ii}{C}_{jj}}}.$$

$$C(a,a^{\prime })\equiv ⟨[w_{\mathrm{DE}}(a)-{\bar{w}}_{\mathrm{DE}}(a)][w_{\mathrm{DE}}(a)-{\bar{w}}_{\mathrm{DE}}(a^{\prime })]⟩.$$

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

$$w_{\mathrm{DE}}(a)\approx w_0+w_a(1-a),$$

$$w_a=-{\frac{d{w}_{\mathrm{DE}}}{da}|}_{a=1}.$$

$$\mathcal{C}(x,y)=\frac{1}{1+{\left(|x-y|/\xi \right)}^n},$$

$$\mathcal{C}(x,y)=\exp \left[-{\left(|x-y|/\xi \right)}^n\right],$$

$$\mathcal{P}(\{w_i\},\vec{p})\propto \mathcal{L}(\{w_i\},\vec{p})\times {\mathcal{P}}_{\mathrm{prior}}(\vec{p})\times {\mathcal{P}}_{\mathrm{prior}}(\{w_i\}),$$

$${\mathcal{P}}_{\mathrm{prior}}(\{w_i\})\propto e^{-{\sum }_{ij}(w_i-{\bar{w}}_i)[C^{-1}{]}_{ij}(w_j-{\bar{w}}_j)/2}.$$

$${\mathcal{H}}^2=$$

$$\dot{\mathcal{H}}=$$

$$\frac{c{a}^2}{m_0^2}=\frac{3}{2}\left(1+\Omega +a{\Omega }^{\prime }\right){\mathcal{H}}^2-\frac{1}{2}\frac{a^2{\rho }_m}{m_0^2}+\frac{1}{2}\frac{\Lambda a^2}{m_0^2}.$$

$$\left(1+\Omega +\frac{1}{2}a{\Omega }^{\prime }\right)\frac{dy}{d\ln a}$$

$$\frac{c_0}{m_0^2}=\frac{3}{2}{\mathcal{H}}_0^2\left(1+{\Omega }_0+{\Omega }_0^{\prime }-{\Omega }_m^0-{\Omega }_{\Lambda }^0\right),$$

$${\Omega }_m^0\equiv \frac{1}{3}\frac{{\rho }_m^0}{{\mathcal{H}}_0^2{m}_0^2};{\Omega }_{\Lambda }^0\equiv -\frac{1}{3}\frac{{\Lambda }_0}{{\mathcal{H}}_0^2{m}_0^2}.$$

$${\mathcal{H}}^2$$

$$\dot{\mathcal{H}}$$

$$w\equiv \frac{P_{DE}}{{\rho }_{DE}}=\frac{-2\dot{\mathcal{H}}-{\mathcal{H}}^2-P_m{a}^2/m_0^2}{3{\mathcal{H}}^2-{\rho }_m{a}^2/m_0^2}.$$

Full text

I Introduction

The cosmological constant $\Lambda$ is the simplest form of Dark Energy (DE) capable of driving the observed accelerating cosmic expansion Riess et al. (1998); Perlmutter et al. (1999). In General Relativity (GR), the bare (geometric) cosmological constant is a free parameter. However, $\Lambda$ also includes a constant vacuum energy density that receives contributions from the zero-point fluctuations of known particles Weinberg (1989). To match the observed value, the geometric contribution to $\Lambda$ needs to be separately fine-tuned to compensate for the contribution of each different particle type individually, making the fine-tuning technically unnatural Burgess (2015).

The extraordinary degree of tuning required to explain the observed value of $\Lambda$ within GR has spurred a broad search for mechanisms that can naturally set $\Lambda$ to a small value. A common proposal involves finding ways to prevent the vacuum from gravitating altogether, as in the so-called “degravitation” Dvali et al. (2007), “self-tuning” Charmousis et al. (2012) or “sequestering” Kaloper and Padilla (2014) scenarios, which would solve the naturalness problem at the cost of invoking some new mechanism – generally a new field or force – to cause the cosmic acceleration. These closely-related issues have motivated studies of alternative models of DE and modifications of gravity (MG). While not necessarily successful in resolving the fine-tuning problem themselves, many of these models can account for the acceleration, often predicting an effective DE that is dynamical. As the quest for a complete understanding of cosmic acceleration continues, it is of significant interest to know the extent to which DE can be dynamical and still remain consistent with observations.

Given the values of the Hubble parameter today, $H_{0}$, and the current matter density fraction, $\Omega_{\rm M}$, an arbitrary expansion history can be reproduced by assuming a flat Friedmann-Robertson-Walker (FRW) universe with a DE component that has an equation of state (EoS) $w_{\rm DE}$ with an appropriately chosen dependence on redshift Caldwell et al. (1998). Constraining $w_{\rm DE}(z)$, which we define here as the effective EoS of all non-dust contributions to the Friedmann equation at late times, is one of the primary science goals of ongoing and future surveys of large scale structure. A detection of $w_{\rm DE}(z)\neq-1$ would be strong evidence of new gravitational physics, and would constitute a potentially vital clue in understanding the source of cosmic acceleration.

Fitting a constant $w_{\rm DE}$ to current data results in good agreement with $-1$ Ade et al. (2016), but such fits would have missed subtle variations in $w_{\rm DE}(z)$, especially if the average happened to be close to $-1$. By the same token, using the Chevallier-Polarski-Linder (CPL) parameterization Chevallier and Polarski (2001); Linder (2003), or any particular parametric form of $w_{\rm DE}(z)$, is prone to biasing the outcome Sahni and Starobinsky (2006). Forecasts using Principle Component Analysis (PCA) Crittenden et al. (2009); Albrecht et al. (2009) have shown that future surveys like Euclid and LSST will be able to constrain several eigenmodes of $w_{\rm DE}(z)$, justifying the use of more flexible parameterizations.

A relatively general way to proceed is to use a piecewise representation, defining $w_{\rm DE}(z)$ in terms of its values in discrete bins in $z$. In the limit of many bins, fitting such a parametrization to data would avoid any bias, as any arbitrary $w_{\rm DE}(z)$ could be reproduced. In practice, however, if more than a few bins are used, many of them will be left unconstrained by the data, with values in neighboring bins effectively becoming degenerate. A possible way to remedy this problem is to introduce correlations between the bins by supplementing the likelihood with a prior covariance Crittenden et al. (2012). In this approach, the bias in the reconstructed $w_{\rm DE}(z)$ can be controlled for a given survey by using a PCA forecast. Specifically, as demonstrated in Crittenden et al. (2012); Zhao et al. (2012), one can tune the strength of the prior to avoid the bias in the reconstructed evolution of $w_{\rm DE}(z)$ on timescales larger than a given correlation scale. A similar effect can be achieved using Gaussian Processes Shafieloo et al. (2012); Seikel et al. (2012), which model correlations in the value of $w_{\rm DE}$ as a smooth function of redshift. Another approach, taken in Shafieloo et al. (2012); Daly and Djorgovski (2003); Shafieloo et al. (2006), attempts to derive $w_{\rm DE}(z)$ directly from data, by taking derivatives of the luminosity distances $d_{\rm L}(z)$. While seemingly model-independent, such methods require smoothing of the (noisy) data before derivatives can be taken, which effectively amounts to adopting a correlation prior.

Since any attempt to constrain $w_{\rm DE}(z)$ will require a prior in one form or the other, one would like to develop reconstruction tools based on theoretically motivated priors that are explicit and easy to interpret. Scalar fields (fundamental or effective) are ubiquitous in theoretical cosmology and particle physics, and are well-suited for representing a dynamical DE. A natural way to construct a theoretically-consistent, yet not-too-model-specific prior, is by generating a representative ensemble of scalar field models. Such approach was taken, for example, in Crittenden et al. (2007); Marsh et al. (2014), where priors on the CPL parameters $w_{0}$ and $w_{a}$ were derived based on an ensemble of very general quintessence DE models, i.e. minimally-coupled scalar fields Ratra and Peebles (1988); Caldwell et al. (1998). In this paper, we go beyond minimally-coupled scalars, and consider the most general scalar-tensor models that are “stable”, in the sense that they lack ghosts, gradient instabilities, and other pathologies. To generate the ensemble of models, we use the so-called “unifying” or “effective field theory” (EFT) approach to DE and MG Gubitosi et al. (2013); Bloomfield et al. (2013); Gleyzes et al. (2013); Bloomfield (2013); Piazza and Vernizzi (2013); Gleyzes et al. (2015a, b). Rather than working with $w_{0}$ and $w_{a}$, we derive the prior covariance for a binned $w_{\rm DE}(z)$, which makes it possible to project onto other parameterizations, and allows for reconstruction of $w_{\rm DE}(z)$ using the methods introduced in Crittenden et al. (2012); Zhao et al. (2012).

We generate the prior covariance of $w_{\rm DE}(z)$ for Horndeski theories, i.e. the general class of scalar-tensor theories with up to second order equations of motion Horndeski (1974); Deffayet et al. (2011). We also separately consider the widely studied sub-classes of the minimally couple scalar field, or quintessence Ratra and Peebles (1988); Caldwell et al. (1998), and models with the canonical form of the scalar field kinetic energy, i.e. the generalized Brans-Dicke (GBD) models Brans and Dicke (1961); Carroll (2004). A non-minimally coupled scalar field mediates a fifth force which is tightly constrained in the Solar System. We will not concern ourselves with satisfying such constraints, as the models may include a screening mechanism that suppresses the fifth force in dense environments. For our purposes, a scalar-tensor theory is viable if it has stable cosmological perturbations.

The paper is organized as follows. In Sec. II we introduce the EFT description of the cosmological background and linear perturbations in scalar-tensor theories and discuss a couple of representative forms of $w_{\rm DE}(z)$. We detail our method for sampling the space of viable Horndeski models in Sec. III. In Sec. IV we present our results for the prior probability distribution for $w_{\rm DE}(z)$ and its covariance. We also provide an analytical form that accurately represents the correlation of $w_{\rm DE}(z)$ at any two points in $z$. We conclude with a summary in Sec. V.

II Modeling viable scalar field

Dark Energy

We focus on the broad class of scalar-tensor models of gravity with second order equations of motion. The corresponding action in $(3+1)$ dimensions was derived by Horndeski Horndeski (1974), and later rediscovered in the context of generalized Galileons Deffayet et al. (2011). For the purpose of our analysis, rather than working with the full covariant action of Horndeski gravity, it is convenient to employ the unifying framework, or effective theory approach (EFT), formulated in Gubitosi et al. (2013); Bloomfield et al. (2013) and further developed in Gleyzes et al. (2013); Bloomfield (2013); Piazza and Vernizzi (2013); Gleyzes et al. (2015a, b). This allows us to model the background evolution and linear perturbations in a model independent way, in terms of a handful of free functions of time. The relevant EFT action reads

[TABLE]

where $R$ is the four-dimensional Ricci scalar, $\delta g^{00}$, $\delta{K}{{}^{\mu}_{\nu}}$, $\delta{K}{{}^{\mu}_{\mu}}$ and $\delta\mathcal{R}$ are respectively the perturbations of the upper time-time component of the metric, the extrinsic curvature, its trace, and the three dimensional spatial Ricci scalar of the constant-time hypersurfaces. The six functions, $\{\Omega,\Lambda,c,M_{2}^{4},\bar{M}_{1}^{3},\bar{M}_{3}^{2}\}$, are arbitrary functions of time allowed by the breaking of the time-diffeomorphism invariance, to which we refer to as the “EFT functions”. Finally, $S_{m}$ is the action for all matter fields $\chi_{m}$ minimally coupled to the metric $g_{\mu\nu}$. Action (1) is built in the unitary gauge, where the additional scalar field is absorbed into the metric, and represents an expansion to quadratic order in perturbations around the flat Friedmann-Robertson-Walker (FRW) universe. A few extra terms were present in the original formulation of EFT of dark energy Gubitosi et al. (2013); Bloomfield et al. (2013), and different generalizations of (1) have recently been studied in the literature Kase and Tsujikawa (2015); Frusciante et al. (2015, 2016) to include models outside the Horndeski class.

Perturbations in any specific Horndeski model can be mapped onto the action (1), with the correspondence between the EFT functions and the functions that appear in the full Lagrangian of the model given in Bloomfield (2013). In the absence of a preferred scalar-tensor theory, one can adopt the agnostic point of view, and treat the EFT functions as free functions of time. We note that the GBD theories are represented by terms in the first line of (1), while quintessence only requires specifying $\Lambda(\tau)$, with all other functions set to zero.

One of the advantages of using the EFT formulation is the ability to separate the terms that affect the background evolution from those that only concern the perturbations. In particular, only three of the EFT functions, $\Omega$, $c$ and $\Lambda$, play a role in determining the dynamics of the background, hence we will refer to them as the background EFT functions. Furthermore, as detailed in Appendix A, only two out of these three functions are sufficient to fix the background dynamics. In our approach, we specify $\{\Omega(a),\Lambda(a)\}$, and find the Hubble parameter, ${\mathcal{H}}(a)$, by solving a differential equation:

[TABLE]

where $y\equiv{\mathcal{H}}^{2}$, and the prime indicates differentiation with respect to the scale factor. Given the solution for ${\mathcal{H}}(a)$, the effective DE EoS is defined via

[TABLE]

where $\rho_{m}$ and $P_{m}$ are the combined energy density and the pressure of all particle species, and the over-dot denotes a conformal time derivative. The full details of this procedure are given in Appendix A. By sampling from a broad range of possible $\{\Omega(a),\Lambda(a)\}$ functions, as described in Sec. III, we are able to generate a representative ensemble of possible expansion histories in single field DE models. We note that this procedure is different from the so-called designer approach used in e.g. Hu et al. (2015), in which one provides ${\mathcal{H}}$ (and $\Omega$), and solves for $\{c,\Lambda\}$.

Since we are interested in reproducing the expansion histories of theoretically viable Horndeski models, we sample the functions $\{M_{2}^{4},\bar{M}_{1}^{3},\bar{M}_{3}^{2}\}$ along with $\{\Omega(a),\Lambda(a)\}$. While the former affect only the perturbations and, hence, would appear not to matter for the prior on $w_{\rm DE}(z)$, they do in fact play a role in determining whether a given background solution corresponds to a stable Horndeski theory. By perturbing the action (1) around a given background, one can derive constraints on combinations of EFT functions and their derivatives that exclude instabilities. Specifically, in our analysis, we impose the no-ghost and no-gradient-type instability conditions for scalar and tensor modes Frusciante et al. (2016). We implement the procedure for solving for the background discussed above in EFTCAMB111http://eftcamb.org Hu et al. (2014a); Raveri et al. (2014) and use its stability module to filter out models with ghost and gradient type instabilities. We find that imposing these conditions leads to the model acceptance rate $\sim 50$% for quintessence, $\sim 10$% for GBD and $\sim 1$% for Horndeski, thus removing a noticeable fraction of the parameter space. The details on technical implementation of the stability conditions can be found in Hu et al. (2014b).

II.1 Two simple case studies

To gain intuition into the expected behavior of $w_{\rm DE}(z)$ in scalar-tensor theories, let us examine the effect of varying $\Lambda$ and $\Omega$.

First, we consider a model in which all EFT functions except $\Lambda$ are set to their $\Lambda$CDM values. This choice corresponds to the minimally coupled quintessence models, with $\Lambda$ representing the Lagrangian of the scalar field, or the difference between its kinetic and potential energy densities. For our illustration, we take $\Lambda(a)=\Lambda_{0}\left(1+\lambda(a)\right)$, with $\lambda(a)=-0.1\,a$, and solve the background equations. Panel (a) of Fig. 1 shows the DE EoS $w_{\rm DE}(z)$, while panel (b) shows the matter, radiation and DE energy fractions as a function of redshift. For this model, we see that $w_{\rm DE}\approx 0$ during the matter and radiation eras, before evolving towards $-1$ at low redshift. This behavior is quite generic, and can be understood as follows. The effective DE EoS is given by

[TABLE]

where $K$ and $P$ are the quintessence kinetic and potential energy densities, and $\rho_{\rm DE}\equiv 3m_{0}^{2}H^{2}-\rho_{M}-\rho_{r}$ is the total energy density excluding matter and radiation. In typical quintessence models, the field is slowly rolling at late times (low redshifts), with $P>K$, but retains a small non-zero kinetic term, hence $w_{\rm DE}\gtrsim-1$. The energy density of the scalar field with $w_{\rm DE}\geq-1$ increases with redshift, unless the kinetic energy of the field is tuned to be exactly zero (which would be the case of the cosmological constant). This drives $w_{\rm DE}\to 0$ at high redshift, as seen in Fig. 1.

Next, let us consider a model in which both $\Omega$ and $\Lambda$ can vary in time. For illustration, we take $\Omega(a)=0.2(a-1)/[1+0.7(a-1)]$ and $\lambda(a)=0.3+1.0(a-1)$. As can be seen from Fig. 1, there are two key differences from the $\Lambda$-only model. Firstly, $w_{\rm DE}$ has a tracking behavior at early times, with $w_{\rm DE}\sim 0$ during the matter era and $w_{\rm DE}\sim 1/3$ in the radiation era. Secondly, $w_{\rm DE}$ crosses $-1$ at late times, when dark energy starts becoming important. Both features are quite representative of non-minimally coupled models. To understand the tracking behavior, observe that changing the value of $\Omega$ changes the gravitational coupling, while keeping the densities of the matter and the radiation species the same. In this case, the effective DE fluid compensates for the difference between the expected and the actual gravitational contributions of matter and radiation. The fluid must, therefore, take on an EoS that tracks the EoS of the dominant species at any given time.

The fact that the effective $w_{\rm DE}$ can cross $-1$ in non-minimally coupled models is well known Das et al. (2006). To see this, consider the form of the effective $w_{\rm DE}$ in the simple case of constant $\Lambda$ and $\Omega$,

[TABLE]

where $P_{m}$ and $\rho_{m}$ are the matter pressure and energy densities. It is clear from the new term in the denominator that the coupling to matter can act to decrease the energy density of the effective dark energy fluid, pushing $w_{\rm DE}$ below $-1$. While the expression for $w_{\rm DE}$ in the case of varying $\Lambda$ and $\Omega$ is more complicated, the physics is essentially the same.

III Methods

This section explains how we build the ensemble of viable Horndeski models, and the way in which we store and present the priors on $w_{\rm DE}(z)$.

III.1 Sampling the ensemble of viable scalar field dark energy models

A particular model realization is defined by specifying the EFT functions that appear in action (1). Ideally, one would adopt the most general functional forms possible for these functions, to avoid unduly focusing on only a few – potentially unrepresentative – corners of the space of possible models. In practise, this means choosing a parametrization that is broad enough to represent many different functions of time, while having sufficiently few parameters that it can be efficiently sampled using a Monte Carlo algorithm. To strike a balance between generality and tractability, we will use a number of different series expansions in the scale factor, $a$.

The first such expansion is a truncated Taylor polynomial, defined by

[TABLE]

where $N$ is the order at which we choose to truncate the expansion, $a_{0}$ is the point around which we expand, and $\{\alpha_{n}\}$ is a set of coefficients to be drawn from a random distribution. Given $\alpha_{n}$ coefficients with identical prior distributions, such a parametrization would favour the lower order terms.

The second expansion that we consider is a polynomial expansion,

[TABLE]

This differs from the Taylor expansion in the absence of the $n!$ term that suppresses higher order terms. Depending on how the $\alpha_{n}$ coefficients are drawn, this will allow rapid variations in $f(a)$ to arise more easily.

Finally, we will also use a Padé expansion,

[TABLE]

where the truncation order is now set by both $N$ and $M$. This expansion is well suited for describing models in which the functions transition from one value at small $a$ to another at large $a$.

We consider these three expansions around $a_{0}=0$ and $a_{0}=1$, to represent models that exhibit “thawing” and “freezing” behaviors respectively. We also explore several truncation orders, ranging from zero-order expansions, where all the three expansions reduce to constants, to ninth-order expansions that should allow very general model behaviors to be captured. We progressively raise the truncation order, obtaining results on the DE EoS prior for each order in turn. We find that the results stabilize beyond the fifth order, changing very little as the expansion order is further increased.

The prior distributions for all the coefficients $\{\alpha_{n}\}$ and $\{\beta_{m}\}$, that define particular realizations of the EFT functions, are chosen to be uniform in the range $[-1,1]$. We checked that extending the ranges of the flat priors does not affect the results.

We have also adopted weak priors on cosmological parameters. The relative density of radiation is fixed by the CMB temperature, while the matter density today (the sum of the baryon and the CDM contributions) was drawn from the range $\Omega_{m}\in[0,1]$. We assume the universe to be flat and the sum of neutrino masses to be $\sum_{\nu}m_{\nu}=0.06$ eV. The present-day DE density (a free parameter that is not fixed by flatness in the models we considered), was allowed to span $\Omega_{\rm DE}\in[0,1]$. The Hubble constant was allowed to vary in the range $H_{0}\in[20,100]\,\mbox{km/s/Mpc}$. In addition, we impose a weak prior on the present value of the gravitational coupling, allowing variations of no more than $10\%$, and the present value of the speed of gravitational waves, allowing variations of no more than $20\%$. The latter two conditions help to exclude models that are in obvious contradiction with laboratory experiments and observations of the nearby universe Perenon et al. (2015); Caves (1980); Moore and Nelson (2001); Beltran Jimenez et al. (2016).

When deriving the probability distributions of $w_{\rm DE}$ values at different $z$, we also take into account mild background data constraints to remove histories that are obviously ruled out. Namely, we use information from Baryon Acoustic Oscillation (BAO) measurements Beutler et al. (2011); Percival et al. (2010); Padmanabhan et al. (2012); Anderson et al. (2013), estimates of the Hubble constant Riess et al. (2011), and supernovae distance measures Suzuki et al. (2012), but with a significantly enlarged covariance (by a factor of four) to avoid the biasing of our results by tensions between these datasets. Let us stress that we do not use this information when deriving the covariance of the $w_{\rm DE}$ bins, so that it can be used in reconstructions of $w_{\rm DE}(z)$ as a purely theoretical prior.

We explore the space of cosmological parameters and $\{\alpha_{n}\}$ coefficients using a Monte Carlo sampling procedure. For each sample, after solving the background equations (with initial conditions defined at the present day; see Appendix A), we check the stability of the corresponding model and, if accepted, compute $w_{\rm DE}(z)$. To ensure good coverage, we enforce a minimum number of $10^{5}$ accepted samples. Depending on the acceptance rate, this results in $\sim 10^{6}-10^{8}$ total samples.

While we have chosen several different ways to parametrize the behavior of the EFT functions, we can marginalize over these choices if desired. Given a set of models $\mathcal{M}$ corresponding to the different available options, we can compute

[TABLE]

where $P\left(\mathcal{M}\right)$ is the prior probability of each of the choices. Note that we are not interested in the normalization of $P\left(\vec{w}\right)$, which can be recovered a posteriori, so there is no need to impose the constraint $\sum_{\mathcal{M}}P(\mathcal{M})=1$. We will assume uninformative priors on the choice of parametrization, so that $P(\mathcal{M})$ is equal for all models. Since we have generated Monte Carlo samples of $P\left(\vec{w}|\mathcal{M}\right)$ with the same number of accepted points, we can then obtain $P(\vec{w})$ simply by merging the Monte Carlo samples from all parametrizations and re-weighting them with their respective acceptance rates.

III.2 The priors on $w_{\rm DE}(z)$

Given an expansion history from a particular model realization, we store the corresponding $w_{\rm DE}(z)$ in $100$ linearly-spaced bins in scale factor that cover the redshift range $z\in[0,6]$. From the Monte Carlo samples of the binned EoS, $w_{\rm DE}(z_{i})=w_{i}$, we compute the mean and the covariance,

[TABLE]

One can also define the normalized correlation matrix as

[TABLE]

Under the assumption of Gaussianity, the prior covariance matrix can be turned into a prior probability distribution and used in reconstructions of $w_{\rm DE}(z)$ from data, following the procedure suggested in Crittenden et al. (2012).

While having a numerically-obtained discrete prior covariance matrix $C_{ij}$ may be sufficient for many practical applications, it is useful to also have an analytical expression characterizing the correlation of $w_{\rm DE}(z)$ between an arbitrary pair of redshifts. To this end, following Crittenden et al. (2009, 2012), we can introduce a two-point correlation function $C(a,a^{\prime})$ defined as

[TABLE]

As in the discrete case, Eq. (11), this prior covariance can be expressed in terms of the normalized correlation, $\mathcal{C}(a,a^{\prime})$, equal to unity for $a=a^{\prime}$, and the auto-correlation (or variance), $C(a)\equiv C(a,a)$, i.e.

[TABLE]

In Sec. IV, we will obtain analytical forms of $C(a,a^{\prime})$ for different classes of models by separately fitting $\mathcal{C}(a_{i},a_{j})$ and $C(a_{i})$ to the numerically-obtained $C_{ij}$ and $C_{ii}$.

As derived quantities, in Sec. IV, we also compute projections of the binned $w_{\rm DE}(z)$ model onto parameters of the CPL parametrization,

[TABLE]

where $w_{a}$ is obtained from

[TABLE]

This, amongst other things, allows us to compare our results to those in Marsh et al. (2014), where priors on ($w_{0}$,$w_{a}$) were derived for quintessence by sampling a broad range of scalar field potentials.

IV Results

In this section, we present the results of our Monte Carlo exploration of allowed DE EoS histories in scalar-tensor theories. We first discuss the PDFs of $w_{\rm DE}$ at different $z$, and then the covariance of $w_{\rm DE}(z)$ and the corresponding analytical fitting formulae.

IV.1 The ensemble of $w_{\rm DE}(z)$ functions

By following the procedure for sampling the ensemble of scalar field theories and calculating the corresponding expansion histories described earlier, we obtain an ensemble of functions $w_{\rm DE}(z)$. Fig. 2 (left column) shows the probability density of $w_{\rm DE}$ as a function of $z$, derived from the Monte Carlo sample of two classes of models: the minimally coupled quintessence field, and the full Horndeski theory. These distributions have been marginalized over cosmological parameters, boundary conditions, EFT function parametrizations, orders of the Taylor expansion, and freezing/thawing behaviors (values of $a_{0}$). In addition, the mild observational constraints discussed in the previous section have been applied. The figures are shaded according to probability density, and show the mean $w_{\rm DE}(z)$, as well as the 68%, 95%, and 99% C.L. contours.

Despite the large amount of freedom in the specification of the models, there is a remarkably well-defined structure to the $w_{\rm DE}(z)$ prior probability density. The overall behavior is similar between the two model classes, despite the significant increase in model freedom in going from the quintessence ($\Lambda$-only) class to the most general Horndeski class. In both cases, the EoS transitions from $w_{\rm DE}\approx 0$ deep in the matter-dominated regime to $\sim-1$ at low redshift. The high-redshift behavior is quite generic, for the reasons given in Sect. II.1: only quintessence models with kinetic energy fine-tuned to be close to zero can remain potential-dominated (and thus have $w\ll 0$) at early times, while the tracking behavior is common in non-minimally coupled models, giving $w_{\rm DE}\approx w_{m}=0$ in the matter-dominated regime. The low-redshift behavior is typical of freezing models, where the kinetic energy of the scalar field always remains sub-dominant, and is significantly driven by the mild data constraints. These down-weight all of the $w_{\rm DE}(z)$ functions that result in non-accelerating cosmologies with $w_{\rm DE}>-1/3$ at $z=0$.

There are, however, some significant qualitative differences between the two model classes. As shown in panel 1a) of Fig. 2, all models in the quintessence class respect $w_{\rm DE}\geq-1$ at all times, as expected. This bound is not respected in non-minimally coupled models, which have a significant probability density of models with $w_{\rm DE}<-1$ at low redshift. These “phantom” models are relatively common, with $40\%$ of viable Horndeski models having $w_{\rm DE}<-1$ at $z=0$.

Another feature worthy of note is the lack of models with $w_{\rm DE}\simeq-1={\rm const.}$, i.e. a cosmological constant-like evolution at all redshifts. While values of $w_{\rm DE}\approx-1$ are very common at low redshift, they are rather unlikely in both model classes at $z\simeq 2$. This can be explained dynamically, by the tracking and kinetic energy tuning arguments above, but has interesting observational implications. While the fractional DE density is low at $z\gtrsim 2$, making measurements of the EoS more difficult, our results suggest that EoS constraints in this regime may actually be the most powerful for distinguishing a cosmological constant from dynamical DE.

The right-hand column of Fig. 2 shows the projection of the $w_{\rm DE}(z)$ curves onto the $w_{0}$ and $w_{a}$ parameters. The results for the quintessence class are shown in panel 1b). This should be compared with Fig. 1 of Marsh et al. (2014), which also shows the ($w_{0},w_{a}$) distribution for a Monte Carlo-generated sample of quintessence models. The orientation of the prior is roughly the same, with most models following a track that goes from $(w_{0},w_{a})=(-1,0)$ towards more positive $w_{0}$ and more negative $w_{a}$. This is in part due to the $w_{\rm DE}\geq-1$ restriction that applies to quintessence models, which excludes a substantial region of the $(w_{0},w_{a})$ plane.

The distribution shown in the panel 1b) of Fig. 2 is quite broad, with a significant fraction of models having $w_{a}>0$. In contrast, the models of Marsh et al. (2014) occupied a much narrower track that is almost entirely at $w_{a}<0$. Part of the reason for this is that the Monte Carlo method in Marsh et al. (2014) used randomly-generated scalar field potentials to define the quintessence models, while we are working directly in terms of randomly-generated EFT functions. This effectively leads to a different weighting of the models; simple random distributions over the parametrized potentials may translate to non-trivial distributions over the EFT functions. In this case, the weighting obtained by parametrizing the theories at the EFT function level seems to make it easier to obtain thawing models (c.f. Huterer and Peiris (2007), who also find relatively more models with $w_{a}>0$). Also, the models in Marsh et al. (2014) were selected to be “physically-motivated,” in the sense that the forms of their potentials were inspired by high-energy theory, like moduli fields and axions. These models form a more restricted class than the general, phenomenological EFT approach considered here. In particular, treating $\Lambda(a)$, which represents a particular evolution history for the quintessence Lagrangian, as being a priori arbitrary, amounts to a significantly more agnostic approach, leading to a broader allowed range of ($w_{0},w_{a}$) values.

The $(w_{0},w_{a})$ distribution for the Horndeski class is broader and much less constrained, primarily as a result of the broadening of the $w_{\rm DE}(z)$ distributions, and the ability of the models to extend to $w_{\rm DE}<-1$. The bulk of the models remain within the vicinity of the $(-1,0)$ point in both classes, although there is a very slight shift towards less negative $w_{0}$ and positive $w_{a}$ for the full Horndeski class. This is despite the fact that very few of the models have equations of state that resemble a cosmological constant out to high redshift. In fact, the mean functional form in the left-hand panels approximately satisfies $w|_{z=0}\to-1$ and $w|_{z\to\infty}\to 0$ in both model classes, which would imply $w_{0}=-1$ and $w_{a}=+1$. These values are within the 68% region for all three classes, but shifted slightly from the mean (shown by the white points in Fig. 2). This is, at least in part, because both $w_{0}$ and $w_{a}$ have been evaluated directly at $z=0$. It is clear from Fig. 2, however, that the derivative of $w_{\rm DE}(z)$ is not constant with redshift, i.e. the $(w_{0},w_{a})$ parametrization cannot fully reproduce its true redshift evolution when defined in this way.

In addition to considering the quintessence sub-class of Horndeski theories, we have separately studied the case of the GBD models, i.e. models with a canonical scalar field kinetic term. Within the EFT framework, they are described by $\Lambda(a)\neq{\rm const.}$ and $\Omega(a)\neq 0$, with $c(a)$ derived from $\{\Lambda,\Omega\}$, and the other functions set to zero. We find the results for the PDF of $w_{\rm DE}(z)$ in this case are practically the same as in the full Horndeski case. On one hand, this is expected, since only $\Lambda(a)$ and $\Omega(a)$ affect the background evolution. On other hand, the higher order functions figure in the stability conditions, which do affect the acceptance rate in the Monte Carlo sampling. We find, therefore, that the stabilization effect of higher order functions has a limited impact on the shape of the prior.

In the analysis presented so far, we have applied mild observational constraints to ensure that the models are broadly consistent with the properties of the real Universe (see Sect. III). The results are qualitatively similar even if these constraints are relaxed, with the only significant difference being the presence of trajectories that remain above $w_{\rm DE}=-1/3$ at low redshift.

IV.2 Theoretical priors on correlations of $w_{\rm DE}(z)$

In this section, we present the covariance between values of $w_{\rm DE}$ at different redshifts, obtained from averaging over the ensemble of $w_{\rm DE}(z)$ functions from each of the model classes. In doing so, we do not employ any observational constraints, as even relatively mild constraints result in correlations induced by the data overwhelming those coming from theory, effectively hiding the information that we want to extract.

We separate the covariance into the normalized correlation matrix, $\mathcal{C}(a_{i},a_{j})$, and the auto-correlation (or variance), $C(a_{i})$, as in Eq. (12). The top left panels in Figs. 3 and 4 show the numerically obtained $\mathcal{C}(a_{i},a_{j})$ at binned values of $a$ for the quintessence and the general Horndeski classes, respectively. The points in Fig. 5 show the variance $C(a_{i})$ for the two models.

As mentioned in Sec. III.2, having analytical expression for the correlation prior can simplify practical applications of our priors to reconstructing $w_{\rm DE}(z)$ from data using methods similar to those in Zhao et al. (2012, 2017). We can obtain them by fitting simple functional forms to the numerically obtained discrete covariance matrices.

First, we perform a least-squares fit to the auto-correlation $C(a_{i})$ using the following functional forms:

• •

a Taylor expansion, $C(x)=\alpha+\beta\,(x-x_{0})$;

• •

an exponential, $C(x)=\alpha+\beta\,\exp[\gamma\,(x-x_{0})]$;

• •

a power law, $C(x)=\alpha+\beta\,(x-x_{0})^{\gamma}$;

where $x$ is the evolution variable for which we try three different choices: the scale factor, $x=a$; the redshift, $x=z$; and the number of e-folds, $x=\ln(a)$. We also try two choices of the reference point, $x_{0}$, in each case: for $x=a$ and $x=\ln(a)$, we either fix $a_{0}=1$ or allow $a_{0}$ to be free, and for $x=z$ we either fix $z_{0}=0$ or allow $z_{0}$ to be free. After fitting all of these different parametrizations to the measured auto-correlation function for each model class, we select the ones that result in the smallest residuals in each case. The choices that give the best fits are summarized in Table 1, and plotted in Fig. 5 as a function of the scale factor. One can see that the exponential fitting formula with $x=\ln a$ and $a_{0}=1$ works very well in the Horndeski case, while the exponential form with $x=a$ and $a_{0}=1$ works best in the quintessence case. The analysis for the GBD class of models gives results that are almost indistinguishable from the Horndeski case.

As expected, the auto-correlation (variance) is generally bigger in the Horndeski case, as these models permit more freedom than quintessence. In both classes, the auto-correlation is small at high redshifts, where the field is expected to have little or no dynamics, and so tends to result in similar values of $w_{\rm DE}$ in all models. The spread is much more significant at low redshifts however, where a greater range of dynamical behaviors is realized. The auto-correlation of the quintessence models is clipped by the hard bound of $w_{\rm DE}\geq-1$, reducing the allowed variation from model to model.

Next, we use a similar procedure to find the analytical form of $\mathcal{C}(a,a^{\prime})$. We consider three fitting forms. One is a generalized version of the CPZ parametrization Crittenden et al. (2009):

[TABLE]

where $x$ and $y$ are either the scale factor, the redshift or $\ln a$, $\xi$ is a parameter defining the correlation time scale, and $n$ is a free parameter. We also separately fit to (16) with $n=2$, which is the CPZ form. In addition, we consider an exponentially decaying correlation,

[TABLE]

where the physical meaning of the parameters is the same as in the previous parametrization (16). We then evaluate the least-squares distance between $\mathcal{C}(a_{i},a_{j})$ and the numerically found $\mathcal{C}_{ij}$, and find the parameters that minimize it, as well as the fitting form that results in the smallest residuals for each class of models.

We present the best fit analytical forms of $\mathcal{C}(a,a^{\prime})$ and the corresponding residuals in Table 1, along with best fit parameters of the CPZ parametrization. The two top right panels of Figs. 3 and 4 show the best fit CPZ and overall best fit $\mathcal{C}(a,a^{\prime})$ for the quintessence and Horndeski models, respectively. To have a more detailed visual check of the goodness of fit, in the bottom panel of Figs. 3 and 4, we also plot $\mathcal{C}(a,a^{\prime})$ as a function of $a^{\prime}$ for several fixed values of $a$. One can see that the exponentially decaying parametrization works marginally better than the CPZ form, as it has more freedom built in. Overall, the analytical expressions work well in reproducing the correlation in different models. We also note that the best fit correlation length, $\xi$, is roughly the same in both parametrizations for each model.

One can clearly see that the correlation $\mathcal{C}(a,a^{\prime})$ is long ranged for quintessence at all $a$, and is slowly decaying with $|a-a^{\prime}|$. In Horndeski models, on the other hand, the correlation is generally shorter ranged, and decays faster at early times compared to late times. This is expected, since Horndeski models have more freedom in the choice of expansion histories compared to quintessence, hence values of $w_{\rm DE}$ at different times are less correlated.

It is also interesting to note that the auto-correlation, $C(a)$, depends linearly on $a$ for quintessence, while it scales as $a^{2}$ in the Horndeski case. Also, the correlation $\mathcal{C}(a,a^{\prime})$ scales as a power of $|a-a^{\prime}|$ for quintessence, while in the case of Horndeski it scales as a power of $|\ln a-\ln a^{\prime}|$. This is due to the fact that, in quintessence, $w_{\rm DE}$ is solely a property of DE, which is only minimally coupled to other fluids. In non-minimally coupled models, on the other hand, $w_{\rm DE}$ is an effective quantity, determined largely by how the other fluids scale with redshift. This is true especially at higher redshifts, where the variance and the correlation of $w_{\rm DE}(z)$ in Horndeski models is set by the dynamics of the matter fluid, which has a fixed time-dependence and uniformly distributed random amplitude set by $\Omega_{M}$.

IV.3 A note on practical applications of the prior

As mentioned earlier, one can use the prior covariance obtained above to aid reconstruction of $w_{\rm DE}(z)$ from data using the method developed in Crittenden et al. (2012) and applied in Zhao et al. (2012, 2017). There, a piece-wise (binned) $w_{\rm DE}(z_{i})=w_{i}$ is fit to data along with other model parameters, $\vec{p}$, using the usual MCMC method for sampling the posterior PDF, ${\cal P}(\{w_{i}\},\vec{p})$. According to Bayes’ theorem, the posterior PDF is the product of the likelihood and the prior PDFs:

[TABLE]

where the prior on $\{w_{i}\}$ is obtained, under the assumption of Gaussianity, from the covariance matrix $C_{ij}$ discussed in the previous subsection:

[TABLE]

This prior probability explicitly depends on the average EoS, $\{\bar{w}_{i}\}$. We note that this is not the mean $w_{\rm DE}(z)$ shown as white lines in panels 1(a) and 2(a) of Fig. 2. The PDFs obtained in that figure used mild data constraints that pushed the mean closer to the observationally favoured $w_{\rm DE}\sim-1$ region.

In practice, any reasonable choice for the fiducial ${\bar{w}}_{i}$, e.g. any choice within the 1$\sigma$ band of the EoS shown in Fig. 2, should be more or less equally acceptable. The primary purpose of the prior is to add curvature to the PDF along otherwise flat directions in the parameter space, thus eliminating the degeneracies between $w_{i}$ and helping the MCMC to converge. A reconstruction that is highly sensitive to the choice of ${\bar{w}}_{i}$ should not be trusted.

Other options, discussed in Crittenden et al. (2012), include using the so-called “running average”, or marginalizing over ${\bar{w}}_{i}$ altogether. We refer the reader to Crittenden et al. (2012); Zhao et al. (2012, 2017) for further details on the method and its application.

V Summary

We have derived theoretical priors on the effective DE EoS within the Horndeski class of scalar-tensor theories, which includes all models with a single scalar field that have second order equations of motion. We separately considered the widely studied sub-classes of the minimally coupled scalar field, or quintessence Ratra and Peebles (1988); Caldwell et al. (1998), and models with the canonical form of the scalar field kinetic energy, i.e. the generalized Brans-Dicke (GBD) models Brans and Dicke (1961); Carroll (2004). Overall, we find that the covariance of $w_{\rm DE}(z)$ in GBD is indistinguishable from that of the general Horndeski case.

Our priors on $w_{\rm DE}(z)$ are stored in the form of a covariance matrix for binned $w_{\rm DE}$, which can be projected onto priors on the parameters of any specific parametrization, such as CPL.

We found that there are notable differences between the case of the minimally coupled quintessence and the non-minimally coupled models, both in the mean values and the covariance of binned $w_{\rm DE}$. We found simple analytical forms for the correlation function, describing the correlation of $w_{\rm DE}$ at different redshifts, that fit our numerical results well. These should simplify the practical application of the priors to reconstructions of $w_{\rm DE}(z)$ from data. Such a reconstruction will be the subject of future work.

Acknowledgements.

We thank Wayne Hu and Gong-Bo Zhao for useful comments. MR is supported by U.S. Dept. of Energy contract DE-FG02-13ER41958. PB is supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, California Institute of Technology, administered by Universities Space Research Association under contract with NASA. AS acknowledges support from the NWO and the Dutch Ministry of Education, Culture and Science (OCW), and also from the D-ITP consortium, a program of the NWO that is funded by the OCW. The work of LP is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).

Appendix A The background solution in EFT

Here, we discuss in detail the procedure for solving for the cosmological background evolution given a choice of EFT functions. In what follows, the over-dot represents a derivative with respect to the conformal time, the accent mark represents a derivative with respect to the scale factor, and the subscript $m$ indicates the sum over all particle species: CDM, baryons, photons, and massless and massive neutrinos. The Friedmann equations for EFT are given by:

[TABLE]

Changing the time coordinate from conformal time to the scale factor, Eq. (20) can be recast as

[TABLE]

Combining this with Eq. (21), and introducing $y\equiv{\mathcal{H}}^{2}$, we obtain

[TABLE]

which is the differential equation we solve to find ${\mathcal{H}}(a)$. The initial conditions for this equation, along with the current values of the model parameters, can be derived from the Friedmann constraint (20) computed today,

[TABLE]

where we have defined

[TABLE]

A few comments regarding this equation are in order. The definition of matter densities follows from what we can directly measure, at least in principle. The quantities $\rho_{m}^{0}$ are in fact the densities that can be measured with a non-gravitational experiment by an observer that is assuming that all these species are moving on geodesics of the metric. These quantities are combined in the equation (24) that is usually used as the flatness constraint. Normally, when assuming flatness, we use this to compute the value of $\Omega_{\Lambda}^{0}$ given the value of $\Omega_{m}$. Here, the situation is slightly different, as both $\Omega_{\Lambda}^{0}$ and $\Omega_{m}^{0}$ can be chosen arbitrarily, and constitute two of the theory parameters that we are going to sample. On the other hand, the flatness constraint is satisfied, once the present day value of the gravitational constant is fixed, by a suitable choice of the present day value of $c$.

Once we have a solution of Eq. (2), we can deduce the effective DE pressure and density from the standard form of the Friedmann equations:

[TABLE]

and compute the effective DE EoS as:

[TABLE]