# Neutrino mass and kinetic gravity braiding degeneracies

Gabriela Garcia-Arroyo  Jorge L. Cervantes-Cota  Ulises Nucamendi

###### Abstract

Modified theories of gravity yield an effective dark energy in the background
dynamics that achieves an accelerated expansion of the universe. In addition,
they present a fifth force that induces gravitational signatures in structure
formation, and therefore in the matter power spectrum and related statistics.
On the other hand, massive neutrinos suppress the power spectrum at scales
that also modified gravity enhances it, so a degeneration of these effects has
been recognized for some gravity models. In the present work, we study both
effects using kinetic gravity braiding ($n$KGB) models to find that in spite
of some degeneracies, the role of the fifth force at very large scales
imprints a bump in the matter power spectrum as a distinctive signature of
this model and, therefore, acts as a smoking gun that seems difficult to match
within the present knowledge of power spectra. These models result
interesting, however, since the $n=1$ presents no $H_{0}$ tension, and all
$n$KGB studied here present no $\sigma_{8}$ tension and, in addition, a null
neutrino mass is excluded.

##  1 Introduction

The standard model of cosmology is based on the existence of dark matter and
dark energy. The first was initially identified as a necessary ingredient to
understand the dynamics in galaxies and clusters [1], and later to encompass
different properties of the Cosmic Microwave Background (CMB), such as its
acoustic peaks [2], and to account for the Baryon Acoustic Oscillations (BAO)
[3]; it is also important for the growth of structure in the Universe at
different scales. The second is important to achieve an accelerated expansion
since the last few Giga years [4, 5], and this impacts different observables,
such as clustering properties as given by the matter power spectrum. However,
none of these dark components have an accepted origin so far, hence different
approaches have been undertaken. As a possible alternative to dark energy are
Modified Gravity (MG) models that apart to generate the accelerated expansion,
a fifth force associated to their scalar degree of freedom yields distinctive
cosmological features [6]. In particular, the linear cosmic structure growth
is modified [7], depending now on time and scale, in contrast to the one in
the $\Lambda$CDM model that is subject only to time evolution.

MG models act mainly in the late universe, when their resulting effective dark
energy accelerates the background expansion and its fifth force influences the
clustering properties. This in turn has an impact in the CMB statistics at
large angles where the cosmic variance makes difficult to distinguish among
models; however, Integrated Sachs-Wolfe (ISW) cross correlated with galaxy
clustering helps to constrain models [8, 9, 10]. On the other hand, structure
formation offers a valuable cosmic probe to test gravity on large scales, and
therefore to be able to distinguish among different MG models. This can be
understood because different fifth forces provoke characteristic effects in
the cosmic, statistical observables, such as the power spectrum (PS) or the
two-point correlation function, and general N-point statistics [11]. However,
MG is not the only factor that imprints special features in these statistics.
Beyond the $\Lambda$CDM model, massive neutrinos have been recognized to
produce a reduction in the structure growth due to its free streaming [12]
that in turn depends on the sum of their masses; CMB constraints pose values
to less than a fraction of 1 eV, where the exact constraints depend on the
data set considered [2]. It is then expected that neutrinos have a range of
masses where MG models also yield structure growth effects. But it has been
found that MG and massive neutrinos work in opposite directions: whereas
neutrinos generate a PS suppression depending on the free streaming scale, MG
models cause a PS enhancement due to the attractive fifth force [13]. In
recent years the interest to study the joint effects of MG and massive
neutrinos has risen, partially because the stage-IV generation of galaxy
surveys and weak lensing probes such as DESI [14], EUCLID [15], or LSST [16]
will reach percent-level accuracy in two-point statistics observables, posing
the interesting challenge to discriminate among different effects due to
gravitational models and neutrino masses.

The effects of neutrinos in the CMB TT power spectrum have been considered in
[17] for the Hu-Sawicki (HS) f(R) gravity model [18]. It is shown that the
presence of neutrinos and a fifth force act in the same direction to lower the
anisotropies for low multipoles; higher multipoles are scarcely affected. But
again the cosmic variance prevent us for discriminating these effects. On the
other hand, the effect of massive neutrinos in the PS for different MG models
shows that a higher neutrino mass (sum of masses) is allowed in MG to render a
similar General Relativity (GR) behavior with a smaller neutrino mass. This
happens because MG rises the matter power spectrum, $P_{\delta\delta}$,
whereas the neutrino free streaming lowers it [13, 19]; also, a similar effect
happens for the weak lensing convergence power spectrum [20]. Whereas an exact
cancellation of these effects does not happen, it seems difficult to
distinguish among them. Fortunately, peculiar velocities help to discriminate:
the density-velocity and velocity-velocity power spectra, $P_{\delta\theta}$
and $P_{\theta\theta}$, are different in GR and MG with or without massive
neutrinos; in ref. [19] this analysis was made for HS model. In the same line
of thought and closer to observations, redshift space distortions (RSD) serve
to sort out these effects. This has been applied again to HS models in refs.
[21, 22], where analyses are made at different redshifts.

Other sources of possible parameter degeneracies, such as baryon effects has
been recently considered for HS models and CPL dark energy parametrizations
[23, 24] to reach the levels of accuracy required to test gravity in the
above-mentioned probes [25]; constraints to the CPL model with massive
neutrinos have also been discussed [26]. Other studies consider the same HS
model to look for degeneracies among neutrinos masses and the halo mass
function and the bias [27], however, disentangling the effects is not an easy
task, again for the level of accuracy required. Here, voids can help in the
analysis [28]. As mentioned, almost all above works considered the HS model.
The influence of neutrinos in other gravity models is less common, but in ref.
[29] the CMB and matter PS for Galileon gravity models are discussed, where
again similar results were found, the mass of neutrinos suppresses the
enhancement of the fifth force.

Clustering effects in different gravity models can be different and they have
to be studied case by case; another approach is to treat these models through
parametrizations [30]. The latter approach serves to characterize general
features of phenomenological models, whereas the former intends to find
consequences of particular models. In the present work we deal with a
particular model, a kinetic gravity braiding (KGB) model that is part of the
full Horndeski gravity, known to yield second order gravitational field
equations. We will consider a model that is much simpler than the full
Horndeski. Its Lagrangian has a kinetic term that depends only on its kinetic
energy ($X$), and not on the scalar field itself, and the coupling to the
D’Alembertian of the scalar field, the braiding, also depends only on $X$.
This model yields a velocity of light for gravitational wave propagation, so
it passes the test of gravity wave speed inferred from the almost simultaneous
arrival of light and gravity waves from the merging of neutron stars, among
other properties considered in refs. [31, 32, 33, 34]; test against cosmic
data for this model are presented in refs. [35, 36].

We will work particularly with $n$KGB models [33] that possess only one or two
free parameters, depending on model choices. They are interesting because
within the same setup one obtains the cubic Galileon model for $n=1$ and one
recovers the $\Lambda$CDM model at background and first-order perturbation
levels for $n\rightarrow\infty$. This model has two parameters, a power law
exponent and its strength. In general, the cosmological background dynamics is
not exactly that of $\Lambda$CDM, but since the models have an attractor
solution, they will tend to a de Sitter phase, and in fact perturbations will
help to distinguish among them. As it happens in other gravity models that
affect large scales, we expect $n$KGB models will also modify linear power
spectra, and therefore we find interesting to consider its gravitational
effects and those of massive neutrinos, to see if they partially cancel each
other, as in the other models mentioned above.

This work is organized as follows: In section 2 we set the theoretical
framework of KGB models, to later understand, in section 3, its effects and
degeneracies with massive neutrinos in the cosmological observables. Then, we
fit these models to current cosmic data in section 4 to understand its role in
the $H_{0}$ and $\sigma_{8}$ tensions. Finally, our final comments and
conclusions are discussed in section 5.

##  2 Kinetic gravity braiding models

KGB models are the most general, minimally-coupled models yielding second
order field equations. Its general action can be written as [31, 32]:

$\displaystyle S_{\text{braiding}}=\int
d^{4}x\sqrt{-g}\left[\frac{M^{2}_{pl}}{2}R+K(\phi,X)-G(\phi,X)\square\phi+{\cal
L}_{\text{mat}}\right],$ (2.1)

where $R$ is the curvature scalar, $K=K(\phi,X)$ and $G=G(\phi,X)$ are
functions of the scalar field $\phi$ and of its kinetic energy density
$X\equiv-g^{\mu\nu}(\nabla_{\mu}\phi)(\nabla_{\nu}\phi)/2$, while
$\square\phi\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\phi$. The term ${\cal
L}_{m}$ stands for the Lagrangian of the matter degrees of freedom. The
gravitational field equations that can be derived from the above action (2.1)
read:

$G_{\mu\nu}=\frac{1}{M^{2}_{pl}}\left[T_{\mu\nu}^{(m)}+T_{\mu\nu}^{(\phi)}\right],$
(2.2)

where $T_{\mu\nu}^{(m)}$ is the total energy momentum tensor, containing dark
matter, baryons, neutrinos, and photons, while $T_{\mu\nu}^{(\phi)}$
represents the energy momentum tensor for the Galileon field $\phi$,

$\displaystyle
T_{\mu\nu}^{(\phi)}=K_{X}(\nabla_{\mu}\phi)(\nabla_{\nu}\phi)+g_{\mu\nu}K+g_{\mu\nu}(\nabla_{\beta}G)(\nabla^{\beta}\phi)$
$\displaystyle-\left[(\nabla_{\mu}G)(\nabla_{\nu}\phi)+(\nabla_{\nu}G)(\nabla_{\mu}\phi)\right]-G_{X}(\square\phi)(\nabla_{\mu}\phi)(\nabla_{\nu}\phi).$
(2.3)

Variation of the action (2.1) with respect to the Galileon scalar field
produces the equation,

$\displaystyle
K_{\phi}+\nabla^{\alpha}\left(K_{X}\nabla_{\alpha}\phi\right)+2G_{\phi\phi}X-2G_{\phi}\square\phi$
$\displaystyle+G_{X\phi}\left[2X\square\phi+2(\nabla^{\mu}\phi)(\nabla_{\mu}\nabla_{\nu}\phi)\nabla^{\nu}\phi\right]$
$\displaystyle+G_{X}\left[(\nabla_{\mu}\nabla_{\nu}\phi)^{2}+R_{\mu\nu}(\nabla^{\mu}\phi)\nabla^{\nu}\phi-(\square\phi)^{2}\right]$
$\displaystyle-
G_{XX}\left[(\nabla_{\nu}X)\nabla^{\nu}X-\nabla^{\mu}\phi(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\nu}\phi)\square\phi\right]$
$\displaystyle=$ $\displaystyle 0\,,$ (2.4)

where in eqs. (2.3) and (2.4) we have defined the first partial derivatives of
the coupling functions as,

$K_{\phi}\equiv\frac{\partial
K(\phi,X)}{\partial\phi},\,K_{X}\equiv\frac{\partial K(\phi,X)}{\partial
X},\,G_{\phi}\equiv\frac{\partial
G(\phi,X)}{\partial\phi},\,G_{X}\equiv\frac{\partial G(\phi,X)}{\partial X},$
(2.5)

and the corresponding second partial derivatives as,

$G_{\phi\phi}\equiv\frac{\partial^{2}G(\phi,X)}{\partial\phi^{2}},\quad
G_{X\phi}\equiv\frac{\partial^{2}G(\phi,X)}{\partial X\partial\phi},\quad
G_{XX}\equiv\frac{\partial^{2}G(\phi,X)}{\partial X^{2}}.$ (2.6)

We note that in the action (2.1) the interaction term $G\square\phi$ contains
the coupling of $G$ with first-order derivatives of both the metric and scalar
field stemming from the Christoffel symbols contained in the D’Alembertian
operator. These couplings produce a class of mixing named kinetic braiding, in
which the essential characteristic is determined by the kinetic dependence of
$G=G(X)$. This implies that the scalar field equation of motion and its energy
momentum tensor contain second derivatives of the metric in an essential way:
there exists no Einstein frame where the kinetic terms are diagonalized [31].
There are many possibilities to choose the $K$ and $G$ functions, we chose
$G=G(X)$ because is precisely the coupling of $X$ with $\square\phi$ that is
responsible of the braiding mixing; we also chose the kinetic term as a
function of the kinetic energy only, $K=K(X)$, to avoid extra complexities, as
it can be seen in equation (2.4). As a consequence, the action (2.1) is
invariant under the change $\phi\longrightarrow\phi+b$ (where b is a
constant); this shift symmetry gives a conserved Noether current $J_{\mu}$
given by,

$J_{\mu}=\left(K_{X}-G_{X}\square\phi\right)\nabla_{\mu}\phi+G_{X}(\nabla_{\mu}\nabla_{\nu}\phi)\nabla^{\nu}\phi\,.$
(2.7)

The conservation equation for the Noether current $\nabla_{\mu}J^{\mu}=0$
gives directly the Galileon field equation (2.4) with $K_{\phi}=G_{\phi}=0$,

$\displaystyle\nabla^{\alpha}\left(K_{X}\nabla_{\alpha}\phi\right)+G_{X}\left[(\nabla_{\mu}\nabla_{\nu}\phi)^{2}+R_{\mu\nu}(\nabla^{\mu}\phi)\nabla^{\nu}\phi-(\square\phi)^{2}\right]$
$\displaystyle-
G_{XX}\left[(\nabla_{\nu}X)\nabla^{\nu}X-\nabla^{\mu}\phi(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\nu}\phi)\square\phi\right]=0\,.$
(2.8)

Before we specify the $K(X)$ and $G(X)$ functions we would like to show the
Friedmann equations. For the spatially flat Friedmann-Robertson-Walker metric
we have:

$\displaystyle 3H^{2}$ $\displaystyle=$ $\displaystyle(8\pi
G)\left[\rho_{m}+\rho_{r}+\rho_{\phi}\right]\quad,$ (2.9) $\displaystyle
2\dot{H}+3H^{2}$ $\displaystyle=$ $\displaystyle-(8\pi
G)\left[P_{r}+P_{\phi}\right]\quad,$ (2.10)

where the dot denotes the derivative with respect to the cosmic time, $H$ is
the Hubble parameter, and the matter and radiation densities obey standard
fluid equations. The density and pressure of the Galileon field are given by:

$\displaystyle\rho_{\phi}$ $\displaystyle=$
$\displaystyle-K+K_{X}\dot{\phi}^{2}+3G_{X}H\dot{\phi}^{3}\quad,$ (2.11)
$\displaystyle P_{\phi}$ $\displaystyle=$ $\displaystyle
K-G_{X}\dot{\phi}^{2}\ddot{\phi}\quad.$ (2.12)

The scalar field equation (2.8) becomes,

$\displaystyle
K_{X}\left(\ddot{\phi}+3H\dot{\phi}\right)+K_{XX}\ddot{\phi}\dot{\phi}^{2}+3G_{X}\left(2H\dot{\phi}\ddot{\phi}+3H^{2}\dot{\phi}^{2}+\dot{H}\dot{\phi}^{2}\right)+3G_{XX}H\dot{\phi}^{3}\ddot{\phi}$
$\displaystyle=$ $\displaystyle 0\,.$ (2.13)

In the rest of this work we consider the $n$KGB model [33]:

$\displaystyle K(X)$ $\displaystyle=$ $\displaystyle-X\quad,$ (2.14)
$\displaystyle G(X)$ $\displaystyle=$ $\displaystyle
g^{(2n-1)/2}\,\Lambda\left(\frac{X}{\Lambda^{4}}\right)^{n}\quad,$ (2.15)

where $g$ is a dimensionless constant that controls the strength of the
braiding and $\Lambda$ has dimensions of energy weighted by the Hubble
constant and the Planck mass such that
$\Lambda^{4n-1}=H_{0}^{2n}M_{pl}^{2n-1}$ [37]. The election $n=1$ corresponds
to the covariant, cubic galileon, that is in big tension with ISW cross
correlated to clustering data [9], but we included it here for completeness.

$n$KGB models are interesting because the sound speed for scalar perturbations
can be smaller or bigger than light velocity, see Eqs. (A.11)-(A.12) and
subsequent paragraph of [31], hence, its clustering properties are nontrivial;
another important characteristic is the existence of an accelerated attractor
exhibiting phantom behaviour in its evolution that tends to a de Sitter phase,
as we will discuss below. In attractor phase and for models with $n>1/2$, it
has been proven that this model does not have ghosts and instabilities for the
cosmological perturbations [31, 35]. Important for our work, their background
and linear cosmological perturbations tend to the $\Lambda$CDM model for
$n\rightarrow\infty$. Also, the sound speed of the tensor modes is equal to
that of light since it is not affected by the braiding coupling, see the
related discussion in Subsection 4.2 of [38]. Additionally, these models have
a well-posed initial value problem and have a Vainshtein mechanism that
recover GR at Solar System scales [31].

In the flat FRW spacetime, the conserved charge density of the Noether current
(2.7) becomes

$J_{0}=\dot{\phi}(3\dot{\phi}G_{X}H-1)\,.$ (2.16)

It has been shown in ref. [31] that the simplest attractor solutions are found
under the condition $J_{0}=0$. From (2.16) we have two branches,

$\displaystyle\dot{\phi}$ $\displaystyle=$ $\displaystyle 0\,,$ (2.17)
$\displaystyle\dot{\phi}$ $\displaystyle=$ $\displaystyle\frac{1}{3G_{X}H}\,,$
(2.18)

the branch $\dot{\phi}=0$ develops ghosts in its cosmological perturbations
[35] and it is not of physical interest. The second branch (2.18) is self-
accelerating and well behaved under cosmological perturbations even though its
equation of state is of phantom type. Let’s mention that an extensive analysis
of scalar cosmological perturbations (without the inclusion of massive
neutrinos) around this background self-accelerating attractor has been done in
[35]. An important property is that its asymptotic evolution describes a de
Sitter attractor that can be analysed through the total equation of state

$w=-1-\frac{2}{3}\frac{\dot{H}}{H^{2}}=-1+\frac{(2n-1)}{3}\frac{(3\Omega_{m}+4\Omega_{r})}{(2n-\Omega_{m}-\Omega_{r})}\,,$
(2.19)

where the $\Omega_{i}$ are functions of time. This expression encodes the
effective evolution of the background dynamics. Evaluated today, $w=-0.82$ for
$n=1$ and $w=-0.7$ for $n\rightarrow\infty$. In this work, we do not force the
background evolution to the self-accelerating attractor, but we take an
arbitrary background setup. However, for some of models we will recover the
attractor solution, when parameters fit to cosmological data, as explained
below.

We now proceed to briefly describe the field equations for the linear
cosmological perturbations of the corresponding $n$KGB model around the FRW
background spacetime. The perturbed metric line element is written in the
Newtonian gauge as

$ds^{2}=-(1+2\Psi)dt^{2}+a^{2}(t)(1+2\Phi)\delta_{ij}dx^{i}dx^{j}\quad,$
(2.20)

while the scalar field is expanded by the background solution $\phi(t)$ and
its perturbation $\delta\phi$ as

$\phi(x^{i},t)=\phi(t)+\delta\phi(x^{i},t)\,.$ (2.21)

For the perturbations of the energy momentum tensor we have,

$\displaystyle\delta T^{i}_{\,\,j}$ $\displaystyle=$
$\displaystyle\delta^{i}_{j}(\delta p)+(\partial_{i}\partial_{j}\Pi)\,,$
(2.22) $\displaystyle\delta T^{i}_{\,\,t}$ $\displaystyle=$
$\displaystyle-\frac{1}{a^{2}}\partial_{i}(\delta q)\,,$ (2.23)
$\displaystyle\delta T^{t}_{\,\,i}$ $\displaystyle=$
$\displaystyle\partial_{i}(\delta q)\,,$ (2.24) $\displaystyle\delta
T^{t}_{\,\,t}$ $\displaystyle=$ $\displaystyle-(\delta\rho)\,,$ (2.25)

where $\delta\rho$, $\delta p$, $\delta q$ denote the perturbations of the
total energy density, pressure, and velocity field, respectively, while $\Pi$
describes the total anisotropic stress.

The perturbed cosmological Einstein equations are given by equations (4.2),
(4.4), (4.6) of reference [33] together with the equation (4.5) modified here
to incorporate the anisotropic stress $\Pi$,

$M^{2}_{pl}(\Psi+\Phi)=-a^{2}\Pi\,,$ (2.26)

while the scalar field equation is described by expression (4.7) in [33] and
it will not be rewritten here. Note that in this work $\Pi\neq 0$, though the
gravitational part has no slip, we consider massive neutrinos that have a
small anisotropic stress $\Pi_{\nu}$, apart from those of photons.

Finally, the perturbed energy momentum conservation equations for every matter
component are,

$\dot{\delta\rho_{i}}+3H\left[\delta\rho_{i}+\delta
p_{i}+\frac{1}{3}\nabla^{2}\Pi_{i}\right]+3(\rho_{i}+p_{i})\dot{\Phi}+\frac{1}{a^{2}}\nabla^{2}(\delta
q_{i})=0\,,$ (2.27) $\dot{\delta q_{i}}+3H\delta
q_{i}+(\rho_{i}+p_{i})\Psi+\delta p_{i}+\nabla^{2}\Pi_{i}=0\,,$ (2.28)

where $i$ denotes CDM, baryons, radiation, and neutrinos.

##  3 MG and neutrino mass degeneracies

The effects of massive neutrinos on the CMB temperature power spectrum are
four [39]: i) A change in the nonrelativistic energy density induces changes
between scales on the last scattering surface and angles on the sky that
affect the angular diameter distance to last scattering. And, on the other
hand, it changes the late ISW effect, that depends on the redshift at matter-
to-DE equality; ii) variation of the metric fluctuations when they become
nonrelativistic; iii) large neutrino masses tend to affect less the CMB
photons due to LSS weak lensing; iv) neutrinos that become nonrelativistic
earlier, due to smaller momenta, are imprinted in high-$\ell$ multipoles. The
first and the third effects are more important to put bounds on the masses. On
the other hand, MG changes the metric potentials involved in the evolution of
the perturbations at large scales and therefore changes the late ISW effect
and the clustering properties of matter that induces differences in the
lensing properties.

Here we show results for $n$KGB models, using hi_class [40] to compute the
spectrum for different $n$ powers, see eqs. (2.14) and (2.15), where we
assumed three species of neutrino with degenerate mass, leaving all other
parameters fixed; for comparison we assume a Planck base-$\Lambda$CDM model
with $M_{\nu}\equiv\Sigma m_{\nu}=0.06$eV. In figure 1 one observes a larger
neutrino mass lowers low multipoles, but braiding gravity enhances them; a
similar behavior is observed for the $n=1$ model in [9]. Due to the point ii)
above, a dip is observed at around $\ell\sim 40-200$ (for the $\Lambda$CDM
version see [39]), that here has two counter effects: whereas massive
neutrinos slightly change potentials wells when they become nonrelativistic,
the gravity effect works in the other direction at these scales; this latter
effect is more prominent for smaller $n$. The oscillatory pattern for
$\ell>200$, known in massive neutrino $\Lambda$CDM cosmologies, see point iii)
above, is also present here, the more pronounced for smaller $n$. One observes
that the $n=3$ curve is closer to $\Lambda$CDM than the $n=1$ curve, since
$n$KGB tends to GR as $n$ grows. In other MG theories the effect induced by
the gravity model could be different to the presented in this work, for
example, in ref. [17] results are shown for the HS model and the deviations of
the CMB TT angular power spectrum looks different from figure 1, but the
effect caused by the sum of neutrino mass is the same.

Figure 1: CMB TT power spectrum for $n$KGB models with $n=1,3$ and different
sum of neutrino masses $M_{\nu}$. Solid lines correspond to
$M_{\nu}=0.06\,\rm{eV}$ and dotted lines to $M_{\nu}=0.23\,\rm{eV}$. Colors
denote different $n$ values. Data points are those from Planck 2018 [2]. The
bottom panel shows the relative difference to $\Lambda$CDM.

We now look for effects on large scale structure observables. Different
combinations of cosmological parameters can cause similar outcomes in an
observable, and this may prevent us to distinguish which parameter is
responsible for the observed effects. In the case of massive neutrinos, it is
well known that they produce a decrease in dark matter perturbations which in
turn is manifested in the matter PS that presents a reduction [12] at
wavenumbers ($k$) between $0.01$ and $1$ $h$ Mpc-1 depending on their masses.
On the other hand, modifications of gravity can also alter the PS, since the
fifth force involved increases the gravitational potential, and the PS grows.

These two effects are shown in figure 2 using linear power spectra at $z=1$
(left panel) and $z=0$ (right panel) for the $n$KGB models ($n=1,3$),
discussed in section 2 employing equations (2.14, 2.15), where the value of
$g$ is fixed by the attractor solution to $0.027$ ($n=1$) and $0.272$ ($n=3$),
and $\Omega_{m}=0.316$ in both cases as the best fit of CMB
TT,TE,EE+lowE+lensing from Planck 2018 [2]. Here we plot the relative matter
PS with respect our base $\Lambda$CDM model.

As in other gravity models, the enhancement produced due to the fifth force is
reduced by the total neutrino mass ($M_{\nu}$). The bigger $M_{\nu}$, the
streaming effect grows, resulting in a more decreased PS. It has been shown in
HS models [13] and in Galileon [29] that one cannot distinguish these models
with a bigger neutrino from GR with a smaller neutrino mass. The effects are
not exactly cancelled out, but the distinction among models results
unpractical for the level of precision we have nowadays or in the near future;
DESI will be providing PS determinations within 1% [14]. Up to intermediate
wavenumbers ($k\gtrsim 0.1$ $h$/Mpc) deviations can be hidden with a selected
$M_{\nu}$ (left panel), that in our case we found $M_{\nu}=0.14$ eV, to
produce tiny differences of less than 1%. However, the bump at linear scales,
$k\sim 0.005$ Mpc/$h$ represents a big deviation (28% for $n=1$ and 25 % for
$n=3$) from $\Lambda$CDM; generic bumps in the PS have been recently studied
in ref. [41], where it is shown that a second smaller bump is induced at
nonlinear scales. As the universe evolves, clustering increases and at $z=0$
(right panel), at wavenumbers $k\gtrsim 0.1$ $h$/Mpc the differences are
bigger than 5% for $M_{\nu}=0.23$ eV, and they can be larger for smaller
neutrino masses. To again force the differences to tiny values (less than 1%)
at $z=0$ for these scales, taking the $n=3$ model (that is the one closest to
GR in our analysis), one would need a neutrino mass as big as $M_{\nu}=0.34$
eV, or $M_{\nu}=0.36$ eV for negative relative PS, values that are still
compatible to terrestrial neutrino mass measurements. We show this in figure
3, but the bump at $k\sim 0.01$ $h$/Mpc is still there and represents a big
deviation (15 %) from $\Lambda$CDM. Therefore, this model cannot get mixed up
with GR provided with lower masses.

Other gravity models such as HS [13, 19] or Galileon [29] for some high range
of neutrino masses are degenerate with GR assuming low masses at the level of
precision required for the matter PS. As demonstrated above, in $n$KGB models
a distinction is possible even for the matter PS. However, more features are
expected, and further ways to distinguish gravity models, when one considers
peculiar velocities. Then, apart from linear overdensity statistics,
$P_{\delta\delta}\equiv P_{L}(k)$, it is useful to consider the linear
density-velocity power spectrum,
$P_{\delta\theta}=\frac{f(k)}{f_{0}}P_{L}(k)$, and the velocity-velocity
spectrum, $P_{\theta\theta}=\left(\frac{f(k)}{f_{0}}\right)^{2}P_{L}(k)$,
where $f_{0}(z)=f(k=0,z)$ is the growth rate in $\Lambda$CDM that is
independent of the wavenumber, for further details we refer to [7]. In ref.
[19] is shown that the power spectra $P_{\delta\theta}$ and $P_{\theta\theta}$
are different for the HS model with massive neutrinos from GR (with or)
without massive neutrinos. In the same line of thought, RSD serve to sort out
these effects for HS models [21, 22]. The main ingredient lies in the growth
rate, that for our models depends on scale and time:

$f(k,z)\equiv\frac{d{\rm ln}D_{+}}{d{\rm
ln}a}=-\frac{1}{2}\,\frac{(1+z)}{P_{L}}\,\frac{dP_{L}}{dz}\,,$ (3.1)

where $D_{+}(k,z)$ is the growth function. Using eq. (3.1) we computed
$f(k,z)$ for $n=1$ and $n=3$ KGB models that are shown in figure 4 evaluated
at $z=1$ (left panel) and $z=0$ (right panel). At large scales (small $k$) the
growth rates of models with different neutrinos masses tend to the same value,
since there massive neutrinos behave simply as cold dark matter. At $z=1$, it
grows around 4% between $k=10^{-3}$ and $k=10^{-2}$ $h$/Mpc for the $n=3$
model and 8% for $n=1$. At $z=0$, changes are more evident, it grows around
35% between $k=10^{-3}$ and $k=10^{-2}$ $h$/Mpc for the $n=3$ model and 40%
for $n=1$. This behavior is expected since the fifth force, that enhances the
gravitational pull, increases the growth of structures, as well. Smaller
deviations of the growth rate are reported for the HS model [42] at $z=0$. We
also plot in figure 5 the resulting $P_{\delta\theta}$ and $P_{\theta\theta}$
with respect to $\Lambda$CDM power spectra, where the differences to the later
model are evident. The most important effect is in $P_{\theta\theta}$ since it
depends quadratically on $f(k,z)$ in comparison to $P_{\delta\theta}$ that
depends on it linearly. Overall the effect of neutrino mass in $f(k,z)$ is
much less important than $n$KGB gravitation.

Nonlinear deviations are computed in ref. [7] for the HS model with zero
neutrino mass, where it is shown that nonlinear terms tend to decrease
$P_{\delta\theta}$ and $P_{\theta\theta}$ in the quasilinear regime. Massive
neutrinos have been recently considered nonlinearly in RSD at 1-loop using
Eulerian perturbation theory for $\Lambda$CDM model [43], resulting in an
excellent agreement with simulation data up to $k\sim 0.25$ $h$/Mpc.

Figure 2: Relative change of the matter PS for $n$KGB models with respect to
the one for the $\Lambda$CDM model. The left panel is computed at redshift
$z=1$ while the right panel at $z=0$. Solid lines correspond to
$M_{\nu}=0.06\rm{eV}$, dash–dotted to $M_{\nu}=0.14\,\rm{eV}$, and dotted
lines to $M_{\nu}=0.23\,\rm{eV}$; colors denote different $n$KGB models.
Figure 3: Relative change of the matter PS for the $n=3$ KGB model with
respect to the one for the $\Lambda$CDM model at redshift $z=0$. Dotted lines
correspond to $M_{\nu}=0.23\rm{eV}$, while small-dashed are for
$M_{\nu}=0.34\,\rm{eV}$ that yields a 1% PS difference and long-dashed lines
correspond to $M_{\nu}=0.36\,\rm{eV}$ that yields a -1% PS difference. A bump
at $k\sim 0.01$ $h$/Mpc is the smoking gun of this model.

Figure 4: Scale–dependent growth rate for $n$KGB models and different sum of
neutrino masses. Left panel is computed at $z=1$ while the right panel is
evaluated at $z=0$.

Figure 5: Relative change of linear power spectra for $n$KGB models and
different sum of neutrino masses evaluated at $z=0$. The left panel shows the
cross-spectrum between the overdensity $\delta$ and the velocity divergence
$\theta$. The right panel is for the velocity divergence PS.

##  4 Fitting nKGB models to cosmological data

We will constrain the base cosmological parameters, the sum of neutrino masses
as well as the parameters of $n$KGB models. This MG model has in principle two
extra parameters, $g$ and $n$ in eq. (2.15). $g$ controls the strength of the
braiding and $n$ is the exponent of the kinetic coupling. In order to
disentangle both contributions, first we chose to explore a free amplitude
$g$, at fixed $n$–values. We consider separately models with $n=1,2,$ and $3$.
And we know that as $n\rightarrow\infty$, the models tend to $\Lambda$CDM.
Later, we will consider both $g$ and $n$ free parameters.

We use the Einstein–Boltzmann solver hi_class [40], which incorporates the
field equations presented in section 2, and the Markov chain Monte Carlo code
MontePython [44] to sample the parameters. For the statistical analysis we
employ the following data sets (Planck + BAO): Planck 2018
TT,TE,EE+lowE+lensing [2] and BAO BOSS DR12 [45]. We include Planck lensing
data, because there is a known dependency between neutrino mass and lensing,
as explained in point iii) of section 3; see however possible issues with
lensing data from Planck and other probes [46]. Since massive neutrinos are
thought to become nonrelativistic after last scattering, they change the
angular diameter distance to this surface, point i) of section 3, and
therefore the neutrino mass ($M_{\nu}$) is degenerate with $H_{0}$. To break
this parameter degeneracy, it is important to include BAO data (that is
measured at low z), to fit better the acoustic scale. On the other side, we
also use correlation function measurements of the KiDS collaboration [47, 48]
set alone to find out the level of the $\sigma_{8}$-tension between weak
lensing and Planck+BAO data in these MG models.As supernovae data have less
constraining power for neutrino masses, we have not included them here.

The contour confidence plots at $68\%$ and $95\%$ for the parameters of
interest are shown in figures 6, 7, and 8 for the $n=1,2,3$ KGB models,
respectively, and the best fit results are reported in table 1 at $95\%$
confidence level (CL.) for the varied parameters
$(\omega_{b},\omega_{c},H_{0},A_{s},n_{s},\tau,M_{\nu},g)$ and derived ones
$(\Omega_{m},\sigma_{8})$. Neutrinos masses are given in $eV$.

Figure 6: Left panel: Contour confidence plots at $68\%$ and $95\%$ CL. of
$\sigma_{8}$ vs $\Omega_{m}$ for the $n=1$ KGB model using Planck+ BAO and
KiDS data; in this model the $\sigma_{8}$ tension between these data is around
$0.28\sigma$, so it essentially disappears. Right panel: $M_{\nu}$ [$eV$] vs
$g$ for the Planck+ BAO data; we did not include KiDS data since this dataset
has much less constraining power for these parameters.

Figure 7: The same contour confidence plots as in figure 6, for the $n=2$ KGB
model.

Figure 8: The same contour confidence plots as in figure 6, for the $n=3$ KGB
model.

It is known that in $\Lambda$CDM an increasing neutrino mass causes a
decreasing value of the Hubble constant, and hence it worsens the tension with
the distance-ladder determination [49]. Some MG models, such as Galileons,
work in the other direction, ameliorating this tension [50]. Comparing the
$H_{0}$ best fit results from table 1 for Planck+BAO data with the reported
value [49] $H_{0}=73.2\pm 1.3\,\rm{km(sMpc)^{-1}}$ we found that the tension
is reduced for these $n$ values in comparison to $\Lambda$CDM. In fact, there
is about $0.43\sigma$ ($n=1$), $1.72\sigma$ ($n=2$), and $2.39\sigma$ ($n=3$)
tension, from which we conclude that the $n=1$ model shows no tension with
distance-ladder values, but bigger $n$’s increase the tension. When a $n$KGB
is in its attractor solution, the asymptotic equation of state for the scalar
field is in the phantom domain ($w_{\phi}<-1$) during the matter era, and
tends to -1 in the scalar field dominated era [33]. This is reason why this
model solves the Hubble tension, as known for phantom Universes [2]. Note that
the KiDS data does not sufficiently constrain the Hubble constant, similar as
in $\Lambda$CDM model [48]. Another aspect is that this model predicts growing
gravitational potentials at low redshift that can be tested using ISW data
cross correlated with clustering [8]. This has been done in [9], where a
conflict to fit data is shown in this model. Another property of general
$n$KGB models is that they exclude the null neutrino mass, see right panels of
figures 6, 7, and 8. These results are in agreement with current neutrino mass
constraints since all the confidence contours (95% CL.) for this parameter are
larger than the lower limit of 60 meV from neutrino oscillations [39] and are
below the upper bound $M_{\nu}<0.8$ eV (90 % CL.) from the tritium decay in
KATRIN [51]. If the mass scale happens to be measured closer to the lower
bound, small $n$ models could be in trouble to explain this small mass. A
normal neutrino hierarchy would put stringer limits than an inverted one.

Another important aspect is to check how the so-called $\sigma_{8}$ tension
behaves in $n$KGB models. It has been previously reported that the $n=1$ model
tends to rise the value of $\sigma_{8}$ [50], but the effect of massive
neutrinos is to lower it. The results when these two modifications are
considered at once are shown in table 1, where we observe no tension for the
considered $n$ values, using Planck+BAO and KiDS data. In small $n$ models
there is no $\sigma_{8}$ tension. This effect is also evident in the left
panels of figures 6, 7, and 8, where we observe that the values from our
Planck+BAO selection are compatible at 1$\sigma$ with the results from the
KiDS collaboration. The results in this work are self consistent when
considering both datasets and the same $n$ value, and also are consistent with
$\Lambda$CDM results from Planck 2015 [47] and from Planck 2018 [2]
collaborations.

To quantify how viable the considered $n$KGB models are compared to
$\Lambda$CDM we use the Akaike Information Criterion (AIC)[52, 53]; recall
that a model with a $\Delta$AIC bigger than 10 is essentially not supported,
and values approaching to zero are strongly supported. Our results in table 1
indicate that models with $n=1,2,3$ are supported by KiDS data, however $n=1$
is discarded by Planck+BAO data; $n=2,3$ values represent models that, while
are not best supported, are not discarded.

Table 1: Constraints at 95% confidence level on the cosmological and $n$KGB parameter models with $n=1,2,3$ and $M_{\nu}$ $[eV]$ as free parameter. To compare the results with those from $\Lambda$CDM we include the $\Delta$AIC values defined as $\Delta\rm{AIC}=\rm{AIC}_{nKGB}-\rm{AIC}_{\Lambda\rm{CDM}}$.  | $n=1$ | $n=2$ | $n=3$
---|---|---|---
|  Planck+BAO |  KiDS |  Planck+BAO |  KiDS |  Planck+BAO |  KiDS
Parameter |  95% limits |  95% limits |  95% limits |  95% limits |  95% limits |  95% limits
$100\leavevmode\nobreak\ \omega_{b}$ |  $2.223^{+0.026}_{-0.026}$ |  $2.23^{+3.6}_{-3.4}$ |  $2.240^{+0.027}_{-0.026}$ |  $2.23^{+0.37}_{-0.34}$ |  $2.246^{+0.017}_{-0.016}$ |  $2.23^{+0.36}_{-0.34}$
$\omega{}_{c}$ |  $0.1204^{+0.0019}_{-0.0019}$ |  $0.132^{+0.080}_{-0.078}$ |  $0.1185^{+0.0020}_{-0.0019}$ |  $0.139^{+0.082}_{-0.081}$ |  $0.1181^{+0.0015}_{-0.0017}$ |  $0.143^{+0.086}_{-0.084}$
$H_{0}$ |  $72.8^{+1.5}_{-1.5}$ |  $75.5^{+6.6}_{-9.4}$ |  $69.9^{+1.4}_{-1.4}$ |  $75.4^{+6.6}_{-9.5}$ |  $68.8^{+1.2}_{-1.2}$ |  $75.4^{+6.7}_{-9.7}$
$10^{9}A_{s}$ |  $2.066^{+0.059}_{-0.064}$ |  $2.8^{+1.3}_{-1.1}$ |  $2.044^{+0.067}_{-0.069}$ |  $2.7^{+1.2}_{-1.0}$ |  $2.033^{+0.067}_{-0.075}$ |  $2.7^{+1.2}_{-1.0}$
$n_{s}$ |  $0.9604^{+0.0073}_{-0.0073}$ |  $1.16^{+0.14}_{-0.19}$ |  $0.9658^{+0.0075}_{-0.0075}$ |  $1.15^{+0.15}_{-0.20}$ |  $0.9669^{+0.0054}_{-0.0054}$ |  $1.15^{+0.15}_{-0.20}$
$\tau{}_{reio}$ |  $0.048^{+0.015}_{-0.016}$ |  – |  $0.045^{+0.017}_{-0.017}$ |  – |  $0.042^{+0.015}_{-0.017}$ |  –
$M_{\nu}$ |  $0.46^{+0.13}_{-0.13}$ |  $<2.27$ |  $0.33^{+0.13}_{-0.13}$ |  $<2.33$ |  $0.31^{+0.12}_{-0.13}$ |  $<2.37$
$g$ |  $1.4^{+2.7}_{-1.2}$ |  $3.2^{+4.5}_{-3.2}$ |  $0.24^{+0.73}_{-0.26}$ |  $12^{+14}_{-12}$ |  $0.9999^{+0.0021}_{-0.0021}$ |  $15^{+17}_{-15}$
$\Omega_{m}$ |  $0.278^{+0.015}_{-0.014}$ |  $0.29^{+0.15}_{-0.14}$ |  $0.296^{+0.016}_{-0.015}$ |  $0.31^{+0.15}_{-0.14}$ |  $0.304^{+0.013}_{-0.013}$ |  $0.31^{+0.15}_{-0.15}$
$\sigma_{8}$ |  $0.817^{+0.033}_{-0.033}$ |  $0.75^{+0.23}_{-0.20}$ |  $0.802^{+0.030}_{-0.031}$ |  $0.75^{+0.23}_{-0.19}$ |  $0.794^{+0.028}_{-0.030}$ |  $0.74^{+0.23}_{-0.19}$
$\Delta$AIC |  $15.4$ |  $3.4$ |  $7.7$ |  $3.8$ |  $8.8$ |  $3.7$

In section 2, we set the conditions for the dynamical system to be in the
attractor solution, following eqs. (2.18, 2.19). This condition in turn
implies a specific determination of our constant $g$ in terms of $\Omega_{i}$
and $H_{0}$ [33]. We checked where in our parameter space this solution is,
founding that for $n=1$ and $n=2$ KGB models, the attractor solution lies
within the 1-$\sigma$ of our resulting parameters; for $n=3$ the best fit is
out of the attractor solution, similar as in other models studied in ref.
[36].

Finally, we have performed a fit to the $n$KGB models letting both $n$ and
$g$, together with the neutrino mass, to be free parameters. We use Planck+BAO
data and assumed a prior $0.5\leq n\leq 30$, the lowest limit stemming from
the physical requirement to avoid perturbation instabilities and the upper
limit was found to be big enough, after performing a few runs. We found the
$n$ is bounded from above at 95% CL. by $22.5$. The contours plots are
depicted in figure 9 and the best fits for the other cosmological parameters
at 95% CL. are: $100\omega_{b}=2.247^{+0.029}_{-0.030}$,
$\omega_{c}=0.1174^{+0.0023}_{-0.0023}$, $H_{0}=68.5^{+1.8}_{-1.8}$,
$10^{9}A_{s}=2.030^{+0.081}_{-0.10}$, $n_{s}=0.9688^{+0.0078}_{-0.0082}$,
$\tau_{reio}=0.042^{+0.021}_{-0.027}$, $M_{\nu}=0.27^{+0.15}_{-0.14}$,
$g<16.6$, $n<22.5$, $\Omega_{m}=0.305^{+0.019}_{-0.018}$,
$\sigma_{8}=0.794^{+0.028}_{-0.030}$. The $\Delta\rm{AIC}$ for this model is
10.2. Its likelihood is a little smaller than the cases in which $n$ is a
fixed parameter, but the model is penalized by the extra parameter, so $n$ as
a free parameter is strongly disfavored. When comparing its neutrino mass with
the one obtained in $\Lambda$CDM for the same data set, given by
$M_{\nu}<0.12$ within 95% CL. [2], this model results in a different interval
that points to higher masses.

Figure 9: Contour confidence plots using Planck+BAO data for $n,g$, and
$M_{\nu}$ at 68 and 95% CL. In these plots $n$ was considered as a free
parameter.

##  5 Conclusions

There is an increasing interest in the scientific community to test gravity at
large scales. In this line of thought one considers models beyond GR, that in
turn imply a larger parameter space. We investigated the cosmological effects
of braiding, $n$KGB models, that possess two extra two parameters $(n,g)$ and
added effects of massive neutrinos to study possible parameter degeneracies.

We first consider the CMB TT power spectrum, founding deviations from
$\Lambda$CDM results mainly at large angles, that are difficult to
discriminate due to large statistical uncertainties of the cosmic variance;
however, cross correlating CMB temperature fluctuations with clustering is an
interesting way to test these theories, as shown in refs. [8, 9, 10].

The fifth force associated to the scalar field induces an extra gravitational
pull that changes the growth history at late times, and this affects the
matter PS. We find that in braiding cosmologies the effect is to enhance PS in
a way that already at linear scales it is possible to distinguish $n$KGB
models from GR. Even in the case that one increases neutrino masses in $n$KGB
models, to lower the PS, this changes in a way that is not possible to mimic
the GR PS with lighter neutrinos masses, as it is is evident in figures 2 and
3. This distinction is not possible to achieve in HS gravity for the level of
precision required, as shown in refs. [13, 19]. It is possible however to
compute the effect of peculiar velocities in the PS that have shown to add
features in MG two-point observables $P_{\delta\theta}$, $P_{\theta\theta}$
[19], or RSD [21, 22]. The key ingredient that distinguishes these statistics
is the growth rate, that for MG models depend on both scale and redshift. We
have computed $f(k,z)$ and added neutrino masses to compare its effects. We
found that the effects of neutrino masses are smaller than the gravitational
effect due to the fifth force of $n$KGB models. In general, the effects of MG
are more evident at lower redshifts, since the growth function has more time
to rise at scales where MG imprints special features. We found that at large
scales (small $k$) the $f(k,z)$ value tends to the $\Lambda$CDM value but at
smaller, still linear scales, the growth rate is larger for $n$KGB models. At
$z=1$, it grows around 4% between $k=10^{-3}$ and $k=10^{-2}$ $h$/Mpc for the
$n=3$ model and 8% for $n=1$. At $z=0$, the changes are bigger, it grows
around 35% between $k=10^{-3}$ and $k=10^{-2}$ $h$/Mpc for the $n=3$ model and
70% for $n=1$. This behavior is expected since the fifth force, that enhances
the gravitational pull, increases the growth of structures, as well. Smaller
deviations of the growth rate are reported for the HS model [42] at $z=0$. For
completeness we also computed $P_{\delta\theta}$ and $P_{\theta\theta}$ and
compared our results to $\Lambda$CDM model, where the differences to the later
model are evident. The most important effect is in $P_{\theta\theta}$ since it
depends quadratically on $f(k,z)$ in comparison to $P_{\delta\theta}$ that
depends on it linearly, as can be appreciated in figure 5. Overall the effect
of neutrino mass in $f(k,z)$ is much less important than the $n$KGB
gravitation.

We constrained the base cosmological parameters, the sum of neutrino masses as
well as the parameters of $n$KGB models. We fitted the models to our selected
Planck+BAO data, and apart to KiDS correlation function data. Best fitted
values are displayed in table 1, showing that: i) the n=1 KGB model reports no
$H_{0}$ tension using Planck+BAO data, whereas bigger $n$ tends to
$\Lambda$CDM and therefore the tension exists. The reason is that KGB models
have a phantom behavior and therefore, as known, these models tend to fix this
tension [2]; ii) All considered $n$KGB models show no $\sigma_{8}$-tension,
when employing Planck+BAO and KiDS data. The reason for this is that although
$n$KGB models increase the PS, as $n$ is bigger, the best fitted values of
neutrino mass go to lower values, that help to relax the tension.

Though some of these models have nice properties such as alleviating the
$H_{0}$ and $\sigma_{8}$ tensions, all predict a bump in the linear matter PS,
as seen in figures 2 and 3, that is also enhanced through $f(k,z)$, and that
surpasses by a significant factor the results of $\Lambda$CDM, cf. plots 5.
Since we know from present data that possible deviations of power spectra are
less than few percent from that of $\Lambda$CDM, see e.g. [54]., then $n$KGB
models are virtually excluded as realistic MG models.

## Acknowledgments

The authors wish to thank A. Aviles for discussing some of the results of this
work. The authors acknowledge partial support by Conacyt project 283151 and
SNI. U.N and G.G.A are grateful to FORDECYT PRONACES-CONACYT for support of
the present research under Grant CF-MG-2558591 and CF-140630-UNAM-UMSNH. U.N.
thanks the Programa para el Desarrollo Profesional Docente of the Secretaría
de Educación Pública (PRODEP-SEP) of México and the Coordinación de la
Investigación Científica of the Universidad Michoacana de San Nicolás de
Hidalgo (CIC-UMSNH).

## References

  * [1] G. Bertone, D. Hooper and J. Silk, _Particle dark matter: Evidence, candidates and constraints_ , _Phys. Rept._ 405 (2005) 279–390, [hep-ph/0404175].
  * [2] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * [3] Z. Zhai et al., _The Clustering of Luminous Red Galaxies at z $\sim$ 0.7 from EBOSS and BOSS Data_, _Astrophys. J._ 848 (2017) 76, [1607.05383].
  * [4] Supernova Search Team collaboration, A. G. Riess et al., _Observational evidence from supernovae for an accelerating universe and a cosmological constant_ , _Astron. J._ 116 (1998) 1009–1038, [astro-ph/9805201].
  * [5] Supernova Cosmology Project collaboration, S. Perlmutter et al., _Measurements of Omega and Lambda from 42 high redshift supernovae_ , _Astrophys. J._ 517 (1999) 565–586, [astro-ph/9812133].
  * [6] K. Koyama, _Cosmological Tests of Modified Gravity_ , _Rept. Prog. Phys._ 79 (2016) 046902, [1504.04623].
  * [7] A. Aviles, G. Valogiannis, M. A. Rodriguez-Meza, J. L. Cervantes-Cota, B. Li and R. Bean, _Redshift space power spectrum beyond Einstein-de Sitter kernels_ , _JCAP_ 04 (2021) 039, [2012.05077].
  * [8] S. Ferraro, B. D. Sherwin and D. N. Spergel, _WISE measurement of the integrated Sachs-Wolfe effect_ , _Phys. Rev. D_ 91 (2015) 083533, [1401.1193].
  * [9] J. Renk, M. Zumalacárregui, F. Montanari and A. Barreira, _Galileon gravity in light of ISW, CMB, BAO and H 0 data_, _JCAP_ 10 (2017) 020, [1707.02263].
  * [10] J. A. Kable, G. Benevento, N. Frusciante, A. De Felice and S. Tsujikawa, _Probing Modified Gravity with Integrated Sachs-Wolfe CMB and Galaxy Cross-correlations_ , 2111.10432.
  * [11] S. Alam et al., _Towards testing the theory of gravity with DESI: summary statistics, model predictions and future simulation requirements_ , _JCAP_ 11 (2021) 050, [2011.05771].
  * [12] J. Lesgourgues and S. Pastor, _Massive neutrinos and cosmology_ , _Phys. Rept._ 429 (2006) 307–379, [astro-ph/0603494].
  * [13] N. Bellomo, E. Bellini, B. Hu, R. Jimenez, C. Pena-Garay and L. Verde, _Hiding neutrino mass in modified gravity cosmologies_ , _JCAP_ 02 (2017) 043, [1612.02598].
  * [14] DESI collaboration, A. Aghamousa et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , 1611.00036.
  * [15] L. Amendola et al., _Cosmology and fundamental physics with the Euclid satellite_ , _Living Rev. Rel._ 21 (2018) 2, [1606.00180].
  * [16] LSST Dark Energy Science collaboration, A. Abate et al., _Large Synoptic Survey Telescope: Dark Energy Science Collaboration_ , 1211.0310.
  * [17] A. Hojjati, L. Pogosian and G.-B. Zhao, _Testing gravity with CAMB and CosmoMC_ , _JCAP_ 08 (2011) 005, [1106.4543].
  * [18] W. Hu and I. Sawicki, _Models of f(R) Cosmic Acceleration that Evade Solar-System Tests_ , _Phys. Rev. D_ 76 (2007) 064004, [0705.1158].
  * [19] S. Hagstotz, M. Gronke, D. Mota and M. Baldi, _Breaking cosmic degeneracies: Disentangling neutrinos and modified gravity with kinematic information_ , _Astron. Astrophys._ 629 (2019) A46, [1902.01868].
  * [20] A. Peel, V. Pettorino, C. Giocoli, J.-L. Starck and M. Baldi, _Breaking degeneracies in modified gravity with higher (than 2nd) order weak-lensing statistics_ , _Astron. Astrophys._ 619 (2018) A38, [1805.05146].
  * [21] B. S. Wright, K. Koyama, H. A. Winther and G.-B. Zhao, _Investigating the degeneracy between modified gravity and massive neutrinos with redshift-space distortions_ , _JCAP_ 06 (2019) 040, [1902.10692].
  * [22] J. E. García-Farieta, F. Marulli, A. Veropalumbo, L. Moscardini, R. A. Casas-Miranda, C. Giocoli et al., _Clustering and redshift-space distortions in modified gravity models with massive neutrinos_ , _Mon. Not. Roy. Astron. Soc._ 488 (2019) 1987–2000, [1903.00154].
  * [23] M. Chevallier and D. Polarski, _Accelerating universes with scaling dark matter_ , _Int. J. Mod. Phys. D_ 10 (2001) 213–224, [gr-qc/0009008].
  * [24] E. V. Linder, _Exploring the expansion history of the universe_ , _Phys. Rev. Lett._ 90 (2003) 091301, [astro-ph/0208512].
  * [25] B. Bose, B. S. Wright, M. Cataneo, A. Pourtsidou, C. Giocoli, L. Lombriser et al., _On the road to per cent accuracy – V. The non-linear power spectrum beyond $\Lambda$CDM with massive neutrinos and baryonic feedback_, _Mon. Not. Roy. Astron. Soc._ 508 (2021) 2479–2491, [2105.12114].
  * [26] M. Zhang, J.-F. Zhang and X. Zhang, _Impacts of dark energy on constraining neutrino mass after Planck 2018_ , _Commun. Theor. Phys._ 72 (2020) 125402, [2005.04647].
  * [27] M. Baldi, F. Villaescusa-Navarro, M. Viel, E. Puchwein, V. Springel and L. Moscardini, _Cosmic degeneracies – I. Joint N-body simulations of modified gravity and massive neutrinos_ , _Mon. Not. Roy. Astron. Soc._ 440 (2014) 75–88, [1311.2588].
  * [28] S. Contarini, F. Marulli, L. Moscardini, A. Veropalumbo, C. Giocoli and M. Baldi, _Cosmic voids in modified gravity models with massive neutrinos_ , 2009.03309.
  * [29] A. Barreira, B. Li, C. Baugh and S. Pascoli, _Modified gravity with massive neutrinos as a testable alternative cosmological model_ , _Phys. Rev. D_ 90 (2014) 023528, [1404.1365].
  * [30] T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, _Modified Gravity and Cosmology_ , _Phys. Rept._ 513 (2012) 1–189, [1106.2476].
  * [31] C. Deffayet, O. Pujolas, I. Sawicki and A. Vikman, _Imperfect Dark Energy from Kinetic Gravity Braiding_ , _JCAP_ 10 (2010) 026, [1008.0048].
  * [32] T. Kobayashi, M. Yamaguchi and J. Yokoyama, _G-inflation: Inflation driven by the Galileon field_ , _Phys. Rev. Lett._ 105 (2010) 231302, [1008.0603].
  * [33] R. Kimura and K. Yamamoto, _Large Scale Structures in Kinetic Gravity Braiding Model That Can Be Unbraided_ , _JCAP_ 04 (2011) 025, [1011.2006].
  * [34] R. Arjona, W. Cardona and S. Nesseris, _Designing Horndeski and the effective fluid approach_ , _Phys. Rev. D_ 100 (2019) 063526, [1904.06294].
  * [35] R. Kimura, T. Kobayashi and K. Yamamoto, _Observational Constraints on Kinetic Gravity Braiding from the Integrated Sachs-Wolfe Effect_ , _Phys. Rev. D_ 85 (2012) 123503, [1110.3598].
  * [36] W. Cardona, R. Arjona, A. Estrada and S. Nesseris, _Cosmological constraints with the Effective Fluid approach for Modified Gravity_ , _JCAP_ 05 (2021) 064, [2012.05282].
  * [37] E. Bellini, I. Sawicki and M. Zumalacárregui, _hi_class: Background Evolution, Initial Conditions and Approximation Schemes_ , _JCAP_ 02 (2020) 008, [1909.01828].
  * [38] T. Kobayashi, _Horndeski theory and beyond: a review_ , _Rept. Prog. Phys._ 82 (2019) 086901, [1901.07183].
  * [39] Particle Data Group collaboration, P. A. Zyla et al., _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * [40] M. Zumalacárregui, E. Bellini, I. Sawicki, J. Lesgourgues and P. G. Ferreira, _hi_class: Horndeski in the Cosmic Linear Anisotropy Solving System_ , _JCAP_ 08 (2017) 019, [1605.06102].
  * [41] D. V. Gomez-Navarro, A. Mead, A. Aviles and A. de la Macorra, _Impact of cosmological signatures in two-point statistics beyond the linear regime_ , _Mon. Not. Roy. Astron. Soc._ 504 (2021) 3284–3297, [2009.12717].
  * [42] Y.-H. Li, J.-F. Zhang and X. Zhang, _Probing $f(R)$ cosmology with sterile neutrinos via measurements of scale-dependent growth rate of structure_, _Phys. Lett. B_ 744 (2015) 213–217, [1502.01136].
  * [43] A. Aviles, A. Banerjee, G. Niz and Z. Slepian, _Clustering in massive neutrino cosmologies via Eulerian Perturbation Theory_ , _JCAP_ 11 (2021) 028, [2106.13771].
  * [44] T. Brinckmann and J. Lesgourgues, _MontePython 3: boosted MCMC sampler and other features_ , _Phys. Dark Univ._ 24 (2019) 100260, [1804.07261].
  * [45] BOSS collaboration, S. 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 (2017) 2617–2652, [1607.03155].
  * [46] E. Di Valentino and A. Melchiorri, _Neutrino Mass Bounds in the era of Tension Cosmology_ , 2112.02993.
  * [47] H. Hildebrandt et al., _KiDS-450: Cosmological parameter constraints from tomographic weak gravitational lensing_ , _Mon. Not. Roy. Astron. Soc._ 465 (2017) 1454, [1606.05338].
  * [48] F. Köhlinger et al., _KiDS-450: The tomographic weak lensing power spectrum and constraints on cosmological parameters_ , _Mon. Not. Roy. Astron. Soc._ 471 (2017) 4412–4435, [1706.02892].
  * [49] A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn et al., _Cosmic Distances Calibrated to 1% Precision with Gaia EDR3 Parallaxes and Hubble Space Telescope Photometry of 75 Milky Way Cepheids Confirm Tension with $\Lambda$CDM_, _Astrophys. J. Lett._ 908 (2021) L6, [2012.08534].
  * [50] A. Barreira, B. Li, C. Baugh and S. Pascoli, _The observational status of Galileon gravity after Planck_ , _JCAP_ 08 (2014) 059, [1406.0485].
  * [51] M. Aker et al., _First direct neutrino-mass measurement with sub-eV sensitivity_ , 2105.08533.
  * [52] H. Akaike, _A new look at the statistical model identification_ , _IEEE Transactions on Automatic Control_ 19 (1974) 716–723.
  * [53] K. P. Burnham and D. R. Anderson, _Multimodel inference: Understanding aic and bic in model selection_ , _Sociological Methods & Research_ 33 (2004) 261–304, [https://doi.org/10.1177/0049124104268644].
  * [54] Y. Kobayashi, T. Nishimichi, M. Takada and H. Miyatake, _Full-shape cosmology analysis of SDSS-III BOSS galaxy power spectrum using emulator-based halo model: a $5\%$ determination of $\sigma_{8}$_, 006110.06969.