Brans–Dicke Gravity with a Cosmological Constant Smoothes Out ΛCDM Tensions

With the above notation, we can write down the BD-field equations in the presence of matter and vacuum components as follows:

$$\begin{eqnarray}&&3{H}^{2}+3H\displaystyle \frac{\dot{\psi }}{\psi }-\displaystyle \frac{{\omega }_{\mathrm{BD}}}{2}{\left(\displaystyle \frac{\dot{\psi }}{\psi }\right)}^{2}=\displaystyle \frac{8\pi }{\psi }{\rho }_{{\rm{T}}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&2\dot{H}+3{H}^{2}+\displaystyle \frac{\ddot{\psi }}{\psi }+2H\displaystyle \frac{\dot{\psi }}{\psi }+\displaystyle \frac{{\omega }_{\mathrm{BD}}}{2}{\left(\displaystyle \frac{\dot{\psi }}{\psi }\right)}^{2}=-\displaystyle \frac{8\pi }{\psi }{p}_{{\rm{T}}}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\ddot{\psi }+3H\dot{\psi }=\displaystyle \frac{8\pi }{2{\omega }_{\mathrm{BD}}+3}({\rho }_{{\rm{T}}}-3{p}_{{\rm{T}}}).\end{eqnarray} \tag{ 3 }$$

For constant ψ, the first two equations reduce to the Friedmann and pressure equations of GR, respectively, with $G\,=1/\psi \,={G}_{N}$, and the third requires $| {\omega }_{\mathrm{BD}}| \to \infty$ for consistency. By combining the above equations one finds a local covariant conservation law that is identical to that of GR since there is no interaction between matter and the BD-field, namely ${\dot{\rho }}_{T}+3H({\rho }_{{\rm{T}}}+{p}_{{\rm{T}}})=0$. Being ${\rho }_{{\rm{\Lambda }}}=$ const. and owing to the equation of state (EoS) of vacuum (${p}_{{\rm{\Lambda }}}=-{\rho }_{{\rm{\Lambda }}}$), the previous equation reduces to $\dot{\rho }+3H(\rho +p)=0$. Taking into account that we assume separate conservation of the different components, such a law can be conveniently split into a conservation law for each component (baryons, dark matter, neutrinos, and photons). Hereafter we use the following definitions:

$$\begin{eqnarray}&&\varphi (t)\equiv {G}_{N}\psi (t)={G}_{N}/G(t),\qquad {\epsilon }_{\mathrm{BD}}\equiv \displaystyle \frac{1}{{\omega }_{\mathrm{BD}}},\end{eqnarray} \tag{ 4 }$$

where φ is the dimensionless BD-field. In this notation, $G(t)\equiv G(\varphi (t))={G}_{N}/\varphi$ plays the role of effective gravitational coupling in the BD context at the level of the action and field equations. We should not expect G(t) to be equal to G_N at the present time (t₀); thus, $\varphi ({t}_{0})$ is in general close, but not exactly equal, to 1. We can recover GR in the limit ${\epsilon }_{\mathrm{BD}}\to 0$, and the ΛCDM model when ${\epsilon }_{\mathrm{BD}}\to 0$ and $\varphi \to 1$.

It proves revealing to analyze Brans–Dicke cosmology from the standpoint of what we may call the "GR-frame." The latter is obtained by parameterizing the BD-field equations as deviations with respect to the usual Friedmann's equations. For instance, in Solà (2018) and de Cruz Pérez & Solà (2018) it was shown, using an approximate treatment, that such kind of approach leads to mimic "running vacuum," even starting from a cosmological constant. Vacuum dynamics, and in general dynamical dark energy (DE), can be phenomenologically favorable, even if not firmly established yet; see, e.g., Solà et al. 2015, 2018, 2019, Zhao et al. 2017, Park & Ratra 2018, Martinelli et al. 2019, and references therein. It can also be a cure for some of the tensions in the ΛCDM (Di Valentino et al. 2017; Solà et al. 2017). Here we undertake a systematic approach and an exact numerical treatment of the BD-field equations in the GR-frame, which will lead to effective dynamical DE. The first step is to rewrite Equation (1) as if it were the usual Friedmann equation, $3{H}^{2}={\kappa }^{2}({\rho }_{{\rm{T}}}+{\rho }_{\mathrm{BD}}).$ It is not difficult to show that the effective energy density associated to the BD-field, ${\rho }_{\mathrm{BD}}$, reads

$$\begin{eqnarray}&&{\kappa }^{2}{\rho }_{\mathrm{BD}}=3{H}^{2}\,{\rm{\Delta }}\varphi -3H\dot{\varphi }+\displaystyle \frac{1}{2{\epsilon }_{\mathrm{BD}}}\displaystyle \frac{{\dot{\varphi }}^{2}}{\varphi },\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Delta }}\varphi \equiv 1-\varphi$. Recall our definitions (4); we have just included all the energy density terms beyond the ΛCDM in the expression for ${\rho }_{\mathrm{BD}}$. Such an energy density can therefore be interpreted as a new DE component within the GR-frame, which must be added to the vacuum part ${\rho }_{{\rm{\Lambda }}}$. The second step is to recast Equation (2) as the usual Friedmann's pressure equation: $2\dot{H}+3{H}^{2}=-{\kappa }^{2}({p}_{{\rm{T}}}+{p}_{\mathrm{BD}})$. This leads to the following expression for the effective pressure associated to the BD-field (playing the role of "extra DE pressure," in addition to ${p}_{{\rm{\Lambda }}}=-{\rho }_{{\rm{\Lambda }}}$):

$$\begin{eqnarray}&&{\kappa }^{2}{p}_{\mathrm{BD}}=-3{H}^{2}\,{\rm{\Delta }}\varphi -2\dot{H}\,{\rm{\Delta }}\varphi +\ddot{\varphi }+2H\dot{\varphi }+\displaystyle \frac{1}{2{\epsilon }_{\mathrm{BD}}}\displaystyle \frac{{\dot{\varphi }}^{2}}{\varphi }.\end{eqnarray} \tag{ 6 }$$

Direct calculation from the three field Equations (1)–(3) shows that such a fluid obeys the additional conservation law

$$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{BD}}+3H({\rho }_{\mathrm{BD}}+{p}_{\mathrm{BD}})=0.\end{eqnarray} \tag{ 7 }$$

One can also cross check it from ${{\rm{\nabla }}}^{\mu }{\tilde{T}}_{\mu \nu }^{\mathrm{BD}}=0$, where the energy–momentum tensor ${\tilde{T}}_{\mu \nu }^{\mathrm{BD}}\equiv -\tfrac{2}{\sqrt{-g}}\tfrac{\delta {\tilde{S}}_{\mathrm{BD}}}{\delta {g}^{\mu \nu }}$ is computed from the part of the full BD-action ${S}_{\mathrm{BD}}$ that remains after we subtract from it the usual Einstein–Hilbert action ${S}_{\mathrm{EH}}$ (with cosmological term) and the matter action, S_m, i.e., ${\tilde{S}}_{\mathrm{BD}}\,={S}_{\mathrm{BD}}-{S}_{\mathrm{EH}}-{S}_{m}$.

From the previous considerations, the BD-fluid can be described by the quantities ${p}_{\mathrm{BD}}$ and ${\rho }_{\mathrm{BD}}$. However, a more useful picture can be obtained by considering the EoS of the full "effective DE" in this context. Obviously, it must receive contributions from the BD-fluid and the vacuum. Using the field equations in the GR-frame as discussed above, we find ${w}_{\mathrm{eff}}=({p}_{{\rm{\Lambda }}}+{p}_{\mathrm{BD}})/({\rho }_{{\rm{\Lambda }}}+{\rho }_{\mathrm{BD}})$. It can be conveniently cast as follows:

$$\begin{eqnarray}&&{w}_{\mathrm{eff}}(t)=-1+\displaystyle \frac{-2\dot{H}{\rm{\Delta }}\varphi +{f}_{1}(\varphi ,\dot{\varphi },\ddot{\varphi })}{{\rm{\Lambda }}+3{H}^{2}\,{\rm{\Delta }}\varphi +{f}_{2}(\varphi ,\dot{\varphi })}.\end{eqnarray} \tag{ 8 }$$

The two functions ${f}_{\mathrm{1,2}}$ need not be specified here, as they can be easily computed from the previous formulae. It suffices to say that they are numerically negligible, in absolute value, as compared to $\dot{H}{\rm{\Delta }}\varphi$ and ${H}^{2}{\rm{\Delta }}\varphi$ since they depend on time derivatives of the slowly varying function φ. The effective EoS (8) is, obviously, a time-evolving quantity that mimics dynamical DE in the above effective GR picture. Since $| {\rm{\Delta }}\varphi |$ is small, as confirmed by our analysis (see Section 4.2), let us remark the following interesting situation near our time (i.e., for cosmological redshift $z\ll 1$):

$$\begin{eqnarray}&&{w}_{\mathrm{eff}}(z)\simeq -1-\displaystyle \frac{2\dot{H}{\rm{\Delta }}\varphi }{{\rm{\Lambda }}}\simeq -1+{\rm{\Delta }}\varphi \,\displaystyle \frac{{{\rm{\Omega }}}_{m0}}{1-{{\rm{\Omega }}}_{m0}}{\left(1+z\right)}^{3},\end{eqnarray} \tag{ 9 }$$

where we have expanded linearly in ${\rm{\Delta }}\varphi$ and expressed the result in terms of the current value of the matter cosmological parameter ${{\rm{\Omega }}}_{m0}={\rho }_{m0}/{\rho }_{c0}$. Here ${\rho }_{c0}=3{H}_{0}^{2}/{\kappa }^{2}$ is the critical density at present and H₀ is the Hubble parameter. Clearly, for ${\rm{\Delta }}\varphi \gt 0$ (resp. < 0) we meet quintessence-like (resp. phantom-like) behavior. We shall further discuss this important issue in Section 4.2.

We use two data sets, which we believe are helpful to better focus on the origin of the main effects. The first one is labeled DS1 and contains the cosmological data Sn Ia+H(z)+BAO+LSS+CMB+R19; the second one, DS2, is just the subset BAO+LSS+CMB+R19 of the first. The CMB part involves the full Planck 2015 TT+lowP+lensing likelihood (Planck XIII 2016). These two data sets are essentially the same as those described in detail in Solà et al. (2019). Here, however, apart from the cosmic chronometer data on the Hubble rate H(z), we have included the latest local value of the Hubble parameter, ${H}_{0}=74.03\pm 1.42$ km s⁻¹ Mpc⁻¹ (Riess et al. 2019), which we have denoted R19, based on distance ladder measurements. This is an important additional input, given the significant existing tension ($\sim 4.4\sigma$) of such a value with the Planck results, as discussed in that reference. In addition, we have added up the independent Sn Ia data from the DES survey (Abbott et al. 2019) to the Pantheon+MCT data set (Riess et al. 2018; Scolnic et al. 2018) already used in the mentioned reference. In the BAO sector, we include a new data point from the DES survey (Camacho et al. 2019) and the combined Lyα-forest data by (Blomqvist et al. 2019). See Tables 4 and 5 of Solà et al. (2019) for the remaining list of BAO and RSD data. Note that, in the last two sets, we include the matter bispectrum data from Gil-Marín et al. (2017; labeled BSP in the aforementioned tables). Finally, we use the WL data from Hildebrandt et al. (2018).

To compare the theoretical predictions of the models under study with the available observational data, we have made use of the Einstein–Boltzmann code CLASS (Blas et al. 2011) in combination with the Monte Carlo Markov Chain sampler MontePython (Audren et al. 2013). Apart from the common set of basic six parameters of the concordance model, the BD-ΛCDM model involves two more, all of them listed in Table 1. One of the characteristic BD parameters is obviously ${\epsilon }_{\mathrm{BD}}$. However, in order to solve numerically the system of differential Equations (1)–(3) we use a second fitting parameter, namely the initial value of the BD-field, ${\varphi }_{\mathrm{ini}}$, which is set at ${z}_{\mathrm{ini}}={10}^{14}$. The evolution of φ proves to be very mild and we naturally take its derivative at that point to be zero. We have checked that the small variation induced in the expansion rate at BBN is within bounds (Uzan 2011). The main fitting results of our analysis are displayed in Table 1 and Figure 1. In the table, we compare the standard ΛCDM model with its BD-ΛCDM counterpart. Remarkably enough, the contour lines in Figure 1 show a preference for relatively high values of H₀, not far from R19, while keeping ${\sigma }_{8}\equiv {\sigma }_{8}(0)$ at an intermediate level between Planck measurements (Planck XIII 2016) and cosmic shear data (Hildebrandt et al. 2018), thus smoothing out this tension as well. In Figure 1, we provide the matrix containing some relevant combinations of two-dimensional marginalized distributions for fitting parameters of the BD model, together with the corresponding one-dimensional marginalized likelihoods for each parameter. We confirm indeed the softening of the two tensions in view of the fitted values of H₀ around 71–72 km s⁻¹ Mpc⁻¹, which coexist peacefully with a sufficiently low ${\sigma }_{8}\simeq 0.80$. The residual H₀-tension with the local measurement (Riess et al. 2019) comes down to $1.8\sigma$ (DS1) and $1.1\sigma$ (DS2) only.

Zoom In Zoom Out Reset image size

Table 1.  The Mean Fit Values and 68.3% Confidence Limits for the Considered Models (ΛCDM and BD-ΛCDM) using Two Data Sets: (1) DS1, i.e., Sn Ia+H(z)+BAO+LSS+CMB+R19 with Full Planck 2015 CMB Likelihood (First Block); and (2) DS2, Based on the Subset BAO+LSS+CMB+R19 (Second Block)

Parameter	DS1		DS2
	ΛCDM	BD-ΛCDM	ΛCDM	BD-ΛCDM
H0 (km/s/Mpc)	${68.65}_{-0.40}^{+0.38}$	${71.03}_{-0.86}^{+0.91}$	${68.69}_{-0.39}^{+0.38}$	${72.00}_{-1.10}^{+1.00}$
${{\rm{\Omega }}}_{m0}$	0.2955 ± 0.0048	0.2742 ± 0.0077	0.2950 ± 0.0047	0.2665 ± 0.0084
${{\rm{\Omega }}}_{b0}$	0.0476 ± 0.0004	0.0453 ± 0.0010	0.0476 ± 0.0004	0.0443 ± 0.0012
τ	${0.063}_{-0.012}^{+0.010}$	${0.081}_{-0.018}^{+0.015}$	${0.063}_{-0.011}^{+0.010}$	0.084 ± 0.018
ns	${0.9700}_{-0.0040}^{+0.0038}$	${0.9891}_{-0.0082}^{+0.0070}$	0.9704 ± 0.0038	${0.9945}_{-0.0086}^{+0.0081}$
${\sigma }_{8}(0)$	${0.804}_{-0.009}^{+0.007}$	0.801 ± 0.010	${0.804}_{-0.008}^{+0.007}$	${0.803}_{-0.010}^{+0.011}$
${\epsilon }_{\mathrm{BD}}$	⋯	$-{0.00277}_{-0.00154}^{+0.00170}$	⋯	$-{0.00315}_{-0.00175}^{+0.00168}$
${\varphi }_{\mathrm{ini}}$	⋯	${0.924}_{-0.023}^{+0.021}$	⋯	${0.901}_{-0.025}^{+0.026}$
$\varphi (0)$	⋯	${0.904}_{-0.029}^{+0.028}$	⋯	0.879 ± 0.032
${w}_{\mathrm{eff}}(0)$	−1	$-{0.961}_{-0.011}^{+0.012}$	−1	$-{0.951}_{-0.013}^{+0.012}$
$\dot{G}(0)/G(0)({10}^{-13}\,{{yr}}^{-1})$	⋯	${3.149}_{-1.924}^{+1.741}$	⋯	${3.625}_{-1.954}^{+1.994}$
ΔDIC (ΔAIC)	⋯	8.34 (7.72)	⋯	9.89 (9.94)

Note. In all cases a massive neutrino of 0.06 eV has been included. First, we display the fitting values for the conventional free parameters: Hubble parameter, H₀, the total nonrelativistic matter parameter ${{\rm{\Omega }}}_{m0}$, and the baryonic part ${{\rm{\Omega }}}_{b0}$, the reionization optical depth τ, the spectral index n_s of the primordial power spectrum, and, for convenience, instead of the amplitude A_s of the spectrum we list the value of ${\sigma }_{8}(0)$. The specific free parameters of the BD model are ${\epsilon }_{\mathrm{BD}}$ and ${\varphi }_{\mathrm{ini}}$; see the text. Owing to their significance for our analysis, we quote their values with the error bars at 1σ, 2σ, and 3σ c.l., to wit: ${\epsilon }_{\mathrm{BD}}=-{0.00277}_{-0.00154-0.00324-0.00484}^{+0.00170+0.00312+0.00458}$ and ${\varphi }_{\mathrm{ini}}={0.924}_{-0.023-0.042-0.064}^{+0.021+0.044+0.066}$ for the DS1 scenario, and ${\epsilon }_{\mathrm{BD}}=-{0.00315}_{-0.00175-0.00341-0.00488}^{+0.00168+0.00338+0.00515}$ and ${\varphi }_{\mathrm{ini}}={0.901}_{-0.025-0.050-0.080}^{+0.026+0.052+0.077}$ for the DS2 one. Finally, we include the computed values of three parameters at present: the value of the BD-field, $\varphi (0)$, the effective EoS ${w}_{\mathrm{eff}}(0)$, Equation (8), and the relative time variation of G, $\dot{G}(0)/G(0)$. In the last row we compare the fit qualities by displaying the differences of the Deviance and Akaike information criteria (DIC and AIC, respectively) between ΛCDM and its BD counterpart for the two data sets under consideration. BD appears to be preferred by the data according to these criteria. See the text for further details.

Download table as:  ASCIITypeset image

How about the global quality of our fits? Occam's razor principle can be implemented rigorously with various Bayesian model selection tools; see, e.g., Burnham & Anderson (2002). In this Letter, we opt for making use of the Deviance (Spiegelhalter et al. 2002) and Akaike (Akaike 1974) information criteria (DIC and AIC, respectively). When a large amount of data is employed, AIC is simply given by $\mathrm{AIC}={\chi }^{2}(\hat{\theta })+2n$, where n is the number of independent fitting parameters and $\hat{\theta }$ is the collection of their mean values. DIC is a more sophisticated generalization of AIC, being defined as

$$\begin{eqnarray}&&\mathrm{DIC}={\chi }^{2}(\hat{\theta })+2{p}_{D}.\end{eqnarray} \tag{ 11 }$$

Here ${p}_{D}=\overline{{\chi }^{2}}-{\chi }^{2}(\hat{\theta })$ is the effective number of parameters of the model and $\overline{{\chi }^{2}}$ the mean of the overall ${\chi }^{2}$ distribution. DIC is particularly suitable for us, since we can easily compute all the quantities involved directly from the Markov Chains and other output generated with MontePython. Both DIC and AIC are reliable provided the posterior distributions are sufficiently Gaussian. This is actually the case here, as reflected in the elliptic shapes of the two-dimensional contours in Figure 1 and also in the normal-like appearance of the one-dimensional distributions shown there. To compare the ability of the ΛCDM and BD-ΛCDM models to fit the data, one has to compute the respective differences of DIC and AIC values between the first and second models. They are denoted ΔDIC and ΔAIC in Table 1, where we provide the results for both data sets DS1 and DS2. Since these differences are positive and both lie in the interval $5\lt {\rm{\Delta }}\mathrm{AIC},{\rm{\Delta }}\mathrm{DIC}\lt 10$ we conclude, following Burnham & Anderson (2002) and Spiegelhalter et al. (2002), that they point to strong evidence in favor of BD-ΛCDM as compared to ΛCDM.

The following observations are also in order. The fitted ${\epsilon }_{\mathrm{BD}}$ in Table 1 entails a large enough value of $| {\omega }_{\mathrm{BD}}| ={ \mathcal O }(300)$ such as to guarantee that the BD-ΛCDM model remains sufficiently close to the concordance one, but not so large as to make the BD approach phenomenologically irrelevant. At the same time, negative values of ${\epsilon }_{\mathrm{BD}}$ are preferred in our numerical analysis. Both features are consistent with previous estimates from analytical power-law solutions found in Solà (2018) and de Cruz Pérez & Solà (2018). Furthermore, the relative time variation of G at present is found to be positive, $\dot{G}(0)/G(0)\simeq +3\times {10}^{-13}\,{\mathrm{yr}}^{-1}$(at roughly 2σ), which indicates "asymptotically free" behavior (i.e., G mildly increasing with the expansion), which is also consistent with the aforementioned power-law solution.

Last but not least, worth noticing is the behavior of the effective EoS (8), which we plot numerically in Figure 2. Being ${\rm{\Delta }}\varphi \gt 0$ throughout our analysis (see Table 1) we can confirm the expected quintessence-like signal, which we had anticipated with the approximate formula (9) near our time ($z\ll 1$)—see the inner plot in that Figure. As can be seen, a rather conspicuous signal in between 3σ and 4σ is obtained, depending on the data set. In the opposite end ($z\gg 1$), i.e., deep into the matter- or radiation-dominated epochs, the behavior of (8) is of the form ${w}_{\mathrm{eff}}\simeq -1-\tfrac{2}{3}\tfrac{\dot{H}}{{H}^{2}}$ (as Λ and ${f}_{\mathrm{1,2}}$ become negligible against $| \dot{H}| {\rm{\Delta }}\varphi$ and ${H}^{2}{\rm{\Delta }}\varphi$). The characteristic EoS of these epochs is then met in sequence (${w}_{\mathrm{eff}}\simeq 0,1/3$).

Zoom In Zoom Out Reset image size