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

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 H\textsc{ii} starburst galaxies, and local measurements of the Hubble parameter $H_0$. 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 [1][2][3][4].Furthermore, observation of baryon acoustic oscillations [5,6], analysis of cosmic microwave background radiation [7][8][9][10], and examination of the power spectrum of matter distributions in the universe [11,12] 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 Λ term into Einstein's equation, known as the "Λ Cold Dark Matter (Λ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 [13] and the finetuning problem [14] 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 [15]) or by means of their kinetic energy (k-essence models [16]).
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 arc-Tanh 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 ΛCDM model), we chose to re-parameterize these models in terms of what is called "deviation/distortion parameter(b)".In this form, with b → 0, a f (R) model tends to the Λ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 [17], 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 [18], 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, [19] 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 [20].Other significant works that investigated Hu-Sawicki and/or exponential models, utilizing various data sets including gravitational lensing data, are [21][22][23].In [24], constraints on Hu-Sawicki, Starobinsky, exponential, and Tsujikawa models were examined in both flat and non-flat spacetimes, utilizing various cosmological data sets.Recently, [25] 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 [17] was extended to include the arcTanh model in [26].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 III and IV, 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 −1 Mpc −1 .
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 where κ = 8πG, g is the determinant of the metric tensor, S m and S 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 where F = ∂f /∂R, □ ≡ g µν ∇ µ ∇ ν is the covariant D'Alambertian and T µν is the energy-momentum tensor of matter and radiation.This field equation can also be expressed as where G µν ≡ R µν − 1 2 Rg µν is Einstein tensor and prime( ′ ) denotes derivative with respect to R. Taking trace of both sides of Eq. 2 we obtain using which in Eq. 3 we rewrite the field equation as where T ≡ T µ µ 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 where a is the FLRW scale factor.
The Hubble parameter(H) is defined as H ≡ ȧ/a where we use a generic notation Ẋ 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 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 ≤ z < 10 4 ) when any interaction between matter and radiation could be ignored, the energymomentum tensor for the fluid can be written as T µ ν = diag(−ρ, p, p, p) with ρ = ρ m + ρ r and p = p r (with radiation pressure p r = ρ 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 ρm + 3Hρ m = 0, ρr + 4Hρ r = 0 . ( The solution to these conservation equations are ρ m = ρ m0 /a 3 , ρ r = ρ 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 ρ and p, Eq. 4 reduces to (9) using which the '00' and 'ii' components of the field Eq. 5 take following respective forms: 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 ρ DE and pressure p DE , i.e.
with the equations 3H 2 = κ (ρ m + ρ r + ρ DE ) and −2 Ḣ = κ (ρ m + ρ r + ρ DE + p m + p r + p DE ), we can deduce the "effective (geometric) dark energy" with density and pressure corresponding to f (R)-theory as and, with the equation-of-state parameter for this effective dark energy defined as Using Eqs. 10 and 11, we may recast Eq. 14 into a more computationally advantageous form (which we use later) as where and we have taken p 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 tot |z = 0 ∼ −0.7, w DE | z=0 ∼ −1, and a transition redshift (zt) ∼ 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 In the ΛCDM model, the value of w DE is fixed at −1 for all redshifts due to the presence of the cosmological constant (Λ) term.However, in f (R) models, the source of w DE | z=0 ∼ −1 arises from the underlying geometry itself.Unlike the Λ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 Λ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 [27,28] for detailed discussions): and, where R 0 denotes the Ricci scalar at the present epoch.
These conditions arise from various considerations.Firstly, the effective gravitational constant (G eff = G/F ) must be positive, ensuring the avoidance of anti-gravity (as expressed in Eq. 18).Besides, any acceptable f (R)model should exhibit stability under perturbations and avoid the instability of Dolgov-Kawasaki type (Eq.19).Similar to the ΛCDM model, a viable f (R) model must also be consistent with local gravity tests (Eqs.19 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 19 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: where the function y(R, b, Λ) serves to measure the deviation from the ΛCDM model, with the parameter b (referred to as the "distortion/deviation parameter") specifically quantifying the extent of the deviation [17].
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.19 and 20 can be alternatively written as lim R→∞ 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 ΛCDM model (which is supported by observations of the cosmic microwave background), we can identify that C 0 = −2Λ, where Λ corresponds to the cosmological constant introduced by Einstein and Hilbert in their action.This amounts to writing symbolically Λ Λ = Λ f (R) where the superscripts Λ and f (R) denote quantities in reference to the ΛCDM model and any viable f (R) -model, respectively.This relation may also be written as where subscript '0' denotes values at present time.Furthermore, considering the definition of the energymomentum tensor and the resulting conservation equation(Eq.8), we can infer that both the ΛCDM model and any viable f (R)-model yield identical matter density and radiation density at the current epoch i.e.
where the subscript i = (m, r) stands for (matter, radiation).Also any viable f (R)-model is expected to exhibit deviations from the standard ΛCDM model predominantly during the late times, in order that the explanation for accelerated expansion at present epoch comes from vanishing Λ.In general terms, this implies 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 (Ω Λ m0 , b, H Λ 0 ) for model-fitting to the data.However, while reporting our findings, we express the results in terms of the parameters (Ω ), which are determined using Eq.24.

A. The Hu-Sawicki Model
The Hu-Sawicki model, initially proposed in [29], is described by the following equation: where c 1 and c 2 are dimensionless parameters, n HS is a positive constant typically assumed to be an integer, and /c 1 ≡ b, we can express Eq. 26 as [17]: We identify the parameter Λ as the usual cosmological constant, and b as the deviation parameter, which indicates the model's deviation from the ΛCDM model.In order to satisfy the viability conditions F > 0 and F ′ > 0 for R ≥ R 0 , it is necessary to consider b > 0 (especially when n HS is an odd integer).Although many researchers have constrained the scenarios where n HS = 1 and/or 2 [17,18,21,22,[24][25][26], they have acknowledged computational challenges as a hindrance to explore the case of n HS = 3 (further elaboration on this point is provided in Sec.V).In this investigation, we have also imposed constraints on the case n HS = 3.

B. The Starobinsky Model
The model proposed by Starobinsky [30] is where n S is a positive constant, λ(> 0) and R S ≈ 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 [17] f

C. The Exponential Model
The exponential model, initially proposed as a viable f (R) model in [31], has been further investigated in [19,25,32,33] (also see references therein), is given by where α and β are the parameters of this model.For large R, an acceptable f (R)-model must approximately resemble the ΛCDM, which is achievable only when α > 0 and β > 0. By substituting Λ = α/2 and b = 2/(αβ) the exponential model can be expressed as from where it becomes evident that when R becomes significantly larger than bΛ (R ≫ bΛ), the function f (R) E approaches R − 2Λ.

D. The Tsujikawa Model
Tsujikawa proposed an alternative model [34] as where ξ(> 0) and R T (> 0) are the model parameters.
With Λ = ξR T /2 and b = 2/ξ, the model can be rewritten as We see clearly that when the parameter b → 0 (which corresponds to ξ → ∞, R T → 0, while ξR T remains finite), the model reduces to f (R) T = R − 2Λ.

E. The ArcTanh Model
In this study we also examined a model proposed in [20] as 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.

A. Type Ia Supernova Data
Type Ia Supernovae (SNIa), known as standard candles [35], have significantly contributed to our comprehension of cosmology.The SNIa observations [2,4] 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 [36].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 HD < 2.26137, where z HD denotes the Hubble diagram redshift.It offers a considerably larger number of low-redshift data compared to the previous Pantheon compilation [37].In our current work, we employ the apparent magnitude at maximum brightness (m b ), heliocentric redshift (z hel ), cosmic microwave background corrected redshift (z cmb ), and the total (statistical + systematic) covariance matrix from this compilation [36].
The theoretical definition of the apparent magnitude involves the luminosity distance, which is given by the integral expression: where the function H(z) is obtained by solving the system of Eqs. 7, 10 and 11.The apparent magnitude is defined as where M represents the absolute magnitude, and 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 [38] where and C represents the total covariance matrix of the data (provided in the Pan-theon+ compilation). 1 is an array of ones of length equal to the number of data points.

B. Cosmic Chronometers Data
By examining the differential age evolution [39][40][41] 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[41][42][43][44][45][46][47][48] (compiled in the Table IV of the Appendix).This so called differential age method, employed in these studies, uses the relation to obtain H(z).The parameters of any model are estimated by minimizing the following residual: where H obs (z i )'s are the observed values of the Hubble parameter function at redshift z i , while σ 2 H,i denotes the corresponding uncertainties associated with the measurements of H 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 th (z i ) in the above Eq.39.

C. 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 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 d ).BAO, therefore, holds the status of standard ruler for length scale in Cosmology [49,50].
We have compiled a collection of 30 data points representing various BAO observables from a range of surveys, as documented in the literature [51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70].These data points, which are used in our current study, are listed in the Table III(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 [71].The BAO observables include the Hubble distance, which is defined as D H = c/H(z), the transverse comoving distance (D M (z)), the angular diameter distance (D A (z)), and the volume-averaged distance (D V (z)) which are defined as and, 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 where A obs (z i ) and σ i respectively denote the observed values and their uncertainties at redshift z i and A th denotes the theoretical prediction from model under consideration.These quantities are given in the columns 2-4 of the Table III.For correlated data points the appropriate residual to be minimized is where C j 's are the covariance matrices of the 4 different data-sets and (A th − A obs ) 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 III.For the whole data-set of BAO, the χ 2 which is to be minimized is D. 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 Hα and Hβ lines, along with a weak continuum.The cosmological significance of HiiG observation comes from the empirically established correlation between the luminosity (L) of the Hβ lines and the velocity dispersion (σ, 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 σ) (see [72] and references therein).
In this study, we have used a total of 181 data points from [72][73][74][75][76] corresponding to the emission of the Balmer Hβ 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 SN > 1.4, while for z HiiG > 1.4 we have 69 data points, and for z HiiG > z 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 σ is given by log where we take β = 5.022 ± 0.058 and α = 33.268± 0.083 from [72,76].From the luminosity distance of HiiG, 1/2 , we can derive the distance modulus as (47) Using Eq. 47, we compute the distance moduli of HiiG from the observables L, σ, and F .For any given cosmological model with theoretical distance moduli µ θ (represented as m th − M in Eq. 36), the parameters associated with that model can be constrained by minimising the following χ 2 function where with and ϵ sys = 0.0.257 as suggested in [72,76,77].The uncertainty contribution of ϵ µ θ,stat in the distance modulus is due to the uncertainties in redshifts(σ z ∼ 10 −4 ) of HiiG and is derived from simple error propagation theory.

E. Local Measurement of H0
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± 0.5 [78], while the locally measured value by the SH0ES collaboration yields H 0 = 73.04 ± 1.04 [79].The discrepancy between these two measurements, known as the Hubble tension, is significant at the 4 − 5σ level.To address this tension, we have included in our study both scenarios -with visa-vis without the SH0ES prior for H 0 .The contribution of this data point to the total χ 2 is given by

F. 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 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 χ 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: (Ω Λ m0 , b, H Λ 0 ).We used the uniform priors Ω Λ m0 ∈ [0, 1] and H Λ 0 ∈ [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 Λ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 HS = 3)(HS3), followed by the Starobinsky(n S = 1)(ST1) and the exponential(EXP) models, and the least similar ones are the Hu-Sawicki model(n HS = 1)(HS1) and the arcTanh(aTanh) models.To allow for the exploration that whether the data supports a model different from the ΛCDM model, an f (R) model which is more similar to the ΛCDM requires a bigger change in b to make them noticeably different from the ΛCDM model.Consequently we use uniform priors b ∈ [0, b max ] where b 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 [80], (ii) to perform MCMC analysis -EMCEE [81], and (iii) for plotting of posterior probability distributions of the parameters -GetDist [82].

V. RESULTS: OBSERVATIONAL CONSTRAINTS ON f (R) MODELS
In this section we will discuss our findings regarding the parameters of the six f (R) models and the Λ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: SCH 0 to denote SNIa+CC+H 0 and SC(H 0 ) to denote SNIa+CC or SNIa+CC+H 0 , depend- ing 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) SCHii(H 0 )-Black, (iv) BCHii(H 0 )-Orange, and (v) SBCHii(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 I and II ing 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) Ω m0 : For all models and data-set combinations, either the median values of Ω m0 fall within the 1-2 sigma interval of the Planck value (Ω m0,Planck = 0.315 ± 0.007 [78]) or, the Planck value is within the 1-2 sigma intervals of the model's median values of Ω 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 ), BCHii(H 0 ), and SBCHii(H 0 ) data-sets, whereas opposite trends are observed for the SC(H 0 ) and SCHii(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,Planck = 67.4± 0.5 [78]) 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,SH0ES = 73.04 ± 1.04 [79].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 Ω m0 or b).Additionally, there is less overall tension among the model-fitted values of H 0 compared to the tension between H 0,Planck and H 0,SH0ES .

A. Constraints on ΛCDM Model
We have also included the results for the ΛCDM model, which is commonly used as a benchmark for assessing f (R) models.We consider a two-parameter ΛCDM model, with the Hubble parameter given by H = H 0 Ω m0 (1 + z) 3 + (1 − Ω m0 ), where the two parameters Ω 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 I and II.The posterior probability distribution of Ω m0 and H 0 , both for the respective cases -without and with the SH0ES prior for H 0 -are depicted in Figs 4 and  TABLE I. Results from the MCMC fitting process(for the cases without SH0ES prior for H0): The median values of the model parameters along with 1-sigma(68.26%)confidence intervals on them, the transition redshift(zt), the current values of EoS parameter for geometric/dark energy component(wDE0) and total EoS parameter(wtot0), the minimum values of χ 2 and other quantities for making statistical inferences such as the reduced chi-square value χ 2 red = χ 2 min /ν(where the number of degrees of freedom ν = 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), ∆AIC = AIC 5.For all data sets (except SC), the values of Ω m0 are compatible with Ω m0,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,Planck or H 0,SH0ES .

B. Constraints on Hu-Sawicki Model and Starobinsky Model
Previous studies (see [18] and references therein) have established that there exists a degeneracy between n HS /n S and Ω m0 .One, therefore , often works with fixed values of n HS and n S to address this degeneracy.In [83], it has been worked out that n HS , n S > 0.9, while TABLE II.Results from the MCMC fitting process(for the cases with SH0ES prior for H0): The median values of the model parameters along with 1-sigma(68.26%)confidence intervals on them, the transition redshift(zt), the current values of EoS parameter for geometric/dark energy component(wDE0) and total EoS parameter(wtot0), the minimum values of χ 2 and other quantities for making statistical inferences such as the reduced chi-square value χ 2 red = χ 2 min /ν(where the number of degrees of freedom ν = 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), ∆AIC = AIC  [84] proposed that n HS and n S should be integers.The difficulty in constraining higher n HS /n 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 HS = 1 and n S = 1.In our study we have also explored the case of n HS = 3.Note that in the re-parameterized version using the deviation parameter b, the Hu-Sawicki model  The Fig. 1 illustrates that as n HS increases (from 1 to 3) the values of Ω m0 approach closer to those from the ΛCDM model (and/or to Ω m0, 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 HS causes the model values of H 0 to approach H 0, Planck , except for the SC dataset where there is almost no change with varying n HS .However, the constraints from the BCHii dataset are not considered reliable due to high χ 2 /ν.When the SH0ES prior is included, there is hardly any change in the model values of H 0 with respect to n HS , and they lie around 69.5−71.5 km s −1 Mpc −1 .As expected, the parameter b increases as n HS increases (see Fig. 2).Also, the constraints on b become weaker.We may infer that for higher values of n HS (say 4, 5, ...) and n S (say 2, 3, ..), the constraints on b would become even more weaker.This implies that with further increase in n HS or n S , one can obtain a model that closely resembles the ΛCDM model in terms of its predictions for the physical parameters Ω m0 and H 0 , despite the constraints on b moving towards larger values.In a sense, this undermines the purpose of explor- ing f (R) models further, and there is a strategic motivation to limit the exploration to n HS = 1, 2 only, putting aside the consideration of computational difficulties in constraining for the higher values of n HS or n 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 ΛCDM model.This finding is also an important result of the present work.

C. Constraints on Exponential Model
In the Tables I and II, 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 BCHii(H 0 ) dataset, the parameter b clearly deviates from zero, with b = 0 barely allowed.Despite this, the values of the parameters Ω m0 and H 0 (as illustrated in Figs. 1 and  3) are reasonably close to the standard values derived from Planck constraints [78].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 ∼ 69.7 to 71.5 km s −1 Mpc −1 .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.

D. Constraints on Tsujikawa Model
The constraints on the parameters of the Tsujikawa model described by Eq. 33 are presented in the Tables I and II, 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 BCHii(H 0 ).We disregard the results from BCHii(H 0 ) on account of very high value of the reduced χ 2 min /ν.Hence, we obtain an f (R) that is clearly distinguishable from the ΛCDM model.The fitted values of Ω m0 and H 0 also fall within reasonable limits compared to Ω m0,Planck and H 0,Planck (or H 0,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 −1 Mpc −1 , which introduces a tension of  2-3 sigma with H 0,Planck or H 0,SH0ES .

E. Constraints on arcTanh Model
The results obtained by constraining the aTanh model are presented in the Tables I and II, 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 BCHii(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 Ω m0 for the model fall within the 1-2 sigma range of Ω m0,Planck .Furthermore, H 0,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 −1 Mpc −1 , which are in 2-3 sigma tension with H 0,Planck or H 0,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 arctanh(Λ/R) ≈ Λ/R for Λ/R ≪ 1 (which holds true in this case since Λ ∼ 0.7H 2 0 and R ≳ 8H 2 0 ).

VI. MODEL COMPARISON
The standard statistical tools commonly used to assess model fitting (and model comparison) in cosmology include the reduced chi-square statistics (χ 2 ν ), the Akaike Information Criterion (AIC) [85], and the Bayesian Information Criterion (BIC) [86] (also see [87][88][89]).The last six columns of Tables I and II contain essential quantities related to the current analysis, specifically pertaining to the utilization of these statistical tools.The reduced chisquare statistics is defined as χ 2 ν = χ 2 min /ν whereas the AIC and the BIC are respectively defined by the following equations The number of degrees of freedom, represented by the symbol ν, is determined by subtracting the number of model parameters (k) from the total number of data points (N ).The minimum value of χ 2 is denoted as χ 2 min and is connected to the maximum likelihood (L max ) through the relation χ 2 min = −2 ln L max .While comparing multiple competing models using a given data set, the model with the lowest values of χ 2 , AIC, and BIC is considered to be more favoured by the data.However, it is insufficient to rely solely on χ 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 χ 2 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 χ 2 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 [87][88][89]).
To perform a comparative analysis of two models, a useful measure ∆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 < ∆X ⩽ 2 (−2 ⩽ ∆X < 0): weak evidence, (ii) 2 < ∆X ⩽ 6 (−6 ⩽ ∆X < −2): positive evidence, (iii) 6 < ∆X ⩽ 10 (−10 ⩽ ∆X < −6): strong evidence, (iv) ∆X > 10 (∆X < −10): very strong (highly pronounced) evidence.
As can be seen from Tables I and II, the reduced χ 2 values for the BCHii(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 I.The ∆AICs values indicate varying degrees of evidence ranging from weak to very strong in favor of f (R) models compared to the ΛCDM model.On the other hand, the values of ∆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 ∆AICs and ∆BICs align, guiding us in the same direction.
When the SH0ES prior for H 0 is included, the results, as shown in Table II, exhibit some differences.The overall support for f (R) models weakens compared to the previous case.Based on consideration of the ∆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 ∆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 Λ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 SBCHiiH 0 data-set, followed by the SBCH 0 data-set.Furthermore, Table II 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 ΛCDM model) is weak or mildly positive according to ∆BIC.Conversely, for models such as HS1, aTanh, and ST (which are less or least similar to the ΛCDM model), the evidence against them is strong when evaluated using ∆BIC.Consequently, we can conclude that with a SH0ES prior for H 0 , the ∆BIC criterion disfavours models that significantly deviate from the Λ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 ΛCDM model (as indicated by ∆AIC and/or ∆BIC) also correspond to the cases where a nonzero value of the deviation parameter b is 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 I and II.

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) ≡ −aä/ ȧ2 = −ä/(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 t ) > 0 [90,91].The quantity z 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 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 eff takes values approximately equal to 1/3, 0, and −1 (for ΛCDM).On the other hand, the EoS parameter (w DE ) associated with the dark energy component remains constant at −1 for all values of redshift (both in the past and future) in ΛCDM model.However, in the case of any viable f (R) model, value of the parameter w DE can cross the so called "phantom-divide" line (w DE (z) = −1) multiple times before finally settling at −1 in the far future [92,93].
In the Tables I and II we present the median values and 1-sigma confidence intervals for the quantities z t , w DE,0 , and w eff,0 for two cases -without and with the SH0ES prior on H 0 .In Fig 18, we display the values of z t for all models and data-sets.This figure shows the median values of z t ∼ 0.5 − 1 for all the models, which is compatible with the model-independent predictions for z t [90,91].In the Figs.19 and 20 we have plotted the evolution of w DE with redshift based on the best constraints obtained for the models, (which happens to be the constraints from the data-set SCHii(H 0 )).In all f (R) models, the values of w 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 DE,0 ) falling within the quintessence region.When the SH0ES prior on H 0 is included, the deviations of w DE (z) for f (R) models from the Λ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 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 [92,93] 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 DE (z) curve up to z → −1(a → ∞), 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 eff (z) is depicted in Figs.21 and  22, revealing a transition from a matter-dominated era (represented by w eff (z) → 0) to a dark energy dominated era around z ∼ 5 for all the f (R) models.Considering the SH0ES prior for H 0 leads to relatively lower values of w eff,0 for all the f (R) models compared to the cases without the prior.Furthermore, the discrepancy between the predictions of w eff (z) obtained from the ΛCDM model and those from the f (R) models reduces when the SH0ES prior is included.

VIII. CONCLUSIONS
In this study, we used the recently released Pantheonplus SNIa data along with cosmic chronometer measure- ments, 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 ≲ z ≲ 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 ΛCDM  The six viable f (R) models examined in this study are the Hu-Sawicki model with n HS = 1, 3, the Starobinsky model (n 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 → 0, these models converge to the standard ΛCDM model.Each of these f (R) models is characterized Regardless of whether the SH0ES prior for H 0 is considered, the Ω m0 values for all the examined models fall within the 1-2 sigma range of the Planck constrained value (Ω m0,Planck = 0.315 ± 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,Planck = 67.4± 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,SH0ES = 73.04 ± 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 ∆AIC and/or ∆BIC based analysis suggest that the f (R) models are (very) strongly favoured over the Λ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 ∆AIC and/or ∆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 t , w eff,0 , and w 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 t ∼ 0.5 − 1.Furthermore, the current values of w DE fall within the quintessential region, having recently crossed over from the phantom region (around z ∼ 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 t , w eff,0 , and w 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.
with Λ = λR S /2 and b = 2/λ.In this study, we have imposed constraints on the cases where n S = 1 and the reason for not exploring higher values of n S is explained later in Sec.V. Note that the Hu-Sawicki model (Eq.27) with n HS = 2 and the Starobinsky model (Eq.29) with n S = 1 are equivalent.Unlike the Hu-Sawicki model (with n 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 ΛCDM model supported by the data, we considered b > 0 without loss of generality.

FIG. 1 .
FIG.1.This figure is helpful in getting a sense of variations of Ω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 H0, 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 2sigma(95.44%,with bigger cap size) confidence intervals for the cases without SH0ES prior.

FIG. 2 .
FIG.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 H0, 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 2sigma(95.44%,with bigger cap size) confidence intervals for the cases without SH0ES prior.

FIG. 3 .
FIG. 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 H0, 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.

FIG. 4 .
FIG. 4. The ΛCDM Model(without SH0ES prior for H0):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.

FIG. 5 .
FIG. 5.The ΛCDM Model(with SH0ES prior for H0):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.

FIG. 6 .
FIG.6.The Hu-Sawicki Model(n HS = 1, without SH0ES prior for H0): 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.

FIG. 7 .
FIG. 7. The Hu-Sawicki Model(n HS = 1, with SH0ES prior for H0): 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.

FIG. 8 .
FIG. 8.The Starobinsky Model(n S = 1, without SH0ES prior for H0): 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.

FIG. 9 .FIG. 10 .
FIG. 9.The Starobinsky Model(n S = 1, with SH0ES prior for H0): 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.

FIG. 11 .
FIG. 11.The Hu-Sawicki Model(n HS = 3, SH0ES prior for H0): 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.

FIG. 12 .
FIG.12.The Exponential Model(without SH0ES prior for H0): 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.

FIG. 13 .
FIG. 13.The Exponential Model(with SH0ES prior for H0):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.

FIG. 14 .
FIG. 14.The Tsujikawa Model(without SH0ES prior for H0): 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.

FIG. 15 .
FIG. 15.The Tsujikawa Model(with SH0ES prior H0):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.

FIG. 16 .
FIG. 16.The arcTanh Model(without SH0ES prior for H0):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.

FIG. 17 .
FIG. 17.The arcTanh Model(with SH0ES prior for H0):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.

FIG. 18 .
FIG.18.Transition Redshift: The blank star and blank circle markers represent median values for the cases with and without SH0ES prior for H0, 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.

FIG. 19 .FIG. 20 .
FIG.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. model

FIG. 21 .
FIG.21.Evolution of the total EoS parameter for all f (R) models and the Λ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.

FIG. 22 .
FIG.22.Evolution of the total EoS parameter for all f (R) models and the ΛCDM model obtained from the data-set SBCHiiH0.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.
f (R) model -AIC ΛCDM model and ∆BIC = BIC f (R)model -BIC ΛCDM model -for all the models and data set combinations considered in this work.
f (R) model -AIC ΛCDM model and ∆BIC = BIC f (R)model -BIC ΛCDM model -for all the models and data set combinations considered in this work.

TABLE IV .
The CC data collected from many cited sources.