# Cosmological effects of the Galileon term in scalar-tensor Theories

Angelo G. Ferrari<EMAIL_ADDRESS>Dipartimento di Fisica e
Astronomia, Università di Bologna, viale Berti Pichat 6/2, 40127 Bologna,
Italy.  INFN, Sezione di Bologna, viale C. Berti Pichat 6/2, 40127 Bologna,
Italy  Mario Ballardini<EMAIL_ADDRESS>Dipartimento di Fisica e
Scienze della Terra, Università degli Studi di Ferrara, via Giuseppe Saragat
1, 44122 Ferrara, Italy  INFN, Sezione di Ferrara, via Giuseppe Saragat 1,
44122 Ferrara, Italy  INAF/OAS Bologna, via Piero Gobetti 101, 40129 Bologna,
Italy  Fabio Finelli<EMAIL_ADDRESS>INAF/OAS Bologna, via Piero
Gobetti 101, 40129 Bologna, Italy  INFN, Sezione di Bologna, viale C. Berti
Pichat 6/2, 40127 Bologna, Italy  Daniela Paoletti<EMAIL_ADDRESS>INAF/OAS Bologna, via Piero Gobetti 101, 40129 Bologna, Italy  INFN, Sezione
di Bologna, viale C. Berti Pichat 6/2, 40127 Bologna, Italy  Nicoletta Mauri
<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia, Università di
Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy.  INFN, Sezione di
Bologna, viale C. Berti Pichat 6/2, 40127 Bologna, Italy

###### Abstract

We study the cosmological effects of a Galileon term in scalar-tensor theories
of gravity. The subset of scalar-tensor theories considered are characterized
by a nonminimal coupling $F(\sigma)R$, a kinetic term with arbitrary sign
$Z(\partial\sigma)^{2}$ with $Z=\pm 1$, a potential $V(\sigma)$, and a
Galileon term $G_{3}(\sigma,(\partial\sigma)^{2})\square\sigma$. In addition
to the modified dynamics, the Galileon term provides a screening mechanism to
potentially reconcile the models with general relativity predictions inside a
Vainshtein radius. Thanks to the Galileon term, the stability conditions,
namely ghost and Laplacian instabilities, in the branch with a negative
kinetic term ($Z=-1$) are fulfilled for a large volume of the parameter space.
Solving numerically the background evolution and linear perturbations, we
derive the constraints on the cosmological parameters in the presence of a
Galileon term for different combination of the cosmic microwave background
data from Planck, baryon acoustic oscillations measurements from BOSS, and
supernovae from the Pantheon compilation. We find that the Galileon term
alters the dynamics of all the studied cases. For a standard kinetic term
($Z=1$), we find that Planck data and a compilation of baryon acoustic
oscillations data constrain the Galileon term to small values that allow
screening very inefficiently. For a negative kinetic term ($Z=-1$), a Galileon
term and a non-zero potential lead to an efficient screening in a physically
viable regime of the theory, with a value for the Hubble constant today which
alleviates the tension between its cosmic microwave background and local
determinations. For a vanishing potential, the case with $Z=-1$ and the
Galileon term driving the late acceleration of the Universe is ruled out by
Planck data.

##  I Introduction

The Lambda cold dark matter $\Lambda$CDM model is the standard cosmological
model, it provides a remarkably good fit to most cosmological observations
such as the measurements of cosmic microwave background (CMB) anisotropies in
temperature and polarization [1], baryon acoustic oscillations (BAO) in the
galaxy and cluster distributions [2, 3], cosmic shear measurements of the
galaxy distribution [4, 5] and of the CMB [6, 7, 8, 9], and measurements of
the luminosity distances of type Ia supernovae (SN Ia) [10, 11, 12].
Nevertheless, in recent years its enduring success has been threatened by the
growing tension between the inference of $H_{0}$ within $\Lambda$CDM from the
Planck measurement of CMB anisotropy [13] and its determination with Supenovae
at low redshift, which has grown to a maximum level of $5\sigma$ with the
latest SH0ES release with SN Ia calibrated using Cepheids [14].

Whereas inferences of $H_{0}$ from different CMB experiments are very
consistent with one another [13, 15, 16], different approaches in calibrating
SN provide different estimates and different uncertainties in the
determination of $H_{0}$, as the latest SH0ES result based on Cepheids is
$H_{0}=73.04\pm 1.04\,\rm km\,s^{-1}\,Mpc^{-1}$ at 68% confidence level (CL)
[14], while the latest Carnegie-Chicago Hubble program estimate, using the tip
of the Red Giant Branch calibration technique is $H_{0}=69.8\pm
0.8\,\text{(stat)}\pm 1.7\,\text{(sys)}\,\rm km\,s^{-1}\,Mpc^{-1}$ at 68% CL
[17]. These differences in the local determinations of $H_{0}$ and their error
obviously affects the tension with the CMB inference in a statistically
significant way and might point to unknown systematics or uncertainties not
properly taken into account. Nevertheless, this tension that emerged with the
first Planck cosmological release [18] has fueled interest in models
explaining a larger value of $H_{0}$ compared to $\Lambda$CDM [19, 20, 21, 22,
23, 24, 25].

Modified gravity (MG) is one of the most promising fundamental physics route
to reconcile the inference of $H_{0}$ from the Planck measurement of CMB
anisotropy and its inferred value determined with the calibration of
supernovae at low redshift within $\Lambda$CDM model. A simple mechanism to
increase $H_{0}$ by an evolving gravitational constant is at play even in the
simplest scalar-tensor models, i.e. induced gravity (IG) (equivalent to the
Brans-Dicke model by a field redefinition). Without the fine-tuning of
parameters, the scalar field regulating the gravitational strength moves
around recombination because of its coupling to pressureless matter. This
mechanism, which can be also seen as a dynamical and self-consistent variation
in the gravitational constant, induces a degeneracy between $H_{0}$ and the
nonminimal coupling of the scalar field to the curvature leading to a larger
$H_{0}$ compared to $\Lambda$CDM [26, 27, 28, 29, 30, 31, 32, 33, 34].

Although this mechanism is completely general, the increase of $H_{0}$ with
respect to $\Lambda$CDM depends on the model and also on the number of
additional parameters with respect to $\Lambda$CDM. With Planck DR3
likelihoods, we have $\Delta H_{0}=H_{0}-H_{0}^{\Lambda{\mathrm{C}DM}}\sim
2\rm\,km\,s^{-1}\,Mpc^{-1}$ for IG [26, 27, 32] and $\Delta H_{0}\sim
4\rm\,km\,s^{-1}\,Mpc^{-1}$ for a nonminimally coupled scalar field with a
quartic potential [33], dubbed early modified gravity (EMG), which is among
the best solutions to the $H_{0}$ tension according to [23].

All the scalar-tensor models of gravity mentioned above do not have Vainshtein
screening and a nearly negligible chameleon screening due to the coupling to
matter. For this reason, one has to compare the volume of the parameter space
of the nonminimal coupling required to increase $H_{0}$ to the Solar System
bounds on the corresponding post-Newtonian parameters [35]. As for the $H_{0}$
degeneracy, also in this case, the comparison is nontrivial and model
dependent. Whereas for the simplest IG models, the Solar System bounds
constrain the nonminimal coupling to values which do not lead to an
appreciable increase of $H_{0}$ compared to $\Lambda$CDM [26, 27, 32], it is
sufficent to consider a nonminimally coupled coupled scalar field to have a
large volume of parameter space where viable models on Solar System scales can
lead to a value of $H_{0}$ larger than in $\Lambda$CDM [28, 31]. Early
modified gravity is a further step in the construction of a scalar-tensor
viable model on Solar System scales with large $H_{0}$ by adding a self-
interaction term which implements an attractor mechanism towards general
relativity (GR) with a cosmological constant.

In this paper, we instead consider the cosmology of IG, with the addition of a
Galileon term [36, 37, 38], whose purpose is to reconcile the theory with GR
on scales below the so-called Vainshtein radius. We therefore would like to
investigate the implications of the addition of an explicit term implementing
a screening mechanism on the cosmology of the simplest scalar-tensor theories
of gravity (see also [39, 40, 41] for previous studies with cosmological
bounds on Galileon gravity).

The paper is organized as follows. In Sec. II we describe the models
considered and their background dynamics together with the stability
conditions and the screening mechanism. In Sec. III we present evolution of
the linear fluctuations, CMB anisotropies and matter power spectrum. Section
IV is devoted to the presentation of the datasets and the details of our
Markov chain Monte-Carlo (MCMC) analysis, joint with the presentation of the
results. We draw our conclusions in Sec. V.

##  II The model

We study the following subclass of Horndeski theories [42, 37, 43, 44, 45]:

${\cal
S}=\int{\mathrm{d}}^{4}x\sqrt{|g|}\left[G_{4}(\sigma)R+G_{2}(\sigma,X)+G_{3}(\sigma,X)\square\sigma+\mathcal{L_{\rm
M}}\right]$ (1)

where $|g|$ is the absolute value of the determinant of the metric
$g_{\mu\nu}$ for which we use a $(-+++)$ signature, $\sigma$ is a scalar
field,
$X\equiv-\nabla_{\mu}\sigma\nabla^{\mu}\sigma/2=-\partial_{\mu}\sigma\partial^{\mu}\sigma/2$,
and $\Box~{}\equiv~{}\nabla_{\mu}\nabla^{\mu}$ is the covariant d’Alambert
operator. $\mathcal{L_{\rm M}}$ is the matter Lagrangian which does not depend
on the scalar field $\sigma$ and it is minimally coupled with the metric $g$.
The action (1) predicts that gravitational waves travel at the speed of light
[45], consistently with the current measurements [46]. We consider the
following form for the $G$ functions:

$\displaystyle\begin{split}G_{4}&=\gamma\sigma^{2}/2\\\ G_{3}&=-2g(\sigma)X\\\
G_{2}&=ZX-V(\sigma)+4\zeta(\sigma)X^{2}.\\\ \end{split}$ (2)

where $Z=\pm 1$ is the sign of the kinetic term. This Lagrangian is equivalent
to the extension of the Brans-Dicke (BD) model [47, 48] with a Galileon term.
In fact with the field redefinition $\phi=\gamma\sigma^{2}/2$, with
$\gamma=Z/(4\omega_{\rm BD})>0\text{ and }Z=\pm 1$, the $G$ functions become

$\displaystyle\begin{split}G_{4}&=\phi\\\ G_{3}&=-2f(\phi)\chi\\\
G_{2}&=2\frac{\omega_{\rm BD}}{\phi}\chi-V(\phi),\\\ \end{split}$ (3)

where $\chi\equiv-\nabla_{\mu}\phi\,\nabla^{\mu}\phi/2$ and the relationship
between $g,\zeta\text{ and}f$ is $g(\sigma)=\gamma\sigma^{3}f(\sigma)$;
$\zeta(\sigma)=\sigma^{-1}g(\sigma)$. Therefore, we refer to the model defined
by Eq. 2 as Brans-Dicke Galileon (BDG), while the model with
$G_{2}=ZX-V(\sigma)$ is the induced gravity Galileon (IGG), as it is the
extension of induced gravity (IG) [49] with a Galileon term
$G_{3}=-2g(\sigma)X$. Therefore even if BD and IG are equivalent, their
extensions with Galileon terms are not. In particular, IGG corresponds to BDG
given by Eq. 2 with, formally, $\zeta=0$, but it is not simply a special case
of BDG when $\zeta=0$. In fact, as described above, in the BDG model the
functions $g$ and $\zeta$ are not independent and setting $\zeta=0$ would mean
setting $g=0$ as well.

We study two different versions of BDG and IGG, relating to the sign of the
kinetic term $Z$:

Standard branch: $Z=+1$ in which the kinetic term has a standard sign. In this
branch we consider both IGG and BDG (IGGst and BDGst, respectively) and we
take $V(\sigma)\text{ and }g(\sigma)$ to be power laws of the scalar field:
$V(\sigma)=\lambda_{n}\sigma^{n}$, $g(\sigma)=\alpha_{m}\sigma^{m}$, for
several combinations of $n\text{ and }m$. In this scenario the potential
dominates the dynamics and provides the acceleration of the expansion of the
Universe, while the $G_{3}$ enters as a small correction to the standard BD
theory. Phantom branch: $Z=-1$ and $\gamma<1/6$, equivalent to $\omega_{\rm
BD}<-3/2$, in which the kinetic term has a nonstandard sign. In this branch we
consider only BDG with $g(\sigma)=\alpha\sigma^{-1}$, either with or without a
potential (cosmological constant). If there is no potential the burden to
provide comic acceleration, a healthy theory and effective screening on small
scales, is on the $G_{3}$ term [50]. Reinserting the potential into the
theory, both the Galileon term and the cosmological constant contribute to the
dynamics and give rise to the late-time acceleration of the Universe. Even if
we have a nonzero potential, the $G_{3}$ is still able to provide a healthy
theory and effective screening on small scales, where we recover (GR).

###  II.1 Covariant equations

By varying the action with respect to the metric we obtain the modified
Einstein field equation

$\displaystyle\begin{split}G_{\mu\nu}=\frac{1}{F(\sigma)}\bigg{[}&T_{\mu\nu}^{\rm(M)}+T_{\mu\nu}^{\rm(G)}+Z\big{(}\partial_{\mu}\sigma\partial_{\nu}\sigma-\frac{1}{2}g_{\mu\nu}\partial^{\rho}\sigma\partial_{\rho}\sigma\big{)}\\\
&-g_{\mu\nu}V(\sigma)+(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\Box)F(\sigma)\bigg{]}\,,\end{split}$
(4)

where $F(\sigma)=\gamma\sigma^{2}$, $G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}R$ is
the Einstein tensor and
$T_{\mu\nu}^{\rm(M)}~{}=~{}-\frac{2}{\sqrt{|g}}\frac{\delta(\sqrt{|g|}\mathcal{L}_{\rm
M})}{\delta g^{\mu\nu}}$ is the energy-momentum tensor of matter.
$T_{\mu\nu}^{(G)}$ is defined as

$\displaystyle\begin{split}T_{\mu\nu}^{\rm(G)}=-2\Big{\\{}\,&g(\sigma)\,\nabla_{\mu}\sigma\nabla_{\nu}\sigma\,\Box\sigma-\nabla_{(\mu}\,\sigma\,\nabla_{\nu)}\big{[}\,g(\sigma)(\partial\sigma)^{2}\,\big{]}\\\
&+\frac{1}{2}g_{\mu\nu}\nabla_{\alpha}\,\sigma\nabla^{\alpha}\big{[}\,g(\sigma)(\partial\sigma)^{2}\,\big{]}-\frac{\zeta(\sigma)}{2}g_{\mu\nu}(\partial\sigma)^{4}\\\
&+2\zeta(\sigma)\nabla_{\mu}\sigma\nabla_{\nu}\sigma\
(\partial\sigma)^{2}\,\Big{\\}},\end{split}$ (5)

with
$\nabla_{(\mu}\,\sigma\,\nabla_{\nu)}=\frac{1}{2}\big{(}\nabla_{\mu}\,\sigma\,\nabla_{\nu}+\nabla_{\nu}\,\sigma\,\nabla_{\mu}\big{)}$.
It is useful to write down the trace of the Einstein equation (4), which is

$R=\frac{1}{F}\Big{[}-T^{(\rm
M)}-T^{\rm(G)}+Z(\partial\sigma)^{2}+4V+3\,\square F(\sigma)\Big{]},$ (6)

with $T^{(G)}$ as

$T^{\rm(G)}=-2\bigg{\\{}g(\sigma)\,(\partial\sigma)^{2}\Box\sigma+\nabla_{\mu}\sigma\nabla^{\mu}\big{[}g(\sigma)\,(\partial\sigma)^{2}\big{]}\bigg{\\}}.$
(7)

Note that $T^{(G)}$ does not depend on $\zeta(\sigma)$.

The equation for the evolution of the scalar field, obtained by varying the
action (1) with respect to the field, is

$\displaystyle\begin{split}&\Box\sigma\big{[}Z-4\zeta\,(\partial\sigma)^{2}\big{]}-2g\,\big{\\{}(\Box\sigma)^{2}-\nabla^{\mu}\nabla^{\nu}\sigma\nabla_{\mu}\nabla_{\nu}\sigma\\\
&-\nabla^{\mu}\sigma\nabla^{\nu}\sigma
R_{\mu\nu}\big{\\}}+4g,_{\sigma}\nabla^{\mu}\sigma\nabla^{\nu}\sigma\nabla_{\mu}\nabla_{\nu}\sigma+g,_{\sigma\sigma}(\partial\sigma)^{4}\\\
&-3\zeta,_{\sigma}(\partial\sigma)^{4}-4\zeta(\sigma)\nabla_{\mu}[(\partial\sigma)^{2}]\nabla^{\mu}\sigma+\gamma\sigma
R-V,_{\sigma}=0.\end{split}$ (8)

###  II.2 Background evolution

Working in a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW)
universe described by the line element

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)\,\mathrm{d}x^{2}\,,$ (9)

the Einstein field equations reduce to:

$\displaystyle\begin{split}3FH^{2}=&\rho+\frac{1}{2}Z\dot{\sigma}^{2}-3H\dot{F}+V(\sigma)\\\
&+\dot{\sigma}^{3}\big{[}6g(\sigma)H-\dot{g}(\sigma)+3\zeta(\sigma)\dot{\sigma}\big{]}\\\
\equiv&\rho+\rho_{\sigma}\,,\end{split}$ (10)
$\displaystyle\begin{split}-2F\dot{H}=&\rho+p+Z\dot{\sigma}^{2}+\ddot{F}-H\dot{F}\\\
&+\dot{\sigma}^{2}\big{[}(6gH-2g,_{\sigma}\dot{\sigma})\dot{\sigma}+4\zeta\dot{\sigma}^{2}-2g\ddot{\sigma}\big{]}\\\
\equiv&\rho+p+\rho_{\sigma}+p_{\sigma}\,,\end{split}$ (11)

where:

$\rho_{\sigma}=\frac{Z}{2}\dot{\sigma}^{2}-3H\dot{F}+V(\sigma)+\dot{\sigma}^{3}\big{[}6g(\sigma)H-\dot{g}(\sigma)+3\zeta(\sigma)\dot{\sigma}\big{]},$
(12)
$p_{\sigma}=\frac{Z}{2}\dot{\sigma}^{2}-V(\sigma)+\ddot{F}+2H\dot{F}-\dot{\sigma}^{4}(g,_{\sigma}-\zeta)-2g\dot{\sigma}^{2}\ddot{\sigma}.$
(13)

The scalar field equation in the FLRW metric takes the form:

$\displaystyle\begin{split}&\ddot{\sigma}\big{(}Z+12gH\dot{\sigma}-4(g,_{\sigma}-3\zeta)\dot{\sigma}^{2}\big{)}-6\gamma\sigma(2H^{2}+\dot{H})+V,_{\sigma}\\\
&+3ZH\dot{\sigma}+6g(3H^{2}+\dot{H})\dot{\sigma}^{2}+12H\zeta\dot{\sigma}^{3}\\\
&-(g,_{\sigma\sigma}-3\zeta,_{\sigma})\dot{\sigma}^{4}=0.\end{split}$ (14)

Because of the nonminimal coupling between the scalar field and the Ricci
scalar in the Lagrangian, the Newton constant in the Friedmann equations is
replaced by $G_{\rm cosm}=(8\pi F)^{-1}$ that now varies with time.

Following the notation of [51] we define the density parameters for radiation
(r), pressureless matter (m) and the scalar field ($\sigma$) as

$\widetilde{\Omega}_{i}=\frac{\rho_{i}}{3FH^{2}}\equiv\frac{\rho_{i}}{\rho_{\rm
crit}}\quad(i={\rm r,m},\sigma).$ (15)

It is also convenient to define dark energy density and pressure parameters in
a framework which mimics Einstein gravity at the present time, by rewriting
the Friedmann equations as [52, 53, 51]

$\displaystyle\begin{split}3F_{0}H^{2}&=\rho+\rho_{\rm DE}\,,\\\
-2F_{0}\dot{H}&=\rho+p+\rho_{\rm DE}+p_{\rm DE},\end{split}$ (16)

which leads to

$\rho_{\rm
DE}=\frac{F_{0}}{F}\rho_{\sigma}+\rho\bigg{(}\frac{F_{0}}{F}-1\bigg{)}\,,$
(17) $p_{\rm
DE}=\frac{F_{0}}{F}p_{\sigma}+p\bigg{(}\frac{F_{0}}{F}-1\bigg{)}\,.$ (18)

In this way, the effective parameter of state for dark energy can be defined
as $w_{\rm DE}\equiv p_{\rm DE}/\rho_{\rm DE}$. In this context the density
parameters mimicking radiation, matter and dark energy (DE) in Einstein
gravity are

$\Omega_{i}=\frac{\rho_{i}}{3F_{0}H^{2}}\quad(i={\rm r,m,DE}).$ (19)

The two definitions in Eqs. 15 and 19 coincide at $z=0$:
$\widetilde{\Omega}_{0,i}=\Omega_{0,i}$.

####  II.2.1 Standard branch $(Z=1)$

In the standard branch we consider IGG: the theory without $\zeta(\sigma)$:
$G_{4}=\gamma\sigma^{2}/2,\,G_{2}=X-V(\sigma),\,G_{3}=-2g(\sigma)X$. In the
region in parameter space that produces a reasonable cosmological background
evolution BDG and IGG are nearly indistinguishable. For this reason we
restrict ourselves to IGG for $Z=1$.

We consider a monomial potential $V(\sigma)=\lambda_{n}\sigma^{n}$ and
$g(\sigma)=\alpha_{m}\sigma^{m}$. For this choice of $V(\sigma)\text{ and
}g(\sigma)$ the IGG model presents exact solutions with accelerated expansion
in the absence of matter where the scale factor is $a(t)\propto t^{p}$, with
$n$ and $m$ related by $m=1-n$, and

$\displaystyle
p=\frac{2(-2+n+4\beta-4\gamma+n^{2}\gamma)}{(24\beta-16\gamma+20n\gamma-8n^{2}\gamma+n^{3}\gamma)}\,,$
(20) $\displaystyle\sigma(t)=c_{0}t^{-2/(n-2)},$ (21)

where $\beta$, defined by $\alpha\equiv\beta c_{0}^{n-2}$, is a
reparametrization of $\alpha$, useful to show that in the limit of
$\alpha\rightarrow 0$ ($\beta\rightarrow 0$) we recover the analogous solution
found in IG [27].

For $n=4$ and $n=2$, there are de Sitter solutions $a(t)\propto{e^{Ht}}$ with
a constant $H>0$. For $n=4$, the solution found for IG [27]: $\sigma=\pm
H\sqrt{3\gamma/\lambda_{4}}$ is still valid in IGG for every possible form of
$g(\sigma)$, as this term does not contribute under the ansatz
$\sigma=\text{constant}$. For $n=2$, when $g(\sigma)\propto\sigma^{-1}$, there
are solutions of the form $\sigma\propto e^{Ht\delta}$, with
$\delta\equiv\delta(\alpha,\gamma,H)$. Its explicit form is

$\delta=-1,\,\text{ or
}\delta=\frac{-1-4\gamma\pm\sqrt{(1+4\gamma)^{2}+48H^{2}\gamma\alpha}}{12H^{2}\alpha}.$
(22)

For $\alpha\rightarrow 0$ these solutions reduce to the ones discussed in
[54]:

$\delta=-1\,\text{ or }\delta=\frac{2\gamma}{1+4\gamma}\,.$ (23)

The evolution of the scalar field is shown in the upper panel of Fig. 1, in
the model with a constant potential $V(\sigma)=\lambda_{0}\equiv\Lambda$ and
$g(\sigma)=\alpha$, compared with the analogous IG model with $\gamma=5\times
10^{-5}$. We consider the rescaled, adimensional quantity
$\tilde{\alpha}=\alpha\times{\rm(Mpc\,[GeV]^{-1})^{-2}\times(M_{\rm
Pl}\,[GeV])^{1+m}}$ and it can be seen that the departure from IG is
significant for larger $\alpha$, i.e. stronger gravity at early times. Deep in
the radiation era the field is nearly at rest but then it grows steeply (the
larger $\alpha$ the steeper the growth) reaching the value expected in IG in
the matter era and evolving until today in the same way as it does in IG.

Figure 1: Time evolution of the scalar field (top panel), the equation of
state parameter for dark energy (middle panel) and the density parameters
(bottom panel) in IGG standard ($Z=1$) with $\gamma=5\times 10^{-5}$, constant
$V(\sigma)=\Lambda$ and $g(\sigma)=\alpha$, for different values of
$\tilde{\alpha}$.

The value of the field at $z=0$ is fixed by requiring that the effective
gravitational constant [55, 56]

$G_{\rm eff}=\frac{1}{16\pi
G_{4}}\Bigg{[}\frac{4G_{4,\sigma}^{2}+G_{4}(G_{2,X}-2G_{3,\sigma})}{3G_{4,\sigma}^{2}+G_{4}(G_{2,X}-2G_{3,\sigma})}\Bigg{]}$
(24)

coincide at $z=0$ with the measured value of the Newton’s constant $G$. For
$m=0$, the $G_{3}$ does not contribute to Eq. 24, as it enters only through
the first derivative with respect to the field, while, for $m\neq 0$ one can
check that the contribution is negligible due to the redshift evolution of the
scalar field in the matter era and we can therefore use the IG approximation
[57]

$8\pi G_{\rm
eff}(z=0)=\frac{1}{\gamma\sigma_{0}^{2}}\frac{1+8\gamma}{1+6\gamma}$ (25)

to fix the value of $\sigma(z=0)=\sigma_{0}$.

In the bottom plot of Fig. 1 we show the evolution of density parameters in
IGG with a constant potential $V(\sigma)=\Lambda$ and $g(\sigma)=\alpha$. The
values of $\alpha$ are chosen large enough to show the effect of the $G_{3}$
term in this scenario, which causes a different evolution of
$\widetilde{\Omega}_{\rm r}$ and $\widetilde{\Omega}_{\sigma}$ in the early
Universe with respect to IG and $\Lambda\rm CDM$.

The middle panel of Fig. 1 presents the evolution of the parameter of state of
dark energy $w_{\rm DE}$ defined above, it tracks the IG behavior at late
times but departs from it at $z\geq 10^{3}$, in correspondence to the
analogous uptick in the evolution of $\widetilde{\Omega}_{\sigma}$. In fact,
$w_{\rm DE}$ follows the dominant component: deep in the radiation epoch it
has a value close to $1/3$ , then in the matter era it decreases towards zero;
finally, at present epoch, it becomes negative, $w_{\rm DE}\simeq-1$,
mimicking a cosmological constant. The bump which approximately occurs in the
radiation era corresponds to the epoch in which the energy density of the
field grows and becomes of the same order of that of radiation. The growth of
$\Omega_{\sigma}$ so early in time is due to inefficient cosmological
Vainshtein screening [58] in this model for the range of parameters
considered.

Figure 2: Time evolution of the scalar field (top panel), the equation of
state parameter for dark energy (middle panel) and the density parameters
(bottom panel) in IGG standard ($Z=1$) with $\gamma=5\times 10^{-5}$, constant
$V(\sigma)=\Lambda$ and $g(\sigma)~{}=~{}\alpha\sigma^{m}$, for different
values of the exponent $m$ and fixed $\tilde{\alpha}~{}=~{}5\times 10^{-6}$.

In Fig. 2 we show the dependence of the background quantities
($\sigma,\,w_{\rm DE},\,\widetilde{\Omega}_{i}$) on the exponent $m$ in the
function $g(\sigma)=\alpha\sigma^{m}$, while keeping $\gamma=5\times
10^{-5}\text{ and }\tilde{\alpha}=5\times 10^{-6}$ fixed, together with the
exponent $n=0$ in $V(\sigma)\propto\sigma^{n}$. We point out that changing $m$
while keeping $\tilde{\alpha}$ fixed modifies the actual value of $\alpha$ and
its dimensions since
$\tilde{\alpha}=\alpha\times{\rm(Mpc\,[GeV]^{-1})^{-2}\times(M_{\rm
Pl}\,[GeV])^{1+m}}$. We can see from the figures that for small values of $m$
the time evolution of the scalar field is suppressed and with that the
equation of state parameter and the density parameters resemble more those of
$\Lambda\rm CDM$. For larger values of the exponent $m$ the scalar field tends
to dominate the dynamics both in the early Universe, reaching equality with
radiation at $z=10^{6}$ for $m=1$, and in the late Universe where it mimics a
cosmological constant. This early epoch of scalar field domination also delays
the matter radiation equality that happens at $z<10^{4}$, with
$\widetilde{\Omega}_{\rm r}=\widetilde{\Omega}_{\rm m}\simeq 0.1$, and leads
to two different epochs in which $\widetilde{\Omega}_{\rm m}$ is equal to the
density parameter of the scalar field. One of these epochs is the beginning of
the matter era, around $z\simeq 10^{-3}$ and one is the onset of the late-time
acceleration of the expansion of the Universe in which the field dominates the
dynamics once again.

####  II.2.2 BDG phantom

In the phantom branch we study BDG with $g(\sigma)=\alpha\sigma^{-1}$ and
$\zeta(\sigma)=\sigma^{-1}g(\sigma)$, because, for this choice, it was shown
in [50] that, for a null potential and in absence of matter and radiation,
$\rho=p=0$, there is a self-accelerating solution with $\dot{H}=\dot{Q}=0$,
with $Q=\dot{\sigma}/\sigma$. This solution satisfies

$y\equiv\frac{Q}{H}=\gamma\frac{-4\pm\sqrt{-32-6Z/\gamma}}{Z+8\gamma};$ (26)

and it is real for $Z=-1\text{ and }\gamma<3/16$ (equivalent to $\omega_{\rm
BD}<-4/3$), for which Eqs. 10 and 14 give

$H^{2}=\frac{\gamma}{\alpha}\,\frac{3+6y+y^{2}/(2\gamma)}{2y^{3}(3+2y)}.$ (27)

For a full study of the theory with $\Lambda=0$ which exactly respects Eq. 27
see appendix A.

Inspired by these analytical solutions, we study an intermediate case always
considering $g(\sigma)=\alpha\sigma^{-1}$ and reinserting the potential. In
doing so we obtain a theory that has a $\Lambda\rm CDM$ limit. In order to
provide the late-time cosmic acceleration the parameters $\alpha$ and
$\Lambda$ should be fine-tuned, the procedure is the following: we pick a
value for $\alpha$ and use a shooting algorithm on the cosmological constant
to obtain, in a flat Universe, $\Omega_{0,{\rm
DE}}(\Lambda,\alpha)=1-\Omega_{0,\rm m}-\Omega_{0,\rm r}$. In this way the
burden to provide the late-time acceleration of the expansion of the Universe
is shared between the potential and the Galileon term. Today’s value of the
scalar field is such that the Planck mass is

$M_{\rm Pl}^{2}(z=0)\equiv\gamma\sigma_{0}^{2}=1.$ (28)

We do not set $\sigma_{0}$ following Eq. (25) or (24) because, due to the
Vainshtein screening mechanism, the theory reduces to GR at small scales, in
the sense that the post-Newtonian parameters are those of GR and the
gravitational constant on small scales is $G_{\rm N}=G_{\rm cosm}$ [59]. In
this way, thanks to the screening mechanism we recover GR with the correct
value of the gravitational constant on small scales today (see Sec. II.4 for
details). In the following, instead of working directly with the parameter
$\alpha$, we use $\tilde{\alpha}$ defined previously, rescaled by a factor
$10^{8}$: $\tilde{\alpha}_{8}\equiv 10^{-8}\tilde{\alpha}$. All the plots and
the MCMC constraints will be expressed as a function of
$1/\tilde{\alpha}_{8}$. The reason for using the inverse of $\tilde{\alpha}$
is that the $\Lambda\rm CDM$ limit is obtained when
$\tilde{\alpha}\rightarrow\infty\,(1/\tilde{\alpha}\rightarrow 0)$.

Figure 3: Time evolution of the scalar field (top panel), the equation of
state parameter for dark energy (middle panel) and the density parameters
(bottom panel) in BDG phantom ($Z=-1$) with $\gamma=5\times 10^{-5}$, constant
$V(\sigma)=\Lambda$ and $g(\sigma)=\alpha\sigma^{-1}$, for different values of
$1/\tilde{\alpha}_{8}$.

In the top panel of Fig. 3 we show the evolution of the scalar field with
redshift as a function of $1/\tilde{\alpha}_{8}$: for larger values of this
parameter the field starts at a lower initial value and it grows more steeply
in the late Universe to reach its value at $z=0$ fixed by the requirement
$\gamma\sigma_{0}^{2}=1$. Contrary to the standard branch where the field
started to evolve deep in the radiation era, here the field is frozen
approximately until the epoch in which the $\sigma$ energy density starts to
grow. This is a manifestation of the cosmological Vainshtein screening
mechanism [58]. The fact that the field is frozen means that for a large
portion of the history of the universe we have a gravitational constant which
is effectively constant but different from the value measured here on Earth
today. In the bottom panel of Fig. 3 we show the evolution of the density
parameters and it can be seen that already for $1/\tilde{\alpha}_{8}=0.1$ the
background expansion in this model is indistinguishable from $\Lambda\rm CDM$.
Only for more extreme values we see departures from $\Lambda\rm CDM$ as the
growth of the dark energy density is delayed and then steeper for
$1/\tilde{\alpha}_{8}=0.4$. This behavior is confirmed also by the time
evolution of $w_{\rm DE}$ that approaches and reaches the value $-1$ earlier
for smaller values of $1/\tilde{\alpha}_{8}$; while for the larger values we
observe a phantom behavior with $w_{\rm DE}<-1$ and $w_{\rm DE}\neq-1$ today,
reaching $-1$ eventually in the future. In this scenario, therefore the
$\Lambda\rm CDM$ limit is approached when $1/\tilde{\alpha}_{8}\rightarrow 0$.
This is the reason why we are considering the inverse of $\alpha$ as a
parameter: in the Markov chain Monte-Carlo analysis we will sample on
$1/\tilde{\alpha}_{8}$ in order to have the $\Lambda\rm CDM$ limit at zero and
not at infinity, as we would if we had sampled on $\alpha$.

###  II.3 Stability conditions in the phantom branch

When the sign of the kinetic term in the Lagrangian is nonstandard, i.e.,
$Z=-1$, and $G_{3}=0$, scalar-tensor theories can suffer from ghost and
Laplacian instabilities. Ghost instabilities occur when the sign of the
kinetic term is negative in the second-order action for perturbations of the
Horndeski action, whereas Laplacian instabilities occur when the speed of
sound squared becomes negative [60]. In general, the theory in Eq. 1 is free
from ghost and Laplacian instabilities if [60, 61]

$\displaystyle\begin{split}q_{s}\equiv\,&4G_{4}\big{\\{}G_{2X}+2G_{3\sigma}+\dot{\sigma}\big{[}(G2_{XX}+G3_{X\sigma})\dot{\sigma}\\\
&-6G_{3X}H\big{]}\big{\\}}+3(2G_{4\sigma}+G_{3X}\dot{\sigma}^{2})^{2}>0,\end{split}$
(29)
$\displaystyle\begin{split}c_{s}^{2}\equiv\,\big{[}&4G_{2X}G_{4}+8G_{3\sigma}G_{4}\\\
&+(6G_{4\sigma}^{2}-G_{3X}\dot{\sigma}^{2})(2G_{4\sigma}^{2}+G_{3X}\dot{\sigma}^{2})\\\
&-8G_{4}(G_{3X}\ddot{\sigma}+2G_{3X}H\dot{\sigma}+G_{3X\sigma}\dot{\sigma}^{2})\big{]}/q_{s}>0,\end{split}$
(30)

where, $G_{i\,\sigma},G_{i\,X}$ indicate the derivative of the $i$th $G$
function with respect to the scalar field and the kinetic term, respectively.
From Eqs. 29 and 30, it can be seen that even IG with $Z=-1$ and
$0<\gamma<1/6$ ($\omega_{\rm BD}<-3/2$) - a region of the parameter space
which would contain otherwise a ghost - can be stable thanks to the addition
of the $G_{3}$ term in the Lagrangian. This is fully confirmed by Fig. 4 where
we show Eqs. 29 and 30 as functions of redshift for the BDG in the phantom
branch. Both conditions are satisfied at all times for the values of
$1/\tilde{\alpha}_{8}$ we have considered.

Note, however, that in BDGph the speed of sound of scalar perturbations can
become temporarily superluminal for large values of $1/\tilde{\alpha}_{8}$ at
redshifts close to the transition to dark energy domination, as can be seen in
the bottom panel of Fig. 4. This was first noticed in [50] for the case
$\Lambda=0$ and it could potentially put a constraint on the value of
$1/\tilde{\alpha}_{8}$. If a theory with superluminal propagation is viable or
not is still a controversial issue: some authors claim it is not problematic
[62, 63], others argue the opposite [64, 65, 66]. We choose conservatively to
not impose any theoretical prior $c_{s}^{2}\leq 1$ in our MCMC analysis. Note
that the stability conditions in the tensor sector are automatically satisfied
by the choice of our Lagrangian and parameters: the regime $\gamma>0$ we
consider ensures that the cosmological Planck mass and the square of the
tensor speed of sound are positive.

Figure 4: Stability conditions: ghost instability (top panel), Laplacian
instability (bottom panel) for BDG in the phantom branch ($Z=-1$) with
$\gamma=5\times 10^{-5}$, constant $V(\sigma)=\Lambda$ and
$g(\sigma)=\alpha\sigma^{-1}$, for different values of $1/\tilde{\alpha}_{8}$.

###  II.4 Screening

The scalar degree of freedom in scalar-tensor theories of gravity can leave
imprints also on scales which have detached from the cosmological expansion
where deviations from GR are strongly constrained. Thus, a theory of MG should
have either small deviations from GR on Solar System scales or be equipped
with a screening mechanism. The Vainshtein screening is a mechanism that
operates in theories with a self-interaction of the form $X\Box\sigma$, like
the models considered in this paper: around local sources the nonlinear
interactions lead to the decoupling of the scalar field from the remaining
degrees of freedom within the so-called Vainshtein radius $r_{V}$. For scales
$r<r_{V}$ the theory of gravity reduces to GR: the post-Newtonian parameters
are those of GR and the gravitational constant on small scales is $G_{\rm
N}=G_{\rm cosm}$ [59].

As it was demonstrated in [59], the Vainshtein mechanism cannot suppress the
cosmological evolution of the scalar field: the mechanism is not able to
determine the coupling constant of the theory (Newton gravitational constant
or $\sigma_{0}$). The cosmological evolution of the scalar field sets the
value of the coupling constant of GR today and on small scales. For this
reason, for models with effective Vainshtein screening, we fix $\sigma_{0}$
according to Eq. 28 in our background evolution.

The Vainshtein radius, which is the radius within which the theory reduces to
GR, in a cosmological background, for a spherical astronomical object of mass
$\delta M$, is given by [59]

$r_{V}=\left(\mathcal{BC}\mu/H^{2}\right)^{1/3},$ (31)

where $\mu=\delta M/(16\pi G_{4})$, and $\mathcal{B}\text{ and }\mathcal{C}$
are functions of the Hordenski functions $G_{i}$s:

${\cal
B}\equiv\frac{4\beta_{0}}{\alpha_{0}+2\alpha_{1}\alpha_{2}+\alpha_{2}^{2}},\quad$
(32) ${\cal
C}\equiv\frac{\alpha_{1}+\alpha_{2}}{\alpha_{0}+2\alpha_{1}\alpha_{2}+\alpha_{2}^{2}},$
(33)

where

$\alpha_{i}(t)\equiv\frac{A_{i}}{{\cal
G}_{T}},\quad\beta_{0}(t)\equiv\frac{B_{0}}{{\cal G}_{T}},\quad$ (34)

with

$\beta_{0}\equiv\frac{\alpha_{1}}{2}+\alpha_{2}\;(\neq 0),$ (35)

and

$\displaystyle{\cal F}_{T}$ $\displaystyle\equiv$ $\displaystyle
2G_{4}\equiv{\cal G}_{T}$ (36) $\displaystyle{\cal E}$ $\displaystyle\equiv$
$\displaystyle 2XG_{2X}-G_{2}-6X\dot{\sigma}HG_{3X}+2XG_{3\sigma}$ (37)
$\displaystyle-6H^{2}G_{4}-6H\dot{\sigma}G_{4\sigma}$ $\displaystyle{\cal P}$
$\displaystyle\equiv$ $\displaystyle
G_{2}+2X\left(G_{3\sigma}+\ddot{\sigma}G_{3X}\right)+2\left(3H^{2}+2\dot{H}\right)G_{4}$
(38)
$\displaystyle+2\left(\ddot{\sigma}+2H\dot{\sigma}\right)G_{4\sigma}+4XG_{4\sigma\sigma}$
$\displaystyle\Theta$ $\displaystyle\coloneqq$
$\displaystyle\dot{\sigma}XG_{3X}+2HG_{4}+\dot{\sigma}G_{4\sigma},$ (39)
$\displaystyle A_{0}$ $\displaystyle\equiv$
$\displaystyle\frac{\dot{\Theta}}{H^{2}}+\frac{\Theta}{H}-{\cal
G}_{T}-2\frac{\dot{\cal G}_{T}}{H}-\frac{{\cal E}+{\cal P}}{2H^{2}},$ (40)
$\displaystyle A_{1}$ $\displaystyle\equiv$ $\displaystyle\frac{1}{H}\dot{\cal
G}_{T}+{\cal G}_{T}-{\cal F}_{T},$ (41) $\displaystyle A_{2}$
$\displaystyle\equiv$ $\displaystyle{\cal G}_{T}-\frac{\Theta}{H},$ (42)
$\displaystyle B_{0}$ $\displaystyle\equiv$
$\displaystyle-\frac{X}{H}\dot{\sigma}G_{3X}.$ (43)

Figure 5: Screening: Vainshtein radius for the sun in BDG phantom ($Z=-1$)
with $\gamma=5\times 10^{-5}$, constant $V(\sigma)=\Lambda$ and
$g(\sigma)=\alpha\sigma^{-1}$, for different values of $1/\tilde{\alpha}_{8}$.

In Fig. 5 we show the cosmological evolution of the Vainshtein radius for an
object with $\delta M=M_{\odot}$, for BDG in the phantom branch. The screening
is effective and the Vainshtein radius starts growing for $z\lesssim 10$,
reaching $r_{V\odot}\gtrsim 90\,{\rm pc}$ today, successfully recovering GR in
the Solar System. We do not present Vainshtein radius in the standard branch
since the screening is not effective and the Vainshtein radius is very small.
This is the reason why $\sigma_{0}$ is fixed using Eq. 25 in the background
evolution in the standard branch.

##  III Linear perturbations and CMB anisotropies

In the synchronous gauge the perturbed FLRW metric is to first order

$\mathrm{d}s^{2}=a^{2}(\tau)\left[-\mathrm{d}^{2}\tau+(\delta_{ij}+h_{ij})\mathrm{d}x^{i}\mathrm{d}x^{j}\right]$
(44)

where $\tau$ is the conformal time, ${\bf x}$ the spatial comoving coordinate,
and $h_{ij}$ is the metric perturbation.

From this point on, we move in Fourier space for the calculation of the
perturbed quantities. The scalar mode of $h_{ij}$, following the conventions
of Ref. [67], can be expressed as a Fourier integral:

$\displaystyle\begin{split}h_{ij}(\tau,{\bf
x})=\int\mathrm{d}^{3}ke^{\imath{\bf k\cdot
x}}\bigg{[}&\hat{k}_{i}\hat{k}_{j}\,h(\tau,{\bf k})\\\
&+\left(\hat{k}_{i}\hat{k}_{j}-\frac{\delta_{ij}}{3}\right)6\,\eta(\tau,{\bf
k})\bigg{]}\end{split}$ (45)

where $\hat{k}_{i}\equiv k_{i}/k$ with $k\equiv|{\bf k}|$ and
$h\equiv\delta^{ij}h_{ij}$ is the Fourier transform of trace of
$h_{ij}(\tau,\,{\bf x})$. The scalar field perturbation $\delta\sigma$ is, in
this notation,

$\delta\sigma(\tau,\,{\bf x})=\int\mathrm{d}^{3}ke^{\imath{\bf k\cdot
x}}\,\delta\sigma(\tau,\,{\bf k}).$ (46)

####  III.0.1 The perturbed Einstein and scalar field equations

Following [67] we use the definitions of the velocity potential $\theta$ and
the anisotropic stress perturbation $\Theta$

$\displaystyle(\bar{\rho}+\bar{p})\theta\equiv\imath k^{i}\delta T^{0}_{\
i}\,,$ (47)
$\displaystyle(\bar{\rho}+\bar{p})\Theta\equiv-\left(\hat{k}_{i}\hat{k}_{j}-\frac{\delta_{ij}}{3}\right)\Sigma^{i}_{\
j}\,,$ (48) $\displaystyle\Sigma^{i}_{\ j}\equiv T^{i}_{\
j}-\delta^{i}_{j}\frac{T^{k}_{\ k}}{3}\,,$ (49)

to write down the perturbed Einstein equations in the following form:

$\displaystyle
k^{2}\eta-\frac{1}{2}\mathcal{H}h^{\prime}=-\frac{a^{2}}{2}\big{(}\delta\tilde{\rho}+\delta\tilde{\rho}^{(G)}\big{)}\,,$
(50) $\displaystyle
k^{2}\eta^{\prime}=\frac{a^{2}}{2}\Big{[}\big{(}\tilde{\rho}+\tilde{P}\big{)}\tilde{\theta}+\big{(}\tilde{\rho}^{(G)}+\tilde{P}^{(G)}\big{)}\tilde{\theta}^{(G)}\Big{]}\,,$
(51) $\displaystyle
h^{\prime\prime}+2\mathcal{H}h^{\prime}-2k^{2}\eta=-3a^{2}\big{(}\delta\tilde{P}+\delta\tilde{P}^{(G)}\big{)}\,,$
(52) $\displaystyle
h^{\prime\prime}+6\eta^{\prime\prime}+2\mathcal{H}(h^{\prime}+6\eta^{\prime})-2k^{2}\eta=-3a^{2}\big{(}\tilde{\rho}+\tilde{P}\big{)}\tilde{\Theta}\,,$
(53)

where a tilde denotes the effective perturbations defined in appendix B. The
equation for the evolution of the scalar field fluctuation in Fourier space is

$\displaystyle\begin{split}\delta\sigma^{\prime\prime}=\frac{\delta\tilde{\cal{G}}}{Z+6\gamma+\frac{2\sigma^{\prime}}{a^{2}}\bigg{[}3\,g\,\Big{(}2\mathcal{H}-\frac{\sigma^{\prime}}{\sigma}\Big{)}+\sigma^{\prime}(6\zeta-2g_{,\sigma})\bigg{]}}\,,\end{split}$
(54)

where the explicit expression for $\delta\tilde{\cal G}$ can be found in
appendix B.

In Eqs. 50, 51, 52 and 53 we distinguished explicitly between the quantities
coming from IG and the ones arising from the Galileon term to make as manifest
as possible the reduction of IGG and BDG to the IG equations of [68, 69] when
$g(\sigma)=\zeta(\sigma)=0$.

####  III.0.2 Imprints on CMB anisotropies and large scale structure

In order to obtain the predictions for the CMB anisotropy angular power
spectra and the matter power spectrum we have used an extension of the code
CLASSig [26] which includes the models studied in our paper. The basis for
this code is the publicly available Einstein-Boltzmann solver
CLASS111https://lesgourg.github.io/class_public/class.html [70].

Figure 6: Relative differences with respect to $\Lambda\rm CDM$ of CMB TT
(left panels) and EE (right panels) power spectrum. Top panel: IGG standard
with $\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda,\;g(\sigma)=\alpha$; bottom
panel: BDG phantom with $\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda\text{ and
}g(\sigma)=\alpha\sigma^{-1}$. Figure 7: Relative differences with respect to
$\Lambda\rm CDM$ in CMB lensing potential power spectrum (left panels) and
matter power spectrum (right panels). Top panel: IGG standard with
$\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda,\;g(\sigma)=\alpha$; bottom panel:
BDG phantom with $\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda\text{ and
}g(\sigma)=\alpha\sigma^{-1}$. Figure 8: Relative differences with respect to
$\Lambda\rm CDM$ of the ISW contribution to the CMB TT angular power spectrum
for BDG phantom with $\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda\text{ and
}g(\sigma)=\alpha\sigma^{-1}$.

In Fig. 6 we show the deviations in the temperature (TT) and and E-mode
polarization (EE) CMB angular power spectra with respect to $\Lambda\rm CDM$,
whereas in Fig. 7 we show how the CMB lensing potential angular power spectrum
and the total linear matter power spectrum at redshift $z=0$ differ in these
theories with respect to the standard model. In the top panels of Figs. 6 and
7 we have IGG in the standard branch, while in the bottom panels we show BDG
phantom with potential $V(\sigma)=\Lambda$.

We note how in the standard branch there are small effects in the TT power
spectrum at at low multipoles, resulting in an ISW effect extremely similar to
the one observed in $\Lambda\rm CDM$. This behavior is common to all the
spectra in IGG at low multipoles, for which the differences with respect to GR
are small. The departure from $\Lambda\rm CDM$ is more evident on the scales
of the acoustic peaks in TT and on smaller scales, where, due to the
modification of gravity around the time of recombination, the contribution of
the scalar field shifts and enhances the peaks of the power spectrum.

On the contrary, in the model with $Z=-1$ we can observe large departures from
$\Lambda\rm CDM$ at low multipoles in TT, due to the interplay of $G_{3}$ and
the ISW effect [38, 40]. We show the relative difference in the ISW
contribution to the TT angular power spectrum in Fig. 8.

For the BDG in the phantom branch, while we see differences on all scales in
all power spectra of Figs. 6 and 7, due to the presence of the cosmological
constant, it is clear that as $1/\tilde{\alpha}_{8}$ gets smaller we recover
the predictions of $\Lambda\rm CDM$ and there is therefore room to provide a
good fit to the data even with $1/\tilde{\alpha}_{8}\neq 0$ but small.

The differences in the matter power spectrum with respect to the standard
model, presented in the right panels of Fig. 7, are interesting targets for
upcoming large scale structure data such as DESI222https://www.desi.lbl.gov,
Euclid333https://www.esa.int/Science_Exploration/Space_Science/Euclid, or Vera
Rubin observatory444https://www.vro.org/, https://www.lsst.org/, which will
significantly tighten constraints on these models.

##  IV Constraints from cosmological observations

In this section, we present the constraints on the cosmological parameters of
the models studied: IGG in the standard branch (IGGst) and BDG in the phantom
branch with a potential (BDGph). In the following results the potential is
always given by a $V(\sigma)=\Lambda$, while the function $g(\sigma)$ is a
constant in the standard branch and $g(\sigma)=\alpha\sigma^{-1}$ for $Z=-1$.
Note that in the following we consider only IGG in the standard branch since,
for the same range of parameters, BDG gives exactly the same results.

###  IV.1 Methodology and datasets

We study these models first with CMB data only and we extend the analysis to
other datasets when necessary in order to constrain a parameter within our
prior or when the model is not ruled out by CMB alone and the results are
interesting enough justify the study with additional likelihoods. We perform a
Markov chain Monte Carlo analysis using the publicly available sampling code
MontePython-v3555https://github.com/brinckmann/montepython_public [71, 72]
connected to the aforementioned extension of CLASSig [26]. For the sampling,
we use the Metropolis-Hastings algorithm with a Gelman-Rubin [73] convergence
criterion $R-1<0.01$. The reported mean values and uncertainties on the
parameters, together with the contour plots have been obtained using
GetDist666https://getdist.readthedocs.io/en/latest/ [74]. We make use of CMB
data in temperature, polarization, and lensing from Planck Data Release $3$
[75, 6]. We consider the Planck baseline likelihood (hereafter P18) which is
composed by the Plik likelihood on high multipoles, $\ell>30$, the commander
likelihood on the lower multipoles for temperature and SimAll for the E-mode
polarization [75], and the conservative multiple range, $8<\ell<400$, for the
CMB lensing. As complementary data to Planck, we use BAO data from to the
post-reconstruction measurements from BOSS DR12 [76]; low-$z$ BAO measurements
from SDSS DR7, 6dF, and MGS [77, 78]; Ly$\alpha$ BAO measurements from eBOSS
DR14, and a combination of those [79, 80, 81]. Hereafter we call this
combination of datasets simply BAO. For the BDGph model, on top of BAO, we
will also consider the combination of P18 with the Pantheon catalog of SN Ia
in the redshift range $0.01<z<2.3$ [82]777https://github.com/dscolnic/Pantheon
(hereafter SN) using a prior on the supernovae peak absolute magnitude $M_{\rm
B}$ [83] [hereafter p(M)] of

$M_{\rm B}=-19.2435\pm 0.0373\;\rm mag.$ (55)

We sample on six standard parameters: $\omega_{\rm b}$, $\omega_{\rm c}$,
$H_{0}$, $\tau_{\rm reio}$, $\ln\left(10^{10}A_{s}\right)$, $n_{s}$, and the
modified gravity parameters, in particular

  * •

For IGGst we sample on $\tilde{\alpha}$ at fixed $\gamma=5\times 10^{-5}$.

  * •

For BDGph we sample on $1/\tilde{\alpha}_{8}$ with fixed
$\gamma~{}=~{}5~{}\times~{}10^{-5}$.

In addition to the cosmological parameters we sample on the nuisance and
foregrounds parameters of the P18 likelihoods and, when considering the
Phanteon dataset, on $M_{\rm B}$. In our analysis we assume two massless
neutrino with $N_{\rm eff}=2.0328$ and a massive one with fixed minimum mass
$m_{\nu}=0.06\rm eV$. As in [27], we fix the primordial ^4He mass fraction
$Y_{\rm p}$ taking into account the different value of the effective
gravitational constant during big bang nucleosynthesis, and the baryon
fraction, $\omega_{\rm b}$, tabulated in the public code PArthENoPE [84, 85].

Moreover, for each run we compute the best-fit values, obtained minimizing the
$\chi^{2}$ following the methods of [23], and quote the difference in the
model $\chi^{2}$ with respect to the $\Lambda\rm CDM$ one, i.e.
$\Delta\chi^{2}\equiv\chi^{2}-\chi^{2}_{\Lambda\rm CDM}$. Thus, negative
values of $\Delta\chi^{2}$ indicate an improvement in the fit with respect to
$\Lambda\rm CDM$. We also compute the Akaike information criterion (AIC) [86]
of the extended model $\mathcal{M}$ relative to that of $\Lambda$CDM:
$\Delta{\rm AIC}=\Delta\chi^{2}+2(N_{\mathcal{M}}-N_{\Lambda\rm CDM})$, where
$N_{\mathcal{M}}$ is the number of free parameters of the model.

###  IV.2 Results

####  IV.2.1 Induced gravity Galileon standard

For IGGst we present the results obtained by considering
$n=m=0\,,\gamma=5\times 10^{-5}$, with $\alpha$ allowed to vary. We show the
P18 results for this model in Figs. 9 and 1, with a comparison to $\Lambda$CDM
and IG with $\gamma=5\times 10^{-5}$. We obtain a tight constrain on the
Galileon term, i.e. $\tilde{\alpha}<2.5\times 10^{-6}$ at 95 % CL with P18
data. This result shows how CMB data are sufficient to strongly constrain IGG
around IG. All the contours for the other parameters overlap with the IG ones.

The marginalized mean and uncertainty for the Hubble constant $H_{0}\,[{\rm
km\,s^{-1}Mpc^{-1}}]$ at $68\%$ CL is $67.72\pm 0.54$ without any hint at a
possible reduction of the Hubble tension, as there is no significant increase
either in the mean or in the uncertainty with respect to the $\Lambda\rm CDM$
value of $67.36\pm 0.54$ [13], or the IGst value $67.64\pm 0.54$ (when
$\gamma=5\times 10^{-5}$ is fixed in the MCMC). When combining BAO with P18 we
obtain tighter constraints but the same qualitative behavior (Tables 1 and 2).

The cosmological constraints from P18 and P18+BAO on $\tilde{\alpha}$ allow a
Vainshtein mechanism to occur only at subparsec scale for an object of a solar
mass.

Figure 9: Marginalized $68\%$ and $95\%$ CL 2D regions using CMB data alone (P18). IGG standard ($Z=1$) with $\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda,\;g(\sigma)=\alpha$ in green; IG standard with $\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda$, in orange and $\Lambda\rm CDM$ in blue. | IGGst | IGst |  $\Lambda$CDM
---|---|---|---
$10^{2}\omega{}_{\rm b}$ | $2.242\pm 0.014$ | $2.240\pm 0.014$ | $2.237\pm 0.015$
$\omega{}_{\rm cdm}$ | $0.1200\pm 0.0012$ | $0.1200\pm 0.0012$ | $0.1200\pm 0.0012$
$H_{0}$ $[{\rm km\,s^{-1}Mpc^{-1}}]$ | $67.72\pm 0.54$ | $67.64\pm 0.54$ | $67.36\pm 0.54$
$\ln(10^{10}A_{\rm s})$ | $3.044\pm 0.014$ | $3.046\pm 0.014$ | $3.044\pm 0.014$
$n_{\rm s}$ | $0.9647\pm 0.0044$ | $0.9662\pm 0.0041$ | $0.9649\pm 0.0042$
$\tau{}_{\rm reio}$ | $0.0545\pm 0.0073$ | $0.0546\pm 0.0075$ | $0.0544\pm 0.0073$
$10^{6}\tilde{\alpha}{}$ |  $<2.5$ (95% CL) | — | —
$\gamma{}$ | $5\times 10^{-5}$ | $5\times 10^{-5}$ | —
$\Delta\chi^{2}$ | 0.2 | 0.2 | —
$\Delta\rm AIC$ | 2.2 | 0.2 | —
Table 1: Constraints on the main parameters (at $68\%$ CL unless otherwise
stated) considering _Planck_ 2018 data alone for IGG standard ($Z=1$) with
$\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda,\;g(\sigma)=\alpha$, IGG standard
with $\gamma=5\times 10^{-5},\;V(\sigma)=\Lambda$ and $\Lambda$CDM.

####  IV.2.2 Brans-Dicke Galileon phantom

For this model in the phantom branch with a cosmological constant, we add BAO
to P18 data in order to constrain the parameter $1/\tilde{\alpha}_{8}$ within
our prior range, $1/\tilde{\alpha}_{8}=10^{8}/\tilde{\alpha}\in[0,0.4]$.
Current CMB data are indeed not enough to constrain this seven parameter model
with $\gamma$ fixed. The presence of the cosmological constant helps in
providing posterior distributions closer to $\Lambda\rm CDM$ and similar
values for most of the cosmological parameters with the exception of $H_{0}$:
due to the shape of the posterior in the plane $1/\tilde{\alpha}_{8}$-$H_{0}$
(see Fig. 10), higher values of the Hubble constant compared to the
$\Lambda\rm CDM$ ones are allowed (as it happens for nonminimally coupled
scalar fields with standard kinetic terms due to the degeneracy
$\gamma$-$H_{0}$). The marginalized mean and uncertainty for the Hubble
constant $H_{0}\,[{\rm km\,s^{-1}Mpc^{-1}}]$ at $68\%$ CL using the
combination P18+BAO (Table 2) is $69.1^{+0.9}_{-1.3}$. For the same dataset,
IG phantom with $\gamma=5\times 10^{-5}$ fixed at the same value considered
here for BDGph, gives $67.62\pm 0.42$, while it gives $67.20^{+0.68}_{-0.55}$
when $\gamma$ is a free parameter. Thus, BDGph can raise the value of the
Hubble constant relatively to IGph and it reduces the Hubble tension to a
significance of $2.5\sigma$.

This result motivates our further analysis using the Pantheon catalog with the
addition of a prior on the peak absolute magnitude of SN Ia as in [87]. We
present the outcomes of this analysis in Tables 3 and 11. When considering
this dataset the Hubble tension has a significance of only $1.7\sigma$ in
BDGph, as the marginalized mean and uncertainty at $68\%$ CL is
$H_{0}=70.58\pm 0.97\,{\rm km\,s^{-1}\,Mpc^{-1}}$. The corresponding values in
IGph, with $\gamma=5\times 10^{-5}$ and $\gamma$ free to vary in the MCMC are,
respectively, $68.15\pm 0.52$ and $68.07\pm 0.56$. This means that the
Galileon term is not only necessary to avoid instabilities present in IGph
(see Sec. II.3) but it also plays an important role in alleviating the Hubble
tension. Using this dataset we obtain $\Delta\chi^{2}=-8.5$ and $\Delta\rm
AIC=-6.5$, showing a weak preference for the model with respect to
$\Lambda$CDM while in all other cases there is no significant improvement, if
any, with respect to $\Lambda$CDM.

We would also like to emphasize another aspect of our findings in terms of the
Vainshtein radius. We obtain $1/\tilde{\alpha}_{8}=0.23^{+0.06}_{-0.05}$ at
$68\%$ CL with P18 + SN + p(M); this means that we see high statistical
significance for $1/\tilde{\alpha}_{8}\neq 0$, and consequently a Vainshtein
radius of $\mathcal{O}(100)\,\rm pc$ for a solar mass. This is in contrast
with the case P18+BAO where $1/\tilde{\alpha}_{8}$ is consistent with zero at
$1\sigma$ and we only have an upper limit $1/\tilde{\alpha}_{8}<0.28$ at
$95\%$ CL, perfectly consistent with the $\Lambda\rm CDM$ limit of the theory.

Figure 10: Marginalized $68\%$ and $95\%$ CL 2D regions, using the combination P18 + BAO. BDG phantom ($Z=-1$) with $\gamma=5\times 10^{-5},\;\;g(\sigma)=\alpha\sigma^{-1}$ in red, IGph and IGph with fixed $\gamma=5\times 10^{-5}$ respectively in green and orange; in all cases $V(\sigma)=\Lambda$. $\Lambda\rm CDM$ in blue. | IGGst | BDGph | IGst |  IGph $\gamma=5\times 10^{-5}$ | IGph |  $\Lambda$CDM
---|---|---|---|---|---|---
$10^{2}\omega{}_{\rm b}$ | $2.245\pm 0.014$ | $2.244\pm 0.014$ | $2.244\pm 0.013$ | $2.245\pm 0.013$ | $2.245\pm 0.014$ | $2.244\pm 0.013$
$\omega{}_{\rm cdm}$ | $0.11940\pm 0.00096$ | $0.11966\pm 0.00097$ | $0.11937\pm 0.00093$ | $0.11910\pm 0.00093$ | $0.11884\pm 0.00099$ | $0.11925\pm 0.00094$
$H_{0}$ $[{\rm km\,s^{-1}Mpc^{-1}}]$ | $67.97\pm 0.44$ | $69.1^{+0.9}_{-1.3}$ | $67.91\pm 0.42$ | $67.62\pm 0.42$ | $67.20^{+0.68}_{-0.55}$ | $67.75\pm 0.43$
$\ln(10^{10}A_{\rm s})$ | $3.047\pm 0.014$ | $3.045\pm 0.014$ | $3.048\pm 0.014$ | $3.048^{+0.014}_{-0.015}$ | $3.048\pm 0.015$ | $3.049\pm 0.014$
$n_{\rm s}$ | $0.9662\pm 0.0041$ | $0.9678\pm 0.0038$ | $0.9676\pm 0.0036$ | $0.9675\pm 0.0037$ | $0.9671\pm 0.0038$ | $0.9675\pm 0.0037$
$\tau{}_{\rm reio}$ | $0.0562\pm 0.0073$ | $0.0543\pm 0.0073$ | $0.0563\pm 0.0071$ | $0.0568^{+0.0068}_{-0.0077}$ | $0.0575\pm 0.0075$ | $0.0568\pm 0.0072$
$10^{6}\tilde{\alpha}{}$ |  $<2.2$ (95% CL) | — | — | — | — | —
$1/\tilde{\alpha}_{8}{}$ | — |  $<0.28$ (95% CL) | — | — | — | —
$\gamma{}$ | $5\times 10^{-5}$ | $5\times 10^{-5}$ | $5\times 10^{-5}$ | $5\times 10^{-5}$ |  $<0.00047$ (95% CL) | —
$\Delta\chi^{2}$ | -0.1 | -1.1 | -0.1 | 0.2 | 0 | —
$\Delta\rm AIC$ | 1.9 | 0.9 | -0.1 | 0.2 | 2 | —
Table 2: Constraints on the main parameters (at $68\%$ CL unless otherwise stated) considering the combination P18+BAO for IGG standard ($Z=1$) with $g(\sigma)=\alpha$, IG standard and phantom ($Z=-1$) with $\gamma=5\times 10^{-5}$ fixed, IG phantom with $\gamma$ free to vary, BDG phantom with $g(\sigma)=\alpha\sigma^{-1}$ and $\Lambda$CDM. In all cases unless otherwise specified $\gamma=5\times 10^{-5}$ and $V(\sigma)=\Lambda$. Figure 11: Marginalized $68\%$ and $95\%$ CL 2D regions, using the combination P18 + SN + p(M). BDG phantom ($Z=-1$) with $\gamma=5\times 10^{-5},\;\;g(\sigma)=\alpha\sigma^{-1}$ in red, IGph and IGph with fixed $\gamma=5\times 10^{-5}$ respectively in green and orange; in all cases $V(\sigma)=\Lambda$. $\Lambda\rm CDM$ in blue. | BDGph | IGph $\gamma=5\times 10^{-5}$ | IGph |  $\Lambda$CDM
---|---|---|---|---
$10^{2}\omega{}_{\rm b}$ | $2.253\pm 0.014$ | $2.256\pm 0.014$ | $2.257\pm 0.014$ | $2.256\pm 0.014$
$\omega{}_{\rm cdm}$ | $0.1190\pm 0.0012$ | $0.1180\pm 0.0011$ | $0.1180\pm 0.0011$ | $0.1181\pm 0.0011$
$H_{0}$ $[{\rm km\,s^{-1}Mpc^{-1}}]$ | $70.58\pm 0.97$ | $68.15\pm 0.52$ | $68.07\pm 0.56$ | $68.31\pm 0.50$
$\ln(10^{10}A_{\rm s})$ | $3.046\pm 0.015$ | $3.053\pm 0.015$ | $3.054\pm 0.015$ | $3.054^{+0.014}_{-0.016}$
$n_{\rm s}$ | $0.9699\pm 0.0040$ | $0.9703\pm 0.0041$ | $0.9702\pm 0.0041$ | $0.9705\pm 0.0040$
$\tau{}_{\rm reio}$ | $0.0555\pm 0.0077$ | $0.0602^{+0.0076}_{-0.0085}$ | $0.0608\pm 0.0079$ | $0.0602^{+0.0072}_{-0.0083}$
$1/\tilde{\alpha}_{8}{}$ | $0.23^{+0.06}_{-0.05}$ | — | — | —
$\gamma{}$ | $5\times 10^{-5}$ | $5\times 10^{-5}$ |  $<0.00022$ (95% CL) | —
$\Delta\chi^{2}$ | -8.5 | 1.1 | 0 | —
$\Delta\rm AIC$ | -6.5 | 1.1 | 2 | —
Table 3: Constraints on the main parameters (at $68\%$ CL unless otherwise
stated) considering the combination P18+SN+p(M) for IG phantom ($Z=-1$) with
$\gamma=5\times 10^{-5}$ fixed and free to vary, BDG phantom with
$\gamma=5\times 10^{-5},\;g(\sigma)=\alpha\sigma^{-1}$ and $\Lambda$CDM. In
all cases $V(\sigma)=\Lambda$.

##  V Conclusions

We have studied the cosmological effects of a model where a scalar field,
$\sigma$, with a coupling to the Ricci scalar of the type
$F(\sigma)=\gamma\sigma^{2}$, has a cubic Galileon term of the form
$\alpha\,\sigma^{m}(\partial\sigma)^{2}\square\sigma$. Since the extensions of
IG and BD with a Galileon term are not equivalent theories up to a field
redefinition, differing by a term $4\zeta(\sigma)X^{2}$ in the Lagrangian, we
have studied the two models separately: IGG and BDG.

This Galileon term leads to the Vaishtein mechanism which can screen gravity,
potentially reconciling the theory with GR inside the so-called Vainshtein
radius also for values of the coupling to the Ricci scalar - or of the Brans-
Dicke parameter $\omega_{\mathrm{BD}}$ \- which evade the Solar System
constraints. Moreover, thanks to the presence of the Galileon term, the theory
is free of ghost and Laplacian instabilities even for a noncanonical sign of
the kinetic term in the Lagrangian. We have therefore considered the theory
(2) in the two branches, called respectively, the standard ($Z=1$) and phantom
($Z=-1$) branches, where $Z$ enters the kinetic term in the Lagrangian as
$ZX$, with $X=-\nabla_{\mu}\sigma\nabla^{\mu}\sigma/2$.

We have shown that the Galileon term leads to important cosmological effects
in all the cases considered. In this paper we have computed the theoretical
predictions for cosmological observables such as the CMB anisotropy and matter
power spectra. We have also compared these predictions with observations by
restricting the analysis to fixed $\gamma$ and $m$, if not otherwise stated.

For a standard kinetic term, we find an instability in the radiation dominated
era due to the presence of the Galileon term. This instability is dissipated
in the subsequent matter dominated era, during which the Galileon term becomes
subleading, and then eventually the scalar field energy density is no more
negligible contributing to the acceleration of the Universe. We find that the
CMB anisotropy pattern is sensitive to the dissipation of the instability in
the matter era. This effect constrains the Galileon term to be small close to
the CMB last scattering surface.

In the branch $Z=-1$ we have considered only BDG for the sake of simplicity.
The phenomenology of a kinetic term with a negative sign - corresponding to a
BD theory in the parameter range $\omega_{\mathrm{BD}}<0$ \- is rather
different. First, the presence of a Galileon term leads to an healthy theory
for all the values of $\gamma$, i.e. for any negative value of
$\omega_{\mathrm{BD}}$, therefore rescuing the range
$\omega_{\mathrm{BD}}<-3/2$ which would contain a ghost in the BD theory with
no Galileon term. The Galileon term leads to a dynamics very different from
the corresponding branch of induced gravity, by freezing the scalar field for
most of the matter dominated era and releasing it at lower redshift.

In the standard branch, Planck 2018 and BAO data constrain the amplitude of
the Galileon parameter as $\tilde{\alpha}<2.2\times 10^{-6}$ at 95 % CL (Table
2) for IGG with $\gamma=5\times 10^{-5}$, $m=0$, and a constant term $\Lambda$
as a potential. The constraints for BDG in the standard branch are identical
to the ones of IGG, we have therefore reported them only once for IGG. The
cosmological constraint on $\tilde{\alpha}$ allows a Vainshtein mechanism to
occur only at subparsec scale for an object of a solar mass. Summarizing, for
a standard kinetic term, we have therefore found tight constraints on the
Galileon term and resulting posterior probabilities for cosmological
parameters very similar to those of IG.

For BDG with $Z=-1$ and a constant term $\Lambda$ as a potential, we have
found that Planck 2018 and BAO data constrain $H_{0}=69.1^{+0.9}_{-1.3}\,{\rm
km\,s^{-1}\,Mpc^{-1}}$ at 68% CL and $\tilde{\alpha}^{-1}<0.281\times 10^{-8}$
at 95% CL, for $\gamma=5\times 10^{-5}$ and $m=-1$ (Tables 2 and 10). The
addition of a Galileon term in the phantom branch is therefore crucial to
obtain a value of $H_{0}$ larger than $\Lambda$CDM. We have therefore a
modified gravity model which leads to a value of $H_{0}$ larger than in
$\Lambda$CDM with a screening of $\mathcal{O}(100)\,\rm pc$ for a solar mass,
as desired. By adding the Pantheon dataset with a prior on the supernovae peak
absolute magnitude $M_{\rm B}=-19.2435\pm 0.0373\;\rm mag$, we find
$H_{0}=70.58\pm 0.97\,{\rm km\,s^{-1}\,Mpc^{-1}}$ and
$\tilde{\alpha}^{-1}=0.23^{+0.06}_{-0.05}\times 10^{-8}$ at 68% CL, always for
$\gamma=5\times 10^{-5}$ and $m=-1$ (Tables 3 and 11). As it happens for the
degeneracy between $H_{0}$ and $\gamma$ in IG [88], the SH0ES measurement
pulls $H_{0}$ along the degeneracy with $\alpha^{-1}$. By adding a Galileon
term, the theory inherits also characteristics of a late model in increasing
$H_{0}$. In addition, the value and the posterior of $S_{8}$ are unchanged
with respect to $\Lambda$CDM, not aggravating the so-called $\sigma_{8}$
tension [89].

We have also analyzed the theoretical predictions of the model introduced by
Silva and Koyama [50] in which there is no potential and the late-time
acceleration is driven exclusively by the Galileon term. This particular model
is physically viable and provides screening on Solar System scales. We have
shown that it leads to CMB predictions which are at odds with the Planck data,
with a $\Delta\chi^{2}=30.6$ with respect to $\Lambda$CDM. Therefore, although
theoretically interesting because the acceleration is not driven by an
effective cosmological constant, the model is ruled out by observations.

In conclusion, we have shown how a Galileon term leads to rather nontrivial
cosmological effects to the simplest scalar-tensor theories of gravity. It can
lead simultaneously to a value of $H_{0}$ larger than in $\Lambda$CDM and to
effective screening for a large volume of parameter space. It is therefore
interesting to find that another term in the Horndeski Lagrangian can lead to
a value larger then the $\Lambda$CDM one, as previously found for the basic
evolving Newton’s constant in [26, 27, 28, 68, 32, 33, 31, 88, 87, 69]. A more
general exploration of the parameter space of this theory in which the speed
of gravitational waves is $c$ without any fine tuning, by allowing to vary
simultaneously the coupling to the Ricci scalar, $\gamma$, and the Galileon
term is in progress.

## Acknowledgments

M.B., F.F., D.P. acknowledge financial support from the ASI/INAF Contract for
the Euclid mission n.2018-23-HH.0, from the INFN InDark initiative and from
the COSMOS network (www.cosmosnet.it) through the ASI (Italian Space Agency)
Grants No. 2016-24-H.0 and 2016-24-H.1-2018, as well as No. 2020-9-HH.0
(participation in LiteBIRD phase A).

This work has made use of computational resources of CNAF HPC cluster in
Bologna.

##  Appendix A BDG in the phantom branch with $\Lambda=0$

In this appendix we study the model, introduced in [50], that exactly respects
the solution (27), considering it as a limiting case of BDGph when
$\Lambda=0$. We discuss the background evolution, including stability
conditions and the Vainshtein screening of the theory, CMB anisotropies and
matter power spectrum and the constraints on the model coming from CMB
observations.

###  A.1 Background evolution, stability conditions and screening

In the absence of the cosmological constant, in order to describe today’s
acceleration of the expansion of the Universe, the parameter $\alpha$ should
be fine-tuned, this is done by solving Eq. 27 for $\alpha$ and using it as an
initial guess for a shooting algorithm that fixes the value of $\alpha$ to
produce the desired $\Omega_{\rm DE}=1-\Omega_{\rm m}-\Omega_{\rm r}$, given
the density parameters of matter and radiation as inputs. With $\alpha$ fixed
by the requirement of cosmic acceleration, the only free parameter of the
theory is $\gamma$; resulting therefore in a theory that can provide cosmic
acceleration without a cosmological constant, with as many free parameters as
the BD model.

Today’s value of the scalar field is such that it satisfies Eq. 28. See Secs.
II.2.2 and II.4 for details on why we do not need to set $\sigma_{0}$
following Eq. 25 or (24) thanks the Vainshtein screening mechanism, but we
still need to respect the condition $M_{\rm Pl}=1$ today. The evolution of the
scalar field with redshift is plotted in the top panel of Fig. 12. The field
is frozen deep in the radiation era and it starts to grow at very late times
reaching the value required by consistency with GR on small scales today. In
the figure we normalized the field to its value at $z=0$ which is different
depending on the chosen value of $\gamma$, since $\gamma\sigma_{0}^{2}=1$.

Figure 12: Time evolution of the scalar field (top panel), the equation of
state parameter for dark energy (middle panel) and the density parameters
(bottom panel) in BDG phantom ($Z=-1$) with no potential and
$g(\sigma)=\alpha\sigma^{-1}$, for different values of $\gamma$.

Departures from $\Lambda\rm CDM$ are evident already at the background level,
especially in the matter era and in the field dominated era. In fact, during
the matter era, before the equality with DE, we observe, from the bottom panel
of Fig. 12, $\widetilde{\Omega}_{\sigma}$ becoming slightly negative. This
effect is more prominent for larger values of $\gamma$. At the time when
$\widetilde{\Omega}_{\sigma}$ is decreasing and becoming negative, the
corresponding $\Omega_{\rm DE}$ in the $\Lambda\rm CDM$ model is already
growing, therefore, a steeper growth of $\widetilde{\Omega}_{\sigma}$ is
needed in order to reach the $\Lambda\rm CDM$ value today and for the matter-
dark energy equality to happen more or less at the same redshift. We emphasize
that $\widetilde{\Omega}_{\sigma}<0$ is not a problem from the physical point
of view: this parameter simply describes the contribution of the scalar field
to the total expansion rate when the Friedmann equations are recast in a form
that resembles Einstein gravity [52, 53].

The evolution of the parameter of state of DE $w_{\rm DE}$ from the middle
panel of Fig. 12 is also interesting, as it presents a phantom behavior
$w_{\rm DE}<-1$ in correspondence to $\widetilde{\Omega}_{\sigma}$ becoming
negative and growing back to dominate the energy content of the Universe. This
behavior is more prominent for a smaller value of $\gamma$ for which the dip
tends to be deeper and $w_{\rm DE}$ reaches lower values. The parameter
$w_{\rm DE}$ today is slightly different from $-1$ and it eventually reaches
$-1$ in the future, as shown in [50].

Figure 13: Stability conditions: ghost instability (top panel), Laplacian
instability (middle panel). Vainshtein radius for a solar mass spherical
object (bottom panel) BDG in the phantom branch ($Z=-1$), without potential,
$g(\sigma)=\alpha/\sigma$, for different values of $\gamma$.

The top two panels of Fig. 13 show the evolution of the stability conditions
(29) and (30) with redshift. The conditions are satisfied for all the values
of $\gamma$ considered thanks to the Galileon term. This term is also
responsible for effective screening on small scales through the Vainshtein
screening mechanism that allows the theory to reduce to GR within the so-
called Vainshtein radius. Thus, within this radius the post-Newtonian
parameters are those of GR. The time evolution of the Vainshtein radius (31)
for a celestial body of a solar mass is shown in the bottom panel of Fig. 13:
the screening is effective and $r_{V\odot}\sim\mathcal{O}(100)\,\rm pc$ today
for all the values of $\gamma$ studied.

###  A.2 CMB anisotropies and matter power spectrum

As we saw in Sec. III the interplay of the $G_{3}$ and the ISW effect can
result in an enhancement of the power spectrum on large scales [38, 40]. For
this reason, in Fig. 14 we can observe very large departures from $\Lambda\rm
CDM$ at the largest scales in TT, due to the enhanced ISW effect. It is worth
pointing out that as there is no potential, the curves do not flatten out
towards zero for smaller $\gamma$, this is due to the fact that without a
cosmological constant there is not a $\Lambda\rm CDM$ limit, as already
discussed in Sec. II. So, while this model provides late-time cosmic
acceleration without a cosmological constant, it shows large differences with
respect to $\Lambda\rm CDM$ at the level of CMB and matter power spectra; this
is a disadvantage for differences $>30\%$ when comparing the predictions
against the data, since $\Lambda\rm CDM$ fits Planck data quite well.

Figure 14: Relative differences with respect to $\Lambda\rm CDM$ for BDG
phantom ($Z=-1$) without potential and $g(\sigma)=\alpha\sigma^{-1}$, for
different values of $\gamma$. CMB TT (top left), EE (top right), lensing
potential (bottom left) and linear matter power spectra (bottom right).

###  A.3 Constraints from cosmological observations

P18
---
$10^{2}\omega{}_{\rm b}$ | $2.266\pm 0.015$
$\omega{}_{\rm cdm}$ | $0.1165\pm 0.0012$
$H_{0}$ $[{\rm kms^{-1}Mpc^{-1}}]$ | $79.57\pm 0.67$
$ln10^{10}A_{\rm s}$ | $2.9979^{+0.0093}_{-0.013}$
$n_{\rm s}$ | $0.9776\pm 0.0041$
$\tau{}_{\rm reio}$ |  $<0.0475\,$ (95% CL)
$\gamma$ |  $<7.70\cdot 10^{-7}\,$ (95% CL)
$\Delta\chi^{2}$ | 30.6
$\Delta\rm AIC$ | 32.6
Table 4: Constraints on the main parameters (at $68\%$ CL unless otherwise
stated) considering Planck 2018 data alone for BDG phantom ($Z=-1$) without
potential and $g(\sigma)=\alpha\sigma^{-1}$.

As we saw, this is the most extreme of the models considered, as it has no
cosmological constant and therefore no $\Lambda\rm CDM$ limit and it cannot
reproduce the CMB and LSS theoretical prediction of $\Lambda\rm CDM$, this put
it at a disadvantage when it is tested against data, and a
$\Delta\chi^{2}=30.6$ when considering P18 confirms it. While this completely
rules out the model, there is a feature that is worth highlighting in addition
to the theoretical appealing characteristic of providing cosmic acceleration
without a cosmological constant: it raises the Hubble constant to
$H_{0}=79.57\pm 0.67\,{\rm km\,s^{-1}\,Mpc^{-1}}$. While this result is still
far from alleviating the tension, the ability to produce a Hubble constant
larger than the SH0ES value [14] is quite interesting and it might be worth to
construct and investigate similar models that, while retaining this ability,
might provide a better fit to the data. The results for the model are shown in
Fig. 15 and Table 4.

Figure 15: Marginalized joint $68\%$ and $95\%$ CL regions 2D parameter space
using P18 data alone. In orange BDG phantom ($Z=-1$) without potential,
$g(\sigma)=\alpha\sigma^{-1}$ and $\Lambda\rm CDM$ in blue.

##  Appendix B Linear perturbed equations

Here we summarize the definition of the coefficients appearing in the
perturbed Einstein and scalar field equations:

$\displaystyle\delta\tilde{\rho}$
$\displaystyle\equiv\frac{\delta\rho_{m}}{\gamma\sigma^{2}}-\frac{h^{\prime}\sigma^{\prime}}{a^{2}\sigma}-\frac{2}{a^{2}}\Bigg{\\{}\frac{\delta\sigma}{\sigma}\bigg{[}\frac{a^{2}\rho_{m}}{\gamma\sigma^{2}}+Z\frac{\sigma^{\prime\,2}}{2\gamma\sigma^{2}}+\frac{a^{2}}{\gamma\sigma}\bigg{(}\frac{V}{\sigma}-\frac{V_{,\sigma}}{2}\bigg{)}-\frac{3\mathcal{H}\sigma^{\prime}}{\sigma}+k^{2}\bigg{]}+\frac{\delta\sigma^{\prime}}{\sigma}\bigg{(}\,3\mathcal{H}-Z\frac{\sigma^{\prime}}{2\gamma\sigma}\bigg{)}\Bigg{\\}}\,,$
(56) $\displaystyle\delta\tilde{\rho}^{(G)}$
$\displaystyle\equiv\,-\frac{2}{\gamma
a^{4}}\frac{\sigma^{\prime\,2}}{\sigma^{2}}\Bigg{\\{}\frac{\delta\sigma}{\sigma}\bigg{[}3\mathcal{H}\sigma^{\prime}\big{(}2\,g-\sigma
g_{,\sigma}\big{)}+\frac{\sigma^{\prime\,2}}{2}\big{(}\sigma
g_{,\sigma,\sigma}-2g_{,\sigma}\big{)}-k^{2}\sigma
g+\frac{3\sigma^{\prime\,2}}{2}\big{(}2\zeta-\zeta_{,\sigma}\sigma\big{)}\bigg{]}$
$\displaystyle+\delta\sigma^{\prime}\Big{(}\,\sigma^{\prime}(2g_{,\sigma}-6\zeta)-9\mathcal{H}g\,\Big{)}-\frac{1}{2}h^{\prime}\sigma^{\prime}g\Bigg{\\}}\,,$
(57) $\displaystyle\big{(}\tilde{\rho}+\tilde{P}\big{)}\tilde{\theta}\,$
$\displaystyle\equiv\,\frac{(\rho_{m}+P_{m})}{\gamma\sigma^{2}}\theta_{m}+\frac{2k^{2}}{a^{2}}\bigg{\\{}\,\frac{\delta\sigma}{\sigma}\bigg{[}\,\frac{\sigma^{\prime}}{2\gamma\sigma}(Z+2\gamma)-\mathcal{H}\bigg{]}+\frac{\delta\sigma^{\prime}}{\sigma}\bigg{\\}}\,,$
(58)
$\displaystyle\big{(}\tilde{\rho}^{(G)}+\tilde{P}^{(G)}\big{)}\tilde{\theta}^{(G)}$
$\displaystyle\equiv\frac{2k^{2}}{\gamma
a^{4}}\frac{\sigma^{\prime\,2}}{\sigma^{2}}\,\Big{[}\delta\sigma(\,3\mathcal{H}g-g_{,\sigma}\sigma^{\prime}+2\zeta\sigma^{\prime}\,)-\delta\sigma^{\prime}g\Big{]}\,,$
(59) $\displaystyle\delta\tilde{P}$ $\displaystyle\equiv\frac{\delta
P_{m}}{\gamma\sigma^{2}}+\frac{1}{a^{2}}\Bigg{\\{}-\frac{\delta\sigma}{\sigma}\bigg{[}\frac{a^{2}}{\gamma\sigma}\bigg{(}V_{,\sigma}-\frac{2V}{\sigma}\bigg{)}+\frac{2a^{2}P_{m}}{\gamma\sigma^{2}}+\frac{\sigma^{\prime\,2}}{\gamma\sigma^{2}}(Z+4\gamma)+\frac{2\sigma^{\prime\prime}}{\sigma}+\frac{2\mathcal{H}\sigma^{\prime}}{\sigma}-\frac{4k^{2}}{3}\bigg{]}$
$\displaystyle+\delta\sigma^{\prime}\bigg{(}\frac{\sigma^{\prime}}{\gamma\sigma^{2}}(Z+4\gamma)+\frac{2\mathcal{H}}{\sigma}\bigg{)}+\frac{2\,\delta\sigma^{\prime\prime}}{\sigma}+\frac{2\sigma^{\prime}}{3\sigma}h^{\prime}\Bigg{\\}}\,,$
(60) $\displaystyle\delta\tilde{P}^{(G)}$ $\displaystyle\equiv\frac{1}{\gamma
a^{4}}\frac{\sigma^{\prime\,2}}{\sigma^{2}}\,\Bigg{\\{}\frac{\delta\sigma}{\sigma}\Big{[}2\,(2\,g-g_{,\sigma}\sigma)(\sigma^{\prime\prime}-\mathcal{H}\sigma^{\prime})+\sigma^{\prime\,2}\big{(}2\,g_{,\sigma}\,-g_{,\sigma,\sigma}\,\sigma-\,2\zeta\,+\zeta_{,\sigma}\,\sigma\big{)}\Big{]}$
$\displaystyle-\delta\sigma^{\prime}\Big{(}\,2g_{,\sigma}\sigma^{\prime}-4\mathcal{H}\,g+2\,g\frac{\sigma^{\prime\prime}}{\sigma^{\prime}}+4\,\zeta\,\sigma^{\prime}\,\Big{)}-2\,g\,\delta\sigma^{\prime\prime}\Bigg{\\}}\,,$
(61) $\displaystyle(\tilde{\rho}+\tilde{P})\tilde{\Theta}$
$\displaystyle\equiv\frac{(\rho_{m}+P_{m})\Theta}{\gamma\sigma^{2}}+\frac{1}{3a^{2}}\bigg{[}\,\frac{4k^{2}\delta\sigma}{\sigma}+2(h^{\prime}+6\eta^{\prime})\frac{\sigma^{\prime}}{\sigma}\,\bigg{]}\,,$
(62) $\displaystyle\delta\tilde{\cal G}$
$\displaystyle\equiv\frac{a^{2}(\delta\rho_{m}-3\,\delta
P_{m})}{(1+6\gamma)\sigma}\,-\frac{h^{\prime}}{2a^{2}}\left[a^{2}\sigma^{\prime}(Z+6\gamma)+2\,g\sigma^{\prime}\left(3\mathcal{H}\sigma^{\prime}-\frac{\sigma^{\prime
2}}{\sigma}+2g\sigma^{\prime\prime}\right)+4\zeta\sigma^{\prime
3}\right]-h^{\prime\prime}g\frac{\sigma^{\prime 2}}{a^{2}}$
$\displaystyle+\frac{1}{a^{2}}\delta\sigma\Bigg{\\{}a^{4}\left[\frac{3P_{m}-\rho_{m}}{\sigma^{2}}-\frac{4V}{\sigma^{2}}+\frac{4V_{,\sigma}}{\sigma}-V_{,\sigma\sigma}+(Z+6\gamma)\left(\frac{\sigma^{\prime
2}}{\sigma^{2}}-k^{2}\right)\right]$
$\displaystyle\;+g\left[2k^{2}\left(\frac{\sigma^{\prime\,2}}{\sigma}-2\mathcal{H}\sigma^{\prime}-2\sigma^{\prime\prime}\right)-6\frac{\sigma^{\prime\,2}}{\sigma^{2}}\sigma^{\prime\prime}\right]-4k^{2}\zeta\,\sigma^{\prime\,2}+g_{,\sigma}\,\sigma^{\prime}\left[6\sigma^{\prime\prime}\left(\frac{\sigma^{\prime}}{\sigma}-2\mathcal{H}\right)-2\sigma^{\prime}\left(3\mathcal{H}^{\prime}+\frac{\sigma^{\prime\,2}}{\sigma^{2}}\right)\right]$
$\displaystyle\;-12\sigma^{\prime\,2}\zeta_{,\sigma}\sigma^{\prime\prime}+2g_{,\sigma\sigma}\sigma^{\prime\,2}\left(\frac{\sigma^{\prime\,2}}{\sigma}-2\mathcal{H}\sigma^{\prime}+2\sigma^{\prime\prime}\right)-3\zeta_{,\sigma\sigma}\sigma^{\prime\,4}-g_{,\sigma\sigma\sigma}\sigma^{\prime\,4}\Bigg{\\}}$
$\displaystyle+\frac{1}{a^{2}}\delta\sigma^{\prime}\Bigg{\\{}a^{2}\left[-2(Z+6\gamma)\left(\mathcal{H}+\frac{\sigma^{\prime}}{\sigma}\right)\right]+12\,g\left[\sigma^{\prime\prime}\left(\frac{\sigma^{\prime}}{\sigma}-\mathcal{H}\right)-\mathcal{H}^{\prime}\sigma^{\prime}\right]-24\zeta\sigma^{\prime}\sigma^{\prime\prime}$
$\displaystyle\qquad+4\,g_{,\sigma}\sigma^{\prime}\left(2\frac{\sigma^{\prime\,2}}{\sigma}-3\mathcal{H}\sigma^{\prime}+2\sigma^{\prime\prime}\right)-12\zeta_{,\sigma}\sigma^{\prime\,3}+4\,g_{,\sigma\sigma}\sigma^{\prime\,3}\Bigg{\\}}\,.$
(63)

## References

  * Aghanim  _et al._ [2020a] N. Aghanim _et al._ (Planck), Planck 2018 results. I. Overview and the cosmological legacy of Planck, Astron. Astrophys. 641, A1 (2020a), arXiv:1807.06205 [astro-ph.CO] .
  * Alam  _et al._ [2021] S. Alam _et al._ (eBOSS), Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory, Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 [astro-ph.CO] .
  * Veropalumbo  _et al._ [2014] A. Veropalumbo, F. Marulli, L. Moscardini, M. Moresco, and A. Cimatti, An improved measurement of baryon acoustic oscillations from the correlation function of galaxy clusters at z $\sim$ 0.3, Mon. Not. Roy. Astron. Soc. 442, 3275 (2014), arXiv:1311.5895 [astro-ph.CO] .
  * Asgari  _et al._ [2021] M. Asgari _et al._ (KiDS), KiDS-1000 Cosmology: Cosmic shear constraints and comparison between two point statistics, Astron. Astrophys. 645, A104 (2021), arXiv:2007.15633 [astro-ph.CO] .
  * Abbott  _et al._ [2022] T. M. C. Abbott _et al._ (DES), Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing, Phys. Rev. D 105, 023520 (2022), arXiv:2105.13549 [astro-ph.CO] .
  * Aghanim  _et al._ [2020b] N. Aghanim _et al._ (Planck), Planck 2018 results. VIII. Gravitational lensing, Astron. Astrophys. 641, A8 (2020b), arXiv:1807.06210 [astro-ph.CO] .
  * Sherwin  _et al._ [2017] B. D. Sherwin, A. van Engelen, N. Sehgal, M. Madhavacheril, G. E. Addison, S. Aiola, _et al._ , Two-season atacama cosmology telescope polarimeter lensing power spectrum, Phys. Rev. D 95, 123529 (2017), arXiv:1611.09753 [astro-ph.CO] .
  * Wu  _et al._ [2019] W. L. K. Wu _et al._ , A Measurement of the Cosmic Microwave Background Lensing Potential and Power Spectrum from 500 deg2 of SPTpol Temperature and Polarization Data, Astrophys. J. 884, 70 (2019), arXiv:1905.05777 [astro-ph.CO] .
  * Darwish  _et al._ [2020] O. Darwish _et al._ , The Atacama Cosmology Telescope: A CMB lensing mass map over 2100 square degrees of sky and its cross-correlation with BOSS-CMASS galaxies, Mon. Not. Roy. Astron. Soc. 500, 2250 (2020), arXiv:2004.01139 [astro-ph.CO] .
  * Riess  _et al._ [1998] A. G. Riess _et al._ (Supernova Search Team), Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ [1999] S. Perlmutter _et al._ (Supernova Cosmology Project), Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae, Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Carr  _et al._ [2022] A. Carr, T. M. Davis, D. Scolnic, D. Scolnic, K. Said, D. Brout, E. R. Peterson, and R. Kessler, The Pantheon+ analysis: Improving the redshifts and peculiar velocities of Type Ia supernovae used in cosmological analyses, Publ. Astron. Soc. Austral. 39, e046 (2022), arXiv:2112.01471 [astro-ph.CO] .
  * Aghanim  _et al._ [2020c] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020c), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Riess  _et al._ [2022] A. G. Riess _et al._ , A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s-1 Mpc-1 Uncertainty from the Hubble Space Telescope and the SH0ES Team, Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Aiola  _et al._ [2020] S. Aiola _et al._ (ACT), The Atacama Cosmology Telescope: DR4 Maps and Cosmological Parameters, JCAP 12, 047, arXiv:2007.07288 [astro-ph.CO] .
  * Balkenhol  _et al._ [2021] L. Balkenhol _et al._ (SPT-3G), Constraints on $\Lambda$CDM extensions from the SPT-3G 2018 EE and TE power spectra, Phys. Rev. D 104, 083509 (2021), arXiv:2103.13618 [astro-ph.CO] .
  * Freedman  _et al._ [2019] W. L. Freedman _et al._ , The Carnegie-Chicago Hubble Program. VIII. An Independent Determination of the Hubble Constant Based on the Tip of the Red Giant Branch, Astrophys. J. 882, 34 (2019), arXiv:1907.05922 [astro-ph.CO] .
  * Ade  _et al._ [2014] P. A. R. Ade _et al._ (Planck), Planck 2013 results. XVI. Cosmological parameters, Astron. Astrophys. 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Knox and Millea [2020] L. Knox and M. Millea, Hubble constant hunter’s guide, Phys. Rev. D 101, 043533 (2020), arXiv:1908.03663 [astro-ph.CO] .
  * Di Valentino  _et al._ [2021a] E. Di Valentino _et al._ , Snowmass2021 - Letter of interest cosmology intertwined II: The hubble constant tension, Astropart. Phys. 131, 102605 (2021a), arXiv:2008.11284 [astro-ph.CO] .
  * Di Valentino  _et al._ [2021b] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, In the realm of the Hubble tension—a review of solutions, Class. Quant. Grav. 38, 153001 (2021b), arXiv:2103.01183 [astro-ph.CO] .
  * Perivolaropoulos and Skara [2022] L. Perivolaropoulos and F. Skara, Challenges for $\Lambda$CDM: An update, New Astron. Rev. 95, 10.1016/j.newar.2022.101659 (2022), arXiv:2105.05208 [astro-ph.CO] .
  * Schöneberg  _et al._ [2022] N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, The H0 Olympics: A fair ranking of proposed models, Phys. Rept. 984, 1 (2022), arXiv:2107.10291 [astro-ph.CO] .
  * Shah  _et al._ [2021] P. Shah, P. Lemos, and O. Lahav, A buyer’s guide to the Hubble constant, Astron. Astrophys. Rev. 29, 9 (2021), arXiv:2109.01161 [astro-ph.CO] .
  * Abdalla  _et al._ [2022] E. Abdalla _et al._ , Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies, JHEAp 34, 49 (2022), arXiv:2203.06142 [astro-ph.CO] .
  * Umiltà  _et al._ [2015] C. Umiltà, M. Ballardini, F. Finelli, and D. Paoletti, CMB and BAO constraints for an induced gravity dark energy model with a quartic potential, JCAP 08, 017, arXiv:1507.00718 [astro-ph.CO] .
  * Ballardini  _et al._ [2016] M. Ballardini, F. Finelli, C. Umiltà, and D. Paoletti, Cosmological constraints on induced gravity dark energy models, JCAP 05, 067, arXiv:1601.03387 [astro-ph.CO] .
  * Rossi  _et al._ [2019] M. Rossi, M. Ballardini, M. Braglia, F. Finelli, D. Paoletti, A. A. Starobinsky, and C. Umiltà, Cosmological constraints on post-Newtonian parameters in effectively massless scalar-tensor theories of gravity, Phys. Rev. D 100, 103524 (2019), arXiv:1906.10218 [astro-ph.CO] .
  * Ballesteros  _et al._ [2020] G. Ballesteros, A. Notari, and F. Rompineve, The $H_{0}$ tension: $\Delta G_{N}$ vs. $\Delta N_{\rm eff}$, JCAP 11, 024, arXiv:2004.05049 [astro-ph.CO] .
  * Zumalacarregui [2020] M. Zumalacarregui, Gravity in the Era of Equality: Towards solutions to the Hubble problem without fine-tuned initial conditions, Phys. Rev. D 102, 023523 (2020), arXiv:2003.06396 [astro-ph.CO] .
  * Braglia  _et al._ [2020] M. Braglia, M. Ballardini, W. T. Emond, F. Finelli, A. E. Gumrukcuoglu, K. Koyama, and D. Paoletti, Larger value for $H_{0}$ by an evolving gravitational constant, Phys. Rev. D 102, 023529 (2020), arXiv:2004.11161 [astro-ph.CO] .
  * Ballardini  _et al._ [2020] M. Ballardini, M. Braglia, F. Finelli, D. Paoletti, A. A. Starobinsky, and C. Umiltà, Scalar-tensor theories of gravity, neutrino physics, and the $H_{0}$ tension, JCAP 10, 044, arXiv:2004.14349 [astro-ph.CO] .
  * Braglia  _et al._ [2021] M. Braglia, M. Ballardini, F. Finelli, and K. Koyama, Early modified gravity in light of the $H_{0}$ tension and LSS data, Phys. Rev. D 103, 043528 (2021), arXiv:2011.12934 [astro-ph.CO] .
  * Lee  _et al._ [2022] B.-H. Lee, W. Lee, E. O. Colgáin, M. M. Sheikh-Jabbari, and S. Thakur, Is local H 0 at odds with dark energy EFT?, JCAP 04 (04), 004, arXiv:2202.03906 [astro-ph.CO] .
  * Bertotti  _et al._ [2003] B. Bertotti, L. Iess, and P. Tortora, A test of general relativity using radio links with the Cassini spacecraft, Nature 425, 374 (2003).
  * Nicolis  _et al._ [2009] A. Nicolis, R. Rattazzi, and E. Trincherini, The Galileon as a local modification of gravity, Phys. Rev. D 79, 064036 (2009), arXiv:0811.2197 [hep-th] .
  * Deffayet  _et al._ [2009] C. Deffayet, G. Esposito-Farese, and A. Vikman, Covariant Galileon, Phys. Rev. D 79, 084003 (2009), arXiv:0901.1314 [hep-th] .
  * Kobayashi  _et al._ [2010] T. Kobayashi, H. Tashiro, and D. Suzuki, Evolution of linear cosmological perturbations and its observational implications in Galileon-type modified gravity, Phys. Rev. D 81, 063513 (2010), arXiv:0912.4641 [astro-ph.CO] .
  * Barreira  _et al._ [2014] A. Barreira, B. Li, C. Baugh, and S. Pascoli, The observational status of Galileon gravity after Planck, JCAP 08, 059, arXiv:1406.0485 [astro-ph.CO] .
  * Renk  _et al._ [2016] J. Renk, M. Zumalacarregui, and F. Montanari, Gravity at the horizon: on relativistic effects, CMB-LSS correlations and ultra-large scales in Horndeski’s theory, JCAP 07, 040, arXiv:1604.03487 [astro-ph.CO] .
  * Frusciante  _et al._ [2020] N. Frusciante, S. Peirone, L. Atayde, and A. De Felice, Phenomenology of the generalized cubic covariant Galileon model and cosmological bounds, Phys. Rev. D 101, 064001 (2020), arXiv:1912.07586 [astro-ph.CO] .
  * Horndeski [1974] G. W. Horndeski, Second-order scalar-tensor field equations in a four-dimensional space, Int. J. Theor. Phys. 10, 363 (1974).
  * Deffayet  _et al._ [2011] C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade, From k-essence to generalised Galileons, Phys. Rev. D 84, 064039 (2011), arXiv:1103.3260 [hep-th] .
  * Kobayashi  _et al._ [2011] T. Kobayashi, M. Yamaguchi, and J. Yokoyama, Generalized G-inflation: Inflation with the most general second-order field equations, Prog. Theor. Phys. 126, 511 (2011), arXiv:1105.5723 [hep-th] .
  * Kase and Tsujikawa [2019] R. Kase and S. Tsujikawa, Dark energy in Horndeski theories after GW170817: A review, Int. J. Mod. Phys. D 28, 1942005 (2019), arXiv:1809.08735 [gr-qc] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A, Astrophys. J. Lett. 848, L13 (2017), arXiv:1710.05834 [astro-ph.HE] .
  * Jordan [1955] P. Jordan,  _Schwerkraft und Weltall: Grundlagen der theoretischen Kosmologie_, Die Wissenschaft (F. Vieweg, 1955).
  * Brans and Dicke [1961] C. Brans and R. H. Dicke, Mach’s principle and a relativistic theory of gravitation, Phys. Rev. 124, 925 (1961).
  * Zee [1979] A. Zee, Broken-symmetric theory of gravity, Phys. Rev. Lett. 42, 417 (1979).
  * Silva and Koyama [2009] F. P. Silva and K. Koyama, Self-accelerating universe in galileon cosmology, Phys. Rev. D 80, 121301 (2009), arXiv:0909.4538 [astro-ph.CO] .
  * Finelli  _et al._ [2008] F. Finelli, A. Tronconi, and G. Venturi, Dark Energy, Induced Gravity and Broken Scale Invariance, Phys. Lett. B 659, 466 (2008), arXiv:0710.2741 [astro-ph] .
  * Torres [2002] D. F. Torres, Quintessence, superquintessence and observable quantities in Brans-Dicke and nonminimally coupled theories, Phys. Rev. D 66, 043522 (2002), arXiv:astro-ph/0204504 .
  * Gannouji  _et al._ [2006] R. Gannouji, D. Polarski, A. Ranquet, and A. A. Starobinsky, Scalar-Tensor Models of Normal and Phantom Dark Energy, JCAP 09, 016, arXiv:astro-ph/0606287 .
  * Cerioni [2011] A. Cerioni,  _Cosmological perturbations in generalized theories of gravity_ , Ph.D. thesis, Alma Mater Studiorum - Università di Bologna (2011).
  * Hohmann [2015] M. Hohmann, Parametrized post-Newtonian limit of Horndeski’s gravity theory, Phys. Rev. D 92, 064019 (2015), arXiv:1506.04253 [gr-qc] .
  * Quiros  _et al._ [2020] I. Quiros, R. De Arcia, R. García-Salcedo, T. Gonzalez, and F. A. Horta-Rangel, An issue with the classification of the scalar–tensor theories of gravity, Int. J. Mod. Phys. D 29, 2050047 (2020), arXiv:1905.08177 [gr-qc] .
  * Boisseau  _et al._ [2000] B. Boisseau, G. Esposito-Farese, D. Polarski, and A. A. Starobinsky, Reconstruction of a scalar tensor theory of gravity in an accelerating universe, Phys. Rev. Lett. 85, 2236 (2000), arXiv:gr-qc/0001066 .
  * Chow and Khoury [2009] N. Chow and J. Khoury, Galileon Cosmology, Phys. Rev. D 80, 024037 (2009), arXiv:0905.1325 [hep-th] .
  * Kimura  _et al._ [2012] R. Kimura, T. Kobayashi, and K. Yamamoto, Vainshtein screening in a cosmological background in the most general second-order scalar-tensor theory, Phys. Rev. D 85, 024023 (2012), arXiv:1111.6749 [astro-ph.CO] .
  * Bellini and Sawicki [2014] E. Bellini and I. Sawicki, Maximal freedom at minimum cost: linear large-scale structure in general modifications of gravity, JCAP 07, 050, arXiv:1404.3713 [astro-ph.CO] .
  * De Felice and Tsujikawa [2012] A. De Felice and S. Tsujikawa, Conditions for the cosmological viability of the most general scalar-tensor theories and their applications to extended Galileon dark energy models, JCAP 02, 007, arXiv:1110.3878 [gr-qc] .
  * Babichev  _et al._ [2008] E. Babichev, V. Mukhanov, and A. Vikman, k-Essence, superluminal propagation, causality and emergent geometry, JHEP 02, 101, arXiv:0708.0561 [hep-th] .
  * Gorini  _et al._ [2008] V. Gorini, A. Y. Kamenshchik, U. Moschella, O. F. Piattella, and A. A. Starobinsky, Gauge-invariant analysis of perturbations in Chaplygin gas unified models of dark matter and dark energy, JCAP 02, 016, arXiv:0711.4242 [astro-ph] .
  * Adams  _et al._ [2006] A. Adams, N. Arkani-Hamed, S. Dubovsky, A. Nicolis, and R. Rattazzi, Causality, analyticity and an IR obstruction to UV completion, JHEP 10, 014, arXiv:hep-th/0602178 .
  * Bonvin  _et al._ [2006] C. Bonvin, C. Caprini, and R. Durrer, A no-go theorem for k-essence dark energy, Phys. Rev. Lett. 97, 081303 (2006), arXiv:astro-ph/0606584 .
  * Ellis  _et al._ [2007] G. Ellis, R. Maartens, and M. A. H. MacCallum, Causality and the speed of sound, Gen. Rel. Grav. 39, 1651 (2007), arXiv:gr-qc/0703121 .
  * Ma and Bertschinger [1995] C.-P. Ma and E. Bertschinger, Cosmological perturbation theory in the synchronous and conformal newtonian gauges, The Astrophysical Journal 455, 7 (1995).
  * Paoletti  _et al._ [2019] D. Paoletti, M. Braglia, F. Finelli, M. Ballardini, and C. Umiltà, Isocurvature fluctuations in the effective Newton’s constant, Phys. Dark Univ. 25, 100307 (2019), arXiv:1809.03201 [astro-ph.CO] .
  * Ballardini  _et al._ [2023] M. Ballardini, A. G. Ferrari, and F. Finelli, Phantom scalar-tensor models and cosmological tensions, JCAP 04, 029, arXiv:2302.05291 [astro-ph.CO] .
  * Lesgourgues [2011] J. Lesgourgues, The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview, arXiv e-prints , arXiv:1104.2932 (2011), arXiv:1104.2932 [astro-ph.IM] .
  * Audren  _et al._ [2013] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, Conservative constraints on early cosmology: an illustration of the monte python cosmological parameter inference code, Journal of Cosmology and Astroparticle Physics 2013 (02), 001.
  * Brinckmann and Lesgourgues [2019] T. Brinckmann and J. Lesgourgues, MontePython 3: boosted MCMC sampler and other features, Phys. Dark Univ. 24, 100260 (2019), arXiv:1804.07261 [astro-ph.CO] .
  * Gelman and Rubin [1992] A. Gelman and D. B. Rubin, Inference from Iterative Simulation Using Multiple Sequences, Statist. Sci. 7, 457 (1992).
  * Lewis [2019] A. Lewis, GetDist: a Python package for analysing Monte Carlo samples, arXiv e-prints , arXiv:1910.13970 (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Aghanim  _et al._ [2020d] N. Aghanim _et al._ (Planck), Planck 2018 results. V. CMB power spectra and likelihoods, Astron. Astrophys. 641, A5 (2020d), arXiv:1907.12875 [astro-ph.CO] .
  * Alam  _et al._ [2017] S. Alam _et al._ (BOSS), The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, Mon. Not. Roy. Astron. Soc. 470, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * Beutler  _et al._ [2011] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, The 6df galaxy survey: baryon acoustic oscillations and the local hubble constant, Monthly Notices of the Royal Astronomical Society 416, 3017 (2011).
  * Ross  _et al._ [2015] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, The clustering of the SDSS DR7 main Galaxy sample – I. A 4 per cent distance measure at $z=0.15$, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * de Sainte Agathe  _et al._ [2019] V. de Sainte Agathe _et al._ , Baryon acoustic oscillations at z = 2.34 from the correlations of Ly$\alpha$ absorption in eBOSS DR14, Astron. Astrophys. 629, A85 (2019), arXiv:1904.03400 [astro-ph.CO] .
  * Blomqvist  _et al._ [2019] M. Blomqvist _et al._ , Baryon acoustic oscillations from the cross-correlation of Ly$\alpha$ absorption and quasars in eBOSS DR14, Astron. Astrophys. 629, A86 (2019), arXiv:1904.03430 [astro-ph.CO] .
  * Cuceu  _et al._ [2019] A. Cuceu, J. Farr, P. Lemos, and A. Font-Ribera, Baryon Acoustic Oscillations and the Hubble Constant: Past, Present and Future, JCAP 10, 044, arXiv:1906.11628 [astro-ph.CO] .
  * Scolnic  _et al._ [2018] D. M. Scolnic _et al._ (Pan-STARRS1), The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample, Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Camarena and Marra [2021] D. Camarena and V. Marra, On the use of the local prior on the absolute magnitude of Type Ia supernovae in cosmological inference, Mon. Not. Roy. Astron. Soc. 504, 5164 (2021), arXiv:2101.08641 [astro-ph.CO] .
  * Pisanti  _et al._ [2008] O. Pisanti, A. Cirillo, S. Esposito, F. Iocco, G. Mangano, G. Miele, and P. D. Serpico, PArthENoPE: Public Algorithm Evaluating the Nucleosynthesis of Primordial Elements, Comput. Phys. Commun. 178, 956 (2008), arXiv:0705.0290 [astro-ph] .
  * Consiglio  _et al._ [2018] R. Consiglio, P. F. de Salas, G. Mangano, G. Miele, S. Pastor, and O. Pisanti, PArthENoPE reloaded, Comput. Phys. Commun. 233, 237 (2018), arXiv:1712.04378 [astro-ph.CO] .
  * Akaike [1974] H. Akaike, A new look at the statistical model identification, IEEE Transactions on Automatic Control 19, 716 (1974).
  * Ballardini  _et al._ [2022] M. Ballardini, F. Finelli, and D. Sapone, Cosmological constraints on the gravitational constant, JCAP 06 (06), 004, arXiv:2111.09168 [astro-ph.CO] .
  * Ballardini and Finelli [2022] M. Ballardini and F. Finelli, Type Ia supernovae data with scalar-tensor gravity, Phys. Rev. D 106, 063531 (2022), arXiv:2112.15126 [astro-ph.CO] .
  * Di Valentino  _et al._ [2021c] E. Di Valentino _et al._ , Cosmology intertwined III: $f\sigma_{8}$ and $S_{8}$, Astropart. Phys. 131, 102604 (2021c), arXiv:2008.11285 [astro-ph.CO] .