This proves our contention that the function $\mathcal{F}$ in (7.32) is of
order $\epsilon_{\rm BD}$ and its effects can be treated as a perturbation to
the above formulas. Incidentally, the relative change of $\varphi$ does not
depend on $\varphi$ itself. From the definition of $\mathcal{F}$ we can now
refine the old Friedmann’s equation (7.15) or (7.17) (only valid for
$\epsilon_{\rm BD}=0$) as follows. Starting from Eq. (1.151), it is easy to
see that it can be cast in the Friedmann-like form:

$H^{2}(a)=\frac{8\pi G_{N}}{3\varphi(a)(1-\mathcal{F}(a))}\,\rho(a).$ (7.35)

Despite $\mathcal{F}(a)$ evolves with the expansion, as shown by (7.32), it is
of order $\epsilon_{\rm BD}$ and evolves very slowly. In this sense, Eq.
(7.35) behaves approximately as an ${\cal O}(\epsilon_{\rm BD})$ correction to
Friedmann’s equation (7.15). Setting $a=1$, the value of the current Hubble
parameter satisfies

$H_{0}^{2}=\frac{8\pi G_{N}}{3\varphi_{0}(1-\mathcal{F}_{0})}\rho^{0},$ (7.36)

where $\varphi_{0}\equiv\varphi(a=1)$ and
$\mathcal{F}_{0}\equiv\mathcal{F}(a=1)$. The above equation implies that
$\rho^{0}=\rho_{c}^{0}\varphi_{0}(1-\mathcal{F}_{0})$. We may now rewrite
(7.35) in the suggestive form:

$H^{2}(a)=H_{0}^{2}\left[\hat{\Omega}_{m}(a)a^{-3}+\hat{\Omega}_{\gamma}(a)a^{-4}+\hat{\Omega}_{\nu}(a)+\hat{\Omega}_{\Lambda}(a)\right]\,,$
(7.37)

provided we introduce the new ‘hatted’ parameters $\hat{\Omega}_{i}(a)$, which
are actually slowly varying functions of the scale factor:

Figure 4: Exact numerical analysis of the evolution of $\varphi$ as a function
of the redshift across the entire cosmic history, starting from the radiation-
dominated epoch up to our time. We use here the values of the BD-$\Lambda$CDM
Baseline+$H_{0}$ dataset indicated in the figure itself (cf. Secs. 7.6 and
7.7). In particular, $\epsilon_{\rm BD}=-0.00199^{+0.00142}_{-0.00147}$. The
band around the central (dotted) curve shows the computed $1\sigma$
uncertainty from the Markov chains of our statistical analysis.

$\hat{\Omega}_{i}(a)=\frac{\Omega_{i}}{\varphi(a)}\frac{1}{1-\mathcal{F}(a)}\simeq\frac{\Omega_{i}}{\varphi(a)}\left(1+\mathcal{F}(a)\right)\,,$
(7.38)

with $\Omega_{i}=\rho^{0}_{i}/\rho_{c}^{0}$ as previously. These functions
also satisfy, exactly, the canonical sum rule at present:
$\sum_{i}\hat{\Omega}_{i}(a=1)=1$.For $\epsilon_{\rm BD}=0$, the hatted
parameters reduce to the old tilded ones (7.18),
$\hat{\Omega}_{i}=\tilde{\Omega}_{i}$, and for typical values of
$|\epsilon_{\rm BD}|\sim\mathcal{O}(10^{-3})$ the two sets of parameter differ
by $\mathcal{O}(\epsilon_{\rm BD})$ only:

$\hat{\Omega}_{i}(a)=\frac{\Omega_{i}}{\varphi}+\mathcal{O}(\epsilon_{\rm
BD})=\tilde{\Omega}_{i}+\mathcal{O}(\epsilon_{\rm BD})\,.$ (7.39)

From (7.31) we obviously recover the previous Eq. (7.25) in the limit
$\epsilon_{\rm BD}\to 0$, and it is easy to see that for non-null values of
$\epsilon_{\rm BD}$ the density contrast acquires a dependence on the ratio
$\varphi^{\prime}/\varphi$ and its derivative, so it is sensitive to the
relative change of $\varphi$ with the expansion. Its time evolution is now
possible by virtue of the third BD-field equation (7.5), which can be expanded
linearly in $\epsilon_{\rm BD}$ in a similar way. After some calculations, we
find

$\varphi^{\prime\prime}+\frac{1}{2a}\left(5-3w(a)\right)\varphi^{\prime}=\frac{3\epsilon_{\rm
BD}}{2a^{2}}(1-3w(a))\varphi\,.$ (7.40)

For $\epsilon_{\rm BD}=0$ we recover the solution $\varphi=$const. In the
radiation dominated epoch (RDE), $w\simeq 1/3$, the r.h.s. vanishes and in
this case $\varphi$ need not be constant. It is easy to see that the exact
solution of this equation in that epoch is

$\varphi(a)=\varphi^{(0)}+\frac{\varphi^{(1)}}{a}\,,$ (7.41)

for arbitrary constants $\varphi^{(i)}$. The variation during the RDE is
therefore very small since the dominant solution is a constant and the
variation comes only through a decaying mode $1/a\sim t^{-1/2}$ ($t$ is the
cosmic time). For the MDE (for which $w=0$) there is some evolution, once more
with a decaying mode but then through a sustained logarithmic term:

$\varphi(a)\sim\varphi^{(0)}\left(1+\epsilon_{\rm BD}\ln
a\right)+\varphi^{(1)}a^{-3/2}\rightarrow\varphi^{(0)}\left(1+\epsilon_{\rm
BD}\ln a\right)\,,$ (7.42)

where coefficient $\varphi^{(0)}$ is to be adjusted from the boundary
conditions between epochs161616Approximate solutions to the BD-field equations
for the main cosmological variables in the different epochs are discussed in
Appendix C.. The dynamics of $\varphi$ for $\epsilon_{\rm BD}\neq 0$ is
actually mild in all epochs since $\epsilon_{\rm BD}$ on the r.h.s. of Eq.
(7.40) is small. However mild it might be, the dynamics of $\varphi$ modifies
both the friction and Poisson terms in Eq. (7.31), and it is therefore of
pivotal importance to understand what are the changes that are induced by
positive and negative values of $\epsilon_{\rm BD}$ on these terms during the
relevant epochs of the structure formation history. An exact (numerical)
solution is displayed in Fig. 4, where we can see that $\varphi$ remains
within the approximate interval $0.918\lesssim\varphi\lesssim 0.932$ for the
entire cosmic history (starting from the RDE up to our time). This plot has
been obtained from the overall numerical fit performed to the observational
data used in this analysis within one of the BD-$\Lambda$CDM baseline datasets
considered (cf. Sec. 7.7). The error band around the main curve includes the
$1\sigma$-error computed from our statistical analysis. Two very important
things are to be noted at this point: on the one hand the variation of
$\varphi$ is indeed small, and on the other hand $\varphi<1$, and hence
$G(\varphi)>G_{N}$ for the whole cosmic span.

Figure 5: Left plot: Relative difference between the friction term $T_{1}$ of
Eq. (7.31),
$\delta_{m}^{\prime}\prime+T_{1}\delta_{m}^{\prime}-T_{2}\delta_{m}=0$, using
the best-fit values of the BD-$\Lambda$CDM model obtained with the
Baseline+$H_{0}$ dataset, i.e. with $\epsilon_{\rm BD}=-0.00199$, and the case
with $\epsilon_{\rm BD}=0$ (cf. Table 3 for the values of the other
parameters); Right plot: The same, but for the last term of Eq. (7.31),
$T_{2}$ (Poisson term). We can clearly appreciate that negatives values of
$\epsilon_{\rm BD}$ produce a higher Hubble friction, i.e. a higher $T_{1}$,
and a lower $T_{2}$, w.r.t. the $\epsilon_{\rm BD}=0$ case. Both things lead
to a decrease of $\delta_{m}$. The two plots have been obtained with our
modified version of CLASS. See the main text for further details.

We can start considering what is the influence of the scalar field dynamics on
the perturbations during the pure MDE. Using the relations (7.33) and (7.34)
in (7.31) we find (setting $w=0$, $\rho\simeq\rho_{m}$ and neglecting
$\rho_{\Lambda}$ in the MDE):

$\delta_{m}^{\prime\prime}+\frac{\delta_{m}^{\prime}}{2a}(3-\epsilon_{\rm
BD})-\frac{3}{2a^{2}}\delta_{m}\left(1+\frac{4\epsilon_{\rm
BD}}{3}\right)=0\,.$ (7.43)

If $\epsilon_{\rm BD}<0$ the Poisson term (the last in the above equation)
decreases and, on top of that, the Hubble friction increases w.r.t. the case
with $\varphi=$const. (or the GR-$\Lambda$CDM model, if we consider the same
energy densities). Both effects help to slow down the structure formation in
the Universe. Of course, if $\epsilon_{\rm BD}>0$ the opposite happens. This
is confirmed by solving explicitly Eq.(7.43). Despite an exact solution to Eq.
(7.43) can be found, it suffices to quote it at $\mathcal{O}(\epsilon_{\rm
BD})$ and neglect the $\mathcal{O}(\epsilon^{2}_{\rm BD})$ corrections. The
growing and decaying modes at leading order read $\delta^{+}_{m}(a)\sim
a^{1+\epsilon_{\rm BD}}$ and $\delta^{-}_{m}(a)\sim
a^{-\frac{1}{2}(3+\epsilon_{\rm BD})}$, respectively. The latter just fades
soon into oblivion and the former explains why negative values of
$\epsilon_{\rm BD}$ are favored by the data on RSD, since $\epsilon_{\rm
BD}<0$ obviously slows down the rate of structure formation and hence acts as
an effective (positive) contribution to the vacuum energy density171717This
feature was already noticed in the preliminary treatment of Ref. [331] for the
BD theory itself, and it was actually pointed out as a general feature of the
class of Running Vacuum Models (RVM), which helps to cure the
$\sigma_{8}$-tension [233, 229]. This is remarkable, since the RVM’s turn out
to mimic BD-gravity, as first noticed in [332]. For a summary, see Appendix
C.. The preference for negative values of $\epsilon_{\rm BD}$ is especially
clear when the RSD data include the matter bispectrum information, which tends
to accentuate the slowing down of the growth function, as noted repeatedly in
a variety of previous works [78, 251, 110, 315]. We may clearly appraise this
feature also in the present study, see e.g. Tables 3-4 (with spectrum and
bispectrum) and 5 (with spectrum but no bispectrum), where the $\sigma_{8}$
value is in general well-behaved ($\sigma_{8}\simeq 0.8$) in the
BD-$\Lambda$CDM framework when $\epsilon_{\rm BD}<0$, but it is clearly
reduced (at a level $\sigma_{8}\simeq 0.78-0.79$) in the presence of
bispectrum data. And in both cases the value of $H_{0}$ is in the range of
$70-71$ km/s/Mpc. Most models trying to explain both tensions usually increase
$\sigma_{8}$ substantially ($0.82-0.85$).

We can also study the pure vacuum-dominated epoch (VDE) in the same way. In
this case $\varphi\sim a^{2\epsilon_{\rm BD}}$ (cf. Appendix D), and hence

$a\frac{\varphi^{\prime}}{\varphi}=2\epsilon_{\rm
BD}+\mathcal{O}(\epsilon^{2}_{\rm BD})\simeq 2\epsilon_{\rm BD}\,,$ (7.44)

again with $\mathcal{F}^{\prime}=\mathcal{O}(\epsilon^{2}_{\rm BD})\simeq 0$.
The Poisson term can be neglected in this case since
$\rho_{m}\ll\rho\simeq\rho_{\Lambda}$, and hence,

$\delta_{m}^{\prime\prime}+\frac{\delta_{m}^{\prime}}{a}(3-\epsilon_{\rm
BD})=0\,.$ (7.45)

When the vacuum energy density rules the expansion of the Universe, there is a
stable constant mode solution $\delta_{m}=$const. and a decaying mode that
decreases faster than in the GR scenario if $\epsilon_{\rm BD}<0$, again due
to the fact that the friction term is in this case larger than in the standard
picture, specifically the latter reads $\delta^{-}_{m}(a)\sim
a^{-2+\epsilon_{\rm BD}}$ in the ${\cal O}(\epsilon_{\rm BD})$ approximation.

Figure 6: As in Fig. 3, but comparing now the $\Lambda$CDM baseline
configuration defined there with the case $\epsilon_{\rm BD}=-0.003$ (and
$\varphi_{ini}=1$) of the BD-$\Lambda$CDM model.

The analytical study of the transition between the matter and vacuum-dominated
epochs is more difficult, but with what we have already seen it is obvious
that the amount of structure generated also in this period of the cosmic
expansion will be lower than in the $\varphi=$const. case if $\epsilon_{\rm
BD}$ takes a negative value. In Fig. 5 we show this explicitly.

From this analysis it should be clear that if $\epsilon_{\rm BD}<0$ there is a
decrease of the matter density contrast for fixed energy densities when
compared with the GR-$\Lambda$CDM scenario, and also with the BD scenario with
$\varphi={\rm const}$. In Fig. 6 we can see this feature directly in the
matter power spectrum, which is seen to be suppressed with respect to the case
$\epsilon_{\rm BD}=0$ for those scales that are relevant for the RSD, i.e.
within the range of wave numbers $0.01h{\rm Mpc^{-1}}\lesssim k\lesssim
0.1h{\rm Mpc^{-1}}$ (corresponding to distance scales roughly between a few
dozen to a few hundred Mpc). However, we still don’t know whether these
negative values of $\epsilon_{\rm BD}$ can also be accommodated by the other
datasets.

Figure 7: As in Fig. 1, but comparing the GR-$\Lambda$CDM baseline scenario
(equivalent to the BD-$\Lambda$CDM one for $\epsilon_{\rm BD}=0$, $\varphi=1$)
with the case $\epsilon_{\rm BD}=\pm 0.003$, using as initial condition
$\varphi_{ini}=1$.

To check this, let us recall from our discussion above (see also Appendices C
and D for more details) that the evolution of the BD-field takes place
basically during the MDE. In the RDE the scalar field is essentially frozen
once the decaying mode becomes irrelevant, and although in the late-time
Universe $\varphi$ evolves faster than in the MDE (compare (7.33) and (7.44))
it remains almost constant in the redshift range $z\lesssim\mathcal{O}(1)$, in
which all the non-CMB data points lie, particularly the LSS data. Let us note
that the typical values of $\epsilon_{\rm BD}$ fitted from the overall set of
data used in our analysis (cf. Sec. 7.6) place that parameter in the ballpark
of $|\epsilon_{\rm BD}|\sim\mathcal{O}(10^{-3})$ (cf. Tables 3-5, for
example). Schematically, we can think that the field takes a value
$\varphi_{ini}$ during the RDE, then if $\epsilon_{\rm BD}<0$ it decreases an
amount $\Delta\varphi_{ini}<0$ with respect to its original value
$\varphi_{ini}$, during the MDE, and finally
$\varphi^{0}\approx\varphi_{ini}+\Delta\varphi_{ini}$, with
$\varphi(z)\approx\varphi^{0}$ at $z\lesssim\mathcal{O}(1)$. Thus, Eq. (7.15)
still applies in good approximation during the RDE and the late-time Universe,
but with two different values of the cosmological Newton’s coupling in the two
(widely separated) cosmological eras. Taking these facts into account it is
easy to understand why SNIa data are not able to tell us anything about
$\epsilon_{\rm BD}$ when used alone, since again $\varphi^{0}$ is fully
degenerated with the absolute magnitude parameter $M$ of the supernovae. This
does not mean, though, that SNIa data cannot tell us anything about the
$H_{0}$-tension when they are considered together with other datasets that do
provide constraints on $H_{0}$, since the constraints that SNIa impose on
$\Omega_{m}$ help to break degeneracies present in the other datasets and to
tighten the allowed region of parameter space. H0LICOW and CCH data, for
instance, will allow us to put constraints again on $\Omega_{m}$ and also on
$H_{0}$ (hence on $\rho^{0}/\varphi^{0}$). BAO will constrain $\rho_{b}^{0}$,
$\Omega_{m}$, and now also $\Delta\varphi_{ini}/\varphi_{ini}$, which is
proportional to $\epsilon_{\rm BD}$. For instance,
$r_{s}H(z)\propto\left(1-\frac{\Delta\varphi_{ini}}{2\varphi_{ini}}\right)$.

The BD effect caused on the temperature spectrum of CMB anisotropies is
presented in the fourfold plot in Fig. 7. Due to the fact that now we have
$\epsilon_{\rm BD}\neq 0$, $\Delta\varphi\neq 0$, a small shift in the
location of the peaks is naturally generated. In the right-bottom plot of such
figure one can see that negative values of $\epsilon_{\rm BD}$ move the peaks
slightly towards lower multipoles, and the other way around for positive
values of this parameter. It is easy to understand why. To start with, let us
remark that we have produced all the curves of this plot using the same
initial condition $\varphi_{ini}=1$ and the same $\rho^{0}_{b}$,
$\rho^{0}_{cdm}$ and $\rho_{\Lambda}$ and, hence, fixing in the same way the
complete evolution of the energy densities for the different plots in that
figure. This means that the differences in the Hubble function can only be due
to differences in the evolution of the BD scalar field. The modified expansion
histories produce changes in the value of $\theta_{*}=r_{s}/D_{A,rec}$ (with
$D_{A,rec}$ being the angular diameter distance to the last scattering
surface), so also in the location of the peaks. If $\epsilon_{\rm BD}<0$,
$\varphi$ decreases with the expansion, so its value at recombination and at
present is lower than when $\epsilon_{\rm BD}=0$, and correspondingly
$G(\varphi)$ will be higher. Because of this, the value of the Hubble function
will be larger, too, and the cosmological distances lower, so the relation
between $\theta_{*}(\epsilon_{\rm BD}\neq 0)$ and $\theta_{*}(\epsilon_{\rm
BD}=0)$ can be written as follows:

$\theta_{*}(\epsilon_{\rm BD}\neq 0)=\frac{r_{s}(\epsilon_{\rm BD}\neq
0)}{D_{A,rec}(\epsilon_{\rm BD}\neq 0)}=\frac{X\cdot r_{s}(\epsilon_{\rm
BD}=0)}{Y\cdot D_{A,rec}(\epsilon_{\rm
BD}=0)}=\frac{X}{Y}\,\theta_{*}(\epsilon_{\rm BD}=0)\,,$ (7.46)

where the rescaling factors satisfy $0<X,Y<1$ for $\epsilon_{\rm BD}<0$. As we
have already mentioned before, most of the variation of $\varphi$ occurs
during the MDE, so the largest length reduction will be in the cosmic stretch
from recombination to the present time, and thereby $Y<X$. Thus, if
$\epsilon_{\rm BD}<0$ we find $\theta_{*}(\epsilon_{\rm
BD}<0)>\theta_{*}(\epsilon_{\rm BD}=0)$ and the peaks of the TT CMB spectrum
shifts towards lower multipoles. Analogously, if $\epsilon_{\rm BD}$ is
positive $X,Y>1$, with $Y>X$, so $\theta_{*}(\epsilon_{\rm
BD}>0)<\theta_{*}(\epsilon_{\rm BD}=0)$ and the peaks move to larger
multipoles. It turns out, however, that these shifts, and also the changes in
the amplitude of the peaks, can be compensated by small changes in the baryon
and DM energy densities, as we will show in Sec. 7.7.
At this point we would like to recall why in the GR-$\Lambda$CDM concordance
model it is not possible to reconcile the local measurements of $H_{0}$ with
its CMB-inferred value. In the concordance model, which we assume spatially
flat, the current value of the cold dark matter density is basically fixed by
the amplitude of the first peak of the CMB temperature anisotropies, and
$\rho_{b}^{0}$ by the relative amplitude of the second and third peaks with
respect to the first one. As a result, even if the cosmological term plays no
role in the early Universe, one finds that in order to explain the precise
location of the CMB peaks the value of $\rho_{\Lambda}$ obtained from the
matching of the predicted $\theta_{*}$ with the measured peak positions
determines $\rho_{\Lambda}$ so precisely that it leaves little margin. This
causes a problem since such narrow range of values is not in the right range
to explain the value of the current Hubble parameter measured with the cosmic
distance ladder technique [280, 321, 84] and the Strong Lensing time delay
angular diameter distances from H0LICOW [97]. In other words, despite the
concordance model fits in a remarkably successful way the CMB and BAO, and
also the SNIa data, there is an irreducible internal discordance in the
parameter needs to explain with precision both the physics of the early and of
the late-time Universe. This is of course the very expression of the
$H_{0}$-tension, to which we have to add the $\sigma_{8}$ one.
Cosmographic analyses based on BAO and SNIa data calibrated with the
GR-$\Lambda$CDM Planck preferred value of $r_{s}$ also lead to low estimates
of $H_{0}$. This is the so-called inverse cosmic distance ladder approach,
adopted for instance in [157, 91, 283, 285]. This has motivated cosmologists
to look for alternative theoretical scenarios (for instance, the generic class
of EDE proposals) able to increase the expansion rate of the Universe before
the decoupling of the CMB photons and, hence, to lower $r_{s}$ down. This, in
principle, demands an increase of the Hubble function at present in order not
to spoil the good fit to the BAO and CMB observable. Nevertheless, not all the
models passing the BAO and CMB constraints and predicting a larger value of
$H_{0}$ satisfy the ‘golden rule’ mentioned in the Introduction, since they
can lead e.g. to a worsening of the $\sigma_{8}$-tension. As an example, we
can mention some early DE models, e.g. those discussed in [333]. In these
scenarios there is a very relevant DE component which accounts for the $\sim
7\%$ of the total matter-energy content of the Universe at redshifts $\sim
3000-5000$, before recombination. This allows of course to enhance the
expansion rate and reduce $r_{s}$. After such epoch, the DE decays into
radiation. In order not to alter the position of the CMB peaks and BAO
relative distances, an increase of the DM energy density is needed. According
to [334], this leads to an excess of density power and an increase of
$\sigma_{8}$ which is not welcome by LSS measurements, including RSD, Weak
Lensing and galaxy clustering data. Another example is the interesting
modified gravity model analyzed in [335], based on changing the cosmological
value of $G$ also in the pre-recombination era, thus mimicking an increase of
the effective number of relativistic degrees of freedom in such epoch. The
additional component gets eventually diluted at a rate faster than radiation
in the MDE and it is not clear if an effect is left at present181818In stark
contrast to the model of [335], in BD-$\Lambda$CDM cosmology the behaviour of
the effective $\rho_{\rm BD}$ (acting as a kind of additional DE component
during the late Universe) mimics pressureless matter during the MDE epoch and
modifies the effective EoS of the DE at present, see the next Sec. 7.3 for
details.. This model also fits the CMB and BAO data well and loosens at some
extent the $H_{0}$-tension, but violates the golden rule of the tension
solver, as it spoils the structure formation owing to the very large values of
$\sigma_{8}\sim 0.84-0.85$ that are predicted (see the discussion in Sec. 7.8
for more details).

Our study shows that a value of the cosmological gravitational coupling about
$\sim 10\%$ larger than $G_{N}$ can ameliorate in a significant way the
$H_{0}$-tension, while keeping the values of all the current energy densities
very similar to those found in the GR-$\Lambda$CDM model. If, apart from that,
we also allow for a very slow running (increase) of the cosmological $G$
triggered by negative values of order $\epsilon_{\rm
BD}\sim-\mathcal{O}(10^{-3})$, we can mitigate at the same time the
$\sigma_{8}$-tension when only the CMB TT+lowE anisotropies are considered.
When the CMB polarizations and lensing are also included in the analysis, then
$\sigma_{8}$ is kept at the GR-$\Lambda$CDM levels, and the sign of
$\epsilon_{\rm BD}$ is not conclusive. In all situations we can preserve the
golden rule. We discuss in detail the numerical results of our analysis in
Sec. 7.7.

###  Effective equation of state of the dark energy in the BD-$\Lambda$CDM
model

Our aim in this section is to write down the Brans-Dicke cosmological
equations (7.3)-(7.5) in the context of what we may call the “effective GR-
picture”. This means to rewrite them in such a way that they can be thought of
as an effective model within the frame of GR, thus providing a parametrized
departure from GR at the background level. We will see that the main outcome
of this task, at least qualitatively, is that the BD-$\Lambda$CDM model
(despite it having a constant vacuum energy density $\rho_{\Lambda}$) appears
as one in the GR class, but with a dynamical DE rather than a CC. The dynamics
of such an effective form of DE is a function of the BD-field $\varphi$. We
wish to compute its effective EoS. In order to proceed, the first step is to
rewrite Eq.(7.3) à la Friedmann:

$3H^{2}=8\pi{G_{N}}(\rho+\rho_{\varphi})\,,$ (7.47)

where $\rho$ is the total energy density as defined previously (coincident
with that of the GR-$\Lambda$CDM model), and $\rho_{\varphi}$ is the
additional ingredient that is needed, which reads

$\rho_{\varphi}\equiv\frac{3}{8\pi{G_{N}}}\left(H^{2}\Delta{\varphi}-H\dot{\varphi}+\frac{\omega_{\rm
BD}}{6}\frac{\dot{\varphi}^{2}}{\varphi}\right).$ (7.48)

Remember the definition $\varphi(t)\equiv G_{N}\psi(t)$ made in (7.1), and we
have now introduced

$\Delta{\varphi}(t)\equiv 1-\varphi(t)\,,$ (7.49)

which tracks the small departure of $\varphi$ from one and hence of
$G(\varphi)$ from $G_{N}$ (cf. Sec. 7.1). Note that $\varphi=\varphi(t)$
evolves in general with the expansion, but very slowly since $\epsilon_{\rm
BD}$ is presumably fairly small.
From the above Eq. (7.48) it is pretty clear that we have absorbed all the
terms beyond the $\Lambda$CDM model into the expression of $\rho_{\varphi}$.
While it is true that we define this quantity as if it were an energy density,
it is important to bear in mind that it is not associated to any kind of
particle, it is just a way to encapsulate those terms that are not present in
the standard model. This quantity, however, satisfies a local conservation law
as if it were a real energy density, as we shall see in a moment. From the
generalized Friedmann equation (7.47) and the explicit expression for
$\rho_{\varphi}$ given above we can write down the generalized cosmic sum rule
verified by the BD-$\Lambda$CDM model in the effective GR-picture:

$\Omega_{m}+\Omega_{r}+\Omega_{\Lambda}+\Omega_{\varphi}=1\,,$ (7.50)

where the $\Omega_{i}$ are the usual (current) cosmological parameters of the
concordance $\Lambda$CDM, whereas $\Omega_{\varphi}$ is the additional one
that parametrizes the departure of the BD-$\Lambda$CDM model from the
GR-$\Lambda$CDM in the context of the GR-picture, and reads

$\Omega_{\varphi}=\frac{\rho_{\varphi}^{0}}{\rho^{0}_{c}}.$ (7.51)

Notice that the above sum rule is exact and it is different from that in Eq.
(7.19) since the latter is only approximate for the case when $\epsilon_{\rm
BD}=0$ or very small. These are two different pictures of the same
BD-$\Lambda$CDM model. The modified cosmological parameters (7.18) depend on
$\varphi$ whereas here the $\varphi$-dependence has been fully concentrated on
$\Omega_{\varphi}$. It is interesting to write down the exact equation (7.50)
in the form

$\Omega_{m}+\Omega_{r}+\Omega_{\Lambda}=1-\Omega_{\varphi}=1-\Delta\varphi_{0}+\frac{\dot{\varphi}_{0}}{H_{0}}-\frac{\omega_{\rm
BD}}{6}\frac{\dot{\varphi}_{0}^{2}}{H_{0}^{2}\varphi_{0}}\,.$ (7.52)

in which $\varphi_{0}=\varphi(z=0)$ and
$\dot{\varphi}_{0}=\dot{\varphi}(z=0)$. For $\epsilon_{\rm BD}\simeq 0$ we
know that $\varphi\simeq$const. and we can neglect the time derivative terms
and then we find the approximate form
$\Omega_{m}+\Omega_{r}+\Omega_{\Lambda}=1-\Omega_{\varphi}\simeq
1-\Delta\varphi_{0}$. This equation suggests that a value of
$\Delta\varphi_{0}\neq 0$ would emulate the presence of a small fictitious
spatial curvature in the GR-picture. See e.g.[50, 54, 53] and references
therein for the study of a variety models explicitly involving spatial
curvature.

The second step in the process of constructing the GR-picture of the BD theory
is to express (7.4) as in the usual pressure equation for GR, and this forces
us to define a new pressure quantity $p_{\varphi}$ associated to
$\rho_{\varphi}$. We find

$2\dot{H}+3H^{2}=-8\pi{G_{N}}(p+p_{\varphi}),$ (7.53)

with

$p_{\varphi}\equiv\frac{1}{8\pi{G_{N}}}\left(-3H^{2}\Delta{\varphi}-2\dot{H}\Delta{\varphi}+\ddot{\varphi}+2H\dot{\varphi}+\frac{\omega_{\rm
BD}}{2}\frac{\dot{\varphi}^{2}}{\varphi}\right).$ (7.54)

On the face of the above definitions (7.48) and (7.54), we can now interpret
the BD theory as an effective theory within the frame of General Relativity,
which deviates from it an amount $\Delta{\varphi}$. Indeed, for
$\Delta{\varphi}=0$ we have $\rho_{\varphi}=p_{\varphi}=0$ and we recover GR.
Mind that $\Delta{\varphi}=0$ means not only that $\varphi=$const. (hence that
$\epsilon_{\rm BD}=0$, equivalently $\omega_{\rm BD}\to\infty$), but also that
that constant is exactly $\varphi=1$. In such case $G_{\rm eff}$ is also
constant and $G_{\rm eff}=G_{N}$ exactly. The price that we have to pay for
such a GR-like description of the BD model is the appearance of the fictitious
BD-fluid with energy density $\rho_{\varphi}$ and pressure $p_{\varphi}$,
which complies with the following conservation equation throughout the
expansion of the Universe 191919The new ‘fluid’ that one has to add to GR to
effectively mimic BD plays a momentous role to explain the $H_{0}$ and
$\sigma_{8}$-tensions. In a way it mimics the effect of the ‘early DE’ models
mentioned in the previous section, except that the BD-fluid persists for the
entire cosmic history and is instrumental both in the early as well as in the
current Universe so as to preserve the golden rule of the tension solver:
namely, it either softens the two tensions of GR or improves one of them
without detriment of the other. :

$\dot{\rho}_{\varphi}+3H(\rho_{\varphi}+p_{\varphi})=0.$ (7.55)

One can check that this equation holds after a straightforward calculation,
which makes use of the three BD-field equations (7.3)-(7.5). Although at first
sight the above conservation equation can be surprising actually it is not,
since it is a direct consequence of the Bianchi identity. Let us now assume
that the effective BD-fluid can be described by an equation of state like
$p_{\varphi}=w_{\varphi}\rho_{\varphi}$, so

$w_{\varphi}=\frac{p_{\varphi}}{\rho_{\varphi}}=\frac{-3H^{2}\Delta{\varphi}-2\dot{H}\Delta{\varphi}+\ddot{\varphi}+2H\dot{\varphi}+\frac{\omega_{\rm
BD}}{2}\frac{\dot{\varphi}^{2}}{\varphi}}{3H^{2}\Delta{\varphi}-3H\dot{\varphi}+\frac{\omega_{\rm
BD}}{2}\frac{\dot{\varphi}^{2}}{\varphi}}.$ (7.56)

The contribution from those terms containing derivatives of the BD-field are
subdominant for the whole cosmic history. We have verified this fact
numerically, see Fig. 4. While the variation of $\varphi$ between the two
opposite ends of the cosmic history is of $\sim 1.5\%$ and is significant for
our analysis, the instantaneous variation is actually negligible. Thus,
$H\dot{\varphi}$ and $\ddot{\varphi}$ are both much smaller than
$\dot{H}\Delta{\varphi}$, and in this limit we can approximate (7.56) very
accurately as

$w_{\varphi}(t)\simeq-1-\frac{2}{3}\frac{\dot{H}}{H^{2}}\,\ \ \ \ \ \ ({\rm
for}\ H\dot{\varphi},\ddot{\varphi}\ll\dot{H}\Delta{\varphi})\,.$ (7.57)

This EoS turns out to be the standard total EoS of the $\Lambda$CDM, which
boils down to the EoS corresponding to the different epochs of the cosmic
evolution (i.e. $w=1/3,0,-1$ for RDE, MDE and VDE). This means that the EoS of
the BD-fluid mimics these epochs. We can go a step further and define not just
the BD-fluid but the combined system of the BD-fluid and the vacuum energy
density represented by the density $\rho_{\Lambda}$ associated to the
cosmological constant. We define the following effective EoS for such combined
fluid:

Figure 8: Left plot: Effective equation of state of the DE in the
BD-$\Lambda$CDM model as a function of the redshift. The inner plot magnifies
the region around our time. We can see that the BD model mimics quintessence
with a significance of more than $3\sigma$; Right plot: It shows once more
Fig. 4 in order to ease the comparison of the EoS evolution, which is
associated to the evolution of the BD-field $\varphi$ – cf. Eq.(7.58). The
shadowed bands in these plots correspond to the 1$\sigma$ regions.

$w_{\rm
eff}\equiv\frac{p_{\Lambda}+p_{\varphi}}{\rho_{\Lambda}+\rho_{\varphi}}=-1+\frac{p_{\varphi}+\rho_{\varphi}}{\rho_{\Lambda}+\rho_{\varphi}}=-1+\frac{-2\dot{H}\Delta{\varphi}+f_{1}(\varphi,\dot{\varphi},\ddot{\varphi})}{\Lambda+3H^{2}\,\Delta{\varphi}+f_{2}(\varphi,\dot{\varphi})}\,,$
(7.58)

where the two functions

$f_{1}(\varphi,\dot{\varphi},\ddot{\varphi})=\ddot{\varphi}-H\dot{\varphi}+\omega_{\rm
BD}\frac{\dot{\varphi}^{2}}{\varphi}\,,\ \ \ \ \ \ \
f_{2}(\varphi,\dot{\varphi})=-3H\dot{\varphi}+\frac{\omega_{\rm
BD}}{2}\frac{\dot{\varphi}^{2}}{\varphi}$ (7.59)

involve differentiations with respect to the slowly varying field $\varphi$
and as before they are negligible, in absolute value, as compared to
$\dot{H}\Delta\varphi$ and ${H^{2}}\Delta\varphi$. The effective EoS (7.58) is
a time-evolving quantity which mimics dynamical DE at low redshifts. At very
high redshifts $z\gg 1$, well beyond the DE dominated epoch, we can neglect
$\Lambda$ in the denominator of the EoS and the dominant term is
$3H^{2}\,\Delta{\varphi}$. Similarly, in the numerator the dominant term is
always $-2\dot{H}\Delta{\varphi}$. Therefore, at high redshifts the effective
EoS (7.58) behaves as (7.57) since $\Delta{\varphi}$ cancels out: $w_{\rm
eff}(z)\simeq w_{\varphi}(z)\,\ (z\gg 1)$, which means that it just reproduces
the standard EoS of the GR-$\Lambda$CDM. The exact EoS (7.58), however, must
be computed numerically, and it is displayed in Fig. 8, together with the
numerical plot of $\varphi$ (which we have already shown in Fig. 4). We have
used the Baseline+$H_{0}$ dataset defined in Sec. 7.6. For a semi-qualitative
discussion of the combined EoS it will suffice an analytical approximation, as
we did before with $w_{\varphi}$. The most relevant part of $w_{\rm eff}(z)$
as to the possibility of disentangling the dynamical DE effects triggered by
the underlying BD model is near the present time ($z<1$). Thus, neglecting the
contribution from the functions $f_{1,2}$, but now keeping the $\Lambda$-term
in the denominator of (7.58) we can use the Hubble function of the concordance
model and we find the following result at linear order in $\Delta{\varphi}$:

$w_{\rm
eff}(z)\simeq-1-\frac{2\dot{H}\Delta{\varphi}}{\Lambda}\simeq-1+\Delta{\varphi}\,\frac{\Omega^{0}_{m}}{\Omega^{0}_{\Lambda}}\,(1+z)^{3}\,,$
(7.60)

where $\Omega^{0}_{m}$ and $\Omega^{0}_{\Lambda}$ are the current values of
the cosmological parameters, which satisfy
$\Omega^{0}_{m}+\Omega^{0}_{\Lambda}=1$ for spatially flat Universe. As has
been stated before, the previous approximate formula is valid only for $z<1$,
but it shows very clearly that for $\Delta{\varphi}>0$ (resp. $<0$) we meet
quintessence-like (resp. phantom-like) behaviour. As we have repeatedly
emphasized, our analysis points to $\epsilon_{\rm BD}<0$ and hence $\varphi$
decreases with the expansion, remaining smaller than one. From Eq. (7.49) this
means $\Delta{\varphi}>0$ and therefore we find that the effective GR
behaviour of the BD-$\Lambda$CDM is quintessence-like. We can be more precise
at this point. We have numerically computed the value of the exact function
(7.58) at $z=0$, taking into account the contribution from all the terms, in
particular the slowly varying functions (7.59), see Tables 3-6. The results
obtained from three of the most prominent datasets defined in Sec. 7.6 read as
follows:

$\displaystyle{\bf Baseline}:\quad$ $\displaystyle w_{\rm eff}(0)=$
$\displaystyle-0.983^{+0.015}_{-0.014}$ (7.61) $\displaystyle{\bf
Baseline+H_{0}}:\quad$ $\displaystyle w_{\rm eff}(0)=$
$\displaystyle-0.966^{+0.012}_{-0.011}$ (7.62) $\displaystyle{\bf
Baseline+H_{0}+SL}:\quad$ $\displaystyle w_{\rm eff}(0)=$
$\displaystyle-0.962\pm 0.011.$ (7.63)

As can be seen, there is a non-negligible departure from the constant EoS
value $-1$ of the GR-$\Lambda$CDM, which reaches the $\sim 3\sigma$ c.l. when
the prior on $H_{0}$ from the local distance ladder measurement by SH0ES [84]
is included in the analysis, and $\sim 3.5\sigma$ c.l. when also the angular
diameter distances from H0LICOW [97] are taken into account. The effective
quintessence EoS $w_{\rm eff}(0)>-1$ is one of the ingredients that allows the
BD-$\Lambda$CDM model to significantly loosen the $H_{0}$-tension, since it is
a direct consequence of having $\varphi<1$ (or, equivalently, $G>G_{N}$) (cf.
Sec. 7.2 for details).

We have obtained the above results from the equations of motion once a metric
was assumed; however, it is possible to obtain all the expressions listed in
this section starting from the BD action (1.147) itself and then considering
the FLRW metric. To show this, let us use the dimensionless field
$\varphi=G_{N}\psi$ and the variable $\Delta{\varphi}$ defined in (7.49).
First of all we split the whole action in three pieces

$S_{\rm BD}[\varphi]=S_{\rm EH}+S_{\rm GR}[\varphi]+S_{m},$ (7.64)

where

$S_{\rm EH}\equiv\int d^{4}x\sqrt{-g}\left[\frac{R}{16\pi
G_{N}}-\rho_{\Lambda}\right],$ (7.65)

is the usual Einstein-Hilbert action, whereas

$S_{\rm GR}[\varphi]\equiv\int d^{4}x\sqrt{-g}\frac{1}{16\pi
G_{N}}\left[-R\Delta{\varphi}-\frac{\omega_{\rm
BD}}{\varphi}g^{\mu\nu}\partial_{\nu}\varphi\partial_{\mu}\varphi\right],$
(7.66)

is the action parametrizing the deviation of the BD theory from the GR-picture
expressed in terms of the scalar field $\varphi$. As expected, for $\varphi=1$
that action vanishes identically. Finally,

$S_{m}\equiv\int d^{4}x\sqrt{-g}\mathcal{L}_{m}(\chi_{i},g_{\mu\nu})$ (7.67)

is the action for the matter fields. Since there is no interaction involving
$\varphi$ with other components, the BD-field $\varphi$ is covariantly
conserved, as remarked in Sec. 7.1. In order to compute the energy-momentum
tensor and find out the effective density and pressure of the BD-field, we
apply the usual definition of that tensor in curved space-time:

$T_{\mu\nu}^{\rm BD}=-\frac{2}{\sqrt{g}}\frac{\delta S_{\rm
GR}[\varphi]}{\delta g^{\mu\nu}}\,.$ (7.68)

After some calculations we arrive at

$\begin{split}T_{\mu\nu}^{\rm BD}=\frac{R_{\mu\nu}}{8\pi
G_{N}}\Delta{\varphi}-\frac{\nabla_{\nu}\nabla_{\mu}\varphi}{8\pi
G_{N}}&+\frac{g_{\mu\nu}\Box\varphi}{8\pi G_{N}}+\frac{\omega_{\rm
BD}}{8\pi\varphi}\partial_{\nu}\varphi\partial_{\mu}\varphi\\\
&-\frac{g_{\mu\nu}}{16\pi G_{N}}\left(R\Delta{\varphi}+\frac{\omega_{\rm
BD}}{\varphi}g^{\alpha\beta}\partial_{\alpha}\varphi\partial_{\beta}\varphi\right).\end{split}$
(7.69)

Since $\varphi$ has no interactions it behaves as any free scalar field, so
its energy-momentum tensor must adopt the perfect fluid form at the background
level:

$T_{\mu\nu}^{\rm
BD}=p_{\varphi}g_{\mu\nu}+(\rho_{\varphi}+p_{\varphi})u_{\mu}u_{\nu}.$ (7.70)

Now we can compare this form with (7.69). It is straightforward to obtain the
energy density as well as the corresponding pressure, we only need to compute
$\rho_{\varphi}=T_{00}^{\rm BD}$ and $p_{\varphi}=(T^{\rm
BD}+\rho_{\varphi})/3$, being $T^{\rm BD}=g^{\mu\nu}T^{\rm BD}_{\mu\nu}$ the
trace of the tensor. Using at this point the spatially flat FLRW metric one
can work out the explicit result for $\rho_{\varphi}$ and $\rho_{\varphi}$ and
reconfirm that it acquires the form previously indicated in the equations
(7.48) and (7.54). This provides perhaps a more formal derivation of these
formulas and serves as a cross-check of them.

###  Structure formation in the linear regime. Perturbations equations

In order to perform a complete analysis of the model, we need to study the
evolution of the perturbed cosmological quantities in the context of BD
theory. For a review of the standard model perturbations equations, see e.g.
[232, 237, 336]. We assume a FLRW metric written in conformal time, denoted by
$\eta$, in which the line element is
$ds^{2}=a^{2}(\eta)[-d\eta^{2}+(\delta_{ij}+h_{ij})dx^{i}dx^{j}]$. Here
$h_{ij}$ is a perturbation on the spatial part of the metric which can be
expressed in momentum space as follows,

$h_{ij}(\eta,\vec{x})=\int
d^{3}k\,e^{-i\vec{k}\cdot\vec{x}}\left[\hat{k}_{i}\hat{k}_{j}h(\eta,\vec{k})+\left(\hat{k}_{i}\hat{k}_{j}-\frac{\delta_{ij}}{3}\right)6\xi(\eta,\vec{k})\right].$
(7.71)

As we see, in momentum space the trace $h\equiv\delta^{ij}h_{ij}$ decouples
from the traceless part of the perturbation, $\xi$. Now, we are going to list
the perturbations equations at late stages of expansion in momentum space at
deep subhorizon scales, that is, we assume $\mathcal{H}^{2}\ll k^{2}$, with
$\mathcal{H}\equiv a^{\prime}/a$. Although primes were used previously for
derivatives with respect to the scale factor, they will henceforth stand for
derivatives with respect to the conformal time: $()^{\prime}\equiv d()/d\eta$.
For example, it is easy to see that $\mathcal{H}=aH$. One may work with the
standard differential equation for the density contrast at deep subhorizon
scales,

$\delta_{m}^{\prime\prime}+\mathcal{H}\delta_{m}^{\prime}-4\pi{G_{N}}\bar{\rho}_{m}a^{2}\delta_{m}=0,$
(7.72)

where $\delta_{m}\equiv\delta\rho_{m}/\bar{\rho}_{m}$ is the density contrast,
the bar over $\bar{\rho}_{m}$ indicates that is a background quantity and the
evolution of $\mathcal{H}$ and $\bar{\rho}_{m}$ is the one expected by the
background equations of the BD theory in Section 7.1. This expression is just
the corresponding one for the $\Lambda$CDM, completely neglecting any possible
perturbation in the BD-field, namely $\delta\varphi$. However, it is possible
to see that a second order differential equation for the density contrast can
be written, even if the perturbation in $\varphi$ is not neglected. This is
done in detail in the Appendix B. In this section, we present the main
perturbations equations in the case of the Synchronous gauge and discuss the
interpretation of the result.

If $\vec{v}_{m}$ is the physical 3-velocity of matter (which is much smaller
than 1 and can be treated as a perturbation), then we can define its
divergence, $\theta_{m}\equiv\vec{\nabla}\cdot(\vec{v}_{m})$. At deep
subhorizon scales it is possible to see that the equation governing its
evolution is

$\theta_{m}^{\prime}+\mathcal{H}\theta_{m}=0.$ (7.73)

Since $da^{-1}/d\eta=-\mathcal{H}/a$, we arrive to a decaying solution
$\theta_{m}\propto a^{-1}$. A common assumption is to set $\theta_{m}\sim 0$
in the last stages of the Universe, which is what we will do in our analysis.
This allows us to simplify the equations. Another simplification occurs if we
take into account that we are basically interested in computing the matter
perturbations only at deep subhorizon scales, namely for
$k^{2}\gg\mathcal{H}^{2}$, which allows us to neglect some terms as well (cf.
Appendix B). Altogether we are led to the following set of perturbations
equations in the synchronous gauge:

$\delta_{m}^{\prime}=-\frac{h^{\prime}}{2}\,.$ (7.74)
$k^{2}\delta\varphi+\frac{h^{\prime}}{2}\bar{\varphi}^{\prime}=\frac{8\pi
G_{N}}{3+2\omega_{\rm BD}}a^{2}\bar{\rho}_{m}\delta_{m}\,,$ (7.75)
$\bar{\varphi}(\mathcal{H}h^{\prime}-2\xi
k^{2})+k^{2}\delta\varphi+\frac{h^{\prime}}{2}\bar{\varphi}^{\prime}=8\pi
G_{N}a^{2}\bar{\rho}_{m}\delta_{m}\,,$ (7.76)
$2k^{2}\delta\varphi+\bar{\varphi}^{\prime}h^{\prime}+\bar{\varphi}\left(h^{\prime\prime}+2h^{\prime}\mathcal{H}-2k^{2}\xi\right)=0\,.$
(7.77)

Combining these four equations simultaneously and without doing any further
approximation, one finally obtains the following compact equation for the
matter density contrast of the BD theory at deep subhorizon scales:

$\delta_{m}^{\prime\prime}+\mathcal{H}\delta_{m}^{\prime}-\frac{4\pi
G_{N}a^{2}}{\bar{\varphi}}\bar{\rho}_{m}\delta_{m}\left(\frac{4+2\omega_{\rm
BD}}{3+2\omega_{\rm BD}}\right)=0\,.$ (7.78)

In other words,

$\delta_{m}^{\prime\prime}+\mathcal{H}\delta_{m}^{\prime}-4\pi G_{\rm
eff}(\bar{\varphi})a^{2}\,\bar{\rho}_{m}\delta_{m}=0\,.$ (7.79)

The quantity

$G_{\rm
eff}(\bar{\varphi})=\frac{G_{N}}{\bar{\varphi}}\left(\frac{4+2\omega_{\rm
BD}}{3+2\omega_{\rm
BD}}\right)=\frac{G_{N}}{\bar{\varphi}}\left(\frac{2+4\epsilon_{\rm
BD}}{2+3\epsilon_{\rm BD}}\right)$ (7.80)

is precisely the effective coupling previously introduced in Eq. (7.24); it
modifies the Poisson term of the perturbations equation with respect to that
of the standard model, Eq. (7.72). There is, in addition, a modification in
the friction term between the two models, which is of course associated to the
change in ${\cal H}$.

The argument of $G_{\rm eff}$ in (7.80) is not $\varphi$ but the background
value $\bar{\varphi}$ since the latter is what remains in first order of
perturbations from the consistent splitting of the field into the background
value and the perturbation: $\varphi=\bar{\varphi}+\delta\varphi$. Notice that
there is no dependence left on the perturbation $\delta\varphi$. As we can see
from (7.80), the very same effective coupling that rules the attraction of two
tests masses in BD-gravity is the coupling strength that governs the formation
of structure in this theory, as it could be expected. But this does not
necessarily mean that the effective gravitational strength governing the
process of structure formation is the same as for two tests masses on Earth.
We shall elaborate further on this point in the next section. At the moment we
note that if we compare the above perturbations equation with the standard
model one (7.72), the former reduce to the latter in the limit $\omega_{\rm
BD}\to\infty$ (i.e. $\epsilon_{\rm BD}\to 0$) and $\bar{\varphi}=1$.

The form of (7.79) in terms of the scale factor variable rather than in
conformal time was given previously in Sec. 7.2.1 when we considered a preview
of the implications of BD-gravity on structure formation data202020Recall,
however, that prime in Eq. (7.79) stands for differentiation w.r.t. conformal
time whereas in Eq. (7.23) denotes differentiation w.r.t. the scale factor.
These equations perfectly agree and represent the same perturbations equation
for the matter density field in BD-gravity in the respective variables. They
are also in accordance with the perturbations equation obeyed by the matter
density field within the context of scalar-tensor theories with the general
action (7.9)[316] of which the form (7.8) is a particular case. . The
transformation of derivatives between the two variables can be easily
performed with the help of the chain rule $d/d\eta=a{\cal H}d/da$.

###  Different BD scenarios and Mach’s Principle

As previously indicated, the relation (7.24), which appears now in the
cosmological context in the manner (7.80), follows from the computation of the
gravitational field felt by two test point-like (or spherical) masses in
interaction in BD-gravity within the weak-field limit[124, 291], see also
[300] and references therein. Such relation shows in a manifest way the
integration of Mach’s principle within the BD context, as it postulates a link
between the measured local value of the gravitational strength, $G_{N}$, as
measured at the Earth surface, and its cosmological value, $G_{\rm
eff}(\varphi)$, which depends on $\varphi$ and $\omega_{\rm BD}$. In
particular, $\varphi$ may be sensitive to the mean energy densities and
pressures of all the matter and energy fields that constitute the Universe. If
there is no mechanism screening the BD-field on Earth, $G_{\rm
eff}(\bar{\varphi})(z=0)=G_{N}$. However, one can still fulfill this condition
if Eq. (7.80) constraints the current value of the cosmological BD-field
$\bar{\varphi}$

$\bar{\varphi}(z=0)=\frac{4+2\omega_{\rm BD}}{3+2\omega_{\rm
BD}}=\frac{2+4\epsilon_{\rm BD}}{2+3\epsilon_{\rm BD}}\simeq
1+\frac{1}{2}\,\epsilon_{\rm BD}+{\cal O}(\epsilon_{\rm BD}^{2})\,.$ (7.81)

That is to say, such constraint permits to reconcile $G_{\rm
eff}(\bar{\varphi})(z=0)$ with $G_{N}$ by still keeping
$\bar{\varphi}(z=0)\neq 1$ and $\epsilon_{\rm BD}\neq 0$. Hence the BD-field
can be dynamical and there can be a departure of $G(\bar{\varphi})$ from
$G_{N}$ even at present. This constraint, however, is much weaker than the one
following from taking the more radical approach in which $G_{\rm
eff}(\varphi)$ and $G_{N}$ are enforced to coincide upon imposing the double
condition $\epsilon_{\rm BD}\to 0$ (i.e. $\omega_{\rm BD}\to\infty$) and
$\bar{\varphi}=1$. It is this last setup which anchors the BD theory to remain
exactly (or very approximately) close to GR at all scales.

However, if we take seriously the stringent constraint imposed by the Cassini
probe on the post-Newtonian parameter $\gamma^{\rm PN}$[317], which leads to a
very large value of $\omega_{\rm BD}\gtrsim 10^{4}$ (equivalently, a very
small value of $\epsilon_{\rm BD}=1/\omega_{\rm BD}$), as we discussed in Sec.
7.1, $\bar{\varphi}$ must remain almost constant throughout the cosmic
expansion, thus essentially equal to $\bar{\varphi}(z=0)$. However, the
Cassini limit leaves $\bar{\varphi}(z=0)$ unconstrained, so this constant is
not restricted to be in principle equal to $1$. In this case the relation
(7.81) may or may not apply; there is in fact no especial reason for it (it
will depend on the effectiveness of the screening mechanism on Earth). If it
does, i.e. if there is no screening, $G_{\rm eff}$ is forced to be very close
to $G_{N}$ $\forall{z}$; if it does not, $\bar{\varphi}$ can freely take
(almost constant) values which do not push $G_{\rm eff}$ to stay so glued to
$G_{N}$. It is interesting to see the extent to which the cosmological
constraints can compete with the local ones given the current status of
precision they can both attain.
So, as it turns out, we find that one of the two interpretations leads to
values of $G_{\rm eff}$ very close to $G_{N}$ $\forall{z}$ on account of the
fact that we are imposing very large values of $\omega_{\rm BD}$ and assuming
(7.81), whereas the other achieves the same aim (viz. $G_{\rm eff}$ can stay
very close to $G_{N}$) for intermediate values of $\omega_{\rm BD}$ provided
they are linked to $\bar{\varphi}(z=0)$ through the constraint (7.81). This
last option, as indicated, is not likely since this would imply the existence
of a screening at the scales probed by Cassini that may become ineffective on
Earth, where the densities are higher. Finally, we may as well have a
situation where the cosmological $G_{\rm eff}$ remains different from $G_{N}$
even if the Cassini limit is enforced. For this to occur we need an
(essentially constant) value of $\bar{\varphi}\neq 1$ (different from the one
associated to that constraint) at the cosmological level. This can still be
compatible with the local constraints provided $\varphi$ is screened on Earth
at $z=0$.
A more open-minded and general approach, which we are going to study in this
chapter, is to take the last mentioned option but without the Cassini limit.
This means that $\epsilon_{\rm BD}$ is not forced to be so small and hence
$\varphi$ can still have some appreciable dynamics. We assume that the pure BD
model applies from the very large cosmological domains to those at which the
matter perturbations remain linear. Equation (7.80) predicts the cosmological
value of $G_{\rm eff}$ once $\omega_{\rm BD}$ and the initial value
$\varphi_{ini}$ are fitted to the data. We can dispense with the Cassini
constraint (which affects $\omega_{\rm BD}$ only) because we assume that some
kind of screening mechanism acts at very low (astrophysical) scales, namely in
the non-linear domain, without altering the pure BD model at the cosmological
level. To construct a concrete screening mechanism would imply to specify some
microscopic interaction properties of $\varphi$ with matter, but these do not
affect the analysis at the cosmological level, where there are no place with
high densities of material particles. But once such mechanism is constructed
(even if not being the primary focus of our work) the value of the BD-field
$\varphi$ is ambient dependent, so to speak, since $\varphi$ becomes sensitive
to the presence of large densities of matter. This possibility is well-known
in the literature through the chameleon mechanism [328] and in the case of the
BD-field was previously considered in [318] without letting the Brans-Dicke
parameter $\omega_{\rm BD}$ to acquire negative values, and using datasets
which now can be considered a bit obsolete. Here we do, instead, allow
negative values for $\omega_{\rm BD}$ (we have seen in Sec. 7.2 the
considerable advantages involved in this possibility), and moreover we are
using a much more complete set of observations from all panels of data taking.
In this scenario we cannot make use of (7.80) to connect the locally measured
value of the gravitational strength $G_{N}$ with the BD-field at cosmological
scales. We just do not need to know how the theory exactly connects these two
values. We reiterate once more: we will not focus on the screening mechanism
itself here but rather on the properties of the BD-field in the Universe in
the large, i.e. at the cosmological level. As it is explained in [318] – see
also [294] – many scalar-tensor theories of gravity belonging to the Horndeski
class [299] could lead to such kind of BD behaviour at cosmological level and,
hence, it deserves a dedicated and updated analysis, which is currently
lacking in the literature.

To summarize, the following interpretations of the BD-gravity framework
considered here are, in principle, possible in the light of the current
observational data:

  * •

BD-Scenario I: Rigid Scenario for both the Local and Cosmological domains. In
it, we have $G_{\rm eff}(\varphi(z))\simeq G(\varphi(z))\simeq G_{N}$, these
three couplings being so close that in practice BD-gravity is
indistinguishable from GR. In this case, the BD-gravity framework is assumed
to hold on equal footing with all the scales of the Universe, local and
cosmological. There are no screening effects from matter. In this context, one
interprets that the limit from the Cassini probe [317] leads to a very large
value of $\omega_{\rm BD}$, which enforces $\varphi$ to become essentially
rigid, and one assumes that such constant value is very close to $1$ owing to
the relation (7.81). Such a rigid scenario is, however, unwarranted. It is
possible, although is not necessary since, strictly speaking, there is no
direct connection between the Cassini bound on $\omega_{\rm BD}$ and the value
that $\varphi$ can take. Thus, in this scenario the relation (7.81) is just
assumed. In point of fact, bounds on $\omega_{\rm BD}$ can only affect the
time evolution of $\varphi$, they do not constraint its value.

  * •

BD-Scenario II (Main): Locally Constrained but Cosmologically Unconstrained
Scenario. It is our main scenario. It assumes a constrained situation in the
local domain, caused by the presence of chameleonic forces, but permits an
unconstrained picture for the entire cosmological range. In other words, the
Cassini limit that holds for the post-Newtonian parameter $\gamma^{\rm
PN}\,$in the local astrophysical level (and hence on $\omega_{\rm BD}$) is
assumed to reflect just the presence of screening effects of matter in that
non-linear domain. These effects are acting on $\varphi$ and produce the
illusion that $\omega_{\rm BD}$ has a very large value (as if ‘dressed’ or
‘renormalized’ by the chamaleonic forces). One expects that the ‘intrusive’
effects of matter are only possible in high density (hence non-relativistic)
environments, and in their presence we cannot actually know the real (‘naked’)
value of $\omega_{\rm BD}$ through local experiments alone. We assume that the
screening disappears as soon as we leave the astrophysical scales and plunge
into the cosmological ones; then, and only then, we can measure the naked or
“bare’ value of $\omega_{\rm BD}$ (stripped from such effects). We may assume
that the screening ceases already at the LSS scales where linear structure
formation occurs, see e.g. [326] for examples of potentials which can help to
realize this mechanism. The bare value of $\omega_{\rm BD}$ can then be fitted
to the overall data, and in particular to the LSS formation data. Since
$\omega_{\rm BD}$ does no longer appear that big (nor it has any a priori
sign) the BD-field $\varphi$ can evolve in an appreciable way at the
cosmological level: it increases with the expansion if $\omega_{\rm BD}>0$,
and decreases with the expansion if $\omega_{\rm BD}<0$. In this context, its
initial value, $\varphi_{ini}$, becomes a relevant cosmological parameter,
which must be taken into account as a fitting parameter on equal footing with
$\omega_{\rm BD}$ and all of the conventional parameters entering the fit.
Using the large wealth of cosmological data, these parameters can be fixed at
the cosmological level without detriment of the observed physics at the local
domain, provided there is a screening mechanism insuring that $G_{\rm
eff}(\varphi)(0)$ stays sufficiently close to $G_{N}$ in the local
neighbourhood. The numerical results for this important scenario are presented
in Tables 3-6 and 9.

  * •

BD-Scenario III: Cassini-constrained Scenario. A more restricted version of
scenario II can appear if the ‘bare value’ of $\omega_{\rm BD}$ is as large as
in the Cassini bound. In such case $\omega_{\rm BD}$ is, obviously, perceived
large in both domains, local and cosmological. Even so, and despite of the
fact that $\varphi$ varies very slowly in this case, one can still exploit the
dependence of the fit on the initial value of the BD-field, $\varphi_{ini}$,
and use it as a relevant cosmological parameter. In practice, this situation
has only one additional degree of freedom as compared to Scenario I (and one
less than in Scenario II), but it is worth exploring – see our results in
Table 10. As these numerical results show, the Cassini bound still leaves
considerable freedom to the BD-$\Lambda$CDM model for improving the $H_{0}$
tension (without aggravating the $\sigma_{8}$ one) since the value of
$\varphi$ is still an active degree of freedom, despite its time evolution is
now more crippled.

  * •

BD-Scenario IV: Variable-$\omega_{\rm BD}(\varphi)$ Scenario. Here one admits
that the parameter $\omega_{\rm BD}$ is actually a function of the BD-field,
$\omega_{\rm BD}=\omega_{\rm BD}(\varphi)$, which can be modeled and adapted
to the constraints of the local and cosmological domains, or even combined
with the screening effects of the local Universe. We have said from the very
beginning that we would assume $\omega_{\rm BD}=$const. throughout our
analysis, and in fact we shall stick to that hypothesis; so here we mention
the variable $\omega_{\rm BD}$ scenario only for completeness. In any case, if
a function $\omega_{\rm BD}(\varphi)$ exists such that it takes very large
values in the local Universe while it takes much more moderate values in the
cosmological scales, that sort of scenario would be in the main tantamount to
Scenario II insofar as concerns its cosmological implications.

In our analysis we basically choose BD- Scenarios II and III (the latter being
a particular case of the former), which represent the most tolerant point of
view within the canonical $\omega_{\rm BD}=$const. option. Scenario II offers
the most powerful framework amenable to provide a cure for the tensions
afflicting the conventional $\Lambda$CDM model based on GR. Thus, we assume
that we can measure the cosmological value of the gravity strength in BD
theory – i.e. the value given in Eq. (7.80) – by using only cosmological data.
We combine the information from the LSS processes involving linear structure
formation with the background information obtained from low, intermediate, and
very high redshift probes, including BAO, CMB, and the distance ladder
measurement of $H_{0}$. The values of $\varphi_{ini}$ and $\omega_{\rm BD}$
are fitted to the data, and with them we obtain not only $\varphi(z=0)$ but we
determine the effective cosmological gravity strength at all epochs from
(7.80). The cosmological value of the gravity coupling can be considered as
the ‘naked’ or ‘bare’ value of the gravitational interaction, stripped from
screening effects of matter, in the same way as $\omega_{\rm BD}$ measured at
cosmological scales is the bare value freed of these effects. Even though
$G_{\rm eff}$ can be different from $G_{N},$ we do not object to that since it
can be ascribed to screening forces caused by the huge amounts of clustered
matter in the astrophysical environments. For this reason we do not adopt the
local constraints for our cosmological analysis presented in this paper, i.e.
we adhere to Scenario II as our main scenario. Remarkably enough, we shall see
that Scenario III still possesses a large fraction of the potentialities
inherent to Scenario II, notwithstanding the Cassini bound. In this sense
Scenarios II and III are both extremely interesting. A smoking gun of such
overarching possible picture is the possible detection of the dynamical dark
energy EoS encoded in the BD theory within the GR-picture (cf. Sec. 7.3),
which reveals itself in the form of effective quintessence, as well as through
the large smoothing achieved of the main tensions afflicting the conventional
$\Lambda$CDM. From here on, we present the bulk of our analysis and detailed
results after we have already discussed to a great extent their possible
implications.

###  Data and methodology

We fit the BD-$\Lambda$CDM together with the concordance GR-$\Lambda$CDM model
and the GR-XCDM (based on the XCDM parametrization of the DE [102]) to the
wealth of cosmological data compiled from distant type Ia supernovae (SNIa),
baryonic acoustic oscillations (BAO), a set of measurements of the Hubble
function at different redshifts, the large scale structure (LSS) formation
data encoded in $f(z_{i})\sigma_{8}(z_{i})$, and the CMB temperature and
low-$l$ polarization data from the Planck satellite. The joint combination of
all these individual datasets will constitute our Baseline Data configuration.
Moreover, we also study the repercussion of some alternative data, by adding
them to the aforementioned baseline setup. These additional datasets are: a
prior on the value of $H_{0}$ (or alternatively an effective calibration prior
on M) provided by the SH0ES collaboration; the CMB high-$l$ polarization and
lensing data from Planck; the Strong Lensing (SL) time delay angular diameter
distances from H0LICOW; and, finally, Weak Lensing (WL) data from KiDS. Below
we provide a brief description of the datasets employed in our analyses,
together with the corresponding references.

SNIa: We use the full Pantheon likelihood, which incorporates the information
from 1048 SNIa [268]. In addition, we also include the 207 SNIa from the DES
survey [337]. These two SNIa samples are uncorrelated, but the correlations
between the points within each sample are non-null and have been duly
incorporated in our analyses through the corresponding covariance matrices.

BAO: We use data on both, isotropic and anisotropic BAO analyses. We provide
the detailed list of data points and corresponding references in Table 1. A
few comments are in order about the use of some of the BAO data points
considered in this chapter. Regarding the Ly$\alpha$-forest data, we opt to
use the auto-correlation information from [338]. Excluding the Ly$\alpha$
cross-correlation data allows us to avoid double counting issues between the
latter and the eBOSS data from [273], due to the partial (although small)
overlap in the list of quasars employed in these two analyses. It is also
important to remark that we consider in our baseline dataset the BOSS data
reported in [147], which contains information from the spectrum (SP) and the
bispectrum (BP). The bispectrum information could capture some details
otherwise missed when only the spectrum is considered, so it is worth to use
it212121See also Ref.[339] for additional comments on the significance of the
bispectrum data as well as its potential implications on the possible
detection of dynamical dark energy.. Therefore, we study the possible
significance of the bispectrum component in the data by carrying out an
explicit comparison of the results obtained with the baseline configuration to
those obtained by substituting the data points from [147] with those from
[278], which only incorporate the SP information. The results are provided in
Tables 3 and 5, respectively. In Tables 4 and 6-8 we use the SP+BSP
combination [147]. In Table 10 we employ both SP and SP+BSP.

Survey | $z$ | Observable | Measurement | References
---|---|---|---|---
6dFGS+SDSS MGS | $0.122$ |  $D_{V}(r_{s}/r_{s,\rm fid})$[Mpc] |  $539\pm 17$[Mpc] | [272]
WiggleZ | $0.44$ |  $D_{V}(r_{s}/r_{s,\rm fid})$[Mpc] |  $1716.4\pm 83.1$[Mpc] | [155]
| $0.60$ |  $D_{V}(r_{s}/r_{s,\rm fid})$[Mpc] |  $2220.8\pm 100.6$[Mpc] |
| $0.73$ |  $D_{V}(r_{s}/r_{s,\rm fid})$[Mpc] |  $2516.1\pm 86.1$[Mpc] |
DR12 BOSS (BSP) | $0.32$ | $Hr_{s}/(10^{3}\rm km/s)$ | $11.549\pm 0.385$ | [147]
|  | $D_{A}/r_{s}$ | $6.5986\pm 0.1337$ |
| $0.57$ | $Hr_{s}/(10^{3}\rm km/s)$ | $14.021\pm 0.225$ |
|  | $D_{A}/r_{s}$ | $9.389\pm 0.1030$ |
DR12 BOSS (SP) | $0.38$ |  $D_{M}(r_{s}/r_{s,\rm fid})$[Mpc] | $1518\pm 22$ | [278]
|  |  $H(r_{s,\rm fid}/r_{s})$[km/s/Mpc] | $81.5\pm 1.9$ |
| $0.51$ |  $D_{M}(r_{s}/r_{s,\rm fid})$[Mpc] | $1977\pm 27$ |
|  |  $H(r_{s,\rm fid}/r_{s})$[km/s/Mpc] | $90.4\pm 1.9$ |
| $0.61$ |  $D_{M}(r_{s}/r_{s,\rm fid})$[Mpc] | $2283\pm 32$ |
|  |  $H(r_{s,\rm fid}/r_{s})$[km/s/Mpc] | $97.3\pm 2.1$ |
DES | $0.81$ | $D_{A}/r_{s}$ | $10.75\pm 0.43$ | [340]
eBOSS DR14 | $1.19$ | $Hr_{s}/(10^{3}\rm km/s)$ | $19.6782\pm 1.5867$ | [273]
|  | $D_{A}/r_{s}$ | $12.6621\pm 0.9876$ |
| $1.50$ | $Hr_{s}/(10^{3}\rm km/s)$ | $19.8637\pm 2.7187$ |
|  | $D_{A}/r_{s}$ | $12.4349\pm 1.0429$ |
| $1.83$ | $Hr_{s}/(10^{3}\rm km/s)$ | $26.7928\pm 3.5632$ |
|  | $D_{A}/r_{s}$ | $13.1305\pm 1.0465$ |
Ly$\alpha$-F eBOSS DR14 | $2.34$ | $D_{H}/r_{s}$ | $8.86\pm 0.29$ | [338]
|  | $D_{M}/r_{s}$ | $37.41\pm 1.86$ |

Table 1: Published values of BAO data, see the quoted references and text in
Sec. 7.6. Although we include in this table the values of
$D_{H}/r_{s}=c/(r_{s}H)$ and $D_{M}/r_{s}$ for the Ly$\alpha$-forest auto-
correlation data from [338], we have performed the fitting analysis with the
full likelihood. The fiducial values of the comoving sound horizon appearing
in the table are $r_{s,{\rm fid}}=147.5$ Mpc for [272], $r_{s,{\rm
fid}}=148.6$ Mpc for [155], and $r_{s,{\rm fid}}=147.78$ Mpc for [278].

Cosmic Chronometers: We use the 31 data points on $H(z_{i})$, at different
redshifts, from [137, 138, 141, 139, 136, 142, 140, 270]. All of them have
been obtained making use of the differential age technique applied to
passively evolving galaxies [341], which provides cosmology-independent
constraints on the Hubble function, but are still subject to systematics
coming from the choice of the stellar population synthesis technique, and also
the potential contamination of young stellar components in the quiescent
galaxies [342, 343, 344]. For this reason we consider a more conservative
dataset that takes into account these additional uncertainties. To be
concrete, we use the processed sample presented in Table 2 of [345]. See
therein for further details.

CMB: In our baseline dataset we consider the full Planck 2018 TT+lowE
likelihood [31]. In order to study the influence of the CMB high-$l$
polarizations and lensing we consider two alternative (non-baseline) datasets,
in which we substitute the TT+lowE likelihood by: (i) the TTTEEE+lowE
likelihood, which incorporates the information of high multipole
polarizations; (ii) the full TTTEEE+lowE+lensing likelihood, in which we also
incorporate the Planck 2018 lensing data. In Tables 6 and 7 these scenarios
are denoted as B+$H_{0}$+pol and B+$H_{0}$+pol+lens, respectively.

Survey | $z$ | $f(z)\sigma_{8}(z)$ | References
---|---|---|---
6dFGS+2MTF | $0.03$ | $0.404^{+0.082}_{-0.081}$ | [346]
SDSS-DR7 | $0.10$ | $0.376\pm 0.038$ | [276]
GAMA | $0.18$ | $0.29\pm 0.10$ | [145]
| $0.38$ | $0.44\pm 0.06$ | [146]
DR12 BOSS (BSP) | $0.32$ | $0.427\pm 0.056$ | [147]
| $0.57$ | $0.426\pm 0.029$ |
WiggleZ | $0.22$ | $0.42\pm 0.07$ | [148]
| $0.41$ | $0.45\pm 0.04$ |
| $0.60$ | $0.43\pm 0.04$ |
| $0.78$ | $0.38\pm 0.04$ |
DR12 BOSS (SP) | $0.38$ | $0.497\pm 0.045$ | [278]
| $0.51$ | $0.458\pm 0.038$ |
| $0.61$ | $0.436\pm 0.034$ |
VIPERS | $0.60$ | $0.49\pm 0.12$ | [277]
| $0.86$ | $0.46\pm 0.09$ |
VVDS | $0.77$ | $0.49\pm 0.18$ |  [151],[152]
FastSound | $1.36$ | $0.482\pm 0.116$ | [274]
eBOSS DR14 | $1.19$ | $0.4736\pm 0.0992$ | [273]
| $1.50$ | $0.3436\pm 0.1104$ |
| $1.83$ | $0.4998\pm 0.1111$ |

Table 2: Published values of $f(z)\sigma_{8}(z)$, see the quoted references
and text in Sec. 7.6.

LSS: In this chapter the LSS dataset contains the data points on the product
of the ordinary growth rate $f(z_{i})$ with $\sigma_{8}(z_{i})$ at different
effective redshifts. They are all listed in Table 2, together with the
references of interest. In order to correct the potential bias introduced by
the particular choice of a fiducial model in the original observational
analyses we apply the rescaling correction explained in [238]. See also Sec.
II.2 of [347]. The internal correlations between the BAO and RSD data from
[147], [278] and [273] have been duly taken into account through the
corresponding covariance matrices provided in these three references.

Prior on $H_{0}$: We include as a prior in almost all the non-baseline
datasets the value of the Hubble parameter measured by the SH0ES
collaboration, $H_{0}=(73.5\pm 1.4)$ km/s/Mpc [84]. It is obtained with the
cosmic distance ladder method using an improved calibration of the Cepheid
period-luminosity relation. As a possible alternative, we also consider the
case in which we use instead the SH0ES effective calibration prior on the
absolute magnitude $M$ of the SNIa, as provided in [348]: $M=-19.2191\pm
0.0405$. It is obtained from the calibration of ‘nearby’ SNIa (at $z\lesssim
0.01$) with Cepheids [321]. It has been recently argued in [349, 348] (and
later on also in [350, 351]) that in cosmological studies it is better to use
this SH0ES constraint rather than the direct prior on $H_{0}$ when combined
with low-redshift SNIa data to avoid double counting issues. We show that the
results obtained with these two methods are compatible and lead to completely
consistent results (see the discussion in Sec. 7.7, and also Tables 3 and 6).

SL: In one of the non-baseline datasets we use in combination with the SH0ES
prior on $H_{0}$ the data extracted from the six gravitational lensed quasars
of variable luminosity reported by the H0LICOW team [97]. In order to know all
the details of the methodology to include these data in the analysis see
Appendix E.

WL: We include the information provided by the survey The Kilo-Degree Survey
(KiDS) [75, 245, 352, 353] in one of our non-baseline scenarios. The reason
why we prefer to keep away this data from the Baseline combination is because
in order to model the non-linear effects it is necessary to choose a
particular cosmological model. As a consequence this data is model-dependent.
See Appendix E for more details.

We study the performance of the BD-$\Lambda$CDM, GR-$\Lambda$CDM and GR-XCDM
models under different datasets. In the following we briefly summarize the
composition of each of them:

  * •

Baseline (B): Here we include the Planck 2018 TT+lowE CMB data, together with
SNIa +BAO+$H(z_{i})$+LSS (see Table 3 and 8). It is important to remark that
for the BOSS BAO+LSS data we consider [147], which includes the information
from the spectrum (SP) as well as from the bispectrum (BSP). We construct some
other datasets using this baseline configuration as the main building block.
See the other items, below.

  * •

Baseline+$H_{0}$ (B+$H_{0}$): Here we add the SH0ES prior on the $H_{0}$
parameter from [84] to the baseline dataset (see again Tables 3 and 8).

  * •

Baseline+$H_{0}$+SL: The inclusion of the Strong Lensing (SL) data from [97]
exacerbates more the $H_{0}$-tension in the context of the GR-$\Lambda$CDM
model (see e.g. [354] and Sec. 7.8), so it is interesting to also study the
ability of the BD-$\Lambda$CDM to fit the SL data when they are combined with
the previous B+$H_{0}$ dataset, and compare the results with those obtained
with the GR-$\Lambda$CDM. The corresponding fitting results are displayed in
Table 4.

  * •

Spectrum: In this dataset we replace the SP+BSP data from [147] used in the
Baseline dataset (see the first item, above) by the data from [278], which
only contains the spectrum (SP) information (i.e. the usual matter power
spectrum).

  * •

Spectrum+$H_{0}$: As in the preceding item, but including the $H_{0}$ prior
from SH0ES [84].

The aforementioned datasets are all based on the BD-Scenario II (cf. Sec. 7.5)
and can be considered as the main ones (cf. Tables 3-5, and 8), nevertheless
we also consider a bunch of alternative datasets (also based on the BD-
Scenario II). We present the corresponding results in Table 6 and the first
five columns of Table 7 for the BD-$\Lambda$-CDM and GR-$\Lambda$CDM models,
respectively.

  * •

B+$M$: In this scenario we replace the prior on $H_{0}$ [84] employed in the
B+$H_{0}$ dataset with the effective SH0ES calibration prior on the absolute
magnitude of SNIa $M$ provided in [348].

  * •

B+$H_{0}$+pol: Here we add the CMB high-$l$ polarization data from Planck 2018
[31] to the B+$H_{0}$ dataset described before, i.e. we consider the Planck
2018 TTTEEE+lowE likelihood for the CMB.

  * •

B+$H_{0}$+pol+lens: In addition to the datasets considered in the above case
we also include the CMB lensing data from Planck 2018 [31], i.e. we use the
Planck 2018 TTTEEE+lowE+ lensing likelihood.

  * •

B+$H_{0}$+WL: In this alternative case we consider the Weak Lensing (WL) data
from KiDS [75, 352], together with the B+$H_{0}$ dataset.

  * •

CMB+BAO+SNIa: By considering only this data combination, we study the
performance of the the BD-$\Lambda$CDM and the GR-$\Lambda$CDM models under a
more limited dataset, obtained upon the removal of the data that trigger the
$H_{0}$ and $\sigma_{8}$ tensions, i.e. the prior on $H_{0}$ from SH0ES and
the LSS data. The use of the BAO+SNIa data helps to break the strong
degeneracies found in parameter space when only the CMB is considered. Here we
use the TT+lowE Planck 2018 likelihood [31].

Finally, in Table 10 we present the results obtained for the BD-$\Lambda$CDM
in the context of the Cassini-constrained scenario, or Scenario III (see Sec.
7.5 for the details). The corresponding results for the GR-$\Lambda$CDM are
shown in the third column of Table 3, and the last three rows of Table 7. In
all these datasets we include the Cassini bound [317] (see Sec. 7.1 for
details). The main purpose of this scenario is to test the ability of the
BD-$\Lambda$CDM to fit the observational data with $\epsilon_{\rm BD}\simeq 0$
and $\varphi\neq 1$.

  * •

B+$H_{0}$+Cassini: It contains the very same datasets as in the
Baseline+$H_{0}$ case, but here we also include the Cassini constraint.

  * •

B+$H_{0}$+Cassini (No LSS): Here we study the impact of the LSS data in the
context of Scenario III, by removing them from the previous B+$H_{0}$+Cassini
dataset.

  * •

Dataset [335]: To ease the comparison with the results obtained in [335], here
we use exactly the same dataset as in that reference, namely: the Planck 2018
TTTEEE+lowE+lensing likelihood [31], the BAO data from [153, 154, 278], and
the SH0ES prior from [321], $H_{0}=74.03\pm 1.42$ km/s/Mpc.

  * •

Dataset[335]+LSS: Here we consider an extension of the previous scenario by
adding the LSS data on top of the data from [335].

We believe that all the datasets and scenarios studied in this work cover a
wide range of possibilities and show in great detail which is the
phenomenological performance of the BD-$\Lambda$CDM, GR-$\Lambda$CDM and GR-
XCDM models.
The speed of gravitational waves at $z\approx 0$, $c_{gw}$, has been recently
constrained to be extremely close to the speed of light, $|c_{gw}/c-1|\lesssim
5\cdot 10^{-16}$ [355]. In the BD-$\Lambda$CDM model $c_{gw}=c\,$
$\forall{z}$, so this constraint is automatically fulfilled. We have also
checked that the BD-$\Lambda$CDM respects the bound on $G(\varphi)$ at the Big
Bang Nucleosynthesis (BBN) epoch, $|G(\varphi_{\rm BBN})/G_{N}-1|\lesssim 0.1$
[168], since $G(\varphi_{\rm BBN})\simeq G(\varphi_{ini})$ and our best-fit
values satisfy $G(\varphi_{ini})>0.9G_{N}$ regardless of the dataset under
consideration, see the fitting results in Tables 3-6 and 10.

To obtain the posterior distributions and corresponding constraints for the
various dataset combinations described above we have run the Monte Carlo
sampler MontePython222222http://baudren.github.io/montepython.html [266]
together with the Einstein-Boltzmann system solver
CLASS232323http://lesgourg.github.io/class_public/class.html [265]. We have
duly modified the latter to implement the background and linear perturbations
equations of the BD-$\Lambda$CDM model. We use adiabatic initial conditions
for all matter species. Let us note that the initial perturbation of the BD-
field and its time derivative can be set to zero, as the DM velocity
divergence when the synchronous gauge is employed, see Appendix C.5 of [95]
for a brief discussion. To get the contour plots and one-dimensional posterior
distributions of the parameters entering the models we have used the Python
package GetDist242424https://getdist.readthedocs.io/en/latest/ [356], and to
compute the full Bayesian evidences for the different models and dataset
combinations, we have employed the code
MCEvidence252525https://github.com/yabebalFantaye/MCEvidence [286]. The
Deviance Information Criterion (DIC) [357] has been computed with our own
numerical code. The results are displayed in Tables 3-10, and also in Figs.
9-11. They are discussed in the next section.

###  Numerical analysis. Results

In the following we put the models under consideration to the test, using the
various datasets described in Sec. 7.6. We perform the statistical analysis of
the models in terms of a joint likelihood function, which is the product of
the individual likelihoods for each data source and includes the corresponding
covariance matrices. For a fairer comparison with the GR-$\Lambda$CDM we use
standard information criteria in which the presence of extra parameters in a
given model is conveniently penalized so as to achieve a balanced comparison
with the model having less parameters. More concretely, we employ the full
Bayesian evidence to duly quantify the fitting ability of the BD-$\Lambda$CDM
model as compared to its GR analogue. Given a dataset $\mathcal{D}$, the
probability of a certain model $M_{i}$ to be the best one among a given set of
models $\\{M\\}$ reads,

Figure 9: The full Bayesian evidence curves for the BD-$\Lambda$CDM as
compared to the GR-$\Lambda$CDM, using different datasets and as a function of
$\Delta_{X}=\Delta\epsilon_{\rm BD}/10^{-2}=\Delta\varphi_{ini}/0.2$, with
$\Delta\epsilon_{\rm BD}$ and $\Delta\varphi_{ini}$ being the (flat) prior
ranges for $\epsilon_{\rm BD}$ and $\varphi_{ini}$, respectively. The curves
are computed using the exact evidence formula, Eq. (7.83). The marked evidence
ranges conform with the conventional Jeffreys’ scale, see the main text in
Sec. 7.7.

$P(M_{i}|\mathcal{D})=\frac{P(M_{i})\mathcal{E}(\mathcal{D}|M_{i})}{P(\mathcal{D})}\,,$
(7.82)

where $P(M_{i})$ is the prior probability of the model $M_{i}$,
$P(\mathcal{D})$ the probability of having the dataset $\mathcal{D}$, and the
normalization condition $\sum_{j\in\\{M\\}}P(M_{j})=1$ is assumed. The
quantity $\mathcal{E}(\mathcal{D}|M_{i})$ is the so-called marginal likelihood
or evidence [55]. If the model $M_{i}$ has $n$ parameters contained in the
vector $\vec{p}^{M_{i}}=(p^{M_{i}}_{1},p^{M_{i}}_{2},...,p^{M_{i}}_{n})$, the
evidence takes the following form:

| Baseline | Baseline+$H_{0}$
---|---|---
Parameter | GR-$\Lambda$CDM | BD-$\Lambda$CDM | GR-$\Lambda$CDM | BD-$\Lambda$CDM
$H_{0}$ (km/s/Mpc) | $68.20^{+0.41}_{-0.40}$ | $68.86^{+1.15}_{-1.24}$ | $68.57^{+0.36}_{-0.42}$ | $70.83^{+0.92}_{-0.95}$
$\omega_{b}$ | $0.02227^{+0.00019}_{-0.00018}$ | $0.02251^{+0.00026}_{-0.00027}$ | $0.02238\pm 0.00019$ | $0.02275^{+0.00024}_{-0.00026}$
$\omega_{cdm}$ | $0.11763^{+0.00090}_{-0.00092}$ | $0.11598^{+0.00159}_{-0.00152}$ | $0.11699^{+0.00092}_{-0.00083}$ | $0.11574^{+0.00164}_{-0.00158}$
$\tau$ | $0.050^{+0.004}_{-0.008}$ | $0.052^{+0.006}_{-0.008}$ | $0.051^{+0.005}_{-0.008}$ | $0.053^{+0.006}_{-0.008}$
$n_{s}$ | $0.9683^{+0.0039}_{-0.0038}$ | $0.9775^{+0.0084}_{-0.0086}$ | $0.9703^{+0.0038}_{-0.0036}$ | $0.9873^{+0.0076}_{-0.0075}$
$\sigma_{8}$ | $0.797^{+0.005}_{-0.006}$ | $0.785\pm 0.013$ | $0.796^{+0.006}_{-0.007}$ | $0.789\pm 0.013$
$r_{s}$ (Mpc) | $147.83^{+0.29}_{-0.30}$ | $145.89^{+2.26}_{-2.49}$ | $147.88\pm 0.31$ | $142.46^{+1.84}_{-1.86}$
$\epsilon_{\rm BD}$ | - | $-0.00184^{+0.00140}_{-0.00142}$ | - | $-0.00199^{+0.00142}_{-0.00147}$
$\varphi_{ini}$ | - | $0.974^{+0.027}_{-0.031}$ | - | $0.932^{+0.022}_{-0.023}$
$\varphi(0)$ | - | $0.960^{+0.032}_{-0.037}$ | - | $0.918^{+0.027}_{-0.029}$
$w_{\rm eff}(0)$ | - | $-0.983^{+0.015}_{-0.014}$ | - | $-0.966^{+0.012}_{-0.011}$
$\dot{G}(0)/G(0)(10^{-13}yr^{-1})$ | - | $2.022^{+1.585}_{-1.518}$ | - | $2.256^{+1.658}_{-1.621}$
$\chi^{2}_{min}$ | 2271.98 | 2271.82 | 2285.50 | 2276.04
$2\ln B$ | - | -2.26 | - | +4.92
$\Delta{\rm DIC}$ | - | -0.54 | - | +4.90

Table 3: The mean fit values and $68.3\%$ confidence limits for the considered
models using our baseline dataset in the first block, i.e.
SNIa+$H(z)$+BAO+LSS+CMB TT data, and baseline+$H_{0}$ in the second one (cf.
Sec. 7.6 for details). These results have been obtained within our main BD
scenario (Scenario II of Sec. 7.5). In all cases a massive neutrino of $0.06$
eV has been included. First we display the fitting results for the six
conventional parameters, namely: $H_{0}$, the reduced density parameters for
baryons ($\omega_{b}=\Omega_{b}h^{2}$) and CDM
($\omega_{cdm}=\Omega_{cdm}h^{2}$), the reionization optical depth $\tau$, the
spectral index $n_{s}$ of the primordial power-law power spectrum, and, for
convenience, instead of the amplitude $A_{s}$ of such spectrum we list the
values of $\sigma_{8}$. We also include the sound horizon at the baryon drag
epoch, $r_{s}$, obtained as a derived parameter. Right after we list the
values of the free parameters that characterize the BD model: $\epsilon_{\rm
BD}$ (7.1) and the initial condition for the BD-field, $\varphi_{ini}$. We
also include the values of the BD-field, the (exact) effective EoS parameter
(7.58), and the ratio between the derivative and the value of Newton’s
coupling, all computed at $z=0$. Finally, we report the values of the minimum
of the $\chi^{2}$-function, $\chi^{2}_{min}$, the exact Bayes ratios (computed
under the conditions explained in the main text of Sec. 7.7), and the DIC. It
is also worth to remark that the baseline dataset employed here includes the
contribution not only of the spectrum, but also the bispectrum information
from BOSS [147], see Sec. 7.6 for details.

$\mathcal{E}(\mathcal{D}|M_{i})=\int\mathcal{L}(\mathcal{D}|\vec{p}^{M_{i}},M_{i})\pi(\vec{p}^{M_{i}})d^{n}p^{M_{i}}\,,$
(7.83)

with $\mathcal{L}(\mathcal{D}|\vec{p}^{M_{i}},M_{i})$ the likelihood and
$\pi(\vec{p}^{M_{i}})$ the prior of the parameters entering the model $M_{i}$.
The evidence is larger for those models that have more overlapping volume
between the likelihood and the prior distributions, but penalizes the use of
additional parameters having a non-null impact on the likelihood. Hence, the
evidence constitutes a good way of quantifying the performance of the model by
implementing in practice the Occam razor principle. We can compare the fitting
performance of BD-$\Lambda$CDM and GR-$\Lambda$CDM models by assuming equal
prior probability for both of them, i.e. $P({\rm BD-}\Lambda{\rm CDM})=P({\rm
GR-}\Lambda{\rm CDM})$ (“Principle of Insufficient Reason”). The ratio of
their associated probabilities can then be directly written as the ratio of
their corresponding evidences, i.e.

| Baseline+$H_{0}$+SL
---|---
Parameter | GR-$\Lambda$CDM | BD-$\Lambda$CDM
$H_{0}$ (km/s/Mpc) | $68.74^{+0.37}_{-0.40}$ | $71.30^{+0.80}_{-0.84}$
$\omega_{b}$ | $0.02242^{+0.00018}_{-0.00019}$ | $0.02281\pm 0.00025$
$\omega_{cdm}$ | $0.11666^{+0.00087}_{-0.00086}$ | $0.11560^{+0.00158}_{-0.00169}$
$\tau$ | $0.051^{+0.005}_{-0.008}$ | $0.053^{+0.006}_{-0.008}$
$n_{s}$ | $0.9708^{+0.0036}_{-0.0038}$ | $0.9901^{+0.0075}_{-0.0070}$
$\sigma_{8}$ | $0.795^{+0.006}_{-0.007}$ | $0.789\pm 0.013$
$r_{s}$ (Mpc) | $147.93^{+0.30}_{-0.31}$ | $141.68^{+1.69}_{-1.73}$
$\epsilon_{\rm BD}$ | - | $-0.00208^{+0.00151}_{-0.00140}$
$\varphi_{ini}$ | - | $0.923^{+0.019}_{-0.021}$
$\varphi(0)$ | - | $0.908^{+0.026}_{-0.028}$
$w_{\rm eff}(0)$ | - | $-0.962\pm 0.011$
$\dot{G}(0)/G(0)(10^{-13}yr^{-1})$ | - | $2.375^{+1.612}_{-1.721}$
$\chi^{2}_{min}$ | 2320.40 | 2305.80
$2\ln B$ | - | $+9.22$
$\Delta{\rm DIC}$ | - | $+9.15$

Table 4: Fitting results as in Table 3, but adding the Strong Lensing data,
i.e. we use here the dataset Baseline+$H_{0}$+SL, for both the GR-$\Lambda$CDM
and the BD-$\Lambda$CDM.

| Spectrum | Spectrum+$H_{0}$
---|---|---
Parameter | GR-$\Lambda$CDM | BD-$\Lambda$CDM | GR-$\Lambda$CDM | BD-$\Lambda$CDM
$H_{0}$ (km/s/Mpc) | $68.00^{+0.47}_{-0.48}$ | $68.86^{+1.27}_{-1.34}$ | $68.61^{+0.46}_{-0.49}$ | $70.94^{+1.00}_{-0.98}$
$\omega_{b}$ | $0.02223^{+0.00020}_{-0.00021}$ | $0.02241\pm 0.00027$ | $0.02239^{+0.00020}_{-0.00019}$ | $0.02264^{+0.00026}_{-0.00025}$
$\omega_{cdm}$ | $0.11809^{+0.00112}_{-0.00095}$ | $0.11743^{+0.00168}_{-0.00170}$ | $0.11695\pm 0.00104$ | $0.11702^{+0.00169}_{-0.00167}$
$\tau$ | $0.051^{+0.005}_{-0.008}$ | $0.053^{+0.006}_{-0.008}$ | $0.053^{+0.006}_{-0.008}$ | $0.054^{+0.007}_{-0.008}$
$n_{s}$ | $0.9673^{+0.0039}_{-0.0044}$ | $0.9742^{+0.0086}_{-0.0090}$ | $0.9705^{+0.0040}_{-0.0041}$ | $0.9845^{+0.0076}_{-0.0077}$
$\sigma_{8}$ | $0.800^{+0.006}_{-0.007}$ | $0.798\pm 0.014$ | $0.798\pm 0.007$ | $0.801^{+0.015}_{-0.014}$
$r_{s}$ (Mpc) | $147.75^{+0.31}_{-0.35}$ | $145.82^{+2.33}_{-2.52}$ | $147.88^{+0.33}_{-0.32}$ | $142.55^{+1.71}_{-1.96}$
$\epsilon_{\rm BD}$ | - | $-0.00079^{+0.00158}_{-0.00157}$ | - | $-0.00081^{+0.00162}_{-0.00165}$
$\varphi_{ini}$ | - | $0.976^{+0.028}_{-0.032}$ | - | $0.935^{+0.020}_{-0.024}$
$\varphi(0)$ | - | $0.970^{+0.034}_{-0.038}$ | - | $0.929^{+0.028}_{-0.030}$
$w_{\rm eff}(0)$ | - | $-0.987^{+0.016}_{-0.014}$ | - | $-0.971^{+0.013}_{-0.011}$
$\dot{G}(0)/G(0)(10^{-13}yr^{-1})$ | - | $0.864^{+1.711}_{-1.734}$ | - | $0.913^{+1.895}_{-1.791}$
$\chi^{2}_{min}$ | 2269.04 | 2268.28 | 2283.66 | 2274.64
$2\ln B$ | - | $-2.94$ | - | $+3.98$
$\Delta{\rm DIC}$ | - | $-3.36$ | - | $+4.76$

Table 5: As in Table 3, but replacing the BOSS BAO+LSS data from [147] with
those from [278], which only includes the spectrum information. See the
discussion of the results in Sec. 7.7.

$\frac{P({\rm BD-}\Lambda{\rm CDM}|\mathcal{D})}{P({\rm GR-}\Lambda{\rm
CDM}|\mathcal{D})}=\frac{\mathcal{E}(\mathcal{D}|{\rm BD-}\Lambda{\rm
CDM})}{\mathcal{E}(\mathcal{D}|{\rm GR-}\Lambda{\rm CDM})}\equiv B\,,$ (7.84)

where $B$ is the so-called Bayes ratio (or Bayes factor) and is the quantity
we are interested in. Notice that when $B>1$ this means that data prefer the
BD-$\Lambda$CDM model over the GR version, but of course depending on how much
larger than 1 it is we will have different levels of statistical significance
for such preference. It is common to adopt in the literature the so-called
Jeffreys’ scale to categorize the level of evidence that one can infer from
the computed value of the Bayes ratio. Jeffrey’s scale actually is usually
written not directly in terms of $B$, but in terms of $2\ln B$. The latter is
sometimes estimated with a simple Schwarz (or Bayesian) information criterion
$\Delta$BIC [186, 242], although $2\ln B$ is a much more rigorous,
sophisticated (and difficult to compute) statistics than just the usual
$\Delta$BIC estimates based on using the minimum value of $\chi^{2}$, the
number of points and the number of independent fitting parameters. If $2\ln B$
lies below $2$ in absolute value, then we can conclude that the evidence in
favor of BD-$\Lambda$CDM (against GR-$\Lambda$CDM) is at most only weak, and
in all cases not conclusive; if $2<2\ln B<6$ the evidence is said to be
positive; if, instead, $6<2\ln B<10$, then it is considered to be strong,
whereas if $2\ln B>10$ one is entitled to speak of very strong evidence in
favor of the BD-$\Lambda$CDM over the GR-$\Lambda$CDM model. For more
technical details related with the evidence and the Bayes ratio we refer the
reader to [242, 187, 55]. Notice that the computation of (7.84) is not easy in
general; in fact, it can be rather cumbersome since we usually work with
models with a high number of parameters, so the multiple integrals that we
need to compute become quite demanding from the computational point of view.
We have calculated the evidences numerically, of course, processing the Markov
chains obtained from the Monte Carlo analyses carried out with
CLASS+MontePython [265, 266] with the numerical code MCEvidence [286], which
is publicly available (cf. Sec. 7.6). We report the values obtained for $2\ln
B$ (7.84) in Tables 3-6, 8, and 10. Table 3 contains the fitting results for
the BD- and GR$-\Lambda$CDM models obtained with the Baseline and
Baseline+$H_{0}$ datasets. In Table 4 we present the results for the
Baseline+$H_{0}$+SL dataset. In Table 5 we show the output of the fitting
analysis for the same models and using the same data as in Table 3, but
changing the BOSS data from [147], which contain both the mater spectrum and
bispectrum information, by the BOSS data from [278], which only incorporate
the spectrum part (i.e. the usual matter power spectrum). In Table 6 we plug
the results obtained for the BD-$\Lambda$CDM with the alternative datasets
described in Sec. 7.6, and in Table 7 we show the corresponding results for
the GR-$\Lambda$CDM. Next, in Table 8 we present the results with the Baseline
and Baseline+$H_{0}$ data configurations obtained using the GR-XCDM
parametrization. In Table 9 we display the values of the parameters
$\sigma_{8}$ and $S_{8}$ for the GR- and BD-$\Lambda$CDM models, as well as
the parameter $\tilde{S}_{8}=S_{8}/\sqrt{\varphi(0)}$ for the BD. Finally, in
Table 10 we present the fitting results for the BD-$\Lambda$CDM model,
considering the (Cassini-constrained) Scenario III described in Sec. 7.5.

Datasets | $H_{0}$ | $\omega_{m}$ | $\sigma_{8}$ |  $r_{s}$ (Mpc) | $\epsilon_{\rm BD}\cdot 10^{3}$ | $w_{\rm eff}(0)$ | $2\ln B$
---|---|---|---|---|---|---|---
B+$M$ | $71.19^{+0.92}_{-1.02}$ | $0.1390^{+0.0014}_{-0.0015}$ | $0.788^{+0.012}_{-0.013}$ | $141.87^{+2.06}_{-1.81}$ | $-2.16^{+1.42}_{-1.36}$ | $-0.963^{+0.012}_{-0.011}$ | $+10.38$
B+$H_{0}$+pol | $69.85^{+0.81}_{-0.85}$ | $0.1409^{+0.0012}_{-0.0011}$ | $0.801\pm 0.011$ | $144.72^{+1.51}_{-1.83}$ | $-0.30^{+1.20}_{-1.23}$ | $-0.985^{+0.012}_{-0.009}$ | $-1.44$
B+$H_{0}$+pol+lens | $69.74^{+0.82}_{-0.77}$ | $0.1416^{+0.0011}_{-0.0010}$ | $0.808\pm 0.009$ | $144.66^{+1.56}_{-1.61}$ | $0.00^{+1.11}_{-1.07}$ | $-0.986\pm 0.010$ | $-1.98$
B+$H_{0}$+WL | $70.69^{+0.91}_{-0.90}$ | $0.1398^{+0.0015}_{-0.0013}$ | $0.794^{+0.011}_{-0.012}$ | $142.76^{+1.79}_{-1.86}$ | $-1.42^{+1.29}_{-1.37}$ | $-0.970^{+0.011}_{-0.010}$ | $+4.34$
CMB+BAO+SNIa | $68.63^{+1.44}_{-1.50}$ | $0.1425\pm 0.0019$ | $0.818^{+0.017}_{-0.018}$ | $146.62^{+2.92}_{-2.93}$ | $1.14^{+1.84}_{-1.68}$ | $-0.999^{+0.020}_{-0.017}$ | $-3.00$

Table 6: Fitting results for the BD-$\Lambda$CDM model obtained with some
alternative datasets and in all cases within the main BD Scenario II. Due to
the lack of space, we employ some abbreviations, namely: B for the Baseline
dataset described in Sec. 7.6; pol for the Planck 2018 (TE+EE) high-$l$
polarization data; and lens for the CMB lensing. The $\omega_{m}$ parameter
contains the contribution of baryons and dark matter. In the last row, CMB
refers to the TT+lowE Planck 2018 likelihood (cf. Sec. 7.6 for details on the
data). $H_{0}$ is given in km/s/Mpc. For a discussion of the results, see Sec.
7.7.

The evidence (7.83) clearly depends on the priors for the parameters entering
the model. For the 6 parameters in common in the BD- and GR-$\Lambda$CDM
models, namely,
$(\omega_{b}=\Omega_{b}h^{2},\omega_{cdm}=\Omega_{cdm}h^{2},H_{0},\tau,A_{s},n_{s})$,
if we use the same flat priors in both models they cancel exactly in the
computation of the Bayes ratio (7.84). Thus, the latter does not depend on the
priors for these parameters if their ranges are big enough so as to not alter
the shape of the likelihood severely. The Bayes ratio is, though, sensitive to
the priors for the two additional parameters introduced in the BD-$\Lambda$CDM
model in our Scenario II (cf. Sec. 7.5), i.e. $\varphi_{ini}$ and
$\epsilon_{\rm BD}$, since they are not canceled in (7.84). We study the
dependence of the evidence on these priors in Fig. 9, where we plot $2\ln B$
obtained for the BD-$\Lambda$CDM model from different datasets, and as a
function of a quantity that we call $\Delta_{X}$, defined as

$\Delta_{X}\equiv\frac{\Delta\epsilon_{\rm
BD}}{10^{-2}}=\frac{\Delta\varphi_{ini}}{0.2}\,,$ (7.85)

with $\Delta\epsilon_{\rm BD}$ and $\Delta\varphi_{ini}$ being the flat prior
ranges of $\epsilon_{\rm BD}$ and $\varphi_{ini}$, centered at $\epsilon_{\rm
BD}=0$ and $\varphi_{ini}=1$, respectively. $\Delta_{X}$ will be equal to one
when $\Delta\epsilon_{\rm BD}=10^{-2}$ and $\Delta\varphi_{ini}=0.2$, which
are natural values for these prior ranges. The former implies $\omega_{\rm
BD}>100$ and the latter could be associated to the range
$0.9<\varphi_{ini}<1.1$. We do not expect $\omega_{\rm BD}\lesssim 100$ since
this would imply an exceedingly large departure from GR, even at cosmological
scales, where this lower bound was already set using the first releases of
WMAP CMB data almost twenty years ago, see e.g. [358, 359]. Regarding the
prior range $0.9<\varphi_{ini}<1.1$, it is also quite natural, since this is
necessary to satisfy the BBN bounds [168]. In all tables we report the values
of $2\ln B$ obtained by setting the natural value $\Delta_{X}=1$, and in Fig.
9, as mentioned before, we also show how this quantity changes with the prior
width, in terms of the variable $\Delta_{X}$ (7.85).
In the Cassini-constrained Scenario III (cf. again Sec. 7.5), we also allow
variations of $\varphi_{ini}$ and $\epsilon_{\rm BD}$ in our Monte Carlo runs,
of course, but the natural prior range for $\epsilon_{\rm BD}$ is now much
smaller than in Scenario II, since now we expect it to be more constrained by
the local observations. It is more natural to take in this case a prior range
$\Delta\epsilon_{\rm BD}=5\cdot 10^{-5}$ (still larger than Cassini’s bound),
and this is what we do in all the analyses of this scenario. See the comments
in Sec. 7.8, and Table 10.
In Tables 3-5 apart from the Bayes ratio, we also include the Deviance
Information Criterion [357], which is strictly speaking an approximation of
the exact Bayesian approach that works fine when the posterior distributions
are sufficiently Gaussian. The DIC is defined as

${\rm DIC}=\chi^{2}(\hat{\theta})+2p_{D}\,.$ (7.86)

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
generated with MontePython. To compare the ability of the BD- and
GR-$\Lambda$CDM models to fit the data, one has to compute the respective
differences of DIC values between the first and second models. They are
denoted $\Delta$DIC in our tables, and this quantity is the analogous to $2\ln
B$.

####  Comparing with the XCDM parametrization

Datasets | $H_{0}$ | $\omega_{m}$ | $\sigma_{8}$ |  $r_{s}$ (Mpc)
---|---|---|---|---
B+$M$ | $68.64^{+0.39}_{-0.38}$ | $0.1398\pm 0.0009$ | $0.796^{+0.005}_{-0.006}$ | $147.87^{+0.29}_{-0.30}$
B+$H_{0}$+pol | $68.50^{+0.33}_{-0.36}$ | $0.1408^{+0.0007}_{-0.0008}$ | $0.799\pm 0.006$ | $147.53\pm 0.022$
B+$H_{0}$+pol+lens | $68.38^{+0.35}_{-0.33}$ | $0.1411\pm 0.0007$ | $0.803\pm 0.006$ | $147.44^{+0.22}_{-0.21}$
B+$H_{0}$+WL | $68.64\pm 0.37$ | $0.1399^{+0.0008}_{-0.0009}$ | $0.795^{+0.005}_{-0.007}$ | $147.89^{+0.31}_{-0.29}$
CMB+BAO+SNIa | $67.91^{+0.39}_{-0.41}$ | $0.1413\pm 0.0009$ | $0.805^{+0.006}_{-0.007}$ | $147.66^{+0.31}_{-0.29}$
B+$H_{0}$ (No LSS) | $68.38^{+0.42}_{-0.37}$ | $0.1405\pm 0.0009$ | $0.802^{+0.007}_{-0.008}$ | $147.76\pm 0.30$
Dataset [335] | $68.17^{+0.43}_{-0.44}$ | $0.1416\pm 0.0009$ | $0.810^{+0.006}_{-0.007}$ | $147.35^{+0.22}_{-0.24}$
Dataset [335] \+ LSS (SP) | $68.36\pm 0.42$ | $0.1412\pm 0.0009$ | $0.806\pm 0.006$ |  $147.43\pm$0.23

Table 7: Different fitting results for the GR-$\Lambda$CDM model. The first
five rows correspond to the different non-baseline datasets explored for the
BD-$\Lambda$CDM in Table 6. The last three rows correspond to other scenarios
tested with the BD-$\Lambda$CDM model in Table 10, see Sec. 7.8 for more
details. $H_{0}$ is given in km/s/Mpc.

As mentioned, in our numerical analysis of the data we also wish to consider
the effect of a simple but powerful DDE parametrization, which is the
traditional XCDM[102]. In this very simple framework, the DE is self-conserved
and is associated to some unspecified entity or fluid (called X) which exists
together with ordinary baryonic and cold dark matter, but it does not have any
interaction with them. The energy density of X is simply given by
$\rho_{X}(a)=\rho_{X0}a^{-3(1+w_{0})}$, $\rho_{X0}=\rho_{\Lambda}$ being the
current DE density value and $w_{0}$ the (constant) EoS parameter of such
fluid. More complex parametrizations of the EoS can be considered, for
instance the CPL one [103, 104], in which there is a time evolution of the EoS
itself. However, we have previously shown its incapability to improve the XCDM
performance in solving the two tensions, see [339]. Thus, in this work we
prefer to stay as closer as possible to the standard cosmological model and we
will limit ourselves to analyze the XCDM only. By setting $w_{0}=-1$ we
retrieve the $\Lambda$CDM model with constant $\rho_{\Lambda}$. For
$w_{0}\gtrsim-1$ the XCDM mimics quintessence, whereas for $w_{0}\lesssim-1$
it mimics phantom DE. The fitting results generated from the XCDM on our
datasets are used in our analysis as a figure of merit or benchmark to compare
with the corresponding fitting efficiency of both the BD-$\Lambda$CDM and the
GR-$\Lambda$CDM models. In the next section, we comment on the comparison.

| GR-XCDM
---|---
Parameter | Baseline | Baseline+$H_{0}$
$H_{0}$ (km/s/Mpc) | $67.34^{+0.63}_{-0.66}$ | $68.40^{+0.60}_{-0.62}$
$\omega_{b}$ | $0.02235^{+0.00021}_{-0.00020}$ | $0.02239^{+0.00019}_{-0.00020}$
$\omega_{cdm}$ | $0.11649^{+0.00108}_{-0.00111}$ | $0.11671^{+0.00117}_{-0.00109}$
$\tau$ | $0.053^{+0.006}_{-0.008}$ | $0.051^{+0.005}_{-0.008}$
$n_{s}$ | $0.9709\pm 0.0043$ | $0.9707^{+0.0042}_{-0.0043}$
$\sigma_{8}$ | $0.782^{+0.011}_{-0.010}$ | $0.792\pm 0.011$
$r_{s}$ (Mpc) | $148.05^{+0.32}_{-0.34}$ | $147.95^{+0.33}_{-0.34}$
$w_{0}$ | $-0.956\pm 0.026$ | $-0.991^{+0.026}_{-0.024}$
$\chi^{2}_{min}$ | 2269.88 | 2285.22
$2\ln B$ | $-2.23$ | $-5.21$

Table 8: As in Table 3, but for the XCDM parametrization (within GR).
Motivated by previous works (see e.g. [339]), we have used the (flat) prior
$-1.1<w_{0}<-0.9$.

###  Discussion and extended considerations

In this chapter, we have dealt with Brans-Dicke (BD) theory in extenso. Our
main goal was to assess if BD-gravity can help to smooth out the main two
tensions besetting the usual concordance $\Lambda$CDM model (based on GR): i)
the $H_{0}$-tension (the most acute existing discordance at present), and ii)
the $\sigma_{8}$-tension, which despite not being so sharp it often occurs
that the (many) models in the market dealing with the former tend to seriously
aggravate the latter. As we have explained at the beginning, the ‘golden rule’
to be preserved by the tension solver should be to find a clue on how to
tackle the main discrepancy (on the local $H_{0}$ parameter) while at the same
time to curb the $\sigma_{8}$ one, or at least not to worsen it. We have found
that BD-gravity could be a key paradigm capable of such achievement.
Specifically, we have considered the original BD model with the only addition
of a cosmological constant (CC), and we have performed a comprehensive
analysis in the light of a rich and updated set of observations. These involve
a large variety of experimental inputs of various kinds, such as the long
chain SNIa+$H(z)$+BAO+LSS+CMB of data sources, which we have considered at
different levels and combinations; and tested with the inclusion of other
potentially important factors such as the influence of the bispectrum
component in the structure formation data (apart from the ordinary power
spectrum); and also assessed the impact of gravitational lensing data of
different sorts (Weak and Strong Lensing).
Although BD-gravity is fundamentally different from GR, we have found very
useful to try to pick out possible measurable signs of the new gravitational
paradigm by considering the two frameworks in the (spatially flat) FLRW metric
and compare the versions of the $\Lambda$CDM model resulting in each case,
which we have called BD-$\Lambda$CDM and GR-$\Lambda$CDM, respectively. We
have parametrized the departure of the former from the latter at the
background level (cf. Sec. 7.3) and we have seen that BD-$\Lambda$CDM can
appear in the form of a dynamical dark energy (DDE) version of the
GR-$\Lambda$CDM, in which the vacuum energy density is evolving through a non-
trivial EoS (cf. Fig. 8). We have called it the ‘GR-picture’ of the BD theory.
The resulting effective behaviour is $\Lambda$CDM-like with, however, a mild
time-evolving quasi-vacuum component. In fact, such behaviour is not ‘pure
vacuum’ – which is why we call it quasi-vacuum – despite of the fact that the
original BD-$\Lambda$CDM theory possesses a rigid cosmological constant.
Specifically, using the numerical fitting results of our analysis we find that
such EoS shows up in effective quintessence-like form at more than $3\sigma$
c.l. (this is perfectly appreciable at naked eye in Fig. 8). Our fit to the
data demonstrates that such an effective representation of BD-gravity can be
competitive with the concordance model with a rigid $\Lambda$-term. It may
actually create the fiction that the DE is dynamical when viewed within the GR
framework, whilst it is just a rigid CC in the underlying BD action. The
practical outcome is that the BD approach with a CC definitely helps to smooth
out some of the tensions afflicting the $\Lambda$CDM in a manner very similar
to the Running Vacuum Model, see e.g. [201, 78, 110, 251, 229, 233], and this
success might ultimately reveal the signs of the BD theory. We conclude that
finding traces of vacuum dynamics, accompanied with apparent deviations from
the standard matter conservation law [181] could be the ‘smoking gun’ pointing
to the possibility that the gravity theory sitting behind these effects is not
GR but BD.

Scenarios | ${\sigma}_{8(\rm GR)}$ | ${\sigma}_{8(\rm BD)}$ | ${S}_{8(\rm GR)}$ | ${S}_{8(\rm BD)}$ | $\tilde{S}_{8(\rm BD)}$
---|---|---|---|---|---
Baseline | $0.797^{+0.005}_{-0.006}$ | $0.785\pm 0.013$ | $0.800^{+0.010}_{-0.011}$ | $0.777^{+0.019}_{-0.020}$ | $0.793\pm 0.012$
Baseline+$H_{0}$ | $0.796^{+0.006}_{-0.007}$ | $0.789\pm 0.013$ | $0.793^{+0.011}_{-0.010}$ | $0.758\pm 0.018$ | $0.792\pm 0.013$
Baseline+$H_{0}$+SL | $0.795^{+0.006}_{-0.007}$ | $0.789\pm 0.013$ | $0.789^{+0.011}_{-0.010}$ | $0.753^{+0.017}_{-0.018}$ | $0.791^{+0.012}_{-0.013}$
Spectrum | $0.800^{+0.006}_{-0.007}$ | $0.798\pm 0.014$ | $0.807\pm 0.013$ | $0.793^{+0.021}_{-0.022}$ | $0.805\pm 0.014$
Spectrum+$H_{0}$ | $0.798\pm 0.007$ | $0.801^{+0.015}_{-0.014}$ | $0.794^{+0.012}_{-0.013}$ | $0.773^{+0.019}_{-0.021}$ | $0.802^{+0.013}_{-0.014}$

Table 9: Fitted values for the $\sigma_{8(M)}$ (here M stands for GR or BD)
obtained under different dataset configurations. We also include the derived
values of $S_{8(M)}$ = $\sigma_{8}\sqrt{\Omega_{m}/0.3}$ and the renormalized
$\tilde{S}_{8}$ = $\sigma_{8}\sqrt{\tilde{\Omega}_{m}/0.3}$, with
$\tilde{\Omega}_{m}$ defined as in Sec. 7.2.1. These results correspond to the
main Scenario II of the BD-$\Lambda$CDM model.

####  Alleviating the $H_{0}$-tension

Our analysis of the BD theory with a CC, taken at face value, suggests that
the reason for the enhancement of $H_{0}$ in the BD model is because the
effective gravitational coupling acting at cosmological scales, $G_{\rm
eff}\sim G_{N}/\varphi$, is higher than the one measured on Earth (see Fig.
4). This possibility allows the best-fit current energy densities of all the
species to remain compatible at $\lesssim 1\sigma$ c.l. with the ones obtained
in the GR-$\Lambda$CDM model (cf. e.g. Tables 3-5). Thus, since $\varphi<1$ we
find that the increase of the Hubble parameter is basically due to the
increase of the effective $G$, and there is no need for a strong modification
of the energy densities of the various species filling the Universe. This is a
welcome feature since the measurable cosmological mass parameter in the
BD-$\Lambda$CDM model is, for sufficiently small $\epsilon_{\rm BD}$, not the
usual $\Omega_{m}$, but precisely the tilded one, related to it through
$\tilde{\Omega}_{m}=\Omega_{m}/\varphi$. The latter is about $\sim 8-9\%$
bigger than its standard model counterpart ($\tilde{\Omega}_{m}>\Omega_{m}$)
as it follows from Fig. 4, where we can read off the current value
$\varphi(z=0)\simeq 0.918$. Now, because at the background level it is
possible to write an approximate Friedmann’s equation (7.17) in terms of the
tilded parameters, these are indeed the ones that are actually measured from
SNIa and BAO observations in the BD context262626Recall that for
$\epsilon_{\rm BD}\neq 0$ the tilded parameters $\tilde{\Omega}_{i}$ (which
were originally defined for $\epsilon_{\rm BD}=0$) receive a correction and
become the hatted parameters $\hat{\Omega}_{i}$ introduced in Eq. (7.38).
However, the difference between them is of ${\cal O}(\epsilon_{\rm BD})$, see
Eq. (7.39), and since $|\epsilon_{\rm BD}|\sim\mathcal{O}(10^{-3})$ it can be
ignored.. The differences, however, as we have just pointed out, are not to be
attributed to a change in the physical energy content of matter but to the
fact that $\varphi<1$ throughout the entire cosmic history, as clearly shown
in Fig. 4. Obviously, the measurement of the parameters $\tilde{\Omega}_{i}$
can be performed through the very same data and procedures well accounted for
in the context of the GR-$\Lambda$CDM framework. This explanation is perfectly
consistent with the fact that when the Friedmann’s law is expressed in terms
of the effective $G$, as indicated in Eq. (7.15), the local value of the
Hubble parameter $H_{0}$ becomes bigger owing to $G_{\rm eff}=G_{N}/\varphi$
being bigger than $G_{N}$. Thus, when we compare the early and local
measurements of $H_{0}$ we do not meet any anomaly in this approach.
We also recall at this point that there is no correction from $\omega_{\rm
BD}$ on the effective coupling $G_{\rm eff}$, Eq. (7.80), in the local domain.
This is because in our context $\omega_{\rm BD}$ appears as being very large
owing to the assumed screening of the BD-field caused by the clustered matter
(cf. Sec. 7.5). From Fig. 4 and Table 3 we find that the BD model leads to a
value of $H_{0}$ a factor $G_{\rm eff}^{1/2}(z=0)/G_{N}^{1/2}\sim
1/\varphi^{1/2}(z=0)=1/\sqrt{0.918}$, i.e. $\sim 4.5\%$, bigger than the one
inferred from the CMB in the GR-$\Lambda$CDM model, in which $G_{\rm
eff}=G_{N}\,(\forall{z})$. It is reassuring to realize that such a
‘renormalization factor’ can enhance the low Planck 2018 CMB measurement of
the Hubble parameter (viz. $H_{0}=67.4\pm 0.5$ km/s/Mpc [31]) up to the range
of $70-71$km/s/Mpc (cf. e.g. Tables 3-6), hence much closer to the local
measurements. For example, SH0ES yields $H_{0}=(73.5\pm 1.4)$ km/s/Mpc [84];
and when the latter is combined with Strong Lensing data from the H$0$LICOW
collab.[97] it leads to $H_{0}=(73.42\pm 1.09)$ km/s/Mpc. This combined value
is $5\sigma$ at odds with the Planck 2018 measurement, a serious tension.

| BD-$\Lambda$CDM (Scenario III: Cassini-constrained)
---|---
Parameter | B+$H_{0}$ (No LSS) | B+$H_{0}$ | Dataset [335] | Dataset [335] \+ LSS (SP)
$H_{0}$ (km/s/Mpc) | $70.99^{+0.94}_{-0.97}$ | $70.80^{+0.81}_{-0.91}$ | $70.01^{+0.86}_{-0.92}$ | $70.03^{+0.90}_{-0.88}$
$\omega_{b}$ | $0.02257\pm 0.00021$ | $0.02256^{+0.00019}_{-0.00020}$ | $0.02271\pm 0.00016$ | $0.02272^{+0.00015}_{-0.00016}$
$\omega_{cdm}$ | $0.11839^{+0.00093}_{-0.00094}$ | $0.11748\pm 0.00089$ | $0.11885^{+0.00092}_{-0.00095}$ | $0.11827^{+0.00089}_{-0.00093}$
$\tau$ | $0.057^{+0.007}_{-0.008}$ | $0.050^{+0.004}_{-0.008}$ | $0.061^{+0.007}_{-0.008}$ | $0.058^{+0.006}_{-0.008}$
$n_{s}$ | $0.9824^{+0.0057}_{-0.0058}$ | $0.9811^{+0.0051}_{-0.0052}$ | $0.9783^{+0.0052}_{-0.0059}$ | $0.9701^{+0.0056}_{-0.0054}$
$\sigma_{8}$ | $0.815^{+0.008}_{-0.009}$ | $0.804^{+0.006}_{-0.007}$ | $0.817\pm 0.007$ | $0.812^{+0.006}_{-0.007}$
$r_{s}$ (Mpc) | $142.14^{+1.91}_{-1.72}$ | $143.31^{+1.72}_{-1.63}$ | $143.58^{+1.62}_{-1.55}$ | $144.10^{+1.62}_{-1.52}$
$\epsilon_{\rm BD}$ | $-0.00002\pm 0.00002$ | $-0.00002\pm 0.00002$ | $-0.00002\pm 0.00002$ | $-0.00002\pm 0.00002$
$\varphi_{ini}$ | $0.933\pm 0.021$ | $0.944\pm 0.020$ | $0.955^{+0.018}_{-0.019}$ | $0.960^{+0.020}_{-0.018}$
$\varphi(0)$ | $0.933^{+0.020}_{-0.021}$ | $0.944^{+0.019}_{-0.020}$ | $0.955^{+0.018}_{-0.019}$ | $0.960^{+0.020}_{-0.017}$
$w_{\rm eff}(0)$ | $-0.972\pm 0.009$ | $-0.977\pm 0.008$ | $-0.981^{+0.008}_{-0.007}$ | $-0.983\pm 0.008$
$\dot{G}(0)/G(0)(10^{-13}yr^{-1})$ | $0.025^{+0.025}_{-0.026}$ | $0.026^{+0.027}_{-0.028}$ | $0.023^{+0.026}_{-0.027}$ | $0.020\pm 0.026$
$\chi^{2}_{min}$ | 2256.14 | 2278.34 | 2797.44 | 2812.68
$2\ln B$ | +9.03 | $+5.21$ | +3.45 | +2.21

Table 10: Fitting results for the BD-$\Lambda$CDM, in the context of the BD-
Scenario III explained in Sec. 7.5 under different datasets. As characteristic
of Scenario III, in all of these datasets the Cassini constraint on the post-
Newtonian parameter $\gamma^{\rm PN}$ has been imposed [317]. In the first two
fitting columns we use the Baseline+$H_{0}$ dataset described in Sec. 7.6.
However, we exclude the LSS data in the first column while it is kept in the
second. In the third and fourth fitting columns we report on the results
obtained using the very same dataset as in Ref. [335], just to ease the
comparison between the BD-$\Lambda$CDM and the variable $G$ model studied in
that reference (cf. their Table 1). This dataset includes the Planck 2018
TTTEEE+lensing likelihood [31], BAO data from [153, 154, 278] and the SH0ES
prior from [321]. In the last fitting column, however, we add the LSS data to
the previous set but with no bispectrum (cf. Table 2 and Sec. 7.6). The
corresponding results for the GR-$\Lambda$CDM can be found in Tables 3 and 7.
Figure 10: Triangular matrix containing the two-dimensional marginalized
distributions for some relevant combinations of parameters in the
BD-$\Lambda$CDM model (at $1\sigma$ and $2\sigma$ c.l.), together with the
corresponding one-dimensional marginalized likelihoods for each of them.
$H_{0}$ is expressed in km/s/Mpc, and $r_{s}$ in Mpc. See Tables 3 and 4 for
the numerical fitting results.

On the other hand, if we compare e.g. the fitting value predicted within the
BD-$\Lambda$CDM model from Table 3 (namely $H_{0}=70.83^{+0.92}_{-0.95}$
km/s/Mpc) with the aforementioned SH0ES determination, we can see that the
difference is of only $1.58\sigma$. If we next compare our fitting result from
Table 4 ($H_{0}=71.30^{+0.80}_{-0.84}$ km/s/Mpc), which incorporates the
H$0$LICOW Strong Lensing data in the fit as well, with the combined SH0ES and
H$0$LICOW result (viz. the one which is in $5\sigma$ tension with the CMB
value) we obtain once more an inconspicuous tension of only $1.55\sigma$. In
either case it is far away from any perturbing discrepancy. In fact, no
discrepancy which is not reaching a significance of at least $3\sigma$ can be
considered sufficiently worrisome.

####  Alleviating the $\sigma_{8}$-tension

Furthermore, the smoothing of the tension applies to the $\sigma_{8}$
parameter as well, with the result that it essentially disappears within a
similar level of inconspicuousness. In fact, values such as
$\sigma_{8}=0.789\pm 0.013$ and $\tilde{S}_{8}=0.792\pm 0.013$ (obtained
within the Baseline$+H_{0}$ dataset, see Table 9) are in good agreement with
weak gravitational lensing observations derived from shear data (cf. the WL
data block mentioned in Sec. 7.6). Let us take the value by Joudaki et al.
2018 of the combined observable $S_{8}=0.742\pm 0.035$, for example, obtained
by KiDS-450, 2dFLenS and BOSS collaborations from a joint analysis of weak
gravitational lensing tomography and overlapping redshift-space galaxy
clustering [245]. These observations can be compared with our prediction for
${S}_{8}=\sigma_{8}\sqrt{\Omega_{m}/0.3}$ and with the ‘renormalized’ form of
that quantity within the BD-$\Lambda$CDM model, namely
$\tilde{S}_{8}=\sigma_{8}\sqrt{\tilde{\Omega}_{m}/0.3}$, which depends on the
modified cosmological parameter $\tilde{\Omega}_{m}$, which is slightly
higher, recall our Eq. (7.18)272727Although we could use the hatted parameter
$\hat{S}_{8}=\sigma_{8}\sqrt{\hat{\Omega}_{m}/0.3}$, instead of
$\tilde{S}_{8}$, we have already pointed out that the difference between
$\hat{\Omega}_{m}$ and $\tilde{\Omega}_{m}$ is negligible for $|\epsilon_{\rm
BD}|\sim\mathcal{O}(10^{-3})$, and so is between $\hat{S}_{8}$ and
$\tilde{S}_{8}$.. Both ${S}_{8}$ and $\tilde{S}_{8}$ are displayed together in
Table 9 for the main scenarios, also in company with $\sigma_{8}$ values for
the GR and BD models. Differences of the mentioned experimental measurements
with respect to e.g. our prediction for the Baseline$+H_{0}$ dataset, are at
the level of $0.5\sigma-1.3\sigma$ depending on whether we use $S_{8}$ or
$\tilde{S}_{8}$, whence statistically irrelevant in any case. More details can
be appraised on some of these observable and their correlation with $H_{0}$ in
Figs. 10 and 11, on which we shall further comment later on.

Figure 11: Constraints obtained for $\sigma_{8}$ and $\tilde{S}_{8}$ versus
$H_{0}$ (in km/s/Mpc) from the fitting analyses of the GR- and BD-$\Lambda$CDM
models, and the GR-XCDM parametrization. We show both, the contour lines in
the corresponding planes of parameter space, and the associated marginalized
one-dimensional posteriors. The centering of the parameters in the ranges
$\sigma_{8}<0.80$ and $H_{0}\gtrsim 71$ km/s/Mpc is a clear sign of the
smoothening of the $\sigma_{8}$-tension and, more conspicuously, of the
$H_{0}$ one within the BD-$\Lambda$CDM model. We can also see that while a
simple XCDM parametrization for the DE can help to diminish $\sigma_{8}$ as
compared to the concordance model, it is however unable to improve the $H_{0}$
tension, which is kept at a similar level as within the concordance model.

####  Comparing different scenarios

We have also tested the performance of the BD- and GR-$\Lambda$CDM models when
we include the CMB high-$l$ (TE+EE) polarization data from Planck 2018, with
and without the CMB lensing, in combination with the baseline dataset and the
SH0ES prior on $H_{0}$ (cf. Tables 6 and 7). The values of the Hubble
parameter in these cases are a little bit lower than when we consider only the
temperature and low-$l$ polarization (TT+lowE) CMB data, but the tension is
nevertheless significantly reduced, being now of only $\sim 2.2\sigma$ c.l.,
whereas in the GR-$\Lambda$CDM model it is kept at the $\sim 3.5\sigma$ level.
The values of $\sigma_{8}$ are still low, $\sim 0.80-0.81$. The information
criteria in these cases, though, have no preference for any of the two models,
they are not conclusive.
We also examine the results that are obtained when we do not include in our
fitting analyses any of the data sources that trigger the tensions. We
consider here the CMB, BAO and SNIa datasets (denoted as CMB+BAO+SNIa in
Tables 6 and 7), but exclude the use of the SH0ES prior on $H_{0}$ and the LSS
data. As expected, the evidence for the BD model decreases since now we do not
give to it the chance of showing its power. Even though the description of the
data improves, it is not enough to compensate for the penalty received owing
to the use of the two additional parameters $(\epsilon_{\rm
BD},\varphi_{ini})$, and in this case we read $2\ln B=-3.00$ (from Table 6).
Thus, there is here a marginal preference for the GR scenario, but as
previously mentioned, this is completely normal, since we are removing
precisely the data sources whose correct description demands for new physics.
Even so, the $H_{0}$-tension is again remarkably reduced from $3.8\sigma$ in
GR-$\Lambda$CDM to only $2.4\sigma$ in the BD-$\Lambda$CDM. The respective
values of $\sigma_{8}$ remain compatible with $0.80$ within $\sim 1\sigma$.
It is also interesting to compare the results that we have obtained within the
BD framework with other approaches in the literature, in which the variation
of $G$ is dealt with as a small departure from GR, namely in a context where
the action still contains a large mass scale $M$ near the Planck mass $m_{\rm
Pl}$, together with some scalar field which parametrizes the deviations from
it. This is of course fundamentally different from the BD paradigm but it
bears relation owing to the variation of the effective $G$, and it has also
been used to try to smooth out the tensions. However, as already mentioned at
the beginning, and also in Sec. 7.2, it is not easy at all for a given model
to fulfill the ‘golden rule’, i.e. to loosen the two tensions at a time, or
just to alleviate one of them without worsening the other. Different proposals
have appeared in the market trying to curb the $H_{0}$-tension, e.g. the so-
called early dark energy models [333, 360, 361], and the model with variable
$G$ recently considered in [335, 362, 363]. Although the physical mechanism of
the EDE and the aforementioned models with variable $G$ is of course very
different, their aim is pretty similar. They reduce the sound horizon $r_{s}$
at recombination in order to force the increase of the Hubble function in the
late Universe. This allows them to generate larger values of $H_{0}$ in order
to keep a good fit to the CMB and BAO data, but this happens only at the
expense of increasing the tension in $\sigma_{8}$, since they do not have any
compensation mechanism able to keep the structure formation in the late
Universe at low enough levels. Some of these models appear not to be
particularly disfavored notwithstanding. But this is simply because they did
not use LSS data in their fits, i.e. they did not put their models to the test
of structure formation and for this reason they have more margin to adjust the
remaining observable without getting any statistical punishment. So the fact
that the significant increase of $\sigma_{8}$ that they find is not
statistically penalized is precisely because they do not use LSS data, such as
e.g. those on the observable $f(z)\sigma_{8}(z)$ displayed in our Table 2. In
this respect, EDE cosmologies are an example; they seem to be unable to
alleviate the $H_{0}$-tension when LSS data are taken into account, as shown
in [334].

####  Imposing the Cassini constraint

To further illustrate the capability of the BD-$\Lambda$CDM model to fit the
data under more severe conditions, in Table 10 we consider four possible
settings to fit our Cassini-constrained BD-Scenario III defined in Sec. 7.5.
Recall that this BD scenario involves the stringent Cassini bound on the post-
Newtonian parameter $\gamma^{\rm PN}$ [317], which we have discussed in Sec.
7.1. The first two fitting columns of Table 10 correspond to our usual
Baseline$+H_{0}$ dataset, in one case (first fitting column in that table) we
omit the LSS data, whereas in the second column we restore it. In this way we
can check the effect of the structure formation data on the goodness of the
fit. The comparison between the results presented in these two columns shows,
first and foremost, that the Cassini bound does not have a drastically
nullifying effect, namely it does not render the BD-$\Lambda$CDM model
irrelevant to the extent of making it indistinguishable from GR-$\Lambda$CDM,
not at all, since the quality of the fits is still fairly high (confer the
Bayes factors in the last row). In truth, the fit quality is still comparable
to that of the Baseline$+H_{0}$ scenario (cf. Table 3). However, the
description of the LSS data is naturally poorer since $\epsilon_{\rm BD}$ is
smaller and the model cannot handle so well the features of the structure
formation epoch, thus yielding slightly higher values of $\sigma_{8}$. Second,
the fact that the scenario without LSS furnishes a higher Bayes factor just
exemplifies the aforementioned circumstance that when cosmological models are
tested without using this kind of data the results may in fact not be
sufficiently reliable. When the LSS data enter the fit (third fitting column
in that table), we observe, interestingly enough, that the BD-$\Lambda$CDM
model is still able to keep $H_{0}$ in the safe range, it does improve the
value of $\sigma_{8}$ as well (i.e. it becomes lower) and, on top of that, it
carries a (‘smoking gun’) signal of almost $2.9\sigma$ c.l. – encoded in the
value of $w_{\rm eff}$ – pointing to quintessence-like behaviour. Overall it
is quite encouraging since it shows that the Cassini bound does not
exceedingly hamper the BD-$\Lambda$CDM model capabilities. Such bound
constraints the time evolution of $\varphi$ (because $|\epsilon_{\rm BD}|$ is
forced to be much smaller) but it does not preclude $\varphi$ from choosing a
suitable value in compensation (cf. BD-Scenario III in Sec. 7.5).
In the last two columns of Table 10, we can further check what are the changes
in the previous fitting results when we use a more restricted dataset, e.g.
the one used in Ref. [335], in which the Cassini bound is also implemented.
These authors study a model which represents a modification of GR through an
effective $G\sim 1/M^{2}$, with $M$ a mass very near the Planck mass, $m_{\rm
Pl}$, which is allowed to change slowly through a scalar field $\phi$ as
follows: $M^{2}\to M^{2}+\beta\phi^{2}$, where $\beta$ is a small
(dimensionless) parameter. The authors assume that the Cassini bound on the
post-Newtonian parameter $\gamma^{\rm PN}$ [317] is in force (see Sec. 7.6 for
details). However, they do not consider LSS data (only CMB and BAO). We may
compare the results they obtain within that variable $G$ model (cf. their
Table 1) with those we obtain within the BD-Scenario III under the very same
dataset as these authors. The results are displayed in the third fitting
column of Table 10. We obtain $H_{0}=(70.01^{+0.86}_{-0.92})$ km/s/Mpc and
$\sigma_{8}=0.817\pm 0.007$, whereas they obtain
$H_{0}=(69.2^{+0.62}_{-0.75})$ km/s/Mpc and
$\sigma_{8}=0.843^{+0.015}_{-0.024}$. Clearly, the BD-$\Lambda$CDM is able to
produce larger central values of $H_{0}$ and lower values of $\sigma_{8}$,
even under the Cassini bound, although the differences are within errors. The
value of $2\ln\,B$ lies around $+3.5$ and hence points to a mild positive
evidence in favor of the BD model. This is consistent with the associated
deviation we find of $G(0)$ from $G_{N}$ at $2.43\sigma$ c.l., and with a
signal of effective quintessence at $2.53\sigma$ c.l. within the GR-picture.
Let us now consider what is obtained if we add up the LSS data to this same
BD-Scenario III, still with the restricted dataset of Ref. [335]. As expected,
the inclusion of the structure formation data pushes the value of
$\sigma_{8}=0.812^{+0.006}_{-0.007}$ down as compared to their absence
($\sigma_{8}=0.817\pm 0.007$). This is the most remarkable difference between
the two cases, as one cannot appreciate significant changes in the other
parameters. Something very similar happens when we compare the values of
$\sigma_{8}$ of the first two fitting columns of Table 10, in which we
consider the B+$H_{0}$ dataset without LSS data (first fitting column) and
with LSS data (second fitting column). The relative improvement w.r.t. the
model of [335] is therefore greater in the presence of LSS data, whose use has
been omitted in that reference. In that variable $G$ model one finds a larger
value of $H(z)$ at recombination thanks to the larger values of $G$ in that
epoch, but at present $G$ is forced to be almost equal to $G_{N}$ (being the
differences not relevant for cosmology). This means that: (i) in that model
$G$ decreases with time, which leads to a kind of effective $\epsilon_{\rm
BD}>0$; (ii) the model cannot increase $H_{0}$ with a large value of $G(z=0)$.
Both facts do limit significantly the effectiveness of the model in loosening
the tensions. In the BD-$\Lambda$CDM model under consideration, instead, we
find that $G$ has to be $\sim 8-9\%$ larger than $G_{N}$ not only in the pre-
recombination Universe, but also at present, and this allows to reduce
significantly the $H_{0}$-tension. Moreover, we find that a mild increase of
the cosmological $G$ with the expansion leads also to an alleviation of the
$\sigma_{8}$-tension.
Under all of the datasets studied in Table 10 we obtain central values of
$|\epsilon_{\rm BD}|\sim\mathcal{O}(10^{-5})$, which are compatible with 0 at
$1\sigma$. Notice that this value is just of order of the Cassini bound on
$\epsilon_{\rm BD}$, as could be expected. Notwithstanding, and remarkably
enough, the stringent bound imposed by the Cassini constraint, which enforces
$\epsilon_{\rm BD}$ to remain two orders of magnitude lower than in the main
Baseline scenarios, is nevertheless insufficient to wipe out the positive
effects from the BD-$\Lambda$CDM model. They are still capable to emerge with
a sizeable part of the genuine BD signal. This is, as anticipated in Sec. 7.5,
mainly due to the fact that the Cassini bound cannot restrict the value of the
BD field $\varphi$, only its time evolution.

####  More on alleviating simultaneously the two tensions

A few additional comments on our results concerning the parameters $H_{0}$ and
$\sigma_{8}$ are now in order. Their overall impact can be better assessed by
examining the triangular matrix of fitted contours involving all the main
parameters, as shown in Fig. 10, in which we offer the numerical results of
several superimposed analyses based on different datasets, all of them within
Baseline scenarios. We project the contour lines in the corresponding planes
of parameter space, and show the associated marginalized one-dimensional
posteriors. The fitted value of the EoS parameter at $z=0$ shown there,
$w_{\rm eff}(0)$, is to be understood, of course, as a derived parameter from
the prime ones of the fit, but we include this information along with the
other parameters in order to further display the significance of the obtained
signal: $\gtrsim 3\sigma$ quintessence-like behaviour. Such signal, therefore,
mimics ‘GR+ DDE’ and hints at something beyond pure GR-$\Lambda$CDM. What we
find in our study is that such time-evolving DE behaviour is actually of
quasi-vacuum type and appears as a kind of signature of the underlying BD
theory282828As we recall in Appendix C, such kind of dynamical behaviour of
the vacuum is characteristic of the Running Vacuum Model (RVM), a version of
the $\Lambda$CDM in which the vacuum energy density is not just a constant but
involves also a dynamical term $\sim H^{2}$. The description of the
BD-$\Lambda$CDM model in the GR-picture mimics a behaviour of this sort [332,
331]..
Figure 10 provides a truly panoramic and graphical view of our main fitting
results, and from where we can comfortably judge the impact of the BD
framework for describing the overall cosmological data. It is fair to say that
it appears at a level highly competitive with GR – in fact, superior to it. In
all datasets involving the local $H_{0}$ input in the fit analyses, the
improvement is substantial and manifest. Let us stand out only three of the
entries in that graphical matrix: i) for the parameters $(\sigma_{8},H_{0})$,
all the contours in the main dataset scenarios are centered around values of
$\sigma_{8}<0.80$ and $H_{0}\gtrsim 71$ km/s/Mpc, which are the coveted ranges
for every model aiming at smoothing the two tensions at a time; ii) for the