# Investigating the accelerated expansion of the Universe through updated
constraints on viable $f(R)$ models within the metric formalism

Kumar Ravi 1<EMAIL_ADDRESS>Anirban Chatterjee 2<EMAIL_ADDRESS>Biswajit Jana 1<EMAIL_ADDRESS>Abhijit Bandyopadhyay 1
<EMAIL_ADDRESS>1Department of Physics, Ramakrishna Mission Vivekananda
Educational and Research Institute, Belur Math 711202, West Bengal, India
2Department of Physics, Indian Institute of Technology Kanpur, Kanpur 208016,
UP, India

###### Abstract

Modified theories of gravity encompass a class of $f(R)$-models that seek to
elucidate the observed late time accelerated expansion of the universe. In
this study, we examine a set of viable $f(R)$ models (Hu-Sawicki: two cases,
Satrobinsky, Tsujikawa, exponential and arcTanh models) in metric formalism,
using recent cosmological data sets: type Ia supernovae data, cosmic
chronometer observations, baryonic acoustic oscillations data, data from Hii
starburst galaxies, and local measurements of the Hubble parameter ($H_{0}$).
We re-parameterize the $f(R)$ models using a distortion/deviation parameter
($b$) which is a measure of their deviation from the standard $\Lambda$CDM
model. The model parameters are constrained using a Bayesian analysis with the
Monte Carlo Markov Chain method. We employ statistical tools such as the
Akaike Information Criterion, Bayesian Information Criterion, and reduced chi-
square statistics to conduct a comparative investigation of these models. We
determine the transition redshift, the evolution of total equation-of-state
(EoS) parameter, and the EoS for the component responsible for current
accelerated expansion to characterize the expansion’s evolution. Taking into
account the “Hubble tension,” we perform the study with and without a Gaussian
prior for $H_{0}$ from local measurements. Our findings are as follows: (i) in
many cases the $f(R)$ models are strongly favored over the standard
$\Lambda$CDM model, (ii) the deviation parameter ($b$) significantly deviates
from zero in several cases, (iii) the inclusion of local $H_{0}$ not only
increases the fitted value of $H_{0}$ (as expected) but also affects the gap
between predictions of $f(R)$ models and the $\Lambda$CDM model, and (iv) the
relevant quantities characterizing the (accelerated) expansion of the universe
obtained in our models are consistent with those obtained in a model-
independent way by others. Our investigation and results present a compelling
case for pursuing further research on $f(R)$ models with future observations
to come.

##  I Introduction

In the last two decades, cosmological investigations have revealed an
accelerated expansion of the current universe and have also identified a
transition from a decelerating phase to the current accelerated phase
occurring during the late time phase of cosmic evolution. The first empirical
aspect of discovering this phenomenon came from the interpretation of
luminosity distance and redshift measurements of type Ia supernovae (SNe Ia)
events Perlmutter et al. (1998, 1999); Schmidt et al. (1998); Riess et al.
(1998). Furthermore, observation of baryon acoustic oscillations Eisenstein et
al. (2005); Percival et al. (2007a), analysis of cosmic microwave background
radiation Spergel et al. (2007); Hinshaw et al. (2009); Ade et al. (2014,
2016), and examination of the power spectrum of matter distributions in the
universe Cole et al. (2005); Percival et al. (2007b) have substantiated
evidence of this late-time cosmic acceleration. A general label describing the
origin of the observed late time cosmic acceleration is “Dark energy” which
refers to a theoretical unclustered form of energy exerting a negative
pressure to counteract the gravitational attraction and thereby causing the
cosmic acceleration.

Despite extensive research over many years, the true nature and origin of dark
energy remains an enigma. Various theoretical perspectives have emerged, each
attempting to construct models that can explain the observed cosmic
acceleration. The introduction of the $\Lambda$ term into Einstein’s equation,
known as the “$\Lambda$ Cold Dark Matter ($\Lambda$CDM) model”, is the
simplest model that can explain the present accelerated expansion of the
universe. However, this model encounters several challenges, such as the
cosmic coincidence issue Zlatev et al. (1999) and the fine-tuning problem
Martin (2012) when considered in the context of particle physics. This
motivates development and exploration of alternative dark energy models from a
range of perspectives. An important class comprises field theoretic models
that involve the incorporation of a scalar field within the energy-momentum
tensor of the Einstein equation. The scalar field plays the role in generating
the necessary negative pressure to drive cosmic acceleration, either through
the slowly varying potentials (quintessence models Tsujikawa (2013)) or by
means of their kinetic energy (k-essence models Armendariz-Picon et al.
(2001)).

Within the context of this study, another crucial category of models under
consideration is modified gravity models which attempt to explain the
acceleration through the geometry itself without modifying the energy-momentum
tensor of the Einstein equation. These models primarily involve modifications
to the geometric component of Einstein’s equation, which may result from
higher-order corrections to the Einstein-Hilbert action. By introducing
suitable modifications, it becomes feasible to induce cosmic acceleration. The
most straightforward types of modifications involve the extension of the Ricci
scalar $R$ to an arbitrary function $f(R)$. The meticulous choice with
appropriate justification of this particular arbitrary function play a crucial
role in all modified gravity models. Theoretical considerations, like the need
for a ghost-free theory with stable perturbations and the presence of Noether
symmetries, impose initial constraints on the forms of the arbitrary function
$f(R)$. However, ability to reproduce the observed features of cosmic
evolution, the behavior of local (solar) systems etc. also serve as the
primary tool to further constrain the $f(R)$-models.

In this study, our objective is to provide updated constraints on viable
$f(R)$ gravity models using the latest cosmological data, including supernovae
type Ia(SNIa) from Pantheonplus compilation, cosmic chronometer(CC)
observations, baryonic acoustic oscillations(BAO) data, Hii starburst
galaxies(HiiG) data, and data from local measurements of the Hubble parameter
($H_{0}$). Specifically, we focus on six viable $f(R)$ models: the Hu-Sawicki
model(two cases), the Starobinsky model, the exponential model, the Tsujikawa
model, and the arcTanh model. By considering these models, we obtained an
update of best-fit values and uncertainties at different confidence levels of
the associated parameters of the models. We have chosen these particular
$f(R)$ gravity models because most of the modified gravity models either fail
to explain the matter-dominated era or have already been ruled out by
observational data sets. The selected models, mentioned above, represent a few
remaining options to study the impact of modified gravity theory within the
framework of metric or Palatini formalism. For our investigation, we adopt the
metric formalism techniques in the context of a homogeneous and isotropic
universe. Rather than working with original forms of these $f(R)$ models(which
for some models give a false impression that they are non-reducible to the
$\Lambda$CDM model), we chose to re-parameterize these models in terms of what
is called “deviation/distortion parameter($b$)”. In this form, with
$b\rightarrow 0$, a $f(R)$ model tends to the $\Lambda$CDM model.

The current state of research in this specific context, which involves
constraining various $f(R)$ models within the metric formalism using
cosmologically relevant data, is experiencing a high level of activity. Here
we refer to several such relevant previous works. In Basilakos et al. (2013),
the authors obtained constraints on the Hu-Sawicki and the Starobinsky models
by utilizing approximate analytical solutions derived from Taylor series
expansions for the Hubble parameter. The constraints were derived using SNIa
data from Union 2.1 compilation, BAO data, cosmic microwave background (CMB)
shift parameter data and growth rate data. In Nunes et al. (2017), the
constraints on Hu-Sawicki, Starobinsky, exponential, and Tsujikawa models were
obtained by utilizing CC data, local measurements of $H_{0}$, SNIa data from
the Joint-Lightcurve-Analysis (JLA) compilation, and BAO data. Furthermore,
Odintsov et al. (2017) constrained the exponential model using SNIa (Union
2.1), CC, BAO, and CMB data, while also discussing the viability of this model
in describing the entire cosmological history. Exploration of various $f(R)$
models, including the Hu-Sawicki and the exponential models, with data sets
such as SNIa (JLA), CC, BAO, CMB, local $H_{0}$, and growth rate, was
conducted in Pérez-Romero and Nesseris (2018). Other significant works that
investigated Hu-Sawicki and/or exponential models, utilizing various data sets
including gravitational lensing data, are D’Agostino and Nunes (2019, 2020);
Odintsov et al. (2021). In Farrugia et al. (2021), constraints on Hu-Sawicki,
Starobinsky, exponential, and Tsujikawa models were examined in both flat and
non-flat spacetimes, utilizing various cosmological data sets. Recently,
Leizerovich et al. (2022) have advocated for the use of quasar X-ray and UV
fluxes data to investigate $f(R)$ models and have constrained Hu-Sawicki and
exponential models accordingly. Moreover, the Taylor series approach from
Basilakos et al. (2013) was extended to include the arcTanh model in Sultana
et al. (2022). Using Gaussian process reconstruction of the Hubble diagram
with CC and HiiG data, the parameters of these three $f(R)$ models were
determined.

This paper is organised as follows. In Section II we derive the modified
Friedmann equations and other relevant equations from the action for $f(R)$
gravity. Section III covers discussions of the conditions that any viable
$f(R)$ models must satisfy, along with brief introductions of the specific
$f(R)$ models investigated in this work. In Section IV, we introduce the
cosmological data sets and the corresponding equations that establish
connection between theory and data. Additionally, we discuss the statistical
procedures employed to obtain constraints on the model parameters. The
obtained constraints on the model parameters for all the models are presented
in Section V. The performance of the different models is assessed using
statistical tools in Section VI. Section VII is dedicated to deriving the
relevant quantities that characterize expansion history of the universe based
on the model constraints. In Section VIII we provide the concluding remarks
for this work. Due to need for careful considerations in solving the modified
Friedmann equations, such as avoiding numerical instabilities and minimizing
computation time, Appendix A is included to discuss the method used for
numerically solving the Friedmann equations. The data for baryonic acoustic
oscillations and cosmic chronometer, collected from different sources, are
tabulated in Appendix in the Tables 3 and 4, respectively.

Unless otherwise mentioned we have set $c=1$ (where, $c$ denotes speed of
light in vacuum), and value of the Hubble parameter is expressed in the unit
km s-1Mpc-1.

##  II $f(R)$ cosmology in Metric formalism

The $f(R)$ theories of gravity involve generalisation of Lagrangian by making
it an arbitrary function of the Ricci scalar $R$ where $f(R)=R$ corresponds to
the standard Einstein theory of gravity. Algebraic expressions for $f(R)$,
different from $f(R)=R$, define different $f(R)$ \- models. The generalized
Einstein-Hilbert action for $f(R)$ theories is

$S=\frac{1}{2\kappa}\int\,d^{4}x\,\sqrt{-g}\,f(R)+S_{\rm m}+S_{\rm r}\,,$ (1)

where $\kappa=8\pi{G}$, $g$ is the determinant of the metric tensor, $S_{\rm
m}$ and $S_{\rm r}$ are the actions for matter fields and radiation fields,
respectively. The variation of this action with respect to the metric gives
the corresponding field equation in metric formalism as

$FR_{{\mu\nu}}-\frac{1}{2}f\textnormal{{g}}_{{\mu\nu}}-\left(\nabla_{\mu}\nabla_{\nu}-\textnormal{{g}}_{{\mu\nu}}\Box\right)F={\kappa}T_{{\mu\nu}}\,,$
(2)

where $F={\partial}f/{\partial}R$,
$\Box\equiv\textnormal{{g}}^{{\mu\nu}}\nabla_{\mu}\nabla_{\nu}$ is the
covariant D’Alambertian and $T_{{\mu\nu}}$ is the energy-momentum tensor of
matter and radiation. This field equation can also be expressed as

$\displaystyle
FG_{{\mu\nu}}-F^{\prime}\nabla_{\mu}\nabla_{\nu}R-F^{\prime\prime}(\nabla_{\mu}R)(\nabla_{\nu}R)$
$\displaystyle+\textnormal{{g}}_{{\mu\nu}}\left[\frac{1}{2}(RF-f)+F^{\prime}\Box{R}+F^{\prime\prime}(\nabla{R})^{2}\right]={\kappa}T_{{\mu\nu}}\,,$
(3)

where $G_{{\mu\nu}}\equiv
R_{{\mu\nu}}-\frac{1}{2}R\textnormal{{g}}_{{\mu\nu}}$ is Einstein tensor and
prime(′) denotes derivative with respect to $R$. Taking trace of both sides of
Eq. 2 we obtain

$\Box{R}=\frac{1}{3F^{\prime}}\left[\kappa{T}-3F^{\prime\prime}(\nabla{R})^{2}+2f-RF\right]\,,$
(4)

using which in Eq. 3 we rewrite the field equation as

$\displaystyle G_{{\mu\nu}}$ $\displaystyle=$
$\displaystyle\frac{1}{F}\bigg{[}F^{\prime\,}\nabla_{\mu}\nabla_{\nu}R+F^{\prime\prime}(\nabla_{\mu}R)(\nabla_{\nu}R)$
(5)
$\displaystyle-\frac{\textnormal{{g}}_{{\mu\nu}}}{6}\left(RF+f+2\kappa{T}\right)+{\kappa}T_{{\mu\nu}}\bigg{]}\,,$

where $T\equiv T^{\mu}_{\mu}$ is the trace of energy-momentum tensor.

In this work we consider the universe to be isotropic and homogeneous at large
scale and described by spatially flat Friedmann-Lemaître-Robertson-Walker
(FLRW) metric

$ds^{2}=-dt^{2}+a(t)^{2}\left[dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta\,d\phi^{2})\right]\,,$
(6)

where $a$ is the FLRW scale factor. The Hubble parameter($H$) is defined as
$H\equiv\dot{a}/a$ where we use a generic notation $\dot{X}$ to denote
derivative of any quantity $X$ with respect to cosmological time ($t$). In
this formalism, for flat FRLW geometry, the Ricci scalar relates to the Hubble
parameter by the relation

$R=6\left(2H^{2}+\dot{H}\right)\,.$ (7)

We consider the content of universe to be a perfect fluid comprising of two
components: radiation and matter in the form of pressureless dust (non-
relativistic). For the epochs($0\leq{z}<10^{4}$) when any interaction between
matter and radiation could be ignored, the energy-momentum tensor for the
fluid can be written as $T^{\mu}_{\nu}=\text{diag}(-\rho,p,p,p)$ with
$\rho=\rho_{\rm m}+\rho_{\rm r}$ and $p=p_{\rm r}$ (with radiation pressure
$p_{\rm r}=\rho_{\rm r}/3$), where the subscripts ‘m’ and ‘r’ stand for matter
and radiation, respectively. Each of the non-interacting components separately
follows the continuity equations

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

The solution to these conservation equations are $\rho_{\rm m}=\rho_{\rm
m0}/a^{3}$, $\rho_{\rm r}=\rho_{\rm r0}/a^{4}$ where subscript ‘0’ denotes
values at present epoch.

In the context of the FLRW universe filled with an ideal perfect fluid
characterised by $\rho$ and $p$, Eq. 4 reduces to

$F^{\prime}\ddot{R}+F^{\prime\prime}\dot{R}^{2}=-\left[\frac{\kappa(3p-\rho)}{3}+\frac{2f-RF}{3}\right]-3HF^{\prime}\dot{R}\,.$
(9)

using which the ‘00’ and ‘$ii$’ components of the field Eq. 5 take following
respective forms:

$\displaystyle-3H^{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{F}\bigg{[}{\kappa}\rho+\frac{RF-f}{2}-3H{F^{\prime}}\dot{R}\bigg{]}\,,$
(10) $\displaystyle-2\dot{H}$ $\displaystyle=$
$\displaystyle\frac{1}{F}\bigg{[}\kappa(\rho+p)+F^{\prime\prime}\dot{R}^{2}+(\ddot{R}-H\dot{R})F^{\prime}\bigg{]}\,.$
(11)

The Eqs. 10 and 11 are modified form of the Friedmann equations for
$f(R)$-models. The temporal profile of the Hubble parameter $H$ is commonly
expressed in the form $H(z)$, $z$ being the redshift related to the FLRW scale
factor by $1+z=1/a$ (where, $a$ is normalised to unity at the present epoch).
We require this profile of $H(z)$ from the cosmological data sets for
obtaining observational constraints on $f(R)$ \- models. For this purpose one
may solve either the system of: (i) Eqs. 9, 10 and 11 or, (ii) Eqs. 7, 10 and
11. We take the path (ii) to solve the system. Eq. 10 serves as a constraint
equation which fixes the initial conditions and must be satisfied at every
integration step during the process of finding solutions. Finding analytical
solutions is almost impossible, and therefore numerical methods are employed.
However, solving this system of ordinary differential equations (ODEs) using a
naive approach often leads to numerical instability. Additional details on how
to solve this system of ODEs are discussed in Appendix A.

If we compare the modified Friedmann Eqs. 10 and 11 to the usual Friedmann
equations with a dark energy component characterised by energy density
$\rho_{\rm DE}$ and pressure $p_{\rm DE}$, i.e. with the equations
$3H^{2}=\kappa\left(\rho_{\rm m}+\rho_{\rm r}+\rho_{\rm DE}\right)$ and
$-2\dot{H}=\kappa\left(\rho_{\rm m}+\rho_{\rm r}+\rho_{\rm DE}+p_{\rm
m}+p_{\rm r}+p_{\rm DE}\right)$, we can deduce the “effective (geometric) dark
energy” with density and pressure corresponding to $f(R)$-theory as

$\rho_{\rm
DE}=\frac{1}{\kappa}\left[\frac{RF-f}{2}-3HF^{\prime}\dot{R}+3(1-F)H^{2}\right]\,,$
(12)

and,

$\begin{split}p_{\rm
DE}&=\frac{1}{\kappa}\bigg{[}\frac{f-RF}{2}+F^{\prime}\ddot{R}+2HF^{\prime}\dot{R}+F^{\prime\prime}\dot{R}^{2}\\\
&\hskip 28.45274pt-(1-F)(2\dot{H}+3H^{2})\bigg{]}\,,\end{split}$ (13)

with the equation-of-state parameter for this effective dark energy defined as

$w_{\rm DE}=\frac{p_{\rm DE}}{\rho_{DE}}\,.$ (14)

Using Eqs. 10 and 11, we may recast Eq. 14 into a more computationally
advantageous form (which we use later) as

$w_{\rm DE}=\frac{w_{\rm tot}-{\kappa}p_{\rm
r}/(3H^{2})}{1-\kappa\left(\rho_{\rm m}+\rho_{\rm r}\right)/(3H^{2})}\,,$ (15)

where

$w_{\rm tot}=-1+\frac{2(1+z)}{3H}\frac{dH}{dz}\,,$ (16)

and we have taken $p_{\rm m}=0$ for the pressureless matter.

The relevant quantities obtained from observations, depending on cosmological
models or even through cosmography (a model-independent kinematical approach),
indicate that the Universe has recently undergone a transition from a phase of
decelerated expansion to accelerated expansion. These quantities include
$w_{\rm tot}|{z=0}\sim-0.7$, $w_{\rm DE}|_{z=0}\sim-1$, and a transition
redshift $(z{\rm t})\sim 0.5-1$. The transition redshift signifies the
redshift at which the transition from decelerated to accelerated expansion
occurred and is determined by the zero-crossing of the deceleration parameter
($q(z)$) given by

$q(z)\equiv-\frac{\ddot{a}}{aH^{2}}=-1+\frac{(1+z)}{H}\frac{dH}{dz}\,.$ (17)

In the $\Lambda$CDM model, the value of $w_{\rm DE}$ is fixed at $-1$ for all
redshifts due to the presence of the cosmological constant ($\Lambda$) term.
However, in $f(R)$ models, the source of $w_{\rm DE}|_{z=0}\sim-1$ arises from
the underlying geometry itself. Unlike the $\Lambda$CDM model, there is no
need to invoke the existence of “dark energy” in $f(R)$ theories of gravity.

##  III The Specific $f(R)$ Models And Their Viability Conditions

In this section we provide a brief introduction to the models examined in this
study, both in their originally proposed forms and their subsequent
transformations into more generalized representations. These transformations
aim to highlight how these models can be more readily reduced to the
$\Lambda$CDM model under appropriate conditions. However, before going into
the details of each model, it is essential to discuss the viability
conditions.

In the metric formalism (unlike Palatini formalism), any viable $f(R)$
cosmological model must satisfy the following set of stringent theoretical
conditions (see De Felice and Tsujikawa (2010); Amendola and Tsujikawa (2015)
for detailed discussions):

$\displaystyle F>0\,,\quad\mathrm{for}\,\,\,R\geq{R}_{0}>0\,,$ (20)
$\displaystyle F^{\prime}>0\,,\quad\mathrm{for}\,\,\,R\geq{R}_{0}>0\,,$
$\displaystyle f(R)\approx R-2\Lambda\,,\quad\mathrm{for}\,\,\,R\gg{R}_{0}\,,$
and, $\displaystyle
0<\frac{RF^{\prime}}{F}(r)<1\quad\mathrm{at}\,\,\,r=-\frac{RF}{f}=-2\,,$ (21)

where $R_{0}$ denotes the Ricci scalar at the present epoch.

These conditions arise from various considerations. Firstly, the effective
gravitational constant ($G_{\rm eff}=G/F$) must be positive, ensuring the
avoidance of anti-gravity (as expressed in Eq. 20). Besides, any acceptable
$f(R)$-model should exhibit stability under perturbations and avoid the
instability of Dolgov-Kawasaki type (Eq. 20). Similar to the $\Lambda$CDM
model, a viable $f(R)$ model must also be consistent with local gravity tests
(Eqs. 20 and 20). Furthermore, the existence of a matter-dominated epoch in
cosmological dynamics necessitates that an $f(R)$ model satisfies the
conditions given by Eqs 20 and 20. Lastly, Eq. 20 ensures the stability of the
late-time de-Sitter solution, from which the late-time accelerated expansion
of the Universe is usually inferred.

The viability conditions mentioned above can be more easily assessed for an
$f(R)$-model if we can somehow transform that model into the following form:

$f(R)=R-2\Lambda{y(R,b,\Lambda)}\,,$ (22)

where the function ${y(R,b,\Lambda)}$ serves to measure the deviation from the
$\Lambda$CDM model, with the parameter $b$ (referred to as the
“distortion/deviation parameter”) specifically quantifying the extent of the
deviation Basilakos et al. (2013). For this reason, all the models examined in
this study are recast into the form of Eq. 22, allowing for a clearer
evaluation of their compliance with the viability conditions.

More specifically, the viability conditions in Eqs. 20 and 20 can be
alternatively written as $\displaystyle\lim_{R\to\infty}f(R)=R+C_{0}$, where
$C_{0}$ represents a constant value. In order for a candidate $f(R)$-model to
exhibit asymptotic behavior similar to the standard $\Lambda$CDM model (which
is supported by observations of the cosmic microwave background), we can
identify that $C_{0}=-2\Lambda$, where $\Lambda$ corresponds to the
cosmological constant introduced by Einstein and Hilbert in their action. This
amounts to writing symbolically $\Lambda^{\Lambda}=\Lambda^{f(R)}$ where the
superscripts $\Lambda$ and $f(R)$ denote quantities in reference to the
$\Lambda$CDM model and any viable $f(R)$ \- model, respectively. This relation
may also be written as

$\Omega_{\Lambda,0}^{\Lambda}(H_{0}^{\Lambda})^{2}=\Omega_{\Lambda,0}^{f(R)}(H_{0}^{f(R)})^{2}\,,$
(23)

where subscript ‘0’ denotes values at present time. Furthermore, considering
the definition of the energy-momentum tensor and the resulting conservation
equation(Eq. 8), we can infer that both the $\Lambda$CDM model and any viable
$f(R)$-model yield identical matter density and radiation density at the
current epoch i.e.

$\Omega_{i,0}^{\Lambda}(H_{0}^{\Lambda})^{2}=\Omega_{i,0}^{f(R)}(H_{0}^{f(R)})^{2}=\frac{8\pi{G}}{3}\rho_{i,0}\,,$
(24)

where the subscript $i=(\mbox{m, r})$ stands for (matter, radiation). Also any
viable $f(R)$-model is expected to exhibit deviations from the standard
$\Lambda$CDM model predominantly during the late times, in order that the
explanation for accelerated expansion at present epoch comes from vanishing
$\Lambda$. In general terms, this implies

$\Omega_{i,0}^{\Lambda}\neq\Omega_{i,0}^{f(R)}\,,\quad H_{0}^{\Lambda}\neq
H_{0}^{f(R)}\,.$ (25)

Based on Eqs. 23-25 and the initial condition requirements for solving the ODE
system described in Eqs. 7, 10, and 11 (see Appendix A), we utilize parameters
($\Omega_{\rm m0}^{\Lambda},\,b,\,H_{0}^{\Lambda}$) for model-fitting to the
data. However, while reporting our findings, we express the results in terms
of the parameters ($\Omega_{\rm m0}^{f(R)},\,b,\,H_{0}^{f(R)}$), which are
determined using Eq. 24.

###  III.1 The Hu-Sawicki Model

The Hu-Sawicki model, initially proposed in Hu and Sawicki (2007), is
described by the following equation:

$f(R)_{\rm HS}=R-\mu^{2}\frac{c{1}(R/\mu^{2})^{n_{{\rm
HS}}}}{1+c{2}(R/\mu^{2})^{n_{{\rm HS}}}}\,,$ (26)

where $c_{1}$ and $c_{2}$ are dimensionless parameters, ${n_{{\rm HS}}}$ is a
positive constant typically assumed to be an integer, and
$\mu^{2}\approx\Omega_{\rm m0}H_{0}^{2}$. By defining
$\mu^{2}c_{1}/2c_{2}\equiv\Lambda$ and $2\left(c_{2}^{1-1/{n_{{\rm
HS}}}}\right)/c_{1}\equiv b$, we can express Eq. 26 as Basilakos et al.
(2013):

$f(R)_{\rm
HS}=R-2\Lambda\left[1-\left\\{1+\left(\frac{R}{b\Lambda}\right)^{{n_{\rm
HS}}}\right\\}^{-1}\right]\,.$ (27)

We identify the parameter $\Lambda$ as the usual cosmological constant, and
$b$ as the deviation parameter, which indicates the model’s deviation from the
$\Lambda$CDM model. In order to satisfy the viability conditions $F>0$ and
$F^{\prime}>0$ for $R\geq{R}_{0}$, it is necessary to consider $b>0$
(especially when ${n_{{\rm HS}}}$ is an odd integer). Although many
researchers have constrained the scenarios where ${n_{{\rm HS}}}=1$ and/or $2$
Basilakos et al. (2013); Nunes et al. (2017); D’Agostino and Nunes (2019,
2020); Farrugia et al. (2021); Leizerovich et al. (2022); Sultana et al.
(2022), they have acknowledged computational challenges as a hindrance to
explore the case of ${n_{{\rm HS}}}=3$ (further elaboration on this point is
provided in Sec. V). In this investigation, we have also imposed constraints
on the case ${n{{}_{\rm HS}}}=3$.

###  III.2 The Starobinsky Model

The model proposed by Starobinsky Starobinsky (2007) is

$f(R)_{\rm ST}=R-{\lambda}R_{\rm S}\left[1-\left(1+\frac{R^{2}}{R_{\rm
S}^{2}}\right)^{-n_{{}_{\rm S}}}\right]\,,$ (28)

where ${n_{{}_{\rm S}}}$ is a positive constant, $\lambda(>0)$ and $R_{\rm
S}\approx R_{0}$ are free parameters (where $R_{0}$ denotes the Ricci scalar
at present epoch). This model too can be reformulated in a more general form
as Basilakos et al. (2013)

$f(R)_{\rm
ST}=R-2\Lambda\left[1-\left\\{1+\left(\frac{R}{b\Lambda}\right)^{2}\right\\}^{-{{n_{{}_{\rm
S}}}}}\right]$ (29)

with $\Lambda={\lambda}R_{\rm S}/2$ and $b=2/\lambda$. In this study, we have
imposed constraints on the cases where ${{n_{{}_{\rm S}}}}=1$ and the reason
for not exploring higher values of ${{n_{{}_{\rm S}}}}$ is explained later in
Sec. V. Note that the Hu-Sawicki model (Eq. 27) with ${{n_{{}_{\rm HS}}}}=2$
and the Starobinsky model (Eq. 29) with ${{n_{{}_{\rm S}}}}=1$ are equivalent.
Unlike the Hu-Sawicki model (with $n_{{}_{\rm HS}}=1$), the viability
condition for the Starobinsky model does not require $b>0$. Based on the
algebraic form of the Starobinsky model (Eq. 29), we can infer that regardless
of the data used, the parameter $b$ must exhibit a symmetric distribution
(centered around $b=0$) from the MCMC fitting procedure. Since our interest
lies in investigating deviations from the $\Lambda$CDM model supported by the
data, we considered $b>0$ without loss of generality.

###  III.3 The Exponential Model

The exponential model, initially proposed as a viable $f(R)$ model in Cognola
et al. (2008), has been further investigated in Linder (2009); Chen et al.
(2015); Odintsov et al. (2017); Leizerovich et al. (2022) (also see references
therein), is given by

$f(R)_{\rm E}=R+\alpha\left[\exp(-\beta{R})-1\right]\,,$ (30)

where $\alpha$ and $\beta$ are the parameters of this model. For large $R$, an
acceptable $f(R)$-model must approximately resemble the $\Lambda$CDM, which is
achievable only when $\alpha>0$ and $\beta>0$. By substituting
$\Lambda=\alpha/2$ and $b=2/(\alpha\beta)$ the exponential model can be
expressed as

$f(R)_{\rm
E}=R-2\Lambda\left[1-\exp\left(-\frac{R}{b\Lambda}\right)\right]\,,$ (31)

from where it becomes evident that when $R$ becomes significantly larger than
$b\Lambda$ ($R\gg{b}\Lambda$), the function $f(R)_{\rm E}$ approaches
$R-2\Lambda$.

###  III.4 The Tsujikawa Model

Tsujikawa proposed an alternative model Tsujikawa (2008) as

$f(R)_{\rm T}=R-{\xi}R_{\rm T}\tanh\left(\frac{R}{R_{\rm T}}\right)\,,$ (32)

where $\xi(>0)$ and $R_{\rm T}(>0)$ are the model parameters. With
$\Lambda={\xi}R_{\rm T}/2$ and $b=2/\xi$, the model can be rewritten as

$f(R)_{\rm T}=R-2\Lambda\tanh\left(\frac{R}{b\Lambda}\right)\,.$ (33)

We see clearly that when the parameter $b\to 0$ (which corresponds to
$\xi\to\infty$, $R_{\rm T}\to 0$, while ${\xi}R_{\rm T}$ remains finite), the
model reduces to $f(R)_{\rm T}=R-2\Lambda$.

###  III.5 The ArcTanh Model

In this study we also examined a model proposed in Pérez-Romero and Nesseris
(2018) as

$f(R)_{\rm
aTanh}=R-\frac{2\Lambda}{1+b\,\text{arctanh}\left(\frac{\Lambda}{R}\right)}\,,$
(34)

where the parameter $b$ is required to be positive, in order to prevent any
occurrence of future singularities.

##  IV Observed Cosmological Data

In this section we present a concise overview of the cosmological data sets
utilized in this study. For any given $f(R)$ model, the system of Eqs. 7, 10
and 11 is solved, primarily yielding the function $H(z)$. So, it becomes
necessary to outline the theoretical equations connecting $H(z)$ with various
observed quantities. Additionally, at the end of this section, a brief
introduction is provided on the statistical techniques employed to extract
model parameters from the data.

###  IV.1 Type Ia Supernova Data

Type Ia Supernovae (SNIa), known as standard candles Phillips (1993), have
significantly contributed to our comprehension of cosmology. The SNIa
observations Perlmutter et al. (1999); Riess et al. (1998) played a pivotal
role in the discovery of the accelerated expansion of the present-day
universe. In this study we used the apparent magnitude data for SNIa obtained
from the recently released PantheonPlus compilation Scolnic et al. (2022).
This compilation comprises of 1701 distinct light curves of 1550 unique
spectroscopically confirmed SNIa, sourced from 18 surveys. This compilation
provides SNIa data within the range $0.00122<z_{\rm HD}<2.26137$, where
$z_{\rm HD}$ denotes the Hubble diagram redshift. It offers a considerably
larger number of low-redshift data compared to the previous Pantheon
compilation Scolnic et al. (2018). In our current work, we employ the apparent
magnitude at maximum brightness ($m_{\rm b}$), heliocentric redshift ($z_{\rm
hel}$), cosmic microwave background corrected redshift ($z_{\rm cmb}$), and
the total (statistical + systematic) covariance matrix from this compilation
Scolnic et al. (2022).

The theoretical definition of the apparent magnitude involves the luminosity
distance, which is given by the integral expression:

$\displaystyle d_{\rm L}(z)$ $\displaystyle=$ $\displaystyle c(1+z_{\rm
hel})\int_{0}^{z_{\rm cmb}}\frac{dz^{\prime}}{H(z^{\prime})}\,,$ (35)

where the function $H(z)$ is obtained by solving the system of Eqs. 7, 10 and
11. The apparent magnitude is defined as

$m_{\rm
th}=M+5\log_{10}\left(\frac{c/H_{0}}{\text{Mpc}}\right)+5\log_{10}\left(D_{\rm
L}(z)\right)+25\,,$ (36)

where $M$ represents the absolute magnitude, and $D_{\rm L}(z)\equiv
H_{0}d_{\rm L}(z)/c$ is the dimensionless Hubble-free luminosity distance. The
computation of apparent magnitude involves a degeneracy between the absolute
magnitude ($M$) and the Hubble parameter at the current epoch ($H_{0}$), as
evident from the Eq. 36. Therefore, after marginalizing over these nuisance
parameters ($M$ and $H_{0}$), the appropriate residual that needs to be
minimized for model fitting is given by Conley et al. (2011)

$\tilde{\chi}^{2}_{\rm sn}=A-\frac{B^{2}}{D}+\log\frac{D}{2\pi}\,,$ (37)

where $A=(m_{\rm b}-m_{\rm th})^{T}C^{-1}(m_{\rm b}-m_{\rm th})$, $B=(m_{\rm
b}-m_{\rm th})^{T}C^{-1}\boldsymbol{1}$, and
$D=\boldsymbol{1}^{T}C^{-1}\boldsymbol{1}$ and $C$ represents the total
covariance matrix of the data (provided in the Pantheon+ compilation).
$\boldsymbol{1}$ is an array of ones of length equal to the number of data
points.

###  IV.2 Cosmic Chronometers Data

By examining the differential age evolution Jimenez and Loeb (2002); Wei
(2018); Simon et al. (2005) of old elliptical galaxies, where star formation
and interactions with other galaxies have ceased, previous studies have
provided 32 data points for $H(z)$ in the redshift range of 0.07-1.965 Simon
et al. (2005); Moresco et al. (2016); Moresco (2015); Zhang et al. (2014);
Stern et al. (2010); Moresco et al. (2012); Ratsimbazafy et al. (2017); Borghi
et al. (2022) (compiled in the Table 4 of the Appendix). This so called
differential age method, employed in these studies, uses the relation

$H(z)=-\frac{1}{1+z}\frac{dz}{dt}\simeq-\frac{1}{1+z}\frac{\Delta{z}}{\Delta{t}}\,,$
(38)

to obtain $H(z)$. The parameters of any model are estimated by minimizing the
following residual:

$\chi^{2}_{{}_{\rm CC}}=\sum_{i=1}^{32}\frac{\left(H_{\rm obs}(z_{i})-H_{\rm
th}(z_{i})\right)^{2}}{\sigma_{H,i}^{2}}\,,$ (39)

where $H_{\rm obs}(z_{i})$’s are the observed values of the Hubble parameter
function at redshift $z_{i}$, while $\sigma_{H,i}^{2}$ denotes the
corresponding uncertainties associated with the measurements of $H_{\rm
obs}(z_{i})$. The theoretical Hubble function at redshift $z_{i}$, which is
model-dependent and obtained from the solutions of Eqs. 7, 10, and 11, is
denoted by $H_{\rm th}(z_{i})$ in the above Eq. 39.

###  IV.3 BAO Data

In the Big Bang Model of the universe, prior to decoupling of matter and
radiation components, the contents of the Universe were evenly distributed,
albeit with small fluctuations. Photons and baryons were strongly coupled
through Thomson scattering. As the universe expanded and cooled, resulting in
a decrease in temperature and density, the fluctuations were amplified by
gravity. The gravitational pull caused the tightly bound photon-baryon mixture
to condense in regions with higher densities, resulting in compressions and
rarefactions in the form of acoustic waves known as Baryonic Acoustic
Oscillations (BAO). Matter and radiation were then decoupled and this epoch
which is marked by the release of baryons from the Compton drag of photons is
known as drag epoch ($z_{\rm d}$), after which the photons travelled freely,
whereas acoustic waves remained frozen in the baryons. The length scale
characterizing the maximum distance traveled by the acoustic wave before
decoupling is known as sound horizon at the epoch of drag ($r_{\rm d}$). BAO,
therefore, holds the status of standard ruler for length scale in Cosmology Hu
and Sugiyama (1996); Eisenstein and Hu (1998).

We have compiled a collection of 30 data points representing various BAO
observables from a range of surveys, as documented in the literature Beutler
et al. (2011); Percival et al. (2010); Padmanabhan et al. (2012); Tojeiro et
al. (2014); Anderson et al. (2013); Bautista et al. (2018); Seo et al. (2012);
Sridhar et al. (2020); de Sainte Agathe et al. (2019); Carter et al. (2018);
Ross et al. (2015); Kazin et al. (2014); Ata et al. (2018); Abbott et al.
(2019); Bautista et al. (2017); du Mas des Bourboux et al. (2017); Font-Ribera
et al. (2014); Alam et al. (2017); Bautista et al. (2020); Neveux et al.
(2020). These data points, which are used in our current study, are listed in
the Table 3(in the Appendix). In our calculations for the drag epoch $z_{d}$
and the sound horizon at the epoch of drag $r_{d}$, we employ improved fits
from Aizpuru et al. (2021). The BAO observables include the Hubble distance,
which is defined as $D_{\rm H}=c/H(z)$, the transverse comoving distance
($D_{\rm M}(z)$), the angular diameter distance ($D_{\rm A}(z)$), and the
volume-averaged distance ($D_{\rm V}(z)$) which are defined as

$D_{A}(z)=\frac{c}{1+z}\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})}\,,$ (40)
$D_{M}(z)=(1+z)D_{A}(z)\,,$ (41)

and,

$D_{V}(z)=\left[(1+z)^{2}D_{A}^{2}(z)\frac{z}{H(z)}\right]^{1/3}\,.$ (42)

Once again $H(z)$ in the above relations is obtained by solving Eqs. 7, 10,
and 11.

For fitting any model to the uncorrelated data points of the BAO data, the
following residual has been used

$\chi^{2}_{{}_{\rm BAO-UC}}=\sum_{i=1}^{20}\left[\frac{A_{\rm
th}(z_{i})-A_{\rm obs}(z_{i})}{\sigma_{i}}\right]^{2}\,,$ (43)

where $A_{\rm obs}(z_{i})$ and $\sigma_{i}$ respectively denote the observed
values and their uncertainties at redshift $z_{i}$ and $A_{\rm th}$ denotes
the theoretical prediction from model under consideration. These quantities
are given in the columns 2-4 of the Table 3. For correlated data points the
appropriate residual to be minimized is

$\chi^{2}_{{}_{\rm BAO-C}}=\sum_{j=1}^{4}\left[\left({\bf A}_{\rm th}-{\bf
A}_{\rm obs}\right)_{j}^{T}{\bf C}_{j}^{-1}\left({\bf A}_{\rm th}-{\bf A}_{\rm
obs}\right)_{j}\right]\,,$ (44)

where ${\bf C}_{j}$’s are the covariance matrices of the 4 different data-sets
and $\left({\bf A}_{\rm th}-{\bf A}_{\rm obs}\right)_{j}$ denotes an array of
difference between observed and theoretical values for each of the 4 data-
sets. These data and the covariance matrix can be found from the corresponding
source papers referred in the last column of Table 3. For the whole data-set
of BAO, the $\chi^{2}$ which is to be minimized is

$\chi^{2}_{{}_{\rm BAO}}=\chi^{2}_{{}_{\rm BAO-UC}}+\chi^{2}_{{}_{\rm
BAO-C}}\,.$ (45)

###  IV.4 H ii starburst Galaxies Data

Hii galaxies (HiiG) are massive and compact starburst structures surrounded by
ionized hydrogen gas. Their optical emission spectra exhibit strong and narrow
Balmer ${\rm H}\alpha$ and ${\rm H}\beta$ lines, along with a weak continuum.
The cosmological significance of HiiG observation comes from the empirically
established correlation between the luminosity ($L$) of the ${\rm H}\beta$
lines and the velocity dispersion ($\sigma$, which is a measure of the width
of the spectral lines). This correlation is attributed to the fact that an
increase in the mass of the starburst component leads to a simultaneous
increase in the number of ionized photons (and thus the luminosity $L$) and
the turbulent velocity (hence, the velocity dispersion $\sigma$) (see
González-Morán et al. (2021) and references therein).

In this study, we have used a total of 181 data points from Erb et al. (2006);
Maseda et al. (2014); Masters et al. (2014); González-Morán et al. (2019,
2021) corresponding to the emission of the Balmer $\rm{H}\beta$ line. These
data points span a redshift range of $0.0088<z<2.5449$. It is important to
note that the redshift coverage provided by the latest SNIa data is as
follows: we have 8 data points for $z_{\rm SN}>1.4$, while for $z_{\rm
H\textsc{ii}G}>1.4$ we have 69 data points, and for $z_{\rm
H\textsc{ii}G}>z_{\rm SN,\,max}$ we have 11 observations of HiiG.
Consequently, the observations of HiiG not only explore hitherto unexplored
higher redshift regions compared to SNIa but also provide a denser coverage of
these higher redshift regions in comparison to SNIa.

The empirical relation between the luminosity $L$ and the dispersion of
velocity $\sigma$ is given by

$\log\left[\frac{L}{\mathrm{erg/s}}\right]=\beta\log\left[\frac{\sigma}{\mathrm{km/s}}\right]+\alpha\,,$
(46)

where we take $\beta=5.022\pm 0.058$ and $\alpha=33.268\pm 0.083$ from
González-Morán et al. (2019, 2021). From the luminosity distance of HiiG,
$d_{\rm L}=\left[L/(4{\pi}F)\right]^{1/2}$, we can derive the distance modulus
as

$\mu_{o}=2.5\left(\alpha+\beta\log\left[\frac{\sigma}{\mathrm{km/s}}\right]-\log\left[\frac{\sigma}{\mathrm{erg/s/cm^{2}}}\right]-40.08\right).$
(47)

Using Eq. 47, we compute the distance moduli of HiiG from the observables $L$,
$\sigma$, and $F$. For any given cosmological model with theoretical distance
moduli $\mu_{\theta}$ (represented as $m_{\rm th}-M$ in Eq. 36), the
parameters associated with that model can be constrained by minimising the
following $\chi^{2}$ function

$\chi^{2}_{{}_{\mathrm{H}\textsc{ii}\mathrm{G}}}=\sum_{i=1}^{181}\frac{\left[\mu_{o}-\mu_{\theta}\right]^{2}}{\epsilon^{2}}\,,$
(48)

where

$\epsilon^{2}=\epsilon_{\mu_{o,\mathrm{stat}}}^{2}+\epsilon_{\mu_{\theta,\mathrm{stat}}}^{2}+\epsilon_{\mathrm{sys}}^{2}\,,$
(49)

with

$\epsilon_{\mu_{o,\mathrm{stat}}}^{2}=6.25\left(\epsilon_{\log{F}}^{2}+\beta^{2}\epsilon_{\log{\sigma}}^{2}+\epsilon_{\beta}^{2}(\log{\sigma})^{2}+\epsilon_{\alpha}^{2}\right)\,,$
(50) $\epsilon_{\mu_{\theta,\mathrm{stat}}}^{2}=\left[\frac{5}{\ln
10}\left(\frac{c(1+z)}{d_{L}(z)H(z)}+\frac{1}{1+z}\right)\sigma_{z}\right]^{2},$
(51)

and $\epsilon_{\mathrm{sys}}=0.0.257$ as suggested in Chávez et al. (2016);
González-Morán et al. (2019, 2021). The uncertainty contribution of
$\epsilon_{\mu_{\theta,\mathrm{stat}}}$ in the distance modulus is due to the
uncertainties in redshifts($\sigma_{z}\sim 10^{-4}$) of HiiG and is derived
from simple error propagation theory.

###  IV.5 Local Measurement of $H_{0}$

A contentious issue remains regarding the current value of the Hubble
parameter $H_{0}$. The Planck constraints provide a value of $H_{0}=67.4\pm
0.5$ Aghanim et al. (2020), while the locally measured value by the SH0ES
collaboration yields $H_{0}=73.04\pm 1.04$ Riess et al. (2022). The
discrepancy between these two measurements, known as the Hubble tension, is
significant at the $4-5\sigma$ level. To address this tension, we have
included in our study both scenarios - with vis-a-vis without the SH0ES prior
for $H_{0}$. The contribution of this data point to the total $\chi^{2}$ is
given by

$\chi^{2}_{\rm SH0ES}=\frac{(H_{0}-73.04)^{2}}{1.04^{2}}\,.$ (52)

###  IV.6 Methodology

We employed the well-known Markov Chain Monte Carlo (MCMC) analysis for
parameter estimations of the models. This involves maximizing the likelihood
function, given by

$\mathcal{L}=\exp\left(-\sum_{i}\chi^{2}_{i}/2\right)\,,$ (53)

where the subscript ‘$i$’, depending on the combination of different data-sets
used, stands for one or more from the set: $\\{$ ‘SN’, ‘CC’, ‘BAO’, ‘HiiG’,
‘SH0ES’$\\}$. These $\chi^{2}_{i}$’s are defined in Eqs. 37, 39, 45, 48 and
52. To get posterior probability distribution of the model parameters we also
need to set priors on them. For the reasons mentioned in the last section and
will further be discussed in the Appendix A, for all the $f(R)$ models in this
work we do data-fitting in terms of the parameters: $(\Omega_{\rm
m0}^{\Lambda},b,H_{0}^{\Lambda})$. We used the uniform priors $\Omega_{\rm
m0}^{\Lambda}\in[0,1]$ and $H_{0}^{\Lambda}\in[50,90]$ for all the models,
whereas priors for $b$ require further considerations. From the algebraic
expressions of these six $f(R)$ \- models we worked out and established a
hierarchy of similarity of these models with the standard $\Lambda$CDM model
using some “measures-of-similarity”(e.g. mean-square-error, correlation,
dynamic time warping etc.). A ballpark hierarchy of similarity we find is: the
most similar $f(R)$ models are the Tsujikawa model(TSUJI) and the Hu-Sawicki
model($n_{{}_{\rm HS}}=3$)(HS3), followed by the Starobinsky($n_{{}_{\rm
S}}=1$)(ST1) and the exponential(EXP) models, and the least similar ones are
the Hu-Sawicki model($n_{{}_{\rm HS}}=1$)(HS1) and the arcTanh(aTanh) models.
To allow for the exploration that whether the data supports a model different
from the $\Lambda$CDM model, an $f(R)$ model which is more similar to the
$\Lambda$CDM requires a bigger change in $b$ to make them noticeably different
from the $\Lambda$CDM model. Consequently we use uniform priors $b\in[0,b_{\rm
max}]$ where $b_{\rm max}$ is 7, 7, 5, 5, 3 and 3 for the models HS3, TSUJI,
ST1, EXP, HS1 and aTanh, respectively.

We developed our own PYTHON codes using the publicly available PYTHON modules:
(i) for solving stiff ODEs Wogan (2021), (ii) to perform MCMC analysis — EMCEE
Foreman-Mackey et al. (2013), and (iii) for plotting of posterior probability
distributions of the parameters — GetDist Lewis (2019).

##  V RESULTS: Observational constraints on $f(R)$ models

Figure 1: This figure is helpful in getting a sense of variations of
$\Omega_{m0}$ from model to model and with data for a given model. Different
colors represent different data-set combinations and these are same as in any
of the parameter distribution plots(e.g. Fig. 6) or see 2nd paragraph of
Sec.V. The blank star and blank circle markers represent median values for the
cases with and without SH0ES prior for $H_{0}$, respectively. The thick and
thin horizontal colored bars represent 1-sigma(68.26%) and 2-sigma(95.44%)
confidence intervals for the cases with SH0ES prior. The colored
continuous/dashed lines represent 1-sigma(68.26%, with smaller cap size) and
2-sigma(95.44%, with bigger cap size) confidence intervals for the cases
without SH0ES prior. Figure 2: This figure is helpful in getting a sense of
variations of the deviation parameter $b$ from model to model and with data
for a given model. Different colors represent different data-set combinations
and these are same as in any of the parameter distribution plots(e.g. Fig. 6)
or see 2nd paragraph of Sec.V. The blank star and blank circle markers
represent median values for the cases with and without SH0ES prior for
$H_{0}$, respectively. The thick and thin horizontal colored bars represent
1-sigma(68.26%) and 2-sigma(95.44%) confidence intervals for the cases with
SH0ES prior. The colored continuous/dashed lines represent 1-sigma(68.26%,
with smaller cap size) and 2-sigma(95.44%, with bigger cap size) confidence
intervals for the cases without SH0ES prior. Figure 3: Hubble Tension:
Different colors represent different data-set combinations and these are same
as in any of the parameter distribution plots(e.g. Fig. 6) or see 2nd
paragraph of Sec.V. The blank star and blank circle markers represent median
values for the cases with and without SH0ES prior for $H_{0}$, respectively.
The thick and thin horizontal colored bars represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals for the cases with SH0ES prior. The
colored continuous/dashed lines represent 1-sigma(68.26%, with smaller cap
size) and 2-sigma(95.44%, with bigger cap size) confidence intervals for the
cases without SH0ES prior.

In this section we will discuss our findings regarding the parameters of the
six $f(R)$ models and the $\Lambda$CDM model. We have considered five
combinations of data sets: SNIa+CC (SC), SNIa+BAO+CC (SBC), SNIa+CC+HiiG
(SCHii), BAO+CC+HIIG (BCHii), and SNIa+BAO+CC+HIIG (SBCHii). We also
separately included a SH0ES prior for the Hubble parameter $H_{0}$ for each of
these combinations. While each data set individually can constrain any model,
they may not necessarily provide tight constraints on the parameters. The
choice of data combinations is also guided by considerations of goodness-of-
fit.

#### Acronym and Color convention:

In addition to the above mentioned acronyms for different data sets, when
including the SH0ES prior for $H_{0}$, we append the notations as follows:
SC$H_{0}$ to denote SNIa+CC+$H_{0}$ and SC($H_{0})$ to denote SNIa+CC or
SNIa+CC+$H_{0}$, depending on the case. Similar notations are used for the
other four data sets. Unless otherwise specified, the following color and data
set correspondences are used in the upcoming figures: (i) SC($H_{0})$—Blue,
(ii) SBC($H_{0})$—Red, (iii) SCH$\textsc{ii}(H_{0})$—Black, (iv)
BCH$\textsc{ii}(H_{0})$—Orange, and (v) SBCH$\textsc{ii}(H_{0})$—Green.

We generated MCMC samples of size 875,000 (with 25 walkers, each taking 35,000
steps) for each parameter in all the cases, where a case refers to a specific
combination of a given model and a data set. These raw MCMC sample chains
underwent tests for convergence and independence. To ensure convergence, an
initial portion of the chain was discarded to obtain a “burned chain.” We
removed the first 125,000 steps (i.e., 5,000 initial steps of each walkers) to
obtain the burned chain. In order to have independent and uncorrelated
samples, the burned chains needed to be properly thinned. We thinned the
burned chains by a factor of 0.75 times the integrated auto-correlation time.
As a result, depending on the case, we obtained convergent and independent
samples of varying sizes, approximately ranging from 15,000 to 20,000. All the
statistical inferences about the model parameters were then obtained from
these burned and thinned subsamples.

The median values of model parameters and 1-sigma confidence intervals on them
are presented in the Tables 1 and 2. The Figs 4-17 display 2D contour plots
depicting the posterior probability distribution of the parameters, which
provide an indication of potential correlations among the parameters, as well
as 1D marginalised distributions of each parameter are also shown there.

Before getting into the detailed results for individual models, we would like
to highlight some overall observations regarding the model parameters in the
light of Figs. 1, 2, and 3. (i) $\Omega_{m0}$: For all models and data-set
combinations, either the median values of $\Omega_{m0}$ fall within the 1-2
sigma interval of the Planck value ($\Omega_{m0,\rm Planck}=0.315\pm 0.007$
Aghanim et al. (2020)) or, the Planck value is within the 1-2 sigma intervals
of the model’s median values of $\Omega_{m0}$ (for cases with or without the
SH0ES prior for $H_{0}$). (ii) $b$: There is a general trend of a shift
towards lower values of $b$ when considering the SH0ES prior for $H_{0}$
compared to the corresponding cases without the SH0ES prior for the
SBC$(H_{0})$, BCH$\textsc{ii}(H_{0})$, and SBCH$\textsc{ii}(H_{0})$ data-sets,
whereas opposite trends are observed for the SC$(H_{0})$ and
SCH$\textsc{ii}(H_{0})$ data-sets. (iii) $H_{0}$: Without the SH0ES prior for
$H_{0}$, for all models and data-set combinations, either the median values of
$H_{0}$ are within a 1-2 sigma interval of the Planck value ($H_{0,\rm
Planck}=67.4\pm 0.5$ Aghanim et al. (2020)) or, the Planck value is within a
1-2 sigma interval of the model’s median values of $H_{0}$. However, with the
SH0ES prior for $H_{0}$, there are slight departures towards higher values of
$H_{0}$ (as expected), but they are not close to $H_{0,\rm SH0ES}=73.04\pm
1.04$ Riess et al. (2022). It is worth noting that the cases with the SH0ES
prior for $H_{0}$ generally have tighter bounds on the model’s $H_{0}$ values
(this observation does not hold true for $\Omega_{m0}$ or $b$). Additionally,
there is less overall tension among the model-fitted values of $H_{0}$
compared to the tension between $H_{0,\rm Planck}$ and $H_{0,\rm SH0ES}$.

Table 1: Results from the MCMC fitting process(for the cases without SH0ES prior for $H_{0}$): The median values of the model parameters along with 1-sigma(68.26%) confidence intervals on them, the transition redshift($z_{\rm t}$), the current values of EoS parameter for geometric/dark energy component($w_{\rm DE0}$) and total EoS parameter($w_{\rm tot0}$), the minimum values of $\chi^{2}$ and other quantities for making statistical inferences such as the reduced chi-square value $\chi^{2}_{\rm red}=\chi^{2}_{\rm min}/\nu$(where the number of degrees of freedom $\nu=N-k$, $N$ is the total number of data points and $k$ is the number of model parameters), the Akaike Information Criterion(AIC), the Bayesian Information Criterion(BIC), $\Delta$AIC = AIC${}_{f(R)\rm\,\,model}$ \- AIC${}_{\Lambda\rm{CDM}\,\,model}$ and $\Delta$BIC = BICf(R)model \- BIC${}_{\Lambda\rm{CDM}\,\,model}$ — for all the models and data set combinations considered in this work. Data | $\Omega_{m0}$ | $b$ | $H_{0}$ | $z_{\rm t}$ | $w_{\rm DE0}$ | $w_{\rm tot0}$ | $\chi^{2}$ | $\chi^{2}_{\rm red}$ | AIC | BIC |  $\Delta$AIC |  $\Delta$BIC
---|---|---|---|---|---|---|---|---|---|---|---|---
$\Lambda$CDM |  |  |  |  |  |  |  |  |  |  |  |
SC | $0.36_{-0.018}^{+0.018}$ | — | $66.42_{-1.674}^{+1.703}$ | $0.53_{-0.039}^{+0.040}$ | $-1$ | $-0.64_{-0.018}^{+0.018}$ | 1777.31 | 1.03 | 1781.31 | 1792.22 | 0.0 | 0.0
SBC | $0.32_{-0.010}^{+0.010}$ | — | $69.20_{-0.669}^{+0.684}$ | $0.63_{-0.025}^{+0.025}$ | $-1$ | $-0.68_{-0.010}^{+0.010}$ | 1815.39 | 1.03 | 1819.39 | 1830.34 | 0.0 | 0.0
SCHii | $0.34_{-0.016}^{+0.017}$ | — | $68.21_{-1.179}^{+1.182}$ | $0.56_{-0.037}^{+0.039}$ | $-1$ | $-0.66_{-0.016}^{+0.017}$ | 2145.62 | 1.12 | 2149.62 | 2160.73 | 0.0 | 0.0
BCHii | $0.30_{-0.011}^{+0.011}$ | — | $68.85_{-0.608}^{+0.619}$ | $0.68_{-0.029}^{+0.029}$ | $-1$ | $-0.70_{-0.011}^{+0.011}$ | 408.22 | 1.69 | 412.22 | 419.19 | 0.0 | 0.0
SBCHii | $0.31_{-0.009}^{+0.010}$ | — | $69.24_{-0.599}^{+0.613}$ | $0.63_{-0.024}^{+0.024}$ | $-1$ | $-0.69_{-0.009}^{+0.010}$ | 2180.25 | 1.12 | 2184.25 | 2195.39 | 0.0 | 0.0
HS1 |  |  |  |  |  |  |  |  |  |  |  |
SC | $0.27_{-0.041}^{+0.039}$ | $1.46_{-0.711}^{+0.662}$ | $67.10_{-1.704}^{+1.721}$ | $0.78_{-0.132}^{+0.188}$ | $-0.78_{-0.078}^{+0.070}$ | $-0.57_{-0.033}^{+0.031}$ | 1773.75 | 1.03 | 1779.75 | 1796.11 | -1.56 | 3.9
SBC | $0.30_{-0.011}^{+0.011}$ | $0.93_{-0.303}^{+0.368}$ | $64.99_{-1.145}^{+1.175}$ | $0.68_{-0.035}^{+0.037}$ | $-0.84_{-0.041}^{+0.040}$ | $-0.59_{-0.027}^{+0.027}$ | 1800.65 | 1.02 | 1806.65 | 1823.07 | -12.74 | -7.26
SCHii | $0.25_{-0.035}^{+0.037}$ | $1.64_{-0.676}^{+0.553}$ | $68.67_{-1.213}^{+1.181}$ | $0.87_{-0.145}^{+0.188}$ | $-0.77_{-0.071}^{+0.061}$ | $-0.57_{-0.031}^{+0.029}$ | 2139.44 | 1.12 | 2145.44 | 2162.1 | -4.18 | 1.37
BCHii | $0.30_{-0.012}^{+0.012}$ | $0.23_{-0.161}^{+0.305}$ | $67.73_{-1.165}^{+0.977}$ | $0.67_{-0.033}^{+0.032}$ | $-0.96_{-0.032}^{+0.065}$ | $-0.67_{-0.028}^{+0.050}$ | 408.19 | 1.7 | 414.19 | 424.64 | 1.97 | 5.44
SBCHii | $0.31_{-0.010}^{+0.010}$ | $0.57_{-0.206}^{+0.255}$ | $66.66_{-0.927}^{+0.973}$ | $0.65_{-0.029}^{+0.030}$ | $-0.89_{-0.039}^{+0.037}$ | $-0.61_{-0.027}^{+0.026}$ | 2171.26 | 1.12 | 2177.26 | 2193.97 | -6.99 | -1.42
ST1/HS2 |  |  |  |  |  |  |  |  |  |  |  |
SC | $0.30_{-0.034}^{+0.031}$ | $2.75_{-1.165}^{+0.786}$ | $67.01_{-1.716}^{+1.740}$ | $0.75_{-0.118}^{+0.155}$ | $-0.81_{-0.071}^{+0.085}$ | $-0.56_{-0.035}^{+0.042}$ | 1773.24 | 1.02 | 1779.24 | 1795.61 | -2.07 | 3.39
SBC | $0.31_{-0.011}^{+0.011}$ | $2.64_{-0.754}^{+0.668}$ | $65.72_{-1.052}^{+1.038}$ | $0.72_{-0.040}^{+0.042}$ | $-0.82_{-0.043}^{+0.047}$ | $-0.57_{-0.030}^{+0.032}$ | 1799.22 | 1.02 | 1805.22 | 1821.64 | -14.17 | -8.7
SCHii | $0.28_{-0.033}^{+0.030}$ | $2.99_{-1.035}^{+0.568}$ | $68.62_{-1.190}^{+1.234}$ | $0.84_{-0.132}^{+0.169}$ | $-0.78_{-0.074}^{+0.087}$ | $-0.56_{-0.038}^{+0.046}$ | 2138.76 | 1.12 | 2144.76 | 2161.43 | -4.86 | 0.7
BCHii | $0.31_{-0.013}^{+0.014}$ | $1.04_{-0.729}^{+1.234}$ | $67.62_{-1.130}^{+1.068}$ | $0.69_{-0.032}^{+0.033}$ | $-0.91_{-0.077}^{+0.074}$ | $-0.64_{-0.059}^{+0.060}$ | 408.18 | 1.7 | 414.18 | 424.62 | 1.95 | 5.43
SBCHii | $0.31_{-0.010}^{+0.010}$ | $2.03_{-0.661}^{+0.691}$ | $66.88_{-0.847}^{+0.840}$ | $0.69_{-0.034}^{+0.035}$ | $-0.85_{-0.038}^{+0.039}$ | $-0.59_{-0.027}^{+0.028}$ | 2168.85 | 1.12 | 2174.85 | 2191.57 | -9.4 | -3.83
HS3 |  |  |  |  |  |  |  |  |  |  |  |
SC | $0.33_{-0.028}^{+0.026}$ | $5.02_{-2.092}^{+1.001}$ | $66.84_{-1.708}^{+1.684}$ | $0.70_{-0.094}^{+0.115}$ | $-0.84_{-0.071}^{+0.105}$ | $-0.56_{-0.038}^{+0.057}$ | 1772.23 | 1.02 | 1778.23 | 1794.6 | -3.07 | 2.38
SBC | $0.31_{-0.011}^{+0.011}$ | $5.62_{-1.130}^{+0.643}$ | $66.54_{-0.993}^{+1.030}$ | $0.76_{-0.040}^{+0.040}$ | $-0.79_{-0.064}^{+0.079}$ | $-0.54_{-0.045}^{+0.055}$ | 1798.35 | 1.02 | 1804.35 | 1820.77 | -15.04 | -9.57
SCHii | $0.31_{-0.026}^{+0.026}$ | $5.24_{-1.657}^{+0.673}$ | $68.53_{-1.162}^{+1.216}$ | $0.78_{-0.106}^{+0.113}$ | $-0.80_{-0.083}^{+0.105}$ | $-0.56_{-0.043}^{+0.062}$ | 2137.24 | 1.12 | 2143.24 | 2159.9 | -6.38 | -0.83
BCHii | $0.31_{-0.013}^{+0.013}$ | $2.96_{-2.019}^{+2.134}$ | $67.93_{-0.874}^{+0.824}$ | $0.72_{-0.038}^{+0.039}$ | $-0.90_{-0.072}^{+0.073}$ | $-0.62_{-0.057}^{+0.056}$ | 408.21 | 1.7 | 414.21 | 424.65 | 1.99 | 5.46
SBCHii | $0.32_{-0.010}^{+0.010}$ | $5.01_{-1.389}^{+0.924}$ | $67.42_{-0.883}^{+0.827}$ | $0.73_{-0.037}^{+0.036}$ | $-0.83_{-0.052}^{+0.065}$ | $-0.57_{-0.037}^{+0.046}$ | 2166.71 | 1.12 | 2172.71 | 2189.42 | -11.54 | -5.97
EXP |  |  |  |  |  |  |  |  |  |  |  |
SC | $0.31_{-0.042}^{+0.036}$ | $2.34_{-0.792}^{+0.597}$ | $66.97_{-1.738}^{+1.729}$ | $0.73_{-0.126}^{+0.160}$ | $-0.83_{-0.090}^{+0.093}$ | $-0.56_{-0.041}^{+0.044}$ | 1773.24 | 1.02 | 1779.24 | 1795.61 | -2.07 | 3.39
SBC | $0.31_{-0.011}^{+0.011}$ | $2.42_{-0.390}^{+0.365}$ | $66.05_{-1.102}^{+1.097}$ | $0.78_{-0.042}^{+0.043}$ | $-0.79_{-0.048}^{+0.050}$ | $-0.54_{-0.033}^{+0.034}$ | 1799.21 | 1.02 | 1805.21 | 1821.63 | -14.18 | -8.71
SCHii | $0.28_{-0.041}^{+0.036}$ | $2.55_{-0.659}^{+0.480}$ | $68.68_{-1.185}^{+1.242}$ | $0.83_{-0.137}^{+0.166}$ | $-0.79_{-0.086}^{+0.086}$ | $-0.55_{-0.041}^{+0.043}$ | 2138.92 | 1.12 | 2144.92 | 2161.58 | -4.7 | 0.85
BCHii | $0.30_{-0.013}^{+0.014}$ | $1.16_{-0.812}^{+0.940}$ | $68.21_{-1.208}^{+0.869}$ | $0.72_{-0.046}^{+0.048}$ | $-0.93_{-0.068}^{+0.102}$ | $-0.65_{-0.053}^{+0.080}$ | 408.22 | 1.7 | 414.22 | 424.67 | 2.0 | 5.47
SBCHii | $0.32_{-0.010}^{+0.010}$ | $2.10_{-0.383}^{+0.372}$ | $67.21_{-0.883}^{+0.899}$ | $0.76_{-0.036}^{+0.037}$ | $-0.83_{-0.044}^{+0.044}$ | $-0.56_{-0.031}^{+0.031}$ | 2167.67 | 1.12 | 2173.67 | 2190.38 | -10.58 | -5.01
TSUJI |  |  |  |  |  |  |  |  |  |  |  |
SC | $0.33_{-0.036}^{+0.032}$ | $4.32_{-1.398}^{+0.828}$ | $66.91_{-1.704}^{+1.743}$ | $0.73_{-0.126}^{+0.160}$ | $-0.83_{-0.090}^{+0.093}$ | $-0.56_{-0.041}^{+0.044}$ | 1772.91 | 1.02 | 1778.91 | 1795.28 | -2.39 | 3.06
SBC | $0.31_{-0.011}^{+0.011}$ | $4.74_{-0.649}^{+0.552}$ | $66.36_{-1.067}^{+1.107}$ | $0.78_{-0.042}^{+0.043}$ | $-0.79_{-0.048}^{+0.050}$ | $-0.54_{-0.033}^{+0.034}$ | 1799.12 | 1.02 | 1805.12 | 1821.54 | -14.27 | -8.8
SCHii | $0.30_{-0.035}^{+0.032}$ | $4.60_{-1.019}^{+0.559}$ | $68.58_{-1.196}^{+1.247}$ | $0.83_{-0.137}^{+0.166}$ | $-0.79_{-0.086}^{+0.086}$ | $-0.55_{-0.041}^{+0.043}$ | 2138.49 | 1.12 | 2144.49 | 2161.15 | -5.13 | 0.42
BCHii | $0.31_{-0.013}^{+0.015}$ | $2.35_{-1.663}^{+1.877}$ | $68.33_{-1.069}^{+0.829}$ | $0.72_{-0.046}^{+0.048}$ | $-0.93_{-0.068}^{+0.102}$ | $-0.65_{-0.053}^{+0.080}$ | 408.22 | 1.7 | 414.22 | 424.66 | 2.0 | 5.47
SBCHii | $0.32_{-0.010}^{+0.010}$ | $4.29_{-0.701}^{+0.590}$ | $67.39_{-0.879}^{+0.886}$ | $0.76_{-0.036}^{+0.037}$ | $-0.83_{-0.044}^{+0.044}$ | $-0.56_{-0.031}^{+0.031}$ | 2167.17 | 1.12 | 2173.17 | 2189.88 | -11.08 | -5.51
aTanh |  |  |  |  |  |  |  |  |  |  |  |
SC | $0.27_{-0.042}^{+0.041}$ | $1.69_{-0.856}^{+0.865}$ | $67.08_{-1.707}^{+1.659}$ | $0.79_{-0.137}^{+0.186}$ | $-0.78_{-0.077}^{+0.065}$ | $-0.57_{-0.031}^{+0.027}$ | 1773.73 | 1.03 | 1779.73 | 1796.1 | -1.57 | 3.88
SBC | $0.30_{-0.011}^{+0.012}$ | $0.99_{-0.339}^{+0.412}$ | $65.02_{-1.169}^{+1.158}$ | $0.68_{-0.035}^{+0.036}$ | $-0.84_{-0.039}^{+0.038}$ | $-0.59_{-0.026}^{+0.026}$ | 1800.58 | 1.02 | 1806.58 | 1823.0 | -12.81 | -7.34
SCHii | $0.25_{-0.035}^{+0.036}$ | $1.87_{-0.768}^{+0.721}$ | $68.65_{-1.151}^{+1.188}$ | $0.88_{-0.142}^{+0.177}$ | $-0.76_{-0.065}^{+0.052}$ | $-0.58_{-0.028}^{+0.025}$ | 2139.45 | 1.12 | 2145.45 | 2162.12 | -4.17 | 1.39
BCHii | $0.30_{-0.011}^{+0.012}$ | $0.24_{-0.169}^{+0.334}$ | $67.68_{-1.190}^{+0.986}$ | $0.67_{-0.033}^{+0.031}$ | $-0.95_{-0.037}^{+0.065}$ | $-0.66_{-0.031}^{+0.051}$ | 408.2 | 1.7 | 414.2 | 424.64 | 1.97 | 5.45
SBCHii | $0.31_{-0.010}^{+0.010}$ | $0.59_{-0.225}^{+0.290}$ | $66.67_{-0.938}^{+0.958}$ | $0.65_{-0.029}^{+0.030}$ | $-0.88_{-0.038}^{+0.036}$ | $-0.61_{-0.027}^{+0.026}$ | 2171.09 | 1.12 | 2177.09 | 2193.81 | -7.16 | -1.59
Table 2: Results from the MCMC fitting process(for the cases with SH0ES prior for $H_{0}$): The median values of the model parameters along with 1-sigma(68.26%) confidence intervals on them, the transition redshift($z_{\rm t}$), the current values of EoS parameter for geometric/dark energy component($w_{\rm DE0}$) and total EoS parameter($w_{\rm tot0}$), the minimum values of $\chi^{2}$ and other quantities for making statistical inferences such as the reduced chi-square value $\chi^{2}_{\rm red}=\chi^{2}_{\rm min}/\nu$(where the number of degrees of freedom $\nu=N-k$, $N$ is the total number of data points and $k$ is the number of model parameters), the Akaike Information Criterion(AIC), the Bayesian Information Criterion(BIC), $\Delta$AIC = AIC${}_{f(R)\rm\,\,model}$ \- AIC${}_{\Lambda\rm{CDM}\,\,model}$ and $\Delta$BIC = BICf(R)model \- BIC${}_{\Lambda\rm{CDM}\,\,model}$ — for all the models and data set combinations considered in this work. Data | $\Omega_{m0}$ | $b$ | $H_{0}$ | $z_{\rm t}$ | $w_{\rm DE0}$ | $w_{\rm tot0}$ | $\chi^{2}$ | $\chi^{2}_{\rm red}$ | AIC | BIC |  $\Delta$AIC |  $\Delta$BIC
---|---|---|---|---|---|---|---|---|---|---|---|---
$\Lambda$CDM |  |  |  |  |  |  |  |  |  |  |  |
SC$H_{0}$ | $0.34_{-0.016}^{+0.016}$ | — | $71.28_{-0.887}^{+0.898}$ | $0.58_{-0.037}^{+0.038}$ | $-1$ | $-0.66_{-0.016}^{+0.016}$ | 1788.16 | 1.03 | 1792.16 | 1803.07 | 0.0 | 0.0
SBC$H_{0}$ | $0.32_{-0.010}^{+0.010}$ | — | $70.40_{-0.589}^{+0.591}$ | $0.61_{-0.024}^{+0.025}$ | $-1$ | $-0.68_{-0.010}^{+0.010}$ | 1824.91 | 1.04 | 1828.91 | 1839.86 | 0.0 | 0.0
SCH$\textsc{ii}H_{0}$ | $0.33_{-0.015}^{+0.015}$ | — | $70.98_{-0.794}^{+0.799}$ | $0.60_{-0.036}^{+0.038}$ | $-1$ | $-0.67_{-0.015}^{+0.015}$ | 2154.83 | 1.13 | 2158.83 | 2169.94 | 0.0 | 0.0
BCH$\textsc{ii}H_{0}$ | $0.30_{-0.011}^{+0.011}$ | — | $69.97_{-0.546}^{+0.547}$ | $0.66_{-0.028}^{+0.028}$ | $-1$ | $-0.70_{-0.011}^{+0.011}$ | 420.25 | 1.74 | 424.25 | 431.23 | 0.0 | 0.0
SBCH$\textsc{ii}H_{0}$ | $0.32_{-0.009}^{+0.010}$ | — | $70.25_{-0.541}^{+0.552}$ | $0.62_{-0.023}^{+0.023}$ | $-1$ | $-0.68_{-0.009}^{+0.010}$ | 2190.08 | 1.13 | 2194.08 | 2205.23 | 0.0 | 0.0
HS1 |  |  |  |  |  |  |  |  |  |  |  |
SC$H_{0}$ | $0.24_{-0.031}^{+0.036}$ | $1.55_{-0.630}^{+0.499}$ | $71.43_{-0.889}^{+0.897}$ | $0.89_{-0.145}^{+0.168}$ | $-0.77_{-0.067}^{+0.054}$ | $-0.59_{-0.029}^{+0.027}$ | 1782.1 | 1.03 | 1788.1 | 1804.47 | -4.06 | 1.4
SBC$H_{0}$ | $0.32_{-0.010}^{+0.010}$ | $0.15_{-0.101}^{+0.134}$ | $69.70_{-0.799}^{+0.735}$ | $0.61_{-0.025}^{+0.026}$ | $-0.97_{-0.019}^{+0.032}$ | $-0.66_{-0.016}^{+0.022}$ | 1824.39 | 1.04 | 1830.39 | 1846.81 | 1.48 | 6.95
SCH$\textsc{ii}H_{0}$ | $0.23_{-0.027}^{+0.033}$ | $1.67_{-0.581}^{+0.432}$ | $71.18_{-0.823}^{+0.780}$ | $0.96_{-0.149}^{+0.158}$ | $-0.76_{-0.061}^{+0.048}$ | $-0.59_{-0.029}^{+0.026}$ | 2146.44 | 1.12 | 2152.44 | 2169.1 | -6.39 | -0.84
BCH$\textsc{ii}H_{0}$ | $0.31_{-0.011}^{+0.011}$ | $0.05_{-0.035}^{+0.068}$ | $69.76_{-0.575}^{+0.578}$ | $0.65_{-0.028}^{+0.029}$ | $-0.99_{-0.005}^{+0.011}$ | $-0.69_{-0.013}^{+0.014}$ | 420.28 | 1.74 | 426.28 | 436.73 | 2.03 | 5.51
SBCH$\textsc{ii}H_{0}$ | $0.32_{-0.010}^{+0.010}$ | $0.14_{-0.091}^{+0.121}$ | $69.69_{-0.688}^{+0.642}$ | $0.62_{-0.024}^{+0.025}$ | $-0.98_{-0.016}^{+0.028}$ | $-0.66_{-0.015}^{+0.020}$ | 2189.75 | 1.13 | 2195.75 | 2212.46 | 1.67 | 7.24
ST1/HS2 |  |  |  |  |  |  |  |  |  |  |  |
SC$H_{0}$ | $0.27_{-0.032}^{+0.029}$ | $2.85_{-0.971}^{+0.601}$ | $71.45_{-0.873}^{+0.901}$ | $0.86_{-0.132}^{+0.168}$ | $-0.79_{-0.069}^{+0.084}$ | $-0.57_{-0.035}^{+0.044}$ | 1781.71 | 1.03 | 1787.71 | 1804.08 | -4.45 | 1.01
SBC$H_{0}$ | $0.32_{-0.010}^{+0.010}$ | $0.83_{-0.445}^{+0.578}$ | $69.51_{-0.799}^{+0.800}$ | $0.63_{-0.029}^{+0.030}$ | $-0.94_{-0.050}^{+0.045}$ | $-0.63_{-0.034}^{+0.033}$ | 1823.2 | 1.04 | 1829.2 | 1845.62 | 0.28 | 5.75
SCH$\textsc{ii}H_{0}$ | $0.26_{-0.031}^{+0.030}$ | $3.03_{-0.889}^{+0.440}$ | $71.19_{-0.820}^{+0.802}$ | $0.94_{-0.144}^{+0.179}$ | $-0.77_{-0.074}^{+0.084}$ | $-0.57_{-0.038}^{+0.047}$ | 2145.93 | 1.12 | 2151.93 | 2168.6 | -6.89 | -1.34
BCH$\textsc{ii}H_{0}$ | $0.31_{-0.011}^{+0.011}$ | $0.21_{-0.147}^{+0.238}$ | $69.85_{-0.579}^{+0.577}$ | $0.66_{-0.028}^{+0.028}$ | $-1.00_{-0.003}^{+0.015}$ | $-0.69_{-0.012}^{+0.016}$ | 420.25 | 1.74 | 426.25 | 436.71 | 2.0 | 5.48
SBCH$\textsc{ii}H_{0}$ | $0.32_{-0.009}^{+0.010}$ | $0.75_{-0.427}^{+0.521}$ | $69.53_{-0.723}^{+0.727}$ | $0.64_{-0.027}^{+0.027}$ | $-0.94_{-0.046}^{+0.046}$ | $-0.64_{-0.031}^{+0.033}$ | 2188.66 | 1.13 | 2194.66 | 2211.38 | 0.58 | 6.15
HS3 |  |  |  |  |  |  |  |  |  |  |  |
SC$H_{0}$ | $0.30_{-0.026}^{+0.025}$ | $4.98_{-1.683}^{+0.715}$ | $71.42_{-0.885}^{+0.901}$ | $0.80_{-0.107}^{+0.114}$ | $-0.82_{-0.077}^{+0.100}$ | $-0.57_{-0.041}^{+0.057}$ | 1780.99 | 1.03 | 1786.99 | 1803.36 | -5.17 | 0.29
SBC$H_{0}$ | $0.33_{-0.010}^{+0.010}$ | $3.38_{-1.363}^{+1.213}$ | $69.55_{-0.671}^{+0.678}$ | $0.67_{-0.033}^{+0.033}$ | $-0.90_{-0.030}^{+0.035}$ | $-0.60_{-0.024}^{+0.026}$ | 1820.05 | 1.03 | 1826.05 | 1842.47 | -2.87 | 2.61
SCH$\textsc{ii}H_{0}$ | $0.28_{-0.024}^{+0.024}$ | $5.12_{-1.303}^{+0.539}$ | $71.14_{-0.779}^{+0.817}$ | $0.85_{-0.108}^{+0.108}$ | $-0.79_{-0.082}^{+0.097}$ | $-0.56_{-0.045}^{+0.059}$ | 2144.92 | 1.12 | 2150.92 | 2167.59 | -7.91 | -2.35
BCH$\textsc{ii}H_{0}$ | $0.31_{-0.012}^{+0.013}$ | $0.57_{-0.402}^{+1.553}$ | $69.78_{-0.579}^{+0.605}$ | $0.67_{-0.029}^{+0.030}$ | $-0.99_{-0.006}^{+0.073}$ | $-0.68_{-0.016}^{+0.055}$ | 420.23 | 1.74 | 426.23 | 436.68 | 1.98 | 5.45
SBCH$\textsc{ii}H_{0}$ | $0.32_{-0.010}^{+0.010}$ | $3.23_{-1.289}^{+1.198}$ | $69.49_{-0.622}^{+0.594}$ | $0.68_{-0.032}^{+0.032}$ | $-0.90_{-0.029}^{+0.034}$ | $-0.61_{-0.023}^{+0.025}$ | 2185.47 | 1.13 | 2191.47 | 2208.19 | -2.61 | 2.96
EXP |  |  |  |  |  |  |  |  |  |  |  |
SC$H_{0}$ | $0.27_{-0.040}^{+0.036}$ | $2.48_{-0.662}^{+0.473}$ | $71.46_{-0.891}^{+0.888}$ | $0.85_{-0.140}^{+0.172}$ | $-0.80_{-0.084}^{+0.085}$ | $-0.56_{-0.040}^{+0.042}$ | 1781.83 | 1.03 | 1787.83 | 1804.2 | -4.32 | 1.13
SBC$H_{0}$ | $0.33_{-0.010}^{+0.010}$ | $1.41_{-0.466}^{+0.381}$ | $69.75_{-0.758}^{+0.732}$ | $0.68_{-0.038}^{+0.034}$ | $-0.90_{-0.042}^{+0.039}$ | $-0.61_{-0.031}^{+0.030}$ | 1820.84 | 1.03 | 1826.84 | 1843.26 | -2.08 | 3.4
SCH$\textsc{ii}H_{0}$ | $0.26_{-0.038}^{+0.036}$ | $2.62_{-0.572}^{+0.392}$ | $71.22_{-0.801}^{+0.807}$ | $0.93_{-0.147}^{+0.189}$ | $-0.77_{-0.084}^{+0.084}$ | $-0.56_{-0.041}^{+0.043}$ | 2146.11 | 1.12 | 2152.11 | 2168.78 | -6.71 | -1.16
BCH$\textsc{ii}H_{0}$ | $0.31_{-0.011}^{+0.011}$ | $0.47_{-0.318}^{+0.412}$ | $69.92_{-0.571}^{+0.562}$ | $0.67_{-0.032}^{+0.033}$ | $-1.00_{-0.005}^{+0.050}$ | $-0.68_{-0.017}^{+0.036}$ | 420.24 | 1.74 | 426.24 | 436.7 | 2.0 | 5.47
SBCH$\textsc{ii}H_{0}$ | $0.32_{-0.010}^{+0.010}$ | $1.36_{-0.449}^{+0.368}$ | $69.68_{-0.662}^{+0.658}$ | $0.69_{-0.036}^{+0.033}$ | $-0.90_{-0.041}^{+0.038}$ | $-0.61_{-0.030}^{+0.029}$ | 2185.98 | 1.13 | 2191.98 | 2208.69 | -2.1 | 3.47
TSUJI |  |  |  |  |  |  |  |  |  |  |  |
SC$H_{0}$ | $0.29_{-0.035}^{+0.032}$ | $4.44_{-1.031}^{+0.569}$ | $71.45_{-0.896}^{+0.892}$ | $0.85_{-0.140}^{+0.172}$ | $-0.80_{-0.084}^{+0.085}$ | $-0.56_{-0.040}^{+0.042}$ | 1781.6 | 1.03 | 1787.6 | 1803.97 | -4.56 | 0.89
SBC$H_{0}$ | $0.33_{-0.010}^{+0.010}$ | $3.13_{-0.873}^{+0.740}$ | $69.82_{-0.723}^{+0.717}$ | $0.68_{-0.038}^{+0.034}$ | $-0.90_{-0.042}^{+0.039}$ | $-0.61_{-0.031}^{+0.030}$ | 1819.49 | 1.03 | 1825.49 | 1841.91 | -3.43 | 2.04
SCH$\textsc{ii}H_{0}$ | $0.27_{-0.034}^{+0.032}$ | $4.58_{-0.856}^{+0.447}$ | $71.19_{-0.811}^{+0.787}$ | $0.93_{-0.147}^{+0.189}$ | $-0.77_{-0.084}^{+0.084}$ | $-0.56_{-0.041}^{+0.043}$ | 2145.82 | 1.12 | 2151.82 | 2168.49 | -7.0 | -1.45
BCH$\textsc{ii}H_{0}$ | $0.31_{-0.011}^{+0.012}$ | $1.06_{-0.731}^{+1.069}$ | $69.91_{-0.570}^{+0.559}$ | $0.67_{-0.032}^{+0.033}$ | $-1.00_{-0.005}^{+0.050}$ | $-0.68_{-0.017}^{+0.036}$ | 420.22 | 1.74 | 426.22 | 436.67 | 1.97 | 5.44
SBCH$\textsc{ii}H_{0}$ | $0.32_{-0.010}^{+0.010}$ | $3.05_{-0.847}^{+0.713}$ | $69.73_{-0.649}^{+0.658}$ | $0.69_{-0.036}^{+0.033}$ | $-0.90_{-0.041}^{+0.038}$ | $-0.61_{-0.030}^{+0.029}$ | 2184.7 | 1.12 | 2190.7 | 2207.41 | -3.38 | 2.19
aTanh |  |  |  |  |  |  |  |  |  |  |  |
SC$H_{0}$ | $0.24_{-0.030}^{+0.035}$ | $1.75_{-0.752}^{+0.627}$ | $71.47_{-0.879}^{+0.868}$ | $0.90_{-0.140}^{+0.155}$ | $-0.77_{-0.066}^{+0.047}$ | $-0.59_{-0.027}^{+0.024}$ | 1782.11 | 1.03 | 1788.11 | 1804.48 | -4.05 | 1.41
SBC$H_{0}$ | $0.32_{-0.010}^{+0.010}$ | $0.16_{-0.104}^{+0.139}$ | $69.65_{-0.781}^{+0.764}$ | $0.61_{-0.026}^{+0.026}$ | $-0.97_{-0.022}^{+0.034}$ | $-0.65_{-0.018}^{+0.023}$ | 1824.32 | 1.04 | 1830.32 | 1846.74 | 1.4 | 6.88
SCH$\textsc{ii}H_{0}$ | $0.23_{-0.027}^{+0.033}$ | $1.86_{-0.705}^{+0.555}$ | $71.17_{-0.770}^{+0.770}$ | $0.96_{-0.144}^{+0.151}$ | $-0.76_{-0.059}^{+0.041}$ | $-0.59_{-0.026}^{+0.023}$ | 2146.48 | 1.12 | 2152.48 | 2169.15 | -6.34 | -0.79
BCH$\textsc{ii}H_{0}$ | $0.31_{-0.011}^{+0.011}$ | $0.05_{-0.034}^{+0.066}$ | $69.76_{-0.587}^{+0.579}$ | $0.66_{-0.029}^{+0.029}$ | $-0.99_{-0.006}^{+0.012}$ | $-0.69_{-0.013}^{+0.015}$ | 420.27 | 1.74 | 426.27 | 436.73 | 2.02 | 5.5
SBCH$\textsc{ii}H_{0}$ | $0.32_{-0.009}^{+0.010}$ | $0.14_{-0.092}^{+0.128}$ | $69.64_{-0.690}^{+0.677}$ | $0.62_{-0.024}^{+0.024}$ | $-0.97_{-0.019}^{+0.031}$ | $-0.66_{-0.016}^{+0.022}$ | 2189.68 | 1.13 | 2195.68 | 2212.4 | 1.6 | 7.17

###  V.1 Constraints on $\Lambda$CDM Model

Figure 4: The $\Lambda$CDM Model(without SH0ES prior for $H_{0}$): The
posterior probability distribution plots of fitted parameters. The color
correspondence for different data-set combinations can be seen in the figure
legends. The darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively. Figure 5: The $\Lambda$CDM
Model(with SH0ES prior for $H_{0}$): The posterior probability distribution
plots of fitted parameters. The color correspondence for different data-set
combinations can be seen in the figure legends. The darker and lighter shades
of colors represent 1-sigma(68.26%) and 2-sigma(95.44%) confidence intervals,
respectively.

We have also included the results for the $\Lambda$CDM model, which is
commonly used as a benchmark for assessing $f(R)$ models. We consider a two-
parameter $\Lambda$CDM model, with the Hubble parameter given by
$H=H_{0}\sqrt{\Omega_{m0}(1+z)^{3}+(1-\Omega_{m0})}$, where the two parameters
$\Omega_{m0}$ and $H_{0}$ stand for the matter density and the Hubble
parameter at the present epoch, respectively. The median values along with
1-sigma confidence intervals for these parameters are presented in Tables 1
and 2. The posterior probability distribution of $\Omega_{m0}$ and $H_{0}$,
both for the respective cases - without and with the SH0ES prior for $H_{0}$
\- are depicted in Figs 4 and 5. For all data sets (except SC), the values of
$\Omega_{m0}$ are compatible with $\Omega_{m0,{\rm Planck}}$ within 1-2(3)
sigma. There is a tension of approximately 1.5-3 sigma between the values of
$H_{0}$ when comparing the cases with and without the SH0ES prior for $H_{0}$.
There exists more or less same order of tensions between the fitted values of
$H_{0}$ and $H_{0,{\rm Planck}}$ or $H_{0,{\rm SH0ES}}$.

###  V.2 Constraints on Hu-Sawicki Model and Starobinsky Model

Figure 6: The Hu-Sawicki Model($n_{{}_{\rm HS}}=1$, without SH0ES prior for
$H_{0}$): The posterior probability distribution plots of fitted parameters.
The color correspondence for different data-set combinations can be seen in
the figure legends. The darker and lighter shades of colors represent
1-sigma(68.26%) and 2-sigma(95.44%) confidence intervals, respectively. Figure
7: The Hu-Sawicki Model($n_{{}_{\rm HS}}=1$, with SH0ES prior for $H_{0}$):
The posterior probability distribution plots of fitted parameters. The color
correspondence for different data-set combinations can be seen in the figure
legends. The darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively. Figure 8: The Starobinsky
Model($n_{{}_{\rm S}}=1$, without SH0ES prior for $H_{0}$): The posterior
probability distribution plots of fitted parameters. The color correspondence
for different data-set combinations can be seen in the figure legends. The
darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively. Figure 9: The Starobinsky
Model($n_{{}_{\rm S}}=1$, with SH0ES prior for $H_{0}$): The posterior
probability distribution plots of fitted parameters. The color correspondence
for different data-set combinations can be seen in the figure legends. The
darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively. Figure 10: The Hu-Sawicki
Model($n_{{}_{\rm HS}}=3$, without SH0ES prior for $H_{0}$): The posterior
probability distribution plots of fitted parameters. The color correspondence
for different data-set combinations can be seen in the figure legends. The
darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively. Figure 11: The Hu-Sawicki
Model($n_{{}_{\rm HS}}=3$, with SH0ES prior for $H_{0}$): The posterior
probability distribution plots of fitted parameters. The color correspondence
for different data-set combinations can be seen in the figure legends. The
darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively.

Previous studies (see Nunes et al. (2017) and references therein) have
established that there exists a degeneracy between $n_{{}_{\rm HS}}/n_{{}_{\rm
S}}$ and $\Omega_{m0}$. One, therefore , often works with fixed values of
$n_{{}_{\rm HS}}$ and $n_{{}_{\rm S}}$ to address this degeneracy. In
Capozziello and Tsujikawa (2008), it has been worked out that $n_{{}_{\rm
HS}},\,n_{{}_{\rm S}}>0.9$, while Fu et al. (2010) proposed that $n_{{}_{\rm
HS}}$ and $n_{{}_{\rm S}}$ should be integers. The difficulty in constraining
higher $n_{{}_{\rm HS}}/n_{{}_{\rm S}}$ values due to numerical instability
issues while solving the modified Hubble equations is also well known.
Consequently, it has become common practice to work with $n_{{}_{\rm HS}}=1$
and $n_{{}_{\rm S}}=1$. In our study we have also explored the case of
$n_{{}_{\rm HS}}=3$. Note that in the re-parameterized version using the
deviation parameter $b$, the Hu-Sawicki model with $n_{{}_{\rm HS}}=2$ is
equivalent to the Starobinsky model with $n_{{}_{\rm S}}=1$.

The quantitative results of fitting the Hu-Sawicki model (Eq. 27) with
$n_{{}_{\rm HS}}=1,3$ and the Starobinsky model (Eq. 29) with $n_{{}_{\rm
S}}=1$, following the procedure discussed in Sec. IV, are presented in the
Tables 1 and 2. The 2D contour plots of the posterior probability distribution
of the model parameters and 1D marginalised distribution of each of the
parameters are shown in the Figs. 6 and 7 (for $n_{{}_{\rm HS}}=1$), Figs. 8
and 9 (for $n_{{}_{\rm HS}}=2$ or, $n_{{}_{\rm S}}=1$), Figs. 10 and 11 (for
$n_{{}_{\rm HS}}=3$), for the respective cases - without and with the SH0ES
prior for $H_{0}$.

The Fig. 1 illustrates that as $n_{{}_{\rm HS}}$ increases (from 1 to 3) the
values of $\Omega_{m0}$ approach closer to those from the $\Lambda$CDM model
(and/or to $\Omega_{m0,\,\rm Planck}$) and the constraints on them become
tighter. This applies to both cases with and without the SH0ES prior for
$H_{0}$. We may observe from Fig. 3 that when the SH0ES prior for $H_{0}$ is
not considered, increasing $n_{{}_{\rm HS}}$ causes the model values of
$H_{0}$ to approach $H_{0,\,{\rm Planck}}$, except for the SC dataset where
there is almost no change with varying $n_{{}_{\rm HS}}$. However, the
constraints from the BCHii dataset are not considered reliable due to high
$\chi^{2}/{\nu}$. When the SH0ES prior is included, there is hardly any change
in the model values of $H_{0}$ with respect to $n_{{}_{\rm HS}}$, and they lie
around $69.5-71.5$ km s-1Mpc-1. As expected, the parameter $b$ increases as
$n_{{}_{\rm HS}}$ increases (see Fig. 2). Also, the constraints on $b$ become
weaker. We may infer that for higher values of $n_{{}_{\rm HS}}$ (say 4, 5, …)
and $n_{{}_{\rm S}}$ (say 2, 3, ..), the constraints on $b$ would become even
more weaker. This implies that with further increase in $n_{{}_{\rm HS}}$ or
$n_{{}_{\rm S}}$ , one can obtain a model that closely resembles the
$\Lambda$CDM model in terms of its predictions for the physical parameters
$\Omega_{m0}$ and $H_{0}$, despite the constraints on $b$ moving towards
larger values. In a sense, this undermines the purpose of exploring $f(R)$
models further, and there is a strategic motivation to limit the exploration
to $n_{{}_{\rm HS}}=1,\,2$ only, putting aside the consideration of
computational difficulties in constraining for the higher values of
$n_{{}_{\rm HS}}$ or $n_{{}_{\rm S}}$. For all three models, there are
instances where $b=0$ is only very marginally allowed(i.e., within a 2-3 sigma
limit), indicating distinguishable models from the $\Lambda$CDM model. This
finding is also an important result of the present work.

###  V.3 Constraints on Exponential Model

Figure 12: The Exponential Model(without SH0ES prior for $H_{0}$): The
posterior probability distribution plots of fitted parameters. The color
correspondence for different data-set combinations can be seen in the figure
legends. The darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively. Figure 13: The Exponential
Model(with SH0ES prior for $H_{0}$): The posterior probability distribution
plots of fitted parameters. The color correspondence for different data-set
combinations can be seen in the figure legends. The darker and lighter shades
of colors represent 1-sigma(68.26%) and 2-sigma(95.44%) confidence intervals,
respectively.

In the Tables 1 and 2, we present the median values and 1-sigma (68.26%)
confidence intervals for the parameters of the exponential model (Eq. 31). The
corresponding 2D contour plots illustrating the posterior probability
distribution of the parameters, along with the 1D marginal distributions for
each parameter, are depicted in Figs 12 and 13 for the cases - without and
with the SH0ES prior for $H_{0}$, respectively. We observe from Figs 12, 13,
and 2 that except for the BCH$\textsc{ii}(H_{0})$ dataset, the parameter $b$
clearly deviates from zero, with $b=0$ barely allowed. Despite this, the
values of the parameters $\Omega_{m0}$ and $H_{0}$ (as illustrated in Figs. 1
and 3) are reasonably close to the standard values derived from Planck
constraints Aghanim et al. (2020). The significance of this result becomes
further clear in the subsequent section on model comparisons. When considering
the SH0ES prior for $H_{0}$, the range of model values of $H_{0}$ is $\sim
69.7$ to $71.5$ km s-1Mpc-1.

###  V.4 Constraints on Tsujikawa Model

Figure 14: The Tsujikawa Model(without SH0ES prior for $H_{0}$): The posterior
probability distribution plots of fitted parameters. The color correspondence
for different data-set combinations can be seen in the figure legends. The
darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively. Figure 15: The Tsujikawa
Model(with SH0ES prior for $H_{0}$): The posterior probability distribution
plots of fitted parameters. The color correspondence for different data-set
combinations can be seen in the figure legends. The darker and lighter shades
of colors represent 1-sigma(68.26%) and 2-sigma(95.44%) confidence intervals,
respectively.

The constraints on the parameters of the Tsujikawa model described by Eq. 33
are presented in the Tables 1 and 2, as well as in Figs. 14 and 15. These
constraints are shown separately for cases without and with the SH0ES prior on
$H_{0}$. Similar to the findings of the exponential model, here also we
observe that the deviation parameter $b$ is mostly nonzero, with only a very
marginal allowance for it to be zero, across all the datasets, except
BCH$\textsc{ii}(H_{0})$. We disregard the results from BCH$\textsc{ii}(H_{0})$
on account of very high value of the reduced $\chi^{2}_{\rm min}/\nu$. Hence,
we obtain an $f(R)$ that is clearly distinguishable from the $\Lambda$CDM
model. The fitted values of $\Omega_{m0}$ and $H_{0}$ also fall within
reasonable limits compared to $\Omega_{m0,{\rm Planck}}$ and $H_{0,{\rm
Planck}}$ (or $H_{0,{\rm SH0ES}}$). With SH0ES prior for $H_{0}$, the model’s
median values for $H_{0}$ range from approximately 69.7 to 71.5 km s-1Mpc-1,
which introduces a tension of 2-3 sigma with $H_{0,{\rm Planck}}$ or
$H_{0,{\rm SH0ES}}$.

###  V.5 Constraints on arcTanh Model

Figure 16: The arcTanh Model(without SH0ES prior for $H_{0}$): The posterior
probability distribution plots of fitted parameters. The color correspondence
for different data-set combinations can be seen in the figure legends. The
darker and lighter shades of colors represent 1-sigma(68.26%) and
2-sigma(95.44%) confidence intervals, respectively. Figure 17: The arcTanh
Model(with SH0ES prior for $H_{0}$): The posterior probability distribution
plots of fitted parameters. The color correspondence for different data-set
combinations can be seen in the figure legends. The darker and lighter shades
of colors represent 1-sigma(68.26%) and 2-sigma(95.44%) confidence intervals,
respectively.

The results obtained by constraining the aTanh model are presented in the
Tables 1 and 2, as well as in the Figs. 16 and 17, separately for cases
without and with the SH0ES prior for $H_{0}$. In several cases, excluding
BCH$\textsc{ii}(H_{0})$, the parameter $b=0$ falls well within the 1-2 sigma
range. However, there are instances where $b=0$ lies outside the 2-sigma
range. In the subsequent section on model comparison, we will observe that
these latter cases are also more favored, making this outcome significant. The
fitted values of $\Omega_{m0}$ for the model fall within the 1-2 sigma range
of $\Omega_{m0,{\rm Planck}}$. Furthermore, $H_{0,{\rm Planck}}$ is
approximately within the 1-2 sigma range of that predicted by the model when
the SH0ES prior is not considered. Upon inclusion of the SH0ES prior for
$H_{0}$, the model’s predicted median values for $H_{0}$ range from
approximately 69.65 to 71.47 km s-1Mpc-1, which are in 2-3 sigma tension with
$H_{0,{\rm Planck}}$ or $H_{0,{\rm SH0ES}}$.

It is important to highlight that the outcomes of this model are not
significantly different from those of the HS1 model. This similarity would not
sound surprising, if we examine Eqs 27 and 34 along with the fact that
$\mathrm{arctanh}(\Lambda/R)\approx{\Lambda/R}$ for $\Lambda/R\ll 1$ (which
holds true in this case since $\Lambda\sim 0.7H_{0}^{2}$ and $R\gtrsim
8H_{0}^{2}$).

##  VI Model Comparison

The standard statistical tools commonly used to assess model fitting (and
model comparison) in cosmology include the reduced chi-square statistics
($\chi^{2}_{\nu}$), the Akaike Information Criterion (AIC) Akaike (1974), and
the Bayesian Information Criterion (BIC) Schwarz (1978) (also see Liddle
(2004, 2007); Burnham and Anderson (2004)). The last six columns of Tables 1
and 2 contain essential quantities related to the current analysis,
specifically pertaining to the utilization of these statistical tools. The
reduced chi-square statistics is defined as $\chi^{2}_{\nu}=\chi^{2}_{\rm
min}/{\nu}$ whereas the AIC and the BIC are respectively defined by the
following equations

$\displaystyle{\rm AIC}$ $\displaystyle=$ $\displaystyle-2\ln\mathcal{L}_{\rm
max}+2k\,,$ (54) $\displaystyle{\rm BIC}$ $\displaystyle=$
$\displaystyle-2\ln\mathcal{L}_{\rm max}+k\ln{N}.$ (55)

The number of degrees of freedom, represented by the symbol $\nu$, is
determined by subtracting the number of model parameters ($k$) from the total
number of data points ($N$). The minimum value of $\chi^{2}$ is denoted as
$\chi^{2}_{\rm min}$ and is connected to the maximum likelihood
($\mathcal{L}_{\rm max}$) through the relation $\chi^{2}_{\rm
min}=-2\ln\mathcal{L}_{\rm max}$. While comparing multiple competing models
using a given data set, the model with the lowest values of $\chi^{2}$, AIC,
and BIC is considered to be more favoured by the data. However, it is
insufficient to rely solely on $\chi^{2}$ due to the principle of Occam’s
razor, which emphasizes the importance of considering the number of model
parameters. Typically, for a nested model, as the number of parameters
increases, the fit improves, leading to a decrease in $\chi^{2}_{\rm min}$ (or
an increase in likelihood), regardless of the relevance of the newly included
parameter(s). Both AIC and BIC incorporate a penalty term ($2k$ and $k\ln{N}$,
respectively) that penalizes models with more parameters, in addition to the
$\chi^{2}_{\rm min}$ term, taking into account any improvement in fit. Thus a
balance between the quality of fit and the complexity of the model is achieved
through AIC and BIC. Comparing the Eqs. 54 and 55 for AIC and BIC,
respectively, it can be observed that the penalty for models with a greater
number of parameters is more severe in BIC than in AIC. Unfortunately, the
conclusions derived from applying these two criteria may sometimes disagree.
In such situations of disagreement, it is necessary to investigate whether any
violations of the assumptions on which AIC and BIC are based have occurred
(for details see Liddle (2004, 2007); Burnham and Anderson (2004)).

To perform a comparative analysis of two models, a useful measure
$\Delta{X}=X_{2}-X_{1}$, computed as the difference between the values of $X$
(= AIC or BIC) for the two models: 1 and 2, facilitates a quantitative
assessment of the evidence supporting the preference of model 1 over model 2.
A rule of thumb which is commonly followed as a guideline to indicate the
degree of strength of the evidence with which the model 2 (model 1) is
favoured over the other, is as follows: (i) $0<\Delta{X}\leqslant 2$
($-2\leqslant\Delta{X}<0$): weak evidence, (ii) $2<\Delta{X}\leqslant 6$
($-6\leqslant\Delta{X}<-2$): positive evidence, (iii) $6<\Delta{X}\leqslant
10$ ($-10\leqslant\Delta{X}<-6$): strong evidence, (iv) $\Delta{X}>10$
($\Delta{X}<-10$): very strong (highly pronounced) evidence.

As can be seen from Tables 1 and 2, the reduced $\chi^{2}$ values for the
BCH$\textsc{ii}(H_{0})$ data-set do not come close to the value 1, for all the
models. Consequently, these cases are not considered in our subsequent
discussions related to model comparisons. The purpose of inclusion of these
cases here is to illustrate that not all possible combinations of data-sets
yield significantly meaningful results.

When we analyze different data-sets without considering the SH0ES prior for
$H_{0}$, we observe interesting patterns in Table 1. The $\Delta$AICs values
indicate varying degrees of evidence ranging from weak to very strong in favor
of $f(R)$ models compared to the $\Lambda$CDM model. On the other hand, the
values of $\Delta$BIC suggest evidence that is either weak or positive against
$f(R)$ models in certain cases, while in other cases, the evidence is weak,
positive, or even strong in favour of $f(R)$ models. The strongest support for
all $f(R)$ models is from the SBC data-set, followed by the SBCHii data-set.
In both cases, the trends of $\Delta$AICs and $\Delta$BICs align, guiding us
in the same direction.

When the SH0ES prior for $H_{0}$ is included, the results, as shown in Table
2, exhibit some differences. The overall support for $f(R)$ models weakens
compared to the previous case. Based on consideration of the $\Delta$AICs, the
maximum support for any $f(R)$ model is now only strong (no longer very
strong), and weak evidence against the $f(R)$ models is observed in a few
cases. On the other hand, analyzing the $\Delta$BICs for data-sets with the
SH0ES prior for $H_{0}$ reveals weak, positive, or strong evidence against the
$f(R)$ models when compared to the $\Lambda$CDM model. Among the data-sets
with the SH0ES prior for $H_{0}$, the strongest evidence against the $f(R)$
models is observed in the SBCH$\textsc{ii}H_{0}$ data-set, followed by the
SBC$H_{0}$ data-set. Furthermore, Table 2 highlights that when a SH0ES prior
for $H_{0}$ is employed, the evidence against models such as HS3, EXP, and
TSUJI (which are more similar to the $\Lambda$CDM model) is weak or mildly
positive according to $\Delta$BIC. Conversely, for models such as HS1, aTanh,
and ST (which are less or least similar to the $\Lambda$CDM model), the
evidence against them is strong when evaluated using $\Delta$BIC.
Consequently, we can conclude that with a SH0ES prior for $H_{0}$, the
$\Delta$BIC criterion disfavours models that significantly deviate from the
$\Lambda$CDM model.

A crucial observation to note is that the cases of $f(R)$ models with strong
or very strong evidence for being preferred over the $\Lambda$CDM model (as
indicated by $\Delta$AIC and/or $\Delta$BIC) also correspond to the cases
where a non-zero value of the deviation parameter $b$ is favored, with the
possibility of a very marginal acceptance of $b=0$. This finding is an
important result of this study, which can be observed from Figs. 6-17 and
Tables 1 and 2.

##  VII Accelerated Expansion of the Universe

In this section we discuss the findings regarding the relevant quantities that
describe the current accelerated expansion of the universe. The behavior of
the deceleration parameter ($q(z)$), which is defined as
$q(z)\equiv-a\ddot{a}/\dot{a}^{2}=-\ddot{a}/(H^{2}a)$, provides valuable
insight into whether the cosmic expansion is accelerating ($q<0$) or
decelerating ($q>0$), serving as an indicator of the respective phases.
Through several model-independent measurements of $q(z)$, it has been observed
that $q(z=0)<0$ and $q(z>z_{\rm t})>0$ Rani et al. (2015); Jesus et al.
(2020). The quantity $z_{\rm t}$, called transition redshift, indicates the
redshift at which the universe has undergone transition from a decelerated
phase to an accelerated phase of expansion.

The total equation-of-state (EoS) parameter ($w_{\rm eff}$) as a function of
redshift provides valuable information about the evolution of the universe
through the eras of radiation domination, matter domination, and the
subsequent late time accelerated phase of expansion dubbed as dark energy era.
During these three phases, $w_{\rm eff}$ takes values approximately equal to
$1/3$, $0$, and $-1$ (for $\Lambda$CDM). On the other hand, the EoS parameter
($w_{\rm DE}$) associated with the dark energy component remains constant at
$-1$ for all values of redshift (both in the past and future) in $\Lambda$CDM
model. However, in the case of any viable $f(R)$ model, value of the parameter
$w_{\rm DE}$ can cross the so called “phantom-divide” line ($w_{\rm
DE}(z)=-1$) multiple times before finally settling at $-1$ in the far future
Jaime et al. (2012); Bamba et al. (2010).

In the Tables 1 and 2 we present the median values and 1-sigma confidence
intervals for the quantities $z_{\rm t}$, $w_{\rm DE,0}$, and $w_{\rm eff,0}$
for two cases - without and with the SH0ES prior on $H_{0}$. In Fig 18, we
display the values of $z_{\rm t}$ for all models and data-sets. This figure
shows the median values of $z_{\rm t}\sim 0.5-1$ for all the models, which is
compatible with the model-independent predictions for $z_{\rm t}$ Rani et al.
(2015); Jesus et al. (2020). In the Figs. 19 and 20 we have plotted the
evolution of $w_{\rm DE}$ with redshift based on the best constraints obtained
for the models, (which happens to be the constraints from the data-set
SCH$\textsc{ii}(H_{0})$). In all $f(R)$ models, the values of $w_{\rm DE}(z)$
mark a crossover from the so-called “phantom regime” in the distant past to
the so-called “quintessence regime” in the more recent past, with the present-
epoch values ($w_{\rm DE,0}$) falling within the quintessence region. When the
SH0ES prior on $H_{0}$ is included, the deviations of $w_{\rm DE}(z)$ for
$f(R)$ models from the $\Lambda$CDM model are reduced and the transition from
the phantom regime to the quintessence regime happens more recently compared
to the situation without the SH0ES prior.

For any viable $f(R)$ model, the $w_{\rm DE}(z)$ curve is expected to
eventually reach a stable, constant value of $-1$ in the distant past \- a
behavior which is necessary to account for the matter-dominated era in the
context of the model. From Figs. 19 and 20, we can readily observe that the
Tsujikawa model and the Exponential model exhibit a clear convergence to $-1$,
while for the other models this convergence takes place in a more remote past,
characterized by higher values of $z$. The trend for possible future crossing
into phantom regime is clearly visible for the HS1 model, the Starobinsky
model and the aTanh model, whereas this trend is not as apparent in the other
models. However, theoretical studies Jaime et al. (2012); Bamba et al. (2010)
have indicated that these models are also expected to cross the phantom divide
in the future, albeit in an oscillatory manner. Although it would be an
intriguing future project to extend the analysis of the $w_{\rm DE}(z)$ curve
up to $z\rightarrow-1(a\rightarrow\infty)$, it is beyond the scope of the
present work. This is due to the potential limitations of the employed ODE
system, which may not accurately capture the oscillatory features in such a
far future regime.

The behavior of $w_{\rm eff}(z)$ is depicted in Figs. 21 and 22, revealing a
transition from a matter-dominated era (represented by $w_{\rm
eff}(z)\rightarrow 0$) to a dark energy dominated era around $z\sim 5$ for all
the $f(R)$ models. Considering the SH0ES prior for $H_{0}$ leads to relatively
lower values of $w_{\rm eff,0}$ for all the $f(R)$ models compared to the
cases without the prior. Furthermore, the discrepancy between the predictions
of $w_{\rm eff}(z)$ obtained from the $\Lambda$CDM model and those from the
$f(R)$ models reduces when the SH0ES prior is included.

Figure 18: Transition Redshift: The blank star and blank circle markers
represent median values for the cases with and without SH0ES prior for
$H_{0}$, respectively. Different colors represent different data-set
combinations and these are same as in any of the parameter distribution
plots(e.g. Fig. 6) or see 2nd paragraph of Sec.V. The thick and thin
horizontal colored bars represent 1-sigma(68.26%) and 2-sigma(95.44%)
confidence intervals for the cases with SH0ES prior. The colored
continuous/dashed lines represent 1-sigma(68.26%, with smaller cap size) and
2-sigma(95.44%, with bigger cap size) confidence intervals for the cases
without SH0ES prior. Figure 19: Evolution of EoS parameter of the “effective
geometric dark energy” for all $f(R)$ models obtained from the data-set
SBCHii. Different colored solid lines show median values for different models
as indicated by legends whereas with same color shaded areas show 1-sigma
confidence interval. Figure 20: Evolution of EoS parameter of the “effective
geometric dark energy” for all $f(R)$ models obtained from the data-set
SBCH$\textsc{ii}H_{0}$. Different colored solid lines show median values for
different models as indicated by legends whereas shaded areas with colors show
1-sigma confidence interval. Figure 21: Evolution of the total EoS parameter
for all $f(R)$ models and the $\Lambda$CDM model obtained from the data-set
SBCHii. Different colored solid lines show median values for different models
as indicated by legends whereas dashed lines(or shaded area) with colors show
1-sigma confidence interval. For better contrast we choose $\log_{10}(z)$ as
independent variable. Figure 22: Evolution of the total EoS parameter for all
$f(R)$ models and the $\Lambda$CDM model obtained from the data-set
SBCH$\textsc{ii}H_{0}$. Different colored solid lines show median values for
different models as indicated by legends whereas dashed lines(or shaded area)
with colors show 1-sigma confidence interval. For better contrast we choose
$\log_{10}(z)$ as independent variable.

##  VIII Conclusions

In this study, we used the recently released Pantheonplus SNIa data along with
cosmic chronometer measurements, baryonic acoustic oscillation measurements,
Hii starburst galaxies data, and the local measurement of the Hubble
parameter, to constrain popularly studied $f(R)$ models in the metric
formalism. The redshift coverage of these data-sets is $0.0012\lesssim
z\lesssim 2.55$, and the sparse region beyond $z>1.4$ is complemented by the
data from Hii galaxies. The inclusion of Hii galaxies data has shown
improvements in parameter constraints, even though its use in constraining
$f(R)$ models is not yet widespread. In light of the so-called “Hubble
tension”, we have also examined cases with vis-a-vis without a SH0ES prior for
$H_{0}$ — a feature which is very often lacking in many earlier studies (see
references in Sec I). It is noteworthy to find that incorporating this prior
for $H_{0}$, based on local measurements, generally leads to reduction in the
differences between predictions of $f(R)$ models and the $\Lambda$CDM model

The six viable $f(R)$ models examined in this study are the Hu-Sawicki model
with $n_{{\rm HS}}=1,3$, the Starobinsky model ($n_{{\rm S}}=1)$, the
exponential, the Tsujikawa and the arcTanh model. Reparametrisation of these
$f(R)$ models using the deviation parameter $b$, shows that in the limit $b\to
0$, these models converge to the standard $\Lambda$CDM model. Each of these
$f(R)$ models is characterized by three parameters: the deviation parameter
$b$, the present-epoch values of the matter density parameter $\Omega_{m0}$,
and the Hubble parameter $H_{0}$.

Regardless of whether the SH0ES prior for $H_{0}$ is considered, the
$\Omega_{m0}$ values for all the examined models fall within the 1-2 sigma
range of the Planck constrained value ($\Omega_{m0,\rm Planck}=0.315\pm
0.007$) or, conversely, the Planck value lies within the 1-2 sigma range of
the model values. When no SH0ES prior for $H_{0}$ is considered, the fitted
values of $H_{0}$ for the models are also within the 1-2 sigma range of that
from the Planck constraint ($H_{0,\rm Planck}=67.4\pm 0.5$), or alternatively,
the Planck value lies within the 1-2 sigma range of the fitted model values.
When incorporating the SH0ES prior on $H_{0}$, the fitted model values of
$H_{0}$ tend to lie within the 2-3 sigma range on the lower side of the SH0ES
measurement of $H_{0,\rm SH0ES}=73.04\pm 1.04$. An additional impact resulting
from the SH0ES prior on $H_{0}$ is the tightening of constraints on all three
model parameters for all the models and data-sets.

In our study, we have found instances of all the investigated $f(R)$ models
where the deviation parameter $b$ significantly deviates from zero, indicating
that $b=0$ is only marginally allowed. It is worth highlighting that these
cases also correspond to the situations where $\Delta$AIC and/or $\Delta$BIC
based analysis suggest that the $f(R)$ models are (very) strongly favoured
over the $\Lambda$CDM model. Based on our knowledge of previous studies
(referred in Sec I), it can be stated that up until now, the support for
$f(R)$ models has mainly been weak or positive (based on $\Delta$AIC and/or
$\Delta$BIC). However, our current research yields (very) strong support for
$f(R)$ models, along with support for non-zero values of the deviation
parameter $b$. This indicates that the viability of $f(R)$ models have not yet
been ruled out by the available cosmological data sets.

We find that the relevant quantities characterizing the (accelerated)
expansion of the universe — $z_{\rm t}$, $w_{\rm eff,0}$, and $w_{\rm DE,0}$ —
estimated in this study are compatible with their model-independent estimates
from earlier works. All the models examined in this study, predict that the
universe underwent a transition from a decelerated phase of expansion to an
accelerated phase of expansion in the recent past, occurring at $z_{\rm t}\sim
0.5-1$. Furthermore, the current values of $w_{\rm DE}$ fall within the
quintessential region, having recently crossed over from the phantom region
(around $z\sim 0.5-2$, depending on the specific models).

The instances of strong evidence in favor of $f(R)$ models accompanied by
clear non-zero $b$(or only very marginal allowance for $b=0$), agreements
between the derived quantities $z_{\rm t}$, $w_{\rm eff,0}$, and $w_{\rm
DE,0}$ here and their expected values from model-independent predictions,
suggest that the analyzed cosmological data sets allow room for the
consideration of viability of $f(R)$ models as an explanation for the observed
late time cosmic acceleration. In fact, it would be worthwhile to conduct even
further investigations of these models using future data sets.

## Data Availability

The observed cosmological data such as SNIa, CC, BAO, Hii starburst galaxies
and local measurement of $H_{0}$ are publicly available — the references to
which are cited in the text. The simulated data generated in this work are
available from the corresponding author, KR, upon reasonable request.

## Acknowledgements

K.R. would like to thank HoD, Dept. of Comp. Sc., RKMVERI, for providing
computational facilities. A.C. would like to thank Indian Institute of
Technology, Kanpur, for supporting this work by means of Institute Post-
Doctoral Fellowship (Ref. No. DF/PDF197/2020-IITK/970).

##  Appendix A The Modified Friedmann Equations as a System of First Order
ODEs

Figure 23: The color projection plot of relation Eq. 67 which show
$z^{i}(\Omega_{\rm m0}^{\Lambda},\,b)$ below which the solutions need to be
switched from the $\Lambda$CDM to $f(R)$(at $z^{i}$ we obtain the initial
conditions for the ODE system Eqs. 60 to 65 from $\Lambda$CDM approximations).
Top panel is plot is for the Tsujikawa model and the bottom panel plot is for
the Hu-Sawicki model($n_{{}_{\rm HS}}=1$).

The strategies employed to solve the system of Eqs. 7, 10 and 11 have been
proposed and utilized in various earlier studiesAmendola et al. (2007);
Carloni et al. (2009); Abdelwahab et al. (2012); de la Cruz-Dombriz et al.
(2016); Odintsov and Oikonomou (2017); Odintsov et al. (2017); Farrugia et al.
(2021). Although there has been proposals for using model-dependent schemes
Odintsov et al. (2017); Leizerovich et al. (2022) too, in our current work, we
were able to obtain stable solutions for all the considered $f(R)$ models
using a single scheme. This scheme involves transforming the ODEs into a
system of first-order ODEs using dynamical variables, which aviods numerical
instability and reduces overall computational cost. Based on the literature
reviewed, we define the following dimensionless variables

$\begin{split}&X_{1}=\frac{R}{6(\eta
H_{0}^{\Lambda})^{2}},\,\,X_{2}=\frac{\dot{R}F^{\prime}}{(\eta
H_{0}^{\Lambda})F},\,\,X_{3}=\frac{f}{6(\eta H_{0}^{\Lambda})^{2}F}\,,\\\
&O_{\rm m}=\frac{\Omega_{\rm m0}^{\Lambda}(1+z)^{3}}{\eta^{2}F},\,\,O_{\rm
r}=\frac{\Omega_{\rm r,0}^{\Lambda}(1+z)^{4}}{\eta^{2}F}\,,\\\ &H=\eta
H_{0}^{\Lambda},\,\,\Gamma=\frac{F}{RF^{\prime}},\,\,r=\frac{R}{R_{\star}}\,,\end{split}$
(56)

where $z$ denotes redshift, $\Omega_{\rm m0}^{\Lambda}$ and $\Omega_{\rm
r,0}^{\Lambda}$ respectively represent the matter and radiation density at
present epoch, and $\eta=H/H_{0}^{\Lambda}$. The superscript $\Lambda$ denotes
quantities inferred in a $\Lambda$CDM paradigm. The parameter $R_{\star}$ has
the dimension of Ricci scalar, and it’s value gets fixed by the choice of a
specific $f(R)$ model (e.g. for the Hu-Sawicki model $R_{\star}\equiv\mu^{2}$,
for the Starobinsky model $R_{\star}\equiv{R}_{\rm S}$, for the exponential
model $R_{\star}\equiv 1/\beta$ and for the Tsujikawa model $R_{\star}\equiv
R_{\rm T}$). The fact that any $f(R)$ model is expected to approach
$\Lambda$CDM cosmology at sufficiently high redshifts (say $z^{i}$, initial
redshift) allows us to set the necessary initial conditions for the ODE
system. The determination of $z^{i}$ will be discussed later in this section.

In terms of the above dynamical/dimensionless variables defined in Eq. 56, the
Eqs. 7 and 10 can be written as

$\frac{d\eta}{dz}=\frac{\eta}{(1+z)}(2-X_{1})\,,$ (57)

and,

$1=O_{\rm m}+O_{\rm r}+X_{1}-X_{2}-X_{3}\,,$ (58)

respectively. Thus it becomes evident from the expression in Eq. 58 that it
serves the purpose of setting the initial conditions and monitoring the
solutions at each steps of integration. We can also express the equation for
effective geometric dark energy(Eq. 15) in terms of above variables as

$w_{\rm DE}=\frac{w_{\rm tot}-O_{\rm r}F/3}{1-(O_{\rm m}+O_{\rm r})F}\,.$ (59)

We further define $N=-\log(1+z)$, which allows us to rewrite the system of
first-order ODEs that need to be solved as follows:

$\displaystyle\frac{dX_{1}}{dN}$ $\displaystyle=$ $\displaystyle
X_{1}\left(X_{2}\Gamma-2X_{1}+4\right)\,,$ (60)
$\displaystyle\frac{dX_{2}}{dN}$ $\displaystyle=$ $\displaystyle(O_{\rm
m}+2X_{1}-X_{2}-4X_{3}-X_{1}X_{2}-X_{2}^{2}),$ (61)
$\displaystyle\frac{dX_{3}}{dN}$ $\displaystyle=$
$\displaystyle(X_{1}X_{2}\Gamma-X_{2}X_{3}+4X_{3}-2X_{1}X_{3})\,,$ (62)
$\displaystyle\frac{dO_{\rm m}}{dN}$ $\displaystyle=$ $\displaystyle-O_{\rm
m}(-1+2X_{1}+X_{2})\,,$ (63) $\displaystyle\frac{dr}{dN}$ $\displaystyle=$
$\displaystyle X_{2}\Gamma{r}\,,$ (64) $\displaystyle{\rm and},$
$\displaystyle\frac{dO_{\rm r}}{dN}$ $\displaystyle=$ $\displaystyle-O_{\rm
r}(2X_{1}+X_{2})\,,$ (65)

where $\Gamma=F/(RF^{\prime})$. We choose $\Omega_{\rm
r,0}^{\Lambda}=2.9656\times 10^{-4}\Omega_{\rm m0}^{\Lambda}$, and with this
choice, Eq. 65 becomes redundant as it reduces to Eq. 63. Therefore, we only
need to solve the system of ODEs consisting of Eqs. 60 to 64 and not Eq. 65.
The reason for considering $\Omega_{\rm r,0}^{\Lambda}$ (and hence $O_{r}$)
will be explained later in this section. Ultimately, by solving the above
system, we obtain $H=\sqrt{rR_{\star}/(6X_{1})}$. Using this expression for
$H(z)$, we can determine other cosmological observables as defined in Sec. IV.
While we perform model-fitting with parameters $\Omega_{\rm m0}^{\Lambda}$,
$b$ and $H^{\Lambda}_{0}$, we eventually obtain the relevant parameters for
the $f(R)$ models, viz. $\Omega_{m0}^{f(R)}$, $b$, and $H_{0}^{f(R)}$ (using
Eq. 24), where $H_{0}^{f(R)}$ is the numerical solution of the aforementioned
$H$ at $z=0$.

The initial conditions needed to solve the system of Eqs. 60 to 64 are given
by

$\begin{split}X_{1}^{i}&=\frac{\Omega_{\rm
m0}^{\Lambda}(1+z^{i})^{3}+4\Omega_{\Lambda,0}^{\Lambda}}{2\eta^{i\,2}}\,,\\\
X_{2}^{i}&=0\,,\\\ X_{3}^{i}&=\frac{\Omega_{\rm
m0}^{\Lambda}(1+z^{i})^{3}+2\Omega_{\Lambda,0}^{\Lambda}}{2\eta^{i\,2}}\,,\\\
O_{\rm m}^{i}&=1-O_{\rm r}^{i}-X_{1}^{i}+X_{2}^{i}+X_{3}^{i}\,,\\\
r^{i}&=\frac{3(H_{0}^{\Lambda})^{2}}{R_{\star}}\left[(\Omega_{\rm
m0}^{\Lambda}(1+z^{i})^{3}+4\Omega_{\Lambda,0}^{\Lambda}\right]\,,\end{split}$
(66)

where $\eta^{i\,2}=\Omega_{\rm m0}^{\Lambda}(1+z^{i})^{3}+\Omega_{\rm
r,0}^{\Lambda}(1+z^{i})^{4}+\Omega_{\Lambda,0}^{\Lambda}$, $O_{\rm
r}^{i}=\frac{\Omega_{\rm r,0}^{\Lambda}(1+z^{i})^{4}}{\eta^{i\,2}}$ and
$\Omega_{\Lambda,0}^{\Lambda}=(1-\Omega_{\rm m0}^{\Lambda}-\Omega_{\rm
r,0}^{\Lambda})$.

In order to estimate $z_{i}$ we set $f(R(z_{i}))=R-2\Lambda(1-\epsilon)$ where
$\epsilon\ll 1$. Here the expressions for $R$ and $\Lambda$ correspond to
those in the $\Lambda$CDM cosmology. With this one obtains the following
expression for $z_{i}$ (see Odintsov et al. (2017); Farrugia et al. (2021) for
derivation):

$z_{i}=\left[\frac{4\Omega_{\Lambda,0}^{\Lambda}}{\Omega_{\rm
m0}^{\Lambda}}\left(\frac{b}{4\nu_{\rm f}}-1\right)\right]^{1/3}-1\,,$ (67)

where $\nu_{\rm f}$ for different models (f) are given by following
expressions:

$\begin{split}&\nu_{{}_{\rm HS}}=\frac{1}{(1/\epsilon-1)^{1/n_{{}_{\rm
HS}}}},\,\,\,\nu_{{}_{\rm ST}}=\sqrt{(1/\epsilon)^{1/n_{{}_{\rm S}}}-1}\,,\\\
&\nu_{{}_{\rm E}}=-1/\ln\epsilon\,,\,\,\,\nu_{{}_{\rm
T}}=1/\mathrm{arctanh}(1-\epsilon)\,,\,\,\,\nu_{{}_{\rm
aTanh}}\approx\epsilon\,.\end{split}$ (68)

If, during MCMC sampling, the value of $z^{i}$ for a particular set of
parameters ($\Omega_{\rm m0}^{\Lambda},\,b$) is found to be smaller than the
maximum redshift ($z_{\rm max}$) of the data, we switch from the $\Lambda$CDM
solution within the range $[z^{i},z_{\rm max}]$ to the $f(R)$ solution for the
range $[0,z^{i}]$ at $z^{i}$. Depending on $f(R)$ models, for some samples of
($\Omega_{\rm m0}^{\Lambda},\,b$) Eq. 67 may yield $z^{i}\leq 0$. This
indicates that for such situations the $\Lambda$CDM cosmology is the solution
for all $z>0$(for example, see the top panel plot in Fig. 23).

Prior to running the MCMC code for a large number of samples, we conducted a
preliminary run with a smaller sample size (approximately 50,000) to determine
suitable values for $\epsilon$. This avoids long computation times without
comporimsing with accuracy of the results. Based on the findings of this pilot
project, we choose $\epsilon=10^{-10}$ for the exponential and the Tsujikawa
models, $\epsilon=10^{-8}$ for the Starobinsky model and the Hu-Sawicki model
with $n_{{}_{\rm HS}}=3$, and $\epsilon=10^{-6}$ for the Hu-Sawicki model with
$n_{{}_{\rm HS}}=1$ and the aTanh model.

We have included a color projection plot in Fig. 23 illustrating the
relationship described in Eq . 67 for the Tsujikawa model and the Hu-Sawicki
model ($n_{{}_{\rm HS}}=1$). From this plot, it can be observed that for some
samples of ($\Omega_{\rm m0}^{\Lambda},\,b$), the value of $z^{i}$ can reach
as high as $~{}200$. In such cases, when we have high $z^{i}$ one must include
$\Omega_{\rm r,0}^{\Lambda}$ in the definition of $\eta^{i}$ even though
$\Omega_{\rm r,0}^{\Lambda}$ is not a free parameter and one has data-set
comprising lower redshifts (as in our case). Through our preliminary
investigations, we observed that excluding $\Omega_{\rm r,0}^{\Lambda}$ led to
numerical instability and inaccuracies when solving the system of ODEs
described by Eqs. 60 to 65. Therefore, based on both theoretical principles
and our practical experience, we made the decision to include $\Omega_{\rm
r,0}^{\Lambda}$ in our analysis.

##  Appendix B BAO and CC Data

The data-sets for BAO and CC used in this work is given in the Tables 3 and 4,
respectively. The covariance matrices, taken from respective cited references,
for BAO data points 21-25, 27-28 and 29-30, respectively, are given in the
following equations:

$C_{1}=\begin{bmatrix}624.707&23.729&325.332&8.34963&157.386&3.57778\\\
23.729&5.60873&11.6429&2.33996&6.39263&0.968056\\\
325.332&11.6429&905.777&29.3392&515.271&14.1013\\\
8.34963&2.33996&29.3392&5.42327&16.1422&2.85334\\\
157.386&6.39263&515.271&16.1422&1375.12&40.4327\\\
3.57778&0.968056&14.1013&2.85334&40.4327&6.25936\\\ \end{bmatrix},$ (69)
$C_{2}=\begin{bmatrix}0.0911&-0.0338\\\ -0.0338&0.22\\\ \end{bmatrix},$ (70)

and,

$C_{3}=\begin{bmatrix}0.3047&0.1707\\\ 0.1707&0.6373\\\ \end{bmatrix}.$ (71)

Table 3: The BAO data collected from many cited sources. The three needed covariance matrices for data-points (i)21-25 (ii)27-28 and (iii)29-30 can be found in the respective references and are given here in Eqs. 69, 70, and 71, respectively. No. | $z_{\rm eff}$ | Value | Observable | Reference
---|---|---|---|---
1 | $0.106$ | $0.336\pm 0.015$ | $r_{d}/D_{V}$ | Beutler et al. (2011)
2 | $0.275$ | $0.1390\pm 0.0037$ | $r_{d}/D_{V}$ | Percival et al. (2010)
3 | 0.122 | $3.65\pm 0.115$ | $D_{V}/r_{d}$ | Carter et al. (2018)
4 | 0.15 | $4.47\pm 0.168$ | $D_{V}/r_{d}$ | Ross et al. (2015)
5 | 0.32 | $8.47\pm 0.167$ | $D_{V}/r_{d}$ | Tojeiro et al. (2014)
6 | 0.35 | $8.88\pm 0.17$ | $D_{V}/r_{d}$ | Padmanabhan et al. (2012)
7 | 0.44 | $11.55\pm 0.559$ | $D_{V}/r_{d}$ | Kazin et al. (2014)
8 | $0.57$ | $13.67\pm 0.22$ | $D_{V}/r_{d}$ | Anderson et al. (2013)
9 | 0.6 | $14.95\pm 0.68$ | $D_{V}/r_{d}$ | Kazin et al. (2014)
10 | 0.72 | $16.08\pm 0.406$ | $D_{V}/r_{d}$ | Bautista et al. (2018)
11 | 0.73 | $16.93\pm 0.579$ | $D_{V}/r_{d}$ | Kazin et al. (2014)
12 | 1.52 | $26.0\pm 0.995$ | $D_{V}/r_{d}$ | Ata et al. (2018)
13 | $0.54$ | $9.212\pm 0.41$ | $D_{A}/r_{d}$ | Seo et al. (2012)
14 | 0.697 | $10.39\pm 0.496$ | $D_{A}/r_{d}$ | Sridhar et al. (2020)
15 | $0.81$ | $10.75\pm 0.43$ | $D_{A}/r_{d}$ | Abbott et al. (2019)
16 | 0.874 | $11.37\pm 0.693$ | $D_{A}/r_{d}$ | Sridhar et al. (2020)
17 | 2.3 | $9.07\pm 0.31$ | $D_{H}/r_{d}$ | Bautista et al. (2017)
18 | $2.34$ | $8.86\pm 0.29$ | $D_{H}/r_{d}$ | de Sainte Agathe et al. (2019)
19 | 2.36 | 9.0$\pm 0.3$ | $c/(r_{d}H)$ | Font-Ribera et al. (2014)
20 | $2.4$ | $36.6\pm 1.2$ | $D_{M}/r_{d}$ | du Mas des Bourboux et al. (2017)
21 | $0.38$ | 1512.39 | $D_{M}\left(r_{d,{\rm fid}}/r_{d}\right)$ | Alam et al. (2017)
22 | $0.38$ | 81.2087 | $H(z)\left(r_{d}/r_{d,{\rm fid}}\right)$ | Alam et al. (2017)
23 | $0.51$ | 1975.22 | $D_{M}\left(r_{d,{\rm fid}}/r_{d}\right)$ | Alam et al. (2017)
24 | $0.51$ | 90.9029 | $H(z)\left(r_{d}/r_{d,{\rm fid}}\right)$ | Alam et al. (2017)
25 | $0.61$ | 2306.68 | $D_{M}\left(r_{d,{\rm fid}}/r_{d}\right)$ | Alam et al. (2017)
26 | $0.61$ | 98.9647 | $H(z)\left(r_{d}/r_{d,{\rm fid}}\right)$ | Alam et al. (2017)
27 | $0.698$ | $17.65\pm 0.3$ | $D_{M}/r_{d}$ | Bautista et al. (2020)
28 | $0.698$ | $19.77\pm 0.47$ | $D_{H}/r_{d}$ | Bautista et al. (2020)
29 | $1.48$ | $13.23\pm 0.47$ | $D_{H}/r_{d}$ | Neveux et al. (2020)
30 | $1.48$ | $30.21\pm 0.79$ | $D_{M}/r_{d}$ | Neveux et al. (2020)
Table 4: The CC data collected from many cited sources. $z$ | $H(z)$ | $\sigma_{H(z)}$ | Reference
---|---|---|---
0.07 | 69.0 | 19.6 | Zhang et al. (2014)
0.09 | 69 | 12 | Simon et al. (2005)
0.12 | 68.6 | 26.2 | Zhang et al. (2014)
0.17 | 83 | 8 | Simon et al. (2005)
0.179 | 75 | 4 | Moresco et al. (2012)
0.199 | 75 | 5 | Moresco et al. (2012)
0.20 | 72.9 | 29.6 | Zhang et al. (2014)
0.27 | 77 | 14 | Simon et al. (2005)
0.28 | 88.8 | 36.6 | Zhang et al. (2014)
0.352 | 83 | 14 | Moresco et al. (2012)
0.38 | 83 | 13.5 | Moresco et al. (2016)
0.4 | 95 | 17 | Simon et al. (2005)
0.4004 | 77 | 10.2 | Moresco et al. (2016)
0.425 | 87.1 | 11.2 | Moresco et al. (2016)
0.445 | 92.8 | 12.9 | Moresco et al. (2016)
0.47 | 89.0 | 49.6 | Ratsimbazafy et al. (2017)
0.4783 | 80.9 | 9 | Moresco et al. (2016)
0.48 | 97 | 62 | Stern et al. (2010)
0.593 | 104 | 13 | Moresco et al. (2012)
0.68 | 92 | 8 | Moresco et al. (2012)
0.75 | 98.8 | 33.6 | Borghi et al. (2022)
0.781 | 105 | 12 | Moresco et al. (2012)
0.875 | 125 | 17 | Moresco et al. (2012)
0.88 | 90 | 40 | Stern et al. (2010)
0.9 | 117 | 23 | Simon et al. (2005)
1.037 | 154 | 20 | Moresco et al. (2012)
1.3 | 168 | 17 | Simon et al. (2005)
1.363 | 160 | 33.6 | Moresco (2015)
1.43 | 177 | 18 | Simon et al. (2005)
1.53 | 140 | 14 | Simon et al. (2005)
1.75 | 202 | 40 | Simon et al. (2005)
1.965 | 186.5 | 50.4 | Moresco (2015)

## References

  * Perlmutter et al. [1998] S. Perlmutter et al.  Discovery of a supernova explosion at half the age of the Universe and its cosmological implications.  _Nature_ , 391:51–54, 1998.  doi: 10.1038/34124.
  * Perlmutter et al. [1999] S. Perlmutter et al.  Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae.  _Astrophys. J._ , 517:565–586, 1999.  doi: 10.1086/307221.
  * Schmidt et al. [1998] Brian P. Schmidt et al.  The High Z supernova search: Measuring cosmic deceleration and global curvature of the universe using type Ia supernovae.  _Astrophys. J._ , 507:46–63, 1998.  doi: 10.1086/306308.
  * Riess et al. [1998] Adam G. Riess et al.  Observational evidence from supernovae for an accelerating universe and a cosmological constant.  _Astron. J._ , 116:1009–1038, 1998.  doi: 10.1086/300499.
  * Eisenstein et al. [2005] Daniel J. Eisenstein et al.  Detection of the Baryon Acoustic Peak in the Large-Scale Correlation Function of SDSS Luminous Red Galaxies.  _Astrophys. J._ , 633:560–574, 2005.  doi: 10.1086/466512.
  * Percival et al. [2007a] Will J. Percival et al.  Measuring the matter density using baryon oscillations in the SDSS.  _Astrophys. J._ , 657:51–55, 2007a.  doi: 10.1086/510772.
  * Spergel et al. [2007] D. N. Spergel et al.  Wilkinson Microwave Anisotropy Probe (WMAP) three year results: implications for cosmology.  _Astrophys. J. Suppl._ , 170:377, 2007.  doi: 10.1086/513700.
  * Hinshaw et al. [2009] G. Hinshaw et al.  Five-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Data Processing, Sky Maps, and Basic Results.  _Astrophys. J. Suppl._ , 180:225–245, 2009.  doi: 10.1088/0067-0049/180/2/225.
  * Ade et al. [2014] P. A. R. Ade et al.  Planck 2013 results. XVI. Cosmological parameters.  _Astron. Astrophys._ , 571:A16, 2014.  doi: 10.1051/0004-6361/201321591.
  * Ade et al. [2016] P. A. R. Ade et al.  Planck 2015 results. XIII. Cosmological parameters.  _Astron. Astrophys._ , 594:A13, 2016.  doi: 10.1051/0004-6361/201525830.
  * Cole et al. [2005] Shaun Cole et al.  The 2dF Galaxy Redshift Survey: Power-spectrum analysis of the final dataset and cosmological implications.  _Mon. Not. Roy. Astron. Soc._ , 362:505–534, 2005.  doi: 10.1111/j.1365-2966.2005.09318.x.
  * Percival et al. [2007b] Will J. Percival et al.  The shape of the SDSS DR5 galaxy power spectrum.  _Astrophys. J._ , 657:645–663, 2007b.  doi: 10.1086/510615.
  * Zlatev et al. [1999] Ivaylo Zlatev, Li-Min Wang, and Paul J. Steinhardt.  Quintessence, cosmic coincidence, and the cosmological constant.  _Phys. Rev. Lett._ , 82:896–899, 1999.  doi: 10.1103/PhysRevLett.82.896.
  * Martin [2012] Jerome Martin.  Everything You Always Wanted To Know About The Cosmological Constant Problem (But Were Afraid To Ask).  _Comptes Rendus Physique_ , 13:566–665, 2012.  doi: 10.1016/j.crhy.2012.04.008.
  * Tsujikawa [2013] Shinji Tsujikawa.  Quintessence: A Review.  _Class. Quant. Grav._ , 30:214003, 2013.  doi: 10.1088/0264-9381/30/21/214003.
  * Armendariz-Picon et al. [2001] C. Armendariz-Picon, Viatcheslav F. Mukhanov, and Paul J. Steinhardt.  Essentials of k essence.  _Phys. Rev. D_ , 63:103510, 2001.  doi: 10.1103/PhysRevD.63.103510.
  * Basilakos et al. [2013] Spyros Basilakos, Savvas Nesseris, and Leandros Perivolaropoulos.  Observational constraints on viable f(R) parametrizations with geometrical and dynamical probes.  _Phys. Rev. D_ , 87(12):123529, 2013.  doi: 10.1103/PhysRevD.87.123529.
  * Nunes et al. [2017] Rafael C. Nunes, Supriya Pan, Emmanuel N. Saridakis, and Everton M. C. Abreu.  New observational constraints on $f(R)$ gravity from cosmic chronometers.  _JCAP_ , 01:005, 2017.  doi: 10.1088/1475-7516/2017/01/005.
  * Odintsov et al. [2017] Sergei D. Odintsov, Diego Sáez-Chillón Gómez, and German S. Sharov.  Is exponential gravity a viable description for the whole cosmological history?  _Eur. Phys. J. C_ , 77(12):862, 2017.  doi: 10.1140/epjc/s10052-017-5419-z.
  * Pérez-Romero and Nesseris [2018] Judit Pérez-Romero and Savvas Nesseris.  Cosmological constraints and comparison of viable $f(R)$ models.  _Phys. Rev. D_ , 97(2):023525, 2018.  doi: 10.1103/PhysRevD.97.023525.
  * D’Agostino and Nunes [2019] Rocco D’Agostino and Rafael C. Nunes.  Probing observational bounds on scalar-tensor theories from standard sirens.  _Phys. Rev. D_ , 100(4):044041, 2019.  doi: 10.1103/PhysRevD.100.044041.
  * D’Agostino and Nunes [2020] Rocco D’Agostino and Rafael C. Nunes.  Measurements of $H_{0}$ in modified gravity theories: The role of lensed quasars in the late-time Universe.  _Phys. Rev. D_ , 101(10):103505, 2020.  doi: 10.1103/PhysRevD.101.103505.
  * Odintsov et al. [2021] Sergei D. Odintsov, Diego Sáez-Chillón Gómez, and German S. Sharov.  Analyzing the $H_{0}$ tension in $F(R)$ gravity models.  _Nucl. Phys. B_ , 966:115377, 2021.  doi: 10.1016/j.nuclphysb.2021.115377.
  * Farrugia et al. [2021] Christine R. Farrugia, Joseph Sultana, and Jurgen Mifsud.  Spatial curvature in f(R) gravity.  _Phys. Rev. D_ , 104(12):123503, 2021.  doi: 10.1103/PhysRevD.104.123503.
  * Leizerovich et al. [2022] Matías Leizerovich, Lucila Kraiselburd, Susana J. Landau, and Claudia G. Scóccola.  Testing f(R) gravity models with quasar x-ray and UV fluxes.  _Phys. Rev. D_ , 105(10):103526, 2022.  doi: 10.1103/PhysRevD.105.103526.
  * Sultana et al. [2022] Joseph Sultana, Manoj K. Yennapureddy, Fulvio Melia, and Demosthenes Kazanas.  Constraining f(R) models with cosmic chronometers and the H ii galaxy Hubble diagram.  _Mon. Not. Roy. Astron. Soc._ , 514(4):5827–5839, 2022.  doi: 10.1093/mnras/stac1713.
  * De Felice and Tsujikawa [2010] Antonio De Felice and Shinji Tsujikawa.  f(R) theories.  _Living Rev. Rel._ , 13:3, 2010.  doi: 10.12942/lrr-2010-3.
  * Amendola and Tsujikawa [2015] Luca Amendola and Shinji Tsujikawa.  _Dark Energy: Theory and Observations_.  Cambridge University Press, 1 2015.  ISBN 978-1-107-45398-2.
  * Hu and Sawicki [2007] Wayne Hu and Ignacy Sawicki.  Models of f(R) Cosmic Acceleration that Evade Solar-System Tests.  _Phys. Rev. D_ , 76:064004, 2007.  doi: 10.1103/PhysRevD.76.064004.
  * Starobinsky [2007] Alexei A. Starobinsky.  Disappearing cosmological constant in f(R) gravity.  _JETP Lett._ , 86:157–163, 2007.  doi: 10.1134/S0021364007150027.
  * Cognola et al. [2008] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani, and S. Zerbini.  A Class of viable modified f(R) gravities describing inflation and the onset of accelerated expansion.  _Phys. Rev. D_ , 77:046009, 2008.  doi: 10.1103/PhysRevD.77.046009.
  * Linder [2009] Eric V. Linder.  Exponential Gravity.  _Phys. Rev. D_ , 80:123528, 2009.  doi: 10.1103/PhysRevD.80.123528.
  * Chen et al. [2015] Yun Chen, Chao-Qiang Geng, Chung-Chi Lee, Ling-Wei Luo, and Zong-Hong Zhu.  Constraints on the exponential $f(R)$ model from latest Hubble parameter measurements.  _Phys. Rev. D_ , 91(4):044019, 2015.  doi: 10.1103/PhysRevD.91.044019.
  * Tsujikawa [2008] Shinji Tsujikawa.  Observational signatures of $f(R)$ dark energy models that satisfy cosmological and local gravity constraints.  _Phys. Rev. D_ , 77:023507, 2008.  doi: 10.1103/PhysRevD.77.023507.
  * Phillips [1993] M. M. Phillips.  The absolute magnitudes of Type IA supernovae.  _Astrophys. J. Lett._ , 413:L105–L108, 1993.  doi: 10.1086/186970.
  * Scolnic et al. [2022] Dan Scolnic et al.  The Pantheon+ Analysis: The Full Data Set and Light-curve Release.  _Astrophys. J._ , 938(2):113, 2022.  doi: 10.3847/1538-4357/ac8b7a.
  * Scolnic et al. [2018] D. M. Scolnic et al.  The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample.  _Astrophys. J._ , 859(2):101, 2018.  doi: 10.3847/1538-4357/aab9bb.
  * Conley et al. [2011] A. Conley et al.  Supernova Constraints and Systematic Uncertainties from the First 3 Years of the Supernova Legacy Survey.  _Astrophys. J. Suppl._ , 192:1, 2011.  doi: 10.1088/0067-0049/192/1/1.
  * Jimenez and Loeb [2002] Raul Jimenez and Abraham Loeb.  Constraining cosmological parameters based on relative galaxy ages.  _Astrophys. J._ , 573:37–42, 2002.  doi: 10.1086/340549.
  * Wei [2018] Jun-Jie Wei.  Model-independent Curvature Determination from Gravitational-Wave Standard Sirens and Cosmic Chronometers.  _Astrophys. J._ , 868(1):29, 2018.  doi: 10.3847/1538-4357/aae696.
  * Simon et al. [2005] Joan Simon, Licia Verde, and Raul Jimenez.  Constraints on the redshift dependence of the dark energy potential.  _Phys. Rev. D_ , 71:123001, 2005.  doi: 10.1103/PhysRevD.71.123001.
  * Moresco et al. [2016] Michele Moresco, Lucia Pozzetti, Andrea Cimatti, Raul Jimenez, Claudia Maraston, Licia Verde, Daniel Thomas, Annalisa Citro, Rita Tojeiro, and David Wilkinson.  A 6% measurement of the Hubble parameter at $z\sim 0.45$: direct evidence of the epoch of cosmic re-acceleration.  _JCAP_ , 05:014, 2016.  doi: 10.1088/1475-7516/2016/05/014.
  * Moresco [2015] Michele Moresco.  Raising the bar: new constraints on the Hubble parameter with cosmic chronometers at z $\sim$ 2\.  _Mon. Not. Roy. Astron. Soc._ , 450(1):L16–L20, 2015.  doi: 10.1093/mnrasl/slv037.
  * Zhang et al. [2014] Cong Zhang, Han Zhang, Shuo Yuan, Tong-Jie Zhang, and Yan-Chun Sun.  Four new observational $H(z)$ data from luminous red galaxies in the Sloan Digital Sky Survey data release seven.  _Res. Astron. Astrophys._ , 14(10):1221–1233, 2014.  doi: 10.1088/1674-4527/14/10/002.
  * Stern et al. [2010] Daniel Stern, Raul Jimenez, Licia Verde, Marc Kamionkowski, and S. Adam Stanford.  Cosmic Chronometers: Constraining the Equation of State of Dark Energy. I: H(z) Measurements.  _JCAP_ , 02:008, 2010.  doi: 10.1088/1475-7516/2010/02/008.
  * Moresco et al. [2012] M. Moresco et al.  Improved constraints on the expansion rate of the Universe up to z~1.1 from the spectroscopic evolution of cosmic chronometers.  _JCAP_ , 08:006, 2012.  doi: 10.1088/1475-7516/2012/08/006.
  * Ratsimbazafy et al. [2017] A. L. Ratsimbazafy, S. I. Loubser, S. M. Crawford, C. M. Cress, B. A. Bassett, R. C. Nichol, and P. Väisänen.  Age-dating Luminous Red Galaxies observed with the Southern African Large Telescope.  _Mon. Not. Roy. Astron. Soc._ , 467(3):3239–3254, 2017.  doi: 10.1093/mnras/stx301.
  * Borghi et al. [2022] Nicola Borghi, Michele Moresco, and Andrea Cimatti.  Toward a Better Understanding of Cosmic Chronometers: A New Measurement of H(z) at z $\sim$ 0.7.  _Astrophys. J. Lett._ , 928(1):L4, 2022.  doi: 10.3847/2041-8213/ac3fb2.
  * Hu and Sugiyama [1996] Wayne Hu and Naoshi Sugiyama.  Small scale cosmological perturbations: An Analytic approach.  _Astrophys. J._ , 471:542–570, 1996.  doi: 10.1086/177989.
  * Eisenstein and Hu [1998] Daniel J. Eisenstein and Wayne Hu.  Baryonic features in the matter transfer function.  _Astrophys. J._ , 496:605, 1998.  doi: 10.1086/305424.
  * Beutler et al. [2011] Florian Beutler, Chris Blake, Matthew Colless, D. Heath Jones, Lister Staveley-Smith, Lachlan Campbell, Quentin Parker, Will Saunders, and Fred Watson.  The 6dF Galaxy Survey: Baryon Acoustic Oscillations and the Local Hubble Constant.  _Mon. Not. Roy. Astron. Soc._ , 416:3017–3032, 2011.  doi: 10.1111/j.1365-2966.2011.19250.x.
  * Percival et al. [2010] Will J. Percival et al.  Baryon Acoustic Oscillations in the Sloan Digital Sky Survey Data Release 7 Galaxy Sample.  _Mon. Not. Roy. Astron. Soc._ , 401:2148–2168, 2010.  doi: 10.1111/j.1365-2966.2009.15812.x.
  * Padmanabhan et al. [2012] Nikhil Padmanabhan, Xiaoying Xu, Daniel J. Eisenstein, Richard Scalzo, Antonio J. Cuesta, Kushal T. Mehta, and Eyal Kazin.  A 2 per cent distance to $z$=0.35 by reconstructing baryon acoustic oscillations - I. Methods and application to the Sloan Digital Sky Survey.  _Mon. Not. Roy. Astron. Soc._ , 427(3):2132–2145, 2012.  doi: 10.1111/j.1365-2966.2012.21888.x.
  * Tojeiro et al. [2014] Rita Tojeiro et al.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: galaxy clustering measurements in the low redshift sample of Data Release 11.  _Mon. Not. Roy. Astron. Soc._ , 440(3):2222–2237, 2014.  doi: 10.1093/mnras/stu371.
  * Anderson et al. [2013] Lauren Anderson et al.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations in the Data Release 9 Spectroscopic Galaxy Sample.  _Mon. Not. Roy. Astron. Soc._ , 427(4):3435–3467, 2013.  doi: 10.1111/j.1365-2966.2012.22066.x.
  * Bautista et al. [2018] Julian E. Bautista et al.  The SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations at redshift of 0.72 with the DR14 Luminous Red Galaxy Sample.  _Astrophys. J._ , 863:110, 2018.  doi: 10.3847/1538-4357/aacea5.
  * Seo et al. [2012] Hee-Jong Seo et al.  Acoustic scale from the angular power spectra of SDSS-III DR8 photometric luminous galaxies.  _Astrophys. J._ , 761:13, 2012.  doi: 10.1088/0004-637X/761/1/13.