# New late-time constraints on $f(R)$ gravity

Suresh Kumar<EMAIL_ADDRESS>Department of Mathematics, Indira Gandhi
University, Meerpur, Haryana 122502, India  Rafael C. Nunes
<EMAIL_ADDRESS>Instituto de Física, Universidade Federal do Rio Grande
do Sul, 91501-970 Porto Alegre RS, Brazil  Divisão de Astrofísica, Instituto
Nacional de Pesquisas Espaciais, Avenida dos Astronautas 1758, São José dos
Campos, 12227-010, SP, Brazil  Supriya Pan<EMAIL_ADDRESS>Department of Mathematics, Presidency University, 86/1 College Street, Kolkata
700073, India  Institute of Systems Science, Durban University of Technology,
PO Box 1334, Durban 4000, Republic of South Africa  Priya Yadav
<EMAIL_ADDRESS>Department of Mathematics, Indira Gandhi University,
Meerpur, Haryana 122502, India

###### Abstract

Modification of general relativity (GR) inspired by theories like $f(R)$
gravity is among the most popular ones to explain the late-time acceleration
of the Universe as an alternative to the $\Lambda$CDM model. In this work, we
use the state-of-the-art BAO+BBN data and the most recent Type Ia supernovae
(SNe Ia) sample namely PantheonPlus, including the Cepheid host distances and
covariance from SH0ES samples, to robustly constrain the $f(R)$ gravity
framework via two of the most popular $f(R)$ models in literature, namely, the
Hu-Sawicki and Starobinsky models. Additionally, we consider how the time
variation of the Newton’s gravitational constant affects the supernovae
distance modulus relation. We find a minor evidence for $f(R)$ gravity under
the Hu-Sawicki dynamics from BAO+BBN and BAO+BBN+uncalibrated supernovae joint
analysis, but the inclusion of Cepheid host distances, makes the model
compatible with GR. Further, we notice tendency of this model to relax the
$H_{0}$ tension. In general, in all the analyses carried out in this study
with the late time probes, we find both the $f(R)$ models to be consistent
with GR at 95% CL.

##  I Introduction

Astronomical data are precious for modern cosmology. From the detection of the
cosmic microwave background anisotropy to the late-time dynamics of the
Universe, we have witnessed the crucial role played by the astronomical data.
For instance, the dynamics of our Universe at its late time got abruptly
changed since 1998 from the observations of Type Ia supernovae (SNe Ia) which
first reported one of the trailblazing results in modern cosmology — the
accelerating expansion of our Universe Riess  _et al._ (1998); Perlmutter  _et
al._ (1999). This late-time accelerating expansion demands that a revision of
the standard cosmology is essential and we need to invoke some exotic type of
fluids into the gravitational equations. This can be done effectively by two
distinct ways, either one can modify the matter sector of the Universe without
touching the gravitational sector described by the Einstein’s General
Relativity (GR) which leads to various Dark Energy (DE) models Peebles and
Ratra (2003); Copeland  _et al._ (2006); Sahni and Starobinsky (2006); Bamba
_et al._ (2012); Li  _et al._ (2013); Motta  _et al._ (2021), or the
Einstein’s GR can be modified in various ways, known as modified gravity (MG)
theories Nojiri and Odintsov (2006a); Sotiriou and Faraoni (2010); De Felice
and Tsujikawa (2010); Clifton  _et al._ (2012); Capozziello and De Laurentis
(2011); Koyama (2016); Cai  _et al._ (2016); Nojiri  _et al._ (2017); Ferreira
(2019); Bahamonde  _et al._ (2021). Following both the approaches, over the
last several years, a cluster of DE and MG models have been tested with the
available astronomical data (see Refs. Copeland  _et al._ (2006); Bamba  _et
al._ (2012); Li  _et al._ (2013); Nojiri and Odintsov (2006a); Sotiriou and
Faraoni (2010); De Felice and Tsujikawa (2010); Clifton  _et al._ (2012);
Capozziello and De Laurentis (2011); Koyama (2016); Cai  _et al._ (2016);
Nojiri  _et al._ (2017); Ferreira (2019); Bahamonde  _et al._ (2021) and the
references therein).

Among the existing DE and MG models, the $\Lambda$-Cold Dark Matter
($\Lambda$CDM) cosmological model, where $\Lambda>0$ acts as a DE candidate in
the context of GR, is an excellent cosmological model that fits to a large
span of astronomical datasets. Nevertheless, $\Lambda$CDM cosmology faces many
theoretical and observational challenges. Recent observations of $\Lambda$CDM-
based Planck 2018 Aghanim  _et al._ (2020a) and the SH0ES (Supernovae and
$H_{0}$ for the Equation of State of dark energy) collaboration Riess  _et
al._ (2022a, b) suggest that the Hubble constant $H_{0}$ from Planck
$\Lambda$CDM is at more than $5\sigma$ tension with the SH0ES measurement
Riess  _et al._ (2022a, b). In addition, measurements of the parameter $S_{8}$
($=\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$; $\sigma_{8}$ is amplitude of the
matter power spectrum and $\Omega_{\rm m}$ is the matter density parameter at
present time) estimated by the Planck 2018 Aghanim  _et al._ (2020a), weak
lensing experiments Heymans  _et al._ (2021); Tröster  _et al._ (2021); Secco
_et al._ (2022); Abbott  _et al._ (2022) and Redshift-Space Distortions
measurements Kazantzidis and Perivolaropoulos (2018); Skara and
Perivolaropoulos (2020); Nunes and Vagnozzi (2021) are in tension at more than
$3\sigma$. These suggest that a revision of the $\Lambda$CDM cosmology is
needed to agree with the observational evidences. As a consequence, several
alternative proposals to the $\Lambda$CDM cosmology appeared to explain such
observational discrepancies Di Valentino  _et al._ (2016); Yang  _et al._
(2019); Vagnozzi (2020); Visinelli  _et al._ (2019); Alestas  _et al._ (2020);
Di Valentino  _et al._ (2021a); Yang  _et al._ (2021a, b); Kumar (2021);
Chudaykin  _et al._ (2022); Ballardini  _et al._ (2023) (see the recent
reviews in this direction Di Valentino  _et al._ (2021b); Perivolaropoulos and
Skara (2022); Schöneberg  _et al._ (2022a); Kamionkowski and Riess (2022)).
However, despite many new and appealing cosmological models, it has been
observed that simultaneous solution to both the tensions are quite difficult
to obtain Abdalla  _et al._ (2022). Thus, understanding the nature of the
cosmological tensions and their solutions demands further attention through
new observational probes and cosmological models.

In this article we focus on one of the viable alternatives to the $\Lambda$CDM
cosmology — the modified gravity theory, where in particular, we consider the
most natural modification to the Einstein’s GR, namely, the $f(R)$ gravity
which has been greatly investigated considering both the theoretical and
observational perspectives Bessa  _et al._ (2022); Leizerovich  _et al._
(2022); Farrugia  _et al._ (2021); Pan  _et al._ (2021); Negrelli  _et al._
(2020); Capozziello  _et al._ (2019); Akarsu  _et al._ (2018); Lazkoz  _et
al._ (2018); Nunes  _et al._ (2017, 2017), as well as to assuage the current
cosmological tensions D’Agostino and Nunes (2020); Wang (2021); Odintsov  _et
al._ (2021). However, unlike in the past works, our approach in this work
significantly differs in the treatment of its observational analysis that
enters through the physics of SNe Ia, and such a difference is caused due to
the consideration of new scalar degree(s) of freedom beyond GR which results
in a time dependent Newton’s gravitational constant $G$. Such a varying $G$
may induce a redshift ($z$-) dependent effect on the peak luminosity of SNe Ia
from the mass of the white dwarf progenitors Amendola  _et al._ (1999);
Gaztanaga  _et al._ (2002); Wright and Li (2018) and this may result in
changes in the cosmological constraints of the modified gravity models. The
revision in the evolution of intrinsic luminosity of SNe Ia due to variation
of $G$ has been considered to constrain several cosmological models
Perivolaropoulos and Skara (2021); Alestas  _et al._ (2022); Sapone  _et al._
(2021); Kumar  _et al._ (2022); Ballardini and Finelli (2022).

This means that for precise understanding of the cosmology of modified gravity
theories, the impact of modified gravity theories on the astrophysics of SNe
Ia should be considered, and through the estimation of the cosmological
parameters using the modified formalism, such impact can be decoded. Following
this, the key aim of this article is to employ the above modifications to
constrain the $f(R)$ gravity models, and study the resulting implications
mainly in light of the $H_{0}$ tension. To test this hypothesis, we use for
the first time the Pantheon+ sample to constrain the free parameters of the
$f(R)$ gravity models. In addition to these perspectives, we also consider for
the first time in this work how the state-of-the-art assumptions on BAO+BBN
joint analysis can constrain the behavior of the $f(R)$ gravity models at late
times.

The article is organized as follows. In Sec. II, we provide a brief
introduction to the cosmology of $f(R)$ gravity and introduce two well known
models that we investigate in this article. In Sec. III, we describe the
observational data-sets and our methodology to constrain the baseline of the
proposed $f(R)$ gravity models. In Sec. IV, we describe the observational
constraints on the $f(R)$ models and discuss our main results. Finally, we
describe our conclusions and perspectives in Sec. V.

##  II $f(R)$ gravity and cosmology

The gravitational action of $f(R)$ gravity, in Jordan frame, is given by

$\mathcal{S}=\frac{1}{16\pi G}\;\int d^{4}x\sqrt{-g}\,\,f(R)+\mathcal{S}_{\rm
m}+\mathcal{S}_{\rm r}\,,$ (1)

where $R$ denotes the Ricci scalar and $G$ is the Newton’s gravitational
constant. Additionally, eqn. (1) includes the actions for the matter sector
($\mathcal{S}_{\rm m}$) and the radiation sector ($\mathcal{S}_{\rm r}$). We
assume that there is no interaction at the non-gravitational level between
matter sector and the radiation sector, that means both these sectors are
independently conserved. Now varying the action (1) with respect to the metric
$g_{\mu\nu}$, we obtain the gravitational equations

$\displaystyle
FG_{\mu\nu}=-\frac{1}{2}g_{\mu\nu}\left(FR-f(R)\right)+\nabla_{\mu}\nabla_{\nu}F-g_{\mu\nu}\Box
F$ $\displaystyle+8\pi G\,\left[T^{(\rm m)}_{\mu\nu}+T^{(\rm
r)}_{\mu\nu}\right]\,,$ (2)

where $G_{\mu\nu}=R_{\mu\nu}-\left(1/2\right)g_{\mu\nu}R$ stands for the
Einstein tensor; ${\nabla}_{\mu}$ is the covariant derivative, $\Box\equiv
g^{\mu\nu}{\nabla}_{\mu}{\nabla}_{\nu}$; $F=F(R)\equiv f_{,R}=df(R)/dR$
(similarly by $f_{,RR}$ we shall mean $d^{2}f(R)/dR^{2}$); $T^{(\rm
m)}_{\mu\nu}$ and $T^{(\rm r)}_{\mu\nu}$ respectively denote the energy-
momentum tensor for the matter sector and the radiation sector. Note that for
$f(R)=R$ in eqn. (1), one recovers the Einstein-Hilbert action for General
Relativity. Now we proceed towards the cosmological evolution in the context
of $f(R)$ gravity theory. As usual, we start with the homogeneous and
isotropic background of our Universe which is well described by the Friedmann-
Lemaître-Robertson-Walker (FLRW) line element

$\displaystyle\text{d}s^{2}=-\text{d}t^{2}+a^{2}(t)\left[\frac{\text{d}r^{2}}{1-kr^{2}}+r^{2}(\text{d}\theta^{2}+\sin^{2}\theta\,\text{d}\phi^{2})\right],$
(3)

where $(t,r,\theta,\phi)$ are the co-moving coordinates; $a(t)$ describes the
expansion scale factor of the Universe and $k$ corresponds to the spatial
geometry of the Universe where $k=0$, $+1$ and $-1$, respectively denote a
spatially flat, closed and open Universe. Now for the spatially flat FLRW line
element ($k=0$), eqn. (2) leads to

$\displaystyle 3FH^{2}=8\pi
G\left(\rho_{\text{m}}+\rho_{\text{r}}\right)+\frac{1}{2}\left(FR-f(R)\right)-3H\dot{F}\,,$
(4) $\displaystyle-2F\dot{H}=8\pi
G\left(\rho_{\text{m}}+p_{\text{m}}+\rho_{\text{r}}+p_{\text{r}}\right)+\ddot{F}-H\dot{F}\,,$
(5)

where an overhead dot denotes the derivative with respect to the cosmic time
$t$; $H=\dot{a}(t)/a(t)$ is the Hubble parameter; ($\rho_{\text{m}}$,
$p_{\text{m}}$), $(\rho_{\text{r}},p_{\text{r}})$ denote the (energy density,
pressure) of the matter sector and the radiation sector respectively. Note
that in the spatially flat FLRW Universe, the Ricci scalar $R$ takes the form
$R=6(2H^{2}+\dot{H})$. It might be interesting to note that the gravitational
equations (4) and (5) can also be expressed as $3H^{2}=8\pi
G\left(\rho_{m}+\rho_{r}+\rho_{\rm eff}\right)$ and $2\dot{H}=-8\pi
G\left(\rho_{m}+p_{m}+\rho_{r}+p_{r}+\rho_{\rm eff}+p_{\rm eff}\right)$
respectively, where through $\rho_{\rm eff}$, $p_{\rm eff}$, defined as

$\displaystyle\rho_{\rm eff}\equiv\frac{1}{8\pi
G}\left[\frac{1}{2}\left(FR-f\right)-3H\dot{F}+3\left(1-F\right)H^{2}\right],$
(6) $\displaystyle p_{\rm eff}\equiv\frac{1}{8\pi
G}\Bigg{[}-\frac{1}{2}\left(FR-f\right)-\left(1-F\right)\left(2\dot{H}+3H^{2}\right)$
$\displaystyle+\ddot{F}+2H\dot{F}\Bigg{]},$ (7)

one introduces an effective dark energy scenario from the modifications of the
gravitational sector where the effective equation-of-state parameter $w_{\rm
eff}\equiv p_{\rm eff}/\rho_{\rm eff}$ reads

$\displaystyle w_{\rm
eff}=-1-\frac{H\dot{F}+2\dot{H}-2F\dot{H}-\ddot{F}}{\frac{1}{2}\left(FR-f\right)-3H\dot{F}+3\left(1-F\right)H^{2}},$
(8)

which represents a deviation from the cosmological constant $w_{\Lambda}=-1$
induced by the gravitational modifications. Thus, for any modified $f(R)$
gravity model, one can calculate the effective equation-of-state and estimate
how far the model is deviating from the cosmological constant. One can notice
that $\rho_{\rm eff}$ and $p_{\rm eff}$ defined in eqns. (6), (7) satisfy the
usual conservation equation $\dot{\rho}_{\rm eff}+3H(\rho_{\rm eff}+p_{\rm
eff})=0$. As the matter and radiation sectors enjoy independent conservation,
therefore, their conservation equations can be expressed as

$\displaystyle\dot{\rho}_{\text{m}}+3H(1+w_{\text{m}})\rho_{\text{m}}=0,$ (9)
$\displaystyle\dot{\rho}_{\text{r}}+3H(1+w_{\text{r}})\rho_{\text{r}}=0,$ (10)

where $w_{\text{m}}=p_{\text{m}}/\rho_{\text{m}}$ and
$w_{\text{r}}=p_{\text{r}}/\rho_{\text{r}}$ are respectively the equation of
state parameters of the matter sector and the radiation sector. We assume the
standard cases where $w_{\text{m}}=0$ (i.e., pressure-less matter) and
$w_{\text{r}}=1/3$. Therefore, from eqns. (9) and (10), one can derive that
$\rho_{\text{m}}\propto a^{-3}$ and $\rho_{\text{r}}\propto a^{-4}$,
respectively.

Now, for a given $f(R)$ model, using the gravitational equations (4) and (5)
together with the conservation equations for the matter and radiation sectors,
in principle, one can determine the cosmological dynamics. However, an
arbitrary $f(R)$ model may suffer from a number of cosmological problems, e.g.
the matter instability Faraoni (2006), instability at the level of
perturbations Bean  _et al._ (2007), absence of matter dominated era Amendola
_et al._ (2007a), inability to satisfy the local gravity constraints Chiba
_et al._ (2007), dark energy oscillations Nojiri and Odintsov (2006b) etc.
Thus, in order to construct viable $f(R)$ models, one needs to impose the
following conditions De Felice and Tsujikawa (2010); Amendola  _et al._
(2007b):

Figure 1: Left panel: Relative difference in the expansion rate of the
Universe, $\Delta H(z)=\left(H^{f(R)\,{\rm Gravity}}(z)/H^{\Lambda{\rm
CDM}}(z)\right)-1$, for the Hu-Sawicki $f(R)$ model considering several values
of $b$. Right panel: Same as in left panel, but for the Starobinsky $f(R)$
model.

$\displaystyle f_{,R}>0\;\text{and}\;f_{,RR}>0\;\mbox{for}\leavevmode\nobreak\
R\geq R_{0}\leavevmode\nobreak\ (>0),$ (11)

where $R_{0}$ is the present value of $R$. The condition $f_{,R}>0$ ensures
that there are no ghosts and $f_{,RR}>0$ ensures the avoidance of tachyonic
instability De Felice and Tsujikawa (2010). Moreover, from the observational
perspectives, a viable $f(R)$ model reproducing the matter dominated era,
satisfying the local gravity constraints plus to be consistent with the
equivalence principle, should behave like

$\displaystyle f(R)\rightarrow R-2\Lambda,\;\text{for}\ \ R\geq R_{0},$ (12)

where $\Lambda$ is a constant and to depict a late-time stable de Sitter
solution De Felice and Tsujikawa (2010), the $f(R)$ model also needs to
satisfy

$\displaystyle 0<\left(\frac{Rf_{,RR}}{f_{,R}}\right)_{r}<1\ \ \text{at}\ \
r=-\frac{Rf_{,R}}{f}=-2.$ (13)

Combining all these conditions altogether, the viable $f(R)$ models up to two
parameters can be recast as

$\displaystyle f(R)=R-2\Lambda y(R,b),$ (14)

where the function $y(R,b)$ gives an idea about the deviation of the
underlying $f(R)$ model from GR in which $b$ is a free parameter. In the
following, we consider two viable $f(R)$ models, namely the Hu-Sawicki $f(R)$
model Hu and Sawicki (2007) and the Starobinsky $f(R)$ model Starobinsky
(2007).

  1. 1.

The Hu-Sawicki $f(R)$ model reads as Hu and Sawicki (2007)

$\displaystyle
f(R)=R-\frac{c_{1}R_{\mathrm{HS}}\left(R/R_{\mathrm{HS}}\right)^{p}}{c_{2}\left(R/R_{\mathrm{HS}}\right)^{p}+1},$
(15)

where $c_{1}$, $c_{2}$, $R_{\mathrm{HS}}$ and $p\leavevmode\nobreak\ (>0)$ are
the free parameters of the model. One can rewrite eqn. (15) to the form of
eqn. (14) where $y(R,b)$ adopts the following expression Bamba  _et al._
(2013); Basilakos  _et al._ (2013); Nunes  _et al._ (2017)

$y(R,b)=1-\frac{1}{1+\Bigl{(}\frac{R}{\Lambda b}\Bigr{)}^{p}},$ (16)

in which $c_{1}R_{\mathrm{HS}}/c_{2}=2\Lambda$ and $2c_{2}^{1-1/p}/c_{1}=b$.
In this article we consider $p=1$ as usual, since this parameter is completely
statistically degenerate. Now, one may notice that for $b\rightarrow 0$ (i.e.,
for $c_{1}\rightarrow\infty$) and $R_{\mathrm{HS}}\rightarrow 0$ with
$c_{1}R_{\mathrm{HS}}\rightarrow 2\Lambda c_{2}$, the Hu-Sawicki $f(R)$ model
can also approximate the $\Lambda$CDM cosmology, i.e., $f(R)\rightarrow
R-2\Lambda$.

  2. 2.

The Starobinsky dark energy $f(R)$ model is given by Starobinsky (2007).

$\displaystyle f(R)=R-\lambda
R_{\mathrm{S}}\left[1-\left(1+\frac{R^{2}}{R_{\mathrm{S}}^{2}}\right)^{-n}\right],$
(17)

where $\lambda\leavevmode\nobreak\ (>0)$, $R_{\mathrm{S}}$ and
$n\leavevmode\nobreak\ (>0)$ are the free parameters of this model. In a
similar fashion, one can rewrite eqn. (17) to the form of eqn. (14) where
$y(R,b)$ takes the form Basilakos  _et al._ (2013); Nunes  _et al._ (2017)

$y(R,b)=1-\frac{1}{\Bigl{[}1+\left(\frac{R}{\Lambda\,b}\right)^{2}\Bigr{]}^{n}},$
(18)

in which $\Lambda=\lambda R_{\mathrm{S}}/2$ and $b=2/\lambda$. Throughout the
article we have considered $n=1$ as usual, since this parameter is completely
statistically degenerate. Now, one may notice that for $b\rightarrow 0$ (i.e.,
for $\lambda\rightarrow\infty$) and $R_{\mathrm{S}}\rightarrow 0$ with
$\lambda R_{\mathrm{S}}\rightarrow 2\Lambda$, the Starobinsky $f(R)$ model can
also approximate the $\Lambda$CDM cosmology, i.e., $f(R)\rightarrow
R-2\Lambda$. Note that this model does not have the same functional structure
as the famous Starobinsky inflation model Starobinsky (1980), but represents a
parametric form that allows to generate an accelerated expansion at late
times, where for $R>>R_{s}$ we have $f(R)=R-2\Lambda$, where the high
curvature value of the effective cosmological constant is $\Lambda=\lambda
R_{s}/2$.

Thus, one can see that the free parameter $b$ quantifies the deviation from GR
($b=0$). In Appendix A, we show the scalar potential for both models in the
Einstein frame. Now, in order to understand the evolution of the Universe for
the proposed $f(R)$ models, one needs to trace the expansion rate of the
Universe. We use the same methodology as in Refs. Basilakos  _et al._ (2013);
Sultana  _et al._ (2022) to derive the expansion rate of the Universe, i.e.,
the $H(z)$ function, for the proposed $f(R)$ models. Fig. 1 shows the
theoretical prediction for the expansion rate of the Universe at late times
for both the models under consideration in this work taking reasonable and
different values of $b$. We quantify the difference from the $\Lambda$CDM
model by using the fixed values of $H_{0}$ and $\Omega_{\rm m}$ to their
canonical values from CMB observations Aghanim  _et al._ (2020b), i.e.,
$H_{0}=67.4$ km/s/Mpc and $\Omega_{\rm m}=0.31$. For the Hu-Sawicki model (see
the left panel of Fig. 1), we note that the expansion of the Universe is very
sensitive to $b$, irrespective of the positive or negative values as clearly
depicted here for $b\in[-0.1,0.1]$. Specifically, for $b>0$ and $z>0.38$, the
expansion rate of the Universe within the Hu-Sawicki $f(R)$ model is greater
than the $\Lambda$CDM model while for $z<0.38$, we notice the inverse
scenario. For $b<0$, the dynamics is opposite to the previous case assuming
$b>0$. On the other hand, for the Starobinsky $f(R)$ model (see the right
panel of Fig. 1), we see that irrespective of the positive and negative values
of $b$ within $[-0.1,0.1]$, the behaviour in the expansion rate within this
$f(R)$ gravity model remains similar. In fact, the expansion rate of the
Universe within this $f(R)$ model is unresponsive to variations in the sign of
parameter $b$ with respect to the $\Lambda$CDM model. Overall, we find that
for $z>0.12$, the expansion rate of the Universe is higher than the
$\Lambda$CDM model while for $z<0.12$, the expansion rate of the Universe is
lower than the $\Lambda$CDM model. In summary, one can see that the Hu-Sawicki
and the Starobinsky $f(R)$ models are quantitatively not the same at late
times.

##  III Data and Methodology

In order to derive constraints on the model baseline, we use the following
datasets.

  * •

BAO: Baryon Acoustic Oscillation (BAO) data consist of isotropic BAO
measurements of $D_{V}(z)/r_{d}$, where $D_{V}(z)$ and $r_{d}$ stand for
spherically averaged volume distance, and sound horizon at baryon drag
respectively and anisotropic BAO measurements of $D_{M}(z)/r_{d}$ and
$D_{H}(z)/r_{d}$ (with $D_{M}(z)$ the comoving angular diameter distance and
$D_{H}(z)=c/H(z)$ the Hubble distance) from the final measurements of the SDSS
collaboration that cover eight distinct redshift intervals, acquired and
ameliorated over the past 20 years Alam  _et al._ (2021). All the above
mentioned BAO-only measurements are compiled in Table 3 of Ref. Alam  _et al._
(2021). We assume that the uncertainties are Gaussian approximations to the
likelihoods for each tracer ignoring the correlations between measurements as
suggested in the SDSS collaboration paper Alam  _et al._ (2021).

  * •

BBN: The Big Bang Nucleosynthesis (BBN) are considered with the state-of-the-
art assumptions, which consist of measurements of the primordial abundances of
helium, $Y_{P}$, from Aver  _et al._ (2015), and the deuterium measurement,
$y_{DP}=10^{5}n_{D}/n_{H}$, obtained in Cooke  _et al._ (2018). This BBN
likelihood is sensitive to the physical baryon density
$\omega_{b}\equiv\Omega_{b}h^{2}$ and the effective number of neutrino species
$N_{\rm eff}$ constraints. In the present work, we fix $N_{\rm eff}=3.046$ .

  * •

Type Ia supernovae and Cepheid: Type Ia supernovae (SNe Ia) have generally
been an important astrophysical tools in establishing the standard
cosmological model. SNe Ia distance moduli measurements constrain the
uncalibrated luminosity distance $H_{0}d_{L}(z)$, or in other words the slope
of the late-time expansion rate, which as a result constrains the matter
density parameter $\Omega_{\rm m}$. For a supernova at redshift $z$, the
theoretical apparent magnitude $m_{B}$ is given by

$\displaystyle m_{B}=5\log_{10}\left[\frac{d_{L}(z)}{1Mpc}\right]+25+M_{B},$
(19)

where $M_{B}$ is the absolute magnitude. The distance modulus reads as
$\mu(z)=m_{B}-M_{B}$. The calibrated SNe Ia absolute magnitude $M_{B}$ is in
general assumed to be truly a constant, i.e., the parameter $M_{B}$ should be
independent of the redshift. It has been argued that a possible variation of
the absolute magnitude $M_{B}$ and equivalently of the absolute luminosity as
$L\sim 10^{-2M_{B}/5}$, could be due to a variation in the value of Newton’s
gravitational constant $G$ Gaztanaga  _et al._ (2002); Wright and Li (2018).
This is due to the fact that the absolute luminosity is proportional to the
Chandrasekhar mass as $L\sim M_{\rm Chandra}$, which depends on $G$ as $L\sim
G^{-3/2}$. Therefore, any modification of gravity will generate an effective
gravitational constant in the form of $G_{\rm eff}$ that will induce a natural
correction to the distance modulus. The presence of a varying effective
gravitational constant leads to rewrite eq. (19) as

$\displaystyle\mu_{\rm th}=m_{B}-M_{B}$ $\displaystyle\;\;\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
=5\log_{10}d_{L}(z)+25+\frac{15}{4}\log_{10}\frac{G_{\rm eff}(z)}{G}.$ (20)

Taking the quasi-static approximation and the modified Poisson equation, it is
well known that in $f(R)$ gravity context, we have Tsujikawa (2007)

$\frac{G_{\rm
eff}(z)}{G}=\frac{1}{f_{R}}\left(\frac{1+4k^{2}m/a^{2}}{1+3k^{2}m/a^{2}}\right),$
(21)

where $m=f_{RR}/f_{R}$ and $G_{\rm eff}(z)$ is the effective gravitational
constant in the $f(R)$ gravity framework.

Figure 2:  The effective gravitational coupling $G_{\rm eff}(z)/G$ as a
function of redshift $z$ for different and reasonable values of $b$ for the
Hu-Sawicki parameterization within the framework adopted in this work.

Note that eq. (• ‣ III) reduces to GR when $f(R)=R-2\Lambda$, i.e., the
$\Lambda$CDM model. We follow Refs. Basilakos  _et al._ (2013); Nesseris  _et
al._ (2017) and set $k=0.1$ h/Mpc, which is necessary as now the Newton’s
gravitational constant depends on the scale $k$ as well. In Fig. 2 we show the
relative difference on the effective gravitational coupling $G_{\rm eff}(z)$
for the reasonable values of the quantifying parameter $b$ under the
perspective of the Hu-Sawicki model. We can note that the effective
gravitational constant has effects of the order of $10^{-4}$ within the
adopted fix range. On the other hand, for the Starobinsky model, we note that
$G_{\rm eff}$ has corrections of the order of $10^{-11}$, and hence this can
be considered to have negligible effects.

In de la Cruz-Dombriz  _et al._ (2008), the authors discuss the validity of
some parameterizations obtained in the quasi-static approximation. In Ref.
Orjuela-Quintana and Nesseris (2023), it is demonstrated that the quasi-static
approximation for the Hu-Sawicki model is valid for the scale 0.01 h/Mpc
$<k<0.2$ h/Mpc, which is in accordance with the values adopted in our
analyses. The validity of the quasi-static approximation using the N-body
simulation has been examined in Bose  _et al._ (2015). The quasi-static regime
for the Starobinsky model was investigated in Gannouji  _et al._ (2009).
Therefore, we can consider these two most popular scenarios in literature as
stable theories within these perspectives. Note that, in principle, the
effective gravitational constant is a function of scale as well, and this
could have a significant effect on the CMB predictions on very large scale,
where the late-time integrated Sachs Wolfe effect has strong contributions.
Then the $f(R)$ gravity could contribute with extra corrections due to that
fact besides only the change on the expansion rate of the Universe and metric
potential corrections on very large scale at CMB level. Also, the adopted
scale $k$ values coincide with the cosmic scales where the oscillation effects
are predominant on the matter power spectrum, which could, in principle, also
affect the model-dependent predictions measurements of BAO. Both points are
not yet addressed in the literature, and we hope to check it out in future
communications.

We use the SNe Ia distance moduli measurements from the Pantheon+ sample Brout
_et al._ (2022), which consists of 1701 light curves of 1550 distinct SNe Ia
ranging in the redshift interval $z\in[0.001,2.26]$, publicly available at
https://pantheonplussh0es.github.io/. We refer to this dataset as
PantheonPlus. We also consider the SH0ES Cepheid host distance anchors, which
facilitate constraints on both $M_{B}$ and $H_{0}$. When utilizing SH0ES
Cepheid host distances, the SNe Ia distance residuals are modified following
the relationship eq.(14) of Ref. Brout  _et al._ (2022). We refer to this
dataset as PantheonPlus&SH0ES.

Thus, it is possible that the modification on the distance moduli induced from
the $f(R)$ gravity framework may carry useful information about the dynamics
of these scenarios.

In our analyses, we allow the parameters $\omega_{b}$, $\omega_{\rm cdm}$
(physical cold dark matter density), $H_{0}$, $b$ and $M_{B}$ (in the analyses
with PantheonPlus data) with wide ranges of flat priors. Note that the extra
correction in eq. (• ‣ III) can be interpreted as a new and time dependent
absolute magnitude in the form of $M_{B}(z)=M_{B}+15/4\;\rm
log_{10}(G_{\rm{eff}}$ $(z)/\rm G)$, where $M_{B}$ is a constant parameter.
Several cosmological tests have recently been done to test the robustness of
the constancy of the Supernova absolute magnitude $M_{B}$ Benisty  _et al._
(2023); Sapone  _et al._ (2021); Tutusaus  _et al._ (2017); Martinelli and
Tutusaus (2019); Di Valentino  _et al._ (2020), where no clear statistical
evidence proves that $M_{B}$ can be time-dependent.

For $G_{\rm eff}/G$ = 1, we recover the default case. As default from a
statistical point of view, we treat $M_{B}$ as a free and nuisance parameter
in all analyses. We ran CLASS+MontePython code Lesgourgues (2011); Blas  _et
al._ (2011); Audren  _et al._ (2013); Brinckmann and Lesgourgues (2018) using
Metropolis-Hastings mode to derive constraints on the cosmological parameters
for the $f(R)$ gravity models defined in Sec. II using several combinations of
the datasets. All of our runs reached a Gelman-Rubin convergence criterion of
$R-1<10^{-2}$. Further, we use
MCEvidence111github.com/yabebalFantaye/MCEvidence algorithm to compute the
Bayesian evidence and perform a model comparison through the Jeffreys’ scale
Vázquez  _et al._ (2012). For model comparison, we use the log-Bayesian
evidence for each of the models relative to the standard $\Lambda$CDM model,
i.e., $\Delta\rm ln\mathcal{Z}=\rm ln\mathcal{Z}_{\Lambda CDM}-\rm
ln\mathcal{Z}_{f(R)gravity}$. For the interpretation of the results, we refer
to the revised Jeffrey’s scale and accordingly the evidence is inconclusive if
$0\leq|\Delta\rm ln\mathcal{Z}|<1$, weak if $1\leq|\Delta\rm
ln\mathcal{Z}|<2.5$, moderate if $2.5\leq|\Delta\rm ln\mathcal{Z}|<5$, strong
if $5\leq|\Delta\rm ln\mathcal{Z}|<10$, and very strong if $|\Delta\rm
ln\mathcal{Z}|\geq 10$ Kass and Raftery (1995); Trotta (2008). In what
follows, we discuss the main results of our analyses.

##  IV Main results and discussions

In Table 1, we report the summary of the statistical analyses considering the
Hu-Sawicki and Starobinsky $f(R)$ models obtained from various observational
datasets, namely, BAO+BBN, BAO+BBN+PantheonPlus, BAO+BBN+PantheonPlus&SH0ES
and PantheonPlus&SH0ES. In addition, we also show the constraints on the
$\Lambda$CDM model using the same datasets in Table 1 in order to compare the
$\Lambda$CDM results with the Hu-Sawicki and Starobinsky $f(R)$ models.
Moreover, in Fig. 3, we display the parametric space at 68% CL and 95% CL for
the Hu-Sawicki (left panel) and Starobinsky (right panel) $f(R)$ models.

Table 1: Constraints at 68% CL on some selected parameters of the Hu-Sawicki, Starobinsky and $\Lambda$CDM models obtained from BAO+BBN, BAO+BBN+PantheonPlus, BAO+BBN+PantheonPlus&SH0ES and PantheonPlus&SH0ES data. Note that here $\Delta\rm ln\mathcal{Z}$ = $\rm ln\mathcal{Z}_{\Lambda CDM}$ \- $\rm ln\mathcal{Z}_{f(R)gravity.}$ Data | BAO+BBN | BAO+BBN+PantheonPlus | BAO+BBN+PantheonPlus&SH0ES | PantheonPlus&SH0ES
---|---|---|---|---
Model | Hu-Sawicki | Hu-Sawicki | Hu-Sawicki | Hu-Sawicki
| Starobinsky | Starobinsky | Starobinsky | Starobinsky
|  $\Lambda$CDM |  $\Lambda$CDM |  $\Lambda$CDM |  $\Lambda$CDM
$b$ | $0.64^{+0.38}_{-0.29}$ | $0.46^{+0.21}_{-0.15}$ | $-0.36\pm 0.21$ | $-0.003\pm 0.029$
| $-0.02\pm 0.91$ | $-0.02\pm 0.80$ | $0.01\pm 0.31$ | $-0.037^{+0.051}_{-0.044}$
| 0 | 0 | 0 | 0
$H_{0}\,[{\rm km}/{\rm s}/{\rm Mpc}]$ | $66.7^{+1.2}_{-1.0}$ | $67.16^{+0.92}_{-1.1}$ | $70.89\pm 0.92$ | $73.76\pm 0.83$
| $65.4^{+2.4}_{-1.6}$ | $66.1\pm 1.5$ | $71.53\pm 0.77$ | $74.0^{+1.0}_{-0.85}$
| $67.5^{+1.1}_{-1.2}$ | $68.3^{+1.0}_{-0.91}$ | $71.47\pm 0.68$ | $73.74\pm 0.98$
$\Omega_{\rm m}$ | $0.273^{+0.024}_{-0.031}$ | $0.268^{+0.024}_{-0.027}$ | $0.360\pm 0.019$ | $0.334\pm 0.017$
| $0.306^{+0.019}_{-0.029}$ | $0.300\pm 0.017$ | $0.335\pm 0.013$ | $0.332\pm 0.019$
| $0.297^{+0.017}_{-0.020}$ | $0.319\pm 0.013$ | $0.336\pm 0.013$ | $0.333\pm 0.018$
$M_{B}$ | $-$ | $-19.525^{+0.054}_{-0.063}$ | $-19.280\pm 0.027$ | $-19.241\pm 0.023$
| $-$ | $-19.477\pm 0.050$ | $-19.304\pm 0.024$ | $-19.235^{+0.029}_{-0.026}$
| $-$ | $-19.409^{+0.035}_{-0.031}$ | $-19.307\pm 0.022$ | $-19.242\pm 0.028$
$\Delta\rm ln\mathcal{Z}$ | $-1.62$ | $-1.84$ | $-0.84$ | $2.33$
| $-1.99$ | $-2.28$ | $0.13$ | $2.04$
| $0$ | $0$ | $0$ | $0$
Table 2: Constraints at 68% CL on some selected parameters of the Hu-Sawicki and Starobinsky models obtained from BAO+BBN+PantheonPlus&SH0ES and PantheonPlus&SH0ES data without the extra corrections on the distance moduli in eq. (• ‣ III). Data | BAO+BBN+PantheonPlus&SH0ES | PantheonPlus&SH0ES
---|---|---
Model | Hu-Sawicki | Hu-Sawicki
| Starobinsky | Starobinsky
$b$ | $-0.36^{+0.39}_{-0.42}$ | $-0.09^{+0.19}_{-0.12}$
| $0.00\pm 0.31$ | $0.02\pm 0.77$
$H_{0}\,[{\rm km}/{\rm s}/{\rm Mpc}]$ | $70.9^{+1.8}_{-1.8}$ | $73.5^{+2.1}_{-2.1}$
| $71.50\pm 0.78$ | $74.10^{+0.97}_{-1.3}$
$\Omega_{\rm m}$ | $0.361^{+0.038}_{-0.037}$ | $0.347^{+0.042}_{-0.039}$
| $0.335\pm 0.013$ | $0.300^{+0.047}_{-0.022}$
$M_{B}$ | $-19.280^{+0.053}_{-0.052}$ | $-19.236^{+0.060}_{-0.059}$
| $-19.304\pm 0.024$ | $-19.240\pm 0.028$

Figure 3: One-dimensional posterior distributions and two-dimensional
marginalized confidence regions (68% CL and 95% CL) for $b$, $\Omega_{\rm m}$
and $H_{0}$ obtained from the BAO+BBN, BAO+BBN+PantheonPlus and
BAO+BBN+PantheonPlus&SH0ES for the Hu-Sawicki model (left panel) and
Starobinsky model (right panel). The parameter $H_{0}$ is in units of
km/s/Mpc.

Figure 4: This figure shows the magnitude-redshift relation of the
PantheonPlus sample in the range $0<z<2.3$ for the best fit values from
BAO+BBN+PantheonPlus and BAO+BBN+PantheonPlus&SH0ES analysis summarized in
Table 1 for Hu-Sawicki model (left panel) and Starobinsky model (right panel).
The $\Lambda$CDM best-fit prediction is also shown in both panels.

The combined dataset BAO+BBN probes the background history of the model
independently of both CMB and supernovae data. As the $H_{0}$ tension directly
invites a straight conflict between the CMB and the local distance ladder
measurements, thus, it will be interesting to find new routes to estimate the
Hubble constant. The joint analysis BAO+BBN has been proved to be a
competitive cosmological test Schöneberg  _et al._ (2019); Cuceu  _et al._
(2019); Schöneberg  _et al._ (2022b) which can provide accurate confidence
limits on the baseline parameters of the models. Thus, we choose BAO+BBN to be
our minimum data set. As well known, the constraints on $H_{0}$ from BAO+BBN
data in the $\Lambda$CDM context fully agree with the CMB data. When applying
BAO+BBN in the context of $f(R)$ gravity, we notice this behavior, being
$H_{0}$ compatible with low values obtained in the CMB measurements, and at
1.6% and 3% accuracy from the Hu-Sawicki and Starobinsky models, respectively.
For the parameter that quantifies the deviation from GR, i.e., the parameter
$b$, we find different results in the two different $f(R)$ gravity models. For
the Hu-Sawicki model, we obtain $b>0$ at more than 68% CL for BAO+BBN
($b=0.64^{+0.38}_{-0.29}$ at 68% CL). On the other hand, for the Starobinsky
$f(R)$ model, we find that $b$ is compatible to zero within 68% CL for
BAO+BBN, that means, no deviation from GR is suggested within this $f(R)$
model for BAO+BBN. This is not unexpected because the models have different
dynamical behaviors at late times, as previously discussed (see Fig. 1). We
further noticed that the free parameter $b$ in the Starobinsky model is not
much sensitive to its sign change with results being largely symmetric under
the change of the sign on the parameter $b$. Thus, the posterior tends to be
bimodal based on the prior adopted in our analysis. The addition of SH0ES
Cepheid host distances tends only to smooth out the bimodal effect. Because of
this bimodality, the estimation of parameter $b$ turns out to be compatible
with the null hypothesis, i.e, $b=0$. Since there are no reasons to impose
some hard prior like $b>0$ and/or $b<0$, we choose to leave this parameter
free within a large prior range.

Now, we move on considering the addition of SNe Ia and Cepheid host distance
measurements from the SH0ES team, while considering the additional corrections
on the distance moduli, i.e., eq. (• ‣ III) due to the $f(R)$ gravity model
being one of the main motivations of this work. As argued in Refs. Efstathiou
(2021); Camarena and Marra (2021); Nunes and Di Valentino (2021), the tension
on $H_{0}$ should be replaced as the tension on the supernova absolute
magnitude $M_{B}$, as the estimate of $H_{0}$ from SH0ES collaboration comes
directly from the estimate of $M_{B}$. So, in our analysis we first consider
the uncalibrated supernovae sample, i.e., the PantheonPlus. When analyzing
with BAO+BBN+PantheonPlus, we find that $b>0$ at 95% CL
($b=0.46^{+0.35}_{-0.38}$) for the Hu-Sawicki $f(R)$ model, while for the
Starobinksy $f(R)$ model, the null hypothesis $b=0$ is fully compatible within
68% CL. The constraints on $H_{0}$ from BAO+BBN+PantheonPlus for both the Hu-
Sawicki and Starobinksy $f(R)$ models are fully compatible with the estimates
of $H_{0}$ obtained from BAO+BBN.

Now, following Ref. Brout  _et al._ (2022) we consider the inclusion of the
Cepheid-host distances measurements in direct combination with the SNe Ia
sample, i.e., the full dataset PantheonPlus&SH0ES. We notice that the
parameter $b$ in both cases becomes fully compatible with GR, i.e., $b=0$.
Thus, the inclusion of the Cepheid host distances and the full covariance
matrix from SH0ES samples, makes the dynamics of the $f(R)$ models similar to
$\Lambda$CDM. On the other hand, $H_{0}$ and $\Omega_{\rm m}$ get larger
values compared to the previous analyses without the inclusion of SH0ES
measurement. Moreover, it is possible to quantify the level of tension between
two estimates $H_{0,i}$ and $H_{0,j}$ of $H_{0}$ by means of the simple
1-dimensional tension metric, which can be constructed as

$T_{H_{0}}\equiv\frac{|H_{0,i}-H_{0,j}|}{\sqrt{\sigma^{2}_{H_{0,i}}+\sigma^{2}_{H_{0,j}}}}\,,$
(22)

measured in equivalent Gaussian standard deviations. In particular, we find
that for the Hu-Sawicki $f(R)$ model (Starobinsky $f(R)$ model), the $H_{0}$
value from BAO+BBN+PantheonPlus&SH0ES is at 2.8$\sigma$ (3.2$\sigma$) and
2.9$\sigma$ (2.9$\sigma$) tensions with the $H_{0}$ values from the
BAO+BBN+PantheonPlus and BAO+BBN analyses, respectively.

We also consider the PantheonPlus&SH0ES data without external probes.
Considering the fact that this sample is at more than 2$\sigma$ tension with
BAO+BBN, we analyze their effects separately. As also shown in Ref. Brout  _et
al._ (2022), for the flat $\Lambda$CDM, the joint analysis from
PantheonPlus&SH0ES tends to generate high values of $H_{0}$ (see our results
in Table 1). When analyzing PantheonPlus&SH0ES for the $f(R)$ gravity model,
we noticed the same behavior, i.e., $H_{0}$ gets high values compatible with
local measurements. That is, without external probes, $H_{0}$ inferred from
PantheonPlus&SH0ES sample for non-standard models of type $f(R)$ gravity,
$H_{0}$ is constrained to high values. We note that Ref. Dhawan  _et al._
(2020) pointed out a possible insensitivity of the local $H_{0}$ constraint
from the Cepheid distance ladder in some model beyond the $\Lambda$CDM
cosmology. In this sense, and based on the present results, we can conclude
that the same is valid for scenarios like $f(R)$ gravity.

It is well known that BAO+SNe Ia joint analysis prefers low values of $H_{0}$
compatible with CMB observations, and on the other hand, this joint analysis
is in tension with SNe Ia+Cepheid sample analysis (see general discussions
introduced in Efstathiou (2021)). Thus, from the analyses presented here for
BAO+BBN+PantheonPlus and PantheonPlus&SH0ES, we find that the $f(R)$ gravity
does not significantly change the local distance ladder value of $H_{0}$, and
therefore these models are not able to solve the $H_{0}$ tension in the light
of the late time probes. In the analyses of the Hu-Sawicki model with the
BAO+BBN and BAO+BBN+PantheonPlus data, we notice negative correlation between
$H_{0}$ and $b$ (see left panel of Fig. 3). So smaller values of $b$
correspond to the larger values of $H_{0}$, while we notice opposite scenario
in the analysis with BAO+BBN+PantheonPlus&SH0ES data. On the other hand, in
all the analyses of the Starobinsky $f(R)$ model, the $b$ parameter is
insensitive to $H_{0}$. Similar and equivalent conclusions can be drawn from
the point of view of the $M_{B}$ estimation.

It has been shown in the literature that a phantom dark energy model or
scenarios that provide an effective phantom behavior at late times, may be the
simplest solutions to the Hubble tension problem Di Valentino  _et al._ (2016,
2021a); Vagnozzi (2020); Yang  _et al._ (2021b, a); Visinelli  _et al._
(2019); Yang  _et al._ (2019); Chudaykin  _et al._ (2022); Ballardini  _et
al._ (2023); Kumar (2021). Also it can be seen in Alestas  _et al._ (2022,
2021, 2020) that the $H_{0}$ tension is alleviated due to the impact of an
effective gravitational coupling evolution. As investigated in Arjona  _et
al._ (2019), the Hu-Sawicki model studied here has an equation of state which
crosses at late times from a phantom to quintessence dynamics for $b>0$, and
on the contrary for $b<0$. Thus, we can note that the joint analysis with BAO
+ BBN has a best-fit with a tendency towards the quintessential behavior, i.e,
$b>0$, while for the BAO + BBN + PantheonPlus&SH0ES combination, the best-fit
indicates a phantom behavior. It is important to emphasize that there is no
clear evidence for $f(R)$ gravity dynamics at more than 1$\sigma$. Thus, the
Hubble tension is alleviated because of the effective phantom dark energy at
late times induced by the $f(R)$ gravity dynamics.

From Table 1, we notice $|\Delta\rm ln\mathcal{Z}|<2.5$ for all the analyses.
Therefore, the Bayesian evidence of $f(R)$ models compared to $\Lambda$CDM is
either weak or inconclusive. So the $f(R)$ models cannot be discriminated from
$\Lambda$CDM statistically in all analysis carried out here. Finally, in Fig.
4 we show the magnitude-redshift relation of the PantheonPlus sample for the
best fit values from BAO+BBN+PantheonPlus and BAO+BBN+PantheonPlus&SH0ES for
both the $f(R)$ gravity models under consideration in the work plus the
reference $\Lambda$CDM model, where we note that the model’s predictions for
low-$z$ are almost indistinguishable from each other, but it may slightly
differ at high-$z$.

Finally, in Table 2 we show the summary of the statistical analyses
considering both the $f(R)$ gravity models investigated in this work, but
without taking into account the extra correction on the distance modulus in
eq. (• ‣ III), i.e, the distance modulus is now only affected due to the
expansion rate of the Universe, $H(z)$, predicted by each $f(R)$ model. Thus,
we can quantify statistically how much the effective time variation of the
Newton’s gravitational constant affects the observational constrains on the
models. From our analyses we find that the parameter $b$ gets affected. For
the Hu-Sawicki model, we noticed an 86% and 18% improvement in the constraint
on $b$ when considering the effective gravitational coupling on eq. (• ‣ III)
from the BAO+BBN+PantheonPlus&SH0ES and PantheonPlus&SH0ES joint analysis,
respectively. Clearly a significant improvement is observed in the parameter
$b$ within the Hu-Sawicki framework. Also, clearly the other baseline
parameters are also affected by these considerations. While comparing the
constraints between the parameters displayed in Tables 1 and 2, a significant
and noticeable improvement on the error bars are clear when the presence of a
varying effective gravitational constant is considered on the distance
modulus. Therefore, one can conclude that the presence of the extra correction
on the distance modulus certainly improves the observational constraints on
the cosmological parameters. As described above, changes on the parameter $b$
in the Starobinsky framework has low impact on the observables. Thus, for this
specific case, within the formalism used in this work, this model does not
show any improvements in this regard. In all cases, there was no change in
conclusions referring to the Bayesian evidence calculation.

##  V Final Remarks

In this work, we have considered new extra degree of freedom of the
gravitational origin by modifying the gravity sector that possesses GR as a
particular limit. The new scalar degree(s) of freedom are proposals of intense
investigation as alternatives to $\Lambda$CDM frameworks in the last two
decades. Certainly, one of the most popular theories to explain the late-time
acceleration in this sense is the $f(R)$ gravity theory. In this work, we have
presented an update of observational constraints with new perspectives on two
well known and widely used $f(R)$ gravity models, viz., Hu-Sawicki and
Starobinsky models. The robustness of the state-of-the-art assumptions on
BAO+BBN data is used for the first time to constrain the dynamics of these
models. Then the most recent SNe Ia data, by taking the time variation of the
Newton’s gravitational constant over cosmic time to correct the supernovae
distance modulus relation predictions, are used in the joint analysis with
BAO+BBN. Finally, the inclusion of the very low-z Cepheid host distances,
including the full covariance of the SH0ES sample, is considered to
investigate the $f(R)$ models under consideration in this work. We have found
a minor evidence for $f(R)$ gravity under the Hu-Sawicki dynamics from BAO+BBN
and BAO+BBN+uncalibrated supernovae joint analysis, but the inclusion of
Cepheid host distances, makes the model compatible with GR. In general, in all
the analyses, we find that $b$ is consistent with 0 at 95% CL for both of the
$f(R)$ models. So we have not found any significant deviation from GR, i.e.,
$b=0$, after the application of late-time data sets. For the Hu-Sawicki model,
we have noticed correlation between $b$ and $H_{0}$ from different
observational data sets, this shows the tendency of the model to relax the
$H_{0}$ tension. Furthermore, the free parameter $b$ of the theories still is
weakly constrained which is clearly observed from its large error bars. The
generalization of the perspectives considered here can be carried out with CMB
data from Planck, the Atacama Cosmology Telescope and with full-shape galaxy
power spectrum sample in the light of the $H_{0}$ tension. The full-shape
galaxy power spectrum has recently been well developed for non-standard models
Philcox  _et al._ (2020); Chudaykin  _et al._ (2021); Nunes  _et al._ (2022);
Simon  _et al._ (2023); Carrilho  _et al._ (2023); Reeves  _et al._ (2023);
Simon (2022), including modified gravity scenarios Piga  _et al._ (2023). We
hope to report the results in this direction in future communications.

##  Appendix A $f(R)$ gravity in the Einstein frame

In this appendix, we write the models under consideration in this work in the
Einstein frame. In what follows, we follow De Felice and Tsujikawa (2010);
Sotiriou and Faraoni (2010). Any $f(R)$ gravity metric can be written under
the action

$\mathcal{S}=\frac{1}{2\kappa}\;\int d^{4}x\sqrt{-g}\,\,f(R)+\mathcal{S}_{\rm
m}(g_{\mu\nu},\psi),$ (23)

where $\kappa=8\pi G$ and $\psi$ represents the matter fields. One can
introduce a new field $\chi$ and write the dynamically equivalent action

$\mathcal{S}=\frac{1}{2\kappa}\;\int
d^{4}x\sqrt{-g}[f(\chi)+f^{\prime}(\chi)(R-\chi)]+\mathcal{S}_{\rm
m}(g_{\mu\nu},\psi).$ (24)

Variation with respect to $\chi$ leads to the equation
$f^{\prime\prime}(\chi)(R-\chi)=0$. Thus, $\chi=R$ if
$f^{\prime\prime}(\chi)=0$, which reproduces the action 23. Re-parameterizing
the field $\chi$ by $\phi=f^{\prime}(\chi)$ and setting
$V(\phi)=\chi(\phi)\phi-f(\chi(\phi))$, the action takes the form

$\mathcal{S}=\frac{1}{2\kappa}\;\int d^{4}x\sqrt{-g}[\phi
R-V(\phi)]+\mathcal{S}_{\rm m}(g_{\mu\nu},\psi).$ (25)

One can perform a conformal transformation and rewrite the action above in the
Einstein frame. Specifically, by performing the conformal transformation
$\tilde{g}_{\mu\nu}=f(R)g_{\mu\nu}=\phi g_{\mu\nu}$, and the scalar field
redefinition $\phi=f^{\prime}(R)$ to $\tilde{\phi}$ with

$d\tilde{\phi}=\sqrt{\frac{3}{2\kappa}}\frac{d\phi}{\phi},$ (26)

a scalar-tensor theory is mapped into the Einstein frame, in which the new
scalar field $\tilde{\phi}$ couples minimally to the Ricci curvature and has
canonical kinetic energy, as described by the following action

$\mathcal{S}=\int
d^{4}x\sqrt{-\tilde{g}}\Big{[}\frac{\tilde{R}}{2\kappa}-\frac{1}{2}\partial^{\alpha}\tilde{\phi}\partial_{\alpha}\tilde{\phi}-U(\tilde{\phi}).\Big{]}$
(27)

For the equivalent of $f(R)$ gravity metric, we have

$\phi=f^{\prime}(R)=e^{\sqrt{2\kappa/3}\tilde{\phi}}$ (28)

and

$U(\tilde{\phi})=\frac{Rf^{\prime}(R)-f(R)}{2\kappa(f^{\prime}(R))^{2}},$ (29)

where $R=R(\tilde{\phi})$. The complete and final action takes the form

$\displaystyle\mathcal{S}$ $\displaystyle=\int
d^{4}x\sqrt{-\tilde{g}}\Big{[}\frac{\tilde{R}}{2\kappa}-\frac{1}{2}\partial^{\alpha}\tilde{\phi}\partial_{\alpha}\tilde{\phi}-U(\tilde{\phi})$
$\displaystyle+S_{M}(e^{-\sqrt{2\kappa/3}\tilde{\phi}}\tilde{g}_{\mu\nu},\psi)\Big{]}$

This represents a direct transformation from the Jordan frame to the Einstein
frame. Now, from eq. (29), we can see the type of potential for the scalar
degree of freedom. For the Hu-Sawicki model, the scalar potential is given by

$U=\dfrac{\Lambda
R^{2}}{\kappa(R+b\Lambda)^{2}}.\dfrac{1}{(f^{\prime}(R))^{2}},$ (30)

where $f^{\prime}(R)=1-\dfrac{2\Lambda^{2}b}{(R+b\Lambda)^{2}}$.
For the Starobinsky model, the scalar potential is given by

$U=-\Lambda R^{2}(2b^{2}-R^{2}-b^{2}\Lambda^{2}).\dfrac{1}{\kappa
f^{\prime}(R)},$ (31)

where
$f^{\prime}(R)=1-\dfrac{4\Lambda^{3}b^{2}R}{(R^{2}+b^{2}\Lambda^{2})^{2}}$.

Note that for $b=0$, we get $U=$ constant and the kinetic term in the action
vanishes. Then the action reduces to $\Lambda$CDM model.

###### Acknowledgements.

S.K. gratefully acknowledges support from the Science and Engineering Research
Board (SERB), Govt. of India (File No. CRG/2021/004658). R.C.N thanks the CNPq
for partial financial support under the project No. 304306/2022-3. S.P.
acknowledges the financial support from the Department of Science and
Technology (DST), Govt. of India under the Scheme “Fund for Improvement of S&T
Infrastructure (FIST)” (File No. SR/FST/MS-I/2019/41). P.Y. is supported by a
Junior Research Fellowship (CSIR/UGC Ref. No. 191620128350) from the
University Grants Commission (UGC), Govt. of India.

## References

  * Riess  _et al._ (1998) A. G. Riess _et al._ (Supernova Search Team), Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ (1999) S. Perlmutter _et al._ (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Peebles and Ratra (2003) P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559 (2003), arXiv:astro-ph/0207347 .
  * Copeland  _et al._ (2006) E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006), arXiv:hep-th/0603057 .
  * Sahni and Starobinsky (2006) V. Sahni and A. Starobinsky, Int. J. Mod. Phys. D 15, 2105 (2006), arXiv:astro-ph/0610026 .
  * Bamba  _et al._ (2012) K. Bamba, S. Capozziello, S. Nojiri, and S. D. Odintsov, Astrophys. Space Sci. 342, 155 (2012), arXiv:1205.3421 [gr-qc] .
  * Li  _et al._ (2013) M. Li, X.-D. Li, S. Wang, and Y. Wang, Front. Phys. (Beijing) 8, 828 (2013), arXiv:1209.0922 [astro-ph.CO] .
  * Motta  _et al._ (2021) V. Motta, M. A. García-Aspeitia, A. Hernández-Almada, J. Magaña, and T. Verdugo, Universe 7, 163 (2021), arXiv:2104.04642 [astro-ph.CO] .
  * Nojiri and Odintsov (2006a) S. Nojiri and S. D. Odintsov, eConf C0602061, 06 (2006a), arXiv:hep-th/0601213 .
  * Sotiriou and Faraoni (2010) T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451 (2010), arXiv:0805.1726 [gr-qc] .
  * De Felice and Tsujikawa (2010) A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Clifton  _et al._ (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Capozziello and De Laurentis (2011) S. Capozziello and M. De Laurentis, Phys. Rept. 509, 167 (2011), arXiv:1108.6266 [gr-qc] .
  * Koyama (2016) K. Koyama, Rept. Prog. Phys. 79, 046902 (2016), arXiv:1504.04623 [astro-ph.CO] .
  * Cai  _et al._ (2016) Y.-F. Cai, S. Capozziello, M. De Laurentis, and E. N. Saridakis, Rept. Prog. Phys. 79, 106901 (2016), arXiv:1511.07586 [gr-qc] .
  * Nojiri  _et al._ (2017) S. Nojiri, S. Odintsov, and V. Oikonomou, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * Ferreira (2019) P. G. Ferreira, Ann. Rev. Astron. Astrophys. 57, 335 (2019), arXiv:1902.10503 [astro-ph.CO] .
  * Bahamonde  _et al._ (2021) S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. L. Said, J. Mifsud, and E. Di Valentino, (2021), arXiv:2106.13793 [gr-qc] .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020a), arXiv:1807.06209 [astro-ph.CO] .
  * Riess  _et al._ (2022a) A. G. Riess _et al._ , Astrophys. J. Lett. 934, L7 (2022a), arXiv:2112.04510 [astro-ph.CO] .
  * Riess  _et al._ (2022b) A. G. Riess, L. Breuval, W. Yuan, S. Casertano, L. M. ~Macri, D. Scolnic, T. Cantat-Gaudin, R. I. Anderson, and M. C. Reyes, (2022b), arXiv:2208.01045 [astro-ph.CO] .
  * Heymans  _et al._ (2021) C. Heymans _et al._ , Astron. Astrophys. 646, A140 (2021), arXiv:2007.15632 [astro-ph.CO] .
  * Tröster  _et al._ (2021) T. Tröster _et al._ (KiDS), Astron. Astrophys. 649, A88 (2021), arXiv:2010.16416 [astro-ph.CO] .
  * Secco  _et al._ (2022) L. F. Secco _et al._ (DES), Phys. Rev. D 105, 023515 (2022), arXiv:2105.13544 [astro-ph.CO] .
  * Abbott  _et al._ (2022) T. M. C. Abbott _et al._ (DES), (2022), arXiv:2207.05766 [astro-ph.CO] .
  * Kazantzidis and Perivolaropoulos (2018) L. Kazantzidis and L. Perivolaropoulos, Phys. Rev. D 97, 103503 (2018), arXiv:1803.01337 [astro-ph.CO] .
  * Skara and Perivolaropoulos (2020) F. Skara and L. Perivolaropoulos, Phys. Rev. D 101, 063521 (2020), arXiv:1911.10609 [astro-ph.CO] .
  * Nunes and Vagnozzi (2021) R. C. Nunes and S. Vagnozzi, Mon. Not. Roy. Astron. Soc. 505, 5427 (2021), arXiv:2106.01208 [astro-ph.CO] .
  * Di Valentino  _et al._ (2016) E. Di Valentino, A. Melchiorri, and J. Silk, Phys. Lett. B 761, 242 (2016), arXiv:1606.00634 [astro-ph.CO] .
  * Yang  _et al._ (2019) W. Yang, S. Pan, E. Di Valentino, E. N. Saridakis, and S. Chakraborty, Phys. Rev. D99, 043543 (2019), arXiv:1810.05141 [astro-ph.CO] .
  * Vagnozzi (2020) S. Vagnozzi, Phys. Rev. D 102, 023518 (2020), arXiv:1907.07569 [astro-ph.CO] .
  * Visinelli  _et al._ (2019) L. Visinelli, S. Vagnozzi, and U. Danielsson, Symmetry 11, 1035 (2019), arXiv:1907.07953 [astro-ph.CO] .
  * Alestas  _et al._ (2020) G. Alestas, L. Kazantzidis, and L. Perivolaropoulos, Phys. Rev. D 101, 123516 (2020), arXiv:2004.08363 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021a) E. Di Valentino, A. Mukherjee, and A. A. Sen, Entropy 23, 404 (2021a), arXiv:2005.12587 [astro-ph.CO] .
  * Yang  _et al._ (2021a) W. Yang, E. Di Valentino, S. Pan, Y. Wu, and J. Lu, Mon. Not. Roy. Astron. Soc. 501, 5845 (2021a), arXiv:2101.02168 [astro-ph.CO] .
  * Yang  _et al._ (2021b) W. Yang, S. Pan, E. Di Valentino, O. Mena, and A. Melchiorri, JCAP 10, 008 (2021b), arXiv:2101.03129 [astro-ph.CO] .
  * Kumar (2021) S. Kumar, Phys. Dark Univ. 33, 100862 (2021), arXiv:2102.12902 [astro-ph.CO] .
  * Chudaykin  _et al._ (2022) A. Chudaykin, D. Gorbunov, and N. Nedelko, (2022), arXiv:2203.03666 [astro-ph.CO] .
  * Ballardini  _et al._ (2023) M. Ballardini, A. G. Ferrari, and F. Finelli, JCAP 04, 029 (2023), arXiv:2302.05291 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021b) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021b), arXiv:2103.01183 [astro-ph.CO] .
  * Perivolaropoulos and Skara (2022) L. Perivolaropoulos and F. Skara, New Astron. Rev. 95, 101659 (2022), arXiv:2105.05208 [astro-ph.CO] .
  * Schöneberg  _et al._ (2022a) N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, Phys. Rept. 984, 1 (2022a), arXiv:2107.10291 [astro-ph.CO] .
  * Kamionkowski and Riess (2022) M. Kamionkowski and A. G. Riess, (2022), arXiv:2211.04492 [astro-ph.CO] .
  * Abdalla  _et al._ (2022) E. Abdalla _et al._ , JHEAp 34, 49 (2022), arXiv:2203.06142 [astro-ph.CO] .
  * Bessa  _et al._ (2022) P. Bessa, M. Campista, and A. Bernui, Eur. Phys. J. C 82, 506 (2022), arXiv:2112.00822 [astro-ph.CO] .
  * Leizerovich  _et al._ (2022) M. Leizerovich, L. Kraiselburd, S. J. Landau, and C. G. Scóccola, Phys. Rev. D 105, 103526 (2022), arXiv:2112.01492 [astro-ph.CO] .
  * Farrugia  _et al._ (2021) C. R. Farrugia, J. Sultana, and J. Mifsud, Phys. Rev. D 104, 123503 (2021), arXiv:2106.04657 [astro-ph.CO] .
  * Pan  _et al._ (2021) Y. Pan, Y. He, J. Qi, J. Li, S. Cao, T. Liu, and J. Wang, Astrophys. J. 911, 135 (2021), arXiv:2103.05212 [astro-ph.CO] .
  * Negrelli  _et al._ (2020) C. Negrelli, L. Kraiselburd, S. J. Landau, and M. Salgado, Phys. Rev. D 101, 064005 (2020), arXiv:2002.12073 [gr-qc] .
  * Capozziello  _et al._ (2019) S. Capozziello, R. D’Agostino, and O. Luongo, Gen. Rel. Grav. 51, 2 (2019), arXiv:1806.06385 [gr-qc] .
  * Akarsu  _et al._ (2018) O. Akarsu, R. Brilenkov, M. Eingorn, V. Shulga, and A. Zhuk, Eur. Phys. J. C 78, 609 (2018), arXiv:1806.02669 [gr-qc] .
  * Lazkoz  _et al._ (2018) R. Lazkoz, M. Ortiz-Baños, and V. Salzano, Eur. Phys. J. C 78, 213 (2018), arXiv:1803.05638 [astro-ph.CO] .
  * Nunes  _et al._ (2017) R. C. Nunes, S. Pan, E. N. Saridakis, and E. M. C. Abreu, JCAP 01, 005 (2017), arXiv:1610.07518 [astro-ph.CO] .
  * D’Agostino and Nunes (2020) R. D’Agostino and R. C. Nunes, Phys. Rev. D 101, 103505 (2020), arXiv:2002.06381 [astro-ph.CO] .
  * Wang (2021) D. Wang, Eur. Phys. J. C 81, 482 (2021), arXiv:2008.03966 [astro-ph.CO] .
  * Odintsov  _et al._ (2021) S. D. Odintsov, D. Sáez-Chillón Gómez, and G. S. Sharov, Nucl. Phys. B 966, 115377 (2021), arXiv:2011.03957 [gr-qc] .
  * Amendola  _et al._ (1999) L. Amendola, P. S. Corasaniti, and F. Occhionero, (1999), arXiv:astro-ph/9907222 .
  * Gaztanaga  _et al._ (2002) E. Gaztanaga, E. Garcia-Berro, J. Isern, E. Bravo, and I. Dominguez, Phys. Rev. D 65, 023506 (2002), arXiv:astro-ph/0109299 .
  * Wright and Li (2018) B. S. Wright and B. Li, Phys. Rev. D 97, 083505 (2018), arXiv:1710.07018 [astro-ph.CO] .
  * Perivolaropoulos and Skara (2021) L. Perivolaropoulos and F. Skara, Phys. Rev. D 104, 123511 (2021), arXiv:2109.04406 [astro-ph.CO] .
  * Alestas  _et al._ (2022) G. Alestas, D. Camarena, E. Di Valentino, L. Kazantzidis, V. Marra, S. Nesseris, and L. Perivolaropoulos, Phys. Rev. D 105, 063538 (2022), arXiv:2110.04336 [astro-ph.CO] .
  * Sapone  _et al._ (2021) D. Sapone, S. Nesseris, and C. A. P. Bengaly, Phys. Dark Univ. 32, 100814 (2021), arXiv:2006.05461 [astro-ph.CO] .
  * Kumar  _et al._ (2022) S. Kumar, R. C. Nunes, and P. Yadav, (2022), arXiv:2209.11131 [astro-ph.CO] .
  * Ballardini and Finelli (2022) M. Ballardini and F. Finelli, Phys. Rev. D 106, 063531 (2022), arXiv:2112.15126 [astro-ph.CO] .
  * Faraoni (2006) V. Faraoni, Phys. Rev. D 74, 104017 (2006), arXiv:astro-ph/0610734 .
  * Bean  _et al._ (2007) R. Bean, D. Bernat, L. Pogosian, A. Silvestri, and M. Trodden, Phys. Rev. D 75, 064020 (2007), arXiv:astro-ph/0611321 .
  * Amendola  _et al._ (2007a) L. Amendola, D. Polarski, and S. Tsujikawa, Phys. Rev. Lett. 98, 131302 (2007a), arXiv:astro-ph/0603703 .
  * Chiba  _et al._ (2007) T. Chiba, T. L. Smith, and A. L. Erickcek, Phys. Rev. D 75, 124014 (2007), arXiv:astro-ph/0611867 .
  * Nojiri and Odintsov (2006b) S. Nojiri and S. D. Odintsov, Phys. Lett. B 637, 139 (2006b), arXiv:hep-th/0603062 .
  * Amendola  _et al._ (2007b) L. Amendola, R. Gannouji, D. Polarski, and S. Tsujikawa, Phys. Rev. D 75, 083504 (2007b), arXiv:gr-qc/0612180 .
  * Hu and Sawicki (2007) W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007), arXiv:0705.1158 [astro-ph] .
  * Starobinsky (2007) A. A. Starobinsky, JETP Lett. 86, 157 (2007), arXiv:0706.2041 [astro-ph] .
  * Bamba  _et al._ (2013) K. Bamba, A. Lopez-Revelles, R. Myrzakulov, S. D. Odintsov, and L. Sebastiani, Class. Quant. Grav. 30, 015008 (2013), arXiv:1207.1009 [gr-qc] .
  * Basilakos  _et al._ (2013) S. Basilakos, S. Nesseris, and L. Perivolaropoulos, Phys. Rev. D 87, 123529 (2013), arXiv:1302.6051 [astro-ph.CO] .
  * Starobinsky (1980) A. A. Starobinsky, Phys. Lett. B 91, 99 (1980).
  * Sultana  _et al._ (2022) J. Sultana, M. K. Yennapureddy, F. Melia, and D. Kazanas, Mon. Not. Roy. Astron. Soc. 514, 5827 (2022), arXiv:2206.10761 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020b), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Alam  _et al._ (2021) S. Alam _et al._ (eBOSS), Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 [astro-ph.CO] .
  * Aver  _et al._ (2015) E. Aver, K. A. Olive, and E. D. Skillman, JCAP 07, 011 (2015), arXiv:1503.08146 [astro-ph.CO] .
  * Cooke  _et al._ (2018) R. J. Cooke, M. Pettini, and C. C. Steidel, Astrophys. J. 855, 102 (2018), arXiv:1710.11129 [astro-ph.CO] .
  * Tsujikawa (2007) S. Tsujikawa, Phys. Rev. D 76, 023514 (2007), arXiv:0705.1032 [astro-ph] .
  * Nesseris  _et al._ (2017) S. Nesseris, G. Pantazis, and L. Perivolaropoulos, Phys. Rev. D 96, 023542 (2017), arXiv:1703.10538 [astro-ph.CO] .
  * de la Cruz-Dombriz  _et al._ (2008) A. de la Cruz-Dombriz, A. Dobado, and A. L. Maroto, Phys. Rev. D 77, 123515 (2008), arXiv:0802.2999 [astro-ph] .
  * Orjuela-Quintana and Nesseris (2023) J. B. Orjuela-Quintana and S. Nesseris, (2023), arXiv:2303.14251 [gr-qc] .
  * Bose  _et al._ (2015) S. Bose, W. A. Hellwing, and B. Li, JCAP 02, 034 (2015), arXiv:1411.6128 [astro-ph.CO] .
  * Gannouji  _et al._ (2009) R. Gannouji, B. Moraes, and D. Polarski, Journal of Cosmology and Astroparticle Physics 2009, 034 (2009).
  * Brout  _et al._ (2022) D. Brout _et al._ , Astrophys. J. 938, 110 (2022), arXiv:2202.04077 [astro-ph.CO] .
  * Benisty  _et al._ (2023) D. Benisty, J. Mifsud, J. Levi Said, and D. Staicova, Phys. Dark Univ. 39, 101160 (2023), arXiv:2202.04677 [astro-ph.CO] .
  * Tutusaus  _et al._ (2017) I. Tutusaus, B. Lamine, A. Dupays, and A. Blanchard, Astron. Astrophys. 602, A73 (2017), arXiv:1706.05036 [astro-ph.CO] .
  * Martinelli and Tutusaus (2019) M. Martinelli and I. Tutusaus, Symmetry 11, 986 (2019), arXiv:1906.09189 [astro-ph.CO] .
  * Di Valentino  _et al._ (2020) E. Di Valentino, S. Gariazzo, O. Mena, and S. Vagnozzi, (2020), arXiv:2005.02062 [astro-ph.CO] .
  * Lesgourgues (2011) J. Lesgourgues, “The cosmic linear anisotropy solving system (class) i: Overview,”  (2011).
  * Blas  _et al._ (2011) D. Blas, J. Lesgourgues, and T. Tram, Journal of Cosmology and Astroparticle Physics 2011, 034 (2011).
  * Audren  _et al._ (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302, 001 (2013), arXiv:1210.7183 [astro-ph.CO] .
  * Brinckmann and Lesgourgues (2018) T. Brinckmann and J. Lesgourgues, (2018), arXiv:1804.07261 [astro-ph.CO] .
  * Vázquez  _et al._ (2012) J. A. Vázquez, A. Lasenby, M. Bridges, and M. Hobson, Monthly Notices of the Royal Astronomical Society 422, 1948 (2012), arXiv:1103.4619 [astro-ph.CO] .
  * Kass and Raftery (1995) R. E. Kass and A. E. Raftery, J. Am. Statist. Assoc. 90, 773 (1995).
  * Trotta (2008) R. Trotta, Contemp. Phys. 49, 71 (2008), arXiv:0803.4089 [astro-ph] .
  * Schöneberg  _et al._ (2019) N. Schöneberg, J. Lesgourgues, and D. C. Hooper, JCAP 10, 029 (2019), arXiv:1907.11594 [astro-ph.CO] .
  * Cuceu  _et al._ (2019) A. Cuceu, J. Farr, P. Lemos, and A. Font-Ribera, JCAP 10, 044 (2019), arXiv:1906.11628 [astro-ph.CO] .
  * Schöneberg  _et al._ (2022b) N. Schöneberg, L. Verde, H. Gil-Marín, and S. Brieden, JCAP 11, 039 (2022b), arXiv:2209.14330 [astro-ph.CO] .
  * Efstathiou (2021) G. Efstathiou, Mon. Not. Roy. Astron. Soc. 505, 3866 (2021), arXiv:2103.08723 [astro-ph.CO] .
  * Camarena and Marra (2021) D. Camarena and V. Marra, Mon. Not. Roy. Astron. Soc. 504, 5164 (2021), arXiv:2101.08641 [astro-ph.CO] .
  * Nunes and Di Valentino (2021) R. C. Nunes and E. Di Valentino, Phys. Rev. D 104, 063529 (2021), arXiv:2107.09151 [astro-ph.CO] .
  * Dhawan  _et al._ (2020) S. Dhawan, D. Brout, D. Scolnic, A. Goobar, A. G. Riess, and V. Miranda, The Astrophysical Journal 894, 54 (2020).
  * Alestas  _et al._ (2021) G. Alestas, L. Kazantzidis, and L. Perivolaropoulos, Phys. Rev. D 103, 083517 (2021), arXiv:2012.13932 [astro-ph.CO] .
  * Arjona  _et al._ (2019) R. Arjona, W. Cardona, and S. Nesseris, Phys. Rev. D 99, 043516 (2019), arXiv:1811.02469 [astro-ph.CO] .
  * Philcox  _et al._ (2020) O. H. E. Philcox, M. M. Ivanov, M. Simonović, and M. Zaldarriaga, JCAP 05, 032 (2020), arXiv:2002.04035 [astro-ph.CO] .
  * Chudaykin  _et al._ (2021) A. Chudaykin, K. Dolgikh, and M. M. Ivanov, Phys. Rev. D 103, 023507 (2021), arXiv:2009.10106 [astro-ph.CO] .
  * Nunes  _et al._ (2022) R. C. Nunes, S. Vagnozzi, S. Kumar, E. Di Valentino, and O. Mena, Phys. Rev. D 105, 123506 (2022), arXiv:2203.08093 [astro-ph.CO] .
  * Simon  _et al._ (2023) T. Simon, P. Zhang, V. Poulin, and T. L. Smith, Phys. Rev. D 107, 063505 (2023), arXiv:2208.05930 [astro-ph.CO] .
  * Carrilho  _et al._ (2023) P. Carrilho, C. Moretti, and A. Pourtsidou, JCAP 01, 028 (2023), arXiv:2207.14784 [astro-ph.CO] .
  * Reeves  _et al._ (2023) A. Reeves, L. Herold, S. Vagnozzi, B. D. Sherwin, and E. G. M. Ferreira, Mon. Not. Roy. Astron. Soc. 520, 3688 (2023), arXiv:2207.01501 [astro-ph.CO] .
  * Simon (2022) T. Simon, in  _33rd Rencontres de Blois: Exploring the Dark Universe_ (2022) arXiv:2212.03004 [astro-ph.CO] .
  * Piga  _et al._ (2023) L. Piga, M. Marinucci, G. D’Amico, M. Pietroni, F. Vernizzi, and B. S. Wright, JCAP 04, 038 (2023), arXiv:2211.12523 [astro-ph.CO] .