#  $f(Q)$-gravity and neutrino physics

Luís Atayde<EMAIL_ADDRESS>Instituto de Astrofisíca e Ciências do
Espaço, Faculdade de Ciências da Universidade de Lisboa, Edificio C8, Campo
Grande, P-1749016, Lisboa, Portugal  Noemi Frusciante
<EMAIL_ADDRESS>Dipartimento di Fisica “E. Pancini”, Università
degli Studi di Napoli “Federico II”, Compl. Univ. di Monte S. Angelo, Edificio
G, Via Cinthia, I-80126, Napoli, Italy

(September 3, 2024)

###### Abstract

Within the $f(Q)$-gravity framework we perform a phenomenological study of the
cosmological observables in light of the degeneracy between neutrinos physics
and the modified gravity parameter and we identify specific patterns which
allow to break such degeneracy. We also provide separately constraints on the
total mass of the neutrinos, $\Sigma m_{\nu}$, and on the effective number of
neutrino species, $N_{\rm eff}$, using cosmic microwave background (CMB),
baryon acoustic oscillation (BAO), redshift space distortion (RSD), supernovae
(SNIa), galaxy clustering (GC) and weak gravitational lensing (WL)
measurements. The strongest upper bound on the total mass of the neutrinos is
found for the combination of CMB+BAO+RSD+SNIa and it is $\Sigma m_{\nu}<0.277$
eV at 95% C.L. For the same combination of data we find $N_{\rm
eff}=2.93^{+0.31}_{-0.34}$ at 95% C.L. We also find that all combinations of
data we consider, prefer a stronger gravitational interaction than
$\Lambda$CDM. Finally, we consider the $\chi^{2}$ and deviance information
criterion statistics and find the $f(Q)+\Sigma m_{\nu}$ model to be
statistically supported by data over the standard scenario. On the contrary
$f(Q)+N_{\rm eff}$ is supported by CMB+BAO+RSD+SNIa but a moderate evidence
against it is found with GC and WL data.

##  I Introduction

Extension to the symmetric teleparallel gravity Nester and Yo (1999); Adak
_et al._ (2013); Adak (2018) recently got attention from the cosmology
community as a possibility to explore new physics beyond the standard
cosmological model ($\Lambda$CDM). In this alternative theory, the
gravitational interaction is attributed to the non-metricity scalar $Q$. While
an action constructed only with this scalar is equivalent to the one of
General Relativity (GR) in flat space Beltrán Jiménez  _et al._ (2019), the
action built with a general function of $Q$, i.e. $f(Q)$ Beltrán Jiménez  _et
al._ (2018); Harko  _et al._ (2018); Xu  _et al._ (2019); Beltrán Jiménez  _et
al._ (2020a); Järv  _et al._ (2018); Rünkla and Vilson (2018); Beltrán Jiménez
_et al._ (2019) can generate a gravitational interaction at cosmological
scales which shows new interesting patterns both at the level of the
background evolution and in the propagation of scalar and tensor perturbations
Beltrán Jiménez  _et al._ (2020a); Jiménez  _et al._ (2019); Lazkoz  _et al._
(2019); Ayuso  _et al._ (2021); Khyllep  _et al._ (2021); Frusciante (2021).
Furthermore, for specific forms of the $f(Q)$ function it has been shown that
it can alleviate the $\sigma_{8}$ tension Barros  _et al._ (2020), while
others allow for a better fit to cosmological data Anagnostopoulos  _et al._
(2021); Atayde and Frusciante (2021a). This theory is then challenging
$\Lambda$CDM from different perspectives.

When exploring modified gravity (MG) theories, such as $f(Q)$-gravity, it is
important to consider that MG effects on cosmological observables are highly
degenerate with neutrino physics Barreira  _et al._ (2014); Shim  _et al._
(2014); Baldi  _et al._ (2014); He (2013); Dossett  _et al._ (2014); Hojjati
_et al._ (2011); Motohashi  _et al._ (2013); Hu  _et al._ (2015); Bellomo  _et
al._ (2017); Frusciante and Benetti (2021); Ballardini  _et al._ (2020).
Indeed, to name a few, a modified gravitational interaction can, for instance,
impact the background expansion history, shape the temperature-temperature
cosmic microwave background radiation (CMB) power spectrum at all multipoles
($\ell$) through a modified early and late integrated Sachs-Wolfe (ISW) effect
Sachs and Wolfe (1967); Kofman and Starobinsky (1985) or a change of the
amplitude and position of the high-$\ell$ peaks. These effects at high-$\ell$
are due to a modified background expansion history Hu and White (1996) or to a
coupling with dark matter Amendola  _et al._ (2012). MG can also impact the
lensing potential Acquaviva and Baccigalupi (2006) and change the growth of
structures Peebles (1984); Barrow and Saich (1993) with a direct effect on the
lensing and matter power spectra. Similarly, neutrinos leave detectable
signatures on cosmological observations that can be used to constrain their
properties, such as the effective number of neutrino species, $N_{\rm eff}$,
and the total neutrino mass, $\Sigma m_{\nu}$ Lattanzi and Gerbino (2018);
Abazajian  _et al._ (2015). Increasing $N_{\rm eff}$ results in a faster
expansion at earlier times altering the radiation to matter equality which in
turn enhances the first acoustic peak due to the early ISW effect, while the
other CMB peaks are moved to higher multipoles. Additionally it also changes
the scale of Silk damping which has the effect of lowering the damping tail of
the CMB spectrum and finally it suppresses the matter power spectrum.
Increasing the sum of the neutrino masses can decrease the low-$\ell$ tail of
the CMB TT power spectrum due to the late ISW effect, change the position of
the high-$\ell$ CMB peaks, reduce the weak lensing effect and dump the growth
of structures on small scales. Therefore, constraints on neutrinos properties
are highly dependent on the kind of gravitational interaction one considers.
Within the standard cosmological model the upper bound on the sum of the
neutrino masses at 95% for the combination of Planck15+lensing+BAO+JLA+$H_{0}$
is $\Sigma m_{\nu}<0.23$ eV Ade  _et al._ (2016a) and the determination on the
effective number of neutrinos is $N_{\rm eff}=3.15\pm 0.23$ at 95% using
Planck15+BAO Ade  _et al._ (2016a), the latter being compatible with the
prediction of the standard neutrino decoupling model, $N_{\rm eff}=3.046$.
Bounds on these two parameters have been found in the MG context as well, for
example, when $N_{\rm eff}$ is fixed to its standard value, one has for Hořava
gravity $\Sigma m_{\nu}<0.165$ eV (95%, Planck18+lensing+DES+JLA+BKP15)
Frusciante and Benetti (2021), for designer $f(R)$ model $\Sigma m_{\nu}<0.32$
eV (95%, Planck13 +WMAP+BAO+lensing+WiggleZ) Hu  _et al._ (2015), for the same
data set when considering a simple effective field theory of MG model, one has
$\Sigma m_{\nu}<0.26$ eV Hu  _et al._ (2015); for the generalized cubic
covariant galileon model instead one is able to find a lower bound $\Sigma
m_{\nu}>0.11$ eV at 1$\sigma$ using Planck15+BAO+RSD+JLA Frusciante  _et al._
(2020); for the simplest model of scalar-tensor theories with a conformally
coupled scalar field it is $\Sigma m_{\nu}<0.13$ (95%, Planck18+BAO+HST),
while if we vary $N_{\rm eff}$ one finds $3.16\pm 0.19$ at 68% and when
varying them simultaneously one has $\Sigma m_{\nu}<0.14$ (95%) and $N_{\rm
eff}=3.14\pm 0.20$ (68%) for the same combinations of data Ballardini  _et
al._ (2020). Finally we mention again the designer $f(R)$ when varying both of
them simultaneously, one finds $N_{\rm eff}=3.58^{+0.72}_{-0.69}$ and $\Sigma
m_{\nu}<0.860$ eV (95%, Planck13+WMAP+BAO+ACT+SPT) He (2013). One can also
notice that when providing bounds on neutrino physics in the presence of
modifications of gravity the constraints might be weaker.

In this paper we consider a specific form of the $f(Q)$ function such that the
background evolution matches exactly the one of $\Lambda$CDM Beltrán Jiménez
_et al._ (2020b). This model is characterized by one extra free parameter and
modification to the gravitational interaction can be appreciated only at
perturbation level Frusciante (2021). Cosmological constraints on the modified
gravity parameter and on the cosmological parameters have been already derived
considering massive neutrinos to have a fixed mass ($\Sigma m_{\nu}=0.06$ eV,
the minimum value allowed by oscillation experiments) and $N_{\rm eff}=3.046$
Atayde and Frusciante (2021b). In this work we open the neutrino sector by
considering the mass of the neutrinos or their effective number as additional
parameters. We will then be able to investigate whether a degeneracy exists
with the modified gravity parameter and if current cosmological data can break
such degeneracy. We will also provide for the first time constraints on
$N_{\rm eff}$ and $\Sigma m_{\nu}$ in the context of $f(Q)$-gravity.

This paper is organized as follows. In Section II we review the theory behind
the $f(Q)$-gravity and we present the specific model we analyze. In particular
we describe how linear perturbation theory is treated in this analysis. In
Sec. III we describe the methodology, the codes and the data sets we adopt to
study the phenomenology associated with the degeneracy between a modified
gravitational interaction and the neutrino sector and to derive the
cosmological constraints for model and cosmological parameters, which are
thoroughly presented and discussed respectively in Secs. IV and V. We then
provide a model selection analysis in Sec. VI and finally we conclude in Sec.
VII.

##  II The model

We adopt the Palatini formalism which is characterized by the metric,
$g_{\mu\nu}$ and the connection $\Gamma^{\alpha}_{\phantom{X}\mu\nu}$ to be
independent fields. The action of the $f(Q)$ theory has the form Beltrán
Jiménez  _et al._ (2020a):

$S=\int d^{4}x\sqrt{-g}\left\\{-\frac{1}{16\pi
G_{N}}\left[Q+f(Q)\right]+\mathcal{L}_{m}(g_{\mu\nu},\chi_{i})\right\\}\,,$
(1)

where $g$ is the determinant of the metric $g_{\mu\nu}$, $G_{N}$ is the
Newtonian constant and $Q$ is the non-metricity scalar which is defined as
$Q=-Q_{\alpha\mu\nu}P^{\alpha\mu\nu}$ where $Q_{\alpha\mu\nu}$ is the non-
metricity tensor, which reads $Q_{\alpha\mu\nu}=\nabla_{\alpha}g_{\mu\nu}$ and
$P^{\alpha}_{\phantom{\alpha}\mu\nu}=-L^{\alpha}_{\phantom{\alpha}\mu\nu}/2+\left(Q^{\alpha}-\tilde{Q}^{\alpha}\right)g_{\mu\nu}/4-\delta^{\alpha}_{(\mu}Q_{\nu)}/4$is
the non-metricity coniugate, where we further define
$Q_{\alpha}=g^{\mu\nu}Q_{\alpha\mu\nu}$,
$\tilde{Q}_{\alpha}=g^{\mu\nu}Q_{\mu\alpha\nu}$ and
$L^{\alpha}_{\phantom{\alpha}\mu\nu}=(Q^{\alpha}_{\phantom{\alpha}\mu\nu}-Q_{(\mu\nu)}^{\phantom{(\mu\nu)}\alpha})/2$.
The deviation from GR is encoded in the general function of the non-metricity
scalar, i.e. $f(Q)$. Indeed the action (1) is equivalent to GR for $f(Q)=0$
Jiménez  _et al._ (2019). Finally $\mathcal{L}_{m}$ is the matter Lagrangian
of standard matter fields, $\chi_{i}$. Let us note that our action (1) when
compared with the one in Ref. Beltrán Jiménez  _et al._ (2020a) includes a
redefinition, i.e. $f(Q)\rightarrow\frac{1}{16\pi G_{N}}\left(Q+f(Q)\right)$.
They are completely equivalent but this form better fits our purpose.

We also note that in the following treatment we adopt a coordinate choice so
that the connection vanishes, following Beltrán Jiménez  _et al._ (2020a).
This choice is called coincident gauge.

In this work we will consider a flat Friedmann-Lemaître-Robertson-Walker
(FLRW) metric defined by

$ds^{2}=-dt^{2}+a(t)^{2}\delta_{i}^{j}dx^{i}dx_{j}\,,$ (2)

where $t$ is the cosmic time, $x_{i}$ are the spatial coordinates and $a(t)$
is the scale factor.

On this background the Friedman equations read Beltrán Jiménez  _et al._
(2018):

$\displaystyle H^{2}(1+2f_{Q})-\frac{1}{6}f=\frac{8\pi G_{N}}{3}\rho,$ (3)
$\displaystyle(12H^{2}f_{QQ}+f_{Q}+1)\dot{H}=-4\pi G_{N}(\rho+p)\,,$ (4)

where $H\equiv\dot{a}/a$ is the Hubble parameter and a dot stands for a
derivative with respect to cosmic time. We note that on a FLRW background
$Q=6H^{2}$ Beltrán Jiménez  _et al._ (2018, 2020b). In Eqs. (3) and (4) we
have also defined: $f_{Q}\equiv\partial f/\partial Q$,
$f_{QQ}\equiv\partial^{2}f/\partial Q^{2}$. Finally $\rho$ and $p$ are
respectively the sum of the energy density and pressure of non-relativistic
matter (baryons and cold dark matter, ”m”), photons (”$\gamma$”) and neutrinos
(both massless and massive species, ”$\nu$”). Each of them obeys the
continuity equation for perfect fluids, i.e.
$\dot{\rho}_{i}+3H(\rho_{i}+p_{i})=0$. Baryons and cold dark matter are
pressurless fluids, the energy density and the pressure for photons and
massless neutrinos are $p_{\gamma,\nu}=1/3\rho_{\gamma,\nu}$. For massive
neutrinos we assume a normal hierarchy obeying the Fermi-Dirac distribution
function.

From the set of background equations, it is clear that either the functional
form of $f(Q)$ or the Hubble expansion $H(a)$ needs to be chosen a priori. The
latter is known as the designer approach, firstly applied to $f(R)$ theory
Song  _et al._ (2007); Pogosian and Silvestri (2008a) and later to $f(Q)$
theory as well Albuquerque and Frusciante (2022). In this investigation we
want to consider a background evolution that matches the one of the standard
cosmological model, $\Lambda$CDM, in this way it is possible to solve
analytically the Friedmann equation for $f(Q)$. Then one gets Beltrán Jiménez
_et al._ (2020b):

$f(Q)=\alpha H_{0}\sqrt{Q}+6H_{0}^{2}\Omega_{\Lambda},$ (5)

where the parameter $\alpha$ is a dimensionless constant, $H_{0}\equiv H(a=1)$
is the Hubble parameter evaluated at present time ($a=1$) and
$\Omega_{\Lambda}$ is the energy density parameter of the cosmological
constant. This construction then allows to investigate deviations from
$\Lambda$CDM that appear at linear perturbation level.

A perturbed flat FLRW metric can be written in Newtonian gauge as follows

$ds^{2}=-(1+2\Psi)dt^{2}+a^{2}(1-2\Phi)\delta_{i}^{j}dx^{i}dx_{j}\,,$ (6)

where $\Phi(t,x_{i})$ and $\Psi(t,x_{i})$ are the two gravitational
potentials. In GR these two potentials are equal but for modified gravity
theories they can be different. Then the relations between the two
gravitational potentials and the gauge-invariant density contrast, $\Delta$,
are altered compared to GR. We adopt the following general way to parameterize
modified gravity in Fourier space, which consists of introducing two coupling
functions, $\mu(a,k)$ and $\gamma(a,k)$, which are defined through the
following equations Amendola  _et al._ (2008); Bean and Tangmatitham (2010);
Silvestri  _et al._ (2013); Pogosian  _et al._ (2010); Amendola  _et al._
(2020):

$\displaystyle-\frac{k^{2}}{a^{2}}\Psi$ $\displaystyle=$ $\displaystyle 4\pi
G_{N}\mu(a,k)[\rho\Delta+3(\rho+p)\sigma],$ (7)
$\displaystyle\frac{k^{2}}{a^{2}}[\phi-\gamma(a,k)\Psi]$ $\displaystyle=$
$\displaystyle 12\pi G_{N}\mu(a,k)(\rho+p)\sigma\,.$ (8)

Here we define the gauge-invariant density contrast as

$\rho\Delta\equiv\rho\delta+3\frac{aH}{k^{2}}(\rho+p)\,,$ (9)

with $\delta$ being the density contrast, $\sigma$ is the matter anisotropic
stress and $k$ is the Fourier mode. In particular $\mu$ is the effective
gravitational coupling and accounts for the modifications of gravity on the
clustering of matter. Having $\mu>1$ implies a stronger gravitational
interaction than in GR ($\mu_{GR}=1$), while a weaker gravity is obtained for
$\mu<1$. Instead $\gamma$ is known as the slip parameter but it does not have
a direct connection to observables. For this reason usually one considers
another coupling function $\Sigma(a,k)$, known as the light deflection
parameter, which describes how light travels on cosmological distances. It is
defined via the following equation Amendola  _et al._ (2008); Bean and
Tangmatitham (2010); Silvestri  _et al._ (2013); Pogosian  _et al._ (2010);
Amendola  _et al._ (2020):

$\displaystyle\frac{k^{2}}{a^{2}}(\Phi+\Psi)=-4\pi
G_{N}\Sigma(a,k)[2\rho\Delta+3(\rho+p)\sigma].$ (10)

In the limiting case of negligible matter anisotropic stress the three
coupling functions are connected:

$\Sigma=\frac{\mu}{2}(1+\gamma)\,.$ (11)

In this approach the coupling functions $\mu$, $\Sigma$ and $\gamma$ in
general can depend on both scale and time and $\gamma\neq 1$, i.e.
$\phi\neq\psi$. However, when computing analytically their forms in specific
theories of MG, by using the quasi static approximation (QSA) Boisseau  _et
al._ (2000); De Felice  _et al._ (2011), they show up to be scale independent,
e.g. Brans Dicke theory Brans and Dicke (1961), DGP Lombriser  _et al._
(2013), K-mouflage Benevento  _et al._ (2019), Galileon/Horndeski theories De
Felice  _et al._ (2011). On the contrary $f(R)$ models show an explicitly
scale dependence Pogosian and Silvestri (2008b). Additionally some of them
have $\gamma=1$, e.g. Cubic Galileon models and DGP, which seems a feature
connected with an unchanged speed of propagation of tensor modes with respect
to GR Pogosian and Silvestri (2016). Let us note that the validity of the QSA
for dark energy and modified gravity models has been extensively discussed in
Ref. Sawicki and Bellini (2015) and numerically proved to be valid for some
specific models for modes with $k\gtrsim 0.001$ h/Mpc in Refs. Pogosian and
Silvestri (2016); Frusciante and Perenon (2020).

In the case of $f(Q)$-gravity, it is possible to obtain the equations to
linear order Beltrán Jiménez  _et al._ (2020b), which assuming the QSA, leads
to $\gamma=1$, similarly to GR and to a modified Poisson equation (7) with
Beltrán Jiménez  _et al._ (2020b):

$\mu(a)=\frac{1}{1+f_{Q}}\,,$ (12)

which is scale independent. We note that due to stability requirement
$1+f_{Q}>0$. For the model we consider in this work (Eq. (5)) it reads
Frusciante (2021):

$\mu(a)=\frac{12H}{12H+\sqrt{6}\alpha H_{0}}\,.$ (13)

Such modification in the strength of the gravitational interaction has been
shown to lead to significant and measurable signatures on cosmological
observables compared to $\Lambda$CDM Frusciante (2021): when $\alpha>0$ the
gravitational interaction is weaker than in GR and this leads to a suppressed
lensing potential auto-correlation power spectrum and matter power spectrum
and to an enhanced CMB ISW-tail; the opposite holds in the case of $\alpha<0$
which corresponds to a stronger gravitational interaction. Let us note that
the model we are analysing does not have early time modifications since the
$f(Q)$ function reduces to its $\Lambda$CDM limit at early times, showing
deviations from standard model only at $z+1\lesssim 20$ Frusciante (2021).

##  III Methodology and data sets

In this work we are interested in studying the degeneracy between a modified
gravitational interaction given by the $f(Q)$ gravity and the neutrino sector.
For this purpose we adopt the Einstein-Boltzmann code MGCAMB Hojjati  _et al._
(2011); Zucca  _et al._ (2019) which implements the framework for linear
perturbations described in Section II and has an accurate modelling of massive
neutrinos Zucca  _et al._ (2019). In detail we use a modified version of it
Frusciante (2021) which includes a new patch to evolve the model in Eq. (5).
For specific details about the explicit modifications of the equations in
MGCAMB on the inclusion of massive neutrinos in the context of MG and the
$f(Q)$ model under consideration we refer the reader to Zucca  _et al._ (2019)
and Frusciante (2021) respectively. We will then use the Monte-Carlo Markov
Chain code, MGCosmoMC Hojjati  _et al._ (2011) to provide cosmological
constraints. The sampled chains from MGCosmoMC are analyzed using GetDist
Lewis (2019).

We will present in the following Sections both the phenomenology analysis and
the cosmological constraints which will be splitted into two cases:

  1. 1.

investigation of the impact of varying massive neutrinos. For this case we
will assume three neutrino mass eigenstates $m_{1}$, $m_{2}$, and $m_{3}$ with
$m_{3}\gg m_{2}>m_{1}$ (normal hierarchy) and a fixed $N_{\rm eff}=3.046$.

  2. 2.

investigation of the impact of varying the neutrino effective number, $N_{\rm
eff}$, due to neutrinos and to any additional massless particles.

Let us now present the sets of data we will use for the cosmological
constraints.

We use the Planck 2018 Aghanim  _et al._ (2020) data of CMB temperature
likelihood for large angular scales ($\ell=[2,29]$ for TT power spectrum) and
for the small angular scales a joint TT, TE and EE likelihoods
($\ell=[30,2508]$ for TT power spectrum, $\ell=[30,1996]$ for TE cross-
correlation and EE power spectra). We refer to this data set as ”Plk18”.

We also consider baryonic acoustic oscillation (BAO) data from the Sloan
Digital Sky Survey (SDSS) DR7 Main Galaxy Sample Ross  _et al._ (2015) and the
6dF Galaxy Survey Beutler  _et al._ (2011) and the joint BAO and RSD datasets
from the SDSS DR12 consensus release Alam  _et al._ (2017).

Additionally we include the Joint Light-curve Array (JLA) of Supernova Type IA
(SNIa) from the Supernova Legacy Survey (SNLS) and SDSS Betoule  _et al._
(2014). We refer to the combined analysis of Planck18, BAO, RSD and SNIa data
as ”PBRS”.

The latest data set we consider includes galaxy clustering (GC) and weak
gravitational lensing (WL) measurements from the Dark Energy Survey Year-One
(DES-1Y) data Abbott  _et al._ (2018). In details these are measurements of
the angular two-point correlation functions of GC, cosmic shear and galaxy-
galaxy lensing. We note that for the present analysis we adopt a standard DES
cutoff (by setting $\Delta\chi^{2}=5$) of the non-linear regime following Ade
_et al._ (2016b); Abbott  _et al._ (2019); Zucca  _et al._ (2019) because for
$f(Q)$-gravity a model of non-linearities is not yet available. We name the
joint analysis of DES data with the previous data sets as ”PBRSD”.

Additional details about the analysis include a flat prior on the model’s
parameter $\alpha\in[-3,3]$ and we vary the base cosmological parameters such
as the physical densities of baryons ($\Omega_{b}h^{2}\in[0.005,0.1]$, where
$h=H_{0}/100$) and cold dark matter ($\Omega_{c}h^{2}\in[0.001,0.99]$), the
primordial amplitude ($\ln(10^{10}A_{s})\in[1.61,3.91]$), the reionization
optical depth ($\tau\in[0.01,0.8]$), the angular size of the sound horizon at
recombination ($\theta_{MC}\in[0.5,10]$) and spectral index of scalar
perturbations ($n_{s}\in[0.8,1.2]$). Depending on which of the two cases we
are considering we also vary either $\Sigma m_{\nu}\in[0.01,0.9]$ or $N_{\rm
eff}\in[1.0,5.0]$. We will show the results in terms of $\sigma_{8}^{0}$
(which is a derived parameter) in place of $A_{s}$.

##  IV Phenomenology: degeneracy with the neutrino sector

In this Section we will investigate the degeneracy between the MG parameter,
$\alpha$, and, alternately, with $\Sigma m_{\nu}$ and $N_{\rm eff}$. To this
purpose we will show the effects of these parameters on some cosmological
observables such as the CMB TT, lensing and matter power spectra. We will use
the following values for the cosmological parameters:
$\Omega_{c}h^{2}=0.12011$, $\Omega_{b}h^{2}=0.02238$, $H_{0}=67.32$ km/s/Mpc,
$A_{s}=3.0448\times 10^{-9}$ and $n_{s}=0.966$. For the MG parameter we
consider two cases: $\alpha=0.5$ and $\alpha=-0.5$. Finally when studying the
impact of the total mass of neutrinos we have $\Omega_{\nu}h^{2}=0.008602$
with $\Sigma m_{\nu}=0.8$ eV and $\Omega_{\nu}h^{2}=0.000645$ with $\Sigma
m_{\nu}=0.06$ eV and we fix $N_{\rm eff}=3.046$. When investigating the
degeneracy with $N_{\rm eff}$ we select $N_{\rm eff}=3.046$ and $N_{\rm
eff}=3.8$ and massive neutrinos is set to zero. Let us stress that these
values are purely illustrative.

Figure 1: CMB TT power spectrum
($D^{TT}_{\ell}=\ell(\ell+1)C^{TT}_{\ell}/(2\pi)$, upper panel), lensing
potential auto-correlation power spectra
($D^{\phi\phi}_{\ell}=\ell(\ell+1)C^{\phi\phi}_{\ell}/(2\pi)$, middle panel)
and matter power spectra $P(k)$ (bottom panel) for $f(Q)$-gravity and
$\Lambda$CDM when two values for $\Sigma m_{\nu}$ are chosen. For each of them
we provide the percentage relative difference ($\Delta$) with respect to the
$\Lambda$CDM model with $\Sigma m_{\nu}=0.06$ eV.

  * •

The neutrino mass

In Figure 1 we show the impact on the cosmological observables of varying the
total mass of neutrinos and the modified gravity parameter. These have been
already discussed in the Introduction and Sec. II respectively. Let us now
focus on their degeneracy.

The specific modification of gravity we are considering impacts the CMB TT
power spectrum at low-$\ell$ (top panels) as an increase of the sum of
neutrino mass does. In particular we can notice that a weaker gravitational
interaction ($\alpha>0$) has the opposite effect of a larger $\Sigma m_{\nu}$:
while the former enhances the CMB ISW tail, the latter suppresses it, implying
that the values of these two parameters can be fine-tuned to cancel their
effects. Thus showing the existence of a degeneracy between them. When the
gravitational interaction is stronger ($\alpha<0$) the ISW-tail is suppressed
by modified gravity and the inclusion of massive neutrinos will further
suppress the power spectrum at low-$\ell$. This means that though the two
effects cannot cancel each other, they can be fine-tuned to obtain a given CMB
TT power spectrum at low-$\ell$. What can make a difference to break the
degeneracy in the CMB TT power spectrum between these two parameter is the
impact $\alpha$ and $\Sigma m_{\nu}$ have on the CMB TT power spectrum at
high-$\ell$. The former does not change its shape because the background is
the same of $\Lambda$CDM, while the higher is the total mass of the neutrinos
the larger is the shift of the acoustic peaks to higher multipoles. This is an
important aspect because data sets which are able to constrain the background
evolution can be used to strongly constrain $\Sigma m_{\nu}$, thus mitigating
its degeneracy with $\alpha$.

In the two middle panels of Figure 1, we show the effects of varying $\Sigma
m_{\nu}$ and $\alpha$ on the lensing potential auto-correlation power spectrum
and we also show the relative percentage difference with respect to the
$\Lambda$CDM model with $\Sigma m_{\nu}=0.06$ eV. An increase of the neutrino
mass reduces the lensing effect as it is well known Lattanzi and Gerbino
(2018); Abazajian  _et al._ (2015), similarly to the effect a positive
$\alpha$ has. On the contrary a negative value of $\alpha$ increases it. Let
us note that differently of what happens in the ISW tail of the CMB spectrum,
here the direction of the degeneracies goes in the opposite direction, i.e.
while in the ISW tail of CMB spectrum an effect due to a positive $\alpha$ can
be compensated by increasing $\Sigma m_{\nu}$, in the lensing power spectrum
they will both contribute in suppressing its power. This is another clear
pattern that can contribute in breaking the degeneracy between these two
parameters.

Finally we discuss the degeneracy that appears in the matter power spectrum
(bottom panels). Increasing the total mass of neutrinos damps the matter power
spectrum on large $k$, same effect generated by a weak gravitational
interaction. Instead a stronger gravitational interaction enhances the matter
power spectrum. We note that a compensation of the MG effect and that of an
increasing $\Sigma m_{\nu}$ might be hard to obtain due to the scale
dependence introduced by the massive neutrinos. A scale dependence that is not
present in our $f(Q)$-gravity. Let us note that while here we focus on massive
neutrinos and $\alpha$ degeneracy, other degeneracies between $\alpha$/$\Sigma
m_{\nu}$ and other cosmological parameters might exist.

In summary we have identified and discussed degeneracy effects between
$\alpha$ and $\Sigma m_{\nu}$ but we have also highlighted those patterns that
can allow us to break their degeneracies.

Figure 2: CMB TT power spectrum
($D^{TT}_{\ell}=\ell(\ell+1)C^{TT}_{\ell}/(2\pi)$, upper panel), lensing
potential auto-correlation power spectra
($D^{\phi\phi}_{\ell}=\ell(\ell+1)C^{\phi\phi}_{\ell}/(2\pi)$, middle panel)
and matter power spectra $P(k)$ (bottom panel) for $f(Q)$-gravity and
$\Lambda$CDM when two values for $N_{\rm eff}$ are chosen. For each of them we
provide the percentage relative difference ($\Delta$) with respect to the
$\Lambda$CDM model with $N_{\rm eff}=3.046$. Figure 3: Marginalised
constraints at 68% (darker) and 95% (lighter) C.L. on the model parameter
$\alpha$ and five cosmological parameters $\Sigma m_{\nu}$, $n_{s}$, $H_{0}$,
$\Omega_{m}^{0}$ and $\sigma_{8}^{0}$ obtained with the CMB data from Planck
2018 (PLK18, red), its combination with BAO, RSD and SNIa data (PBRS, green)
and with DES data (PBRSD, blue).

  * •

The neutrino effective number

In Figure 2 we show the effects of varying $N_{\rm eff}$ along with $\alpha$
on some cosmological observables.

In the top two panels we can notice that modified gravity and an increasing
$N_{\rm eff}$ have completely different effects on the CMB spectrum: while the
former affects the amplitude of the ISW-tail at low-$\ell$ due to a late time
ISW effect, the latter impacts the shape and positions of the acoustic peaks
at high-$\ell$ due to a combined effect of the early ISW effect and background
cosmology.

In the two middle panels we show the effects on the lensing power spectrum. A
higher value of $N_{\rm eff}$ suppresses the lensing spectrum and, as we
already saw before, this is a similar effect produced by a weaker
gravitational interaction, contrary to what happen when the gravitational
interaction is stronger. While here some degeneracy might exist we need to
consider that the dependence of $N_{\rm eff}$ and MG from the angular scale is
different as it is evident in the plot of the percentage relative difference.

Finally in the bottom two panels we show the matter power spectrum and the
percentage relative difference. Here we find the same situation discussed for
the case of massive neutrinos: a degeneracy exists but the scale dependence
and the damping of the acoustic scale introduced by a larger value of $N_{\rm
eff}$ make the difference.

In summary the degeneracy between MG and $N_{\rm eff}$ appears to be less
relevant than the one identified for the case of $\Sigma m_{\nu}$.

Figure 4: Top panel: marginalized cosmological parameters in the case of
$f(Q)$-gravity with varying massive neutrinos (solid lines, results in this
paper) and the case with fixed total mass of neutrinos $\Sigma m_{\nu}=0.06$
eV in Ref. Atayde and Frusciante (2021b) (dashed lines). Bottom panel:
marginalized cosmological parameters in the case of $f(Q)$-gravity with
varying $N_{\rm eff}$ (solid lines, results in this paper) and the case of
fixed $N_{\rm eff}=3.046$ in Ref. Atayde and Frusciante (2021b) (dashed
lines).

##  V Cosmological constraints

In this Section we present the cosmological constraints on the model parameter
and cosmological parameters for the cases of varying the total mass of
neutrinos and the neutrino effective number. They can be found respectively in
Tables 1 and 2 which also include the constraints for the respective cases in
the $\Lambda$CDM scenario for comparison.

  * •

The neutrino mass

We now discuss the cosmological constraints for the $f(Q)+\Sigma m_{\nu}$
case. From Fig. 3 we can notice that the stronger upper bound on $\Sigma
m_{\nu}$ is found for the PBRS data ($\Sigma m_{\nu}<0.277$ eV at 95% C.L.)
which is driven by the power in constraining background data have on
$\Omega_{m}^{0}$. Regardless of the combination of data we note that for MG
the upper bounds are weaker than the corresponding cases in $\Lambda$CDM. In
particular we note that for the full combination of data the bound degrades a
bit because DES data prefer a weaker gravitational interaction, while CMB data
for the $f(Q)$ model support a stronger gravitational interaction. Though the
value of $\alpha$ is higher for this data combination with respect to the one
for Plk18, it is not sufficient to obtain a weaker gravitational interaction
which is then supplied by a large total mass of neutrinos. Looking at the
constraints on the cosmological parameters we note they are compatible within
the error with those of the standard scenario. The constraints on the $\alpha$
parameter are in line with what is found in Ref. Atayde and Frusciante
(2021a). In detail a larger negative mean value of $\alpha$ is preferred by
the Plk18 because it allows for a better fit of the low-$\ell$ tail of TT
power spectrum. On the contrary a higher value is preferred by DES data. In
Fig. 4 top panel we compare our results with those obtained for the same model
but fixed $\Sigma m_{\nu}=0.06$ eV Atayde and Frusciante (2021a). We note that
they are compatible within the errors but the $f(Q)+\Sigma m_{\nu}$ case in
general prefers smaller mean values for $\alpha$, this is because the effects
of $\alpha$ are then compensated by a larger $\Sigma m_{\nu}$. We also notice
that the mean values of $\sigma_{8}^{0}$ and $H_{0}$ are slightly lower for
the varying massive neutrinos case.

Varying massive neutrinos
---
Model | $\alpha$ | ${\rm{ln}}(10^{10}A_{s})$ | $n_{s}$ | $H_{0}$ | $\Omega_{m}^{0}$ | $\sigma_{8}^{0}$ | $\Sigma m_{\nu}$
$\Lambda$CDM (Plk18) | - | $3.17^{+0.10}_{-0.12}$ | $0.97\pm 0.01$ | $66.6^{+2.9}_{-4.1}$ | $0.325^{+0.054}_{-0.038}$ | $0.823^{+0.064}_{-0.080}$ | $<0.622$
$\Lambda$CDM (PBRS) | - | $3.14^{+0.09}_{-0.10}$ | $0.970\pm 0.008$ | $67.97^{+0.87}_{-0.91}$ | $0.306\pm 0.011$ | $0.839\pm 0.035$ | $<0.219$
$\Lambda$CDM (PBRSD) | - | $3.13^{+0.10}_{-0.10}$ | $0.973^{+0.010}_{-0.0093}$ | $68.08^{+0.93}_{-1.0}$ | $0.304^{+0.012}_{-0.011}$ | $0.822^{+0.035}_{-0.034}$ | $<0.277$
$f(Q)$ (Plk18) | $-0.70^{+0.65}_{-0.60}$ | $3.09\pm 0.13$ | $0.970\pm 0.010$ | $66.0^{+3.4}_{-4.2}$ | $0.332^{+0.057}_{-0.044}$ | $0.806^{+0.070}_{-0.084}$ | $<0.731$
$f(Q)$ (PBRS) | $-0.66^{+0.66}_{-0.60}$ | $3.06^{+0.12}_{-0.11}$ | $0.972\pm 0.009$ | $67.92^{+0.93}_{-1.0}$ | $0.307^{+0.012}_{-0.011}$ | $0.830^{+0.036}_{-0.037}$ | $<0.277$
$f(Q)$ (PBRSD) | $-0.09^{+0.36}_{-0.36}$ | $3.12^{+0.11}_{-0.11}$ | $0.9731^{+0.0098}_{-0.0089}$ | $68.04^{+0.95}_{-1.0}$ | $0.305^{+0.012}_{-0.011}$ | $0.819^{+0.036}_{-0.037}$ | $<0.293$
Table 1: Marginalized constraints on cosmological and model parameters at 95%
C.L. for the $\Lambda$CDM and $f(Q)$ models when the total mass of the
neutrinos is considered as free parameter.

  * •

The neutrino effective number

In Fig. 5 we show the cosmological constraints for the $f(Q)+N_{\rm eff}$
case. The bounds on $N_{\rm eff}$ are always compatible within the errors with
the prediction of the standard neutrino decoupling model. We note that
regardless of the combination of data the $f(Q)+N_{\rm eff}$ case has slightly
lower mean values for $N_{\rm eff}$ compared to the ones for the
$\Lambda$CDM$+N_{\rm eff}$ case. The constraints on the other cosmological
parameters are very similar to those of the standard scenario. When varying
$N_{\rm eff}$ we notice that the mean values of $\alpha$ are higher compared
to those we obtain for $f(Q)+\Sigma m_{\nu}$ or those discussed in Ref. Atayde
and Frusciante (2021a) for fixed $N_{\rm eff}$ and $\Sigma m_{\nu}$ (see Fig.
4 bottom panel for comparison) because higher values of $\alpha$ are necessary
to compensate the lower value of $N_{\rm eff}$.

Varying the neutrino effective number
---
Model | $\alpha$ | ${\rm{ln}}(10^{10}A_{s})$ | $n_{s}$ | $H_{0}$ | $\Omega_{m}^{0}$ | $\sigma_{8}^{0}$ | $N_{eff}$
$\Lambda$CDM (Plk18) | - | $3.14^{+0.11}_{-0.11}$ | $0.971^{+0.022}_{-0.020}$ | $68.7^{+3.6}_{-3.3}$ | $0.300^{+0.024}_{-0.025}$ | $0.862^{+0.052}_{-0.048}$ | $3.07^{+0.43}_{-0.40}$
$\Lambda$CDM (PBRS) | - | $3.106^{+0.082}_{-0.084}$ | $0.966^{+0.013}_{-0.013}$ | $67.9^{+2.0}_{-2.0}$ | $0.304^{+0.011}_{-0.011}$ | $0.846^{+0.039}_{-0.037}$ | $2.98^{+0.32}_{-0.33}$
$\Lambda$CDM (PBRSD) | - | $3.072^{+0.080}_{-0.080}$ | $0.964^{+0.012}_{-0.013}$ | $67.8^{+1.9}_{-2.1}$ | $0.302^{+0.011}_{-0.011}$ | $0.828^{+0.035}_{-0.034}$ | $2.91^{+0.29}_{-0.34}$
$f(Q)$ (Plk18) | $-0.62^{+0.59}_{-0.59}$ | $3.06^{+0.13}_{-0.11}$ | $0.970^{+0.021}_{-0.019}$ | $68.8^{+3.5}_{-3.2}$ | $0.297^{+0.023}_{-0.023}$ | $0.860^{+0.047}_{-0.044}$ | $3.03^{+0.42}_{-0.39}$
$f(Q)$ (PBRS) | $-0.50^{+0.59}_{-0.54}$ | $3.03^{+0.10}_{-0.091}$ | $0.964^{+0.013}_{-0.013}$ | $67.7^{+2.0}_{-2.1}$ | $0.304^{+0.011}_{-0.011}$ | $0.841^{+0.036}_{-0.032}$ | $2.93^{+0.31}_{-0.34}$
$f(Q)$ (PBRSD) | $-0.06^{+0.36}_{-0.37}$ | $3.061^{+0.099}_{-0.099}$ | $0.964^{+0.012}_{-0.013}$ | $67.8^{+2.0}_{-2.0}$ | $0.302^{+0.011}_{-0.011}$ | $0.827^{+0.036}_{-0.034}$ | $2.91^{+0.31}_{-0.33}$
Table 2: Marginalized constraints on cosmological and model parameters at 95%
C.L. for the $\Lambda$CDM and $f(Q)$ model when $N_{\rm eff}$ is considered as
free parameter.  Figure 5: Marginalised constraints at 68% (darker) and 95%
(lighter) C.L. on the model parameter $\alpha$ and five cosmological
parameters $N_{\rm eff}$, $n_{s}$, $H_{0}$, $\Omega_{m}^{0}$ and
$\sigma_{8}^{0}$ obtained with the CMB data from Planck 2018 (PLK18, red), its
combination with BAO, RSD and SNIa data (PBRS, green) and with DES data
(PBRSD, blue).

##  VI Model selection analysis

We provide a model selection analysis based on the $\chi^{2}$-statistics and
the Deviance Information Criterion (DIC) Spiegelhalter  _et al._ (2014). While
the former allows us to measure how a model compares to actual observed data,
the latter will allow us to quantify the preference of one model over the
other and it is particulary useful when the posterior distributions of the
models have been obtained by MCMC simulation. In particular in both cases we
will compute the difference of the $f(Q)$ model with respect to $\Lambda$CDM:
$\Delta\chi^{2}=\chi^{2}_{\rm f(Q)}-\chi^{2}_{\Lambda{\rm CDM}}$ and
$\Delta\text{DIC}=\text{DIC}_{\text{f(Q)}}-\text{DIC}_{\text{$\Lambda$CDM}}$.
Before moving to the actual discussion of results, let us introduce the
definition we use to compute the DIC, which is

$\text{DIC}:=\chi_{\text{eff}}^{2}+2p_{\text{D}},$ (14)

where $\chi_{\text{eff}}^{2}$ is the value of the effective $\chi^{2}$
corresponding to the maximum likelihood and
$p_{\text{D}}=\overline{\chi}_{\text{eff}}^{2}-\chi_{\text{eff}}^{2}$, with
the bar being the average of the posterior distribution. This definition
allows us to take into account both the goodness of fit and the bayesian
complexity of the model, according to which more complex models are
disfavored. Therefore models with smaller DIC are preferred to models with
larger DIC Liddle (2009); Peirone  _et al._ (2019a, b); Frusciante  _et al._
(2020); Frusciante and Benetti (2021); Anagnostopoulos  _et al._ (2021);
Rezaei and Malekjani (2021). The values for both the $\Delta\chi^{2}$ and
$\Delta{\rm DIC}$ for the case of varying the mass of neutrinos are shown in
Table 3 and for the case of varying the effective number are Table 4.

In the case of $f(Q)+\Sigma m_{\nu}$ we find that, regardless of the
combination of data, the model can fit the data better than
$\Lambda$CDM$+\Sigma m_{\nu}$ according to the $\chi^{2}$-statistics. This is
mostly driven by a better fit to CMB data. Additionally according to the DIC
the model is also statistically preferred by data specifically for the Plk18
and PBRS data sets, while no evidence is found when GC and WL data are
included. However let us note that the $\Delta\text{DIC}$ between $f(Q)+\Sigma
m_{\nu}$ and $\Lambda$CDM$+\Sigma m_{\nu}$, when computed for the full
combination, gets better if compared to the one found between $f(Q)$ and
$\Lambda$CDM for a fixed $\Sigma m_{\nu}$, $\Delta\text{DIC}=4.8203$ Atayde
and Frusciante (2021a). In the case under consideration a large mass of the
neutrinos can suppress the matter power spectrum and compensate the
enhancement due to the stronger gravitational interaction which is preferred
by CMB data.

In the case of $f(Q)+N_{\rm eff}$, the $\chi^{2}$-statistics tells us that the
model fits better the data with respect to $\Lambda$CDM$+N_{\rm eff}$ for all
combinations. The $\Delta\text{DIC}$ supports the model over the standard
scenario when considering CMB, BAO, RSD and SNIa, however this is not the case
when we include DES data. That is because, though the value of $N_{\rm eff}$
is lowered, it cannot compensate the higher value of $\alpha$ preferred by DES
data which degrades the CMB fit.

In summary the model selection analysis we performed shows a support for the
$f(Q)$ model over $\Lambda$CDM for the CMB, BAO, RSD and SNIa data. When DES
data are included we do not find a definite conclusion when we consider
$f(Q)+\Sigma m_{\nu}$ which needs to be further explored, possibly including
non-linear scales, and a moderate evidence against $f(Q)+N_{\rm eff}$
($\Delta\text{DIC}>5$ according to Jeffrey’s scale Kass and Raftery (1995)).

Varying massive neutrinos
---
Data | $\Delta\chi^{2}$ | $\Delta DIC$
Plk18 | -2.284 | -6.547
PBRS | -2.938 | -3.987
PBRSD | -2.470 | 0.732
Table 3:  Results for the $\Delta\chi^{2}=\chi^{2}_{\rm
f(Q)}-\chi^{2}_{\Lambda{\rm CDM}}$ and $\Delta{\rm DIC}={\rm DIC}_{\rm
f(Q)}-{\rm DIC}_{\Lambda{\rm CDM}}$.  Varying $N_{\rm eff}$
---
Data | $\Delta\chi^{2}$ | $\Delta DIC$
Plk18 | -5.282 | -1.975
PBRS | -3.204 | -3.153
PBRSD | -2.468 | 4.985
Table 4:  Results for the $\Delta\chi^{2}=\chi^{2}_{\rm
f(Q)}-\chi^{2}_{\Lambda{\rm CDM}}$ and $\Delta{\rm DIC}={\rm DIC}_{\rm
f(Q)}-{\rm DIC}_{\Lambda{\rm CDM}}$.

##  VII Conclusions

We studied the degeneracy in the cosmological observables between either the
total mass of neutrinos or the effective number of neutrino species with the
modification of gravity introduced by an $f(Q)$ model which is characterized
by a background evolution identical to the $\Lambda$CDM one. In this way we
limited the degeneracy at the linear perturbation level only. We found that
degeneracy exists and that the one between MG and $\Sigma m_{\nu}$ is more
relevant than the one present with $N_{\rm eff}$. In particular we noticed
that such degeneracies can be mitigated by a combination with background data
which will strongly constrain either $N_{\rm eff}$ and $\Sigma m_{\nu}$ since
the MG effects enter only at linear perturbation level given the choice we
made for the form of the $f(Q)$ function.

We also put observational constraints on the model parameter $\alpha$ and the
other cosmological parameters, separately considering the case of varying
$\Sigma m_{\nu}$ and $N_{\rm eff}$ by adopting MCMC simulation with CMB, BAO,
RSD, SNIa, GC and WL measurements. The strongest upper bound on $\Sigma
m_{\nu}$ is obtained with the combination of CMB, BAO, RSD, SNIa with the
constraining power coming from BAO, RSD, SNIa. $N_{\rm eff}$ is always found
to be compatible with the standard value of $N_{\rm eff}=3.046$ but with lower
mean values. The best-fit values of the $f(Q)$ model gave for all combinations
of data a better effective $\chi^{2}$-statistic than the standard cosmological
scenario. When we considered the DIC statistic we noted that the combinations
of Plk18 alone and the one with BAO, RSD, SNIa were those for which we noted a
statistical support by data. When the combination includes DES data we did not
find an evidence in support or against $f(Q)$ in the case of varying $\Sigma
m_{\nu}$ but we found a moderate evidence against $f(Q)+N_{\rm eff}$. Thus the
$f(Q)+\Sigma m_{\nu}$ model can be a compelling and viable candidate for
explaining the late time acceleration.

In the future we will perform further investigations which will consider a
background for $f(Q)$ which is not exactly $\Lambda$CDM. Of interest will be
also the inclusion of non-linear scales in modeling the GC and WL so that we
can fully consider DES data.

###### Acknowledgements.

We thank Bruno J. Barros and Francesco Pace for useful comments on the
results. L.A. is supported by Fundação para a Ciência e a Tecnologia (FCT)
through the research grants UIDB/04434/2020, UIDP/04434/2020 and from the FCT
PhD fellowship grant with ref. number 2022.11152.BD . N.F. is supported by the
Italian Ministry of University and Research (MUR) through the Rita Levi
Montalcini project “Tests of gravity on cosmic scales” with reference
PGR19ILFGP. The authors also acknowledge the FCT project with ref. number
PTDC/FIS-AST/0054/2021 and the COST Action CosmoVerse, CA21136, supported by
COST (European Cooperation in Science and Technology).

## References

  * Nester and Yo (1999) J. M. Nester and H.-J. Yo, Chin. J. Phys. 37, 113 (1999), arXiv:gr-qc/9809049 .
  * Adak  _et al._ (2013) M. Adak, O. Sert, M. Kalay, and M. Sari, Int. J. Mod. Phys. A 28, 1350167 (2013), arXiv:0810.2388 [gr-qc] .
  * Adak (2018) M. Adak, Int. J. Geom. Meth. Mod. Phys. 15, 1850198 (2018), arXiv:1809.01385 [gr-qc] .
  * Beltrán Jiménez  _et al._ (2019) J. Beltrán Jiménez, L. Heisenberg, and T. S. Koivisto, Universe 5, 173 (2019), arXiv:1903.06830 [hep-th] .
  * Beltrán Jiménez  _et al._ (2018) J. Beltrán Jiménez, L. Heisenberg, and T. Koivisto, Phys. Rev. D 98, 044048 (2018), arXiv:1710.03116 [gr-qc] .
  * Harko  _et al._ (2018) T. Harko, T. S. Koivisto, F. S. N. Lobo, G. J. Olmo, and D. Rubiera-Garcia, Phys. Rev. D 98, 084043 (2018), arXiv:1806.10437 [gr-qc] .
  * Xu  _et al._ (2019) Y. Xu, G. Li, T. Harko, and S.-D. Liang, Eur. Phys. J. C 79, 708 (2019), arXiv:1908.04760 [gr-qc] .
  * Beltrán Jiménez  _et al._ (2020a) J. Beltrán Jiménez, L. Heisenberg, T. S. Koivisto, and S. Pekar, Phys. Rev. D 101, 103507 (2020a), arXiv:1906.10027 [gr-qc] .
  * Järv  _et al._ (2018) L. Järv, M. Rünkla, M. Saal, and O. Vilson, Phys. Rev. D 97, 124025 (2018), arXiv:1802.00492 [gr-qc] .
  * Rünkla and Vilson (2018) M. Rünkla and O. Vilson, Phys. Rev. D 98, 084034 (2018), arXiv:1805.12197 [gr-qc] .
  * Jiménez  _et al._ (2019) J. B. Jiménez, L. Heisenberg, and T. S. Koivisto, Universe 5, 173 (2019), arXiv:1903.06830 [hep-th] .
  * Lazkoz  _et al._ (2019) R. Lazkoz, F. S. N. Lobo, M. Ortiz-Baños, and V. Salzano, Phys. Rev. D 100, 104027 (2019), arXiv:1907.13219 [gr-qc] .
  * Ayuso  _et al._ (2021) I. Ayuso, R. Lazkoz, and V. Salzano, Phys. Rev. D 103, 063505 (2021), arXiv:2012.00046 [astro-ph.CO] .
  * Khyllep  _et al._ (2021) W. Khyllep, A. Paliathanasis, and J. Dutta, Phys. Rev. D 103, 103521 (2021), arXiv:2103.08372 [gr-qc] .
  * Frusciante (2021) N. Frusciante, Phys. Rev. D 103, 044021 (2021), arXiv:2101.09242 [astro-ph.CO] .
  * Barros  _et al._ (2020) B. J. Barros, T. Barreiro, T. Koivisto, and N. J. Nunes, Phys. Dark Univ. 30, 100616 (2020), arXiv:2004.07867 [gr-qc] .
  * Anagnostopoulos  _et al._ (2021) F. K. Anagnostopoulos, S. Basilakos, and E. N. Saridakis, (2021), arXiv:2104.15123 [gr-qc] .
  * Atayde and Frusciante (2021a) L. Atayde and N. Frusciante, Phys. Rev. D 104, 064052 (2021a), arXiv:2108.10832 [astro-ph.CO] .
  * Barreira  _et al._ (2014) A. Barreira, B. Li, C. Baugh, and S. Pascoli, Phys. Rev. D 90, 023528 (2014), arXiv:1404.1365 [astro-ph.CO] .
  * Shim  _et al._ (2014) J. Shim, J. Lee, and M. Baldi, (2014), arXiv:1404.3639 [astro-ph.CO] .
  * Baldi  _et al._ (2014) M. Baldi, F. Villaescusa-Navarro, M. Viel, E. Puchwein, V. Springel, and L. Moscardini, Mon. Not. Roy. Astron. Soc. 440, 75 (2014), arXiv:1311.2588 [astro-ph.CO] .
  * He (2013) J.-h. He, Phys. Rev. D 88, 103523 (2013), arXiv:1307.4876 [astro-ph.CO] .
  * Dossett  _et al._ (2014) J. Dossett, B. Hu, and D. Parkinson, JCAP 03, 046 (2014), arXiv:1401.3980 [astro-ph.CO] .
  * Hojjati  _et al._ (2011) A. Hojjati, L. Pogosian, and G.-B. Zhao, JCAP 08, 005 (2011), arXiv:1106.4543 [astro-ph.CO] .
  * Motohashi  _et al._ (2013) H. Motohashi, A. A. Starobinsky, and J. Yokoyama, Phys. Rev. Lett. 110, 121302 (2013), arXiv:1203.6828 [astro-ph.CO] .
  * Hu  _et al._ (2015) B. Hu, M. Raveri, A. Silvestri, and N. Frusciante, Phys. Rev. D 91, 063524 (2015), arXiv:1410.5807 [astro-ph.CO] .
  * Bellomo  _et al._ (2017) N. Bellomo, E. Bellini, B. Hu, R. Jimenez, C. Pena-Garay, and L. Verde, JCAP 02, 043 (2017), arXiv:1612.02598 [astro-ph.CO] .
  * Frusciante and Benetti (2021) N. Frusciante and M. Benetti, Phys. Rev. D 103, 104060 (2021), arXiv:2005.14705 [astro-ph.CO] .
  * Ballardini  _et al._ (2020) M. Ballardini, M. Braglia, F. Finelli, D. Paoletti, A. A. Starobinsky, and C. Umiltà, JCAP 10, 044 (2020), arXiv:2004.14349 [astro-ph.CO] .
  * Sachs and Wolfe (1967) R. K. Sachs and A. M. Wolfe, Astrophys. J. 147, 73 (1967).
  * Kofman and Starobinsky (1985) L. Kofman and A. A. Starobinsky, Sov. Astron. Lett. 11, 271 (1985).
  * Hu and White (1996) W. Hu and M. J. White, Astrophys. J. 471, 30 (1996), arXiv:astro-ph/9602019 .
  * Amendola  _et al._ (2012) L. Amendola, V. Pettorino, C. Quercellini, and A. Vollmer, Phys. Rev. D 85, 103008 (2012), arXiv:1111.1404 [astro-ph.CO] .
  * Acquaviva and Baccigalupi (2006) V. Acquaviva and C. Baccigalupi, Phys. Rev. D 74, 103510 (2006), arXiv:astro-ph/0507644 .
  * Peebles (1984) P. J. E. Peebles, Astrophys. J. 284, 439 (1984).
  * Barrow and Saich (1993) J. D. Barrow and P. Saich, MNRAS 262, 717 (1993).
  * Lattanzi and Gerbino (2018) M. Lattanzi and M. Gerbino, Front. in Phys. 5, 70 (2018), arXiv:1712.07109 [astro-ph.CO] .
  * Abazajian  _et al._ (2015) K. N. Abazajian _et al._ (Topical Conveners: K.N. Abazajian, J.E. Carlstrom, A.T. Lee), Astropart. Phys. 63, 66 (2015), arXiv:1309.5383 [astro-ph.CO] .
  * Ade  _et al._ (2016a) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016a), arXiv:1502.01589 [astro-ph.CO] .
  * Frusciante  _et al._ (2020) N. Frusciante, S. Peirone, L. Atayde, and A. De Felice, Phys. Rev. D 101, 064001 (2020), arXiv:1912.07586 [astro-ph.CO] .
  * Beltrán Jiménez  _et al._ (2020b) J. Beltrán Jiménez, L. Heisenberg, T. S. Koivisto, and S. Pekar, Phys. Rev. D 101, 103507 (2020b), arXiv:1906.10027 [gr-qc] .
  * Atayde and Frusciante (2021b) L. Atayde and N. Frusciante, Phys. Rev. D 104, 064052 (2021b), arXiv:2108.10832 [astro-ph.CO] .
  * Song  _et al._ (2007) Y.-S. Song, W. Hu, and I. Sawicki, Phys. Rev. D 75, 044004 (2007), arXiv:astro-ph/0610532 .
  * Pogosian and Silvestri (2008a) L. Pogosian and A. Silvestri, Phys. Rev. D 77, 023503 (2008a), [Erratum: Phys.Rev.D 81, 049901 (2010)], arXiv:0709.0296 [astro-ph] .
  * Albuquerque and Frusciante (2022) I. S. Albuquerque and N. Frusciante, Phys. Dark Univ. 35, 100980 (2022), arXiv:2202.04637 [astro-ph.CO] .
  * Amendola  _et al._ (2008) L. Amendola, M. Kunz, and D. Sapone, JCAP 04, 013 (2008), arXiv:0704.2421 [astro-ph] .
  * Bean and Tangmatitham (2010) R. Bean and M. Tangmatitham, Phys. Rev. D 81, 083534 (2010).
  * Silvestri  _et al._ (2013) A. Silvestri, L. Pogosian, and R. V. Buniy, Phys. Rev. D 87, 104015 (2013), arXiv:1302.1193 [astro-ph.CO] .
  * Pogosian  _et al._ (2010) L. Pogosian, A. Silvestri, K. Koyama, and G.-B. Zhao, Phys. Rev. D 81, 104023 (2010), arXiv:1002.2382 [astro-ph.CO] .
  * Amendola  _et al._ (2020) L. Amendola, D. Bettoni, A. M. Pinho, and S. Casas, Universe 6, 20 (2020), arXiv:1902.06978 [astro-ph.CO] .
  * Boisseau  _et al._ (2000) B. Boisseau, G. Esposito-Farese, D. Polarski, and A. A. Starobinsky, Phys. Rev. Lett. 85, 2236 (2000), arXiv:gr-qc/0001066 .
  * De Felice  _et al._ (2011) A. De Felice, T. Kobayashi, and S. Tsujikawa, Phys. Lett. B 706, 123 (2011), arXiv:1108.4242 [gr-qc] .
  * Brans and Dicke (1961) C. Brans and R. H. Dicke, Phys. Rev. 124, 925 (1961).
  * Lombriser  _et al._ (2013) L. Lombriser, J. Yoo, and K. Koyama, Phys. Rev. D 87, 104019 (2013), arXiv:1301.3132 [astro-ph.CO] .
  * Benevento  _et al._ (2019) G. Benevento, M. Raveri, A. Lazanu, N. Bartolo, M. Liguori, P. Brax, and P. Valageas, JCAP 05, 027 (2019), arXiv:1809.09958 [astro-ph.CO] .
  * Pogosian and Silvestri (2008b) L. Pogosian and A. Silvestri, Phys. Rev. D 77, 023503 (2008b), [Erratum: Phys.Rev.D 81, 049901 (2010)], arXiv:0709.0296 [astro-ph] .
  * Pogosian and Silvestri (2016) L. Pogosian and A. Silvestri, Phys. Rev. D 94, 104014 (2016), arXiv:1606.05339 [astro-ph.CO] .
  * Sawicki and Bellini (2015) I. Sawicki and E. Bellini, Phys. Rev. D92, 084061 (2015), arXiv:1503.06831 [astro-ph.CO] .
  * Frusciante and Perenon (2020) N. Frusciante and L. Perenon, Phys. Rept. 857, 1 (2020), arXiv:1907.03150 [astro-ph.CO] .
  * Zucca  _et al._ (2019) A. Zucca, L. Pogosian, A. Silvestri, and G.-B. Zhao, JCAP 05, 001 (2019), arXiv:1901.05956 [astro-ph.CO] .
  * Lewis (2019) A. Lewis, (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A5 (2020), arXiv:1907.12875 [astro-ph.CO] .
  * Ross  _et al._ (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Beutler  _et al._ (2011) F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. Roy. Astron. Soc. 416, 3017 (2011), arXiv:1106.3366 [astro-ph.CO] .
  * Alam  _et al._ (2017) S. Alam _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 470, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * Betoule  _et al._ (2014) M. Betoule _et al._ (SDSS), Astron. Astrophys. 568, A22 (2014), arXiv:1401.4064 [astro-ph.CO] .
  * Abbott  _et al._ (2018) T. M. C. Abbott _et al._ (DES), Phys. Rev. D 98, 043526 (2018), arXiv:1708.01530 [astro-ph.CO] .
  * Ade  _et al._ (2016b) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A14 (2016b), arXiv:1502.01590 [astro-ph.CO] .
  * Abbott  _et al._ (2019) T. M. C. Abbott _et al._ (DES), Phys. Rev. D 99, 123505 (2019), arXiv:1810.02499 [astro-ph.CO] .
  * Spiegelhalter  _et al._ (2014) D. J. Spiegelhalter, N. G. Best, B. P. Carlin, and A. van der Linde, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76, 485 (2014).
  * Liddle (2009) A. R. Liddle, Ann. Rev. Nucl. Part. Sci. 59, 95 (2009), arXiv:0903.4210 [hep-th] .
  * Peirone  _et al._ (2019a) S. Peirone, G. Benevento, N. Frusciante, and S. Tsujikawa, Phys. Rev. D 100, 063540 (2019a), arXiv:1905.05166 [astro-ph.CO] .
  * Peirone  _et al._ (2019b) S. Peirone, G. Benevento, N. Frusciante, and S. Tsujikawa, Phys. Rev. D 100, 063509 (2019b), arXiv:1905.11364 [astro-ph.CO] .
  * Rezaei and Malekjani (2021) M. Rezaei and M. Malekjani, Eur. Phys. J. Plus 136, 219 (2021), arXiv:2102.10671 [astro-ph.CO] .
  * Kass and Raftery (1995) R. E. Kass and A. E. Raftery, J. Am. Statist. Assoc. 90, 773 (1995).