pair $(H_{0},r_{s})$, the contours are centered around the same range of
relevant $H_{0}$ values as before, and also around values of the comoving
sound horizon (at the baryon drag epoch) $r_{s}\lesssim 142$ Mpc, these being
significantly smaller than those of the concordance $\Lambda$CDM (cf. Table 7)
and hence consistent with larger values of the expansion rate at that epoch;
iii) and for $(w_{\rm eff}(0),H_{0})$, the relevant range $H_{0}\gtrsim 71$
km/s/Mpc is once more picked out, together with an effective quintessence
signal $w_{\rm eff}(0)>-1$ at more than $3\sigma$ (specifically $3.45\sigma$
for the scenario of Table 4, in which strong lensing data are included in the
fit).
It is also interesting to focus once more our attention on Fig. 11, where we
provide devoted contours involving both $\tilde{S}_{8}$ and $\sigma_{8}$
versus $H_{0}$. On top of the observations we have previously made on these
observable, we can compare here our basic dataset scenarios for the
BD-$\Lambda$CDM model with the yield of a simple XCDM parametrization of the
DDE. In previous studies we had already shown that such parametrization can
help to deal with the $\sigma_{8}$ tension [339]. Nonetheless, as we can see
here, it proves completely impotent for solving or minimally helping to
alleviate the $H_{0}$-tension since the values predicted for this parameter
stay as low as in the concordance model. This shows, once more, that in order
to address a possible solution to the two tensions simultaneously, it is not
enough to have just some form of dynamics in the DE sector; one really needs a
truly specific one, e.g. the one provided (in an effective way) by the
BD-$\Lambda$CDM model.

####  Predicted relative variation of the effective gravitational strength

In the context of the BD framework it is imperative, in fact mandatory, to
discuss the current values of the relative variation of the effective
gravitational strength, viz. of $\dot{G}(0)/G(0)$, which follow from our
fitting analyses (see the main Tables 3-5). The possible time evolution of
that quantity hinges directly on $\epsilon_{\rm BD}$, of course, since the
latter is the parameter that controls the (cosmological) evolution of the
gravitational coupling in the BD theory. It is easy to see from Eq. (7.1) that
$\dot{G}(0)/G(0)=-\dot{\varphi}(0)/\varphi(0)\simeq-\epsilon_{\rm BD}H_{0}$,
where we use the fact that $\varphi\sim a^{\epsilon_{\rm BD}}$ in the matter-
dominated epoch (cf. Appendix C). Recalling that $H_{0}\simeq 7\times
10^{-11}$ yr-1 (for $h\simeq 0.70$), we find
$\dot{G}(0)/G(0)\simeq-\epsilon_{\rm BD}\cdot 10^{-10}\,yr^{-1}$. Under our
main BD-Scenario II (cf. Sec. 7.5) we obtain values for $\dot{G}(0)/G(0)$ of
order $\mathcal{O}(10^{-13})\,yr^{-1}$ (and positive), just because
$\epsilon_{\rm BD}\sim\mathcal{O}(10^{-3})$ (and negative). Being
$\dot{G}(0)/G(0)>0$ it means that the effective gravitational coupling
obtained by our global cosmological fit increases with the expansion, and
hence it was smaller in the past. This suggests that the sign $\epsilon_{\rm
BD}<0$, which is directly picked out by the data, prefers a kind of
asymptotically free behaviour for the gravitational coupling since the epochs
in the past are more energetic, in fact characterized by larger values of $H$
(with natural dimension of energy). The central values show a mild time
variation at present, at a level of $1.3\sigma$, both for $\epsilon_{\rm BD}$
and $\dot{G}(0)/G(0)$, when the bispectrum data from BOSS is also included in
the analysis (cf. Tables 3 and 4). Such departure goes below $1\sigma$ level
when only the spectrum is considered (see Table 5). In the context of the BD-
Scenario III, in which $\epsilon_{\rm BD}$ is very tightly constrained by the
Cassini bound [317], namely at a level of $\mathcal{O}(10^{-5})$, we find
$\dot{G}(0)/G(0)\sim 10^{-15}\,yr^{-1}$, which is compatible with 0 at
$1\sigma$. All that said, we should emphasize once more that the fitting
values that we obtain for $\dot{G}(0)/G(0)$ refer to the cosmological time
variation of $G$ and, therefore, cannot be directly compared with constraints
existing in the literature based on strict local gravity measurements, such as
e.g. those from the lunar laser ranging experiment – $\dot{G}(0)/G(0)=(2\pm
7)\cdot 10^{-13}yr^{-1}$ [364] – (see e.g. the review [168] for a detailed
presentation of many other local constraints on $\dot{G}(0)/G(0)$). Even
though this bound turns out to be preserved within our analysis, it is not in
force at the cosmological level provided an screening mechanism acting at
these scales is assumed, as in our case. Thus, the local measurements have no
bearing a priori on the BD-$\Lambda$CDM cosmology. The opposite may not be
true, for despite the fact that the values reported in our tables are model-
dependent, they prove to be quite efficient and show that the cosmological
observations can compete in precision with the local measurements.

Figure 12: CMB temperature power spectrum for the GR-$\Lambda$CDM (in black) and BD-$\Lambda$CDM (in orange), obtained from the fitting results within the Baseline+$H_{0}$ dataset (cf. Table 3 and Sec. 7.6). We plot the central curves together with the corresponding $1\sigma$ bands. In the inner plot we zoom in the multipole range $l\in[0,30]$ and include the Planck 2018 [31] data error bars (in green). At low multipoles ($l\lesssim 30$) the BD-$\Lambda$CDM model produces less power than the concordance GR-$\Lambda$CDM model owing to the suppression of the Integrated Sachs-Wolfe effect that we have discussed in Sec. 7.2. The differences are $\gtrsim 1\sigma$, and allow to soften a well-known low-multipole CMB anomaly. See the main text for further discussion. | Baseline | Baseline+$H_{0}$
---|---|---
Parameter | GR-$\Lambda$CDM | BD-$\Lambda$CDM | GR-$\Lambda$CDM | BD-$\Lambda$CDM
$H_{0}$ (km/s/Mpc) | $67.75^{+0.46}_{-0.48}$ | $68.86^{+1.26}_{-1.22}$ | $68.35^{+0.49}_{-0.46}$ | $70.81^{+0.95}_{-0.92}$
$\omega_{b}$ | $0.02231^{+0.00018}_{-0.00020}$ | $0.02214\pm 0.00031$ | $0.02239\pm 0.00019$ | $0.02239^{+0.00029}_{-0.00032}$
$\omega_{cdm}$ | $0.11655^{+0.00110}_{-0.00111}$ | $0.11860^{+0.00194}_{-0.00197}$ | $0.11625^{+0.00099}_{-0.00106}$ | $0.11821^{+0.00214}_{-0.00200}$
$\tau$ | $0.053^{+0.006}_{-0.008}$ | $0.054^{+0.006}_{-0.008}$ | $0.053^{+0.006}_{-0.008}$ | $0.054^{+0.006}_{-0.009}$
$n_{s}$ | $0.9706^{+0.0041}_{-0.0040}$ | $0.9691^{+0.0104}_{-0.0093}$ | $0.9717^{+0.0041}_{-0.0042}$ | $0.9791^{+0.0097}_{-0.0090}$
$\sigma_{8}$ | $0.770^{+0.016}_{-0.018}$ | $0.759^{+0.0184}_{-0.0163}$ | $0.780^{+0.018}_{-0.015}$ | $0.762^{+0.021}_{-0.018}$
$r_{s}$ (Mpc) | $148.03\pm 0.33$ | $146.16^{+2.40}_{-2.58}$ | $148.04^{+0.32}_{-0.34}$ | $142.81^{+1.91}_{-2.03}$
$m_{\nu}$ (eV) | $0.161^{+0.058}_{-0.059}$ | $0.409^{+0.139}_{-0.198}$ | $0.118^{+0.053}_{-0.058}$ | $0.409^{+0.183}_{-0.228}$
$\epsilon_{\rm BD}$ | - | $0.00459^{+0.00316}_{-0.00319}$ | - | $0.00433^{+0.00395}_{-0.00336}$
$\varphi_{ini}$ | - | $0.979^{+0.028}_{-0.032}$ | - | $0.938^{+0.023}_{-0.024}$
$\varphi(0)$ | - | $1.016^{+0.041}_{-0.049}$ | - | $0.972^{+0.042}_{-0.043}$
$w_{\rm eff}(0)$ | - | $-1.005^{+0.021}_{-0.017}$ | - | $-0.986^{+0.017}_{-0.016}$
$\dot{G}(0)/G(0)(10^{-13}yr^{-1})$ | - | $-5.048^{+3.485}_{-3.434}$ | - | $-4.915^{+3.749}_{-4.501}$
$\chi^{2}_{min}$ | 2270.54 | 2268.44 | 2285.02 | 2274.76
$2\ln B$ | - | -0.12 | - | +9.34
$\Delta{\rm DIC}$ | - | -0.81 | - | +8.22
Table 11: As in Table 3, but here we allow the variation of the mass of the
massive neutrino ($m_{\nu}$) in the Monte Carlo routine, instead of setting it
to $0.06$ eV. The other two neutrinos remain massless. We have used the same
conservative prior range for $m_{\nu}\in[0,1]$ eV in both, the GR- and
BD-$\Lambda$CDM models.

####  One more bonus: suppressing the power at low multipoles

An additional bonus from the BD cosmology is worth mentioning before we close
this lengthy study. It is found in the description of the CMB temperature
anisotropies. As we have discussed in Sec. 7.2, the BD-$\Lambda$CDM model is,
in principle, able to suppress the power at low multipoles ($l\lesssim 30$),
thereby softening one of the so-called CMB anomalies that are encountered in
the context of the GR-$\Lambda$CDM model. This is basically due to the low
values of $\varphi<1$ preferred by the data, which in turn produce a
suppression of the Integrated Sachs-Wolfe (ISW) effect [323, 324]. We have
confirmed that this suppression actually occurs for the best-fit values of the
parameters in our analysis, cf. Fig. 12. The aforementioned anomaly is not
very severe, since the power at low multipoles is affected by a large cosmic
variance and cannot be measured very precisely. Nevertheless, it is a subtle
anomaly which has been there unaccounted for a long time and could not be
improved in a consistent way: that is to say, usually models ameliorating the
low tail of the spectrum do spoil the high part of it. However, here the
suppression of power with respect to the GR-$\Lambda$CDM at $l\lesssim 30$ is
fully consistent and is another very welcome feature of the BD-$\Lambda$CDM
model, which is not easy at all to attain. It is interesting to mention that
the ISW effect can be probed by cross correlating the CMB temperature maps
with the LSS data, e.g. with foreground galaxies number counts, especially if
using future surveys which should have much smaller uncertainties. This can be
a useful probe for DE and a possible distinctive signature for DDE theories
[365, 366] 292929We thank L. Pogosian for interesting comments along these
lines.. The upshot of our investigation is that once more the ‘golden rule’
mentioned at the beginning is preserved here and the curing effects from the
BD-$\Lambda$CDM stay aligned: the three tensions of the GR-$\Lambda$CDM
($H_{0}$, $\sigma_{8}$ and the exceeding CMB power at low multipoles) can be
improved at a time.

####  The effect of massive neutrinos

Finally, it is worth assessing the impact of leaving the sum of the three
neutrino masses ($M_{\nu}\equiv\sum m_{\nu}$) as a free parameter, rather than
fixing it. We model this scenario in two different ways:

  * •

Considering one massive ($m_{\nu}$) and two massless neutrinos, so
$M_{\nu}=m_{\nu}$. Here we try to mimic the physics encountered when the
neutrino masses follow the normal hierarchy, in which one of the neutrinos is
much heavier than the other two, for reasonable values of $M_{\nu}$.

  * •

Considering three massive neutrinos of equal mass $m_{\nu}$, such that
$M_{\nu}=3m_{\nu}$. This is the degenerated case, utilized also in the
analysis by Planck 2018 [31] (cf. Sec. 7.5.1 therein).

In both cases, the BD-$\Lambda$CDM has one additional parameter w.r.t. the
same model with fixed mass $M_{\nu}=m_{\nu}=0.06$ eV that we have discussed in
the previous sections, and three more than the vanilla $\Lambda$CDM model:
$(\epsilon_{\rm BD},\varphi_{ini},M_{\nu})$. Upon taking into account the
constraints obtained from neutrino oscillation experiments on the mass-squared
splittings through the corresponding likelihoods we would perhaps find more
precise bounds on the sum of the neutrino masses [367], but here we want to
carry out a more qualitative analysis to study how massive neutrinos may
impact on our results, considering only cosmological data. The results of the
first mass scenario are shown in Table 11. Those for the second scenario are
not tabulated, but are commented below.
As can be seen from Table 11, the effect of having a neutrino with an
adjustable mass picked out by the fitting process is non-negligible. It
produces a significant lowering of the value of $\sigma_{8}$ while preserving
the value of $H_{0}$ at a level comparable to the previous tables in which the
neutrino had a fixed mass of $0.06$ eV, therefore preserving what we have
called the ‘golden rule’. Let us note, however, that the sign of
$\epsilon_{\rm BD}$ has changed now with respect to the situation with a fixed
light mass, it is no longer negative but positive and implies a value of the
BD parameter of $\omega_{\rm BD}\simeq 230$. Last but not least, the fitted
value of the neutrino mass is $m_{\nu}=0.409$ eV, which is significantly
higher than the upper bound placed by Planck 2018 under the combination
TTTEEE+lowE+lensing+BAO: $M_{\nu}\equiv\sum m_{\nu}<0.120$ eV ($95\%$ c.l.)
[31].
Similar results are obtained with the second mass scenario mentioned above.
Let us summarize them. For the Baseline+$H_{0}$ dataset, the common fitted
mass value obtained for the three neutrinos is
$m_{\nu}=0.120^{+0.054}_{-0.068}$ eV, which means that $M_{\nu}\simeq 0.360$
eV, slightly lower than in the first scenario. The corresponding BD parameter
reads comparable, $\omega_{\rm BD}\simeq 261$, and the values for $\sigma_{8}$
and $H_{0}$ are also very close to the previous case, so the two scenarios
share similar advantages. The fact that the Planck 2018 upper limit for
$M_{\nu}$ is overshooted in both neutrino mass scenarios does not necessarily
exclude them, as the limits on $M_{\nu}$ are model-dependent, see e.g. [367]
and [363]. In particular, Planck 2018 obviously used GR-$\Lambda$CDM. We can
check in Table 11 that the neutrino mass values for BD-$\Lambda$CDM are
substantially different. The results that we have obtain for GR-$\Lambda$CDM
are fully compatible with those from Planck 2018. We conclude that the
influence of neutrino masses on the BD-$\Lambda$CDM fitting results is
potentially significant, but it cannot be settled at this point. The subject
obviously deserves further consideration in the future.

###  Conclusions

To summarize, we have presented a rather comprehensive work on the current
status of the Brans-Dicke theory with a cosmological constant in the light of
the modern observations. Such framework constitutes a new version of the
concordance $\Lambda$CDM model in the context of a distinct gravity paradigm,
in which the gravitational constant is no longer a fundamental constant of
Nature but a dynamical field. We have called this framework BD-$\Lambda$CDM
model to distinguish it from the conventional one, the GR-$\Lambda$CDM model,
based on General Relativity. Our work is a highly expanded and fully updated
analysis of our previous and much shorter presentation [83], in which we have
replaced the Planck 2015 data by the Planck 2018 data, and we have included
additional sets of modern cosmological observations. We reconfirm the results
of [83] and provide now a bunch of new results which fit in with the
conclusions of our previous work and reinforce its theses. To wit: in the
light of the figures and tables that we have presented in the current study we
may assert that the BD-$\Lambda$CDM model fares better not only as compared to
the GR-$\Lambda$CDM with a rigid cosmological constant (cf. Tables 3-7 and 10)
but also when the CC term is replaced with a dynamical parametrization of the
DE, such as the traditional XCDM, which acts as a benchmark (cf. Table 8). We
find that the GR-XCDM is completely unable to enhance the value of $H_{0}$
beyond that of the concordance model. In particular, in Tables 3-4 we can see
that the information criteria (Bayes factor and Deviance Information
Criterion) do favor significantly and consistently the BD-$\Lambda$CDM model
as compared to GR-$\Lambda$CDM and GR-XCDM. There is a very good resonance
between the Bayesian evidence criterion and the DIC differences, which
definitely uphold the BD framework at a level of $+5$ units for the
Baseline+$H_{0}$ dataset scenario, meaning that the degree of support of BD
versus GR is in between positive to strong (cf. Sec. 7.7). This support is
further enhanced up to more than $+9$ units, hence in between strong to very
strong, for the case when we include the Strong Lensing data in the fit (see
Table 4). The exact Bayesian evidence curves computed in Fig. 9 reconfirm
these results in a graphical way. The pure baseline dataset scenario (in which
the local $H_{0}$ value is not included) shows weak evidence; however, as soon
as the local $H_{0}$ value is fitted along with the remaining parameters the
evidence increases rapidly and steadily, reaching the status of positive,
strong and almost very strong depending on the datasets.
Another dataset scenario which is particularly favored in our analysis is the
one based on considering the effective calibration prior on the absolute
magnitude $M$ of the nearer SNIa data in the distance ladder (as defined in
Sec. 7.6), instead of the local value $H_{0}$ from SH0ES. The results for the
Baseline dataset in combination with $M$ (denoted B+$M$) can be read off from
Table 6 (first row). We can see it yields a tantalizing overall output, with
values of $H_{0}$ and $\sigma_{8}$ in the correct ranges for solving the two
tensions, and fully compatible with the results obtained using the prior on
$H_{0}$, as expected. In addition, the corresponding Bayes factor for this
scenario points to a remarkably high value $2\ln B>+10$, thereby carrying a
very strong Bayesian evidence, in fact comparable to the Baseline+$H_{0}$+SL
scenario of Table 4. In all of the mentioned cases in our summary the
information criteria definitely endorse the BD-cosmology versus the GR one.
Finally, we have also assessed the influence of the neutrino masses in the
context of the BD-$\Lambda$CDM model, see Table 11. We have found that massive
neutrinos can help to further reduce the predicted value of $\sigma_{8}$ to
the level of what is precisely needed to describe the weak gravitational
lensing observations derived from direct shear data. This would completely
dissolve the $\sigma_{8}$-tension without detriment of the positive results
obtained to loosen the $H_{0}$-tension, i.e. by preserving the ‘golden rule’.
Such a conclusion is, however, provisional as it requires a devoted study of
the neutrino sector (extending the analysis of Sec. 7.8.8) which is beyond the
scope of the current work.
Overall, the statistical support in favor of the BD-$\Lambda$CDM model against
the concordance GR-$\Lambda$CDM model is rather significant. It is not only
that the $H_{0}$ and $\sigma_{8}$ tensions are simultaneously dwarfed to a
level where they are both rendered inessential ($\lesssim 1.5\sigma$), but
also that all tested BD scenarios involving the local $H_{0}$ value provide a
much better global fit than the concordance model on the basis of a rich and
updated set of modern observations from all the main cosmological data sources
available at present. If we take into account that the BD-$\Lambda$CDM
framework is not just some ad hoc phenomenological toy-model, or some last-
minute smart parametrization just concocted to solve or mitigate the two
tensions, but the next-to-leading fundamental theory candidate directly
competing with GR, it may give us a sense of the potential significance of
these results.

## Conclusions

Throughout this thesis a wide variety of models, beyond the standard one, have
been studied. The corresponding analyses have been carried out by studying in
detail the theoretical predictions at the background and perturbation level,
with the purpose of testing them with the large amount of cosmological data to
which we currently have access. The ultimate goal is to see if we can detect
signs of new physics that help to alleviate some of the tensions that affect
the $\Lambda$CDM. As we have stated before, the concordance model, has
remained robust and unbeaten for a long time since it is roughly consistent
with a large body of cosmological data. Because of this fact, it is not
reasonable to look for models with a very different behaviour than the
$\Lambda$CDM, but to study models that exhibit small departures with respect
to the standard model in key aspects.
We have studied the Running Vacuum Models (RVM’s) in depth. They are
characterized by having a time-evolving vacuum energy density, whose
functional expression is motivated in the context of quantum field theory
(QFT) in curved space-time. Neglecting the terms that are not relevant for the
late-time Universe, the vacuum energy density, can be written as:

$\rho_{\rm
vac}(H,\dot{H})=\frac{3}{8\pi{G_{N}}}(c_{0}+\nu{H^{2}}+\alpha{\dot{H}}).$

The constant term $c_{0}$ turns out to be a key ingredient of the above
expression. Without it, it is not possible to generate the expected transition
from a decelerated to an accelerated Universe and, what is more, the fit of
the structure formation data is ruined. This is in total agreement with what
we have said above, since the small dimensionless parameters $\nu$ and
$\alpha$ introduce a mild evolution of the vacuum energy density with respect
to the standard behaviour, represented by the constant term $c_{0}$.

We have also studied the Peebles & Ratra model, which is a particularly
successful scalar field model $\phi$CDM for which the potential takes the form
$V(\phi)\sim\phi^{-\alpha}$. The dimensionless parameter $\alpha$ encodes the
extra degree of freedom that this model has with respect to the standard
model. It is found to be small and positive, therefore $V(\phi)$ can mimic and
approximate cosmological constant that is decreasing slowly with time. In the
late Universe the contribution of the scalar field, $\phi$, surpasses the
matter density, thus becoming the dominant component.

Not all the models studied are motivated within a theoretical framework. For
instance we have analyzed the performance of a couple of dark energy
parameterizations, namely, the XCDM and the CPL. They were introduced as the
simplest way to track a possible dynamics for dark energy. The rigid
cosmological constant is replaced by an unspecified quantity $X$ whose energy
density coincides with $\rho^{0}_{\Lambda}$ at present. For the XCDM
parameterization we have a constant equation of state (EoS) parameter, but
different from $-1$, whereas CPL is characterized by having a time-evolving
EoS.
We have also contrasted the performance of the RVM against two
phenomenological models, denoted by: $Q_{dm}$ and $Q_{\Lambda}$, respectively.
For both a totally ad-hoc expression for $\dot{\rho}_{\rm vac}$ is chosen, and
then in order to allow the time evolution of the vacuum energy density an
interaction with dark matter (DM) is considered.

Last but not least, at the end of the thesis the Brans & Dicke (BD) gravity
model was studied in detail. The main feature of this model is that the
Newtonian constant coupling $G_{N}$ is replaced by a dynamical scalar field
$G(t)=1/\psi(t)$, coupled to the curvature. As a consequence the gravitational
interaction is not only mediated by the metric field, as in the General
Relativity (GR) case but also for the aforementioned scalar field $\psi$. The
theory contains a dimensionless parameter, named BD-parameter and denoted by
$\omega_{\rm BD}$. Unlike in the original formulation of the model we do
consider the presence of the constant $\rho_{\Lambda}$ in the action. Defined
in this way, the model contains two extra d.o.f. with respect to the standard
model, to wit: the before mentioned parameter $\omega_{\rm BD}$ and the
initial value of the BD scalar field. Therefore, in order to recover the
$\Lambda$CDM we have to enforce both limits, $\omega_{\rm
BD}\rightarrow\infty$ and $\psi\rightarrow{G_{N}}$. As we have seen the fact
of having these two degrees of freedom turns out to be fundamental to solve
the $H_{0}$-tension and $\sigma_{8}$-tension at a time.

Let us summarize the main conclusions obtained in this dissertation:

  * •

In the light of the results obtained we can confirm that the general class of
RVM appears as a serious candidate for the description of the current state of
the Universe in accelerated expansion. Very clear signs of dynamical vacuum
energy have been found, being the evidence in the range 3-4$\sigma$ depending
on the dataset employed. The dimensionless parameter $\nu$, which can be
related with the $\beta$-function of the running, takes a value of order
$\sim\mathcal{O}(10^{-3})$, which is in complete agreement with the
theoretical predictions. A decreasing vacuum energy density (due to $\nu>0$)
(which is more natural from a thermodynamic point of view) throughout the
cosmic history is favoured by the data, meaning that the amount of vacuum
energy grows in the past. As a result the repulsive force exerted by this
component is greater in the past that it is in the $\Lambda$CDM, therefore
there is less structure formation predicted in this model. This important
feature is revealed as essential to loosen the $\sigma_{8}$-tension that
affects the standard model of cosmology.

  * •

By looking at the results obtained for the $\phi$CDM, with the Peebles & Ratra
potential, we can reconfirm that the hypothesis $\Lambda={\rm const.}$ is once
again, disfavoured. As in the case of the RVM’s the level of significance
achieved is between 3-4$\sigma$ and the dimensionless parameter $\alpha$,
which is responsible for the slow evolution of the scalar field, takes a value
of order $\sim\mathcal{O}(0.1)$. It is reassuring the fact of finding the same
level of evidence in favor of dynamical DE even when two different models like
$\phi$CDM and the general class of RVM’s, are test against the very same
dataset. Both type of models seem to point out in the same direction, namely,
that the DE is decreasing with the expansion and therefore that it behaves
effectively as quintessence.

  * •

Even for the XCDM and CPL, which are simple DE parameterizations, it is
possible to detect significant signs of evolving dark energy. This may be
pointing out that perhaps the signature of such dynamical signature is sitting
in the cosmological data and therefore is not exclusive of a particular model.
On the other hand the level of significance does indeed vary from one model to
another. In fact, this is what happens if we compare the performance of the
phenomenological models $Q_{dm}$ and $Q_{\Lambda}$ with the performance of the
RVM. Neither of the first two models reaches the level of significance
achieved by the Running Vacuum Model.

  * •

The quality fits obtained in the different chapters of this thesis, that point
out the preference of the cosmological data for the option $\Lambda\neq{\rm
const.}$ have been reassessed with the help of the time-honored Akaike and
Bayesian information criteria. These two tests, conveniently penalize the fact
of having a greater number of degrees of freedom. According to the statistical
standards, the evidence in favour of the dynamical dark energy goes from
“positive evidence” to “very strong evidence” depending on the model under
consideration. This supports the conclusion that there exists an improvement
of the description of the cosmological data, when we consider a dynamical DE
model, in comparison to the $\Lambda$CDM.

  * •

By studying the effect that each one of the datasets has on the results, we
have identified which are the most important data ingredients responsible for
the dynamical dark energy signal. When the triad BAO+LSS+CMB is not included
in the analysis the signal is completely lost. The impact of the $H(z)$ and
SNIa data turns out to be more moderate, however they are still important to
consolidate the results. We have also verified that the inclusion of the
Gravitational Lensing data, either the Weak Lensing one or the Strong Lensing
data, has a significant impact on the results.

  * •

While the dynamical vacuum energy models and also the quasi-vacuum models do a
great job of reducing the $\sigma_{8}$-tension, they are not able to alleviate
by any means the $H_{0}$-tension. When the LSS data is included in the dataset
the values obtained for the $H_{0}$ are in full agreement with the one
obtained by the Planck collaboration, but in tension with the value reported
by the SH0ES team obtained from the application of the cosmic distance ladder
method.

  * •

In regards to the BD model, or BD-$\Lambda$CDM as we have called it in Chapter
7, the results clearly show that even with a $G_{\rm eff}$ being higher than
in the GR-$\Lambda$CDM, the model is capable to fit the wealth of the
cosmological data employed in the analysis. The main consequence of having
higher value of $G_{\rm eff}$ is that the value of the $H_{0}$ parameter can
be in the range 70-72 km/s/Mpc (depending on the dataset utilized),
alleviating then the aforementioned $H_{0}$-tension. What is more, due to the
role played by the $\omega_{\rm BD}$ parameter in helping to suppress the
structure formation processes in the Universe, the model is also able to
loosen the $\sigma_{8}$-tension. Therefore, the main outcome that can be
extracted from the comprehensive analysis of the BD-$\Lambda$CDM model is that
it is able to alleviate the main two tensions affecting the GR-$\Lambda$CDM at
a time.

As it can be appreciated the results obtained in this thesis, after years of
studying the different cosmological models in depth, provide very encouraging
signs in favor of dynamical dark energy. It is very remarkable the fact that
this signal has been perceived even if the models are based on totally
different assumptions which clearly reinforces the conclusions collected. Such
promising results may be reconfirmed thanks to the incoming new generation of
cosmological data. Perhaps, in a very short time we will obtain more clues
that help us to decipher some of the existing mysteries in cosmology, either
from the observational side or from the theoretical perspective. We are
entering in an era that promises to be fascinating for cosmology !

##  Appendix A Sign conventions in General Relativity

Due to the existence of different sign conventions in cosmology we deem it is
necessary to dedicate this appendix to see how the different elements that we
work with are affected by these conventions. 303030This appendix is based on
the notes that the supervisor of this thesis, the Professor Joan Solà
Peracaula, prepared for his students.

Let us define the following four signs, $S1$ (metric), $S2$ (Riemann), $S3$
(Einstein) and $S_{4}$ (Ricci). Only three of them are independent. We
arbitrarily take to be the first three the independent ones, and we call the
sign sequence

$(S1,S2,S3)$ (A.1)

the “MTW sign convention” for a given author. The name comes from the time-
honored truly encyclopedic GR book by Misner-Thorn-Wheeler [368]. Each GR book
uses a given one of the 8 possible conventions. For this book the convention
is, of course, $(+,+,+)$. These three fundamental signs are defined as
follows:

$g_{\mu\nu}=[S1]\times\left(-,+,+,+\right)\,,$ (A.2)
$R^{\alpha}_{\mu\beta\nu}=[S2]\times\left(\partial_{\beta}\Gamma_{\mu\nu}^{\alpha}-\partial_{\nu}\Gamma^{\alpha}_{\mu\beta}+\Gamma^{\alpha}_{\sigma\beta}\Gamma^{\sigma}_{\nu\mu}-\Gamma^{\alpha}_{\sigma\nu}\Gamma^{\sigma}_{\beta\mu}\right),$
(A.3) $G_{\mu\nu}=[S3]\times 8\pi G_{N}T_{\mu\nu}\,.$ (A.4)

The fourth sign is for the Ricci tensor, but is no longer independent if the
first three ones are given:

$R_{\mu\nu}=[S4]\times R^{\alpha}_{\mu\alpha\nu}=[S2]\times[S3]\times
R^{\alpha}_{\mu\alpha\nu}=[S3]\times\left(\partial_{\alpha}\Gamma_{\mu\nu}^{\alpha}-...\right)\,.$
(A.5)

There is no need to show the rest of the terms in the Riemann tensor since it
all depends on the starting ones and then the rest is unambiguously fixed.

The line element squared for a metric widely employed in this thesis, the FLRW
one, therefore reads

$ds^{2}=[S1]\times\left(-dt^{2}+a^{2}(t)\delta_{ij}dx^{i}dx^{j}\right).$ (A.6)

Sometimes the MTW convention is expressed as

$(S1,S2,S4)\,.$ (A.7)

In that case $[S3]$ is fixed from $[S3]=[S2]\times[S4]$.

The following observations should be bear in mind for the rest of the
appendix:

  * •

The Christoffel symbols do not change under metric conventions since the
metric appears twice, derivately and non-derivatively.

  * •

$R_{\mu\nu}$ does not flip sign under metric flip changes. The reason is that
the Riemann tensor is made of Christoffel symbols and hence does not flip sign
under metric flips, and then $R_{\mu\nu}$ involves one more bilinear
dependence on the metric that leads to a $\delta$-Kronecker, and hence no
additional sign flip:

$R_{\mu\nu}=[S4]\times R^{\alpha}_{\mu\alpha\nu}=[S4]\times
g^{\alpha\beta}R_{\alpha\mu\beta\nu}=[S4]\times
g^{\alpha\beta}g_{\alpha\sigma}R^{\sigma}_{\mu\beta\nu}=[S4]\times\delta^{\beta}_{\alpha}R^{\alpha}_{\mu\beta\nu}.$
(A.8)

  * •

In contrast, the sign of the Ricci scalar

$R=g^{\mu\nu}R_{\mu\nu}\,,$ (A.9)

does flip with the metric since $R_{\mu\nu}$ does not.

  * •

It follows that $g_{\mu\nu}R$ does not flip with the change of the metric
signature and hence the entire Einstein tensor
$G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R$ does not flip either.

  * •

Let us repeat here that if we take the first three signs to be independent,
then the Ricci sign (i.e. the sign $[S4]$ on the r.h.s. of Einstein equations,
see Eq. (A.5)) is fixed as follows:

$[S4]=[S2]\times[S3]\,,$ (A.10)

as it can be seen from (A.8) and (A.4) and the definition of $G_{\mu\nu}$.
This can be shown by substituting the last piece of the r.h.s. of (A.5), which
already used Eq.(A.10), into $G_{\mu\nu}$ on the l.h.s. of (A.4) and observe
that all signs cancel out in Einstein’s equations, once these are expressed in
terms of $\partial_{\alpha}\Gamma_{\mu\nu}^{\alpha}-...$, as these terms are
no longer dependent on any sign, not even on $[S1]$. Conversely, if we use the
alternative signs (A.7), then the Einstein sign is fixed through
$[S3]=[S2]\times[S4].$

Let us now derive the signs of Einstein’s equations from the signs in the EH-
action, which in any MTW convention must have the form

$\displaystyle S_{\rm EH}$ $\displaystyle=$
$\displaystyle\frac{1}{16\pi\,G_{N}}\,\int
d^{4}x\sqrt{-g}\,\left([S1]\times[S3]\times\,R-2\,\Lambda\,\right)$ (A.11)
$\displaystyle=$ $\displaystyle\ \ \ \int
d^{4}x\sqrt{-g}\,\left(\frac{[S1]\times[S3]}{16\pi\,G_{N}}\,R-\rho_{\Lambda}\right)\,.\
\ $

Notice that the sign in front of $\Lambda$ or $\rho_{\Lambda}$ must _always_
be negative since for $\Lambda>0$ it must give a positive contribution to the
potential energy (density) of the Lagrangian. As for the sign of $R$, recall
it flips sign with the metric, hence $[S1]\times R$ does not, and finally
there must be $[S3]$ because this controls the sign of $G_{\mu\nu}$ and hence
of $R_{\mu\nu}$, see Eq. (A.5), which appears in the variation of $R$. Indeed,
let us check that. Recall the formulae:

$\delta\sqrt{-g}=-\frac{1}{2}\,\sqrt{-g}\,g_{\mu\nu}\delta g^{\mu\nu}\,,\ \ \
\ \ \ \ \delta\left(\sqrt{-g}\,R\right)=\sqrt{-g}\,G_{\mu\nu}\,\delta
g^{\mu\nu}\,.$ (A.12)

Now we must variate the total action $S=S_{\rm EH}+S_{m}$, where $S_{m}$ is
the matter action. The latter is related to the energy-momentum tensor as
follows:

$T^{\mu\nu}=+[S1]\times\frac{2}{\sqrt{-g}}\frac{\delta S_{m}}{\delta
g_{\mu\nu}}\ \ \Longleftrightarrow\ \ \delta
S_{m}=+\frac{1}{2}\,[S1]\times\int d^{4}x\sqrt{-g}\,T^{\mu\nu}\,\delta
g_{\mu\nu}\,.$ (A.13)

From these formulas it is now straightforward to compute the total variation
of the action $\delta S=\delta S_{\rm EH}+\delta S_{m}$ and equate it to zero:

$\displaystyle\delta S$ $\displaystyle=$ $\displaystyle\int
d^{4}x\,\sqrt{-g}\,\left\\{\frac{1}{16\pi\,G_{N}}\left([S1]\times[S3]\times
G_{\mu\nu}+\Lambda g_{\mu\nu}\right)-\frac{1}{2}[S1]\times
T_{\mu\nu}\right\\}\delta g^{\mu\nu}=0.$

Equating the braced expression to zero and multiplying it all with $16\pi
G_{N}$ we get

$[S1]\times[S3]\times G_{\mu\nu}+\Lambda g_{\mu\nu}=[S1]\times 8\pi
G_{N}T_{\mu\nu}.$ (A.15)

Multiplying this equation by $[S1]\times[S3]$ we find Einstein’s equations in
all MTW conventions:

$G_{\mu\nu}+[S1]\times[S3]\times\Lambda\,g_{\mu\nu}=[S3]\times 8\pi G_{N}\
{T}_{\mu\nu}.$ (A.16)

So far we have been dealing with the geometric part of Einstein’s equations,
but now it is time to consider some of the particular forms that the energy-
momentum tensor can take and how they change under the sign conventions. Let
us consider two of the most common cases the energy-momentum tensor of a
scalar field and the energy-momentum tensor of a perfect fluid.

### Signs for the energy-momentum tensor of a scalar field

Notice that since the coordinates are defined carrying upper indices, the
metric tensor to which they couple is $g_{\mu\nu}$ rather than its inverse
$g^{\mu\nu}$. So the primary gravitational field is $g_{\mu\nu}$, not
$g^{\mu\nu}$, and with respect to the primary one we have a $+$ sign in the
definition of $T^{\mu\nu}$ – see (A.13). However, owing to $\delta
g^{\mu\nu}=-g^{\mu\alpha}g^{\beta\nu}\delta g_{\alpha\beta}$ an extra sign
appears, which comes from the fact that
$g^{\mu\alpha}g_{\alpha\nu}=\delta^{\mu}_{\nu}$. Incidentally, the one-
dimensional analog is $\delta x^{-1}=-x^{-2}\delta x$, which appears from
differentiating $x^{-1}\,x=1$. This allows to extend the formulae (A.12) as
follows:

$\delta\sqrt{-g}=-\frac{1}{2}\,\sqrt{-g}\,g_{\mu\nu}\delta
g^{\mu\nu}=+\frac{1}{2}\,\sqrt{-g}\,g^{\mu\nu}\delta g_{\mu\nu}\,,\ \ \ \
\delta\left(\sqrt{-g}\,R\right)=\sqrt{-g}\,G_{\mu\nu}\,\delta
g^{\mu\nu}=-\sqrt{-g}\,G^{\mu\nu}\,\delta g_{\mu\nu}\,.$ (A.17)

The sign difference between using $\delta g^{\mu\nu}$ and $\delta g_{\mu\nu}$
is transferred to the corresponding formula for $T_{\mu\nu}$ with lower
indices:

$T_{\mu\nu}=-[S1]\times\frac{2}{\sqrt{-g}}\frac{\delta S_{m}}{\delta
g^{\mu\nu}}\ \ \Longleftrightarrow\ \ \delta
S_{m}=-\frac{1}{2}\,[S1]\times\int d^{4}x\sqrt{-g}\,T_{\mu\nu}\,\delta
g^{\mu\nu}\,.$ (A.18)

We are going to work out the expression for the energy-momentum tensor of a
scalar field from the corresponding action

$\displaystyle S_{\Phi}=\int
d^{4}x\,\sqrt{-g}\,\left[-[S1]\times\frac{1}{2}\,g^{\mu\nu}\,\partial_{\mu}\Phi\,\partial_{\nu}\Phi-U(\Phi)\right]\,,$
(A.19)

we obtain the generic form for $T^{\mu\nu}$

$T^{\mu\nu}=\partial^{\mu}\Phi\,\partial^{\nu}\Phi+[S1]\times
g^{\mu\nu}\left[-[S1]\times\frac{1}{2}g^{\alpha\beta}\partial_{\alpha}\Phi\,\partial_{\beta}\Phi-U(\Phi)\right]=\partial^{\mu}\Phi\,\partial^{\nu}\Phi+[S1]\times
g^{\mu\nu}{\cal L}$ (A.20)

or

$T_{\mu\nu}=\partial_{\mu}\Phi\,\partial_{\nu}\Phi+[S1]\times
g_{\mu\nu}\left[-[S1]\times\frac{1}{2}g^{\alpha\beta}\partial_{\alpha}\Phi\,\partial_{\beta}\Phi-U(\Phi)\right]=\partial_{\mu}\Phi\,\partial_{\nu}\Phi+[S1]\times
g_{\mu\nu}{\cal L}\,.$ (A.21)

The previous calculations follow from the standard rule for functional
differentiation, but it is possible to just use normal partial differentiation
for any action of the form

$S_{\Phi}=\int d^{4}x\,\sqrt{-g}\,{\cal L}(\Phi,g_{\mu\nu},\partial\Phi)$
(A.22)

in which the Lagrangian (density) ${\cal L}$ does not depend on the
derivatives of the metric, as in (A.19). Using (A.13) we have (see below):

${T^{\mu\nu}=[S1]\times\left[-2\frac{\partial\cal L}{\partial
g_{\mu\nu}}+g^{\mu\nu}{\cal L}\right]}.$ (A.23)

The same formula is valid with lower indices and no sign change anywhere, as
can be easily checked. Using these formulas we can easily recover equations
(A.20) and (A.21) using normal partial differentiation.

To check (A.23) we can either use variational calculus or functional
differentiation. Let us start with the latter. Let us recall that from the
point of view of functional differentiation the formulas (A.17) must be
interpreted conveniently. For instance,

$\frac{\delta\sqrt{-g(x)}}{\delta
g^{\mu\nu}(x^{\prime})}=-\frac{1}{2}\,\sqrt{-g(x)}\,g_{\mu\nu}\,\delta^{(4)}(x-x^{\prime})\,,\
\ \ \ \frac{\delta\sqrt{-g}}{\delta
g_{\mu\nu}(x^{\prime})}=\frac{1}{2}\,\sqrt{-g}\,g^{\mu\nu}\,\delta^{(4)}(x-x^{\prime})\,.$
(A.24)

Similarly, from $\delta g^{\mu\nu}=-g^{\mu\alpha}g^{\beta\nu}\delta
g_{\alpha\beta}$ we obtain

$\frac{\delta g^{\mu\nu}(x)}{\delta
g_{\alpha\beta}(x^{\prime})}=-g^{\mu\alpha}g^{\beta\nu}\delta^{(4)}(x-x^{\prime})\,.$
(A.25)

Using these formulas we may functionally differentiate (A.22):

$\frac{\delta S_{\Phi}}{\delta g_{\mu\nu}(x^{\prime})}=\int
d^{4}x\left(\frac{\delta\sqrt{-g(x)}}{\delta g_{\mu\nu}(x^{\prime})}\,{\cal
L}+\frac{\delta{\cal L}}{\delta g_{\mu\nu}(x^{\prime})}\sqrt{-g(x)}\right)$
(A.26)

where the first term in the parenthesis is given by the second in (A.24). And
the second in the parenthesis is computed with the help of (A.25), rendering:

$\frac{\delta{\cal L}}{\delta g_{\mu\nu}(x^{\prime})}=\frac{\delta}{\delta
g_{\mu\nu}(x^{\prime})}\left[-[S1]\times\frac{1}{2}\,g^{\alpha\beta}\,\partial_{\alpha}\Phi\,\partial_{\beta}\Phi-U(\Phi)\right]$
$=-[S1]\times\frac{1}{2}\,\frac{\delta g^{\alpha\beta}(x)}{\delta
g_{\mu\nu}(x^{\prime})}\,\partial_{\alpha}\Phi\,\partial_{\beta}\Phi=+[S1]\times\frac{1}{2}\,\partial^{\mu}\Phi\,\partial^{\nu}\Phi\,\delta^{(4)}(x-x^{\prime}).$
(A.27)

Integrating out the $\delta$-functions, and trivially noting that
$[S1]^{2}=1$, we find from (A.13) and the previous results:

$T^{\mu\nu}=+[S1]\times\frac{2}{\sqrt{-g}}\frac{\delta S_{\Phi}}{\delta
g_{\mu\nu}}=\partial^{\mu}\Phi\,\partial^{\nu}\Phi+[S1]\times g^{\mu\nu}{\cal
L}\,.$ (A.28)

This result is exactly equivalent to (A.23) because the second term is already
matched, whereas the first term is just also right so:

$[S1]\times\left[-2\frac{\partial\cal L}{\partial
g_{\mu\nu}}\right]=[S1]\times\left[-2[S1]\times\frac{1}{2}\frac{\partial
g^{\alpha\beta}}{\partial
g_{\mu\nu}}\,\partial_{\alpha}\Phi\partial_{\beta}\Phi\right]=\partial^{\mu}\Phi\,\partial^{\nu}\Phi\,,$
(A.29)

where from $g^{\alpha\beta}=g^{\alpha\mu}g^{\beta\nu}g_{\mu\nu}$ we get
$\partial g^{\alpha\beta}/\partial g_{\mu\nu}=g^{\alpha\mu}g^{\beta\nu}$.
Notice that upper and lower indices variables are treated as independent
variables in the differentiation process. We could equally proceed
differentiating with respect to $g^{\mu\nu}$ (i.e. with upper indices) and
compute

$\frac{\delta S_{\Phi}}{\delta g^{\mu\nu}(x^{\prime})}=\int
d^{4}x\left(\frac{\delta\sqrt{-g(x)}}{\delta g^{\mu\nu}(x^{\prime})}\,{\cal
L}+\frac{\delta{\cal L}}{\delta g^{\mu\nu}(x^{\prime})}\sqrt{-g(x)}\right)$
(A.30)

where the first term in the parenthesis is given by the first of (A.24) and
has now a minus sign, and the second is computed from

$\frac{\delta{\cal L}}{\delta g^{\mu\nu}(x^{\prime})}=\frac{\delta}{\delta
g^{\mu\nu}(x^{\prime})}\left[-[S1]\times\frac{1}{2}\,g^{\alpha\beta}\,\partial_{\alpha}\Phi\,\partial_{\beta}\Phi-U(\Phi)\right]$
$=-[S1]\times\frac{1}{2}\,\frac{\delta g^{\alpha\beta}(x)}{\delta
g^{\mu\nu}(x^{\prime})}\,\partial_{\alpha}\Phi\,\partial_{\beta}\Phi=-[S1]\times\frac{1}{2}\,\partial_{\mu}\Phi\,\partial_{\nu}\Phi\,\delta^{(4)}(x-x^{\prime})$
(A.31)

where we have used

$\frac{\delta g^{\mu\nu}(x)}{\delta
g^{\alpha\beta}(x^{\prime})}=+\delta^{\mu}_{\alpha}\delta^{\beta}_{\nu}\delta^{(4)}(x-x^{\prime})\,.$
(A.32)

Note the following: since all metric indices are in upper positions, the
previous relation obviously does not introduce an extra sign – in contrast to
(A.25). These sign changes in each term are now compensated with the extra
minus sign that the definition of $T_{\mu\nu}$ (with lower indices) involves,
and we obtain the desired result:

$T_{\mu\nu}=-[S1]\times\frac{2}{\sqrt{-g}}\frac{\delta S_{\Phi}}{\delta
g^{\mu\nu}}=\partial_{\mu}\Phi\partial_{\nu}\Phi+[S1]\times g_{\mu\nu}{\cal
L}\,.$ (A.33)

This result is consistent with (A.28) because they must be related of course
through

$T^{\mu\nu}=g^{\mu\alpha}\,g^{\nu\beta}\,T_{\alpha\beta}\,.$ (A.34)

One can entirely avoid using functional differentiation and just use
variational calculus. In this case one has to vary the action and then
identify $T^{\mu\nu}$ using the equivalent expression on the r.h.s. of Eqs.
(A.13) or (A.18) etc. The results are of course the same. For example, let us
variate (A.22) with respect to the metric with _upper_ indices, $\delta
g^{\mu\nu}$:

$\delta S_{\Phi}=\int d^{4}x\left[\delta\sqrt{-g}\,{\cal
L}+\sqrt{-g}\left(-[S1]\times\frac{1}{2}\delta
g^{\mu\nu}\partial_{\mu}\Phi\partial_{\nu}\Phi\right)\right]$ $=\int
d^{4}x\sqrt{-g}\left[-\frac{1}{2}\,g_{\mu\nu}{\cal
L}-[S1]\times\frac{1}{2}\partial_{\mu}\Phi\partial_{\nu}\Phi\right]\delta
g^{\mu\nu}\,.$ (A.35)

Comparing with the second expression of (A.18) we must have

$-\frac{1}{2}[S1]\times T_{\mu\nu}=-\frac{1}{2}\,g_{\mu\nu}{\cal
L}-[S1]\times\frac{1}{2}\partial_{\mu}\Phi\partial_{\nu}\Phi.$ (A.36)

Multiplying both sides for $-2[S1]$ we find the desired result:

$T_{\mu\nu}=\partial_{\mu}\Phi\,\partial_{\nu}\Phi+[S1]\times g_{\mu\nu}{\cal
L}\,.$ (A.37)

Had we varied the action (A.22) with respect to the metric with _lower_
indices, $\delta g_{\mu\nu}$, we would have started again from

$\delta S_{\Phi}=\int d^{4}x\delta\sqrt{-g}\,{\cal
L}+\sqrt{-g}\left[-[S1]\times\frac{1}{2}\delta
g^{\mu\nu}\partial_{\mu}\Phi\partial_{\nu}\Phi\right]\,,$ (A.38)

but now we use $\delta\sqrt{-g}=+\frac{1}{2}\,\sqrt{-g}\,g^{\mu\nu}\delta
g_{\mu\nu}$ for the first term, whereas we use

$\delta
g^{\mu\nu}\partial_{\mu}\Phi\,\partial_{\nu}\Phi=-\,g^{\mu\alpha}\,g^{\nu\beta}\delta
g_{\alpha\beta}\,\partial_{\mu}\Phi\,\partial_{\nu}\Phi=-\partial^{\alpha}\Phi\,\partial^{\beta}\Phi\,\delta
g_{\alpha\beta}=-\partial^{\mu}\Phi\,\partial^{\nu}\Phi\,\delta g_{\mu\nu}\,,$
(A.39)

for the second term. Therefore (A.38) becomes

$\delta S_{\Phi}=\int d^{4}x\sqrt{-g}\left[+\frac{1}{2}\,g^{\mu\nu}{\cal
L}+[S1]\times\frac{1}{2}\partial^{\mu}\Phi\partial^{\nu}\Phi\right]\delta
g_{\mu\nu}\,.$ (A.40)

Comparing with the second expression of (A.13) we must have

$+\frac{1}{2}[S1]\times T_{\mu\nu}=+\frac{1}{2}\,g_{\mu\nu}{\cal
L}+[S1]\times\frac{1}{2}\partial^{\mu}\Phi\partial^{\nu}\Phi.$ (A.41)

Multiplying both sides for $2[S1]$ we find

$T^{\mu\nu}=\partial^{\mu}\Phi\,\partial^{\nu}\Phi+[S1]\times g^{\mu\nu}{\cal
L}\,,$ (A.42)

which is consistent once more with the previous form (A.37) because both must
be related through Eq. (A.34). It is perhaps faster with variational calculus,
but functional differentiation is a very useful tool in many other situations
and it is very convenient to practice it as well. It is actually the most
rigorous procedure, see below.

Let us notice that while the functional differentiation was unambiguous all
the time and provided the right answer whether differentiating with respect to
$g_{\mu\nu}$ or $g^{\mu\nu}$, the variational calculation was ambiguous at the
point (A.38). In fact, we could have written the kinetic Lagrangian as ${\cal
L}=-[S1]\times(1/2)g_{\mu\nu}\partial^{\mu}\Phi\,\partial_{\nu}\Phi$ where now
the metric indices are down and the partial derivative indices up. The
variation with respect to the metric could have been interpreted simply as as
$\delta{\cal L}=-[S1]\times(1/2)\delta
g_{\mu\nu}\partial^{\mu}\Phi\,\partial_{\nu}\Phi$, with no sign change as
compared to (A.39). If so, the sign of the first term of (A.42) would be
negative , and hence wrong because it would very evidently violate the
fundamental tensor relation (A.42).

### Signs for the energy-momentum tensor of a perfect fluid

Let us also note that the energy-momentum tensor for a perfect fluid is given
by

$T_{\mu\nu}=[S1]\times pg_{\mu\nu}+(\rho+p)U_{\mu}U_{\nu}$ (A.43)

where we generically consider a fluid characterized by its energy density
$\rho$ and the corresponding pressure $p$. The square of the 4-velocity vector
satisfies the closure relation

$U_{\mu}U^{\mu}=g_{\mu\nu}U^{\mu}U^{\nu}=-[S1].$ (A.44)

The above relation is fulfilled in any frame since it is an invariant. Let us
now consider the components of the energy-momentum tensor in the proper frame,
where we simply have $U^{\mu}=(-1,0,0,0)$. It is important to note that the
components have different form depending on whether we use the mixed tensor
$T^{\mu}_{\nu}$, with upper and lower indices, or the original tensor
$T_{\mu\nu}$ given in (A.43) with all indices in the same position. As we
shall see, the first is easy to operate because the $\delta$-Kronecker appears
in it;

$T^{\mu}_{\nu}=[S1]\times p\delta^{\mu}_{\nu}+(\rho+p)U^{\mu}U_{\nu}\,,$
(A.45)

but on the other hand the components of this tensor do depend on the signature
of the metric, i.e. on $[S1]$, whereas $T_{\mu\nu}$ does not. This is expected
since $T_{00}=T^{00}$ gives the sign of the energy density in _all_
conventions and therefore it cannot depend on the convention, whereas
$T_{0}^{0}$ is convention-dependent.

Let us check it. From the definition of $U^{\mu}$ given above and the property
$U_{\mu}U^{\mu}=-[S1]$ we have

$T^{0}_{0}=[S1]\times{p}-(\rho+p)\times[S1],\quad T^{i}_{i}=[S1]\times p$
(A.46)

consequently

$T^{\mu}_{\nu}=[S1]\times{\rm diag}(-\rho,p,p,p)\,.$ (A.47)

From here we can compute the trace:

$T^{\mu}_{\mu}=[S1]\times(3p-\rho)\,,$ (A.48)

which does indeed depend on the metric signature since it is not an
observable. Notice, however, from (A.2) and (A.47) that

$T^{00}=g^{\mu
0}T^{0}_{\mu}=g^{00}T^{0}_{0}=[S1]\times\left(-T^{0}_{0}\right)=[S1]\times[S1]\times\rho=\rho\,,$

and similarly

$T^{ii}=g^{\mu i}T^{i}_{\mu}=g^{ii}T^{i}_{i}=[S1]\times
T^{i}_{i}=[S1]\times[S1]\times p=p\,.$

Therefore the energy-momentum tensor with upper indices in the proper frame
reads

$T^{\mu\nu}={\rm diag}(\rho,p,p,p)\,,$ (A.49)

which is fully independent of the metric convention. These are the observable
quantities. Similarly

$T_{\mu\nu}={\rm diag}(\rho,p,p,p)\,.$ (A.50)

owing to the tensor property Eq. (A.34), which involves twice the metric
tensor.

123

123

##  Appendix B Cosmological perturbations

It is a fact that our Universe is really far from being completely homogeneous
and isotropic. Only on large scales ($\sim$100 Mpc) the description of the
expanding Universe in terms of quantities that only depend on time, remains
valid. As we have explained in the Introduction the aforementioned scenario is
well described in terms of the Friedmann-Lemaître-Robertson-Walker (FLRW)
metric and the corresponding Friedmann’s and pressure equations. However if we
want to go a step further, the previous description is no longer valid and we
need to employ the perturbation theory in order to explain the observed
structures in the Universe. So, how were these structures formed ? Let us
answer this question with a brief summary, without going into the mathematical
details, of this fascinating process based on a simple concept, the
gravitational instability. It is commonly believed that after the Big Bang,
the Universe undergone a phase of rapidly expansion called inflation, where
the space expanded itself at an approximately exponential rate. While it is
true that we do not have a solid evidence of this theory, the absence of
alternatives and the fact that it provides beautiful solutions for some
worrisome problems, like the observed flatness, the homogeneity and isotropy
of the current Universe plus the absence of exotic relics, makes it very
appealing for cosmologists. One of the most important predictions of inflation
is the generation of quantum-mechanical fluctuations when the relevant scales
were causally connected. These tiny fluctuations were amplified by inflation
and seeds the structure formation that we can observe. The end of inflation is
called reheating. This process is poorly understood, nevertheless the most
accepted theory stands that along this period the particle responsible for the
inflationary period (whatever it is) decays into a hot thermal plasma of other
particles. Therefore, at this stage the Universe was made up of a dense plasma
of photons, electrons, baryons and dark matter (DM). Because of the radiation
pressure, structure formation in the baryonic sector was inhibited at that
time, nevertheless since DM does not interact except gravitationally remained
in the center of the overdensities. The electromagnetic interaction forced
baryons and photons to move together outward from those overdense regions,
creating spherical sound waves which are the origin of the baryonic acoustic
oscillations (BAO). As the Universe expanded it cooled, so at some point the
temperature was cool enough to allow protons to capture electrons and form
neutral hydrogen atoms. This is known as the recombination time and as can be
imagined it was of tremendous importance in the cosmic history. The
interaction between photons and baryons ceased and the first ones started to
stream freely making up, what we know as, the cosmic microwave background
radiation. Baryons fell into the potential wells, previously created by DM,
forming the first structures in the Universe. The Universe continued to expand
and cool and eventually (well within the matter dominated epoch) the
gravitational attraction exerted by baryons and DM allows the structure
formation started to grow. This is the starting point where we will focus on
this appendix.

### Cosmological perturbation equations for dark energy models

In this section we present the basic steps that must be followed to obtain the
set of perturbation equations at first order. We will keep the FLRW metric as
the background metric (it will be denoted as $\bar{g}_{\mu\nu}$) and for the
rest of the chapter we stick to the convention (+,+,+) explained in the
previous appendix. We introduce a tiny perturbation in the metric in the
following way:

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=(\bar{g}_{\mu\nu}+\delta{g_{\mu\nu}})dx^{\mu}dx^{\nu}=a^{2}(\eta_{\mu\nu}+{h_{\mu\nu}})dx^{\mu}dx^{\nu},$
(B.1)

where $\eta_{\mu\nu}={\rm diag}(-1,+1,+1,+1)$ is the Minkowski metric. We have
to bear in mind that the elements of the perturbed metric $\delta{g_{\mu\nu}}$
(or equivalently ${a^{2}}h_{\mu\nu}$) are small in comparison with the
background components $\bar{g}_{\mu\nu}$. The metric field $g_{\mu\nu}$ must
be symmetric, therefore $\delta{g_{\mu\nu}}$ also has to be. This means that
the total number of d.o.f. is reduced from 16 to 10. The General Relativity
(GR) field equations are invariant under a general coordinate change, as a
consequence, we can fix 4 of the remaining d.o.f. by fixing the gauge. It can
be shown that in a completely general way the perturbed metric can be
decomposed as

$ds^{2}=a^{2}\left[-(1+2\Phi)d\eta^{2}+\omega_{i}{d\eta{dx^{i}}}+(\delta_{ij}+h_{ij})dx^{i}{dx^{j}}\right],$
(B.2)

where $\Phi$ is a scalar function, $\omega_{i}$ is a 3-vector and $h_{ij}$ is
a traceless tensor. Remember that $\eta$ represents the conformal time and it
can be related with the cosmic time through the expression
$\eta=\int^{t}_{0}dt^{\prime}/a(t^{\prime})$. There are three types of
perturbation, to wit: scalar, vector and tensor. Here we are only interested
in the first type. The vector perturbations are suppressed in the early
Universe whereas the tensor ones contribute as a gravitational waves to the
B-mode of the CMB polarizations. The scalar perturbations are the only ones
coupled to the matter perturbations and this is why here we only focus on this
particular type. The perturbed metric can be rewritten as:

$ds^{2}=a^{2}\left[-(1+2\Phi)d\eta^{2}+\partial_{i}E{d\eta{dx^{i}}}+((1+2\Psi)\delta_{ij}+D_{ij}B)dx^{i}{dx^{j}}\right]$
(B.3)

and now $\Phi$, $E$, $\Psi$ and $B$ are all scalar functions. We can now fix
the gauge by imposing some conditions over these functions. For instance, if
we want to work in the synchronous gauge we must set up $\Phi$ = $E$ = 0. On
the other hand if we want to work in the conformal Newtonian gauge then we
impose $E$ = $B$ = 0. Let us focus on the last one, for which the final form
of the perturbed metric can be written as:

$ds^{2}=a^{2}\left[-(1+2\Phi)d\eta^{2}+(1+2\Psi)\delta_{ij}dx^{i}{dx^{j}}\right].$
(B.4)

Once the perturbed metric, which is a key element in the perturbation theory,
is clearly determined we can start the process that will lead us to obtain the
perturbed cosmological equations.

The following comment about the notation employed in this appendix is in
order. When we write a quantity (could be a scalar, a vector or a tensor)
generally denoted as $X$ and we do not add nor a bar neither a $\delta$ right
in front of the term, we are referring to a quantity which is the sum of a
background part $\bar{X}$ plus a perturbed part $\delta{X}$. For instance in
the case of the metric we have (as it can be observed from (B.1))
$g_{\mu\nu}=\bar{g}_{\mu\nu}+\delta{g_{\mu\nu}}$.

To derive the first-order Einstein field equations we need to introduce a
perturbation not only in the Einstein tensor $G_{\mu\nu}$ but also in the
energy-momentum tensor $T_{\mu\nu}$. We end up with the set of equations

$\bar{G}_{\mu\nu}+\delta{G}_{\mu\nu}=8\pi{G_{N}}\left(\bar{T}_{\mu\nu}+\delta{T_{\mu\nu}}\right).$
(B.5)

From now on, we are going to drop the pure background terms since in this part
we are only interested in the perturbation equations. It turns out really
convenient to write down the covariant elements of the perturbed metric and
the corresponding contravariant elements:

$\begin{split}g_{00}=-a^{2}(1+2\Phi)\qquad
g_{ij}=a^{2}(1+2\Psi)\delta_{ij}\quad\quad
g^{00}=-\frac{1}{a^{2}}(1-2\Phi)\qquad
g^{ij}=\frac{1}{a^{2}}(1-2\Psi)\delta^{ij}.\end{split}$ (B.6)

Once we have the entries of the perturbed metric we can move forward and
compute the perturbed Christoffel symbols $\delta\Gamma^{\mu}_{\nu\lambda}$ by
using the expression

$\delta\Gamma^{\mu}_{\nu\lambda}=\frac{1}{2}\delta{g^{\mu\alpha}}(\partial_{\lambda}\bar{g}_{\alpha\nu}+\partial_{\nu}\bar{g}_{\alpha\lambda}-\partial_{\alpha}\bar{g}_{\nu\lambda})+\frac{1}{2}\bar{g}^{\mu\alpha}(\partial_{\lambda}\delta
g_{\alpha\nu}+\partial_{\nu}\delta g_{\alpha\lambda}-\partial_{\alpha}\delta
g_{\nu\lambda}).$ (B.7)

The Christoffel symbols different from zero, in the Newtonian gauge, are:

$\Gamma^{0}_{00}=\mathcal{H}+\Phi^{\prime}\qquad\Gamma^{0}_{0i}=\Gamma^{i}_{00}=\partial_{i}\Phi\qquad\Gamma^{0}_{ij}=\delta_{ij}[\mathcal{H}(1+2\Psi-2\Phi)+\Psi^{\prime}]$
(B.8)
$\Gamma^{i}_{j0}=\delta^{i}_{j}(\mathcal{H}+\Psi^{\prime})\qquad\Gamma^{i}_{jl}=\delta^{i}_{j}\partial_{l}\Psi+\delta^{i}_{l}\partial_{j}\Psi-\delta_{jl}\partial_{i}\Psi.$

Where $\mathcal{H}=\frac{1}{a}\frac{da}{d\eta}$ represents the Hubble function
in conformal time and the primes denote derivatives with respect to that
variable. We can repeat exactly the same procedure to obtain from

$\delta{R_{\mu\nu}}=\partial_{\alpha}\delta{\Gamma^{\alpha}_{\mu\nu}}-\partial_{\nu}\delta{\Gamma^{\alpha}_{\mu\alpha}}+\delta\Gamma^{\alpha}_{\mu\nu}\bar{\Gamma}^{\beta}_{\alpha\beta}+\bar{\Gamma}^{\alpha}_{\mu\nu}\delta\Gamma^{\beta}_{\alpha\beta}-\delta\Gamma^{\alpha}_{\mu\beta}\bar{\Gamma}^{\beta}_{\alpha\nu}-\bar{\Gamma}^{\alpha}_{\mu\beta}\delta\Gamma^{\beta}_{\alpha\nu}$
(B.9)

the non-null components of the Ricci tensor:

$R_{00}=-3\mathcal{H}^{\prime}+\nabla^{2}\Phi-3\Psi^{\prime\prime}+3\mathcal{H}(\Phi^{\prime}-\Psi^{\prime}),\qquad
R_{0i}=-2\partial_{i}\Psi^{\prime}+2\mathcal{H}\partial_{i}\Phi,$ (B.10)
$R_{ij}=-\partial_{i}\partial_{j}(\Psi+\Phi)+\delta_{ij}\left[(2\mathcal{H}^{2}+\mathcal{H}^{\prime})(1+2\Psi-2\Phi)-\nabla^{2}\Psi+\Psi^{\prime\prime}+5\mathcal{H}\Psi^{\prime}-\mathcal{H}\Phi^{\prime}\right].$

Contracting the indices of the previous tensor we are able to compute the
Ricci scalar

${a^{2}}R=6(\mathcal{H}^{2}+\mathcal{H}^{\prime})(1-2\Phi)-2\nabla^{2}(\Phi+2\Psi)+6\Psi^{\prime\prime}-6\mathcal{H}\Phi^{\prime}+18\mathcal{H}\Psi^{\prime}.$
(B.11)

Finally we are in a position to calculate the components of Einstein’s tensor

$G_{00}=3\mathcal{H}^{2}+6\mathcal{H}\Psi^{\prime}-2\nabla^{2}\Psi,$ (B.12)
$G_{ij}=-\partial_{i}\partial_{j}(\Psi+\Phi)+\delta_{ij}\left[-(\mathcal{H}^{2}+2\mathcal{H}^{\prime})(1+2\Psi-2\Phi)+2\mathcal{H}(\Phi^{\prime}-2\Psi^{\prime})+\nabla^{2}(\Psi+\Phi)-2\Psi^{\prime\prime}\right],$
$G_{0i}=-2\partial_{i}\Psi^{\prime}+2\mathcal{H}\partial_{i}\Phi.$

It is time now to obtain the perturbed expression for the energy-momentum
tensor of a perfect fluid, which contains the contribution of the different
species considered. The starting point is the expression valid at the
background level

$\bar{T}_{\mu\nu}=(\bar{\rho}+\bar{p})\bar{U}_{\mu}\bar{U}_{\nu}+\bar{p}{\bar{g}_{\mu\nu}},$
(B.13)

by plugging the corresponding perturbations we can get

$\delta{T}_{\mu\nu}=(\bar{\rho}+\bar{p})(\delta{U}_{\mu}\bar{U}_{\nu}+\bar{U}_{\mu}\delta{U}_{\nu})+(\delta\rho+\delta{p})\bar{U}_{\mu}\bar{U}_{\nu}+\bar{p}\delta{g}_{\mu\nu}+\delta{p}\bar{g}_{\mu\nu}.$
(B.14)

Where we denote the total energy density and pressure, at the background
level, by $\bar{\rho}$ and $\bar{p}$ respectively. We need to obtain the
components of the perturbed 4-velocity. To do so, we only require the
definition of the total 4-velocity $U^{\mu}=dx^{\mu}/d\tau$, where $\tau$ is
the proper time and the closure relation $U^{\mu}U_{\mu}=-1$

$g_{\mu\nu}\frac{dx^{\mu}}{d\tau}\frac{dx^{\nu}}{d\tau}=\left(\frac{d\eta}{d\tau}\right)^{2}\left[g_{00}+g_{ij}\frac{dx^{i}}{d\eta}\frac{dx^{j}}{d\eta}\right]=-1,$
(B.15)

where $v^{i}\equiv dx^{i}/d\eta$ is defined as the peculiar velocity with
respect to the comoving frame. If we assume that the peculiar velocities are
small quantities, we can get the perturbed expression for the first component
of the 4-velocity vector

$\frac{d\eta}{d\tau}=\frac{1}{a}(1-\Phi)+\mathcal{O}(2).$ (B.16)

Regarding the spatial components we only need to employ the chain rule
together with the aforementioned expression for the peculiar velocities

$\frac{dx^{i}}{d\tau}=\frac{dx^{i}}{d\eta}\frac{d\eta}{d\tau}=\frac{v^{i}}{a}+\mathcal{O}(2).$
(B.17)

Once we have computed the different components we can write the expression for
the 4-velocity, at first order, in the covariant and contravariant form:

$U^{\mu}=\left[\frac{1}{a}(1-\Phi),\frac{v^{i}}{a}\right]\quad
U_{\mu}=[-a(1+\Phi),av^{i}].$ (B.18)

For convenience we are going to rewrite the 3-vector of the peculiar
velocities using the Helmholtz decomposition
$\vec{v}=\vec{v}_{\parallel}+\vec{v}_{\perp}$, where
$\vec{\nabla}\times\vec{v}_{\parallel}=0$ and
$\vec{\nabla}\cdot\vec{v}_{\perp}=0$. As we stated at the beginning of this
section we are only interested in scalar perturbations, therefore, only
$\vec{v}_{\parallel}=\vec{\nabla}{v}$ plays a role in the calculations. For
the sake of convenience we list the non-null components of the total energy-
momentum tensor:

$\displaystyle T_{00}$
$\displaystyle=a^{2}\sum_{n}\left[(1+2\Phi)\bar{\rho}_{n}+\delta\rho_{n}\right]$
(B.19) $\displaystyle T_{ij}$
$\displaystyle=a^{2}\sum_{n}\left[(1+2\Psi)\bar{p}_{n}+\delta{p}_{n}\right]\delta_{ij}$
(B.20) $\displaystyle T_{0i}$
$\displaystyle=-a^{2}\partial_{i}\sum_{n}(\bar{\rho}_{n}+\bar{p}_{n})v_{n},$
(B.21)

where the subindex $n$ runs for all components considered. Once we have
computed all the required elements we can write down the perturbation
equations for the different DE models under study. Before doing so, we are
going to perform a Fourier expansion of the perturbed quantities, following
the transformation between the position-space and the momentum-space:

$f(\eta,\vec{x})=\frac{1}{(2\pi)^{3}}\int
f_{k}(\eta,\vec{k})e^{i\vec{k}\cdot\vec{x}}d^{3}k.$ (B.22)

The above expression means that we decompose all of the perturbed quantities
considered in the equations, as a sum of plane waves. Due to the fact that we
have restricted the analysis at first order, we do not have to worry about the
mixing between the different modes. They all $f_{k}(\eta,\vec{k})$ evolve
independently.

It is finally time to write down the perturbed equations in momentum space
derived from Einstein’s field equations:

  * •

$\delta{G_{00}}=8\pi{G_{N}}\delta{T_{00}}$

$k^{2}\Psi+3\mathcal{H}\left(\Psi^{\prime}-\mathcal{H}\Phi\right)=4\pi{G_{N}}a^{2}\sum_{n}\delta\rho_{n}$
(B.23)

  * •

$\delta{G_{0i}}=8\pi{G_{N}}\delta{T_{0i}}$

$\Psi^{\prime}-\mathcal{H}\Phi=4\pi{G_{N}}a^{2}\sum_{n}\left(\bar{\rho}_{n}+\bar{p}_{n}\right)v_{n}$
(B.24)

  * •

$\delta{G_{ij}}=8\pi{G_{N}}\delta{T_{ij}}$ ($i\neq j$)

$\Psi=-\Phi$ (B.25)

  * •

$\delta{G_{ij}}=8\pi{G_{N}}\delta{T_{ij}}$ ($i=j$)

$\Psi^{\prime\prime}-\Phi\left(\mathcal{H}^{2}+2\mathcal{H}^{\prime}\right)-\mathcal{H}\left(\Phi^{\prime}-2\Psi^{\prime}\right)=-4\pi{G_{N}}a^{2}\sum_{n}\delta{p}_{n}$
(B.26)

  * •

$\nabla^{\mu}T_{\mu{0}}=0$

$\sum_{n}\left[\delta\rho^{\prime}_{n}+3\mathcal{H}\left(\delta\rho_{n}+\delta{p_{n}}\right)+\left(\bar{\rho}_{n}+\bar{p}_{n}\right)\left(3\Psi^{\prime}-k^{2}v_{n}\right)\right]=0$
(B.27)

  * •

$\nabla^{\mu}T_{\mu{i}}=0$

$\sum_{n}\left[\left(\bar{\rho}^{\prime}_{n}+\bar{p}^{\prime}_{n}\right)v_{n}+\left(\bar{\rho}_{n}+\bar{p}_{n}\right)\left(4\mathcal{H}v_{n}+v^{\prime}_{n}+\Phi\right)+\delta{p_{n}}\right]=0.$
(B.28)

Although not all the equations listed are independent, they are all useful in
one way or another. If we look at (B.25) we can see that it is only valid if
we consider that the energy budget is composed just by perfect fluids. We do
know that this is not accurate if we consider, either the interaction between
photons and baryons or we deal with the full treatment of the massive
neutrinos. Taking into account that we are considering equations valid at the
matter dominated epoch, well after the decoupling time when the interaction
between photons and baryons was relevant, as a good approximation we can work
with the assumption that both of them evolve independently. As far as
neutrinos are concerned, as a first approximation we can consider that they
behave as if they were a perfect fluid. As a consequence the relation
$\Psi+\Phi\simeq 0$ remains valid.

Unfortunately it is not possible to find an analytical solution for this set
of differential equations, nevertheless, we can find solutions under some
approximations. In this section we are going to consider physical modes much
more shorter than the Hubble radius, $1/\mathcal{H}$, since they are the only
ones that can be observed. What is more, since the structure formation started
within the matter dominated epoch, well after the recombination time, we drop
the radiation contributions, namely we make the approximation
$\bar{\rho}_{r}/\bar{\rho}_{m}\simeq 0$. So, considering the before mentioned
conditions for the modes, $k^{2}\gg\mathcal{H}^{2}$, we can find the following
set of equations:

$\displaystyle
k^{2}\Psi-4\pi{G_{N}}a^{2}\left(\delta\rho_{m}+\delta\rho_{\Lambda}\right)=0$
(B.29)
$\displaystyle\sum_{i=m,\Lambda}\left[\delta\rho^{\prime}_{i}+3\mathcal{H}\left(\delta\rho_{i}+\delta{p_{i}}\right)+\left(\bar{\rho}_{i}+\bar{p}_{i}\right)\left(3\Psi^{\prime}-k^{2}v_{i}\right)\right]=0$
(B.30)
$\displaystyle\sum_{i=m,\Lambda}\left[\left(\bar{\rho}^{\prime}_{i}+\bar{p}^{\prime}_{i}\right)v_{i}+\left(\bar{\rho}_{i}+\bar{p}_{i}\right)\left(4\mathcal{H}v_{i}+v^{\prime}_{i}+\Phi\right)+\delta{p_{i}}\right]=0.$
(B.31)

As usual $\bar{\rho}_{m}=\bar{\rho}_{b}+\bar{\rho}_{dm}+\bar{\rho}^{\rm
NR}_{\nu}$ (and of course
$\delta\rho_{m}=\delta\rho_{b}+\delta\rho_{dm}+\delta\rho^{\rm NR}_{\nu}$),
consequently we get
$v_{m}=\left(\bar{\rho}_{b}{v_{b}}+\bar{\rho}_{dm}v_{dm}+\bar{\rho}^{\rm
NR}_{\nu}v^{\rm NR}_{\nu}\right)/\bar{\rho}_{m}$. Where $\rho^{\rm NR}_{\nu}$
stands for the contribution of massive neutrinos in the non-relativistic
regime. It is worthwhile to note that we make no assumptions about DE, since
we do not consider any particular value for
$w_{\Lambda}=\bar{p}_{\Lambda}/\bar{\rho}_{\Lambda}$ and we do not neglect the
contribution of $\delta\rho_{\Lambda}$. In order to compare the theoretical
predictions for a given model with the large scale structure (LSS)
observational data it turns out essential to compute the total matter density
contrast, defined as:

$\delta_{m}=\frac{\delta\rho_{b}+\delta\rho_{dm}+\delta\rho^{\rm
NR}_{\nu}}{\bar{\rho}_{m}}.$ (B.33)

It goes without saying that the particularities of the models under study,
like the equation of sate (EoS) parameter considered or some sort of possible
interactions between the different species, will affect the value of (B.33).
Therefore when the background data is not enough to distinguish between two
given models, the large scale structure (LSS) data can help by pointing out
departures in the theoretical predictions of the models at the perturbation
level.

### Perturbation equations for the scalar field models

In this section our aim is to write down the most important perturbation
equations for the non-interactive scalar field models. Throughout this thesis
we have denoted the 1-dimensional scalar field as $\Phi$, however, due to the
fact that here we are going to deal with the perturbed metric, which in the
Newtonian gauge already includes the term $\Phi$, we are forced to change
slightly the notation. Thus we refer to the scalar field and the corresponding
perturbation as:

$\Phi_{S}=\bar{\Phi}_{S}+\delta\Phi_{S}.$ (B.34)

Consequently the action at the background level must be rewritten as

$S[\bar{\Phi}_{S}]=-\int
d^{4}x\sqrt{-g}\left(\frac{1}{2}\bar{g}^{\alpha\beta}\partial_{\alpha}\bar{\Phi}_{S}\partial_{\beta}\bar{\Phi}_{S}+U(\bar{\Phi}_{S})\right).$
(B.35)

Considering (B.35), employing (B.34) and taking into account the perturbations
of the metric $g_{\mu\nu}=\bar{g}_{\mu\nu}+\delta{g}_{\mu\nu}$, we can work
out the expression for the symmetric energy-momentum tensor of a scalar field:

$T^{\Phi_{S}}_{\mu\nu}=\bar{T}^{{\Phi}_{S}}_{\mu\nu}+\delta{T}^{\Phi_{S}}_{\mu\nu}$
(B.36)

where

$\displaystyle\bar{T}^{{\Phi}_{S}}_{\mu\nu}$
$\displaystyle=\partial_{\mu}\bar{\Phi}_{S}\partial_{\nu}\bar{\Phi}_{S}+\bar{g}_{\mu\nu}\left[\frac{\bar{\Phi}^{\prime{2}}_{S}}{2a^{2}}-U(\bar{\Phi}_{S})\right]$
(B.37) $\displaystyle\delta{T}^{\Phi_{S}}_{\mu\nu}$
$\displaystyle=\partial_{\mu}\bar{\Phi}_{S}\partial_{\nu}\delta\Phi_{S}+\partial_{\mu}\delta\Phi_{S}\partial_{\nu}\bar{\Phi}_{S}+\bar{g}_{\mu\nu}\left[\frac{\bar{\Phi}^{\prime}_{S}\delta\Phi^{\prime}_{S}}{a^{2}}-\frac{\Psi\bar{\Phi}^{\prime{2}}_{S}}{a^{2}}-\frac{\partial{U(\bar{\Phi}_{S})}}{\partial\bar{\Phi}_{S}}\delta\Phi_{S}\right]$
(B.38)
$\displaystyle+\delta{g}_{\mu\nu}\left[\frac{\bar{\Phi}^{\prime{2}}_{S}}{2a^{2}}-U(\bar{\Phi}_{S})\right].$

It will be useful to write down the components of the total energy-momentum
tensor $T^{\Phi_{S}}_{\mu\nu}$:

$\displaystyle T^{\Phi_{S}}_{00}$
$\displaystyle=\frac{\bar{\Phi}^{\prime{2}}_{S}}{2}+a^{2}{U(\bar{\Psi}_{S})}+\bar{\Phi}^{\prime}_{S}\delta\Phi^{\prime}_{S}+a^{2}\frac{\partial{U(\bar{\Phi}_{S})}}{\partial\bar{\Phi}_{S}}\delta\Phi_{S}+2a^{2}\Psi{U(\bar{\Psi}_{S})}$
(B.39) $\displaystyle T^{\Phi_{S}}_{ij}$
$\displaystyle=\left[\frac{\bar{\Phi}^{\prime{2}}_{S}}{2}-a^{2}{U(\bar{\Psi}_{S})}+\bar{\Phi}^{\prime}_{S}\delta\Phi^{\prime}_{S}+a^{2}\frac{\partial{U(\bar{\Phi}_{S})}}{\partial\bar{\Phi}_{S}}\delta\Phi_{S}\right.$
(B.40)
$\displaystyle\left.+(\Psi-\Phi)\bar{\Phi}^{\prime{2}}_{S}-a^{2}\frac{\partial{U(\bar{\Phi}_{S})}}{\partial\bar{\Phi}_{S}}\delta\Phi_{S}-2a^{2}\Psi{U(\bar{\Psi}_{S})}\right]\delta_{ij}$
$\displaystyle T^{\Phi_{S}}_{0i}$
$\displaystyle=\partial_{i}\left(\bar{\Phi}^{\prime}_{S}\delta\Phi^{\prime}_{S}\right).$
(B.41)

Once we have the non-null components of the energy-momentum tensor in addition
to the components of the Einstein tensor (B.12) and the matter energy momentum
tensor $T^{m}_{\mu\nu}$ (B.19)-(B.21) (considering that here the sum over the
different components only contains the contribution of the non-relativistic
matter) we can list the Einstein equations in momentum space:

  * •

$\delta{G}_{00}=8\pi{G_{N}}\left(\delta
T^{m}_{00}+\delta{T^{\Phi_{S}}_{00}}\right)$

$k^{2}\Psi+3\mathcal{H}\left(\Psi^{\prime}-\mathcal{H}\Phi\right)=4\pi{G_{N}}a^{2}\left[\delta\rho_{m}+\frac{\bar{\Phi}^{\prime}_{S}\delta\Phi^{\prime}_{S}}{a^{2}}+\frac{\partial{U(\bar{\Phi}_{S})}}{\partial\bar{\Phi}_{S}}\delta\Phi_{S}-\frac{\bar{\Phi}^{\prime{2}}_{S}\Phi}{a^{2}}\right]$
(B.42)

  * •

$\delta{G}_{0i}=8\pi{G_{N}}\left(\delta
T^{m}_{0i}+\delta{T^{\Phi_{S}}_{0i}}\right)$

$\Psi^{\prime}-\mathcal{H}\Phi=4\pi{G_{N}}a^{2}\left((\bar{\rho}_{m}+\bar{p}_{m})v_{m}-\frac{\bar{\Phi}^{\prime}_{S}\delta\Phi_{S}}{a^{2}}\right)$
(B.43)

  * •

$\delta{G}_{ij}=8\pi{G_{N}}\left(\delta
T^{m}_{ij}+\delta{T^{\Phi_{S}}_{ij}}\right)$ ($i\neq j$)

$\Psi=-\Phi$ (B.44)

  * •

$\delta{G}_{ij}=8\pi{G_{N}}\left(\delta
T^{m}_{ij}+\delta{T^{\Phi_{S}}_{ij}}\right)$ ($i=j$)

$\Psi^{\prime\prime}-\Phi\left(\mathcal{H}^{2}+2\mathcal{H}^{\prime}\right)-\mathcal{H}\left(\Phi^{\prime}-2\Psi^{\prime}\right)=-4\pi{G_{N}}a^{2}\left(\delta{p_{m}}+\frac{\bar{\Phi}^{\prime}_{S}\delta\Phi^{\prime}_{S}}{a^{2}}-\frac{\partial{U(\bar{\Phi}_{S})}}{\partial\bar{\Phi}_{S}}\delta\Phi_{S}-\frac{\bar{\Phi}^{\prime{2}}_{S}\Phi}{a^{2}}\right)$
(B.45)

  * •

$\nabla^{\mu}{T^{\Phi_{S}}}_{\mu{0}}=0$

$\displaystyle\frac{\bar{\Phi}^{\prime{2}}_{S}}{a^{2}}\left(3\Psi^{\prime}-\Phi^{\prime}-4\mathcal{H}\Phi\right)$
$\displaystyle+\frac{\bar{\Phi}^{\prime}_{S}}{a^{2}}\left(4\mathcal{H}\delta\Phi^{\prime}_{S}-2\Phi\bar{\Phi}^{\prime\prime}_{S}+k^{2}\delta\Phi_{S}+\delta\Phi^{\prime\prime}_{S}\right)$
(B.46) $\displaystyle+$
$\displaystyle\frac{d}{d\eta}\left(\frac{\partial{U(\bar{\Phi}_{S})}}{\partial\bar{\Phi}_{S}}\right)\delta\Phi_{S}+\delta\Phi^{\prime}_{S}\left(\frac{\partial{U(\bar{\Phi}_{S})}}{\partial\bar{\Phi}_{S}}+\frac{\bar{\Phi}^{\prime\prime}_{S}}{a^{2}}\right)=0.$

It is important to note that at scales deep inside the horizon, where remember
the condition $k^{2}\gg\mathcal{H}^{2}$ is fulfilled, the perturbation of the
scalar field is strongly suppressed, $k^{2}\delta\Phi_{S}\simeq 0$ as it can
be deduced from (B.46). As a consequence the differential equation for
$\delta_{m}$ at subhorizon scales turns out to be the same as in the
$\Lambda$CDM case, with of course the replacement
$\rho^{0}_{\Lambda}\rightarrow\rho_{\Phi_{S}}$.

### Perturbation equations for the Brans-Dicke model with cosmological
constant

We are now interested in obtaining the perturbation equations for the Brans-
Dicke model when a cosmological constant is included in the action. As it can
be easily understood by taking a quick look at the equations (7.3) and (7.4)
the process will be a little longer than the previous cases, however, the
philosophy remains untouched, namely: we introduce the usual perturbations
$\delta{g}_{\mu\nu}$, $\delta\rho_{n}$, $\delta{p}_{n}$ in addition to the
$\delta\varphi$, which corresponds to the perturbation of the BD scalar field,
in the field equations and we obtain the perturbed equations up to the first
order. First of all, we provide the formulas of some expressions depending on
$\varphi$ that will be useful in subsequent computations:

$\displaystyle a^{2}\Box\varphi$
$\displaystyle=-\bar{\varphi}^{\prime\prime}-2\mathcal{H}\bar{\varphi}^{\prime}-\delta\varphi^{\prime\prime}+2\bar{\varphi}^{\prime\prime}\Phi+\nabla^{2}\delta\varphi-2\mathcal{H}\delta\varphi^{\prime}+\bar{\varphi}^{\prime}(\Phi^{\prime}-3\Psi^{\prime})+4\mathcal{H}\Phi\bar{\varphi}^{\prime},$
(B.47) $\displaystyle\partial_{\alpha}\varphi\partial^{\alpha}\varphi$
$\displaystyle=-\frac{(\bar{\varphi}^{\prime})^{2}}{a^{2}}(1-2\Phi)-\frac{2}{a^{2}}\bar{\varphi}^{\prime}\delta\varphi^{\prime},$
(B.48) $\displaystyle\nabla_{0}\nabla_{0}\varphi$
$\displaystyle=\bar{\varphi}^{\prime\prime}+\delta\varphi^{\prime\prime}-\bar{\varphi}^{\prime}(\mathcal{H}+\Phi^{\prime})-\mathcal{H}\delta\varphi^{\prime},$
(B.49) $\displaystyle\nabla_{i}\nabla_{j}\varphi$
$\displaystyle=-k_{i}{k}_{j}\delta\varphi-\bar{\varphi}^{\prime}\delta_{ij}[\mathcal{H}(1+2\Psi-2\Phi)+\Psi^{\prime}]-\mathcal{H}\delta\varphi^{\prime}\delta_{ij},$
(B.50) $\displaystyle\nabla_{i}\nabla_{0}\varphi$
$\displaystyle=i{k_{i}}(\delta\varphi^{\prime}-\bar{\varphi}^{\prime}\Phi-\mathcal{H}\delta\varphi).$
(B.51)

Since the local energy conservation equation (1.154) has exactly the same form
than the one of the standard model, the equations $\nabla^{\mu}T_{\mu{0}}=0$
and $\nabla^{\mu}T_{\mu{i}}=0$, do not change with respect to the ones
displayed in (B.27) and (B.28), respectively. The field equations, written in
the momentum space, that have been modified because of the extra d.o.f., are:

  * •

For $\mu=\nu=0$

$\begin{split}\bar{\varphi}(6\mathcal{H}\Psi^{\prime}+2k^{2}\Psi)+3\mathcal{H}^{2}\delta\varphi+&k^{2}\delta\varphi+3\mathcal{H}\delta\varphi^{\prime}+3\Psi^{\prime}\bar{\varphi}^{\prime}\\\
{}&+\frac{\omega_{\rm
BD}}{2\bar{\varphi}}\left[\frac{\delta\varphi}{\bar{\varphi}}(\bar{\varphi}^{\prime})^{2}-2\bar{\varphi}^{\prime}\delta\varphi^{\prime}\right]=8\pi
G_{N}a^{2}\sum_{n}(2\Phi\bar{\rho}_{n}+\delta\rho_{n})\,.\end{split}$ (B.52)

  * •

$\mu=i$ and $\nu=j$ ($i\neq j$)

$\Psi+\Phi=-\frac{\delta\varphi}{\bar{\varphi}}.$ (B.53)

  * •

$\mu=i$ and $\nu=j$ ($i=j$)

$\begin{split}&(\Phi^{\prime}-2\Psi^{\prime})(2\mathcal{H}\bar{\varphi}+\bar{\varphi}^{\prime})+(\Phi-\Psi)\left[2\bar{\varphi}(\mathcal{H}^{2}+2\mathcal{H}^{\prime})+2\bar{\varphi}^{\prime\prime}+2\mathcal{H}\bar{\varphi}^{\prime}+\frac{\omega_{\rm
BD}}{\bar{\varphi}}(\bar{\varphi}^{\prime})^{2}\right]\\\
{}&-2\Psi^{\prime\prime}\bar{\varphi}-\delta\varphi^{\prime\prime}-\mathcal{H}\delta\varphi^{\prime}-\delta\varphi(\mathcal{H}^{2}+2\mathcal{H}^{\prime})+\frac{\omega_{\rm
BD}}{2\bar{\varphi}}\left[\frac{\delta\varphi}{\bar{\varphi}}(\bar{\varphi}^{\prime})^{2}-2\bar{\varphi}^{\prime}\delta\varphi^{\prime}\right]=8\pi
G_{N}a^{2}\sum_{n}(2\bar{p}_{n}\Psi+\delta p_{n})\,.\end{split}$ (B.54)

  * •

The perturbed Klein-Gordon equation

$-\delta\varphi^{\prime\prime}+2\bar{\varphi}^{\prime\prime}\Phi-k^{2}\delta\varphi-2\mathcal{H}\delta\varphi^{\prime}+\bar{\varphi}^{\prime}(\Phi^{\prime}-3\Psi^{\prime})+4\mathcal{H}\Phi\bar{\varphi}^{\prime}=\frac{8\pi
G_{N}}{3+2\omega_{\rm BD}}a^{2}\sum_{n}(3\delta p_{n}-\delta\rho_{n}).$ (B.55)

As has been done for the other models, we study now the evolution of matter
perturbations at deep subhorizon scales, i.e. $k^{2}\gg\mathcal{H}^{2}$. In
this limit, Eq. (B.27) boils down to

$\delta_{m}^{\prime}+\theta_{m}+3\Psi^{\prime}=0\,,$ (B.56)

and (B.28) can be written as

$\theta^{\prime}_{m}+\mathcal{H}\theta_{m}-k^{2}\Phi=0,$ (B.57)

where we have used the definition $\theta_{m}\equiv\partial_{i}{v^{i}_{m}}$.
These equations can be easily combined to make disappear the dependence on
$\theta_{m}$. If we do that, we obtain the following approximate second order
differential equation for the matter density contrast,

$\delta_{m}^{\prime\prime}+\mathcal{H}\delta_{m}^{\prime}+k^{2}\Phi=0.$ (B.58)

The problem is now reduced to find an expression for $k^{2}\Phi$ in terms of
background quantities and $\delta_{m}$. Collecting (B.52), (B.53) and (B.55),

$-\frac{\delta\varphi}{\bar{\varphi}}=\Psi+\Phi,$ (B.59)
$2k^{2}\Psi+k^{2}\frac{\delta\varphi}{\bar{\varphi}}=\frac{8\pi
G_{N}a^{2}}{\bar{\varphi}}\bar{\rho}_{m}\delta_{m},$ (B.60)
$k^{2}\delta\varphi=\frac{8\pi G_{N}a^{2}}{3+2\omega_{\rm
BD}}\bar{\rho}_{m}\delta_{m},$ (B.61)

respectively. We can see, as expected, that for $\omega_{\rm BD}\to\infty$ and
$\bar{\varphi}=1$ the first equation above gives $\Phi=-\Psi$ (no anisotropy
stress) and the third one renders a trivial equality ($0=0$), whereas the
second equation yields Poisson equation $k^{2}\Psi=-k^{2}\Phi=4\pi
G_{N}a^{2}\sum_{n}\delta\rho_{n}$. In the general case, by combining the above
equations one finds:

$k^{2}\Phi=-\frac{4\pi
G_{N}a^{2}\bar{\rho}_{m}\delta_{m}}{\bar{\varphi}}\left(\frac{4+2\omega_{\rm
BD}}{3+2\omega_{\rm BD}}\right)\,.$ (B.62)

So, finally, inserting the previous relation in (B.58) we are led to the
desired equation for the density contrast at deep subhorizon scales:

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

or, alternatively, in terms of the cosmic time $t$,

$\ddot{\delta}_{m}+2H\dot{\delta}_{m}-\frac{4\pi
G_{N}}{\bar{\varphi}}\bar{\rho}_{m}\delta_{m}\left(\frac{4+2\omega_{\rm
BD}}{3+2\omega_{\rm BD}}\right)=0\,,$ (B.64)

with the dots denoting derivatives with respect to $t$. As expected, we can
recover the standard $\Lambda$CDM result for $\omega_{\rm BD}\to\infty$ and
$\bar{\varphi}\to 1$.
As already mentioned above, because of (B.59) non-null scalar field
perturbations induce a deviation of the anisotropic stress, $-\Psi/\Phi$, from
$1$, i.e. the GR-$\Lambda$CDM (see the notation employed in Chapter 7) value.
At scales well below the horizon,

$-\frac{\Psi}{\Phi}=\frac{1+\omega_{\rm BD}}{2+\omega_{\rm
BD}}=1-\epsilon_{\rm BD}+\mathcal{O}(\epsilon^{2}_{\rm BD})\,,$ (B.65)

so constraints on the anisotropic stress directly translate into constraints
on $\epsilon_{\rm BD}=1/\omega_{\rm BD}$, and a deviation of this quantity
from $1$ at the linear regime would be a clear signature of non-standard
gravitational physics. A model-independent reconstruction of the anisotropic
stress from observations has been recently done in [369]. Unfortunately, the
error bars are still of order $\mathcal{O}(1)$, so these model-independent
results cannot put tight constraints on $\epsilon_{\rm BD}$ (yet).

##  Appendix C Semi-analytical solutions for the Brans-Dicke model in
different epochs

Our aim in this appendix is to provide semi-analytical solutions for the
Brans-Dicke (BD) equations in the various epochs of the cosmic history, by
using a perturbative approach. The nomenclature and the notation employed here
are the same that can be found in Chapter 7. We express the BD-field (7.1) and
the scale factor up to linear order in $\epsilon_{\rm BD}$ as follows,

$\displaystyle\varphi(t)=\varphi^{(0)}+\epsilon_{\rm
BD}\varphi^{(1)}(t)+\mathcal{O}(\epsilon_{\rm BD}^{2})$ (C.1) $\displaystyle
a(t)=a^{(0)}(t)+\epsilon_{\rm BD}a^{(1)}(t)+\mathcal{O}(\epsilon_{\rm
BD}^{2})\,,$ (C.2)

with the functions with superscript $(0)$ denoting the solutions of the
background equations in standard GR with a constant Newtonian coupling that
can be in general different from $G_{N}$, and the functions with superscript
$(1)$ denoting the first-order corrections induced by a non-null
$\epsilon_{\rm BD}$. Neglecting the higher-order terms is a very good
approximation for all the relevant epochs of the expansion history due to the
small values of $\epsilon_{\rm BD}$ allowed by the data. Plugging these
expressions in Eqs. (7.3)-(7.6) we can solve the system and obtain the
dominant energy density at each epoch and the Hubble function, which of course
can also be written as

$\displaystyle\rho_{N}(t)=\rho_{N}^{(0)}(t)+\epsilon_{\rm
BD}\rho_{N}^{(1)}(t)+\mathcal{O}(\epsilon_{\rm BD}^{2})\,$ (C.3)
$\displaystyle H(t)=H^{(0)}(t)+\epsilon_{\rm
BD}H^{(1)}(t)+\mathcal{O}(\epsilon_{\rm BD}^{2})\,,$ (C.4)

respectively, where $N=R,M,\Lambda$ denotes the solution at the radiation,
matter, and $\Lambda$-dominated epochs313131For the sake of simplicity and to
ease the obtaining of analytical expressions, in this appendix we consider
three massless neutrinos.. We make use of the following relations,

$\frac{1}{2\omega_{\rm BD}+3}=\frac{\epsilon_{\rm
BD}}{2}+\mathcal{O}(\epsilon_{\rm
BD}^{2})\,,\qquad\frac{\dot{\varphi}}{\varphi}=\epsilon_{\rm
BD}\frac{\dot{\varphi}^{(1)}}{\varphi^{(0)}}+\mathcal{O}(\epsilon_{\rm
BD}^{2})\,,\qquad\frac{\omega_{\rm
BD}}{2}\left(\frac{\dot{\varphi}}{\varphi}\right)^{2}=\frac{\epsilon_{\rm
BD}}{2}\left(\frac{\dot{\varphi}^{(1)}}{\varphi^{(0)}}\right)^{2}+\mathcal{O}(\epsilon_{\rm
BD}^{2})\,.$ (C.5)

The Klein-Gordon equation is already of first order in $\epsilon_{\rm BD}$,

$\ddot{\varphi}^{(1)}+3H^{(0)}\dot{\varphi}^{(1)}=4\pi
G_{N}(\rho^{(0)}_{N}-3p^{(0)}_{N})\,,$ (C.6)

and this allows us to find $\varphi^{(1)}$ without knowing $H^{(1)}$ nor the
linear corrections for the energy densities and pressures. The Friedmann
equation leads to

$H^{(0)}=\left(\frac{8\pi{G_{N}}}{3\varphi^{(0)}}\rho_{N}^{(0)}\right)^{1/2}\,,$
(C.7)

at zeroth order, and

$6H^{(0)}H^{(1)}+3H^{(0)}\frac{\dot{\varphi}^{(1)}}{\varphi^{(0)}}-\frac{1}{2}\left(\frac{\dot{\varphi}^{(1)}}{\varphi^{(0)}}\right)^{2}=3\left(H^{(0)}\right)^{2}\left(\frac{\rho^{(1)}_{N}}{\rho_{N}^{(0)}}-\frac{\varphi^{(1)}}{\varphi^{(0)}}\right)\,,$
(C.8)

at first order, which let us to compute $a^{(1)}(t)$ once we get
$\rho^{(1)}_{N}(t,a^{(1)}(t))$ from the conservation equation and substitute
it in the r.h.s. $H^{(1)}$ is then trivially obtained using the computed scale
factor. Alternatively, one can also combine (C.8) with the pressure equation
to obtain a differential equation for $H^{(1)}$ and directly solve it.
Proceeding in this way one obtains, though, an additional integration constant
that must be fixed using the Friedmann and conservation equations.

### Radiation dominated epoch (RDE)

In the RDE the trace of the energy-momentum tensor is negligible when compared
with the total energy density in the Universe, so we can set the right-hand
side of (C.6) to zero. The leading order of the scale factor and the Hubble
function take the following form, respectively: $a^{(0)}(t)=At^{1/2}$,
$H^{(0)}(t)=1/(2t)$, with
$A\equiv(32\pi{G_{N}}\rho^{0}_{r}/3\varphi^{(0)})^{1/2}$ and $\rho^{0}_{r}$
the current value of the radiation energy density. Using these relations we
can find the BD-field as well as the scale factor at first order. They read,

$\varphi(t)=\varphi^{(0)}+\epsilon_{\rm
BD}\left(C_{R1}+C_{R2}\,t^{-1/2}\right)+\mathcal{O}(\epsilon_{\rm BD}^{2})\,,$
(C.9)

and

$a(t)=At^{1/2}\left(1+\epsilon_{\rm
BD}\left[\frac{C_{R2}^{2}\ln(t)}{24(\varphi^{(0)})^{2}{t}}-\frac{C_{R1}}{4\varphi^{(0)}}+\frac{C_{R3}}{t}\right]+\mathcal{O}(\epsilon_{\rm
BD}^{2})\right)\,,$ (C.10)

where $C_{R1}$, $C_{R2}$ and $C_{R3}$ are integration constants. Let us note
that in the RDE the evolution of the BD-field $\varphi$ is essentially frozen,
since there is no growing mode. The term evolving as $\sim t^{-1/2}$ is the
decaying mode, and after some time we are eventually left with a constant
contribution.

Finally, it is easy to find the corresponding Hubble function

$H(t)=\frac{1}{2t}\left(1+\epsilon_{\rm
BD}\left[-\frac{C_{R2}^{2}}{12(\varphi^{(0)})^{2}}\frac{\ln(t)}{t}+\frac{1}{t}\left(\frac{C_{R2}^{2}}{12(\varphi^{(0)})^{2}}-2C_{R3}\right)\right]+\mathcal{O}(\epsilon_{\rm
BD}^{2})\right)\,.$ (C.11)

At late enough times it is natural to consider that the decaying mode is
already negligible, and this allows us to simplify a lot the expressions, by
setting $C_{R2}=0$. The scalar field remains then constant in very good
approximation when radiation rules the expansion of the Universe, and the
other cosmological functions are the same as in GR ($a\sim t^{1/2}$,
$H\sim\frac{1}{2t}$), but with an effective gravitational coupling
$G=G_{N}/\varphi$.

### Matter dominated epoch (MDE)

When non-relativistic matter is dominant in the Universe the scalar field
evolves as

$\varphi(t)=\varphi^{(0)}+\epsilon_{\rm
BD}\left[\frac{C_{M1}}{t}+\frac{2\varphi^{(0)}}{3}\ln(t)+C_{M2}\right]+\mathcal{O}(\epsilon_{\rm
BD}^{2})\,,$ (C.12)

where $C_{M1}$ and $C_{M2}$ are integration constants. At leading order the
scale factor and the Hubble function take the following form,
$a^{(0)}(t)=Bt^{2/3}$ and $H^{(0)}(t)=2/3t$ respectively, with
$B\equiv\left(6\pi{G_{N}}\rho^{0}_{m}/\varphi^{(0)}\right)^{1/3}$ and
$\rho^{0}_{m}$ the current value of the matter energy density. If we neglect
the decaying mode in (C.12) we find $\varphi(a)=\varphi^{(0)}(1+\epsilon_{\rm
BD}\ln a+\mathcal{O}(\epsilon^{2}_{\rm BD}))$. This solution is also found
from the analysis of fixed points at leading order in $\epsilon_{\rm BD}$ that
will be presented in Appendix D. The scale factor reads in this case,

$a(t)=Bt^{2/3}\left(1+\epsilon_{\rm
BD}\left[-\left(\frac{C_{M1}}{\varphi^{(0)}}\right)^{2}\frac{1}{8t^{2}}-\frac{C_{M2}}{3\varphi^{(0)}}-\frac{1}{18}-\frac{2}{9}\ln(t)+\frac{C_{M3}}{t}\right]+\mathcal{O}(\epsilon_{\rm
BD}^{2})\right)\,,$ (C.13)

where $C_{M3}$ is another integration factor, and the Hubble function,

$H(t)=\frac{2}{3t}\left(1+\epsilon_{\rm
BD}\left[\left(\frac{C_{M1}}{\varphi^{(0)}}\right)^{2}\frac{3}{8t^{2}}-\frac{1}{3}-\frac{3C_{M3}}{2t}\right]+\mathcal{O}(\epsilon_{\rm
BD}^{2})\right)\,.$ (C.14)

Once more, after we neglect the contribution from the decaying modes, the
usual cosmological functions are as in GR in the MDE ($a\sim t^{2/3}$,
$H\sim\frac{2}{3t}$). However, in contrast to the RDE there is some mild
evolution (a logarithmic one with the cosmic time) of the BD-field. Since
$\varphi^{(0)}$ is obviously positive, it follows from Eq. (C.12) that the
sign of such evolution (implying growing or decreasing behaviour) is entirely
defined by the sign of the BD-parameter $\epsilon_{\rm BD}$. The fit to the
overall data clearly shows that $\epsilon_{\rm BD}<0$ and hence $\varphi$
decreases with the expansion during the MDE, which means that the effective
gravitation coupling $G=G_{N}/\varphi$ increases with the expansion.

###  $\Lambda$-dominated or VDE

Here we assume that the energy density of the Universe is completely dominated
by a vacuum fluid, with constant energy density $\rho_{\Lambda}$. This period
occurs in the very early Universe during inflation, and in the very late one
when matter has diluted significantly and a new period of inflation occurs.
The usual solution for the Hubble function in that epoch is
$H^{(0)}(t)=H_{\Lambda}$, where $H_{\Lambda}$ is a constant, fixed by (C.7).
Taking this into account we find

$\varphi(t)=\varphi^{(0)}+\epsilon_{\rm
BD}\left(2H_{\Lambda}\varphi^{(0)}t+C_{\Lambda 1}e^{-3H_{\Lambda}t}+C_{\Lambda
2}\right)+\mathcal{O}(\epsilon_{\rm BD}^{2}),$ (C.15)

and for the scale factor

$\begin{split}&a(t)=C_{\Lambda 4}e^{H_{\Lambda}t}\left(1+\epsilon_{\rm
BD}\left[-\frac{H_{\Lambda}^{2}t^{2}}{2}-\left(\frac{C_{\Lambda
1}}{\varphi^{(0)}}\right)^{2}\frac{e^{-6H_{\Lambda}t}}{8}+C_{\Lambda
3}H_{\Lambda}t\left(\frac{C_{\Lambda
2}}{2\varphi^{(0)}}+\frac{2}{3}\right)\right]+\mathcal{O}(\epsilon_{\rm
BD}^{2})\right),\end{split}$ (C.16)

where $C_{\Lambda 1},$ $C_{\Lambda 2}$, $C_{\Lambda 3}$, and $C_{\Lambda 4}$
are integration constants. Note that by performing the limit $\epsilon_{\rm
BD}\rightarrow 0$ we recover the usual GR expressions, as in all the previous
formulas. The Hubble function at first order takes the form

$\begin{split}&H(t)=H_{\Lambda}\left(1+\epsilon_{\rm
BD}\left[-\frac{2}{3}+\frac{3}{4}\left(\frac{C_{\Lambda
1}}{\varphi^{(0)}}\right)^{2}e^{-6H_{\Lambda}t}-H_{\Lambda}t-\frac{C_{\Lambda
2}}{2\varphi^{(0)}}\right]+\mathcal{O}(\epsilon_{\rm
BD}^{2})\right).\end{split}$ (C.17)

The term accompanied by the constant $C_{\Lambda 1}$ in (C.15) is a decaying
mode, which we considered to be negligible already in the RDE. Thus, we can
safely remove it, and $C_{\Lambda 2}$ can be fixed by the condition
$\varphi(t_{*})=\varphi_{*}$, at some $t_{*}$ deeply in the VDE. The scalar
field evolves then as $\varphi(t)\sim 2\epsilon_{\rm BD}H_{\Lambda}t\sim
2\epsilon_{\rm BD}\ln a$, which is the behaviour we obtain also from the
analysis of fixed points (cf. Appendix D). We can also see that during this
epoch the BD-field decreases with the expansion because $\epsilon_{\rm BD}<0$,
as indicated before.

### Mixture of matter and vacuum energy

In a Universe with a non-negligible amount of vacuum and matter energy
densities, it is also possible to obtain an analytical expression for the BD
scalar field at leading order in $\epsilon_{\rm BD}$. Unfortunately, this is
not the case for the scale factor and the Hubble function, so we will present
here only the formula for $\varphi$. In the $\Lambda$CDM the scale factor is
given by

$a^{(0)}(t)=\left(\frac{\tilde{\Omega}_{m}}{\tilde{\Omega}_{\Lambda}}\right)^{1/3}\sinh^{2/3}\left(\frac{3}{2}\sqrt{\tilde{\Omega}_{\Lambda}}H_{0}t\right)\,,$
(C.18)

so

$H^{(0)}(t)=H_{0}\sqrt{\tilde{\Omega}_{\Lambda}}\coth\left(\frac{3}{2}\sqrt{\tilde{\Omega}_{\Lambda}}H_{0}t\right)\,.$
(C.19)

Solving the Klein-Gordon equation, we obtain

$\varphi^{(1)}(t)=\sqrt{\tilde{\Omega}_{\Lambda}}H_{0}t\coth\left(\frac{3}{2}\sqrt{\tilde{\Omega}_{\Lambda}}H_{0}t\right)+\frac{2}{3}\ln\left(\sinh\left(\frac{3}{2}\sqrt{\tilde{\Omega}_{\Lambda}}H_{0}t\right)\right)\,.$
(C.20)

One can easily check that in the limits $H_{0}t\ll 1$ and $H_{0}t\gg 1$ we
recover the behaviour that we have found in previous sections for the matter
and $\Lambda$-dominated Universes, respectively. When we substitute the
previous expression in Eq. (C.1), we confirm once more that $\varphi$
decreases with the expansion (since $\epsilon_{\rm BD}<0$ and
$\dot{\varphi}^{(1)}(t)>0$ $\forall{t}$). This is the period when the Universe
is composed of a mixture of matter and vacuum energy at comparable
proportions, and corresponds to the current Universe. Thus $G=G_{N}/\varphi$
increases with the expansion in the present Universe as it was also the case
in the preceding MDE period, which is of course an important feature that
helps to solve the $H_{0}$-tension, as explained in different parts of the
Chapter 7.

### Connection of the BD-$\Lambda$CDM model with the Running Vacuum Model

As we have just seen it is possible to search for approximate solutions in the
different epochs, which can help to better understand the numerical results
and the qualitative behaviour of the BD-$\Lambda$CDM model. Actually, a first
attempt in this direction trying to show that BD-$\Lambda$CDM can mimic the
Running Vacuum Model (RVM) was done in [332, 331] and we refer the reader to
these references for details. Here we just summarize the results and adapt
them to the notation employed in this thesis. It is based on searching for
solutions in the MDE in the form of a power-law ansatz in which the BD-field
$\varphi$ evolves very slowly:

$\varphi(a)=\varphi_{0}\,a^{-\epsilon}\ \qquad(|\epsilon|\ll 1)\,.$ (C.21)

Obviously $\epsilon$ must be a very small parameter in absolute value since
$G(a)\equiv G(\varphi(a))$ cannot depart too much from $G_{N}$. On comparing
with the analysis of fixed points given in Appendix D – cf. Eq. (D.23) – we
can anticipate that $\epsilon\propto-\epsilon_{\rm BD}$, although we do not
expect perfect identification since (C.21) is a mere ansatz solution in the
MDE whereas (D.23) is an exact phase trajectory in that epoch. For
$\epsilon>0$ (hence $\epsilon_{\rm BD}<0$), the effective coupling increases
with the expansion and therefore is asymptotically free since $G(a)$ is
smaller in the past, which is the epoch when the Hubble rate (with natural
dimension of energy) is bigger. For $\epsilon<0$ ($\epsilon_{\rm BD}>0$),
instead, $G(a)$ decreases with the expansion. Using the power-law ansatz
(C.21) we find:

$\frac{\dot{\varphi}}{\varphi}=-\epsilon{H}\,,\ \ \ \ \ \ \ \ \ \ \ \
\frac{\ddot{\varphi}}{\varphi}=-\epsilon\dot{H}+\epsilon^{2}{H^{2}}.$ (C.22)

Plugging these relations into the system of equations (7.3)-(7.5) and after
some calculation it is possible to arrive at the following pair of Friedmann-
like equations to ${\cal O}(\epsilon)$[332, 331]:

$H^{2}=\frac{8\pi G}{3}\left(\rho^{0}_{m}a^{-3+\epsilon}+\rho_{\rm
DE}(H)\right)$ (C.23)

and

$\frac{\ddot{a}}{a}=-\frac{4\pi
G}{3}\,\left(\rho_{m}^{0}a^{-3+\epsilon}+\rho_{\rm
DE}(H)+3p_{\Lambda}\right)\,,$ (C.24)

with $G=G_{N}/\varphi_{0}$. The first equation emulates an effective
Friedmann’s equation with time-evolving cosmological term, in which the DE
appears as dynamical:

$\rho_{\rm DE}(H)=\rho_{\Lambda}+\frac{3\,\nu_{\rm eff}}{8\pi G}H^{2}\,.$
(C.25)

Here

$\ \nu_{\rm eff}\equiv\epsilon\left(1+\frac{1}{6}\,\omega_{\rm
BD}\epsilon\right)$ (C.26)

is the coefficient controlling the dynamical character of the dark energy
(C.25). The structure of this dynamical dark energy (DDE) is reminiscent of
the Running Vacuum Model (RVM), see [33, 370, 37] and references therein. In
the language used in Chapter 7, the analogue of the above RVM form of the
Friedmann equation can be derived from Eq. (7.35) upon taking into account
that the function ${\cal F}$ is of ${\cal O}(\epsilon_{\rm BD})$. In fact,
${\cal F}$ is the precise analogue of $\nu_{\rm eff}$, for if we set
$\epsilon\to-\epsilon_{\rm BD}$ in (C.26) it boils down to the value quoted in
Eq.(7.34). The two languages are similar, but not identical, for the reasons
explained above. Notice from (C.23) that, to ${\cal O}(\epsilon)$:

$\Omega_{m}+\Omega_{\Lambda}=1-\nu_{\rm eff}\,,$ (C.27)

so the usual sum rule of GR is slightly violated by the BD model when
parametrized as a deviation with respect to GR. Only for $\epsilon=0$ we have
$\nu_{\rm eff}=0$ and then we recover the usual cosmic sum rule. The parameter
$\nu_{\rm eff}$ becomes associated to the dynamics of the DE. Worth noticing,
the above expression adopts the form of the RVM, see [33, 370, 37] and
references therein, in particular [371, 111, 39] – where the running parameter
is usually denoted $\nu$ and is associated to the $\beta$-function of the
running vacuum. Recently, the parameter $\nu$ (and in general the structure of
the RVM) has been elucidated from direct calculations in QFT in curved space-
time within GR [42]. For additional discussions on the running of the CC term,
see e.g. [372, 373, 374]. The RVM has been shown to be phenomenologically
promising to alleviate some of the existing tensions within the $\Lambda$CDM,
particularly the $\sigma_{8}$-tension [78, 251, 233, 229, 110, 315, 201, 171,
375, 376, 377]. It is therefore not surprising that the mimicking of the RVM
by the BD-$\Lambda$CDM model enjoys of the same virtues. In actual fact, the
particular RVM form obtained in BD-gravity (we may call it “BD-RVM” for short)
is even more successful since it can cure both tensions, the $H_{0}$ and
$\sigma_{8}$ one. The reason why the BD-RVM can cure also the $H_{0}$-tension
is because we need the evolution of the effective gravitational coupling
$G_{\rm eff}$ to achieve that, as we have seen in the preview Sec. 7.2,
whereas the $\sigma_{8}$-tension can be cured with $\nu_{\rm eff}$, which is
associated to $\epsilon\propto-\epsilon_{\rm BD}$ (the second ingredient
characteristic of BD-gravity), and hence the two key elements are there to
make a successful phenomenological performance.

On the other hand, from (C.24) it follows that the EoS for the effective DDE
is

$w_{\rm eff}(z)=\frac{p_{\Lambda}}{\rho_{\rm DE}(H)}\simeq-1+\frac{3\nu_{\rm
eff}}{8\pi G\rho_{\Lambda}}\,H^{2}(z)=-1+\frac{\nu_{\rm
eff}}{\Omega_{\Lambda}}\,\frac{H^{2}(z)}{H_{0}^{2}}\,,$ (C.28)

where use has been made of (C.25). It follows that the BD-RVM, in contrast to
the original RVM, does not describe a DE of pure vacuum form
($p_{\Lambda}=-\rho_{\Lambda}$) but a DE whose EoS departs slightly from the
pure vacuum. In fact, for $\epsilon>0\ (\epsilon<0)$ we have $\nu_{\rm eff}>0\
(\nu_{\rm eff}<0)$ and the effective DDE behaves quintessence (phantom)-like.
For $\epsilon\to 0$ (hence $\nu_{\rm eff}\to 0$) we have $w_{\rm eff}\to-1$
($\Lambda$CDM). As could be expected, Eq.(C.28) is the BD-RVM version of the
effective EoS that we obtained in Sec. 7.3 – see Eq. (7.60). The two languages
are consistent. Indeed, by comparison we see that $\nu_{\rm eff}$ here plays
the role of $\Delta{\varphi}$ there. We know that
$\Delta{\varphi}=1-\varphi>0$, i.e. $\varphi<1$, for $\epsilon_{\rm BD}<0$, as
we have shown previously, which is consistent with the fact that $\nu_{\rm
eff}\propto\epsilon\propto-\epsilon_{\rm BD}>0$. Finally, since $w_{\rm eff}$
approaches $-1$ from above it corresponds to an effective quintessence
behaviour, which is more pronounced the more we explore the EoS into our past.

##  Appendix D Fixed Points in BD-$\Lambda$CDM cosmology

In order to study the fixed points of the BD theory we must define new
variables such that the system of differential equations that characterize the
model becomes of first order in the derivatives when it is rewritten in terms
of the new variables. It is useful though to firstly carry out the change
$t\to N\equiv\ln(a)$. This preliminary step will help us to identify in an
easier way how we must define the new variables. When written in terms of $N$
the system takes the following form323232Primes in this appendix stand for
derivatives w.r.t. to the variable $N=\ln a$, i.e. $()^{\prime}\equiv d()/dN$.
We consider, as in Appendix C, three massless neutrinos.:

$\frac{\psi^{\prime\prime}}{\psi}+\frac{H^{\prime}}{H}\frac{\psi^{\prime}}{\psi}+3\frac{\psi^{\prime}}{\psi}=\frac{8\pi}{3+2\omega_{\rm
BD}}\left(\frac{\rho_{m}+4\rho_{\Lambda}}{\psi H^{2}}\right)\,,$ (D.1)
$3+3\frac{\psi^{\prime}}{\psi}-\frac{\omega_{\rm
BD}}{2}\left(\frac{\psi^{\prime}}{\psi}\right)^{2}=\frac{8\pi}{\psi
H^{2}}(\rho_{r}+\rho_{m}+\rho_{\Lambda})\,,$ (D.2)
$3+\frac{H^{\prime}}{H}\left(2+\frac{\psi^{\prime}}{\psi}\right)+\frac{\psi^{\prime\prime}}{\psi}+2\frac{\psi^{\prime}}{\psi}+\frac{\omega_{\rm
BD}}{2}\left(\frac{\psi^{\prime}}{\psi}\right)^{2}=\frac{8\pi}{\psi
H^{2}}\left(\rho_{\Lambda}-\frac{\rho_{r}}{3}\right)\,,$ (D.3)
$\rho^{\prime}_{r}+4\rho_{r}=0\,,$ (D.4) $\rho^{\prime}_{m}+3\rho_{m}=0\,.$
(D.5)

Now one can define the following quantities:

$x_{\psi}\equiv\frac{\psi^{\prime}}{\psi}\quad;\quad
x_{i}^{2}\equiv\frac{8\pi\rho_{i}}{H^{2}\psi}\,,$ (D.6)

where $i=r,m,\Lambda$. In terms of these variables the system of equations can
be easily written as follows:

$x_{\psi}^{\prime}+x_{\psi}^{2}+\frac{H^{\prime}}{H}x_{\psi}+3x_{\psi}=\frac{x_{m}^{2}+4x_{\Lambda}^{2}}{3+2\omega_{\rm
BD}}\,,$ (D.7) $3+3x_{\psi}-\frac{\omega_{\rm
BD}}{2}x_{\psi}^{2}=x_{r}^{2}+x_{m}^{2}+x_{\Lambda}^{2}\,,$ (D.8)
$3+\frac{H^{\prime}}{H}\left(2+x_{\psi}\right)+x_{\psi}^{\prime}+x_{\psi}^{2}+2x_{\psi}+\frac{\omega_{\rm
BD}}{2}x_{\psi}^{2}=x_{\Lambda}^{2}-\frac{x_{r}^{2}}{3}\,,$ (D.9)
$x_{r}^{\prime}=-x_{r}\left(\frac{H^{\prime}}{H}+2+\frac{x_{\psi}}{2}\right)\,,$
(D.10)
$x_{m}^{\prime}=-x_{m}\left(\frac{H^{\prime}}{H}+\frac{3}{2}+\frac{x_{\psi}}{2}\right)\,.$
(D.11)

This system is of first order, as wanted. We have five equations and five
unknowns, namely: $x_{r},x_{m},x_{\psi},x_{\Lambda},H^{\prime}/H$. We can
reduce significantly the complexity of the system if we just isolate
$H^{\prime}/H$ from (D.7) and $x_{\Lambda}$ from (D.8),

$\frac{H^{\prime}}{H}=\frac{1}{x_{\psi}}\left[\frac{x_{m}^{2}+4x_{\Lambda}^{2}}{3+2\omega_{\rm
BD}}-x_{\psi}^{\prime}-3x_{\psi}-x_{\psi}^{2}\right]\,,$ (D.12)
$x_{\Lambda}^{2}=3+3x_{\psi}-\frac{\omega_{\rm
BD}}{2}x_{\psi}^{2}-x_{r}^{2}-x_{m}^{2}\,,$ (D.13)

and substitute the resulting expressions in the other equations, i.e. in
(D.9), (D.10), and (D.11). Doing this, and after a little bit of algebra, one
finally obtains three equations written only in terms of
$x_{r},x_{m},x_{\psi}$:

$x_{\psi}^{\prime}=-x_{\psi}\left[3+3x_{\psi}-\frac{1}{2}\omega_{\rm
BD}x_{\psi}^{2}-\frac{2}{3}x_{r}^{2}-\frac{x_{m}^{2}}{2}-\left(\frac{1}{x_{\psi}}+\frac{1}{2}\right)\left(\frac{12+12x_{\psi}-2\omega_{\rm
BD}x_{\psi}^{2}-4x_{r}^{2}-3x_{m}^{2}}{3+2\omega_{\rm BD}}\right)\right]\,,$
(D.14) $x_{r}^{\prime}=-x_{r}\left[2+\frac{5}{2}x_{\psi}-\frac{\omega_{\rm
BD}}{2}x_{\psi}^{2}-\frac{2}{3}x_{r}^{2}-\frac{x_{m}^{2}}{2}-\frac{1}{2}\left(\frac{12+12x_{\psi}-2\omega_{\rm
BD}x_{\psi}^{2}-4x_{r}^{2}-3x_{m}^{2}}{3+2\omega_{\rm BD}}\right)\right]\,,$
(D.15)
$x_{m}^{\prime}=-x_{m}\left[\frac{3}{2}+\frac{5}{2}x_{\psi}-\frac{\omega_{\rm
BD}}{2}x_{\psi}^{2}-\frac{2}{3}x_{r}^{2}-\frac{x_{m}^{2}}{2}-\frac{1}{2}\left(\frac{12+12x_{\psi}-2\omega_{\rm
BD}x_{\psi}^{2}-4x_{r}^{2}-3x_{m}^{2}}{3+2\omega_{\rm BD}}\right)\right]\,.$
(D.16)

They allow us to search for the fixed points of the system. There is an
important restriction produced by (D.10) and (D.11). Supposing that $x_{r}\neq
0$ and $x_{m}\neq 0$, we see that the mentioned equations impose that

$\frac{x_{r}^{\prime}}{x_{r}}=-\frac{x_{m}^{\prime}}{x_{m}}+\frac{1}{2}.$
(D.17)

This equation is not compatible with the conditions of fixed point, so that we
should assume that $x_{r}=0$, $x_{m}=0$ or both conditions at the same time.
The fixed points are:

RDE

$(x_{r},x_{m},x_{\Lambda},x_{\psi})_{\rm RD}=\left(\sqrt{3},0,0,0\right)\,,$
(D.18)

the Jacobian of the non-linear system (D.14), (D.15), (D.16) has eigenvalues
$\lambda^{\rm RD}_{1}=-1$, $\lambda^{\rm RD}_{2}=1/2$, $\lambda^{\rm
RD}_{3}=4$ so that is an unstable point.

MDE

$(x_{r},x_{m},x_{\Lambda},x_{\psi})_{\rm
MD}=\left(0,\frac{\sqrt{12+17\omega_{\rm BD}+6\omega_{\rm
BD}^{2}}}{\sqrt{2}|1+\omega_{\rm BD}|},0,\frac{1}{1+\omega_{\rm
BD}}\right)\,,$ (D.19)

the Jacobian of the non-linear system (D.14), (D.15), (D.16) has eigenvalues
$\lambda^{\rm MD}_{1}=-1/2,$ $\lambda^{\rm MD}_{2}\simeq-3/2,$ $\lambda^{\rm
MD}_{3}\simeq 3$ so that is also an unstable point.

$\Lambda$-dominated or VDE

$(x_{r},x_{m},x_{\Lambda},x_{\psi})_{\rm\Lambda
D}=\left(0,0,\frac{\sqrt{15+28\omega_{\rm BD}+12\omega_{\rm
BD}^{2}}}{|1+2\omega_{\rm BD}|},\frac{4}{1+2\omega_{\rm BD}}\right)\,$ (D.20)

the Jacobian of the non-linear system (D.14), (D.15), (D.16) has eigenvalues
$\lambda^{\Lambda}_{1}=-2,$ $\lambda^{\Lambda}_{2}=-3/2,$
$\lambda^{\Lambda}_{3}\simeq-3$ so that is a stable point.

The first and second fixed points are unstable, whereas the latter is stable.
We have assumed that $|\omega_{\rm BD}|\gg 1$ for approximating the
eigenvalues. The first one is very well-known, since regardless of the initial
conditions for the BD scalar field we already know that the velocity of the
scalar field decays during the RDE, and the solution tends to the attractor
with $\psi^{\prime}=0$ ($x_{\psi}=0$) and full domination of radiation, i.e.
$x_{r}=\sqrt{3}$. The Hubble function and BD scalar field are:

$\psi_{\rm RD}(a)=\psi_{r}\,,$ (D.21) $H^{2}_{\rm
RD}(a)=\frac{8\pi}{3\psi_{r}}\rho^{0}_{r}a^{-4}\,,$ (D.22)

where $\psi_{r}$ is an arbitrary constant and $\rho^{0}_{r}$ is the value of
the radiation energy density at present. This fixed point is unstable because
at some moment non-relativistic matter starts to dominate the expansion. When
this happens the solution starts to look for the new attractor, the one of the
MDE. We stress that this is an exact solution. The BD scalar field and the
Hubble function during the MDE take the following form:

$\psi_{\rm MD}(a)=Ca^{\frac{1}{1+\omega_{\rm BD}}}\,,$ (D.23) $H^{2}_{\rm
MD}(a)=\frac{16\pi\rho_{m}^{0}(1+\omega_{\rm
BD})^{2}a^{-\left(\frac{4+3\omega_{\rm BD}}{1+\omega_{\rm
BD}}\right)}}{C(12+17\omega_{\rm BD}+6\omega_{\rm BD}^{2})}\,,$ (D.24)

where $C$ is an arbitrary constant and $\rho_{m}^{0}$ is the value of the
matter energy density at present. The MD fixed point is, again, unstable,
because the MDE finishes and the VDE starts. The solution searches now for the
last fixed point, which is stable. In this last case, we have obtained the
following Hubble function and BD field,

$\psi_{\Lambda}(a)=Da^{4/(1+2\omega_{\rm BD})}\,,$ (D.25)
$H^{2}_{\Lambda}(a)=\frac{8\pi\rho_{\Lambda}}{D}\frac{(1+2\omega_{\rm
BD})^{2}}{15+28\omega_{\rm BD}+12\omega_{\rm BD}^{2}}a^{-4/(1+2\omega_{\rm
BD})}\,,$ (D.26)

where $D$ is an arbitrary constant and $\rho_{\Lambda}$ is the constant value
of vacuum energy.

One can easily check that the solutions computed in Appendix C of first order
in $\epsilon_{\rm BD}=1/\omega_{\rm BD}$, coincide (once the decaying modes
become irrelevant) with the ones presented here when the latter are Taylor-
expanded up to first order in this parameter as well. The results we have
found here are consistent with previous studies on fixed points in
cosmological dynamical systems, see e.g. [378] and references therein.

##  Appendix E Cosmological data and statistical analysis

If one desires to carry out a serious study of a given cosmological model it
is fundamental to go beyond the theoretical predictions, and put them in the
light of the increasing amount of observational data at our disposal, allowing
then to test their performance and make a model comparison, with the ultimate
purpose of discerning which is the most favoured model. In this appendix we
explicitly show the theoretical definitions for the various observable used in
the analysis, necessary to compare with the experimental data. Additionally,
we also collect some basic (but technical) information about the fitting
procedure that has been employed to get the results displayed in Chapters 2-7.

Let us start by defining the $\chi^{2}$ function which represents the building
block of the fitting procedure:

$\chi^{2}_{\rm
tot}(\vec{p})=\sum_{s}[\vec{x}_{s}(\vec{p})-\vec{d}_{s}]^{T}C_{s}^{-1}[\vec{x}_{s}(\vec{p})-\vec{d}_{s}]\,,$
(E.1)

where the subindex $s$ runs over all the data sets considered in the analysis,
for instance the CMB, BAO or SNIa 333333We must remark that actually the form
of (E.1) does not apply for the Strong Lensing data, as it will be explained
in the corresponding section. . Later on, we will go into the details of the
different data considered. The $\vec{x}_{s}$ are the vectors containing the
theoretical predictions extracted from the different models under study,
$\vec{p}$ is the vector containing the free parameters, they represent the
different d.o.f. that characterize a certain model. Finally $C_{s}$ is the
covariance matrix for each data set $s$, which can be constructed from the
correlation matrix $\rho_{s}$ and the $1\sigma$ uncertainties as follows:
$C_{s,ij}=\rho_{s,ij}\sigma_{s,i}\sigma_{s,j}$.

As we shall see, it is fundamental to properly build the function
$\chi^{2}_{\rm tot}(\vec{p})$. One must be very careful to only pick up data
for which either the correlations with other data are well under control or
there are no correlations at all. Another issue that must be considered it is
the possibility of the double counting. Throughout the last years, many data
releases from the different surveys have appeared and as it can be understood
they are not independent since they have been obtained considering at least an
overlapping sample. Let us specify how to construct the $\chi^{2}$ function by
treating all the data sets employed in this thesis one by one.

### Cosmological data

In the following we comment on the most important details of the observational
data employed in the different chapters, as well as, the cosmological
observable that allow us to contrast them with the theoretical predictions.
While it is true that we do not take into consideration all the possible
sources of cosmological data, we do use, to the best of our knowledge, the
most reliable ones. It is important to remark that the extraction of data,
through cosmological observations, is a task far from be a piece of cake. A
great effort has been made, over the years, to perfect the employed
techniques. This is the reason why we prefer to be conservative when adding
new data to the dataset.

#### SNIa

As we have mentioned in the introductory chapter supernovae are extremely
luminous explosions that can be divided in different types according to their
spectrum. They played a fundamental role in the confirmation of the
accelerated expansion of the Universe and they are one of the most
consolidated sources of cosmological data at our disposal. For the SNIa
sector, the fitted quantity is the distance modulus, defined as:

$\mu(z,\vec{p})=5\log{D_{L}(z,\vec{p})}+M$ (E.2)

where $M$ is a normalization parameter and $D_{L}(z,\vec{p})$ is the
luminosity distance:

$D_{L}(z,\vec{p})={c}{(1+z)}\int_{0}^{z}\frac{{d}z^{\prime}}{H(z^{\prime},\vec{p})}\;,$
(E.3)

with $c$ being the speed of light (now included explicitly in this formula for
better clarity). The previous formula (E.3) for the luminosity distance
applies only for a spatially flat Universe, which is something that we have
assumed throughout this thesis.

There are two different options to deal with $M$ usually referred to as
nuisance parameter: it can be included as if it were just another free
parameter in the fitting procedure or it can be integrated out through the
corresponding marginalization of the SNIa likelihood. This integration can be
performed analytically, so this also helps to improve the computational
efficiency by reducing in one dimension the fitting parameter space. Let us
provide the details of the analytical integration.

After this marginalization, the resulting effective $\chi^{2}$ function for
this data set reads,

$\chi^{2}_{\rm
SNIa}(\vec{p})=\vec{y}(\vec{p})^{T}J\vec{y}(\vec{p})-\frac{[\sum_{i,j=1}^{N}J_{ij}y_{i}(\vec{p})]^{2}}{\sum_{i,j=1}^{N}J_{ij}}\,,$
(E.4)

where $J\equiv C_{\rm SNIa}^{-1}$, $y_{i}(\vec{p})\equiv
5\log{D_{L}(z_{i},\vec{p})}-{\mu}_{\rm obs,i}$ and $N$ is the number of data
points for each case. (E.4) is the direct generalization of formula (13.16) of
[55], since we are now considering a non-diagonal covariance matrix, instead
of a diagonal one. Let us prove (E.4). The original SNIa likelihood can be
written as follows,

$\mathcal{L}_{\rm
SNIa}(\vec{p},M)=\mathcal{L}_{0}\,e^{-\frac{1}{2}[y_{i}(\vec{p})+M]J_{ij}[y_{j}(\vec{p})+M]}\,,$
(E.5)

where we have assumed the Einstein summation notation in order to write the
result more concisely and $\mathcal{L}_{0}$ is the normalization constant. Let
us now marginalize this distribution over the parameter $M$,

$\tilde{\mathcal{L}}(\vec{p})\equiv\int_{-\infty}^{+\infty}\mathcal{L}_{\rm
SNIa}(\vec{p},M)\,dM\,.$ (E.6)

This can be expressed like

$\tilde{\mathcal{L}}(\vec{p})={\mathcal{L}}_{0}e^{-\frac{1}{2}B_{0}(\vec{p})}F(\vec{p})\,,$
(E.7)

with

$F(\vec{p})\equiv\int_{-\infty}^{+\infty}e^{-\frac{B_{1}}{2}M^{2}-MB_{2}(\vec{p})}\,dM\,,$
(E.8)

and

$\begin{array}[]{rcl}B_{0}(\vec{p})&\equiv&\sum_{i,j=1}^{31}y_{i}(\vec{p})J_{ij}y_{j}(\vec{p})\,,\\\
B_{1}&\equiv&\sum_{i,j=1}^{31}J_{ij}\,,\\\
B_{2}(\vec{p})&\equiv&\sum_{i,j=1}^{31}J_{ij}y_{i}(\vec{p})\,.\end{array}$
(E.9)

The problem is basically reduced to compute the function $F(\vec{p})$
appearing in (E.7). Let us perform the derivative of the latter with respect
to $M$ and the ith component of the fitting vector $\vec{p}$,

$\frac{\partial F}{\partial
M}=0=-\int_{-\infty}^{+\infty}\left[MB_{1}+B_{2}({\vec{p}})\right]e^{-\frac{B_{1}}{2}M^{2}-MB_{2}(\vec{p})}\,dM\,,$
(E.10) $\frac{\partial F}{\partial p_{i}}=-\frac{\partial
B_{2}(\vec{p})}{\partial
p_{i}}\int_{-\infty}^{+\infty}M\,e^{-\frac{B_{1}}{2}M^{2}-MB_{2}(\vec{p})}\,dM\,.$
(E.11)

Combining these two expressions we can obtain a differential equation for
$F(\vec{p})$, which is very easy to solve,

$\frac{1}{F(\vec{p})}\frac{\partial F}{\partial p_{i}}=\frac{\partial
B_{2}}{\partial p_{i}}\frac{B_{2}(\vec{p})}{B_{1}}\,,$ $F(\vec{p})\propto
e^{\frac{B_{2}^{2}(\vec{p})}{2B_{1}}}\,.$ (E.12)

The integration constant is irrelevant since it can be absorbed by the
normalization factor $\mathcal{L}_{0}$ in (E.7), giving rise to a new one,
$\tilde{\mathcal{L}}_{0}$. Thus, the marginalized likelihood reads,

$\tilde{\mathcal{L}}(\vec{p})=\tilde{\mathcal{L}}_{0}\,e^{-\frac{\chi^{2}_{\rm
SNIa}(\vec{p})}{2}}\,,$ (E.13)

with $\chi^{2}_{\rm SNIa}(\vec{p})$ given by (E.4).

#### BAO

Galaxy surveys measure angular and redshift distributions of galaxies. Using
these data, they are able to obtain the matter power spectrum
$P_{m}(z,\vec{k})$ in redshift space upon modeling the bias factor and
transforming these angles and redshifts into distances. However, the distance
of two sources along the line of sight, being given by the difference of
redshift $\Delta z$, depends on the cosmological model. Similarly, the
distance of two sources in the direction perpendicular to the line of sight is
given by the angular separation, $\Delta\theta$, and also depends on the
cosmological model. Therefore, it is convenient to first define a fiducial
model to which one refers the predictions of the new models being studied.
Usually such model is the $\Lambda$CDM with appropriately chosen values of the
free parameters. One also has to disentangle the effect of the redshift space
distortions that are due to the peculiar velocities of galaxies falling into
the gravitational potential wells. The comoving wavevector $\vec{k}$ is
usually decomposed in a parallel to the line of sight component,
$k_{\parallel}$, and another one which is perpendicular to it, $k_{\perp}$.
These two directions determine the two-dimensional (2D) space in which the
distribution of galaxies is studied. The analysis of the BAO signal in the 2D
power spectrum carries precious information about the angular and longitudinal
size of the BAO standard ruler at the measured redshift $z$, namely

$\Delta\theta_{s}(z)=\frac{ar_{s}(a_{d})}{D_{A}(a)}=\frac{r_{s}(z_{d})}{(1+z)D_{A}(z)}\,,$
(E.14)

where $a=1/(1+z)$ is the scale factor, and

$\Delta z_{s}(z)=\frac{r_{s}(z_{d})H(z)}{c}\,,$ (E.15)

$r_{s}(z_{d})$ being the sound horizon at the redshift of the drag epoch (see
below), with $D_{A}$ the proper diameter angular distance,

$D_{A}(z,\vec{p})\equiv\frac{c}{(1+z)}\int^{z}_{0}\frac{dz^{\prime}}{H(z^{\prime},\vec{p})}\,.$
(E.16)

In principle, galaxy surveys are able to extract anisotropic BAO information
from their analyses. That is to say, they can grant constraints on the
quantities $r_{s}(z_{d})/D_{A}(z)$ and $H(z)r_{s}(z_{d})$ through the
perpendicular and parallel dilation scale factors

$\displaystyle\alpha_{\perp}$ $\displaystyle\equiv$ $\displaystyle\frac{r^{\rm
fid}_{s}(z_{d})/D^{\rm
fid}_{A}(z)}{r_{s}(z_{d})/D_{A}(z)}=\frac{D_{A}(z)r^{\rm
fid}_{s}(z_{d})}{D^{\rm fid}_{A}(z)r_{s}(z_{d})}$
$\displaystyle\alpha_{\parallel}$ $\displaystyle\equiv$
$\displaystyle\frac{H^{\rm fid}(z)r^{\rm
fid}_{s}(z_{d})}{H(z)r_{s}(z_{d})}\,,$ (E.17)

where $r^{\rm fid}_{s}(z_{d})$, $D^{\rm fid}_{A}(z)$ and $H^{\rm fid}(z)$ are
obtained in the fiducial cosmology that the survey has used to convert
redshift and angles into distances. The quantity $r_{s}(z_{d})/D_{A}(z)$ is
the sound horizon at the drag epoch measured in units of the angular diameter
distance at redshift $z$, and similarly $r^{\rm fid}_{s}(z_{d})/D^{\rm
fid}_{A}(z)$ is the corresponding value in the fiducial model. These ratios
are just the comoving angular sizes (E.14), i.e. the angular sizes divided by
the scale factor. Similarly, $H(z)r_{s}(z_{d})$ and $H^{\rm fid}(z)r^{\rm
fid}_{s}(z_{d})$ are, respectively, the sound horizon in units of the Hubble
horizon in the given model and in the fiducial model.

Sometimes the surveys are not equipped with enough BAO data to obtain (E.14)
and (E.15) separately. In this case they limit themselves to compute the
volume-averaged spectrum in a conveniently defined volume $V$ and obtain a
constraint which is usually encapsulated in the isotropic BAO estimator
$r_{s}(z_{d})/D_{V}(z)$ (sometimes denoted $d_{z}$ [194]). Here $D_{V}(z)$ is
an effective distance or dilation scale [379] obtained from the cubic root of
the volume $V$ (whose value is defined from the square of the transverse
dilation scale times the radial dilation scale):

$D_{V}(z,\vec{p})=\left[zD_{H}(z,\vec{p})D^{2}_{M}(z,\vec{p})\right]^{1/3}\,.$
(E.18)

In this expression, $D_{M}=(1+z)D_{A}$ is the comoving angular diameter
distance (playing the role of the transverse dilation scale at redshift $z$)
and $D_{H}$ is the Hubble radius

$D_{H}(z,\vec{p})\equiv\frac{c}{H(z,\vec{p})},$ (E.19)

with $zD_{H}$ acting as the radial dilation scale at redshift $z$. Therefore,
the distilled BAO estimator $d_{z}$ measures the sound horizon distance at the
drag epoch in units of the effective dilation scale (E.18).

Again a comparison with a fiducial value is necessary. The isotropic BAO
parameter relating the fiducial model value of $D_{V}$ with the actual value
of a given model is defined in [153] as

$\alpha\equiv\frac{D_{V}(z)}{D^{\rm fid}_{V}(z)}\,.$ (E.20)

In other cases $\alpha$ is defined in terms of $d_{z}$,

$\alpha\equiv\frac{r^{\rm fid}_{s}(z_{d})/D^{\rm
fid}_{V}(z)}{r_{s}(z_{d})/D_{V}(z)}=\frac{D_{V}(z)r^{\rm
fid}_{s}(z_{d})}{D^{\rm fid}_{V}(z)r_{s}(z_{d})}\,,$ (E.21)

as in [154, 155]. These isotropic BAO estimators (E.20 and E.21) are akin to
the Alcock-Paczynski (AP) test [380], in which the ratio of the observed
radial/redshift to the angular size at different redshifts, $\Delta
z_{s}/\Delta\theta_{s}$, is used to obtain cosmological constraints on the
product $H(z)D_{M}(z)$ and hence of $H(z)D_{A}(z)$. Both for the AP test and
the isotropic BAO measurement, the value of the Hubble parameter $H(z)$ cannot
be disentangled from $D_{A}(z)$, there is a degeneracy. This is in
contradistinction to the situation with the anisotropic BAO parameters (E), in
which a measurement of both $\alpha_{\perp}$ and $\alpha_{\parallel}$ permits
to extract the individual values of $H(z)$ and $D_{A}(z)$.

In our fitting analysis, carried out in the different chapters of this thesis,
we use data of both types of BAO estimators: isotropic and anisotropic ones.
The anisotropic BAO data contains more information than the isotropic one
because the aforementioned degeneracy between $H(z)$ and $D_{A}(z)$ has been
broken. Thus, anisotropic BAO is richer and yields stronger cosmological
constraints [381, 382].

To compute the theoretical predictions for the above BAO estimators, we need
some extra formulas. For instance, the sound horizon at redshift $z$, i.e.
$r_{s}(z)$, is given by the expression

$r_{s}(z)=\int^{\infty}_{z}\frac{c_{s}(z^{\prime})dz^{\prime}}{H(z^{\prime})}\,,$
(E.22)

where

$c_{s}(z)=\frac{c}{\sqrt{3(1+\mathcal{R}(z)})}\,$ (E.23)

is the sound speed in the photo-baryon plasma and
$\mathcal{R}(z)=\frac{\delta\rho_{b}(z)}{\delta\rho_{\gamma}(z)}$. If the
expressions for the baryon and photon energy densities take the usual form,
then we have $\mathcal{R}(z)=(3\Omega_{b}/4\Omega_{\gamma})/(1+z)$, being
$\Omega_{b}$ and $\Omega_{\gamma}$ the baryon and photon density parameter at
present time, respectively. The redshift at the drag epoch is
$z_{d}\sim\mathcal{O}(10^{3})$. Below such redshift value the baryon
perturbations effectively decouple from the photon ones and start to grow with
dark matter perturbations. The precise value of $z_{d}$ depends in a
complicated way on the different cosmological parameters. The fitted formula
obtained in [177], can help us get an idea of how $z_{d}$ depends on the
standard parameters:

$\displaystyle
z_{d}=\frac{1291\,\omega_{m}^{0.251}}{1+0.659\,{\omega_{m}}^{0.828}}\left[1+\beta_{1}{\omega_{b}}^{\beta_{2}}\right]$
$\displaystyle\beta_{1}=0.313\,{\omega_{m}}^{-0.419}\left[1+0.607\,{\omega_{m}}^{0.674}\right]$
(E.24) $\displaystyle\beta_{2}=0.238\,{\omega_{m}}^{0.223}\,.$

Let us point out that many experimental groups do not make use of these fitted
formulas (denoted with a subindex $ff$ below), but of the complete set of
Boltzmann (Boltz) equations, which must be solved numerically. Both approaches
can give rise to differences that are around $2-3\%$. Consequently, if the
observational values are computed with a Boltzmann code and one is interested
in using the fitted formulas to get the theoretical observable, then it is
necessary to apply a re-scaling in order to perform the comparison with our
theoretical predictions. In these cases we follow the same procedure applied
in [155] upon re-scaling the observational data by $f_{r}\equiv
r_{s}(z_{d})^{\rm fid}_{ff}/r_{s}(z_{d})^{\rm fid}_{{\rm Boltz}}$. The two
quantities involved in this ratio are computed using the same fiducial
$\Lambda$CDM cosmology chosen by the observational teams, but while the first
is obtained with the fitted formulas (E.22)-(E), the second is obtained with a
Boltzmann code. One can also compute the latter using the existing
approximated formulas for the sound horizon at $z_{d}$ that are obtained by
fitting the data extracted from e.g. the CAMB code, since they are very
accurate, with errors less than the 0.1% [382],

$r_{s}(z_{d})_{\rm Boltz}=\frac{55.234(1+\Omega_{\nu}h^{2})^{-0.3794}\,{\rm
Mpc}}{(\Omega_{dm}h^{2}+\Omega_{b}h^{2})^{0.2538}(\Omega_{b}h^{2})^{0.1278}}\,,$
(E.25)

where $\Omega_{\nu}h^{2}$ is the neutrino reduced density parameter and
depends, of course, on the effective number of neutrino species and the photon
CMB temperature, which is taken from [56].

We point out that the aforementioned re-scaling does not change the values of
the dilation scale factors (E), (E.20) and (E.21) furnished by the
experimental teams. If they, for instance, hand over their results through the
product $\alpha_{\perp}(D^{\rm fid}_{A}/r_{s}(z_{d})^{\rm fid}_{{\rm
Boltz}})$, then we just multiply it by $f_{r}^{-1}$ in order to keep our
fitting procedure consistent. For the same reason, we re-scale $(H^{\rm
fid}r_{s}(z_{d})^{\rm fid}_{{\rm Boltz}})/\alpha_{\parallel}$ multiplying it
by $f_{r}$. But we insist that the dilation scale factors, which are the
fundamental outputs extracted from the analysis of the BAO signal, are not
modified by this re-scaling.

In the different analysis considered in this thesis the BAO data has been
updated and corrected (when necessary) as the experimental teams have made it
public. Take a look at the Chapters 2-7 in order to know the particular data
points employed in each case.

#### Cosmic Chronometers

Data on the Hubble parameter at different redshift points. We use $31$ (in
Chapters 2-4 we employ 30 but when the data point from [270] was released we
incorporated it into our table) data points on $H(z_{i})$ at different
redshifts. We use only $H(z_{i})$ values obtained by the so-called
differential-age (or cosmic chronometer) techniques applied to passively
evolving galaxies. That is to say, we use data where one estimates $H(z)$ from
$(1+z)H(z)=-dz/dt$, where $dz/dt$ is extracted from a sample of passive
galaxies (i.e. with essentially no active star formation) of different
redshifts and ages. See the corresponding references for more details. The
important point to remark here is that these $H(z_{i})$ values, although
relying on the theory of spectral evolution of galaxies, are uncorrelated with
the BAO data points.

The covariance matrix for the $H(z_{i})$ data is diagonal and therefore,
$C_{H,ij}=\sigma_{H,i}^{2}\delta_{ij}$. The possibility of having non-zero
off-diagonal terms in the correlation matrix is not excluded, but these
coefficients are not found in the literature. Thus, the corresponding
$\chi^{2}$ function adopts the following simple form,

$\chi^{2}_{H}(\vec{p})=\sum_{i=1}^{31}\left[\frac{H(z_{i},\vec{p})-H_{\rm
obs}(z_{i})}{\sigma_{H,i}}\right]^{2}\,.$ (E.26)

Where $H(z_{i},\vec{p})$ represents the theoretical value, for a given model
which depends on the parameters encapsulated in $\vec{p}$, evaluated at each
redshift $z_{i}$ and $H_{\rm obs}(z_{i})$ are the observational values.

#### Large Scale Structure

By observing the existing anisotropies in the clustering of galaxies (in the
redshift space) produced by the peculiar velocities, it is possible to obtain
measurements that will allow us to test the growth rate of structure of a
given model. An observable that turns out to be particularly useful, since it
does not contain the bias factor which is model dependent, is the combination
$f(a,\vec{k})\sigma_{8}(a)$. The term $f(a,\vec{k})$ is known as the growth
rate and it can be computed from the following expression:

$f(a,\vec{k})=\frac{a}{\delta_{m}(a,\vec{k})}\frac{d\delta_{m}(a,\vec{k})}{da},$
(E.27)

where $\delta_{m}=\delta_{m}/\rho_{m}$ is the density contrast for non-
relativistic matter. If we are well within the horizon ($k\gg aH$) and it is
valid the linear regime it is possible to express the matter density contrast
as:

$\delta_{m}(a,\vec{k})=D(a)F(\vec{k})$ (E.28)

where as it can be appreciated the dependence on $\vec{k}$ has been factored
out. The properties of $F(\vec{k})$ are determined by the initial conditions
and $D(a)$ is called the growth function. The above relation is going to be
useful due to the fact that we want to make use of the matter power spectrum
(which is related with the matter density contrast as we are about to see) to
obtain the value of $f(a)$. The relation between the matter power spectrum and
the density contrast reads as follows:

$P_{m}(a,\vec{k})=C\langle\delta_{m}(a,\vec{k})\delta^{*}_{m}(a,\vec{k})\rangle=CD^{2}(a)\langle
F(\vec{k})F^{*}(\vec{k})\rangle\equiv D^{2}(a)P_{0}(\vec{k})$ (E.29)

where $C$ is a constant and $P_{0}(\vec{k})=C\langle
F(\vec{k})F^{*}(\vec{k})\rangle$ is the primordial power spectrum which is