Generalized Brans-Dicke theories in light of evolving dark energy

TL;DR

This paper reconstructs generalized Brans-Dicke scalar-tensor theories from recent dark energy data, analyzing their effects on cosmological observables and showing that certain signatures can rule out these theories as explanations for dark energy dynamics.

Contribution

It introduces a method to reconstruct GBD Lagrangians from data and investigates their observable signatures, providing a way to test scalar-tensor theories against cosmological observations.

Findings

Reconstructed GBD models can mimic observed dark energy behavior.

The ISW effect can distinguish GBD theories from other dark energy models.

Certain GBD models are ruled out by CMB and galaxy cross-correlation data.

Abstract

The expansion history of the Universe reconstructed from a combination of recent data indicates a preference for a changing Dark Energy (DE) density. Moreover, the DE density appears to be increasing with cosmic time, with its equation of state being below -1 on average, and possibly crossing the so-called phantom divide. Scalar-tensor theories, in which the scalar field mediates a force between matter particles, offer a natural framework in which the effective DE equation of state can be less than -1 and cross the phantom barrier. We consider the generalized Brans-Dicke (GBD) class of scalar-tensor theories and reconstruct their Lagrangian given the effective DE density extracted from recent data. Then, given the reconstructed Lagrangian, we solve for the linear perturbations and investigate the characteristic signatures of these reconstructed GBD in the cosmological observables, such as the cosmic microwave background (CMB) anisotropy, the galaxy number counts, and their cross-correlations. In particular, we demonstrate that the Integrated Sachs-Wolfe (ISW) effect probed by the cross-correlation of CMB with the matter distribution can rule out scalar-tensor theories as the explanation of the observed DE dynamics independently from the laboratory and solar system fifth force constraints.

Equations

$$H^2=H_0^2\left[E_m(a)+{\Omega }_{\Lambda }X(a)\right].$$

$$\xi (|\Delta a|)=\frac{\xi (0)}{1+(|\Delta a|/a_c{)}^2}.$$

$$X_{\mathrm{fit}}(a)=A\tanh [B(a-C)]+D,$$

$$\begin{array}{rl}S& =\int d^{4}x\sqrt{-g}\left[\frac{m_0^2}{2}F(\varphi )R-\frac{1}{2}K(\varphi )(\partial \varphi {)}^2-U(\varphi )\right]\\ & +S_m[g_{\mu \nu },{\chi }_i],\end{array}$$

$$F{G}_{\mu \nu }=\frac{1}{m_0^2}\left(T_{\mu \nu }^{\mathrm{m}}+T_{\mu \nu }^{\varphi }\right)+{\nabla }_{\mu }{\nabla }_{\nu }F-g_{\mu \nu }\square F,$$

$$T_{\mu \nu }^{\varphi }\equiv {\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -g_{\mu \nu }\left[\frac{1}{2}{\partial }_{\alpha }\varphi {\partial }^{\alpha }\varphi +U(\varphi )\right].$$

$$\square \varphi =U_{\varphi }-\frac{m_0^2}{2}{F}_{\varphi }R,$$

$$G_{\mu \nu }$$

$$(T_{\mathrm{DE}}^{\mathrm{eff}}{)}_{\mu \nu }\equiv F^{-1}\left\{T_{\mu \nu }^{\varphi }+{\nabla }_{\mu }{\nabla }_{\nu }F-g_{\mu \nu }\square F+(1-F)T_{\mu \nu }^m\right\}.$$

$${\rho }_{\mathrm{DE}}^{\mathrm{eff}}=F^{-1}\left\{{\dot{\varphi }}^2/(2a^2)+U(\varphi )-3\mathcal{H}\dot{F}/a^2+(1-F){\rho }_m\right\},$$

$$p_{\mathrm{DE}}^{\mathrm{eff}}=F^{-1}\left\{{\dot{\varphi }}^2/(2a^2)-U(\varphi )+\mathcal{H}\dot{F}/a^2+\ddot{F}/a^2\right\}.$$

$${\mathcal{H}}^2={\left(\frac{\dot{a}}{a}\right)}^2=\frac{a^2}{3m_0^2}[{\rho }_m(a)+{\rho }_{\mathrm{DE}}^{\mathrm{eff}}(a)],$$

$$w_{\mathrm{DE}}^{\mathrm{eff}}=\frac{{\dot{\varphi }}^2/(2a^2)-U(\varphi )+\mathcal{H}\dot{F}/a^2+\ddot{F}/a^2}{{\dot{\varphi }}^2/(2a^2)+U(\varphi )-3\mathcal{H}\dot{F}/a^2+(1-F){\rho }_M},$$

$$F(\varphi )=\exp (\xi \varphi ),$$

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

$$\begin{array}{rl}\mathcal{G}\frac{\ddot{a}}{a}& =2\mathcal{D}{\mathcal{H}}^2-\frac{a^2{\mathcal{H}}^2}{2}\left(\frac{1}{3}+F_{\varphi \varphi }\right)({\varphi }^{\prime }{)}^2\\ & -\frac{1}{2}\frac{(\rho +P)a^2}{m_0^2}-\frac{1}{2}{F}_{\varphi }{\mathcal{H}}^2\left({\varphi }^{\prime \prime }-{\varphi }^{\prime }\right),\end{array}$$

$$\mathcal{D}$$

$$\mathcal{G}$$

$$\begin{array}{rl}{\varphi }^{\prime \prime }& =-\frac{1+F_{\varphi \varphi }}{F_{\varphi }}({\varphi }^{\prime }{)}^2\\ & +\left(1+\frac{1}{2}\frac{3E_m+4E_r-E_{\nu }^{\prime }-{\Omega }_{\Lambda }{X}^{\prime }}{E_m+E_r+E_{\nu }+{\Omega }_{\Lambda }X}\right){\varphi }^{\prime }\\ & +\frac{1}{F_{\varphi }}\frac{(F-1)(3E_m+4E_r-E_{\nu }^{\prime })-F{\Omega }_{\Lambda }{X}^{\prime }}{E_m+E_r+E_{\nu }+{\Omega }_{\Lambda }X},\end{array}$$

$$\begin{array}{rl}U{a}^2& =3\mathcal{D}{H}_0^2{a}^2(E_m+E_r+E_{\nu }+{\Omega }_{\Lambda }X)\\ & -3H_0^2{a}^2(E_m+E_r+E_{\nu }).\end{array}$$

$$\frac{{\varphi }_{\mathrm{ini}}^{\prime }}{m_0}\in [-10^{-6},10^{-6}],{\log }_{10}{a}_{\mathrm{ini}}\in [-3,-1],\xi \in [0.1,10],$$

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

$$\displaystyle\begin{split}\dot{\mathcal{H}}&=\frac{1}{F+\frac{1}{2}aF^{\prime}}\biggl{\{}\left[F+2aF^{\prime}+a^{2}F^{\prime\prime}\right]\mathcal{H}^{2}\\ &-\frac{1}{2m_{0}^{2}}\left[Pa^{2}+\frac{1}{2}\dot{\phi}^{2}-Ua^{2}\right]\biggr{\}},\end{split}$$

$$\begin{array}{rl}\dot{\mathcal{H}}& =\left\{\left[\frac{5}{2}F+\frac{7}{2}a{F}^{\prime }+a^2{F}^{\prime \prime }\right]{\mathcal{H}}^2-\frac{(\rho +P)a^2}{2m_0^2}-\frac{1}{2m_0^2}{\dot{\varphi }}^2\right\}\\ & \times {\left(F+\frac{1}{2}a{F}^{\prime }\right)}^{-1}.\end{array}$$

$$\begin{array}{rl}\frac{U{a}^2}{m_0^2}& ={\mathcal{H}}^2\left[\frac{1}{2}F-\frac{1}{2}a{F}^{\prime }-a^2{F}^{\prime \prime }\right]\\ & +\frac{(P-\rho )a^2}{2m_0^2}+\dot{\mathcal{H}}\left[F+\frac{1}{2}a{F}^{\prime }\right].\end{array}$$

$$F(a)$$

$${\alpha }_i\in [-1,1].$$

$$S=\int d^{4}x\sqrt{-g}\left\{\frac{m_0^2}{2}\left[1+\Omega (\tau )\right]R+\Lambda (\tau )+c(\tau )a^2\delta g^{00}\right\}$$

$$\Omega (a)$$

$$\frac{c{a}^2}{m_0^2}$$

$$\frac{\Lambda a^2}{m_0^2}$$

$$C_{\ell }^{Tg}=\frac{2}{\pi }\int dk{k}^2{\Delta }_{\ell }^{\mathrm{ISW}}(k,{\tau }_0){\Delta }_{\ell }^{\mathrm{GNC}}(k,{\tau }_0){\mathcal{P}}_{\mathcal{R}}(k),$$

$${\Delta }_{\ell }^{\mathrm{ISW}}(k,{\tau }_0)=-{\int }_{{\tau }^{\ast }}^{{\tau }_0}d\tau (\dot{\Phi }+\dot{\Psi })j_{\ell }[k({\tau }_0-\tau )],$$

Full text

I Introduction

The observed accelerated expansion of the Universe has been puzzling cosmologists since its discovery two decades ago (Perlmutter et al., 1999; Riess et al., 1998). Within the context of General Relativity (GR), it implies the existence of an energy-momentum component with a negative equation of state (EOS), referred to as Dark Energy (DE). The standard cosmological model, $\Lambda$CDM, in which DE is the constant energy of the vacuum, provides a good fit to a plethora of cosmological observations such as the cosmic microwave background (CMB) anisotropies (Bennett et al., 2013; Ade et al., 2016a), baryon acoustic oscillations (BAO) (Percival et al., 2002; Beutler et al., 2011; Alam et al., 2017), type Ia supernovae (Conley et al., 2011; Suzuki et al., 2012), galaxy clustering (Zehavi et al., 2011) and galaxy lensing (Heymans et al., 2012; Hildebrandt et al., 2017). However, $\Lambda$CDM is not fully satisfactory from the theoretical perspective, as the observed value of the vacuum energy requires an extreme fine tuning of the cosmological constant $\Lambda$ in the context of the present understanding of particle interactions (Burgess, 2017). Also, with the data becoming more accurate, several “tensions” between different datasets have begun to arise when interpreting observations within the $\Lambda$CDM model Ade et al. (2016a); Riess et al. (2016, 2019); Abbott et al. (2018); Wong et al. (2019). Although these tensions might just be due to unaccounted systematic effects or rare statistical fluctuations (Scott, 2018), they generated significant interest in possible extensions of $\Lambda$CDM capable of relieving the tensions Dvorkin et al. (2014); Di Valentino et al. (2016, 2018a, 2018b); Adhikari and Huterer (2019); Lin et al. (2019); Peirone et al. (2019), including the possibility of the DE density evolving with time Zhao et al. (2017a); Wang et al. (2018); Yang et al. (2018); Di Valentino et al. (2019).

Using a combination of available observations, non-parametric reconstructions of the DE dynamics were performed in (Zhao et al., 2017a; Wang et al., 2018). Interestingly, they show a preference for an increasing effective DE density, i.e. one with an EOS, $\smash{w_{\rm DE}^{\rm eff}<-1}$. The reconstruction also shows crossing of the so-called phantom divide Caldwell (2002); Carroll et al. (2003); Vikman (2005) of $\smash{w_{\rm DE}^{\rm eff}=-1}$. Such dynamics cannot be explained by a minimally coupled quintessence field DE, but could be realized in scalar-tensor extensions of GR where the additional scalar field $\phi$ mediates a force between particles Carroll et al. (2003); Das et al. (2006); Amendola and Tsujikawa (2008). In fact, scalar-tensor theories possess enough freedom to reproduce any expansion history.

The aim of this paper is to investigate scalar-tensor theories of the generalized Brans-Dicke (GBD) type capable of realizing the expansion histories reconstructed in (Zhao et al., 2017a; Wang et al., 2018). Using the observed $H(a)$ as input, and making certain assumptions about the scalar field coupling function, we systematically scan the parameter space to reconstruct the GBD Lagrangians consistent with that $H(a)$. We then solve for the cosmological perturbations and calculate predictions for the CMB and galaxy power spectra and other observables to isolate the theories that are in agreement with current data.

Late time deviations from $\Lambda$CDM are mainly encoded in the CMB temperature through the integrated Sachs-Wolfe (ISW) effect. Although too small to be detected from the CMB temperature auto-correlation, the ISW contribution can be probed by cross correlating the CMB temperature maps with the foreground galaxies number counts (Crittenden and Turok, 1996; Afshordi, 2004; Giannantonio et al., 2008), which can be a useful probe for DE Pogosian (2005); Pogosian et al. (2005). In $\Lambda$CDM, the accelerating expansion results in decaying gravitational potentials, yielding a strictly positive ISW effect. In scalar-tensor theories however, the ISW effect can have a positive or negative sign depending on whether the enhanced clustering due to the fifth force, which yields a negative ISW, dominates over the effect of the accelerating expansion Song et al. (2007); Renk et al. (2016). We find that most of the GBD theories reconstructed in this work predict CMB-matter cross-correlations that are significantly different from those in $\Lambda$CDM and, therefore, can be ruled out or confirmed with the next generation galaxy surveys, such as DESI DESI Collaboration et al. (2016), LSST Ivezić et al. (2019); Abell et al. (2009) and Euclid Refregier (2009); Cimatti et al. (2009); Laureijs (2009).

The paper is structured as follows. In Sec. II we discuss the DE reconstruction results from (Wang et al., 2018) and introduce expansion models used in this work, along with a brief review of GBD theories. In Sec. III we introduce the GBD theories and describe two different ways in which we have reconstructed their Lagrangians from a given expansion history. Then, in Sec. IV, we evaluate the key cosmological observables predicted by the viable GBD theories that includes the CMB-galaxy cross-correlation spectra at three representative redshifts, showing how the latter can help to rule out or confirm GBD as an alternative to $\Lambda$CDM. We summarize our findings in Sec. V.

II The reconstructed dark energy density

A Bayesian, non-parametric reconstruction of the time-evolution of the DE density was performed in (Wang et al., 2018) using the correlated prior method introduced in (Crittenden et al., 2012, 2009). The effective DE energy density is modelled through the dimensionless function $X(a)$ that enters the Friedmann equation via

[TABLE]

where $\smash{E_{m}\equiv\sum_{i}\rho_{i}(a)/\rho_{\rm crit}^{0}}$ includes contributions of all matter and radiation fields, i.e. baryons, cold dark matter (CDM), photons and neutrinos, and $\smash{X(a)\equiv\rho_{\rm DE}^{\rm eff}(a)/\rho_{\rm DE}^{\rm eff}(a=1)}$ is due to any contribution to the standard Friedmann equation from terms other than the matter and radiation. Solving for the cosmological perturbations would require making additional assumptions regarding the underlying DE or modified gravity (MG) theory Silvestri and Trodden (2009); Clifton et al. (2012), hence, only observables probing the background expansion were used in (Wang et al., 2018) to keep the reconstruction model-independent. The datasets included the CMB distance priors, the “Joint Light-curve Analysis” sample of supernovae type Ia (SNe Ia) (Betoule et al., 2014), the Hubble parameter $H_{0}$ from (Riess et al., 2016), the Observational Hubble parameters Data (Moresco et al., 2016), and the BAO distance measurements from i) the 6dF Galaxy Survey (Beutler et al., 2011), ii) SDSS DR7 Main Galaxy Sample (Ross et al., 2015), iii) the tomographic BOSS DR12 (Zhao et al., 2017b; Wang et al., 2017), iv) eBOSS DR14 quasar sample (Ata et al., 2018) and the Lyman-$\alpha$ forest of BOSS DR11 quasars (Font-Ribera et al., 2014; Delubac et al., 2015).

The correlated prior method assumes that, at each scale factor $a$, $X(a)$ is a Gaussian random variable with values at different $a$ correlated according to a specified correlation function, taken to be (Crittenden et al., 2012, 2009),

[TABLE]

As demonstrated in (Crittenden et al., 2012), the details of the particular functional form are not essential as long as it has the rough shape that interpolates from 1 at $\Delta a=0$ to 0 at $\Delta a\gg a_{c}$. Here, $a_{c}$ determines the correlation length and $\xi(0)$ sets the strength of the prior and is related to the expected variance of the mean $\sigma_{\bar{X}}^{2}$ through $\sigma_{\bar{X}}^{2}\simeq\pi\xi(0)a_{c}/(a_{\rm max}-a_{\rm min})$. The Gaussian prior effectively acts as an extra term in the total $\chi^{2}$, that is used to constrain the values of $X(a)$ in $40$ bins in the interval $a\in[0.001,1]$. The advantage of this approach is that it allows one to control the strength of the prior and find the Bayesian evidence for each choice of the prior parameters. If the evidence for DE dynamics is larger than that for $\Lambda$CDM for a broad range of values of $a_{c}$ and $\sigma_{\bar{X}}$, i.e. does not require one to optimize them to improve the evidence, then one could say that dynamical DE is favoured by observations. One can also define the evidence-weighted reconstruction, in which departures from $X(a)=1$ with low evidence get suppressed (see (Wang et al., 2018) for details).

Fig. 1 shows the DE density reconstruction performed with the “standard” choice of the prior, $\sigma_{\bar{X}}\equiv 0.04$, $a_{c}=0.06$ (in green) along with the evidence-weighted reconstruction (in blue). They show two apparent trends: an overall increase in the effective DE density and an oscillatory behaviour at $a\gtrsim 0.6$. The increase is driven by the local measurements of the Hubble constant $H_{0}$, whose larger value could be interpreted as an increase in DE density. The measurement of the BAO scale from the Lyman-$\alpha$ forest, which prefers a lower $H(z)$ at $z\sim 2.3$, further contributes to the same trend111The reconstruction in (Wang et al., 2018) was based on the Ly$\alpha$BAO analysis of BOSS DR11 quasars (Delubac et al., 2015) that showed a $2.5\sigma$ deviation from the best fit $\Lambda$CDM. The tension has since been reduced to $1.7\sigma$ with the eBOSS DR14 Ly$\alpha$BAO analysis performed in de Sainte Agathe et al. (2019); Blomqvist et al. (2019). Note that, although the tension is lower for the eBOSS sample, which has roughly 20% more Lyman-$\alpha$ sources than BOSS, it does not necessarily mean that the tension will further decrease with a larger sample. The source of the tension, which could be new physics or a yet unknown systematic error, remains unknown, and forthcoming experiments, such as DESI, will help to clarify the issue.. Oscillations, on the other hand, are caused by the combination of the tomographic BAO and the JLA SNe Ia data which happen to have matching oscillatory patterns.

One can see that the apparently large deviation from $X(a)=1$ at high redshifts, seen in the “standard” reconstruction in Fig. 1, is not present in the evidence-weighted curve. The ability of data to constrain DE at $z>3$ is very weak and the reconstruction there is almost completely determined by the prior. This implies no Bayesian evidence for large deviations at high $z$, although the data still prefers a modest increase in DE density.

The lower panels in Fig. 1 show the corresponding effective DE EOS $w_{\rm DE}^{\rm eff}(a)$. They are obtained by generating an ensemble of $X(a)$ from its mean and the covariance matrix and, for each realization, evaluate $w_{\rm DE}^{\rm eff}(a)$ from the conservation of the effective DE fluid. Averaging over the ensemble gives the mean and the uncertainty in $w_{\rm DE}^{\rm eff}(a)$ shown in the plots. If a sampled $X(a)$ happens to have a $|X(a)|<10^{-5}$, we replace it with $X(a)=10^{-5}$ to prevent a singularity in $w_{\rm DE}^{\rm eff}$. As expected, the uncertainty in $w_{\rm DE}^{\rm eff}$ is very large at high redshifts in the case of the standard prior (left panel). This is because $w_{\rm DE}^{\rm eff}(a)$ is determined by the derivative of $X(a)$ and each sampled $X(a)$ can fluctuate within the range allowed by the variance. In the case of the evidence-weighted reconstruction (right panel), $X(a)$ is a linear superposition of many reconstructions obtained with different priors. The different priors all prefer $X(a)$ to be constant, without enforcing any particular value of the constant. Thus, while there is $\sim$15% uncertainty around the value of $X(a)$ at high $z$, its derivative is zero with a much higher certainty, which explains why the uncertainty in $w_{\rm DE}^{\rm eff}(a)$ is so small near $a=0$.

Interestingly, as shown in (Wang et al., 2018), the Bayesian evidence ($\Delta\ln E$) for the oscillatory features is positive at $2.8\sigma$, and they appear equally prominently in both reconstructions in Fig. 1. We also note that, although the Bayesian evidence for dynamical DE is weak, it has increased over the years, with the dynamical pattern being largely consistent with the reconstruction performed in 2012 Zhao et al. (2012).

As we will see later, the oscillatory features in the reconstructions can, in certain circumstances, trigger fast-growing instabilities in cosmological perturbations. Also, the oscillatory pattern and the overall increase in DE density are driven by entirely different datasets. For this reason, we have also considered $X(a)$ obtained by fitting a monotonic function to the reconstructions, which capture the overall increase but do not allow for oscillations. We take the form to be

[TABLE]

where the parameter $D$ is chosen such that $X_{\rm fit}(a=1)=1$. The fitted functions and the corresponding DE EOS are shown with black solid and dashed lines in Fig. 1. Thus, in what follows, we will consider four $X(a)$ histories:

• H1: using the standard prior (green line, Fig. 1);
• H1F: the monotonic fit to H1 (black solid line);
• H2: evidence-weighted reconstruction (blue line);
• H2F: the monotonic fit to H2 (black dashed line).

As an increasing effective DE density cannot be realized in simple quintessence models, one is prompted to consider more complex gravity theories. In the next section we explore the GBD theories as a possible framework for explaining the observed DE dynamics.

III Generalized Brans-Dicke theories and ways to reconstruct them

The non-minimal coupling of the scalar field in the GBD theories could explain the observed “ghostly” behaviour of the effective DE density. We stress that, in this context, the phantom dynamics is only an apparent phenomenon perceived by a cosmologist fitting the conventional Friedmann equation to data while being unaware of the non-minimal coupling.

The GBD action can be written as (Bergmann, 1968; Nordtvedt, 1970; Wagoner, 1970)

[TABLE]

where $m_{0}\equiv(8\pi G)^{-1/2}$ is the Planck mass in terms of the Newton’s constant $G$ measured on Earth, $\phi$ is the extra scalar degree of freedom, $(\partial\phi)^{2}\equiv g^{\mu\nu}\partial_{\mu}\partial_{\nu}\phi$, $U(\phi)$ is the GBD potential and $S_{m}$ denotes the action for the matter fields $\chi_{i}$ minimally coupled to the (Jordan frame) metric $g_{\mu\nu}$. We set $K(\phi)=1$, as one can always do so by a redefinition of $\phi$. The modified Einstein equations are obtained by varying the action with respect to (w.r.t.) the metric $g_{\mu\nu}$,

[TABLE]

where $\nabla_{\mu}$ denotes the covariant derivative w.r.t. the coordinate $x^{\mu}$, $\square\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$, $T_{\mu\nu}^{\rm m}$ is the matter energy-momentum tensor and

[TABLE]

The equation of motion for the scalar field $\phi$ is then obtained by extremizing the action (4) w.r.t. variations of the field $\phi$,

[TABLE]

where the subscript $\phi$ denotes a derivative w.r.t. $\phi$. For convenience we redefine the field, $\smash{\phi\to\phi/m_{0}}$, to make it dimensionless, and the potential, $\smash{U\to U/m_{0}^{2}}$, with the latter measured in ${\rm Mpc}^{-2}$ in agreement with the units convention in CAMB (Lewis et al., 2000).

The freedom in choosing the two functions $F(\phi)$ and $U(\phi)$ translates into the ability of GBD theories to reproduce any expansion history. In particular, the effective DE density $X(a)$ defined in Eq. (1) can increase, and the effective DE EOS can cross $-1$. To see this, we re-write the modified Einstein equation (5) as

[TABLE]

where, in the second line, we have defined the effective DE stress-energy by absorbing into it all the terms on the right hand side other than the usual matter term, i.e.,

[TABLE]

In a flat FRW Universe, the effective DE density is

[TABLE]

with the dot standing for a derivative w.r.t the conformal time, while the effective DE pressure is

[TABLE]

The $\mu=\nu=0$ component of Eq. (8) gives the Friedmann equation,

[TABLE]

which can be recast in the form of Eq. (1).

Note that, by construction, the effective DE “fluid” is conserved, but its EOS,

[TABLE]

is not always well-defined because $\rho^{\rm eff}_{\rm DE}$ in the denominator is allowed to change sign due to the new terms generated by the non-minimal coupling $F(\phi)$. Thus, as previously noted in (Carroll et al., 2005; Das et al., 2006), observing $\smash{w^{\rm eff}_{\rm DE}<-1}$, or finding that $\rho^{\rm eff}_{\rm DE}$ changes its sign, could be a smoking gun of interactions in the dark sector.

The idea of reconstructing the GBD Lagrangian from a given expansion history was previously explored in (Boisseau et al., 2000; Esposito-Farese and Polarski, 2001; Perivolaropoulos, 2005), motivated by the fact that the Hubble function $H(a)$ inferred from the supernovae data available at that time showed a preference for an effective phantom DE equation of state, $\smash{w^{\rm eff}_{\rm DE}<-1}$. As they have shown, one can, in principle, reconstruct both functions $F(\phi)$ and $U(\phi)$ if, in addition to $H(a)$, one knows the evolution of the growth of the matter density contrast $\delta(a)$. Another interesting example is the $f(R)$ gravity where the only unknown function is the function $f$ itself and the full reconstruction can be done with the sole knowledge of the expansion history $H(a)$ (Song et al., 2007; Pogosian and Silvestri, 2008).

In the present work we adopt a slightly different approach. Since the growth of perturbations is rather complicated to extract in a model-independent way because of the redshift-space distortions, non-linearities, bias, etc., we attempt to reconstruct only one of the functions, namely $U(\phi)$, while the other is chosen to either have a given functional form $F(\phi)$ (Model 1), or a given parameterized time-dependence $F(a)$ (Model 2). We will analyze these two cases separately.

While exploring the parameter space, which includes the initial conditions for the scalar field, we restrict to solutions in which the net change in $F(\phi)$ is under 10%, to satisfy the BBN constraints on the variation of the Newton’s constant. We also check for various types of instabilities using the procedure implemented in EFTCAMB (Hu et al., 2014a; Raveri et al., 2014). Specifically, we check for ghost, gradient and mass instabilities discussed in detail in Frusciante et al. (2019) and briefly reviewed below.

After expanding the action up to the second order in perturbations of the metric and matter fields, and removing spurious degrees of freedom, one can isolate the action for the propagating scalar and tensor degrees of freedom De Felice et al. (2017). The conditions for avoiding instabilities can then be formulated in Fourier space in terms of the corresponding kinetic, gradient and mass matrices as follows:

no-ghost: a ghost instability develops when the kinetic term of a field is negative. In the presence of multiple propagating degrees of freedom, a positive definite kinetic matrix guarantees that no ghosts will develop. In practice, this requirement needs to be imposed only at high energies, i.e.  in the high-$k$ limit, since an infrared ghost does not lead to catastrophic instabilities Gumrukcuoglu et al. (2016). 2. 2.

no-gradient: gradient instabilities arise in the high-$k$ regime when the speed of propagation is imaginary. The sound speeds of the propagating degrees of freedom can be identified from the dispersion relations that result from the quadratic action after diagonalizing the kinetic matrix. In order to avoid gradient instabilities we impose $c_{s}^{2}>0$ for all the degrees of freedom. 3. 3.

no-tachyon: whenever the mass matrix of the Hamiltonian contains a negative eigenvalue, the mass instability plagues the low-$k$ regime with the development of a tachyon De Felice et al. (2017). The rate of the instability needs to be taken into account. We will assume that the GBD Hamiltonian exhibits a tachyonic instability when at low momenta a mass eigenvalue $\mu_{i}$ becomes negative and evolves rapidly, i.e. $|\mu_{i}|\gg H^{2}$. Thus, for a theory to be viable, we require $\mu_{i}>0$ or, alternatively, $|\mu_{i}|\lesssim H^{2}$.

This set of conditions was shown to guarantee stability over the full range of linear scales Frusciante et al. (2019) and was implemented in a private version of EFTCAMB.

Public versions of EFTCAMB, as well as other Einstein-Boltzmann solvers like HiClass Zumalacrregui et al. (2016), do not contain the mass condition. Instead, in addition to checking for the no-ghost and no-gradient instabilities, they impose a set of mathematical conditions that prevent the development of exponentially diverging solutions. The latter are worked out from the linear order equation for the scalar field perturbation and are meant to protect against the mass instabilities as well as the ghost and gradient instabilities that could have possibly evaded the checks based on some approximations necessary in setting the conditions. When a mathematical condition is violated, one cannot easily tell which of the three types of instabilities was responsible. In our analysis we used both methods. Namely, we checked for the ghost, gradient and mass stability conditions, as well as using the publicly available stability check that combines the ghost, the gradient and the mathematical conditions.

III.1 Model 1: Reconstructing GBD for a given $F(\phi)$

Given the functional form of $F(\phi)$, we can reconstruct $U(\phi)$ from a given expansion history. We take

[TABLE]

which is a form motivated by high energy theories, e.g. a non-minimally coupled dilaton field representing compactified extra dimensions with the dimensionless parameter $\xi$ controlling the coupling strength.

We begin by writing the two Friedmann equations as

[TABLE]

where

[TABLE]

and ′ denotes derivatives w.r.t. $N\equiv\ln a$. Eq. (15) can be rewritten as an equation for the background evolution of $\phi$:

[TABLE]

where $\smash{E_{m}\equiv\rho_{m}/\rho_{\rm crit}^{0}}$ includes CDM and baryons, $\smash{E_{r}\equiv\rho_{r}/\rho_{\rm crit}^{0}}$ includes photons and massless neutrinos and $\smash{E_{\nu}\equiv\rho_{\nu}/\rho_{\rm crit}^{0}}$ includes massive neutrinos species only. Eq (18) can be solved given the functional form (13) of $F(\phi)$ and the DE density evolution $X(a)$. Given the solution, $\phi(a)$, one can find the potential $U(a)$ from Eq. (14), namely:

[TABLE]

If $\phi(a)$ is monotonic, it can be inverted to obtain $a(\phi)$ and, thus, $U(\phi)$ for the range of $\phi$ covered by the evolution.

Solving Eq. (18) requires setting the value of the field $\phi_{\rm ini}$ and its derivative $\phi_{\rm ini}^{\prime}$ at some initial time $a_{\rm ini}$. To preserve the success of $\Lambda$CDM in explaining the Big Bang Nucleo-synthesis (BBN) and the peak structure of the CMB spectrum, we assume that gravity was close to GR at early times, so that $\smash{F(\phi)=1}$ for $\smash{a\leq a_{\rm ini}}$, but could start deviating from unity at later times. For $\smash{F(\phi)=\exp(\xi\phi)}$ this means $\smash{\phi_{\rm ini}=0}$ and, to explain the features in the reconstructed DE density discussed in the previous section, we will need $\smash{a_{\rm ini}\lesssim 0.1}$. Thus, in addition to providing $X(a)$, we have to specify three parameters: $\xi$, $a_{\rm ini}$ and $\phi_{\rm ini}^{\prime}$.

For the H1, H1F and H2 background histories, reconstructed Model 1 theories contain fast-growing mode instabilities for all choices of initial conditions. For the perturbations around these backgrounds we find that both the mass and the mathematical conditions are not satisfied. It appears that the large rapid increases in $X(a)$ present in H1, H1F and H2 drive the solution towards instability, which could, in principle, be prevented by an appropriate choice of $F(\phi)$. However, the Model 1 coupling function $F(\phi)=\exp(\xi\phi)$ is monotonic and is unable to prevent the onset of instability.

In the case of H2F, which is monotonic and with a relatively small change in $X(a)$, we are able to find viable solutions despite finding negative mass eigenvalues. There, the tachyonic instability corresponding to the negative mass eigenvalues develops on time scales comparable to the Hubble rate allowing for growth of cosmic structure that is in reasonable agreement with observations.

To summarize, we find that Model 1 reconstructions from all four expansion histories are plagued by a mass instabilities. For H1, H1F and H2, this instability develops on time scales small enough for the mathematical condition to detect diverging solutions. In the case of H2F the characteristic time scale is longer and, while our approximate bound of $|\mu_{i}|<H^{2}$ is not satisfied, the instability does not develop to the point of giving diverging solutions.

Fig. 2 shows the allowed range of the non-minimal coupling function $\Omega(a)\equiv F(\phi(a))-1$ for GBD theories reconstructed from H2F. It is obtained by uniformly sampling parameters $(\phi^{\prime}_{\rm ini},\log_{10}a_{\rm ini},\xi)$ from the intervals

[TABLE]

solving for the evolution of $\phi$ and selecting solutions that have $\Omega$ within the allowed range and satisfy the stability condition. The shaded regions in Fig. 2 indicate the confidence level (CL) for having a particular value of $\Omega$ at a given $a$, while the dark line in the middle shows the mean. Examining the numerical solutions, we find that, as expected, the increasing effective DE density drives the field to negative values, resulting in $F(\phi)<1$ and a larger $G_{\rm eff}\propto G/F(\phi)$.

For illustration, in Fig. 3 we show the potential $U(\phi)$ for four GBD theories reconstructed from the H2F DE density with $\phi^{\prime}_{\rm ini}/m_{0}=0$ and $a_{\rm ini}=10^{-2}$ and $\xi=0.5$, $1.5$, $3$ and $10$ respectively. We can see how in all four cases the potential has a cusp at the origin. Stronger couplings leads to steeper potentials. Their shapes resemble the potential in chameleon-like models Khoury and Weltman (2004); Brax et al. (2004), $V(\phi)\propto|\phi|^{-n}$, although the dynamics here is completely different. In the chameleon model, the field tracks the minimum of the effective potential, with the coupling function $F(\phi)$ slowly increasing with the evolution. In our reconstructed theories, the field $\phi$ starts at the top of the cusp and rolls down the potential, with $F(\phi)$ decreasing as it rolls.

III.2 Model 2: Reconstructing GBD with a parameterized $F(a)$

We now change the approach and, instead of working with a given $F(\phi)$, we directly specify the time-dependence of $F$, i.e. $F(a)$. A similar approach was used in (Esposito-Farese and Polarski, 2001; Perivolaropoulos, 2005) to reconstruct the GBD Lagrangian from the expansion history inferred from an early SNe Ia dataset. In this case, we start by writing the modified Friedmann equations as

[TABLE]

where primes denote derivatives w.r.t. the scale factor and overdots denote derivatives w.r.t. the conformal time $\tau$. We can then use (21) to eliminate the potential $U$ in (22) to write

[TABLE]

One can then solve the above equation for $\dot{\phi}$ and use it in (21) to obtain a solution for $U(a)$:

[TABLE]

With the known $U(a)$, one can solve for the kinetic energy $\dot{\phi}^{2}$ from (21) and complete the solution by solving the differential equation to find $\phi(a)$. With the field $\phi(a)$ known, one can convert $U(a)$ and $F(a)$ into $U(\phi)$ and $F(\phi)$, thus reconstructing the functional form of the theory for the range of $\phi(a)$ covered by the solution.

To explore a broad range of possible $F(a)$ histories we adopt a polynomial parametric form

[TABLE]

with coefficients $\alpha_{i}$ sampled uniformly from

[TABLE]

This range is chosen to favour positive values of $F(a)$ close to unity as required by existing bounds.

With $H(a)$ and $F(a)$ specified, one can use EFTCAMB (Hu et al., 2014b), as described in the next Section, to compute the cosmological observables.

We have generated samples of $F(a)$ using the parameterized form (25) and performed reconstructions of the GBD theories for each of the four $X(a)$ histories shown in Fig. 1. The viable ranges of $F(a)$ functions in each case are shown in Fig. 4. One can see that for H1, H1F and H2F, in which $X(a)$ has a large increase or is non-monotonic, $F(a)$ must be non-monotonic to avoid instabilities. In the case of H2F, which has a gently increasing monotonic $X(a)$, a monotonic $F(a)$ has a small probability, but is not excluded.

To reconstruct $U(\phi)$ and $F(\phi)$ one needs to solve a first order ODE for $\phi$, which requires specifying a boundary condition, such as the value of the field at an initial time $a_{\rm ini}$. This means we can only reconstruct $U(\phi)$ and $F(\phi)$ up to an arbitrary shift in the value of $\phi$. The shift has no physical significance, as all the observables are already fully determined. Hence, without loss of generality, we take $a_{\rm ini}=0.001$ and $\phi(a_{\rm ini})=0$.

Fig. 5 shows the non-minimal coupling function $\Omega(\phi)=F(\phi)-1$ and the potential $U(\phi)$ for four representative theories reconstructed from H1, H1F, H2 and HF2. We see that $\Omega(\phi)$ is non-monotonic in these representative cases. The potentials have a runaway shape, being seemingly unbounded from below for large values of the field, although one should keep in mind that the shape is only known over the range covered by the evolution of the field. One can also see small bumps in the potentials derived from H1 and H2, needed to accommodate oscillations in $X(a)$.

IV Cosmological Observables in reconstructed GBD theories

Next we investigate the cosmological implications of the reconstructed GBD theories by computing the CMB anisotropy and the matter power spectra, along with the cross correlation of the CMB temperature and Galaxy Number Counts (GNC), which probes the ISW effect. It is relatively straightforward to calculate these observables using EFTCAMB (Hu et al., 2014a; Raveri et al., 2014)222https://github.com/EFTCAMB/EFTCAMB, which is an implementation of the effective field theory of dark energy (EFTofDE) (Gubitosi et al., 2013; Bloomfield et al., 2013) in the popular Boltzmann solver CAMB (Lewis et al., 2000). In the EFTofDE approach, the most general action for the cosmological background and perturbations in scalar-tensor theories can be written in the unitary gauge, in which the scalar field is uniform on hypersurfaces of constant time, as an expansion in increasing rank-ordered operators invariant under spatial diffeomorphisms. The time-dependent expansion coefficients are referred to as the EFT functions. The part of the EFT action of relevance to the GBD theories is

[TABLE]

where $\tau$ is the conformal time, $\delta g^{00}=g^{00}+1$ is the metric tensor perturbation, and $\Omega$, $\Lambda$ and $c$ are the EFT functions. The GBD theories reconstructed in the previous section can be mapped onto the EFT formalism via

[TABLE]

With this mapping we can use EFTCAMB to compute the CMB spectra and other cosmological observables.

As the initial time of the reconstruction $a_{\rm ini}$ is after the epoch of recombination, the only differences in the CMB anisotropy spectrum compared to the $\Lambda$CDM model can be due to the change in the expansion history, that modifies the distance to last scattering and shifts the positions of the peaks and troughs, and due to the different evolution of the gravitational potentials which affects the late-time ISW contribution to anisotropy. In the latter case, the phenomenology of GBD theories has three competing effects. First, just like in the case of $\Lambda$CDM, the accelerating expansion causes a decay of the metric potentials. In our reconstructed expansion histories the matter dominated era lasts longer, thus delaying the decay of the potentials. Secondly, the background value of the effective Newtons constant that determines the rate of gravitational clustering is $G/F(\phi)$, and can be larger or smaller than $G$, depending on the dynamics of the coupling function, correspondingly increasing or decreasing the rate at which the metric potentials evolve. Thirdly, the scalar field mediates a fifth force on scales smaller that the Compton wavelength of the field, which enhances the growth of the potentials. It is practically impossible to isolate these effects in the CMB anisotropy spectrum, since it only probes the square of the overall integral of the ISW signal. However, one can learn more by studying the correlation of CMB temperature with galaxy distribution at different redshifts Crittenden and Turok (1996); Afshordi (2004). In particular, a characteristic signature of the fifth force would be a negative galaxy-CMB correlation at high redshifts, where one normally expects no ISW signal. A change in the background value of the gravitational coupling could show as either a positive or negative signal, depending on its evolution.

The CMB temperature and GNC cross correlation angular power spectrum can be written as

[TABLE]

where the ISW transfer function is given by

[TABLE]

and the GNC transfer function is given by

[TABLE]

In the above, $\Phi$ and $\Psi$ are the Newtonian gauge metric potentials in Fourier space, $\delta(k,\tau)$ is the matter density contrast, $W(z)$ is the window function that selects galaxies in the given redshift range, and $b_{g}$ is the galaxy bias. The term “corrections” in Eq. (33) denotes collectively the redshift-space-distortion corrections, lensing terms, and other contributions suppressed by $\mathcal{H}/k$ (Challinor and Lewis, 2011). The cross-correlation spectra are then computed using the EFTCAMB patch for CAMB sources 333The latest EFTCAMB patch is not yet compatible with the latest CAMB. In its last update, CAMB and CAMB sources have been merged, so we used the last available iteration of CAMB sources at https://github.com/cmbant/CAMB/tree/CAMB_sources. (Challinor and Lewis, 2011; Lewis and Challinor, 2007).

Since we are not interested in fitting the parameters of the GBD theories to data, but rather in investigating the qualitative features of the ISW effect we choose to show the cross-correlation in three Gaussian windows $W_{1}$, $W_{2}$ and $W_{3}$ centred at redshifts $z_{1}=0.5$, $z_{2}=1$ and $z_{3}=3$. The widths of the window functions are $\sigma_{1}=0.05,\sigma_{2}=0.1$ and $\sigma_{3}=0.5$. The galaxy bias $b_{g}$ is, in general, time and scale dependent. On large scales, relevant for the cross-correlation with CMB, one expects the scale dependence of the bias to be weak and the time dependence to have a simple polynomial dependence (see (Ade et al., 2016b)). The bias is degenerate with the ISW amplitude, but one can calibrate it by jointly studying the GNC auto-correlations and the cross-correlations between GNC and galaxy lensing. As we are only interested in demonstrating the general features of the ISW signal, we fix the galaxy bias to $b_{g}=1$.

IV.1 Observables for Model 1

In the left panel of Fig. 6 we show the distribution of CMB temperature anisotropy spectrum ${\cal D}_{\ell}\equiv\ell(\ell+1)(2\pi)^{-1}C_{\ell}$ for Model 1 theories reconstructed from the H2F DE density obtained by sampling the parameter space as described in Sec. III.1. The shaded regions represent the CL regions to find ${\cal D}_{\ell}$ in the corresponding range, while the white lines show the mean values. In this sampling procedure we used the cosmological parameters obtained in the reconstruction of $X(a)$ in Wang et al. (2018), except for the parameters setting the primordial power spectrum which were not constrained in Wang et al. (2018), and for which we used the best fit $\Lambda$CDM values. The light green band shows the irreducible statistical uncertainty in ${\cal D}_{\ell}$ due to cosmic variance based on the $\Lambda$CDM model. As the ${\cal D}$ measured by Planck are cosmic variance limited over most of the cosmologically relevant $\ell$ Akrami et al. (2018), the shown uncertainty is representative of current data.

As expected, we observe a modified ISW effect at small $\ell$. The small differences in the high-$\ell$ part of the spectra are mainly due to the different distance to the last scattering surface (because of the different expansion history) which causes a shift in the peaks, and also because of the different baryon and CDM densities $\Omega_{b}h^{2}$ and $\Omega_{c}h^{2}$ in the $X(a)$ vs $\Lambda$CDM cases. These high-$\ell$ differences are well-within the cosmic variance band and would likely be accommodated by adjusting other parameters in a comprehensive MCMC parameter estimation.

The right panel of Fig. 6 shows the linear matter power spectrum. First of all, one can note an overall shift upwards for the GBD theories. At early times, before DE begins to dominate the background dynamics, the Planck best-fit $\Lambda$CDM model has more DE density than the GBD models with the reconstructed DE. This means that in the GBD models the matter dominated era (MDE) lasts slightly longer than in the $\Lambda$CDM model, allowing matter to cluster more, hence the overall shift upwards of the matter power spectrum. As in the case of the CMB spectrum, we expect that this difference can accommodated by adjusting other parameters in a comprehensive fit which, however, is beyond the scope of this work. In addition to the change in the matter-DE equality, $P(k)$ is also effected by the larger $G_{\rm eff}$ and the fifth force mediated by the scalar field. This is encoded in the way the deviations from $\Lambda$CDM increase on smaller scales. Finally the oscillations that we note at $k\approx 0.1$ $h/{\rm Mpc}$ are due to the different position of the BAO scale due to a slightly different expansion history of the GBD models.

In Fig. 7 we show the theoretical prediction of the cross correlations for the two classes of reconstructed GBD theories from H2F. Also shown is the cosmic variance statistical uncertainty in the cross-correlation predicted by the $\Lambda$CDM model. As one can see, in some Model 1theories, the ISW effect can become negative signalling a growing gravitational potential due to the fifth force mediated by the extra scalar field. At lower redshifts, when the effective DE becomes larger and the growth of the gravitational potential is overcome by the decay induced by the accelerated expansion, the ISW term is mainly positive.

Some of the Model 1 theories reconstructed from the H2F expansion history are cosmologically viable, at least from the perspective of fitting the CMB spectra. Fig. 8 shows the cosmological observables for the four representative models whose reconstructed potentials $U(\phi)$ were shown in Fig. 3. While the CMB anisotropies are almost the same for each model, they differ considerably in the clustering of matter and this is also noticeable in the cross correlations $C_{\ell}^{Tg}$ at the bottom panels. In the higher redshift windows, the larger values of the couplings constant $\xi$ drive a growth of the gravitational potentials $\Psi$ and $\Phi$ due to the fifth force mediated by the scalar field, causing a negative ISW effect. When DE eventually starts dominating the potentials stop growing and instead decay, turning the sign of the ISW effect.

IV.2 Observables for Model 2

Fig. 9 shows the distribution of the CMB anisotropy spectra corresponding to the sampled Model 2 theories reconstructed from the H1, H1F, H2 and H2F expansion histories. We see in all cases there is a preference for a large ISW contribution to ${\cal D}_{\ell}$. This is especially the case for H1, in which $X(a)$ is non-monotonic and has a large increase. However, since cosmic variance results in large statistical error bars at small $\ell$, there are models on the fringe of the allowed range for H1F, H2 and H2F that can be compatible with the current data.

The Model 2 distribution of $P(k)$ is in good agreement with the data, with the best fit $\Lambda$CDM prediction being well inside the 68% CL, as shown in Fig. 10 in the case of H2F. The matter power spectra in the cases of H1, H1F and H2F are very similar.

A very distinctive observational feature of Model 2 models is a large positive ISW signal at high redshifts, as seen in Fig. 11. This is caused by $F>1$ at $z\gtrsim 1$, which decreases the effective Newton’s constant appearing in the Poisson equation ($G_{\rm eff}\propto G/F$) resulting in a suppression of gravitational potentials during the matter dominated epoch. The enhancement in the high redshift cross-correlation is well in excess of the cosmic variance uncertainty around the $\Lambda$CDM model and would be detectable with the next generation large scale structure surveys such as DESI, LSST and Euclid. We note that recent redshift space distortion measurements slightly favour a lower value of $G_{\rm eff}$ Zhao et al. (2019).

In Fig. 12 we show observables corresponding to the four models in Fig. 5 representing reconstructions using H1, H1F, H2 and H2F. We can see the similarity in general trends, with features being the most pronounces in the case of H1 and less so for H2F. However, in all cases, there is a large positive ISW signal at high redshifts which would be a smoking gun of GBD models with a non-monotonic $F(\phi)$.

V Conclusions

Current observations favour an increasing effective DE density, corresponding to an effective DE EOS that is less than $-1$ (Zhao et al., 2017a; Wang et al., 2018). Such apparently phantom behaviour of DE can also occur in GBD theories, as a manifestation of the additional interaction mediated by the scalar field.

We have set up a reconstruction method to design the Lagrangians of GBD type scalar-tensor theories corresponding to expansion histories extracted using the latest data probing the background (Wang et al., 2018). We then examined the viability of such designer GBD theories, both in terms of their stability and their ability to predict acceptable cosmological observables.

We found that a large increase in the effective DE density, or the apparent oscillatory dynamics also favoured by the data, are difficult to accommodate within a GBD theory with a monotonically evolving coupling function, such as $F(\phi)\propto\exp({\xi\phi})$. However, allowing for an arbitrary $F(a)$, parametrized in terms of a polynomial expansion, results in GBD theories capable of fitting current CMB and matter power spectra.

We find that, in viable models, $F(a)$ increases at high redshifts before decreasing at more recent epochs, leading to a smaller effective gravitational coupling $G_{\rm eff}$ at redshifts $z\gtrsim 1$, and a larger $G_{\rm eff}$ at $z<1$. This leads to a robust prediction of a large positive ISW signal at $z>1$, which would be readily detectable through CMB-galaxy cross-correlation using high redshift sources from DESI, LSST and Euclid.

In our analysis, we opted to provide functions $F(\phi)$ or $F(a)$, and reconstruct the potential $U(\phi)$. One could, alternatively, opt to find $F(\phi)$ for a given $U(\phi)$. We expect that, regardless the choice, the main conclusion about the key role of the ISW effect in falsifying GBD theories will remain the same.

The method developed here is complementary to the reconstruction of the EFT functions (including the GBD subset of the EFTofDE) from cosmological observations performed in Raveri (2019). In that work, the expansion history was reconstructed in conjunction with the scalar-tensor Lagrangian, thus only producing expansion histories that are consistent with the GBD. Our approach is different – we start with an expansion history obtained from the data in a largely model-independent way and checked if there can be GBD theories producing it. The difference is that a joint reconstruction within the GBD framework could miss expansion histories that are difficult to accommodate with smooth EFT functions, potentially missing a hint for dynamics that would correspond to a rare realization of GBD.

Our results show that one could rule out scalar-tensor theories as the explanation of departures from the $\Lambda$CDM background expansion history using purely cosmological data sets. This is particularly important for testing theories in which the scalar field couples only to dark matter, to which the tight laboratory and solar system tests of gravity do not apply.

Acknowledgements.

We acknowledge G. Papadomanolakis and S. Peirone for useful discussions and feedback on our findings. A.Z. and L.P. are supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada. A.S. 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. GBZ is supported by the National Key Basic Research and Development Program of China (No. 2018YFA0404503), and by NSFC Grants 11720101004 and 11673025. Y.W. is supported by the Nebula Talents Program of NAOC. This research was enabled in part by support provided by WestGrid and Compute Canada.