# Phantom scalar-tensor models and cosmological tensions

Mario Ballardini  Angelo G. Ferrari  Fabio Finelli

###### Abstract

We study three different extended scalar-tensor theories of gravity by also
allowing a negative sign for the kinetic term for the scalar field in the
Jordan frame. Our scope is to understand how the observational constraints for
these models cope with the volume of the parameter space in which the theory
is healthy. Models with a negative kinetic term lead to decreasing effective
gravitational constant with redshift and behave as an effective relativistic
component with a negative energy density as opposite to their corresponding
version with a standard kinetic term. As a consequence, we find that the
extended branch with a negative sign for the kinetic term correspond in
general to lower $H_{0}$ and $\sigma_{8}$ compared to $\Lambda$CDM. We find
that in all the cases with a negative sign for the kinetic term studied here,
cosmological observations constrain these models around GR and prefer a volume
of the parameter space in which the theory is not healthy since the scalar
field behave as a ghost also in the related Einstein frame. We show that also
in the phantom branch early modify gravity with a quartic coupling can
substantially reduce the $H_{0}$ tension fitting the combination of cosmic
microwave background data from Planck, baryon acoustic oscillations from BOSS
and eBOSS, and Supernovae from the Pantheon sample with calibration
information by SH0ES.

##  1 Introduction

The $\Lambda$CDM model represents the current standard cosmological model
providing an excellent fit to most of cosmological observations: measurements
of luminosity distances of Type Ia Supernovae (SN Ia) [1, 2, 3], measurements
of cosmic microwave background (CMB) anisotropies in temperature and
polarization [4], measurements of the baryon acoustic oscillations (BAO) in
galaxy and cluster distribution [5, 6], cosmic shear measurements of the CMB
[7, 8, 9] and of the galaxy distribution [10, 11], and the predicted abundance
of light elements [12]. While the $\Lambda$CDM model provides an accurate
description to most of cosmological observations, it relies on a number of
assumptions and unknown ingredients such dark matter, dark energy, and a
suitable mechanism to produce its initial condition.

In addition to the interest in testing at which extent the validity of the
$\Lambda$CDM model holds with better and more data, the theoretical search for
extended models [13, 14, 15, 16, 17, 18, 19, 20, 21] has been fueled by the
persisting cosmological tensions or rather intriguing inconsistencies between
different measurements under the framework of the minimal $\Lambda$CDM model;
see Refs. [22, 23, 24, 25, 26, 27, 28, 29] for reviews on the topic.

Among the many proposed models, there are still difficulties in finding a
candidate able to completely solve the discrepancy between the value of the
Hubble parameter inferred in $\Lambda$CDM using CMB data from the Planck DR3,
i.e. $H_{0}=(67.36\pm 0.54)\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$ at 68%
confidence level (CL) [30], with the measurement from the SH0ES team [31]
obtained with cosmic distance ladder calibration of SN Ia from the revised
Pantheon+ compilation [3], i.e. $H_{0}=(73.0\pm 1.0)\,{\rm km}\,{\rm
s}^{-1}\,{\rm Mpc}^{-1}$ at 68% CL, once all cosmological data are combined.

It is even more difficult to reconcile the value of the Hubble parameter
together with the persistent but less significant tension between Planck and
galaxy shear experiments, quantified through the value of
$S_{8}\equiv\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$, see [27]. Adopting a flat
$\Lambda$CDM model, cosmic shear analysis of the fourth data release of the
Kilo-Degree Survey (KiDS-1000) reported $S_{8}=0.759^{+0.024}_{-0.021}$ [10]
and $S_{8}=0.776\pm 0.017$ from Dark Energy Survey (DES) Year 3 (Y3)
combination of three large-scale structures (LSS) two-point correlation
functions ($3\times 2$ pt) [11], while the value measured by Planck
corresponds to $S_{8}=0.832\pm 0.013$ [30]. Indeed, while minimally and
nonminimally coupled scalar field have been extensively studied as possible
solutions to the Hubble tension, they usually lead to a higher value of the
Hubble constant together with a larger growth of structures on small scales,
i.e. a higher value of $\sigma_{8}$, see Refs. [13, 32, 14, 33, 34, 35, 36].
However, they generally predict a value of $S_{8}$ compatible to the one
obtained in $\Lambda$CDM avoiding to exacerbate the tension on the growth of
structure amplitude since the larger $\sigma_{8}$ is compensated by a larger
value of $H_{0}$ and a lower value of $\Omega_{\rm m}$ [37, 20, 38, 39].

One possibility is to extend the dynamics of the scalar field to behave
differently at early- and late-time in order to solve both tensions at the
same time. The possibility to have models with phenomenology in both the early
and late universe has been tried in the context of modified gravity [20, 38],
early dark sector [40], and combining modified gravity or early dark energy to
extended neutrino physics [19, 41, 42].

In this paper, we study modified gravity models with a nonminimally coupled
scalar field with negative kinetic energy, so-called phantom field. Note that
this is not strictly related to the the phantom dark energy models for which
the dark energy (DE) equation of state can cross the phantom divide line
$w_{\rm DE}=-1$. Moreover, scalar-tensor models can be realized with no
necessity to introduce a ghost field [43] avoiding the problems in ghost
phantom DE [44] to be plagued by classical and quantum instabilities [45].
Such a the non-canonical kinetic energy term can occur in supergravity models
[46] and in higher derivative theories of gravity [47]. We show how a
nonminimally coupled scalar field with the negative sign of its kinetic term
(phantom branch) behaves differently compared to the case with standard
kinetic term (standard branch) and we derive the constraints on these models
combining the information from Planck 2018 DR3 CMB temperature, polarization
and lensing, together with a compilation of BAO measurements from the releases
DR7 and DR12 of the Baryon Oscillation Spectroscopic Survey (BOSS) and
Ly$\alpha$ measurements from the extended Baryon Oscillation Spectroscopic
Survey (eBOSS), and uncalibrated SN Ia from the Pantheon sample.

The paper is organized as follows. After this introduction, we describe the
implementation of the various basic quantities in the context of scalar-tensor
theories in Section 2. In Section 3, we describe the datasets and prior
considered and we discuss our results for the three models studies: induced
gravity, non-minimal coupling, and early modified gravity. In Section 4 we
draw our conclusions. In Appendix A, we collect the tables with the
constraints on all the cosmological parameters obtained with our MCMC
analysis. Background equations, linear perturbations, and initial conditions
for background and cosmological fluctuations are collected in Appendices
B-C-D. In Appendix E, we present a comparison of the results by using CMB data
plus different combinations of LSS measurements.

##  2 Theory and cosmological background dynamics

We study the action for the scalar-tensor theory in Jordan frame [48] which is
given by

${\cal
S}=\int{\mathrm{d}}^{4}x\sqrt{|g|}\left[\frac{F(\sigma)R}{2}-\frac{Z(\sigma)}{2}(\partial\sigma)^{2}-V(\sigma)+{\cal
L}_{m}\right]$ (2.1)

where $|g|$ is the absolute value of the determinant of the metric
$g_{\mu\nu}$, $\sigma$ is the scalar field, $F(\sigma)$ is the non-minimal
coupling function, $R$ is the Ricci scalar, $V(\sigma)$ is the potential for
$\sigma$, and ${\cal L}_{m}$ the Lagrangian density of matter minimally
coupled to the metric (without introducing any direct coupling between the
scalar field and the matter content we guarantee that the weak equivalence
principle is exactly satisfied). The function $Z(\sigma)$ in front of the
kinetic term can be set to $\pm 1$ by a redefinition of the scalar field.

In this model, the effective gravitational constant $G_{\rm eff}$ for the
attraction between two test masses [having the same physical meaning as the
Newton gravitational constant in general relativity (GR)] is given by

$G_{\rm eff}=\frac{1}{8\pi
F}\frac{ZF+2F^{2}_{\sigma}}{ZF+\frac{3}{2}F^{2}_{\sigma}}$ (2.2)

on all scales for which the scalar field is effectively massless [43], i.e.
$V_{\sigma}\simeq 0$ and $V_{\sigma\sigma}\simeq 0$.

The current values of the time derivative and field derivative of coupling $F$
in these theories - assuming a homogeneous evolution of the scalar field for
all the scales - are strongly constrained by Solar System tests of post-
Newtonian parameters (for these quantities, we drop here the subscript 0)

$\displaystyle\gamma_{\rm PN}=1-\frac{F_{\sigma}^{2}}{ZF+2F^{2}_{\sigma}}$
(2.3) $\displaystyle\beta_{\rm
PN}=1+\frac{1}{4}\frac{FF_{\sigma}}{2ZF+3F^{2}_{\sigma}}\frac{\mathrm{d}\gamma_{\rm
PN}}{\mathrm{d}\sigma}$ (2.4)

as well as the time variation of the effective cosmological constant. Current
constraints [49, 50, 51, 52] correspond to

$\displaystyle\gamma_{\rm PN}-1=\left(2.1\pm 2.3\right)\cdot 10^{-5}$ (2.5)
$\displaystyle\beta_{\rm PN}-1=\left(-4.1\pm 7.8\right)\cdot 10^{-4}$ (2.6)
$\displaystyle\dot{G}/G=\left(7.1\pm 7.6\right)\cdot 10^{-14}\,{\rm
yr}^{-1}\,.$ (2.7)

On cosmological scales, post-Newtonian parameters are weakly constrained from
current cosmological data, see Ref. [19], with the perspective to reach the
Solar System accuracy with the combination of future cosmological surveys [53,
54, 38].

There are essentially two stability conditions which impact on these scalar-
tensor theories. The condition

$G_{\rm eff}>0$ (2.8)

is one of the stability conditions of this theory meaning that the graviton is
not a ghost. Moreover, we have the inequality

$\frac{ZF}{F_{\sigma}^{2}}>-\frac{3}{2}$ (2.9)

requiring the positivity of the kinetic energy of the scalar field in the
Einstein frame [55]. Eqs. (2.8)-(2.9) reduce to $ZFF_{\sigma}^{-2}>0$ for
$Z=+1$ or $-3/2<ZFF_{\sigma}^{-2}<0$ for $Z=-1$, and to $F>0$. This condition
can be mapped to a range of allowed parameter space for the parameters
modelling $F(\sigma)$. However, we will consider a larger parameter space in
the following analysis testing the models also for parameters violating the
stability conditions in an agnostic way.

##  3 Constraints and results

In this section, we present our constraints on the cosmological parameters of
the models studied. In particular, we study the nonminimally coupling
$F=\xi\sigma^{2}$, i.e. induced gravity (IG) [56, 57], and $F=N_{\rm
Pl}^{2}+\xi\sigma^{2}$ (hereafter NMC) both with a phantom scalar field, i.e.
$Z=-1$; in both cases we consider $V(\sigma)=\lambda F^{2}(\sigma)/4$ which
yields to an effectively massless dynamic [58, 59]. We study also the early
modified gravity (EMG) model proposed in Ref. [20] extended to $Z=-1$. This
last case is given by $F=M_{\rm Pl}^{2}+\xi\sigma^{2}$ and
$V=\Lambda+\lambda\sigma^{4}/4$ with a negative amplitude $\lambda$ of the
self-interaction term in order to produce the peculiar evolution of the scalar
field damped into coherent oscillations within the phantom branch. We perform
a Markov-chain Monte Carlo (MCMC) analysis using a modified version of the
CLASSig code [13], based on the Einstein-Boltzmann code
CLASS111https://github.com/lesgourg/class_public [60, 61], interfaced to the
publicly sampling code
MontePython222https://github.com/brinckmann/montepython_public [62, 63]. The
datasets used in this work include

  * •

P18 refers to the CMB temperature, polarization, and lensing from Planck DR3
[64, 65].

  * •

FS refers to the combination of pre-reconstructed full-shape monopole and
quadrupole galaxy power spectra for three different sky-cuts CMASS NGC, CMASS
SGC and LOWZ NGC [66] based on the publicly available code
PyBird333https://github.com/pierrexyz/pybird/[67].

  * •

BAO refers to the post-reconstruction measurements from BOSS DR12 [68],
low-$z$ BAO measurements from SDSS DR7 6dF and MGS [69, 70], Ly$\alpha$ BAO
measurements from eBOSS, and combination of those [71, 72, 73].

  * •

SN refers to the Pantheon catalogue of high-redshift supernovae, spanning the
redshift range $0.01<z<2.3$ [74]444https://github.com/dscolnic/Pantheon.

  * •

Additional constraints of a Gaussian prior on the density of baryons
(hereafter BBN) motivated from Big Bang nucleosynthesis (BBN) constraints
corresponding to $\omega_{\rm b}=0.02235\pm 0.0005$ [12], used in combination
to FS and SN for the CMB-independent analysis.

  * •

Additional constraints that include a Gaussian prior on the Hubble constant
[hereafter $p(H_{0})$], $H_{0}=(73.04\pm 1.04)\,{\rm km}\,{\rm s}^{-1}\,{\rm
Mpc}^{-1}$ at 68% CL, from Ref. [31].

We vary 6 standard parameters, i.e. $\omega_{\rm b}$, $\omega_{\rm c}$,
$H_{0}$, $\tau$, $\ln\left(10^{10}A_{\rm s}\right)$, $n_{\rm s}$, and the
modified gravity parameters. We assume 2 massless neutrino with $N_{\rm
ur}=2.0328$, and a massive one with fixed minimum mass $m_{\nu}=0.06\,{\rm
eV}$. We fix the primordial ^4He mass fraction $Y_{\rm p}$ according to the
prediction from PArthENoPE [75, 76], by taking into account the relation with
the baryon fraction $\omega_{\rm b}$ and the varying gravitational constant
which enters in the Friedman equation during nucleosynthesis.

Following the minimization method of Ref. [27], we report for each combination
of datasets the $\Delta\chi^{2}$ values calculated with respect to the
$\Lambda$CDM model where negative values correspond to a better fit of the
dataset.

###  3.1 Phantom induced gravity

For IG (or equivalently extended Jordan-Brans-Dicke), with coupling
$F(\sigma)=\xi\sigma^{2}$, we sample on the quantity $\zeta_{\rm
IG}\equiv\ln\left(1+4\xi\right)$ which corresponds to a linear prior on the
coupling to the parameter $\xi$ for $\xi\ll 1$. Here we impose the following
boundary condition on the current value of the effective gravitational
constant

$G_{\rm eff}(z=0)=G$ (3.1)

which fixes the final value of the scalar field.

Scalar-tensor theories of gravity such as extended Jordan-Brans-Dicke models,
lead to a modification of the Hubble parameter $H_{0}$ due to the time
evolution of $\sigma$ and due to the redshift evolution of the gravitational
strength. Indeed, a variation of the strength of gravity can be connected to a
change of the expansion rate of the universe as

$\frac{H(\xi\neq 0)}{H(\xi=0)}\approx\frac{M_{\rm Pl}^{2}}{F(\sigma)}\,.$
(3.2)

For a fixed matter content, reducing the Planck mass $F(\sigma)<M_{\rm
Pl}^{2}$ with respect to the GR prediction increases the expansion rate at a
given time and consequently reduces the comoving sound horizon at
recombination

$r_{s}=\int_{z_{*}}^{\infty}\mathrm{d}z^{\prime}\frac{c_{\rm
s}(z^{\prime})}{H(z^{\prime})}$ (3.3)

where $z_{*}$ is the redshift parameter at recombination and $c_{\rm s}$ is
the speed of sound in the photon-baryon fluid. We show in Fig. 1 that the
coupling function increases in the branch with standard kinetic term (solid
lines) while decreasing in the phantom branch (dashed lines). This different
behaviour is connected with a different late-time evolution of the Hubble
parameter (when the scalar field starts to evolve driven by the non-
relativistic matter) which is larger than the $\Lambda$CDM case in the
standard branch and lower in the phantom branch.

Figure 1: Time evolution of the coupling to the Ricci scalar $F=\xi\sigma^{2}$
(left panel) and of the Hubble parameter (right panel) for different values of
the coupling parameter $\xi$ in the standard branch (solid lines) and in the
phantom one (dashed lines) for IG ($F=\xi\sigma^{2},\ V=\lambda F^{2}/4$).

This effect induces also a modification on the comoving angular diameter
distance

$D_{\rm M}(z)=\int_{0}^{z}\frac{\mathrm{d}z^{\prime}}{H(z^{\prime})}$ (3.4)

and does not cancel out on the angular size of horizon at the last-scattering
surface $\theta_{*}$

$\theta_{*}=\frac{r_{s}}{D_{\rm M}(z_{*})}$ (3.5)

driving a shift on the acoustic peaks of the CMB connected to the evolution of
the coupling $F$ [77, 19, 78]. In Fig. 2, we show the shift of the acoustic
peaks of the CMB temperature anisotropies angular power spectrum imprinted by
the evolving effective Planck mass. The peaks move to the right in the
positive branch and to the left in the phantom one. Indeed, in order to
compensate this shift (keeping nearly unchanged the value of the CDM density
parameter) the two branches go in the direction of a larger or a smaller value
of $\Omega_{\rm m}$ once the CMB data are included in the analysis, see Fig.
3.

Figure 2: CMB temperature anisotropies power spectrum (left panel) and
relative differences with respect to the $\Lambda$CDM case (right panel) for
different values of the coupling parameter $\xi$ in the standard branch (solid
lines) and in the phantom one (dashed lines) for IG ($F=\xi\sigma^{2},\
V=\lambda F^{2}/4$).

Therefore it is possible to break the degeneracy at background level between
the scalar field $\sigma$ and the density parameters by combining early- and
late-time probes.

Figure 3: Marginalized joint 68% and 95% CL regions 2D parameter space using
the CMB alone data for IG ($F=\xi\sigma^{2},\ V=\lambda F^{2}/4$) with $Z=1$
(orange) and for IG with $Z=-1$ (blue). Figure 4: Marginalized joint 68% and
95% CL regions 2D parameter space using the combination CMB+BAO+SN (orange)
and CMB+BAO+SN+$p(H_{0})$ (blue) for IG ($F=\xi\sigma^{2},\ V=\lambda
F^{2}/4$) in the phantom branch ($Z=-1$).

In Fig. 3, we compare the CMB-only constraints for IG in the phantom branch to
the standard case. We see that the degeneracy direction in the $\xi$-$H_{0}$
plane changes orientation going from one case to the other according to Eq.
(3.2). It turns out that the phantom branch allows much larger values of the
coupling $\xi$ and predicts a lower value for the Hubble constant without any
prospect to reduce the $H_{0}$ tension. It is interesting to note that the
extension of our study to the phantom case strengthen the correspondence
between the kinetic term and the spatial curvature: the standard (phantom)
kinetic term shifts the position to the right (left) as a negative (positive)
spatial curvature.

Figure 5: Time evolution of the amplitude of matter perturbation within
spheres of radius $8\,h^{-1}\,{\rm Mpc}$ (left panel) and relative differences
of the linear matter power spectrum at $z=0$ with respect to $\Lambda$CDM
(right panel) for different values of the coupling parameter $\xi$ in the
standard branch (solid lines) and in the phantom one (dashed lines) for IG
($F=\xi\sigma^{2},\ V=\lambda F^{2}/4$).

Finally, it is interesting to note that the matter density root mean square
fluctuations $\sigma_{8}$ goes toward lower values in the phantom branch
compared to both the standard branch and the $\Lambda$CDM model predictions,
see Fig. 5. This behaviour can be understood studying the late-time solution
of the perturbation equation for the matter density contrast in the linear
regime, on sub-horizon scales

$\delta_{\rm
m}^{\prime\prime}+\left(\frac{3}{a}+\frac{H^{\prime}}{H}\right)\delta_{\rm
m}^{\prime}-\frac{G_{\rm eff}}{2GH^{2}}\frac{\rho_{\rm m}}{a^{2}}\delta_{\rm
m}\simeq 0$ (3.6)

where primes are derivatives with respect to the scale factor $a$. By
rewriting the Friedmann equations (B.4)-(B.5)

$\displaystyle H^{2}=\frac{\rho+V}{3F(1+f)}$ (3.7)
$\displaystyle\frac{H^{\prime}}{H}=-\frac{3}{2a}-\frac{3}{a}w-\frac{F^{\prime}}{2F}+\frac{f^{\prime}}{2}$
(3.8)

introducing the quantity

$f=-a\frac{F_{\sigma}}{F}\sigma^{\prime}-\frac{a^{2}Z}{6F}\sigma^{\prime 2}$
(3.9)

and where we used $\rho^{\prime}/\rho=-3(1+w)/a$, we can write Eq. (3.6) as

$\delta_{\rm
m}^{\prime\prime}+\left[\frac{3}{2a}(1-w)-\frac{F^{\prime}}{2F}+\frac{f^{\prime}}{2}\right]\delta_{\rm
m}^{\prime}-\frac{3}{2a^{2}}\frac{2ZF+4F_{\sigma}^{2}}{2ZF+3F_{\sigma}^{2}}(1-f)\frac{\rho_{\rm
m}}{\rho+V}\delta_{\rm m}\simeq 0\,.$ (3.10)

During the matter-dominated era, the scalar field evolves as $\sigma\sim
a^{2Z\xi}$ [59] leading for IG to

$f\sim-\frac{10Z\xi}{3}$ (3.11)

and consequently

$\delta_{\rm
m}^{\prime\prime}+\frac{3}{2a}\left(1-\frac{4Z\xi}{3}\right)\delta_{\rm
m}^{\prime}-\frac{3}{2a^{2}}\left(1+\frac{16Z\xi}{3}\right)\delta_{\rm
m}\simeq 0\,.$ (3.12)

In the weak coupling regime for $\xi\ll 1$, which turns out to be the range
allowed from observations, the leading-order growing solution of Eq. (3.12)
goes as $\delta_{\rm m}\sim a^{1+4Z\xi}$ showing a slower (faster) growth of
structures compared to the $\Lambda$CDM case for $Z<0$ ($Z>0$) during the
matter dominated era.

We show the results for most of the combination of datasets on Fig. 6 (in
Appendix E, we show a comparison between P18+BAO and P18+BAO+FS). The
marginalized upper bound on the coupling parameter $\xi$ at 95% CL corresponds
to $<0.0024$ for FS+SN, $<0.0018$ for P18, $<0.00046$ for P18+BAO, and
$<0.00040$ for P18+BAO+SN; see Tab. 1 for the constraints on all the
parameters.

Figure 6: Marginalized joint 68% and 95% CL regions 2D parameter space using
the CMB-independent combination FS+SN (red), P18 (green), the combination
P18+BAO (orange), and the combination P18+BAO+SN (blue) for IG
($F=\xi\sigma^{2},\ V=\lambda F^{2}/4$) in the phantom branch ($Z=-1$).

In Fig. 6, we see the larger marginalized uncertainties for the analysis
without CMB information, i.e. combining FS with SN and a Gaussian prior on
$\omega_{\rm b}$ motivated from BBN, and the analysis with CMB alone. In these
cases larger value of $\xi$ can be accommodated by changes in the density
parameters and in the scalar spectral index.

The marginalized means and uncertainties for the Hubble constant
$H_{0}\,\left[{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}\right]$ at 68% CL
correspond to $67.4\pm 1.8$ for FS+SN, $63.6^{+2.7}_{-1.9}$ for P18,
$67.17^{+0.64}_{-0.50}$ for P18+BAO, and $67.29^{+0.60}_{-0.47}$ for
P18+BAO+SN; all of them are lower than the corresponding results we found in
the standard branch, see Refs. [19, 79]. The upper bound on $\xi$ becomes much
tighter, i.e. $\xi<0.00016$ at 95% CL, once we add a prior on $H_{0}$ in order
to reproduce a larger value of the Hubble parameter, i.e. $68.34\pm 0.41$ at
68% CL, see Fig. 4 and Tab. 4.

The marginalized constraint on the present value of $\sigma_{8}$ at 68% CL
corresponds to $0.717\pm 0.049$ for FS+SN, $0.784^{+0.021}_{-0.015}$ for P18,
$0.799^{+0.010}_{-0.009}$ for P18+FS, and $0.8059\pm 0.0058$ for P18+FS+SN.
However, the combination $S_{8}\equiv\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$,
commonly used to quantify the tension between Planck and weak lensing of
galaxies measurements, moves to the wrong direction. In order not to spoil the
fit to CMB and galaxy measurements, an increase of the matter energy density
is needed to compensate for the shifted position of the CMB acoustic peaks and
of the BAO. Indeed, we find for $S_{8}$ $0.744\pm 0.050$ for FS+SN,
$0.850^{+0.016}_{-0.019}$ for P18, $0.831^{+0.011}_{-0.012}$ for P18+FS, and
$0.825\pm 0.012$ for P18+FS+SN at 68% CL, resulting to be larger than in the
standard branch as shown in Fig. 3.

The constraints found are at odds with the parameter space free from ghost
which corresponds for IG to $\xi>1/6$ according to Eqs. (2.8)-(2.9). Note
that, this condition for $\xi$ can be relaxed if one considers a more general
Lagrangian with respect to the one introduced in Eq. (2.1). Higher order terms
in the kinetic energy $X\equiv-\partial_{\mu}\sigma\partial^{\mu}\sigma/2$
[80] appear in low-energy effective string theory [81] or in tachyon
condensation [82]. In a more general Lagrangian containing also a Galileon
term $G_{3}$ [83, 84] as in ${\cal
L}=G_{4}(\sigma)R+G_{2}(\sigma,X)+G_{3}(\sigma,X)\square\sigma$, the
conditions for the avoidance of such instabilities are

$q_{s}\equiv\,4G_{4}\big{\\{}G_{2X}+2G_{3\sigma}+\dot{\sigma}\big{[}(G_{2XX}+G_{3X\sigma})\dot{\sigma}-6G_{3X}H\big{]}\big{\\}}+3(2G_{4\sigma}+G_{3X}\dot{\sigma}^{2})^{2}>0$
(3.13)
$\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}$
(3.14)

which reduce to Eqs. (2.8)-(2.9) when $G_{3}=0$; but, in general, depending on
the functional form of the cubic interaction term and its magnitude, can allow
for $\xi<1/6$ while maintaining the theory free of ghost and Laplacian
instabilities [85, 86, 77, 87].

###  3.2 Phantom non-minimal coupling

For NMC+ (NMC-) [14], with coupling $F(\sigma)=N_{\rm Pl}^{2}+\xi\sigma^{2}$
and coupling $V(\sigma)=\lambda F^{2}(\sigma)/4$, we sample on the
dimensionless parameter $\Delta\tilde{N}_{\rm Pl}\equiv N_{\rm Pl}/M_{\rm
Pl}-1$ and $\xi$.

We show the results for all the combinations of datasets on Fig. 7 for NMC+
and for NMC- in Fig. 8. Marginalized constraints on cosmological parameters
are consistent to the results obtained for IG for each combination of
datasets. The marginalized limits on the coupling parameters for NMC+ (NMC-)
at 95% CL correspond to $\xi<0.0015$ ($>-0.039$) and $N_{\rm Pl}>0.91$
($<1.18$) for P18+BAO, and to $\xi<0.0019$ ($>-0.027$) and $N_{\rm Pl}>0.83$
($<1.21$) for P18+BAO+SN; see Tabs. 2 and 3 for the constraints on all the
parameters. As observed in Refs. [14, 19], there is a strong degeneracy
between the coupling parameters $N_{\rm Pl}$ and $\xi$ for the form of non-
minimal coupling $F(\sigma)=N_{\rm Pl}^{2}+\xi\sigma^{2}$. Since cosmological
observables are affected by contributions ${\cal O}\left(\xi\sigma^{2}/N_{\rm
Pl}^{2}\right)$, it is possible to compensate the effects due to a large value
of $|\xi|$ increasing $|\tilde{N}_{\rm Pl}-1|$ and vice versa.

Figure 7: Marginalized joint 68% and 95% CL regions 2D parameter space using
using the P18 (green), the combination P18+BAO (orange), and the combination
P18+BAO+SN (blue) for NMC+ ($F=N_{\rm Pl}^{2}+\xi\sigma^{2},\ V=\lambda
F^{2}/4$) in the phantom branch ($Z=-1$). Figure 8: Marginalized joint 68% and
95% CL regions 2D parameter space using using the P18 (red), the combination
P18+BAO (green), the combination P18+BAO+SN (orange), and the combination
P18+BAO+SN+$p(H_{0})$ (blue) for NMC- ($F=N_{\rm Pl}^{2}+\xi\sigma^{2},\
V=\lambda F^{2}/4$) in the phantom branch ($Z=-1$).

In this case Eqs. (2.8)-(2.9) reduce to

$\left(\frac{N_{\rm Pl}}{\xi\sigma}\right)^{2}+\frac{1}{\xi}<6\,.$ (3.15)

Also in this case, a large portion of the allowed parameter space is at odds
with Eq. (3.15) despite the larger number of degrees of freedom.

###  3.3 Phantom early modified gravity

For EMG [20], with coupling $F(\sigma)=M_{\rm Pl}^{2}+\xi\sigma^{2}$ and
potential $V(\sigma)=\Lambda+\lambda\sigma^{4}/4$, we sample on the quantity
$\xi$ and $V_{0}$ where $\lambda\equiv-10^{2V_{0}}/M_{\rm Pl}^{4}$. In this
case the scalar field decays into the minimum of the potential, i.e.
$\sigma=0$, and we do not have to impose the boundary condition (3.1) being
automatically satisfied for each initial value of scalar field. This leads to
a third free parameter which we identify with the initial value of the scalar
field $\sigma_{\rm ini}$ $\left[M_{\rm Pl}\right]$.

In Fig. 9, we compare the background evolution and spectra of EMG in the
standard branch ($Z=1$) to the phantom branch case ($Z=-1$). The evolution of
the scalar field $\sigma$ is very similar in the two cases. Starting with the
scalar field at rest in the radiation era, it starts to grow around the
recombination driven by the coupling to the non-relativistic matter component
and it is subsequently driven into damped coherent oscillations around the
minimum of the quartic potential. The different evolution of the coupling
function, which increase in the branch with standard kinetic term while
decreasing in the phantom branch before the scalar field starts to decay, due
to the different sign of the coupling parameter $\xi$ induces different
effects on the spectra. The acoustic peaks of the CMB temperature anisotropies
angular power spectrum are shifted to right in the standard branch ($Z=1$)
when the scalar field starts to move before recombination ($V_{0}=2$) and in
the phantom branch ($Z=-1$) if the scalar field decays after recombination
($V_{0}=-1$), vice versa they shift to the left with respect the $\Lambda$CDM
case. The situation is different on the linear matter power spectrum were we
observe a suppression of power in the standard branch ($Z=1$) despite the
value of $V_{0}$ and an increase of power at small scales in the phantom
branch ($Z=-1$). This highlights the importance of the combination of
combining early- and late-time probes in order to break the degeneracy between
the extra parameters of the model and also to discriminate between the two
different branches.

Figure 9: Time evolution of the scalar field $\sigma$ (upper left panel) and
of the coupling to the Ricci scalar $F(\sigma)$ (upper right panel). Relative
differences of the CMB temperature anisotropies power spectrum with respect to
the $\Lambda$CDM case (bottom left panel) and of the linear matter power
spectrum at $z=0$ (bottom right panel). Different lines correspond to
different value of the amplitude of the effective potential $V_{0}$ for
$|\xi|=0.1$ in the standard branch (solid lines) and in the phantom one
(dashed lines) for EMG ($F=M_{\rm
Pl}^{2}+\xi\sigma^{2},\,V=\Lambda+\lambda\sigma^{4}/4$).

We show the results for the combinations of datasets P18+BAO+SN and P18
+BAO+SN+$p(H_{0})$ on Fig. 10. In this case we find that $\xi$ is not
constrained by data on the prior considered $[-0.1,\,0]$. For this reason, we
show constraints on the combination $\xi\sigma_{\rm ini}^{2}$ $[M_{\rm
Pl}^{2}]$ (connected to the additional contribution to the expansion rate
evolution (3.2) before recombination). The marginalized upper bound on the
coupling combination $\xi\sigma_{\rm ini}^{2}$ at 95% CL corresponds to
$>-0.0026$ for P18+BAO+SN and when we include the Gaussian prior on the Hubble
parameter we obtain at 95% CL $-0.006\pm 0.005$. Analogously, for the initial
value of the scalar field $\sigma_{\rm ini}$ $[M_{\rm Pl}]$ we find $<0.45$
for P18+BAO+SN and $0.35^{+0.17}_{-0.15}$ for P18+BAO+SN+$p(H_{0})$ both at
95% CL. Also $V_{0}$ is not well constrained. We find a 95% CL upper bound
only when we include the Gaussian prior on the Hubble parameter corresponding
to $V_{0}<0.81$.

The marginalized means and uncertainties for the Hubble constant
$H_{0}\,\left[{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}\right]$ at 68% CL
correspond to $68.44^{+0.62}_{-0.79}$ for P18+BAO+SN and
$70.18^{+0.59}_{-0.68}$ for P18+BAO+SN +$p(H_{0})$. The marginalized
constraints on $S_{8}$ correspond to $S_{8}=0.827\pm 0.011$ for P18+BAO+SN and
$S_{8}=0.822\pm 0.010$ for P18+BAO+SN+$p(H_{0})$. See Tab. 5 for the
constraints on all the parameters.

Figure 10: Marginalized joint 68% and 95% CL regions 2D parameter space using
P18+BAO+SN (orange) and the combination P18+BAO+SN+$p(H_{0})$ (blue) for EMG
($F=M_{\rm Pl}^{2}+\xi\sigma^{2},\ V=\Lambda+\lambda\sigma^{4}/4$) in the
phantom branch ($Z=-1$).

##  4 Conclusions

We have studied the dynamics and inferred the cosmological constraints for
modified gravity models with a nonminimally coupled scalar with a kinetic term
which can also have a negative sign. For stable models with an effectively
massless scalar field $\sigma$, like IG and NMC, the change of sign in front
of the kinetic term of the scalar field modifies the evolution of the scalar
field which is at rest during the radiation-dominated epoch and evolves like
$\sigma\sim a^{2Z\xi}$ during the matter-dominated era.

We have shown the effect of the sign of the kinetic term on cosmological
observables. We have computed the marginalized constraints for different
combination of cosmological datasets by allowing the coupling to the Ricci
scalar and the rest of cosmology (standard cosmological parameters and
nuisance ones) to vary. Combining Planck 2018 DR3 measurements with BAO from
BOSS and eBOSS, and uncalibrated SN I$\alpha$ from the Pantheon sample we
constrain the coupling parameters at 95% CL to $\xi<0.00040$ for
$F(\sigma)=\xi\sigma^{2}$ and for $F(\sigma)=N_{\rm Pl}^{2}+\xi\sigma^{2}$ to
$\xi<0.0019$ ($>-0.027$) and $N_{\rm Pl}>0.83$ ($<1.21$).

Nonminimally coupled scalar-tensor theories with early-time deviation from GR
predictions usually lead to higher values of the Hubble parameter $H_{0}$, a
lower value of the matter density parameter $\Omega_{\rm m}$, and a larger
value of the $\sigma_{8}$ [78, 77, 19]. In their phantom construction, the
modified evolution of the scalar field, connected to a different time
evolution of the effective gravitational constant, inverts the degeneracy
between these parameters and the coupling ones. Indeed, we find a lower values
of both $\sigma_{8}$ and $H_{0}$ compared to the branch with standard kinetic
term.

We have also studied the phantom version of the EMG model introduced in Ref.
[20]. While the evolution of the scalar field is very similar, with the
quartic potential leading the scalar field to decay into damped coherent
oscillations, different signatures appear on the cosmological observables.
Compared to the $\Lambda$CDM model, the CMB acoustic peaks of the CMB are
shifted to right in the standard branch ($Z=1$) when the scalar field starts
to move before recombination ($V_{0}=2$) and in the phantom branch ($Z=-1$) if
the scalar field decays after recombination ($V_{0}=-1$), vice versa they
shift to the left. Matter perturbations on sub-horizon scales are suppressed
in the standard branch ($Z=1$) and enhanced in the phantom branch ($Z=-1$)
despite the value of amplitude of the self-interaction term parameterized by
$V_{0}$.

The allowed parameter space for the coupling parameters by our analysis is at
odds with the parameter space free from ghosts and Laplacian instabilities. It
would be interesting to understand if instead there are healthy scalar-tensor
theories which retain the possibility to alleviate the current tensions
between different cosmological observations.

## Acknowledgments

MB and FF acknowledge financial support from the INFN InDark initiative and
from the COSMOS network (www.cosmosnet.it) through the ASI (Italian Space
Agency) Grants 2016-24-H.0 and 2016-24-H.1-2018, as well as 2020-9-HH.0
(participation in LiteBIRD phase A). This work has made use of computational
resources of INAF OAS Bologna and of the CNAF HPC cluster in Bologna.

##  Appendix A Tables

| P18 | P18 + BAO | P18 + BAO + SN
---|---|---|---
$\omega_{\rm b}$ | $0.02223\pm 0.00017$ | $0.02244\pm 0.00013$ | $0.02245\pm 0.00013$
$\omega_{\rm c}$ | $0.1204\pm 0.0012$ | $0.11896\pm 0.00099$ | $0.11882\pm 0.00096$
$H_{0}$ [km s-1Mpc-1] | $63.6^{+2.7}_{-1.9}$ | $67.17^{+0.65}_{-0.48}$ | $67.29^{+0.60}_{-0.47}$
$\tau$ | $0.0523\pm 0.0071$ | $0.0584^{+0.0070}_{-0.0083}$ | $0.0584^{+0.0068}_{-0.0076}$
$\ln\left(10^{10}A_{\rm s}\right)$ | $3.037\pm 0.015$ | $3.051^{+0.014}_{-0.016}$ | $3.051\pm 0.014$
$n_{\rm s}$ | $0.9574^{+0.0067}_{-0.0057}$ | $0.9668\pm 0.0038$ | $0.9671\pm 0.0036$
$\xi$ |  $<0.0018$ (95% CL) |  $<0.00046$ (95% CL) |  $<0.00040$ (95% CL)
$\gamma_{\rm PN}$ |  $>0.9928$ (95% CL) |  $>0.9982$ (95% CL) |  $>0.9984$ (95% CL)
$\delta G_{\mathrm{N}}/G_{\mathrm{N}}$ (z=0) |  $<0.057$ (95% CL) |  $<0.014$ (95% CL) |  $<0.012$ (95% CL)
$10^{13}\dot{G}_{\mathrm{N}}/G_{\rm N}$ (z=0) [yr-1] |  $<2.43$ (95% CL) |  $<0.57$ (95% CL) |  $<0.50$ (95% CL)
$G_{\mathrm{N}}/G$ (z=0) | $0.9982^{+0.0013}_{-0.00074}$ | $0.99965^{+0.00033}_{-0.00013}$ | $0.99968^{+0.00030}_{-0.00012}$
$\Omega_{\rm m}$ | $0.354^{+0.020}_{-0.032}$ | $0.3135^{+0.0059}_{-0.0068}$ | $0.3120^{+0.0056}_{-0.0065}$
$\sigma_{8}$ | $0.784^{+0.021}_{-0.015}$ | $0.8053^{+0.0083}_{-0.0069}$ | $0.8053^{+0.0078}_{-0.0066}$
$S_{8}$ | $0.850^{+0.016}_{-0.019}$ | $0.823\pm 0.011$ | $0.821\pm 0.010$
$r_{s}$ [Mpc] | $148.89^{+0.8}_{-1.4}$ | $147.62^{+0.31}_{-0.50}$ | $147.62^{+0.29}_{-0.44}$
$\Delta\chi^{2}$ | $-2.8$ | $0$ | $-0.5$
Table 1:  Constraints on the main and derived parameters (at 68% CL if not otherwise stated) considering P18 in combination with BAO and BAO+SN for the IG model. | P18 | P18 + BAO | P18 + BAO + SN
---|---|---|---
$\omega_{\rm b}$ | $0.02224^{+0.00018}_{-0.00016}$ | $0.02246\pm 0.00013$ | $0.02246^{+0.00011}_{-0.00014}$
$\omega_{\rm c}$ | $0.1206\pm 0.0012$ | $0.1190^{+0.0014}_{-0.0011}$ | $0.1189\pm 0.0010$
$H_{0}$ [km s-1Mpc-1] | $64.1^{+2.6}_{-1.7}$ | $67.28\pm 0.59$ | $67.42\pm 0.52$
$\tau$ | $0.0514\pm 0.0081$ | $0.0590\pm 0.0052$ | $0.0583\pm 0.0071$
$\ln\left(10^{10}A_{\rm s}\right)$ | $3.037^{+0.015}_{-0.023}$ | $3.0517^{+0.0074}_{-0.015}$ | $3.051\pm 0.014$
$n_{\rm s}$ | $0.9580^{+0.0058}_{-0.0047}$ | $0.9673\pm 0.0042$ | $0.9674\pm 0.0039$
$\xi$ |  $<0.030$ (95% CL) |  $<0.015$ (95% CL) |  $<0.019$ (95% CL)
$N_{\rm Pl}$ [$M_{\rm Pl}$] | $-$ |  $>0.91$ (95% CL) |  $>0.83$ (95% CL)
$\gamma_{\rm PN}$ |  $>0.9941$ (95% CL) |  $>0.9986$ (95% CL) |  $>0.9987$ (95% CL)
$\beta_{\rm PN}$ |  $>0.999965$ (95% CL) |  $>0.999994$ (95% CL) |  $>0.999994$ (95% CL)
$\delta G_{\mathrm{N}}/G_{\mathrm{N}}$ (z=0) |  $<0.052$ (95% CL) |  $<0.011$ (95% CL) |  $<0.012$ (95% CL)
$10^{13}\dot{G}_{\mathrm{N}}/G_{\rm N}$ (z=0) [yr-1] |  $<1.94$ (95% CL) |  $<0.42$ (95% CL) |  $<0.40$ (95% CL)
$G_{\mathrm{N}}/G$ (z=0) | $1.00134^{+0.00063}_{-0.0011}$ | $1.000255^{+0.000097}_{-0.00024}$ | $1.000221^{+0.000084}_{-0.00023}$
$\Omega_{\rm m}$ | $0.349^{+0.017}_{-0.030}$ | $0.3121^{+0.0068}_{-0.0056}$ | $0.3110\pm 0.0061$
$\sigma_{8}$ | $0.788^{+0.021}_{-0.013}$ | $0.8069\pm 0.0069$ | $0.8064^{+0.0079}_{-0.0065}$
$S_{8}$ | $0.849^{+0.013}_{-0.019}$ | $0.823^{+0.014}_{-0.009}$ | $0.821^{+0.012}_{-0.011}$
$r_{s}$ [Mpc] | $148.56^{+0.90}_{-1.3}$ | $147.54^{+0.30}_{-0.48}$ | $147.52^{+0.27}_{-0.44}$
$\Delta\chi^{2}$ | $-1.5$ | $0$ | $-0.5$
Table 2:  Constraints on the main and derived parameters (at 68% CL if not otherwise stated) considering P18 in combination with BAO and BAO+SN for the NMC+ model. | P18 | P18 + BAO | P18 + BAO + SN
---|---|---|---
$\omega_{\rm b}$ | $0.02230\pm 0.00014$ | $0.02245\pm 0.00013$ | $0.02247\pm 0.00013$
$\omega_{\rm c}$ | $0.11982^{+0.00068}_{-0.0011}$ | $0.11891\pm 0.00094$ | $0.11875^{+0.00078}_{-0.0010}$
$H_{0}$ [km s-1Mpc-1] | $64.1^{+3.1}_{-2.1}$ | $67.26^{+0.59}_{-0.45}$ | $67.44^{+0.57}_{-0.45}$
$\tau$ | $0.0548^{+0.0072}_{-0.0059}$ | $0.0573^{+0.0061}_{-0.0074}$ | $0.0590\pm 0.0068$
$\ln\left(10^{10}A_{\rm s}\right)$ | $3.041^{+0.017}_{-0.013}$ | $3.049\pm 0.014$ | $3.052^{+0.014}_{-0.012}$
$n_{\rm s}$ | $0.9604^{+0.0067}_{-0.0045}$ | $0.9669^{+0.0043}_{-0.0035}$ | $0.9675\pm 0.0036$
$\xi$ |  $>-0.036$ (95% CL) |  $>-0.039$ (95% CL) |  $>-0.027$ (95% CL)
$N_{\rm Pl}$ [$M_{\rm Pl}$] |  $<1.13$ (95% CL) |  $<1.18$ (95% CL) |  $<1.21$ (95% CL)
$\gamma_{\rm PN}$ |  $>0.988$ (95% CL) |  $>0.998$ (95% CL) |  $>0.998$ (95% CL)
$\beta_{\rm PN}$ |  $<1.00018$ (95% CL) |  $<1.000022$ (95% CL) |  $<1.000017$ (95% CL)
$\delta G_{\mathrm{N}}/G_{\mathrm{N}}$ (z=0) |  $<0.060$ (95% CL) |  $<0.012$ (95% CL) |  $<0.010$ (95% CL)
$10^{13}\dot{G}_{\mathrm{N}}/G_{\rm N}$ (z=0) [yr-1] |  $<3.85$ (95% CL) |  $<0.62$ (95% CL) |  $<0.50$ (95% CL)
$G_{\mathrm{N}}/G$ (z=0) | $1.00224^{+0.00080}_{-0.0021}$ | $1.00037^{+0.00012}_{-0.00037}$ | $1.00030^{+0.00012}_{-0.00030}$
$\Omega_{\rm m}$ | $0.348^{+0.021}_{-0.033}$ | $0.3125^{+0.0052}_{-0.0065}$ | $0.3106^{+0.0050}_{-0.0068}$
$\sigma_{8}$ | $0.786^{+0.025}_{-0.011}$ | $0.8047\pm 0.0076$ | $0.8064\pm 0.0062$
$S_{8}$ | $0.844^{+0.011}_{-0.018}$ | $0.821\pm 0.011$ | $0.8204^{+0.0091}_{-0.012}$
$r_{s}$ [Mpc] | $148.91^{+0.77}_{-1.6}$ | $147.60^{+0.28}_{-0.43}$ | $147.57^{+0.30}_{-0.40}$
$\Delta\chi^{2}$ | $-2.8$ | $0$ | $-0.3$
Table 3:  Constraints on the main and derived parameters (at 68% CL if not otherwise stated) considering P18 in combination with BAO and BAO+SN for the NMC- model. | IG | NMC+ | NMC-
---|---|---|---
$\omega_{\rm b}$ | $0.02260^{+0.00012}_{-0.00014}$ | $0.02262\pm 0.00013$ | $0.02260^{+0.00014}_{-0.00012}$
$\omega_{\rm c}$ | $0.11747\pm 0.00086$ | $0.11752\pm 0.00086$ | $0.11753^{+0.00095}_{-0.00069}$
$H_{0}$ [km s-1Mpc-1] | $68.34\pm 0.41$ | $68.42^{+0.44}_{-0.36}$ | $68.38\pm 0.41$
$\tau$ | $0.0617^{+0.0067}_{-0.0085}$ | $0.0636^{+0.0071}_{-0.0081}$ | $0.0644^{+0.0068}_{-0.0090}$
$\ln\left(10^{10}A_{\rm s}\right)$ | $3.055^{+0.013}_{-0.017}$ | $3.059^{+0.014}_{-0.016}$ | $3.061^{+0.014}_{-0.017}$
$n_{\rm s}$ | $0.9711\pm 0.0036$ | $0.9716\pm 0.0037$ | $0.9712\pm 0.0035$
$\xi$ |  $<0.000075$ (95% CL) |  $<0.0096$ (95% CL) |  $>-0.022$ (95% CL)
$N_{\rm Pl}$ [$M_{\rm Pl}$] | $0$ |  $>0.82$ (95% CL) |  $<1.24$ (95% CL)
$\gamma_{\rm PN}$ |  $>0.9993$ (95% CL) |  $>0.9995$ (95% CL) |  $>0.9994$ (95% CL)
$\beta_{\rm PN}$ | $1$ |  $>0.999999$ (95% CL) |  $<1.000004$ (95% CL)
$\delta G_{\mathrm{N}}/G_{\mathrm{N}}$ (z=0) |  $<0.0056$ (95% CL) |  $<0.0041$ (95% CL) |  $<0.0042$ (95% CL)
$10^{13}\dot{G}_{\mathrm{N}}/G_{\rm N}$ (z=0) [yr-1] |  $<0.23$ (95% CL) |  $<0.16$ (95% CL) |  $<0.19$ (95% CL)
$G_{\mathrm{N}}/G$ (z=0) | $0.99987^{+0.00013}_{-0.000042}$ | $1.000080^{+0.000028}_{-0.000090}$ | $1.000102^{+0.000043}_{-0.00011}$
$\Omega_{\rm m}$ | $0.2999\pm 0.0050$ | $0.2994^{+0.0046}_{-0.0052}$ | $0.2997^{+0.0053}_{-0.0047}$
$\sigma_{8}$ | $0.8055^{+0.0059}_{-0.0069}$ | $0.8078\pm 0.0066$ | $0.8086\pm 0.0063$
$S_{8}$ | $0.805\pm 0.010$ | $0.807\pm 0.010$ | $0.8082\pm 0.0094$
$r_{s}$ [Mpc] | $147.63^{+0.23}_{-0.29}$ | $147.56\pm 0.25$ | $147.58^{+0.21}_{-0.24}$
Table 4:  Constraints on the main and derived parameters (at 68% CL if not otherwise stated) considering the combination with P18+BAO+SN+$p(H_{0})$ for IG, NMC+, and NMC-. | P18 + BAO + SN |  P18 + BAO + SN + $p(H_{0})$
---|---|---
$\omega_{\rm b}$ | $0.02246\pm 0.00014$ | $0.02255\pm 0.00014$
$\omega_{\rm c}$ | $0.1194\pm 0.0010$ | $0.11900\pm 0.00099$
$H_{0}$ [km s-1Mpc-1] | $68.44^{+0.62}_{-0.79}$ | $70.18^{+0.59}_{-0.68}$
$\tau$ | $0.0536\pm 0.0080$ | $0.0503^{+0.0085}_{-0.0073}$
$\ln\left(10^{10}A_{\rm s}\right)$ | $3.043\pm 0.016$ | $3.035^{+0.017}_{-0.015}$
$n_{\rm s}$ | $0.9671^{+0.0036}_{-0.0042}$ | $0.9687\pm 0.0038$
$\xi\sigma_{\rm ini}^{2}$ [$M_{\rm Pl}^{2}$] |  $>-0.0057$ (95% CL) | $-0.0062^{+0.0028}_{-0.0023}$
$V_{0}$ | $-$ |  $<0.81$ (95% CL)
$\sigma_{\rm ini}$ [$M_{\rm Pl}$] |  $<0.446$ (95% CL) | $0.348^{+0.062}_{-0.097}$
$\Omega_{\rm m}$ | $0.3028\pm 0.0068$ | $0.2875\pm 0.0056$
$\sigma_{8}$ | $0.823^{+0.010}_{-0.013}$ | $0.840\pm 0.011$
$S_{8}$ | $0.827\pm 0.011$ | $0.822\pm 0.010$
$r_{s}$ [Mpc] | $147.00\pm 0.40$ | $146.56\pm 0.46$
Table 5:  Constraints on the main and derived parameters (at 68% CL if not
otherwise stated) considering P18 in combination with BAO+SN and
BAO+SN+$p(H_{0})$ for the EMG model.

P18+BAO+SN | IG |  NMC$+$ |  NMC$-$ |  EMG ($Z=1$) |  EMG $(Z=-1)$
---|---|---|---|---|---
Planck high-$\ell$ TTTEEE | -1.2 | -1.3 | -1.6 | -0.8 | -0.6
Planck low-$\ell$ EE | 0.4 | 0.4 | 1 | 0.1 | -0.3
Planck low-$\ell$ TT | 0.2 | 0.3 | 0.3 | -0.2 | 0.5
Planck lensing | -0.2 | -0.2 | -0.3 | 0 | 0.8
BAO | 0.2 | 0.2 | 0.2 | 0 | -0.5
Pantheon | 0.1 | 0.1 | 0.1 | 0 | -0.2
Total | -0.5 | -0.5 | -0.3 | -0.9 | -0.3

P18+BAO+SN+$p(H_{0})$ | IG |  NMC$+$ |  NMC$-$ |  EMG ($Z=1$) |  EMG $(Z=-1)$
---|---|---|---|---|---
Planck high-$\ell$ TTTEEE | -1.2 | -0.8 | -1.4 | -0.7 | -6.5
Planck low-$\ell$ EE | 0.8 | 0.4 | 0.7 | 0.4 | -1.8
Planck low-$\ell$ TT | -0.1 | -0.2 | -0.2 | -0.2 | 1.1
Planck lensing | -0.1 | -0.1 | -0.2 | -0.1 | 0.6
BAO | 0.1 | 0.1 | 0 | 0 | 5
Pantheon | 0 | 0 | 0 | 0 | 0.4
$H_{0}$ | 0 | 0.2 | 0.6 | -17.8 | -13.7
Total | -0.5 | -0.4 | -0.5 | -18.4 | -14.9

Table 6:  Best-fit $\Delta\chi^{2}$ with respect to the $\Lambda$CDM model for
each dataset for the combination P18+BAO+SN (upper table) and
P18+BAO+SN+$p(H_{0})$ (lower table) for IG, NMC+, NMC-, and EMG.

##  Appendix B Background equations

Starting from Eq. 2.1, it is possible to write down the equations governing
the background evolution. Specializing to a spatially flat Friedmann-Lemaître-
Robertson-Walker (FLRW) universe described by the line element

$\mathrm{d}s^{2}=a^{2}(\eta)\left(-\mathrm{d}\eta^{2}+\mathrm{d}{\bf
x}^{2}\right)\,,$ (B.1)

where $\eta$ is the conformal time and ${\bf x}$ the spatial comoving
coordinate. The Einstein equations are obtained by varying the action 2.1 with
respect to the metric, they correspond to

$G^{\mu}_{\ \nu}=\frac{1}{F}\left[T^{\mu}_{\
\nu}-\frac{Z}{2}\nabla^{\mu}\sigma\nabla_{\nu}\sigma-g^{\mu}_{\
\nu}V+\left(\nabla^{\mu}\nabla_{\nu}-g^{\mu}_{\ \nu}\square\right)F\right]$
(B.2)

where the energy-momentum tensor for a perfect fluid is given by

$T_{\mu\nu}=pg_{\mu\nu}+(\rho+p)u_{\mu}u_{\nu}$ (B.3)

where a sum over all the species in the Universe is taken for granted, i.e.
$\rho\equiv\sum_{i}\rho_{i}$ and $p\equiv\sum_{i}p_{i}$.

From Eq. (B.2), the Friedmann equations in Jordan frame are as follows

$\displaystyle
3F\,\mathcal{H}^{2}=a^{2}\left(\rho+V\right)+\frac{Z\sigma^{\prime\,2}}{2}-3\mathcal{H}F^{\prime}$
(B.4)
$\displaystyle-2F\,\mathcal{H}^{\prime}=\frac{a^{2}}{3}\left(\rho+3p-2V\right)+\frac{2}{3}Z\sigma^{\prime\,2}+F^{\prime\prime}$
(B.5)

and the Einstein trace (the Ricci scalar) equation results

$a^{2}FR=a^{2}(\rho-3p)+4a^{2}V-3F_{\sigma}\left(\sigma^{\prime\prime}+2\mathcal{H}\sigma^{\prime}\right)-(Z+3F_{\sigma\sigma})\sigma^{\prime\,2}\,.$
(B.6)

Finally, the evolution equation of the scalar field $\sigma$ is governed by
the modified Klein-Gordon equation

$\displaystyle\sigma^{\prime\prime}+2\mathcal{H}\sigma^{\prime}-\frac{F_{\sigma}}{2ZF+3F_{\sigma}^{2}}\Bigl{[}$
$\displaystyle
a^{2}(\rho-3p)+4a^{2}\left(V-\frac{F}{2F_{\sigma}}V_{\sigma}\right)$
$\displaystyle-\left(Z+3F_{\sigma\sigma}\right)\sigma^{\prime\,2}\Bigr{]}=0\,.$
(B.7)

##  Appendix C Linear perturbed equations

In the synchronous gauge, up to linear order in perturbed quantities, the
perturbed FLRW metric is

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

where $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}$ can be express as a Fourier
integral as

$h_{ij}(\eta,\,{\bf x})=\int\mathrm{d}^{3}ke^{\imath{\bf k\cdot
x}}\left[\hat{k}_{i}\hat{k}_{j}h(\eta,\,{\bf
k})+\left(\hat{k}_{i}\hat{k}_{j}-\frac{\delta_{ij}}{3}\right)6\xi(\eta,\,{\bf
k})\right]$ (C.2)

where $\hat{k}_{i}=k_{i}/k$ with $k=|{\bf k}|$ and $h\equiv\delta^{ij}h_{ij}$
is the Fourier transform of trace of $h_{ij}(\eta,\,{\bf x})$. We follow the
conventions of Ref. [88].

###  C.1 The perturbed Einstein field equations

Splitting the Einstein tensor as the sum of the background (mean) part
$\bar{G}_{\mu\nu}$ and its corresponding perturbation $\delta G_{\mu\nu}$,
i.e. $G_{\mu\nu}=\bar{G}_{\mu\nu}+\delta G_{\mu\nu}$, the scalar perturbations
in synchronous gauge are usually presented as time-time, longitudinal time-
space, trace space-space, and longitudinal traceless space-space parts of the
Einstein equations in Fourier space as follow

$\displaystyle\xi k^{2}-\frac{{\cal H}h^{\prime}}{2}$
$\displaystyle=a^{2}\frac{\delta G^{0}_{\ 0}}{2}$ $\displaystyle=4\pi
Ga^{2}\delta T^{0}_{\ 0}\,,$ (C.3) $\displaystyle k^{2}\xi^{\prime}$
$\displaystyle=a^{2}\frac{\nabla^{i}\delta G^{0}_{\ i}}{2}$
$\displaystyle=4\pi Ga^{2}(\bar{\rho}+\bar{p})\theta\,,$ (C.4) $\displaystyle
h^{\prime\prime}+2{\cal H}h^{\prime}-2\xi k^{2}$ $\displaystyle=-a^{2}\delta
G^{i}_{\ i}$ $\displaystyle=-8\pi Ga^{2}\delta T^{i}_{\ i}\,,$ (C.5)
$\displaystyle h^{\prime\prime}+6\xi^{\prime\prime}+2{\cal
H}(h^{\prime}+6\xi^{\prime})-2k^{2}\xi$
$\displaystyle=-3a^{2}\left(\hat{k}_{i}\hat{k}_{j}-\frac{\delta_{ij}}{3}\right)\delta
G^{i}_{\ i}$ $\displaystyle=-24\pi Ga^{2}(\bar{\rho}+\bar{p})\Theta$ (C.6)

where we used the definition

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

and splitting the total energy density and pressure in a background and
perturbed parts, we obtain the following elements

$\displaystyle\bar{T}_{00}=a^{2}(\bar{\rho}+\delta\rho)\,,$ (C.10)
$\displaystyle\bar{T}_{0i}=-a^{2}(\bar{\rho}+\bar{p})v^{i}\,,$ (C.11)
$\displaystyle\bar{T}_{ij}=a^{2}\delta_{ij}\left(\bar{p}+\delta
p\right)+a^{2}\Sigma_{ij}\,.$ (C.12)

The first perturbed equation involving total density fluctuations reads

$\xi k^{2}-\frac{{\cal H}h^{\prime}}{2}=-(8\pi
G)a^{2}\frac{\delta\tilde{\rho}}{2F}$ (C.13)

with

$\displaystyle\delta\tilde{\rho}=$
$\displaystyle\delta\rho-\frac{h^{\prime}\bar{\sigma}^{\prime}\bar{F}_{\sigma}}{2a^{2}}+\frac{\delta\sigma^{\prime}}{a^{2}}\left(Z\bar{\sigma}^{\prime}-3{\cal
H}\bar{F}_{\sigma}\right)-\frac{\delta\sigma\bar{F}_{\sigma}}{a^{2}\bar{F}}\left[a^{2}\bar{\rho}+\frac{Z}{2}\bar{\sigma}^{\prime
2}\right.$
$\displaystyle\left.+a^{2}\left(V-V_{\sigma}\frac{F}{F_{\sigma}}\right)-3{\cal
H}\bar{F}_{\sigma}\bar{\sigma}^{\prime}+3{\cal
H}\frac{F_{\sigma\sigma}F}{F_{\sigma}}\bar{\sigma}^{\prime}+k^{2}F\right]\,.$
(C.14)

The second perturbed equation involving total velocity reads

$k^{2}\xi^{\prime}=(8\pi
G)a^{2}\frac{(\tilde{\rho}+\tilde{p})\tilde{\theta}}{2F}$ (C.15)

with

$\left(\tilde{\rho}+\tilde{p}\right)\tilde{\theta}=\left(\bar{\rho}+\bar{p}\right)\theta+\frac{k^{2}}{a^{2}}\left[\left(Z\bar{\sigma}^{\prime}-{\cal
H}F_{\sigma}+F_{\sigma\sigma}\bar{\sigma}^{\prime}\right)\delta\sigma+F_{\sigma}\delta\sigma^{\prime}\right]\,.$
(C.16)

The third perturbed equation involving total pressure reads

$h^{\prime\prime}+2{\cal H}h^{\prime}-2k^{2}\xi=-3(8\pi
G)a^{2}\frac{\delta\tilde{p}}{F}$ (C.17)

with

$\displaystyle\delta\tilde{p}=$ $\displaystyle\delta
p+\frac{h^{\prime}F^{\prime}}{3a^{2}}+\frac{Z}{a^{2}}\bar{\sigma}^{\prime}\delta\sigma^{\prime}-\delta
V+\frac{2}{3a^{2}}k^{2}\delta F+\frac{{\cal H}}{a^{2}}\delta
F^{\prime}+\frac{\delta F^{\prime\prime}}{a^{2}}$
$\displaystyle-\delta\sigma\frac{\bar{F}_{\sigma}}{a^{2}\bar{F}}\left(a^{2}\bar{p}+\frac{Z}{2}\bar{\sigma}^{\prime
2}-a^{2}\bar{V}+{\cal H}F^{\prime}+F^{\prime\prime}\right)\,.$ (C.18)

The fourth perturbed equation involving total shear reads

$h^{\prime\prime}+6\xi^{\prime\prime}+2{\cal
H}(h^{\prime}+6\xi^{\prime})-2k^{2}\xi=-3(8\pi
G)a^{2}\frac{(\tilde{\rho}+\tilde{p})\tilde{\Theta}}{F}$ (C.19)

with

$(\tilde{\rho}+\tilde{p})\tilde{\Theta}=(\bar{\rho}+\bar{p})\Theta+\frac{2k^{2}}{3a^{2}}\left(F_{\sigma}\delta\sigma+F^{\prime}\frac{h^{\prime}+6\xi^{\prime}}{2k^{2}}\right)\,.$
(C.20)

The perturbed Ricci scalar is given by

$\displaystyle a^{2}F\delta R=\,$ $\displaystyle
a^{2}(\delta\bar{\rho}-3\delta\bar{p})-\frac{3h^{\prime}F^{\prime}}{2}-2Z\bar{\sigma}^{\prime}\delta\sigma^{\prime}+4a^{2}V_{\sigma}\delta\sigma-6{\cal
H}\delta F^{\prime}-3\delta F^{\prime\prime}$ $\displaystyle-3k^{2}\delta
F-\frac{\delta\sigma
F_{\sigma}}{F}\left[a^{2}(\bar{\rho}-3\bar{p})-Z\bar{\sigma}^{\prime
2}+4a^{2}\bar{V}-3F^{\prime\prime}-6{\cal H}F^{\prime}\right]\,.$ (C.21)

###  C.2 The perturbed Klein-Gordon equation

The perturbed equation for the evolution of the scalar field perturbation
$\delta\sigma$ is

$Z\delta\sigma^{\prime\prime}+2{\cal
H}Z\delta\sigma^{\prime}+\left[Zk^{2}+a^{2}\left(V_{\sigma\sigma}-\frac{RF_{\sigma\sigma}}{2}\right)\right]\delta\sigma+Z\frac{h^{\prime}\bar{\sigma}^{\prime}}{2}-\frac{a^{2}F_{\sigma}}{2}\delta
R=0\,.$ (C.22)

##  Appendix D Initial conditions

The adiabatic initial condition for the background correspond to

$\displaystyle a(\tau)$ $\displaystyle=\sqrt{\frac{\rho_{\rm r\,0}}{3F_{\rm
ini}}}\tau\left[1+\frac{Z}{4}\omega\tau-\frac{5ZF_{\rm
ini,\sigma}^{2}\left(Z+3F_{\rm ini,\sigma\sigma}\right)}{64ZF_{\rm
ini}+96F_{\rm ini,\sigma}^{2}}\left(\omega\tau\right)^{2}\right]\,,$
$\displaystyle\mathcal{H}(\tau)$
$\displaystyle=\frac{1}{\tau}\left[1+\frac{Z}{4}\omega\tau-Z\frac{2F_{\rm
ini}+F_{\rm ini,\sigma}^{2}\left(8Z+15F_{\rm
ini,\sigma\sigma}\right)}{32ZF_{\rm ini}+48F_{\rm
ini,\sigma}^{2}}\left(\omega\tau\right)^{2}\right]\,,$
$\displaystyle\sigma(\tau)$ $\displaystyle=\sigma_{\rm ini}+\frac{3F_{\rm
ini,\sigma}}{4}\omega\tau-F_{\rm ini,\sigma}\frac{4ZF_{\rm ini}(2Z-3F_{\rm
ini,\sigma\sigma})+27F_{\rm ini,\sigma}^{2}\left(Z+F_{\rm
ini,\sigma\sigma}\right)}{32(2ZF_{\rm ini}+3F_{\rm
ini,\sigma}^{2})}\left(\omega\tau\right)^{2}$ (D.1)

where

$\omega=\frac{\rho_{\rm m,\,0}}{\sqrt{3\rho_{\rm r,\,0}}}\frac{2\sqrt{F_{\rm
ini}}}{2ZF_{\rm ini}+3F_{\rm ini,\sigma}^{2}}\,.$ (D.2)

For the cosmological fluctuations in the synchronous gauge, we have as
adiabatic initial conditions

$\displaystyle\delta_{\gamma}(k\,,\tau)=\delta_{\nu}(k\,,\tau)=\frac{4}{3}\delta_{\rm
b}(k\,,\tau)=\frac{4}{3}\delta_{\rm
c}(k\,,\tau)=-\frac{(k\tau)^{2}}{3}\left(1-\frac{Z\omega\tau}{5}\right)$ (D.3)
$\displaystyle\theta_{\nu}(k\,,\tau)=-\frac{k^{4}\tau^{3}}{36}\frac{23+4R_{\nu}}{15+4R_{\nu}}\left[1-\frac{3}{20}\frac{Z\left(275+50R_{\nu}+8R_{\nu}^{2}\right)F_{\rm
ini}+15(5-4R_{\nu})F_{\rm ini,\sigma}^{2}}{(15+2R_{\nu})(23+4R_{\nu})F_{\rm
ini}}\omega\tau\right]$ (D.4)
$\displaystyle\theta_{\gamma}(k\,,\tau)=\theta_{\rm
b}(k\,,\tau)=-\frac{k^{4}\tau^{3}}{36}\left[1-\frac{3}{20}\frac{Z(1-R_{\nu}+5R_{\rm
b})F_{\rm ini}+\frac{15}{2}R_{\rm b}F_{\rm
ini,\sigma}^{2}}{\left(1-R_{\nu}\right)F_{\rm ini}}\omega\tau\right]$ (D.5)
$\displaystyle\theta_{\rm c}(k\,,\tau)=0$ (D.6)
$\displaystyle\sigma_{\nu}(k\,,\tau)=\frac{2(k\tau)^{2}}{3(15+4R_{\nu})}\left[1+\frac{(-5+4R_{\nu})\left(2ZF_{\rm
ini}+3F_{\rm ini,\sigma}^{2}\right)}{8(15+2R_{\nu})F_{\rm
ini}}\omega\tau\right]$ (D.7) $\displaystyle
h(k\,,\tau)=\frac{\left(k\tau\right)^{2}}{2}\left(1-\frac{Z\omega\tau}{5}\right)$
(D.8)
$\displaystyle\eta(k\,,\tau)=1-\frac{\left(k\tau\right)^{2}}{12(15+4R_{\nu})}\left[5+4R_{\nu}-\frac{2Z(5+4R_{\nu})(65+4R_{\nu})F_{\rm
ini}+75(-5+4R_{\nu})F_{\rm ini,\sigma}^{2}}{20(15+2R_{\nu})F_{\rm
ini}}\omega\tau\right]$ (D.9)
$\displaystyle\delta\sigma(k\,,\tau)=-\frac{1}{16}F_{\rm
ini,\sigma}k^{2}\tau^{3}\omega\left[1-\frac{2Z(8Z-9F_{\rm
ini,\sigma\sigma})F_{\rm ini}+(45F_{\rm ini,\sigma\sigma}+48Z)F_{\rm
ini,\sigma}^{2}}{40ZF_{\rm ini}+60F_{\rm ini,\sigma}^{2}}\omega\tau\right]$
(D.10)

where $R_{\nu}=\rho_{\nu,\,0}/\rho_{\rm r,\,0}$ and $R_{\rm b}=\rho_{\rm
b,\,0}/\rho_{\rm m,\,0}$. These quantities reduce to induced gravity for $Z=1$
and $F=\xi\sigma^{2}$ [89], to a non-minimally coupled scalar field with
standard kinetic term for $Z=1$ and $F=N_{\rm Pl}^{2}+\xi\sigma^{2}$ [14], and
General Relativity for $Z=1$ and $F=M_{\rm Pl}^{2}$ [88].

##  Appendix E Comparison between BAO, and FS + BAO joint analysis

We compare here the results adding different datasets of galaxy information to
the Planck DR3 data such as the full shape (FS) of BOSS DR12 pre-reconstructed
power spectrum measurements [66, 67], BAO of BOSS DR12 post-reconstruction
power spectrum measurements [68], low-$z$ BAO measurements from SDSS DR7 6dF
and MGS [69, 70], Ly$\alpha$ BAO measurements from eBOSS [71, 72, 73], and
combination of those including the covariance among the DR12 datasets.

Figure 11: Marginalized joint 68% and 95% CL regions 2D parameter space using
the CMB P18 data for IG ($F=\xi\sigma^{2}$) in the phantom branch ($Z=-1$) in
combination with BAO from BOSS DR12 (red), BAO from BOSS DR12/SDSS DR7 6dF-
MGS/eBOSS (green), FS plus BAO from BOSS DR12 (orange), and FS plus BAO from
BOSS DR12/SDSS DR7 6dF-MGS/eBOSS (blue).

In Fig. 11, we show the marginalized posterior distributions of the
cosmological parameters for IG with P18 plus different combinations of the FS
and BAO measurements. We see that the posterior distributions are very robust
among the combination considered and that the addition of FS information to
the combination P18+BAO does not change the marginalized constraints for the
models studied here.

## References

  * Riess et al. [1998] Adam G. Riess et al.  Observational evidence from supernovae for an accelerating universe and a cosmological constant.  _Astron. J._ , 116:1009–1038, 1998.  doi: 10.1086/300499.
  * Perlmutter et al. [1999] S. Perlmutter et al.  Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae.  _Astrophys. J._ , 517:565–586, 1999.  doi: 10.1086/307221.
  * Carr et al. [2021] Anthony Carr, Tamara M. Davis, Daniel Scolnic, Khaled Said, Dillon Brout, Erik R. Peterson, and Richard Kessler.  The Pantheon+ Analysis: Improving the Redshifts and Peculiar Velocities of Type Ia Supernovae Used in Cosmological Analyses.  12 2021.
  * Aghanim et al. [2020a] N. Aghanim et al.  Planck 2018 results. I. Overview and the cosmological legacy of Planck.  _Astron. Astrophys._ , 641:A1, 2020a.  doi: 10.1051/0004-6361/201833880.
  * Alam et al. [2021] Shadab Alam et al.  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(8):083533, 2021.  doi: 10.1103/PhysRevD.103.083533.
  * 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(4):3275–3283, 2014.  doi: 10.1093/mnras/stu1050.
  * Sherwin et al. [2017] Blake D. Sherwin et al.  Two-season Atacama Cosmology Telescope polarimeter lensing power spectrum.  _Phys. Rev. D_ , 95(12):123529, 2017.  doi: 10.1103/PhysRevD.95.123529.
  * 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.  doi: 10.3847/1538-4357/ab4186.
  * Darwish et al. [2020] Omar 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(2):2250–2263, 2020.  doi: 10.1093/mnras/staa3438.
  * Asgari et al. [2021] Marika Asgari et al.  KiDS-1000 Cosmology: Cosmic shear constraints and comparison between two point statistics.  _Astron. Astrophys._ , 645:A104, 2021.  doi: 10.1051/0004-6361/202039070.
  * Abbott et al. [2022] T. M. C. Abbott et al.  Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing.  _Phys. Rev. D_ , 105(2):023520, 2022.  doi: 10.1103/PhysRevD.105.023520.
  * Cooke et al. [2018] Ryan J. Cooke, Max Pettini, and Charles C. Steidel.  One Percent Determination of the Primordial Deuterium Abundance.  _Astrophys. J._ , 855(2):102, 2018.  doi: 10.3847/1538-4357/aaab53.
  * 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, 2015.  doi: 10.1088/1475-7516/2015/08/017.
  * Rossi et al. [2019] Massimo Rossi, Mario Ballardini, Matteo Braglia, Fabio Finelli, Daniela Paoletti, Alexei A. Starobinsky, and Caterina Umiltà.  Cosmological constraints on post-Newtonian parameters in effectively massless scalar-tensor theories of gravity.  _Phys. Rev. D_ , 100(10):103524, 2019.  doi: 10.1103/PhysRevD.100.103524.
  * Poulin et al. [2019] Vivian Poulin, Tristan L. Smith, Tanvi Karwal, and Marc Kamionkowski.  Early Dark Energy Can Resolve The Hubble Tension.  _Phys. Rev. Lett._ , 122(22):221301, 2019.  doi: 10.1103/PhysRevLett.122.221301.
  * Agrawal et al. [2019] Prateek Agrawal, Francis-Yan Cyr-Racine, David Pinner, and Lisa Randall.  Rock ’n’ Roll Solutions to the Hubble Tension.  4 2019.  arXiv: 1904.01016.
  * Hart and Chluba [2020] Luke Hart and Jens Chluba.  Updated fundamental constant constraints from Planck 2018 data and possible relations to the Hubble tension.  _Mon. Not. Roy. Astron. Soc._ , 493(3):3255–3263, 2020.  doi: 10.1093/mnras/staa412.
  * Jedamzik and Pogosian [2020] Karsten Jedamzik and Levon Pogosian.  Relieving the Hubble tension with primordial magnetic fields.  _Phys. Rev. Lett._ , 125(18):181302, 2020.  doi: 10.1103/PhysRevLett.125.181302.
  * Ballardini et al. [2020] Mario Ballardini, Matteo Braglia, Fabio Finelli, Daniela Paoletti, Alexei A. Starobinsky, and Caterina Umiltà.  Scalar-tensor theories of gravity, neutrino physics, and the $H_{0}$ tension.  _JCAP_ , 10:044, 2020.  doi: 10.1088/1475-7516/2020/10/044.
  * Braglia et al. [2021] Matteo Braglia, Mario Ballardini, Fabio Finelli, and Kazuya Koyama.  Early modified gravity in light of the $H_{0}$ tension and LSS data.  _Phys. Rev. D_ , 103(4):043528, 2021.  doi: 10.1103/PhysRevD.103.043528.
  * Antony et al. [2022] Akhil Antony, Fabio Finelli, Dhiraj Kumar Hazra, and Arman Shafieloo.  Discordances in cosmology and the violation of slow-roll inflationary dynamics.  2022\.
  * Knox and Millea [2020] Lloyd Knox and Marius Millea.  Hubble constant hunter’s guide.  _Phys. Rev. D_ , 101(4):043533, 2020.  doi: 10.1103/PhysRevD.101.043533.
  * Di Valentino et al. [2021a] Eleonora Di Valentino et al.  Snowmass2021 - Letter of interest cosmology intertwined II: The hubble constant tension.  _Astropart. Phys._ , 131:102605, 2021a.  doi: 10.1016/j.astropartphys.2021.102605.
  * Di Valentino et al. [2021b] Eleonora Di Valentino et al.  Cosmology intertwined III: $f\sigma_{8}$ and $S_{8}$.  _Astropart. Phys._ , 131:102604, 2021b.  doi: 10.1016/j.astropartphys.2021.102604.
  * Di Valentino et al. [2021c] Eleonora Di Valentino, Olga Mena, Supriya Pan, Luca Visinelli, Weiqiang Yang, Alessandro Melchiorri, David F. Mota, Adam G. Riess, and Joseph Silk.  In the realm of the Hubble tension—a review of solutions.  _Class. Quant. Grav._ , 38(15):153001, 2021c.  doi: 10.1088/1361-6382/ac086d.
  * Perivolaropoulos and Skara [2022] Leandros Perivolaropoulos and Foteini Skara.  Challenges for $\Lambda$CDM: An update.  _New Astron. Rev._ , 95, 2022.  doi: 10.1016/j.newar.2022.101659.
  * Schöneberg et al. [2021] Nils Schöneberg, Guillermo Franco Abellán, Andrea Pérez Sánchez, Samuel J. Witte, Vivian Poulin, and Julien Lesgourgues.  The $H_{0}$ Olympics: A fair ranking of proposed models.  7 2021.  arXiv: 2107.10291.
  * Shah et al. [2021] Paul Shah, Pablo Lemos, and Ofer Lahav.  A buyer’s guide to the Hubble constant.  _Astron. Astrophys. Rev._ , 29(1):9, 2021.  doi: 10.1007/s00159-021-00137-4.
  * Abdalla et al. [2022] Elcio Abdalla et al.  Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies.  _JHEAp_ , 34:49–211, 2022.  doi: 10.1016/j.jheap.2022.04.002.
  * Aghanim et al. [2020b] N. Aghanim et al.  Planck 2018 results. VI. Cosmological parameters.  _Astron. Astrophys._ , 641:A6, 2020b.  doi: 10.1051/0004-6361/201833910.  [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * Riess et al. [2022] Adam 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(1):L7, 2022.  doi: 10.3847/2041-8213/ac5c5b.
  * Ballardini et al. [2016] Mario Ballardini, Fabio Finelli, Caterina Umiltà, and Daniela Paoletti.  Cosmological constraints on induced gravity dark energy models.  _JCAP_ , 05:067, 2016.  doi: 10.1088/1475-7516/2016/05/067.
  * Hill et al. [2020] J. Colin Hill, Evan McDonough, Michael W. Toomey, and Stephon Alexander.  Early dark energy does not restore cosmological concordance.  _Phys. Rev. D_ , 102(4):043507, 2020.  doi: 10.1103/PhysRevD.102.043507.
  * Braglia et al. [2020] Matteo Braglia, Mario Ballardini, William T. Emond, Fabio Finelli, A. Emir Gumrukcuoglu, Kazuya Koyama, and Daniela Paoletti.  Larger value for $H_{0}$ by an evolving gravitational constant.  _Phys. Rev. D_ , 102(2):023529, 2020.  doi: 10.1103/PhysRevD.102.023529.
  * Ivanov et al. [2020] Mikhail M. Ivanov, Evan McDonough, J. Colin Hill, Marko Simonović, Michael W. Toomey, Stephon Alexander, and Matias Zaldarriaga.  Constraining Early Dark Energy with Large-Scale Structure.  _Phys. Rev. D_ , 102(10):103502, 2020.  doi: 10.1103/PhysRevD.102.103502.
  * Rezazadeh et al. [2022] K. Rezazadeh, A. Ashoorioon, and D. Grin.  Cascading Dark Energy.  8 2022.
  * Smith et al. [2021] Tristan L. Smith, Vivian Poulin, José Luis Bernal, Kimberly K. Boddy, Marc Kamionkowski, and Riccardo Murgia.  Early dark energy is not excluded by current large-scale structure data.  _Phys. Rev. D_ , 103(12):123542, 2021.  doi: 10.1103/PhysRevD.103.123542.
  * Ballardini et al. [2022] Mario Ballardini, Fabio Finelli, and Domenico Sapone.  Cosmological constraints on the gravitational constant.  _JCAP_ , 06(06):004, 2022.  doi: 10.1088/1475-7516/2022/06/004.
  * Simon et al. [2022] Théo Simon, Pierre Zhang, Vivian Poulin, and Tristan L. Smith.  Updated constraints from the effective field theory analysis of BOSS power spectrum on Early Dark Energy.  8 2022.  arXiv: 2208.05930.
  * McDonough et al. [2022] Evan McDonough, Meng-Xiang Lin, J. Colin Hill, Wayne Hu, and Shengjia Zhou.  Early dark sector, the Hubble tension, and the swampland.  _Phys. Rev. D_ , 106(4):043525, 2022.  doi: 10.1103/PhysRevD.106.043525.
  * Gómez-Valent et al. [2022] Adrià Gómez-Valent, Ziyang Zheng, Luca Amendola, Christof Wetterich, and Valeria Pettorino.  Coupled and uncoupled early dark energy, massive neutrinos and the cosmological tensions.  7 2022.  arXiv: 2207.14487.
  * Reeves et al. [2022] Alexander Reeves, Laura Herold, Sunny Vagnozzi, Blake D. Sherwin, and Elisa G. M. Ferreira.  Restoring cosmological concordance with early dark energy and massive neutrinos?  7 2022.  arXiv: 2207.01501.
  * Boisseau et al. [2000] B. Boisseau, Gilles Esposito-Farese, D. Polarski, and Alexei A. Starobinsky.  Reconstruction of a scalar tensor theory of gravity in an accelerating universe.  _Phys. Rev. Lett._ , 85:2236, 2000.  doi: 10.1103/PhysRevLett.85.2236.
  * Caldwell [2002] R. R. Caldwell.  A Phantom menace?  _Phys. Lett. B_ , 545:23–29, 2002.  doi: 10.1016/S0370-2693(02)02589-3.
  * Cline et al. [2004] James M. Cline, Sangyong Jeon, and Guy D. Moore.  The Phantom menaced: Constraints on low-energy effective ghosts.  _Phys. Rev. D_ , 70:043543, 2004.  doi: 10.1103/PhysRevD.70.043543.
  * Nilles [1984] Hans Peter Nilles.  Supersymmetry, Supergravity and Particle Physics.  _Phys. Rept._ , 110:1–162, 1984.  doi: 10.1016/0370-1573(84)90008-5.
  * Pollock [1988] M. D. Pollock.  On the Initial Conditions for Superexponential Inflation.  _Phys. Lett. B_ , 215:635–641, 1988.  doi: 10.1016/0370-2693(88)90034-2.
  * Gannouji et al. [2006] Radouane Gannouji, David Polarski, Andre Ranquet, and Alexei A. Starobinsky.  Scalar-Tensor Models of Normal and Phantom Dark Energy.  _JCAP_ , 09:016, 2006.  doi: 10.1088/1475-7516/2006/09/016.
  * 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–376, 2003.  doi: 10.1038/nature01997.
  * Clifton et al. [2005] Timothy Clifton, John D. Barrow, and Robert J. Scherrer.  Constraints on the variation of G from primordial nucleosynthesis.  _Phys. Rev. D_ , 71:123526, 2005.  doi: 10.1103/PhysRevD.71.123526.
  * Muller and Biskupek [2007] Jurgen Muller and Liliane Biskupek.  Variations of the gravitational constant from lunar laser ranging data.  _Class. Quant. Grav._ , 24:4533–4538, 2007.  doi: 10.1088/0264-9381/24/17/017.
  * Belgacem et al. [2019] Enis Belgacem, Andreas Finke, Antonia Frassino, and Michele Maggiore.  Testing nonlocal gravity with Lunar Laser Ranging.  _JCAP_ , 02:035, 2019.  doi: 10.1088/1475-7516/2019/02/035.
  * Alonso et al. [2017] David Alonso, Emilio Bellini, Pedro G. Ferreira, and Miguel Zumalacárregui.  Observational future of cosmological scalar-tensor theories.  _Phys. Rev. D_ , 95(6):063502, 2017.  doi: 10.1103/PhysRevD.95.063502.
  * Ballardini et al. [2019] M. Ballardini, D. Sapone, C. Umiltà, F. Finelli, and D. Paoletti.  Testing extended Jordan-Brans-Dicke theories with future cosmological observations.  _JCAP_ , 05:049, 2019.  doi: 10.1088/1475-7516/2019/05/049.
  * Esposito-Farese and Polarski [2001] Gilles Esposito-Farese and D. Polarski.  Scalar tensor gravity in an accelerating universe.  _Phys. Rev. D_ , 63:063504, 2001.  doi: 10.1103/PhysRevD.63.063504.
  * Cooper and Venturi [1981] Fred Cooper and Giovanni Venturi.  Cosmology and Broken Scale Invariance.  _Phys. Rev. D_ , 24:3338, 1981.  doi: 10.1103/PhysRevD.24.3338.
  * Wetterich [1988] C. Wetterich.  Cosmologies With Variable Newton’s ’Constant’.  _Nucl. Phys. B_ , 302:645–667, 1988.  doi: 10.1016/0550-3213(88)90192-7.
  * Amendola [1999] Luca Amendola.  Scaling solutions in general nonminimal coupling theories.  _Phys. Rev. D_ , 60:043501, 1999.  doi: 10.1103/PhysRevD.60.043501.
  * Finelli et al. [2008] F. Finelli, A. Tronconi, and Giovanni Venturi.  Dark Energy, Induced Gravity and Broken Scale Invariance.  _Phys. Lett. B_ , 659:466–470, 2008.  doi: 10.1016/j.physletb.2007.11.053.
  * Lesgourgues [2011] Julien Lesgourgues.  The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview.  4 2011.  arXiv: 1104.2932.
  * Blas et al. [2011] Diego Blas, Julien Lesgourgues, and Thomas Tram.  The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes.  _JCAP_ , 07:034, 2011.  doi: 10.1088/1475-7516/2011/07/034.
  * Audren et al. [2013] Benjamin Audren, Julien Lesgourgues, Karim Benabed, and Simon Prunet.  Conservative Constraints on Early Cosmology: an illustration of the Monte Python cosmological parameter inference code.  _JCAP_ , 02:001, 2013.  doi: 10.1088/1475-7516/2013/02/001.
  * Brinckmann and Lesgourgues [2019] Thejs Brinckmann and Julien Lesgourgues.  MontePython 3: boosted MCMC sampler and other features.  _Phys. Dark Univ._ , 24:100260, 2019.  doi: 10.1016/j.dark.2018.100260.
  * Aghanim et al. [2020c] N. Aghanim et al.  Planck 2018 results. V. CMB power spectra and likelihoods.  _Astron. Astrophys._ , 641:A5, 2020c.  doi: 10.1051/0004-6361/201936386.
  * Aghanim et al. [2020d] N. Aghanim et al.  Planck 2018 results. VIII. Gravitational lensing.  _Astron. Astrophys._ , 641:A8, 2020d.  doi: 10.1051/0004-6361/201833886.
  * Gil-Marín et al. [2016] Héctor Gil-Marín et al.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: RSD measurement from the LOS-dependent power spectrum of DR12 BOSS galaxies.  _Mon. Not. Roy. Astron. Soc._ , 460(4):4188–4209, 2016.  doi: 10.1093/mnras/stw1096.
  * D’Amico et al. [2021] Guido D’Amico, Leonardo Senatore, and Pierre Zhang.  Limits on $w$CDM from the EFTofLSS with the PyBird code.  _JCAP_ , 01:006, 2021.  doi: 10.1088/1475-7516/2021/01/006.
  * Alam et al. [2017] Shadab Alam et al.  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(3):2617–2652, 2017.  doi: 10.1093/mnras/stx721.
  * Beutler et al. [2011] Florian Beutler, Chris Blake, Matthew Colless, D. Heath Jones, Lister Staveley-Smith, Lachlan Campbell, Quentin Parker, Will Saunders, and Fred Watson.  The 6dF Galaxy Survey: Baryon Acoustic Oscillations and the Local Hubble Constant.  _Mon. Not. Roy. Astron. Soc._ , 416:3017–3032, 2011.  doi: 10.1111/j.1365-2966.2011.19250.x.
  * Ross et al. [2015] Ashley J. Ross, Lado Samushia, Cullan Howlett, Will J. Percival, Angela Burden, and Marc 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(1):835–847, 2015.  doi: 10.1093/mnras/stv154.
  * de Sainte Agathe et al. [2019] Victoria 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.  doi: 10.1051/0004-6361/201935638.
  * Blomqvist et al. [2019] Michael Blomqvist et al.  Baryon acoustic oscillations from the cross-correlation of Ly$\alpha$ absorption and quasars in eBOSS DR14.  _Astron. Astrophys._ , 629:A86, 2019.  doi: 10.1051/0004-6361/201935641.
  * Cuceu et al. [2019] Andrei Cuceu, James Farr, Pablo Lemos, and Andreu Font-Ribera.  Baryon Acoustic Oscillations and the Hubble Constant: Past, Present and Future.  _JCAP_ , 10:044, 2019.  doi: 10.1088/1475-7516/2019/10/044.
  * Scolnic et al. [2018] D. M. Scolnic et al.  The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample.  _Astrophys. J._ , 859(2):101, 2018.  doi: 10.3847/1538-4357/aab9bb.
  * 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–971, 2008.  doi: 10.1016/j.cpc.2008.02.015.
  * 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–242, 2018.  doi: 10.1016/j.cpc.2018.06.022.
  * Zumalacarregui [2020] Miguel Zumalacarregui.  Gravity in the Era of Equality: Towards solutions to the Hubble problem without fine-tuned initial conditions.  _Phys. Rev. D_ , 102(2):023523, 2020.  doi: 10.1103/PhysRevD.102.023523.
  * Solà Peracaula et al. [2020] Joan Solà Peracaula, Adrià Gómez-Valent, Javier de Cruz Pérez, and Cristian Moreno-Pulido.  Brans–Dicke cosmology with a $\Lambda$-term: a possible solution to $\Lambda$CDM tensions.  _Class. Quant. Grav._ , 37(24):245003, 2020.  doi: 10.1088/1361-6382/abbc43.
  * Ballardini and Finelli [2022] Mario Ballardini and Fabio Finelli.  Type Ia supernovae data with scalar-tensor gravity.  _Phys. Rev. D_ , 106(6):063531, 2022.  doi: 10.1103/PhysRevD.106.063531.
  * Garriga and Mukhanov [1999] Jaume Garriga and Viatcheslav F. Mukhanov.  Perturbations in k-inflation.  _Phys. Lett. B_ , 458:219–225, 1999.  doi: 10.1016/S0370-2693(99)00602-4.
  * Gasperini and Veneziano [2003] M. Gasperini and G. Veneziano.  The Pre - big bang scenario in string cosmology.  _Phys. Rept._ , 373:1–212, 2003.  doi: 10.1016/S0370-1573(02)00389-7.
  * Arkani-Hamed et al. [2004] Nima Arkani-Hamed, Hsin-Chia Cheng, Markus A. Luty, and Shinji Mukohyama.  Ghost condensation and a consistent infrared modification of gravity.  _JHEP_ , 05:074, 2004.  doi: 10.1088/1126-6708/2004/05/074.
  * Nicolis et al. [2009] Alberto Nicolis, Riccardo Rattazzi, and Enrico Trincherini.  The Galileon as a local modification of gravity.  _Phys. Rev. D_ , 79:064036, 2009.  doi: 10.1103/PhysRevD.79.064036.
  * Deffayet et al. [2009] C. Deffayet, Gilles Esposito-Farese, and A. Vikman.  Covariant Galileon.  _Phys. Rev. D_ , 79:084003, 2009.  doi: 10.1103/PhysRevD.79.084003.
  * Silva and Koyama [2009] Fabio P Silva and Kazuya Koyama.  Self-Accelerating Universe in Galileon Cosmology.  _Phys. Rev. D_ , 80:121301, 2009.  doi: 10.1103/PhysRevD.80.121301.
  * Kobayashi [2010] Tsutomu Kobayashi.  Cosmic expansion and growth histories in Galileon scalar-tensor models of dark energy.  _Phys. Rev. D_ , 81:103533, 2010.  doi: 10.1103/PhysRevD.81.103533.
  * Ferrari and et al. [in preparation (2023] Angelo Ferrari and et al.  in preparation (2023).
  * Ma and Bertschinger [1995] Chung-Pei Ma and Edmund Bertschinger.  Cosmological perturbation theory in the synchronous and conformal Newtonian gauges.  _Astrophys. J._ , 455:7–25, 1995.  doi: 10.1086/176550.
  * 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.  doi: 10.1016/j.dark.2019.100307.