# Observational constraints on a modified-gravity model with an exponential
function of the curvature using the expansion history, the RSD, and the
Pantheon$+$SH0ES data

Mario A. Acero<EMAIL_ADDRESS>A. Oliveros
<EMAIL_ADDRESS>Programa de Física, Universidad del
Atlántico, Carrera 30 8-49, Puerto Colombia, Atlántico, Colombia

###### Abstract

Considering a well motivated $f(R)$ modified-gravity model, in which an
exponential function of the curvature is included, in this paper we implement
a statistical data analysis to set constraints on the parameters of the model,
taking into account an analytic approximate solution for the expansion rate,
$H(z)$. From the Monte Carlo Markov Chain-based analysis of the expansion rate
evolution, the standardized SN distance modulus and the redshift space
distortion observational data, we find that the preferred value for the
perturbative parameter, $b$, quantifying the deviation of the $f(R)$ model
from $\Lambda$CDM, lives in a region which excludes $b=0$ at $\gtrsim
4.5\sigma$ C.L., and that the predicted current value of the Hubble parameter,
$H_{0}$, locates in between the two observational results currently under
scrutiny from Planck and SH0ES collaborations, indicating that the proposed
model would alleviate the apparent tension. Under the implemented approximate
solution, and with the constraints obtained for the parameters, the proposed
$f(R)$ model successfully reproduce the observational data and the predicted
evolution of interesting cosmological parameters resemble the results of
$\Lambda$CDM, as expected, while an oscillatory behavior of the dark energy
equation of state is observed, pointing to deviation from the concordance
cosmological model. The results presented here reinforces the conclusion that
the $f(R)$ modified-gravity model represents a viable alternative to describe
the evolution of the Universe, evading the challenges faced by $\Lambda$CDM.

###### keywords:

Modified gravity, Dark energy, $f(R)$ gravity, Parameter constraints.

###### pacs:

[

PACS]04.50.Kd, 98.80.-k

##  1 Introduction

Although Einstein’s General Relativity (GR) has been an enormous success in
explaining many observations at the astrophysical and cosmological levels,
there are phenomena that cannot be adequately explained within this framework.
For example, the current observed accelerating expansion of the Universe poses
a challenge [1, 2]. A first attempt to explain this late-time cosmic
acceleration is the introduction of a new energy component in the Universe
known as dark energy (DE), characterized by a negative pressure. However, this
proposal has proven to be very difficult to incorporate within the known
theories of physics (for a comprehensive review about this topic see Refs. [3,
4, 5, 6]). DE is commonly associated with a cosmological constant ($\Lambda$),
which drives the late-time cosmic evolution and whose origins are traced back
to early quantum fluctuations of the vacuum. However, this model (known as
$\Lambda$CDM) faces challenges such as the coincidence and cosmological
constant problems, as well as tensions that have arisen among recent
cosmological measurements.

In order to circumvent the above issues, an interesting proposal is the $f(R)$
gravity theories, in which the Einstein-Hilbert action is modified with a
general function of the Ricci scalar $R$, $f(R)$. Nevertheless, the selection
of a specific $f(R)$ function is not arbitrary: it must adhere to several
consistency requirements and various constraints that impose conditions for
the cosmological viability of $f(R)$ models. One crucial requirement is that a
given $f(R)$ adequately describes the different cosmic eras, including the
radiation, matter, and dark energy eras, and probably the inflationary period.
Moreover, it is imperative also that the selected $f(R)$ function satisfies
both cosmological and local gravity constraints, in addition to other relevant
considerations [7]. Numerous works have been undertaken in this context,
exploring various aspects at both the astrophysical and cosmological levels
[8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31, 32]. In general, for a more extensive review about $f(R)$
theories, interested readers are invited to see Refs. [33, 34, 35].

Within this wide plethora of viable $f(R)$ gravity models, two of those have
been highlighted in the literature: the Hu and Sawicki (HS) model [36] and,
the Starobinsky model [37]. Although these models were originally advertised
as models that do not contain the cosmological constant as part of $f(R)$
being distinct from the $\Lambda$CDM form $f(R)=R-2\Lambda$, in Ref. [38] it
has been demonstrated that these models can be arbitrarily close to
$\Lambda$CDM (where the deviation from $\Lambda$CDM is characterized by a
parameter $b$), and they provide predictions that are similar to those of the
usual (scalar field) DE models, particularly concerning the cosmic history,
including the presence of the matter era, the stability of cosmological
perturbations, the stability of the late de Sitter point, etc. Also, they
found that the parameter $b$ is of the order of $b\sim 0.1$ for the HS model,
thus making it practically indistinguishable from the $\Lambda$CDM at the
background level.

In a related investigation, in Ref. [39], the authors introduce a new class of
models that are variants of the HS model that interpolate between the
cosmological constant model and a matter dominated universe for different
values of the parameter $b$, which is usually expected to be small for viable
models and which, in practice, measures the deviation from GR. Recently, in
Ref. [40], the state-of-the-art BAO+BBN data and the most recent Type Ia
supernovae (SNe Ia) sample, PantheonPlus, including the Cepheid host distances
and covariance from SH0ES samples, were used to robustly constrain the HS and
Starobinsky models, and found that both models are consistent with GR at a 95%
CL.

As a contribution to the research in this matter, in this work, we analyze the
parameters governing the $f(R)$ model proposed in [41] with the approximate
analytical solution found by one of the authors [42] (Section 2), setting
constraints on the characteristic parameter, $b$, and on some cosmological
parameters as predicted by the studied model. The constrains are obtained by
performing a statistical analysis based on the Markov Chain Monte Carlo (MCMC)
method (Section 3), and considering three sets of observational data: the
Hubble parameter ($H(z)$), the Type Ia Supernova sample (Pantheon+SH0ES), and
the redshift distortion sample ($f\sigma_{8}$). We find that, although posing
constraints on the model parameters presents some difficulties when individual
datasets are considered, the joint statistical analyses allows to set strong
constraints on the parameters, such that the model fits the data accurately,
within the observational uncertainties. In addition, our model predictions for
the considered cosmological parameters are found to be consistent with those
reported by Planck or DESI (within the $\sim 1\sigma$ C.L.). Remarkably, we
obtain a present value of the Hubble parameter, $H_{0}$, laying between the
values reported by Planck [43] and SH0ES [44, 45], alleviating the tension
between these observations.

Our results also indicate that the value of deviation parameter $b$ that best
fit the data (Section 3) is larger that expected, considering the perturbative
approach implemented in [42] to find the approximate solution. The impact of
such a large value reflects on the redshift-evolution of the cosmological
parameters $w_{\rm{eff}}$, $q$, $w_{\rm{DE}}$ and $\Omega_{\rm{DE}}$ presented
in Section 4, with particular impact on $w_{\rm{DE}}$, from which an
oscillatory evolution at late times is obtained.

We also present results from statistical (Section 3.6) and dynamics (4) tests
to verify the validity of the model under consideration, and address the
conclusions in Section 5.

##  2 $f(R)$ Gravity: Preliminaries

In general, the gravitational action of $f(R)$ gravity in the presence of
matter components is given by

$S=\int{d^{4}x\sqrt{-g}\left(\frac{f(R)}{2\kappa^{2}}+\mathcal{L}_{\rm{M}}\right)},$
(1)

where $g$ denotes the determinant of the metric tensor $g^{\mu\nu}$,
$\kappa^{2}=8\pi G=1/M_{\rm{p}}^{2}$, with $G$ being the Newton’s constant and
$M_{\rm{p}}$ the reduced Planck mass. $\mathcal{L}_{\rm{M}}$ represents the
Lagrangian density for the matter components (relativistic and non-
relativistic perfect matter fluids). The term $f(R)$ is for now an arbitrary
function of the Ricci scalar $R$. Variation with respect to the metric gives
the equation of motion

$f_{R}(R)R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}f(R)+(g_{\mu\nu}\square-\nabla_{\mu}\nabla_{\nu})f_{R}(R)=\kappa^{2}T_{\mu\nu}^{(\rm{M)}},$
(2)

where $f_{R}\equiv\frac{df}{dR}$, $\nabla_{\mu}$ is the covariant derivative
associated with the Levi-Civita connection of the metric, and
$\square\equiv\nabla^{\mu}\nabla_{\mu}$. Plus, $T_{\mu\nu}^{(\rm{M})}$ is the
matter energy–momentum tensor which is assumed to be a perfect fluid.
Considering the flat Friedman-Robertson-Walker (FRW) metric,

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

with $a(t)$ representing the scale factor, the time and spatial components of
Eq. (2) are given, respectively, by

$3H^{2}f_{R}=\kappa^{2}(\rho_{\rm{m}}+\rho_{\rm{r}})+\frac{1}{2}(Rf_{R}-f)-3H\dot{f}_{R},$
(4)

and

$-2\dot{H}f_{R}=\kappa^{2}\left(\rho_{\rm{m}}+\frac{4}{3}\rho_{\rm{r}}\right)+\ddot{f}_{R}-H\dot{f}_{R},$
(5)

where $\rho_{\rm{m}}$ is the matter density and $\rho_{\rm{r}}$ denotes the
density of radiation. The over-dot denotes a derivative with respect to the
cosmic time $t$ and $H\equiv\dot{a}/a$ is the Hubble parameter. Note that in
the spatially flat FLRW Universe, the Ricci scalar $R$ takes the form

$R=6(2H^{2}+\dot{H}).$ (6)

If there is no interaction between non-relativistic matter and radiation, then
these components obey separately the conservation laws:

$\dot{\rho}_{\rm{m}}+3H\rho_{\rm{m}}=0,\quad\dot{\rho}_{\rm{r}}+4H\rho_{\rm{r}}=0.$
(7)

As usual in the literature, it is possible to rewrite the field equations (4)
and (5) in the Einstein-Hilbert form:

$3H^{2}=\kappa^{2}\rho,$ (8) $-2\dot{H}^{2}=\kappa^{2}(\rho+p),$ (9)

where $\rho=\rho_{\rm{m}}+\rho_{\rm{r}}+\rho_{\rm{DE}}$ and
$p=p_{\rm{m}}+p_{\rm{r}}+p_{\rm{DE}}$ correspond to the total effective energy
density and total effective pressure density of the cosmological fluid. In
this case, the dark energy component has a geometric origin, and after a some
manipulation in Eqs. (4) and (5), we obtain the effective dark energy and
pressure corresponding to $f(R)$-theory given by

$\rho_{\rm{DE}}=\frac{1}{\kappa^{2}}\left[\frac{Rf_{R}-f}{2}+3H^{2}(1-f_{R})-3H\dot{f}_{R}\right],$
(10)

and

$p_{\rm{DE}}=\frac{1}{\kappa^{2}}[\ddot{f}_{R}-H\dot{f}_{R}+2\dot{H}(f_{R}-1)-\kappa^{2}\rho_{\rm{DE}}],$
(11)

it is easy to show that $\rho_{\rm{DE}}$ and $p_{\rm{DE}}$ defined in this way
satisfy the usual energy conservation equation

$\dot{\rho}_{\rm{DE}}+3H(\rho_{\rm{DE}}+p_{\rm{DE}})=0,$ (12)

in this case, we assume that the equation of state (EoS) parameter for this
effective dark energy satisfies the relation
$w_{\rm{DE}}=p_{\rm{DE}}/\rho_{\rm{DE}}$, and in explicit form it is given by

$w_{\rm{DE}}=-1-\frac{H\dot{f}_{R}+2\dot{H}(1-f_{R})-\ddot{f}_{R}}{\frac{1}{2}(f_{R}R-f)-3H\dot{f}_{R}+3(1-f_{R})H^{2}},$
(13)

In the following sections, our analysis will be focused on the $f(R)$ gravity
model, defined by

$f(R)=R-2\,\Lambda\,e^{-(b\Lambda/R)^{n}},$ (14)

where $b$ and $n$ are positive real dimensionless parameters, and $\Lambda$ is
the cosmological constant. This model was introduced in Ref. [42], and it is a
reparametrization of a specific viable $f(R)$ gravity model studied in Refs.
[46, 41]. Furthermore, it is shown that the HS model is a limiting case of
this model. In the literature, other authors have studied some $f(R)$ gravity
models with exponential functions of the scalar curvature (see for example
Refs. [47, 48, 49, 50, 51]).

From the specific form of this model, and as has been demonstrated in Ref.
[38], it is possible to derive an analytic approximation for the expansion
rate $H(z)$. This approximate analytical expression was found by one of the
authors in Ref. [42], and it is given by

$\displaystyle E^{2}(z)$
$\displaystyle\equiv\frac{H^{2}(z)}{H_{0}^{2}}=1-\Omega_{m0}+(1+z)^{3}\Omega_{m0}$
(15)
$\displaystyle+\frac{6b(\Omega_{m0}-1)^{2}\left(-4+\Omega_{m0}(9-3\Omega_{m0}+z(3+z(z+3))(1+(3+2z(3+z(z+3)))\Omega_{m0}))\right)}{(4+(-3+z(3+z(z+3)))\Omega_{m0})^{3}}$
$\displaystyle+\frac{b^{2}(\Omega_{m0}-1)^{3}}{(1+z)^{24}\left(\frac{4(1-\Omega_{m0})}{(1+z)^{3}}+\Omega_{m0}\right)^{8}}\bigg{[}5120(\Omega_{m0}-1)^{6}+9216(1+z)^{3}(\Omega_{m0}-1)^{5}\Omega_{m0}$
$\displaystyle-30144(1+z)^{6}(\Omega_{m0}-1)^{4}\Omega_{m0}^{2}+31424(1+z)^{9}(\Omega_{m0}-1)^{3}\Omega_{m0}^{3}-9468(1+z)^{12}$
$\displaystyle\times(\Omega_{m0}-1)^{2}\Omega_{m0}^{4}-4344(1+z)^{15}(\Omega_{m0}-1)\Omega_{m0}^{5}+\frac{37}{2}(1+z)^{18}\Omega_{m0}^{6}\bigg{]},$

where for simplicity, it has been assumed that $\Omega_{r0}=0$, $n=1$ and made
the substitution $N=-\ln{(1+z)}$.

##  3 Cosmological constraints

This section is devoted to the description of the performed statistical
analysis and the considered observational data, to evaluate constraints on the
free parameter of the model, $b$, as well as on some of the relevant
cosmological parameters, as predicted from the $f(R)$ model. We also present a
comparison with the predictions of the $\Lambda$CDM model when the same
statistical analysis and datasets are considered.

The statistical analysis used here is based on the well known Markov Chain
Monte Carlo (MCMC) method implemented with the emcee package [52] to find the
parameters that maximize a user-defined likelihood function

$\mathcal{L}(D|z;\bm{\theta})=-\ln
p(D|z;\bm{\theta})=-\frac{1}{2}\chi^{2}(z;\bm{\theta}),$ (16)

where $D$ refers to the analyzed dataset(s), $\bm{\theta}$ is the vector of
free the parameters to fit (the actual elements of this vector depend on the
dataset under consideration, as explained in the following sections), and $z$
is the independent variable which, for our case corresponds to the redshift.

For each dataset considered in this paper (the Hubble parameter observational,
the Type Ia supernova –Pantheon+– and the redshift space distortion (RSD)
data), a suitable $\chi^{2}$ function is defined, considering the particular
number of data and the observed uncertainties. In addition, combinations of
the different datasets are also considered, looking for strengthen the
constraints on the relevant parameters, in which case the corresponding
$\chi^{2}$ function would be the sum of the individual functions for each
dataset, i.e.,

$\chi_{\rm{tot}}^{2}(z;\bm{\theta})=\sum_{i}\chi_{i}^{2}(z;\bm{\theta}),\quad
i=(H(z),{\rm Pantheon},f\sigma_{8}).$ (17)

Using the MCMC method benefits by the inclusion of any information previously
known about the parameters. This is done by adding suitable _priors_ which
make the emcee package to explore the parameters inside a defined range, with
an specified probability distribution. In order to avoid any possible bias on
the analysis, flat priors are enforced for all the parameters, with the
corresponding ranges shown in table 1.

Table 1: Defined ranges for the parameters to fit, included as flat priors in the MCMC analysis. Parameter |  $H_{0}$111$H_{0}$ is measured in km/s/Mpc. | $\Omega_{m0}$ | $\sigma_{8}$ | $b$
---|---|---|---|---
Range | (60, 80) | (0.1, 0.4) | (0.6, 1.0) | (0, 2)
|  |  |  |

As can be noticed from Table 1, we allow $H_{0}$ to float in an interval which
includes the two values that are currently under discussion given the
independent measurements by Planck, $H_{0}^{\rm{Planck}}=[67.36\pm 0.54]$ km
s-1 Mpc-1 [43], and SH0ES, $H_{0}^{\rm{SH0ES}}=[73.30\pm 1.04]$ km s-1 Mpc-1
[44] (see also Ref. [45] for a reported value with reduced uncertainty). By
doing so, we are able to test if our model shows indications for alleviating
the tension between the two observations. Moreover, it is noteworthy that we
have considered a range of positive values for $b$, as this choice aligns with
the imposed conditions of the current model, i.e. $f_{R}>0$ and $f_{RR}>0$ for
$R>R_{0}$ $(>0)$ (where $R_{0}$ is the Ricci scalar at the present time), and
also $0<\frac{Rf_{RR}}{f_{R}}(r)<1$ at the de Sitter point,
$r=-\frac{Rf_{R}}{f}=-2$.

###  3.1 The Hubble parameter data

For the observational Hubble parameter, we consider the data reported in [53]
based on cosmic chronometers (CC) and radial Baryon Acoustic Oscillations
(BAO) methods (see Ref. [53] for additional details). In this case, the
$\chi^{2}$ used for the likelihood maximization is defined as

$\chi^{2}_{Hz}(z;\bm{\theta})=\sum_{k=1}^{N_{d}}\frac{\left(H_{{\rm
fR},k}(z;\bm{\theta})-H_{{\rm obs},k}(z)\right)^{2}}{\sigma_{k}^{2}},$ (18)

with $N_{d}=40$, and $\bm{\theta}=(H_{0},\Omega_{m0},b)$. Here, $H_{{\rm
fR},k}$ and $H_{{\rm obs},k}$ are the $f(R)$-model prediction (Eq. (15)), and
the observed values of the expansion rate, respectively, and $\sigma_{k}$ is
the corresponding observational error. An analogous $\chi^{2}$ functions is
used to test the $\Lambda$CDM model prediction for $H(z)$, for which the
vector of parameters reduces to
$\bm{\theta}_{\Lambda\rm{CDM}}=(H_{0},\Omega_{m0})$.

The results of this fit are shown in Fig. 1, where the $1\sigma$– and
$2\sigma$–confidence contour plots and posterior probabilities for the
considered parameters are exhibited (purple contours and lines). The best fit
(BF) values of the parameters are also presented in Table 2, where we also
write our results obtained from the fitting this dataset to the $\Lambda$CDM
model, and the reported values by the Planck [43] and DESI [54]
Collaborations.

Figure 1: Contour plots and 1-D posterior probabilities obtained from the MCMC
analysis of the $H(z)$ (purple) and the Pantheon (dark blue) observational
data, as well as the combination those two datasets (red), for the parameters
$\left(H_{0},\Omega_{m0},b\right)$. The numbers over the 1-D posteriors
correspond to the joint analysis.

Regarding the present value of the Hubble parameter, $H_{0}$, within the
$1\sigma$ interval, our model shows an indication of reducing the tension
between the Planck and SH0ES observations, being also consistent with our
$\Lambda$CDM model analysis. Our model prediction for $\Omega_{m0}$ (which is
consistent with our fit to the $\Lambda$CDM model) agrees with the Planck
value within a ${\sim}1.4\sigma$ C.L. On the other hand, although the
parameter $b$ is not strongly constrained by this dataset, the prediction is
consistent with $b=0$, and the allowed interval expands up to $b\lesssim 0.8$
at $\,\sim 2\sigma$ C.L.

###  3.2 The standardized distance modulus - Type Ia Supernova Sample

For the Ia Supernova distance modulus we consider the Pantheon+SH0ES (referred
to as Pantheon here) database described in Refs. [55, 44], comprising 1701
data points in a range of $0.001\leq z\leq 2.3$. The analysis was performed
with a suited $\chi^{2}$ function, considering both statistical and systematic
uncertainties through a covariance matrix, $\bm{C}_{\rm{cov}}$:

$\chi^{2}_{\rm{Pantheon}}=\left[\mu_{\rm{fR}}(z;\bm{\theta})-\mu_{\rm{obs}}(z)\right]^{T}\bm{C}_{\rm{cov}}^{-1}\,\left[\mu_{\rm{fR}}(z;\bm{\theta})-\mu_{\rm{obs}}(z)\right].$
(19)

Both, the covariance matrix and the observed distance modulus $\mu_{\rm{obs}}$
were obtained from the Pantheon+SH0ES data release [56]. For the model
prediction, we have

$\mu_{\rm{fR}}(z;\bm{\theta})=M+25+5\log_{10}D_{L}(z;\bm{\theta}),$ (20)

where $M$ is the absolute magnitude, the parameters vector is
$\bm{\theta}=(H_{0},\Omega_{m0},b)$, and $D_{L}(z;\bm{\theta})$ is the
luminosity distance given by

$D_{L}(z;\bm{\theta})=c\,(1+z)\int_{0}^{z}\frac{dz^{\prime}}{H_{\rm{fR}}(z^{\prime};\bm{\theta})},$
(21)

where $c$ is the speed of light and, as before,
$H_{\rm{fR}}(z^{\prime};\bm{\theta})$ is calculated using Eq. (15). As for the
analysis of the expansion rate in the previous section, here we also carry out
the fit of the $\Lambda$CDM model prediction for the distance modulus with the
same $\chi^{2}$ function, Eq. (19), reducing the parameters vector to
$\bm{\theta}=(H_{0},\Omega_{m0})$.

The contour plots, together with the posterior probabilities for the fitted
parameters for the $f(R)$ model are shown in Fig. 1 (dark blue contours and
lines), while the corresponding best fit values are presented in Table 2. As
clearly exhibited by the (dark blue) contours, there is a strong correlation
among the parameters, producing long allowed areas, covering most of studied
range of values.

In particular, for the present value of the Hubble expansion rate, $H_{0}$, a
statistically weak preference for lower values (in a better agreement with
Planck) is obtained, pushing also $b$ to be small, but different from zero up
to almost $\,\sim 2\sigma$ C.L. (see Table 2). In this scenario (low $H_{0}$
and small $b$), the allowed region for $\Omega_{\rm{m}0}$ extends to largest
explored values, in good agreement with the Planck and DESI observations.

###  3.3 Combining the Hubble expansion rate and the Pantheon+SH0ES datasets

As anticipated at the beginning of this section, a joint analysis considering
the two previously described datasets was also implemented, adding the
corresponding $\chi^{2}$-functions, Eqs. (18) and (19). The combined fit
produces the expected results, shown as (red) contours and 1-D posterior
probabilities in Fig. 1, and in Table 2 (see the row for $H(z)+$SN). First, it
is clear that, despite the larger number of data available from the Pantheon
sample, the constraints on the parameters are dominated by the $H(z)$ dataset;
this is particularly apparent from the 1-D histograms for $H_{0}$ and
$\Omega_{\rm{m}0}$. For these quantities, the joint fit keeps the model
prediction close to that obtained from the $H(z)$ dataset alone, but enhancing
the corresponding limits (i.e. reducing the allowed regions).

This last feature is exceptionally noticeable for the $b$ parameter, for which
not only we obtain a rather large preferred value, but also $b=0$ is excluded
at more that $3\sigma$ C.L. Though one would expect the deviation parameter
$b$ to be close to zero, bringing our model close to $\Lambda$CDM, as it will
be clear later (Sec. 3.5), the proposed $f(R)$ model with a large deviation
parameter successfully fit the considered observational data. In addition,
this results agree with earlier studies [38], where values of $b$ of order
$\mathcal{O}(1)$ are also obtained.

###  3.4 The redshift space distortion, $f\sigma_{8}$ \- The growth Sample

The last dataset considered here is the value of the growth rate $f(z)$
multiplied by the amplitude of the matter power spectrum on the scale of
$8h^{-1}\rm{Mpc}$, $\sigma_{8}(z)$, usually written as $f\sigma_{8}(z)$. This
quantity is considered the best observable to discriminate between modified
gravity theories (such as $f(R)$ gravity models) and $\Lambda$CDM, given that
many $f(R)$ gravity models are virtually indistinguishable from the
$\Lambda$CDM model at the background level [57]. We consider a total of
$N_{d}=26$ data points for different redshifts, $0.013\leq z\leq 1.944$ [39,
58, 59], with the $\chi^{2}$ function defined as

$\chi^{2}_{f\sigma_{8}}=\sum_{k=1}^{N_{d}}\frac{\left[(f\sigma_{8})_{\rm{fR},k}(z;\bm{\theta})-(f\sigma_{8})_{\rm{obs},k}(z)\right]^{2}}{\sigma_{k}^{2}}.$
(22)

As for the previous cases, $\sigma_{k}$ is the corresponding error for each
observational value $(f\sigma_{8})_{\rm{obs},k}(z)$, which is compared against
the model prediction, $(f\sigma_{8})_{\rm{fR},k}(z;\bm{\theta})$, with
$\bm{\theta}=(\sigma_{8},H_{0},\Omega_{m0},b)$. The predicted growth rate is
computed though the following relation [39, 60]:

$(f\sigma_{8})_{\rm{fR}}(z;\bm{\theta})=\sigma_{8}\frac{\delta_{\rm{m}}^{\prime}(z;\bm{\theta})}{\delta_{\rm{m}}(z=0)},$
(23)

where $\sigma_{8}=\sigma_{8}(z=0)$,
$\delta_{\rm{m}}\equiv\delta\rho_{\rm{m}}/\rho_{\rm{m}}$ is the gauge-
invariant matter density perturbation (the density contrast), and the prime
stands for the derivative with respect to the redshift, $z$. To obtain the
theoretical prediction from Eq. (23), it is necessary to calculate
$\delta_{\rm{m}}$. The equation governing the evolution of this quantity for
the $f(R)$ gravity has been derived previously in the literature, considering
the subhorizon approximation ($k^{2}/a^{2}\gg H^{2}$) [61, 62], and it has the
following form:

$\ddot{\delta}_{\rm{m}}+2H\dot{\delta}_{\rm{m}}-4\pi
G_{\rm{eff}}(a,k)\rho_{\rm{m}}\delta_{\rm{m}}=0;$ (24)

here the dot denotes the differentiation with respect to the cosmic time,
$G_{\rm{eff}}(a,k)$ is the effective gravitational “constant”, $k$ is the
comoving wave number, $a$ is the scale factor normalized to unity at present
epoch, and $\rho_{\rm{m}}$ is the background matter density. In order to
facilitate our calculations, we rewrite Eq. (24) in terms of $z$, as follows:

$\delta_{\rm{m}}^{\prime\prime}(z)+\left(\frac{E^{2\,\prime}(z)}{2E^{2}(z)}-\frac{1}{1+z}\right)\delta_{\rm{m}}^{\prime}(z)-\frac{3\Omega_{\rm{m0}}}{2E^{2}(z)}(1+z)\frac{G_{\rm{eff}}(z,k)}{G_{\rm{N}}}\delta_{\rm{m}}(z)=0;$
(25)

in this case, the explicit form for $G_{\rm{eff}}(z,k)$ is

$G_{\rm{eff}}(z,k)=\frac{G_{\rm{N}}}{f_{R}}\left[1+\frac{k^{2}(1+z)^{2}(f_{RR}/f_{R})}{1+3k^{2}(1+z)^{2}(f_{RR}/f_{R})}\right],$
(26)

where $G_{\rm{N}}$ is the Newton constant. Eq. (25) has been expressed in
terms of $E^{2}(z)$, since this function is known in explicit form in our
case. To solve Eq. (25) numerically, we adopt initial conditions for the
density contrast, and its first derivative that are consistent with those
observed at very high redshifts (matter era), matching that of the
$\Lambda$CDM model.

The statistical analysis allows us to set constraint to the parameters,
$\bm{\theta}=(\sigma_{8},H_{0},\Omega_{m0},b)$. However, as observed in Table
2 (see the row for $f(R)-f\sigma_{8}$), the allowed interval obtained for
$H_{0}$ is considerably large (and larger than the result of the other fits);
in fact, the constraints on this parameter are well weaker than from the other
cases, indicating that these data alone are not enough to provide a robust
fit.

To overcome this situation, we performed a series of joint fits combining the
growth sample with the Hubble expansion rate sample, with the Pantheon sample
only, and with both samples at the same time. The corresponding posterior
distributions for the considered parameters resulting from these three
different analyses are shown in Fig. 2, where the numbers on top of each
column correspond to the inferred values from the combination of the three
datasets (gray 2D-contours and histograms). Also, Table 2 exhibit the
intervals at a 68% C.L. for all of the fits.

Figure 2: Contour plots and 1-D posterior probabilities obtained from the MCMC
analysis, for the combination of $H(z)$, $f\sigma_{8}(z)$, and Pantheon
observational data for the parameters
$\left(\sigma_{8},H_{0},\Omega_{m0},b\right)$.

Looking at Fig. 2 one can notice the apparent effect of combining the three
datasets. As expected, all the parameters are better bounded in this case, and
the most evident impact is on $H_{0}$ (which, as mentioned earlier, is not
constrained by the $f\sigma_{8}$ data alone), for which the combined data set
a well constrained allowed region around $H_{0}=70.7$ km/s/Mpc. Interestingly,
bounds on $\sigma_{8}$ are also improved with the combination, showing that,
although this parameter is undoubtedly not limited by the $H(z)$ data, the
joint fit with Pantheon and the growth rate samples makes the model to predict
a smaller allowed region for $\sigma_{8}$, rejecting values outside
$(0.81,0.87)$ (see lower-left panel of Fig. 2).

Table 2: Resulting $1\sigma$ allowed intervals for the fitted parameters from the MCMC analysis. The values in the first (second) row correspond to those for the $\Lambda$CDM model from Planck 2018 [43] (DESI 2024 [54]); the following rows are the result of fitting the $\Lambda$CDM model predictions to the corresponding dataset, while the last part of the table shows the corresponding predictions from our $f(R)$ model. Note that, here, $\rm{SN}$ refers to the Pantheon dataset. Model |  $H_{0}$111$H_{0}$ is measured in km/s/Mpc. | $\Omega_{m0}$ | $\sigma_{8}$ | $b$
---|---|---|---|---
$\Lambda$CDM | Planck [43] | $67.4\pm 0.5$ | $0.315\pm 0.007$ | $0.811\pm 0.006$ | –
$\Lambda$CDM | DESI [54] |  $68.52\pm 0.62$222DESI BAO. |  $0.295\pm 0.015$333DESI BAO + CMB. |  $0.8135\pm 0.0053$444DESI BAO + Planck[plik] + CMB lensing. | –
$\Lambda$CDM | $H(z)$ | $73.1_{-2.8}^{+2.7}$ | $0.237_{-0.025}^{+0.029}$ | – | –
$\rm{SN}$ | $72.8\pm 0.2$ | $0.361\pm 0.018$ | – | –
$H(z)+\rm{SN}$ | $73.9\pm 0.2$ | $0.263\pm 0.009$ | – | –
$f\sigma_{8}$ | – | $0.294_{-0.042}^{+0.048}$ | $0.805\pm 0.033$ | –
$H(z)+f\sigma_{8}$ | $71.6\pm 2.5$ | $0.254_{-0.023}^{+0.026}$ | $0.832\pm 0.028$ | –
${\rm SN}+f\sigma_{8}$ | $72.9\pm 0.2$ | $0.353\pm 0.017$ | $0.777\pm 0.022$ | –
$H(z)+{\rm SN}+f\sigma_{8}$ | $73.8\pm 0.2$ | $0.264\pm 0.009$ | $0.824_{-0.022}^{+0.023}$ | –
$f(R)$ | $H(z)$ | $71.0_{-2.6}^{+2.3}$ | $0.246_{-0.026}^{+0.027}$ | – | $0.342_{-0.234}^{+0.315}$
SN | $67.2_{-2.6}^{+5.4}$ | $0.254_{-0.074}^{+0.068}$ | – | $0.629_{-0.337}^{+0.315}$
$H(z)+{\rm SN}$ | $70.4_{-1.2}^{+1.6}$ | $0.251_{-0.024}^{+0.023}$ | – | $0.615_{-0.144}^{+0.136}$
$f\sigma_{8}$ | $70.1_{-6.8}^{+6.7}$ | $0.218_{-0.069}^{+0.071}$ | $0.863_{-0.052}^{+0.064}$ | $0.801_{-0.441}^{+0.320}$
$H(z)+f\sigma_{8}$ | $71.1_{-2.1}^{+3.3}$ | $0.244_{-0.026}^{+0.025}$ | $0.841_{-0.028}^{+0.030}$ | $0.440\pm 0.268$
${\rm SN}+f\sigma_{8}$ | $77.0_{-2.1}^{+2.0}$ | $0.256_{-0.030}^{+0.045}$ | $0.834_{-0.035}^{+0.031}$ | $0.609_{-0.226}^{+0.137}$
$H(z)+{\rm SN}+f\sigma_{8}$ | $70.7_{-1.2}^{+1.3}$ | $0.245_{-0.023}^{+0.022}$ | $0.844\pm 0.028$ | $0.648_{-0.137}^{+0.125}$
|  |  |  |  |

Notice how the combination of the Pantheon and the growth rate samples (red
2D-contorus and histograms) result in a strong correlation among the
parameters, which was also evident before, from the analysis of the Pantheon
sample alone (see Fig. 1). Adding the growth rate sample improves the
constraints by reducing the allowed regions, but also pushes $H_{0}$ towards
larger values, closer to the SH0ES [44] observation than to the Planck [43]
one, as one would actually expect. Importantly, the inclusion of the Hubble
parameter sample (gray 2D-contorus and histograms in Fig. 2) further improves
the constraints over all the parameter, particularly pushing $H_{0}$ back to
low values, alleviating the tension between Planck and SH0ES. Also, although
still evident in the 2D-contours, the correlation between the parameters is
reduced. We note here that, in the case of the $(\sigma_{8},\Omega_{\rm{m}0})$
space, a similar result was obtained by the authors of Ref. [39], in a context
were variants of the Hu-Sawicki model were studied. Comparing this
2D-parameter space with a more recent study [63], were constraints from the
redshift-space galaxy skew spectra are set for some cosmological parameters
(although not in the context of $f(R)$ gravity models), we see that we obtain
compatible results both, at the 2D-contour level, and the allowed region for
each parameter. This provides an interesting insight about the possibility of
strengthen our constrains even further by the inclusion of the non-Gaussian
information of the cosmic large-scale structure, a task which might be
considered in a future work.

Another relevant aspect to remark here is the results for the perturbative
parameter, $b$. As noticed in the previous section (see Fig. 1, and Tab. 2),
when only the $H(z)$ sample is considered, $b$ allowed region goes to lower
values, including $b=0$; however, when the other two samples are considered
into the analysis, not only the constraints are strengthened, but also $b$
preferred values are moved to rather large values, now excluding $b=0$ at
$\gtrsim 4.5\sigma$. These large values for $b$ might not be awaited,
considering its perturbative nature, but it is not totally unexpected since
similar results have been observed before [38, 64] (see also [39], were a
particular degenerate hypergeometric model was considered, obtaining a
$b\approx 6$ best fit.)

Let us conclude this section by pointing that even if it would be natural to
expect $b$ to be close to zero, our statistical analysis indicate that this is
not quite the case for the considered model and data samples. One must notice,
though, that the large value obtained for the perturbative parameter is
statistically strong and the fit to the data is compensated (see Sec. 3.5,
bellow) by the other relevant cosmological parameters considered in the
analysis, $(H_{0},\Omega_{\rm{m}0},\sigma_{8})$, which are considerably
different to those reported by Planck and DESI (see Table 2), within the
$\Lambda$CDM model. Indeed, when we perform the same statistical analysis
implemented for our $f(R)$ model to the $\Lambda$CDM predictions, the results
appear to be more compatible between these two models (see the $\Lambda$CDM
and $f(R)$ set of rows in Table 2), with a remarkable distinction coming from
the fact that our proposed $f(R)$ model (i.e., with $b\neq 0$) predicted
$H_{0}$ value indicates an alleviation of the Planck-SH0ES measurement
tension.

###  3.5 Model predictions vs. observational data

As an important evaluation of the results presented in the previous sections,
the obtained values for the constrained parameters
($H_{0},\Omega_{\rm{m}0},\sigma_{8},b$) are used to draw the evolution of the
Hubble expansion rate, the distance modulus Ia-SN, and the space distortion,
in terms of the redshift, $z$, as predicted by our $f(R)$ model.

This is shown in Figs. 3-5, where the $f(R)$ predictions (red dash-dotted
line) are compared with each sample published data (black dots with the
vertical line indicating the corresponding data uncertainty), as well as with
the predictions from the $\Lambda$CDM model, considering both, Planck [43]
(blue line) and DESI [44] (black dashed line) reported observations.

Figure 3: Evolution of the Hubble parameter with the redshift, $z$, as
predicted for the $f(R)$ model presented in this work (dot-dashed red),
compared against observational data (black dots with the vertical lines
indicating the uncertainty). The prediction by the $\Lambda$CDM model (full
blue for Planck and dashed black for DESI) is also shown for comparison.
Figure 4: As in Fig. 3, but for the corrected/standardized magnitude as a
function of the redshift. Figure 5: As in Fig. 3, but for the growth rate as a
function of the redshift.

For a complete and more consistent comparison, in Figs. 3-5 we also included
our own fit of the $\Lambda$CDM model predictions to the considered data
samples (blue dotted line). Although not easily visible in Fig. 4, the three
figures also exhibit a light red band obtained thought allowint the parameters
float up to the $1\sigma$ allowed intervals. It is evident that our model very
well reproduces the observations and that, despite the large value of the
perturbative parameter $b$, the proposed model does not deviates considerably
from $\Lambda$CDM.

###  3.6 Information Criteria

In this section we implement a different evaluation of the fits described in
the previous sections, by the using two standard information criteria (IC):
the Akaike Information Criterion (AIC) and the Bayesian Information Criterion
(BIC). This procedure provides a way to compare a set of model with their
predictions given by datasets (see Ref. [65] and references therein for a
complete description, and Refs. [66, 64], where this analysis is also
implemented). This analysis is useful to compare models with different number
of parameters and the number of data points for the different data samples
under consideration.

Specifically, the AIC estimator is given by [65]

$\rm{AIC}=-2\ln(\mathcal{L}_{\rm{max}})+2k+\frac{2k(k+1)}{N_{\rm{tot}}-k-1},$
(27)

while the BIC evidence estimator is computed through

$\rm{BIC}=-2\ln(\mathcal{L}_{\rm{max}})+k\log(N_{\rm{tot}}),$ (28)

where $k$ is the number of free parameters in the proposed model,
$\mathcal{L}_{\rm{max}}$ is the maximum likelihood value of the dataset(s)
considered for analysis, and $N_{\rm{tot}}$ is the number of data points.
Then, to compare the models, we compute the relative differences between the
IC,

$\Delta IC_{\rm{model}}=IC_{\rm{model}}-IC_{\rm{min}},$ (29)

where $IC_{\rm{min}}$ is the minimum value of IC of the set of competing
models [65]. According to the authors of Ref. [65], a value $\Delta IC\leq 2$
indicates the statistical compatibility of the compared models; obtaining
$2<\Delta IC<6$ points to a moderate tension between the models, and $\Delta
IC\geq 10$ hints towards a strong tension. In general, the larger the $\Delta
IC_{\rm{model}}$, the stronger the evidence against the model as compared with
the model with $IC_{\rm{min}}$.

Table 3: Results of the information criteria analyses comparing the $\Lambda$CDM model with the $f(R)$ model proposed in this work, with the different sets of considered data samples. Dataset | Model | $\chi^{2}_{\rm{min}}$ | AIC |  $|\Delta$AIC$|$ | BIC |  $|\Delta$BIC$|$
---|---|---|---|---|---|---
$H(z)$ |  $\Lambda$CDM | 18.69 | 23.02 | 0 | 26.07 | 0
$f(R)$ | 19.09 | 25.75 | 2.73 | 30.15 | 4.08
Pantheon |  $\Lambda$CDM | 1752.51 | 1756.52 | 0 | 1767.39 | 0
$f(R)$ | 1750.48 | 1756.49 | 0.03 | 1772.79 | 5.41
$H(z)$ \+ Pantheon |  $\Lambda$CDM | 1816.37 | 1820.38 | 0 | 1831.29 | 0
$f(R)$ | 1771.20 | 1777.21 | 43.17 | 1793.58 | 37.71
$f\sigma_{8}$ |  $\Lambda$CDM | 14.92 | 19.45 | 0 | 21.44 | 0
$f(R)$ | 13.22 | 20.31 | 0.86 | 22.99 | 1.55
$H(z)+f\sigma_{8}$ |  $\Lambda$CDM | 34.85 | 41.24 | 0 | 47.42 | 0
$f(R)$ | 33.41 | 42.07 | 0.83 | 50.17 | 2.75
$f\sigma_{8}$ \+ Pantheon |  $\Lambda$CDM | 1769.20 | 1775.22 | 0 | 1791.57 | 0
$f(R)$ | 1764.29 | 1772.32 | 2.90 | 1794.11 | 2.54
$H(z)$ \+ $f\sigma_{8}$ \+ Pantheon |  $\Lambda$CDM | 1831.71 | 1837.72 | 0 | 1854.14 | 0
$f(R)$ | 1784.84 | 1792.87 | 44.86 | 1814.75 | 39.39
|  |  |  |  |  |

The results of the IC analysis are presented in Table 3. Though we are using
two models ($\Lambda$CDM vs. $f(R)$) only, the comparison is performed from
the results of the statistical analyses of the different datasets (separately
and jointly), as described above. For each case, in Table 3 we report the
values of $\chi^{2}_{\rm{min}}$, AIC (Eq. (27)), BIC (Eq. (28)), and $|\Delta
IC|$ (computed for both criteria using Eq. (29)).

If only $\chi^{2}_{\rm{min}}$ is considered, we see that the proposed $f(R)$
model provides a better fit to the considered data (excluding ($H(z)$ alone)
than $\Lambda$CDM. Then, by looking at the values of AIC and BIC, for which
the number of parameters is considered, the situation is more convoluted,
since the IC is lower for $\Lambda$CDM in some cases and larger in others. In
particular, notice that the largest differences ($|\Delta IC|$) are obtained
when the Pantheon sample is used in the joint statistical analysis (together
with $H(z)$ and/or $f\sigma_{8}$). In these cases, we observe what appears to
be a strong preference of our proposed model over $\Lambda$CDM. Nevertheless,
notice that the results of the other data samples does not provide an
indication in favor of any of the models, but points to the compatibility
between them, and to the fact that both of the models are equally likely to
reproduce the data. Also, we have to consider the fact that the proposed
$f(R)$ model is originates as a perturbation from $\Lambda$CDM, so the results
are not astonishing 111Let us point out that, as marked in [39] and detailed
in [67], this kind of analysis should not be taken as a final word when
comparing different models, but as a complementary tool., and additional tests
might be performed.

##  4 Cosmological dynamics in late-time

Finally, setting the parameters of the model to the BF values obtained from
the joint fit (Table 2, last row), we can take a look at the cosmological
dynamics at late time as described by the $f(R)$ model studied here.

###  4.1 Om($z$) Diagnostic

An interesting tool to study the dynamics of a particular model is the Om
diagnostic proposed in Ref. [5], which relies on the Hubble parameter, $H(z)$.
With this diagnostic, it is also possible to analyse difference between the
proposed model and $\Lambda$CDM. The diagnostic is performed by computing

$\rm{Om}(z)=\frac{E^{2}(z)-1}{(1+z)^{3}-1},$ (30)

where $E^{2}(z)=H(z)/H_{0}$. Looking at the evolution of $\rm{Om}(z)$, one can
obtain information about the nature of DE as predicted by the considered
model: if the model predicts a quintessence behaviour, $\rm{Om}(z)$ would
exhibit a negative slope (decreasing evolution); if, instead, the prediction
favors a phantom DE, $\rm{Om}(z)$ increases with $z$, showing a positive
slope; finally, $\rm{Om}(z)$ remains constant, corresponds to a cosmological
constant DE, _i.e._ , the standard $\Lambda$CDM model.

For the $f(R)$ model studied in this work, we can compute $\rm{Om}(z)$ by
means of the analytical solution, Eq. (15), considering the BF values of the
relevant parameters $(\Omega_{\rm{m}0},b)$, obtained from the
$H(z)$+Pantheon+$f\sigma_{8}$-data joint statistical analysis (last row of
Table 2).

Figure 6: Evolution of Om terms of the redshift, as predicted by the $f(R)$
model proposed in this work, considering the constraints from $H(z)$, Pantheon
and $f\sigma_{8}$ samples.

The resulting evolution of $\rm{Om}(z)$ is shown in Figure 6. Notice how, for
$z<0$ and $z\gtrsim 3.5$, $\rm{Om}(z)$ presents a negligible variation (zero
slope), indicating that the effective DE would behave like a cosmological
constant. For the region in between, and for $z=0$, on particular,
$\rm{Om}(z)$ decreases (negative slope), implying that the effective DE of our
model displays a quintessence-like behaviour, which is consistent with the
evolution of the DE EoS, $w_{\rm{DE}}$ at most of the corresponding $z$
interval (see left-lower panel of Fig. 7).

###  4.2 Cosmological parameters

We now consider interesting cosmological parameters as given by the proposed
$f(R)$ model, which provide insights on the model predictions and evolution,
as well as a suitable way to compare with $\Lambda$CDM. In particular, here we
examine the effective EoS,

$w_{\rm{eff}}=-1+\frac{1}{3}(1+z)\frac{(E^{2}(z))^{\prime}}{E^{2}(z)},$ (31)

the deceleration parameter,

$q=-1+\frac{1}{2}(1+z)\frac{(E^{2}(z))^{\prime}}{E^{2}(z)},$ (32)

the DE EoS,

$w_{\rm{DE}}=-1+\frac{1}{3}(1+z)\frac{(\rho_{\rm{DE}}(z))^{\prime}}{\rho_{\rm{DE}}(z)},$
(33)

and the DE density,

$\Omega_{\rm{DE}}=1-\frac{\Omega_{\rm{m}0}(1+z)^{3}}{E^{2}(z)};$ (34)

in all the above expressions, the prime indicates derivative with respect to
$z$.

The $f(R)$-model predicted evolution of these quantities is shown in Fig. 7
(red dotted lines) in terms of the redshift, $z$, where we also include the
$\Lambda$CDM prediction (blue lines), for comparison.

Figure 7: Comparison of the evolution of some cosmological parameters in terms
of the redshift, as predicted by the $\Lambda$CDM (blue-full line) model and
the $f(R)$ (red-dotted line) model proposed in this work.

In spite of the fact that the statistical analysis showed a preference towards
$b\sim\mathcal{O}(10^{-1})$, the cosmological evolution of $w_{\rm{eff}}$,
$q$, and $\Omega_{\rm{DE}}$ predicted by the $f(R)$ model closely resembles
the prediction of $\Lambda$CDM. The largest deviation appears in the range
$0.5\lesssim z\lesssim 3$, most certainly due to the fact that the
approximated solution implemented in this analysis considers a perturbative
expansion up to a second order; if additional terms (proportional to $b^{n}$,
$n>2$) had been considered, the $f(R)$ model would have resulted to be much
closer to $\Lambda$CDM, and the red dotted lines in Fig. 7 would be almost
indistinguishable from the blue ones. This is expected since it has already
been shown in Refs. [46, 41] that using the exact (numerical) solution for the
Hubble expansion rate, the $f(R)$ model is essentially indistinguishable from
$\Lambda$CDM at the background level.

As observed in the bottom-left panel of Fig. 7, $w_{\rm{DE}}$ shows a
considerable deviation from $w_{\rm{DE}}=-1$ along the depicted range,
specially for $z\lesssim 4$, where an oscillatory behaviour is observed. This
discrepancy (although with a lower amplitude) was already anticipated in Ref.
[42] for a smaller value of the deviation parameter $b$; it is also apparent
that $w_{\rm{DE}}\to-1$ at the early stages of the Universe ($z\gtrsim 4$), as
also observed in [42]. This is another indication of the effect of using the
perturbative expansion up to the second order. In fact, it is reasonable to
think that, as a consequence of $b\sim\mathcal{O}(10^{-1})$ as resulted from
the observational data analysis, additional terms in the expansion of Eq.
(15), proportional to $b^{3}$ and larger powers, might contribute
substantially to the solution, likely mitigating the oscillations of
$w_{\rm{DE}}$. It is also important to notice that this oscillatory evolution
has already been observed by other authors, for instance, in the context of
modified gravity models [68, 69], or considering dynamical dark energy models
[70, 71].

##  5 Conclusions

In this paper we have performed a statistical analysis of a viable, known
$f(R)$ gravity model which includes an exponential function of the scalar
curvature, Eq. (14), with a specific parameter $b$ governing the deviation of
the $f(R)$ model from $\Lambda$CDM. Within this context, we implemented the
analytical approximate solution for the expansion rate, $H(z)$, shown in Eq.
15, from which some observational quantities can be computed, allowing to
investigate the impact of the truncating the perturbative expansion with
respect to $b$.

In addition to $H(z)$ itself, we considered the distance modulus, $\mu(z)$,
and the growth rate multiplied by the amplitude of the matter power spectrum
at $8h^{-1}\rm{Mpc}$, $f\sigma_{8}(z)$. Hence, for the statistical analysis we
used observational data from cosmic chronometers and radial Baryon Acoustic
Oscillations methods (Section 3.1), the SN Ia Pantheon+SH0ES sample (Section
3.2), and the growth sample (Section 3.4. We analysed these data samples
individually to set constraints on the model parameters,
$\left(\sigma_{8},H_{0},\Omega_{m0},b\right)$ and found that, in particular,
the growth sample alone does not provide reasonable bounds on $H_{0}$, and
that, from the three analyses, the value of the deviation parameter $b$ that
best fit each data set is $\sim\mathcal{O}(10^{-1})$, and that $b=0$ is
excluded at $>1\sigma$ (see Fig. 1 and Table 2).

Strengthened constraints on the parameters were obtained by performing joint
analyses. By only combining $H(z)$ and the Pantheon samples, the bounds on
$\left(H_{0},\Omega_{m0},b\right)$ are considerably improved (Fig. 1),
remarkably locating $H_{0}$ in a region well in between the observations made
by Planck, on one side, and SH0ES, on the other. Again, for the deviation
parameter $b$ the allowed region is such that $b=0$ is excluded even further.
Similar results are obtained when combining $H(z)$ with the growth sample, and
Pantheon with the growth sample (Fig. 2). From the former (blue contours and
lines), a null deviation parameter (i.e., $b=0$) is not prohibited, while from
the later, the best fit and the allowed region of $H_{0}$ is pushed to larger
values, implying a better agreement with SH0ES than with Planck.

As expected, the joint fit of the three data samples delivers the strongest
constraints on the considered parameters (gray contours and lines in (Fig. 2),
and last row of Table 2). Up to the second order of perturbative expansion on
the deviation parameter $b$, the proposed $f(R)$ model appropriately
reproduces the data (Figs. 3 \- 5) with

$\displaystyle H_{0}$ $\displaystyle=70.7^{+1.3}_{-1.2}\,\,\,\rm{km/s/Mpc},\,$
$\displaystyle\Omega_{m0}$ $\displaystyle=0.245^{+0.022}_{-0.023},$
$\displaystyle\sigma_{8}$ $\displaystyle=0.844\pm 0.028,\,$ $\displaystyle b$
$\displaystyle=0.648^{+0.125}_{-0.137},$

results that indicate that our model alleviate the tension between Planck and
SH0ES regarding $H_{0}$, on one side, and that the preferred value of $b$
turns out to be larger than initially expected, and certainly $b\neq 0$ at
$\gtrsim 4.5\sigma$. Notwithstanding, this is not entirely stunning, since
this has also been obtained by different authors previously. Furthermore, we
also looked at the predicted evolution of some interesting cosmological
parameters (Section 4), noticing that the effective equation of state, the
deceleration parameter and the DE density exhibit the expected behaviour,
deviating slightly from $\Lambda$CDM. With regards to the DE EoS, although the
deviation is more evident, its oscillatory evolution is not unexpected (it has
been observed by other authors, e.g., [68, 69, 70, 71]), and leads us to the
conclusion that additional terms in the perturbative expansion should diminish
the observable difference with $\Lambda$CDM.

By performing the Om diagnostic, and using the BF values of the constrained
parameters, we have observed that the proposed $f(R)$ model predicts a DE that
behaves like a cosmological constant at early times $(z\gtrsim 3.5)$ and for
the near future $(z<0)$, while at current and late time, the DE exhibits a
quintessence-like evolution, in agreement with the results discussed above
regarding $w_{\rm{DE}}$.

Finally, as an evaluation of the statistical analysis performed in this study,
and a tool to compare different models, we implemented the AIC and BIC
information criteria (Section 3.6), which results are presented in Table 3. We
found that, depending on the analysed data sample, the IC is lower for
$\Lambda$CDM or larger than for the $f(R)$ model proposed here, and that the
largest differences ($|\Delta IC|$) are obtained when the Pantheon sample is
used in a joint statistical analysis, pointing to a preference of our proposed
model (lower IC) over $\Lambda$CDM. However, as in the other cases
$2\lesssim|\Delta IC|\lesssim 6$, the results indicate that the preference
over one model or the other is modest and the two models are essentially
compatible.

Acknowledgements A. O. is supported by _Patrimonio Autónomo–Fondo Nacional de
Financiamiento para la Ciencia, la Tecnología y la Innovación Francisco José
de Caldas (MINCIENCIAS–COLOMBIA)_ Grant No. 110685269447 RC-80740-465-2020,
projects 69723 and 69553. A. O. and M. A. A. express their gratitude to
Ricardo Vega (Universidad del Atlántico) for allowing them to use the computer
resources of his Laboratory for the MCMC analyses.

## Declarations

  * •

Availability of data and materials. All data used in the present analysis have
been released and published by the corresponding research teams, and we
properly include all necessary References to those works, and hence no further
data deposit is needed.

## References

  *   * Riess et al. [1998] Riess, A.G., et al.: Observational evidence from supernovae for an accelerating universe and a cosmological constant. Astron. J. 116, 1009–1038 (1998) https://doi.org/10.1086/300499 arXiv:astro-ph/9805201
  * Perlmutter et al. [1999] Perlmutter, S., et al.: Measurements of $\Omega$ and $\Lambda$ from 42 High Redshift Supernovae. Astrophys. J. 517, 565–586 (1999) https://doi.org/10.1086/307221 arXiv:astro-ph/9812133
  * Peebles and Ratra [2003] Peebles, P.J.E., Ratra, B.: The Cosmological Constant and Dark Energy. Rev. Mod. Phys. 75, 559–606 (2003) https://doi.org/10.1103/RevModPhys.75.559 arXiv:astro-ph/0207347
  * Copeland et al. [2006] Copeland, E.J., Sami, M., Tsujikawa, S.: Dynamics of dark energy. Int. J. Mod. Phys. D 15, 1753–1936 (2006) https://doi.org/10.1142/S021827180600942X arXiv:hep-th/0603057
  * Sahni and Starobinsky [2006] Sahni, V., Starobinsky, A.: Reconstructing Dark Energy. Int. J. Mod. Phys. D 15, 2105–2132 (2006) https://doi.org/10.1142/S0218271806009704 arXiv:astro-ph/0610026
  * Bamba et al. [2012] Bamba, K., Capozziello, S., Nojiri, S., Odintsov, S.D.: Dark energy cosmology: the equivalent description via different theoretical models and cosmography tests. Astrophys. Space Sci. 342, 155–228 (2012) https://doi.org/10.1007/s10509-012-1181-8 arXiv:1205.3421 [gr-qc]
  * Amendola et al. [2007] Amendola, L., Gannouji, R., Polarski, D., Tsujikawa, S.: Conditions for the cosmological viability of f(R) dark energy models. Phys. Rev. D 75, 083504 (2007) https://doi.org/10.1103/PhysRevD.75.083504 arXiv:gr-qc/0612180
  * Hwang and Noh [2001] Hwang, J.-c., Noh, H.: f(R) gravity theory and CMBR constraints. Phys. Lett. B 506, 13–19 (2001) https://doi.org/10.1016/S0370-2693(01)00404-X arXiv:astro-ph/0102423
  * Nojiri and Odintsov [2003] Nojiri, S., Odintsov, S.D.: Modified gravity with negative and positive powers of the curvature: Unification of the inflation and of the cosmic acceleration. Phys. Rev. D 68, 123512 (2003) https://doi.org/10.1103/PhysRevD.68.123512 arXiv:hep-th/0307288
  * Capozziello et al. [2005] Capozziello, S., Cardone, V.F., Troisi, A.: Reconciling dark energy models with f(R) theories. Phys. Rev. D 71, 043503 (2005) https://doi.org/10.1103/PhysRevD.71.043503 arXiv:astro-ph/0501426
  * Cognola et al. [2005] Cognola, G., Elizalde, E., Nojiri, S., Odintsov, S.D., Zerbini, S.: One-loop f(R) gravity in de Sitter universe. JCAP 02, 010 (2005) https://doi.org/10.1088/1475-7516/2005/02/010 arXiv:hep-th/0501096
  * Capozziello et al. [2006] Capozziello, S., Cardone, V.F., Francaviglia, M.: f(R) Theories of gravity in Palatini approach matched with observations. Gen. Rel. Grav. 38, 711–734 (2006) https://doi.org/10.1007/s10714-006-0261-x arXiv:astro-ph/0410135
  * Nojiri and Odintsov [2006] Nojiri, S., Odintsov, S.D.: Modified f(R) gravity consistent with realistic cosmology: From matter dominated epoch to dark energy universe. Phys. Rev. D 74, 086005 (2006) https://doi.org/10.1103/PhysRevD.74.086005 arXiv:hep-th/0608008
  * Song et al. [2007] Song, Y.-S., Hu, W., Sawicki, I.: The Large Scale Structure of f(R) Gravity. Phys. Rev. D 75, 044004 (2007) https://doi.org/10.1103/PhysRevD.75.044004 arXiv:astro-ph/0610532
  * Faulkner et al. [2007] Faulkner, T., Tegmark, M., Bunn, E.F., Mao, Y.: Constraining f(R) Gravity as a Scalar Tensor Theory. Phys. Rev. D 76, 063505 (2007) https://doi.org/10.1103/PhysRevD.76.063505 arXiv:astro-ph/0612569
  * Olmo [2007] Olmo, G.J.: Limit to general relativity in f(R) theories of gravity. Phys. Rev. D 75, 023511 (2007) https://doi.org/10.1103/PhysRevD.75.023511 arXiv:gr-qc/0612047
  * Sawicki and Hu [2007] Sawicki, I., Hu, W.: Stability of Cosmological Solution in f(R) Models of Gravity. Phys. Rev. D 75, 127502 (2007) https://doi.org/10.1103/PhysRevD.75.127502 arXiv:astro-ph/0702278
  * Faraoni [2007] Faraoni, V.: de Sitter space and the equivalence between f(R) and scalar-tensor gravity. Phys. Rev. D 75, 067302 (2007) https://doi.org/10.1103/PhysRevD.75.067302 arXiv:gr-qc/0703044
  * Bean et al. [2007] Bean, R., Bernat, D., Pogosian, L., Silvestri, A., Trodden, M.: Dynamics of Linear Perturbations in f(R) Gravity. Phys. Rev. D 75, 064020 (2007) https://doi.org/10.1103/PhysRevD.75.064020 arXiv:astro-ph/0611321
  * Nojiri and Odintsov [2007] Nojiri, S., Odintsov, S.D.: Unifying inflation with LambdaCDM epoch in modified f(R) gravity consistent with Solar System tests. Phys. Lett. B 657, 238–245 (2007) https://doi.org/10.1016/j.physletb.2007.10.027 arXiv:0707.1941 [hep-th]
  * Capozziello et al. [2007] Capozziello, S., Stabile, A., Troisi, A.: The Newtonian Limit of f(R) gravity. Phys. Rev. D 76, 104019 (2007) https://doi.org/10.1103/PhysRevD.76.104019 arXiv:0708.0723 [gr-qc]
  * Deruelle et al. [2008] Deruelle, N., Sasaki, M., Sendouda, Y.: Junction conditions in f(R) theories of gravity. Prog. Theor. Phys. 119, 237–251 (2008) https://doi.org/10.1143/PTP.119.237 arXiv:0711.1150 [gr-qc]
  * Appleby and Battye [2008] Appleby, S.A., Battye, R.A.: Aspects of cosmological expansion in F(R) gravity models. JCAP 05, 019 (2008) https://doi.org/10.1088/1475-7516/2008/05/019 arXiv:0803.1081 [astro-ph]
  * Carloni et al. [2008] Carloni, S., Dunsby, P.K.S., Troisi, A.: The Evolution of density perturbations in f(R) gravity. Phys. Rev. D 77, 024024 (2008) https://doi.org/10.1103/PhysRevD.77.024024 arXiv:0707.0106 [gr-qc]
  * Capozziello et al. [2008] Capozziello, S., Cardone, V.F., Salzano, V.: Cosmography of f(R) gravity. Phys. Rev. D 78, 063504 (2008) https://doi.org/10.1103/PhysRevD.78.063504 arXiv:0802.1583 [astro-ph]
  * Dunsby et al. [2010] Dunsby, P.K.S., Elizalde, E., Goswami, R., Odintsov, S., Gomez, D.S.: On the LCDM Universe in f(R) gravity. Phys. Rev. D 82, 023519 (2010) https://doi.org/10.1103/PhysRevD.82.023519 arXiv:1005.2205 [gr-qc]
  * Capozziello and De Laurentis [2012] Capozziello, S., De Laurentis, M.: The dark matter problem from f(R) gravity viewpoint. Annalen Phys. 524, 545–578 (2012) https://doi.org/10.1002/andp.201200109
  * Odintsov and Oikonomou [2019a] Odintsov, S.D., Oikonomou, V.K.: $f(R)$ Gravity Inflation with String-Corrected Axion Dark Matter. Phys. Rev. D 99(6), 064049 (2019) https://doi.org/%****␣AOliveros-MAcero_fR_Constraints_epjc.bbl␣Line␣500␣****10.1103/PhysRevD.99.064049 arXiv:1901.05363 [gr-qc]
  * Odintsov and Oikonomou [2019b] Odintsov, S.D., Oikonomou, V.K.: Unification of Inflation with Dark Energy in $f(R)$ Gravity and Axion Dark Matter. Phys. Rev. D 99(10), 104070 (2019) https://doi.org/10.1103/PhysRevD.99.104070 arXiv:1905.03496 [gr-qc]
  * Odintsov and Oikonomou [2020] Odintsov, S.D., Oikonomou, V.K.: Geometric Inflation and Dark Energy with Axion $F(R)$ Gravity. Phys. Rev. D 101(4), 044009 (2020) https://doi.org/10.1103/PhysRevD.101.044009 arXiv:2001.06830 [gr-qc]
  * Oikonomou [2021a] Oikonomou, V.K.: Rescaled Einstein-Hilbert Gravity from $f(R)$ Gravity: Inflation, Dark Energy and the Swampland Criteria. Phys. Rev. D 103(12), 124028 (2021) https://doi.org/10.1103/PhysRevD.103.124028 arXiv:2012.01312 [gr-qc]
  * Oikonomou [2021b] Oikonomou, V.K.: Unifying inflation with early and late dark energy epochs in axion $F(R)$ gravity. Phys. Rev. D 103(4), 044036 (2021) https://doi.org/10.1103/PhysRevD.103.044036 arXiv:2012.00586 [astro-ph.CO]
  * Nojiri and Odintsov [2011] Nojiri, S., Odintsov, S.D.: Unified cosmic history in modified gravity: from F(R) theory to Lorentz non-invariant models. Phys. Rept. 505, 59–144 (2011) https://doi.org/10.1016/j.physrep.2011.04.001 arXiv:1011.0544 [gr-qc]
  * Clifton et al. [2012] Clifton, T., Ferreira, P.G., Padilla, A., Skordis, C.: Modified Gravity and Cosmology. Phys. Rept. 513, 1–189 (2012) https://doi.org/10.1016/j.physrep.2012.01.001 arXiv:1106.2476 [astro-ph.CO]
  * Nojiri et al. [2017] Nojiri, S., Odintsov, S.D., Oikonomou, V.K.: Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution. Phys. Rept. 692, 1–104 (2017) https://doi.org/10.1016/j.physrep.2017.06.001 arXiv:1705.11098 [gr-qc]
  * Hu and Sawicki [2007] Hu, W., Sawicki, I.: Models of f(R) Cosmic Acceleration that Evade Solar-System Tests. Phys. Rev. D 76, 064004 (2007) https://doi.org/10.1103/PhysRevD.76.064004 arXiv:0705.1158 [astro-ph]
  * Starobinsky [2007] Starobinsky, A.A.: Disappearing cosmological constant in f(R) gravity. JETP Lett. 86, 157–163 (2007) https://doi.org/10.1134/S0021364007150027 arXiv:0706.2041 [astro-ph]
  * Basilakos et al. [2013] Basilakos, S., Nesseris, S., Perivolaropoulos, L.: Observational constraints on viable f(R) parametrizations with geometrical and dynamical probes. Phys. Rev. D 87(12), 123529 (2013) https://doi.org/10.1103/PhysRevD.87.123529 arXiv:1302.6051 [astro-ph.CO]
  * Pérez-Romero and Nesseris [2018] Pérez-Romero, J., Nesseris, S.: Cosmological constraints and comparison of viable $f(R)$ models. Phys. Rev. D 97(2), 023525 (2018) https://doi.org/10.1103/PhysRevD.97.023525 arXiv:1710.05634 [astro-ph.CO]
  * Kumar et al. [2023] Kumar, S., Nunes, R.C., Pan, S., Yadav, P.: New late-time constraints on $f(R)$ gravity. Phys. Dark Univ. 42, 101281 (2023) https://doi.org/10.1016/j.dark.2023.101281 arXiv:2301.07897 [astro-ph.CO]
  * Oliveros and Acero [2023] Oliveros, A., Acero, M.A.: Late-time cosmology in a model of modified gravity with an exponential function of the curvature. Phys. Dark Univ. 40, 101207 (2023) https://doi.org/10.1016/j.dark.2023.101207 arXiv:2302.07022 [gr-qc]
  * Oliveros [2023] Oliveros, A.: A viable $f(R)$ gravity model without oscillations in the effective dark energy. Int. J. Mod. Phys. D 32(12), 2350086 (2023) https://doi.org/10.1142/S0218271823500864 arXiv:2307.11896 [astro-ph.CO]
  * Aghanim et al. [2020] Aghanim, N., et al.: Planck 2018 results. VI. Cosmological parameters. Astron. Astrophys. 641, 6 (2020) https://doi.org/10.1051/0004-6361/201833910 arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)]
  * Riess et al. [2022a] Riess, A.G., et al.: A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s-1 Mpc-1 Uncertainty from the Hubble Space Telescope and the SH0ES Team. Astrophys. J. Lett. 934(1), 7 (2022) https://doi.org/10.3847/2041-8213/ac5c5b arXiv:2112.04510 [astro-ph.CO]
  * Riess et al. [2022b] Riess, A.G., Breuval, L., Yuan, W., Casertano, S., Macri, L.M., Bowers, J.B., Scolnic, D., Cantat-Gaudin, T., Anderson, R.I., Reyes, M.C.: Cluster Cepheids with High Precision Gaia Parallaxes, Low Zero-point Uncertainties, and Hubble Space Telescope Photometry. Astrophys. J. 938(1), 36 (2022) https://doi.org/10.3847/1538-4357/ac8f24 arXiv:2208.01045 [astro-ph.CO]
  * Granda [2020] Granda, L.N.: Modified gravity with an exponential function of curvature. Eur. Phys. J. C 80(6), 539 (2020) https://doi.org/10.1140/epjc/s10052-020-8114-4 arXiv:2003.09006 [gr-qc]
  * Linder [2009] Linder, E.V.: Exponential Gravity. Phys. Rev. D 80, 123528 (2009) https://doi.org/%****␣AOliveros-MAcero_fR_Constraints_epjc.bbl␣Line␣825␣****10.1103/PhysRevD.80.123528 arXiv:0905.2962 [astro-ph.CO]
  * Cognola et al. [2008] Cognola, G., Elizalde, E., Nojiri, S., Odintsov, S.D., Sebastiani, L., Zerbini, S.: A Class of viable modified f(R) gravities describing inflation and the onset of accelerated expansion. Phys. Rev. D 77, 046009 (2008) https://doi.org/10.1103/PhysRevD.77.046009 arXiv:0712.4017 [hep-th]
  * Odintsov et al. [2017] Odintsov, S.D., Sáez-Chillón Gómez, D., Sharov, G.S.: Is exponential gravity a viable description for the whole cosmological history? Eur. Phys. J. C 77(12), 862 (2017) https://doi.org/10.1140/epjc/s10052-017-5419-z arXiv:1709.06800 [gr-qc]
  * Odintsov et al. [2019] Odintsov, S.D., Saez-Chillon Gomez, D., Sharov, G.S.: Testing logarithmic corrections on $R^{2}$-exponential gravity by observational data. Phys. Rev. D 99(2), 024003 (2019) https://doi.org/10.1103/PhysRevD.99.024003 arXiv:1807.02163 [gr-qc]
  * Odintsov et al. [2023] Odintsov, S.D., Sáez-Chillón Gómez, D., Sharov, G.S.: Exponential F(R) gravity with axion dark matter. Phys. Dark Univ. 42, 101369 (2023) https://doi.org/10.1016/j.dark.2023.101369 arXiv:2310.20302 [gr-qc]
  * Foreman-Mackey et al. [2013] Foreman-Mackey, D., Hogg, D.W., Lang, D., Goodman, J.: emcee: The MCMC Hammer. Publ. Astron. Soc. Pac. 125, 306–312 (2013) https://doi.org/10.1086/670067 arXiv:1202.3665 [astro-ph.IM]
  * Cao et al. [2021] Cao, S., Zhang, T.-J., Wang, X., Zhang, T.: Cosmological Constraints on the Coupling Model from Observational Hubble Parameter and Baryon Acoustic Oscillation Measurements. Universe 7(3), 57 (2021) https://doi.org/10.3390/universe7030057 arXiv:2103.03670 [astro-ph.CO]
  * Adame et al. [2024] Adame, A.G., et al.: DESI 2024 VI: Cosmological Constraints from the Measurements of Baryon Acoustic Oscillations (2024) arXiv:2404.03002 [astro-ph.CO]
  * Brout et al. [2022] Brout, D., et al.: The Pantheon+ Analysis: Cosmological Constraints. Astrophys. J. 938(2), 110 (2022) https://doi.org/10.3847/1538-4357/ac8e04 arXiv:2202.04077 [astro-ph.CO]
  * PantheonPlusSH0ES [2022] PantheonPlusSH0ES: PantheonPlusSH0ES/datarelease (2022). https://github.com/PantheonPlusSH0ES/DataRelease
  * Cardone et al. [2012] Cardone, V.F., Camera, S., Diaferio, A.: An updated analysis of two classes of f(R) theories of gravity. JCAP 02, 030 (2012) https://doi.org/10.1088/1475-7516/2012/02/030 arXiv:1201.3272 [astro-ph.CO]
  * Saridakis et al. [2023] Saridakis, E.N., Yang, W., Pan, S., Anagnostopoulos, F.K., Basilakos, S.: Observational constraints on soft dark energy and soft dark matter: Challenging $\Lambda$CDM cosmology. Nucl. Phys. B 986, 116042 (2023) https://doi.org/10.1016/j.nuclphysb.2022.116042 arXiv:2112.08330 [astro-ph.CO]
  * Alestas et al. [2022] Alestas, G., Kazantzidis, L., Nesseris, S.: Machine learning constraints on deviations from general relativity from the large scale structure of the Universe. Phys. Rev. D 106(10), 103519 (2022) https://doi.org/10.1103/PhysRevD.106.103519 arXiv:2209.12799 [astro-ph.CO]
  * Mhamdi et al. [2024] Mhamdi, D., Bouali, A., Dahmani, S., Errahmani, A., Ouali, T.: Cosmological constraints on $f(Q)$ gravity with redshift space distortion data. Eur. Phys. J. C 84(3), 310 (2024) https://doi.org/10.1140/epjc/s10052-024-12549-4
  * Tsujikawa [2007] Tsujikawa, S.: Matter density perturbations and effective gravitational constant in modified gravity models of dark energy. Phys. Rev. D 76, 023514 (2007) https://doi.org/10.1103/PhysRevD.76.023514 arXiv:0705.1032 [astro-ph]
  * Tsujikawa et al. [2009] Tsujikawa, S., Gannouji, R., Moraes, B., Polarski, D.: The dispersion of growth of matter perturbations in f(R) gravity. Phys. Rev. D 80, 084044 (2009) https://doi.org/10.1103/PhysRevD.80.084044 arXiv:0908.2669 [astro-ph.CO]
  * Hou et al. [2024] Hou, J., Moradinezhad Dizgah, A., Hahn, C., Eickenberg, M., Ho, S., Lemos, P., Massara, E., Modi, C., Parker, L., Blancard, B.R.-S.: Cosmological constraints from the redshift-space galaxy skew spectra. Phys. Rev. D 109(10), 103528 (2024) https://doi.org/10.1103/PhysRevD.109.103528 arXiv:2401.15074 [astro-ph.CO]
  * Sultana et al. [2022] Sultana, J., Yennapureddy, M.K., Melia, F., Kazanas, D.: Constraining f(R) models with cosmic chronometers and the H ii galaxy Hubble diagram. Mon. Not. Roy. Astron. Soc. 514(4), 5827–5839 (2022) https://doi.org/%****␣AOliveros-MAcero_fR_Constraints_epjc.bbl␣Line␣1125␣****10.1093/mnras/stac1713 arXiv:2206.10761 [astro-ph.CO]
  * Mandal et al. [2023] Mandal, S., Pradhan, S., Sahoo, P.K., Harko, T.: Cosmological observational constraints on the power law f(Q) type modified gravity theory. Eur. Phys. J. C 83(12), 1141 (2023) https://doi.org/10.1140/epjc/s10052-023-12339-4 arXiv:2310.00030 [gr-qc]
  * Nesseris et al. [2017] Nesseris, S., Pantazis, G., Perivolaropoulos, L.: Tension and constraints on modified gravity parametrizations of $G_{\textrm{eff}}(z)$ from growth rate and Planck data. Phys. Rev. D 96(2), 023542 (2017) https://doi.org/10.1103/PhysRevD.96.023542 arXiv:1703.10538 [astro-ph.CO]
  * Nesseris and Garcia-Bellido [2013] Nesseris, S., Garcia-Bellido, J.: Is the Jeffreys’ scale a reliable tool for Bayesian model comparison in cosmology? JCAP 08, 036 (2013) https://doi.org/10.1088/1475-7516/2013/08/036 arXiv:1210.7652 [astro-ph.CO]
  * Granda [2021] Granda, L.N.: Modified gravity with disappearing cosmological constant. JHEP 12, 205 (2021) https://doi.org/10.1007/JHEP12(2021)205 arXiv:2007.13956 [gr-qc]
  * Granda [2022] Granda, L.N.: Two-parametric f(R) dark energy model and some observational constraints. Mod. Phys. Lett. A 37(15), 2250089 (2022) https://doi.org/10.1142/S0217732322500894
  * Zhao et al. [2017] Zhao, G.-B., et al.: Dynamical dark energy in light of the latest observations. Nature Astron. 1(9), 627–632 (2017) https://doi.org/10.1038/s41550-017-0216-z arXiv:1701.08165 [astro-ph.CO]
  * Escamilla et al. [2024] Escamilla, L.A., Pan, S., Di Valentino, E., Paliathanasis, A., Vázquez, J.A., Yang, W.: Oscillations in the Dark? (2024) arXiv:2404.00181 [astro-ph.CO]