Analyzing the $H_0$ tension in $F(R)$ gravity models

The Hubble constant tension problem is analysed in the framework of a class of modified gravity, the so-called $F(R)$ gravity. To do so, we explore two models: an exponential and a power-law $F(R)$ gravities, which includes an early dark energy (EDE) term in the latter. These models can describe both an early time inflationary epoch and the late time accelerating expansion of the universe. We confront both models with recent observational data including the Pantheon Type Ia supernovae sample, the latest measurements of the Hubble parameter $H(z)$ from differential ages of galaxies (cosmic chronometers) and separately from baryon acoustic oscillations. Standard rulers data set from the Cosmic Microwave Background radiation are also included in our analysis. The estimations of the Hubble constant appear to be essentially depending on the set of observational data and vary in the range from 68 to 70.3 km/(s$\cdot$Mpc). The fits of the other free parameters of the models are also obtained, revealing interesting consequences.


I. INTRODUCTION
Among current problems in modern cosmology, the tension among estimations of the Hubble constant H 0 is one of the most striking and irritating for researchers.Over the last years such discrepancies among H 0 measurements have been revealed through two different methods: by the Planck collaboration after collecting and analyzing data from the cosmic microwave background radiation (CMB) over the last 7 years [1][2][3], which provides an estimation of H 0 = 67.4± 0.5 km s −1 Mpc −1 (Planck18), and on the other hand by the SH0ES group of Hubble Space Telescope (HST) [4,5] with the last estimate given by H 0 = 74.03± 1.42 km s −1 Mpc −1 (SH0ES19).The HST method includes measurements of the local distance ladder by combining photometry from Cepheids (and their period luminosity relation) with other local distance anchors, Milky Way parallaxes and calibration distances to Cepheids in the nearest galaxies which are hosts of Type Ia Supernovae (SNe Ia).In particular, the above estimation for the Hubble constant by the HST group includes observations of 70 Cepheids in the Large Magellanic Cloud [5].
Currently the mismatch among H 0 estimations by Planck [3] and HST [5] collaborations exceeds 4σ, as this tension has grown over the last years, as shown in Table I.This problem may be dealt as the discrepancy between observations at early and late cosmological time of our Universe [6], since HST group works with late time data while Planck collaboration combines observations from redshifts in a wide range 0 < z < 1100 and uses the standard ΛCDM model as fiducial model, but the issue may be approached through a theoretical way.For the former, some researchers have suggested some different ways for solving the H 0 tension problem.Several groups have analysed the estimations of H 0 by using several approaches independent of the Cepheid distance scale and CMB anisotropies (for a review see [6,7]).Among these methods with new observational results for H 0 , the following approaches may be highlighted: the tip of the red giant branch (TRGB) method used by the Carnegie-Chicago Hubble Program (CCHP) [8,9], lensing objects with strong time delays between multiple images (H0LiCOW project and others) [10][11][12], CMB-lensing data [13], maser (megamaser) hosting galaxies [14], oxygen-rich variable stars (Miras) [15].Some other researchers tried to explain the tension by assuming that Planck or HST measurements might suffer from systematic errors [16], but these analysis did not led to convincing solutions of the problem.I, one can note that most of H 0 estimations lie on the range among Planck18 and SH0ES19 values, while the local (late-time) H 0 measurements are close to the SH0ES19 value, exceed the early-Universe estimations.Only the CCHP estimation obtained by TRGB method violates the latter tendency, which have led to some discussions [6,9].By comparing these facts, many cosmologists over the recent years have considered the H 0 tension as a hint for new physics beyond the standard ΛCDM model with different phenomena in early and/or late times of the universe evolution [17] - [23] (see also the extended list of literature in Ref. [7]).These analysis suggest several ways for solving the problem, which can be generally summarized as a mechanism that shifts the effective H 0 value from early to late time Universe under different factors.The best fit H 0 appears to be essentially depending on the mentioned factors.Following this idea, some scenarios have been studied:

As shown in Table
• Dark energy models with a varying equation of state (EoS), via a varying EoS parameter w or via the dark energy density ρ DE [17].
• Scenarios with an early dark energy (EDE) component, reproduced in different frameworks (scalar fields, axions), which becomes important before the epoch of matter-radiation equality z ≃ 3000 and then decays after faster than radiation [18].
• Models with evolving or decaying dark matter into dark radiation or other species [19].
• Interacting dark energy and dark matter models [20].
• Models with extra relativistic species that modifies the effective number N eff at the recombination era [21].
• Modified gravity models that emerge at an intermediate epoch, as scalar-tensor theories, F (R) gravity, F (T ) and others [22][23][24].In this sense, modified gravities have been widely studied in the literature in the framework of cosmology (for a review see [25]).Particularly, F (R) gravity is very well known by the scientific community, with an extensive literature where many aspects of the theory have been analysed.This modification of GR assumes a generic function of the Ricci scalar for the gravitational action instead of the linear term of the Hilbert-Einstein action.Such modification leads to interesting properties and a rich phenomenology that can solve some of the most important problems in cosmology, as the origin of dark energy.In this sense, F (R) gravity can reproduce well the late-time acceleration with no need of additional fields and alleviate the cosmological constant problem by compensating the large value for the vacuum energy density predicted by quantum field theories (see [26][27][28][29][30][31]). In addition, the same type of modifications of GR have been studied in the framework of inflation where as in the case of dark energy, F (R) gravities can lead to successful scenarios that fit perfectly well the constraints on the spectral index of perturbations as given by the analysis of the CMB ( [32]).With this in mind, models that unify the dark energy epoch and the inflationary paradigm through the corrections introduced in the gravitational action have been proposed with a great success [33][34][35][36][37][38][39][40].In addition, some F (R) gravity models that reproduce late-time acceleration can also recover GR at local scales where this one is very well tested, leading to the so-called viable F (R) gravity models [28,29].Hence, F (R) gravity is particularly of great interest in cosmology.
Hence, supported on the the great knowledge of F (R) gravity and its success on trying to solve some of the most important problems in cosmology, the possibility of alleviating the H 0 tension problem in this framework may be promising, despite has not been studied yet exhaustively.Although some efforts are being done, as the analysis to solve this tension by viable F (R) models, and particularly through the Hu-Sawicki F (R) model [28] in Ref. [23], where concluded that the Hu-Sawicki gravity cannot reduce the H 0 tension.In the present paper, we analyse two F (R) models and the possibility of alleviating the H 0 tension.We confront the models with observational data and estimate the Hubble constant H 0 and other model parameters by using approaches developed in some previous papers [34,35,[41][42][43][44].Here we include in our analysis the following observations: the Type Ia supernovae data (SNe Ia) from the Pantheon sample survey [45], data connected with cosmic microwave background radiation (CMB) and extracted from Planck 2018 observations [3,46], and estimations of the Hubble parameter H(z) for different redshifts z from two different sources: (a) measured from differential ages of galaxies (in other words, from cosmic chronometers, these 31 data points are analyzed separately) and (b) H(z) obtained as observable effect of baryon acoustic oscillations (BAO).We obtain the best fit parameters and compare to the ones from ΛCDM.
The paper is organized as follows.In section II, we introduce F (R) gravity and the two models we analyse along the paper.Section III is devoted to SNe Ia, H(z) and CMB observational data.In section IV we analyze the results, estimations for the Hubble constant H 0 and other model parameters.Finally, section V gathers the conclusions of the paper.

II. F (R) GRAVITY MODELS
We can start by reviewing the basics of what is called F (R) gravity, a generalization of the Einstein-Hilbert action that assumes a more complex Lagrangian in terms of the Ricci scalar R: Here κ 2 = 8πG and S matter is the matter action.In the present paper, we are interested in analyzing the H 0 tension problem in the framework of F (R) gravity with the following general form for the action [37]: The first term here is the Einstein-Hilbert action, F inf is assumed to describe the early-time inflation [37][38][39] becoming negligible at late times z < 3000 (only this epoch is visible in our observational data), whereas F DE plays the role of dark energy, dominating at late times, and is the object under study in this paper.Finally, we have added an extra term, F EDE , that behaves as an early dark energy (EDE) term, i.e. mimics an effective cosmological constant at intermediate times but then dilutes along the expansion, helping to suppress some inadequate behaviors during the intermediate phases between the matter-radiation equality and recombination [18].The general field equations for F (R) gravity are obtained by varying the action (1) with respect to the metric g µν , leading to: where R µν and T µν are the Ricci and energy-momentum tensors respectively.By assuming a spatially-flat Friedman-Lemaître-Robertson-Walker (FLRW) space-time with the scale factor a(t), the FLRW equations in F (R) gravity are obtained: The continuity equation ( 4) can be easily solved for dust matter ρ m and radiation ρ r and yields Here a = 1, ρ 0 m and ρ 0 r are the present time values of the scale factor and the matter densities, while we assume the following estimation for the ratio among densities as provided by Planck [1]: The aim of this paper is to explore and compare two clases of F (R) models of the type described by (2).In both cases, we neglect the inflationary term given by F inf and assume some initial conditions that mimic ΛCDM model at large redshifts, namely [34]: Here the index * refers to parameters as given in the ΛCDM model.In particular, Ω * Λ = Λ 3(H * 0 ) 2 and H * 0 is the Hubble constant in the ΛCDM scenario as measured today.However, the late-time evolution for the F (R) models deviates from these initial conditions and consequently from ΛCDM model, such that the above parameters measured today for our models will be different: Nevertheless, these parameters are connected among themselves [28,34]: It is also convenient as shown below, to redefine the Hubble parameter and the Ricci scalar as dimensionless functions: The first model is given by the following exponential function [31,34,36,40]: Note that this exponential model turns out ΛCDM model at the limit β → ∞.Moreover, at large redshifts, the model also recovers ΛCDM as the curvature becomes large enough R >> Λ/β.Hence, physical solutions for this F (R) action tend asymptotically to ΛCDM solutions at large redshifts, such that the above initial conditions (7) results convenient for the equations.By using the dimensionless variables defined in (9), the corresponding system of equations (3) can be rewritten as: This system of equations can be solved by integrating over the independent variable x = log a = − log(z + 1) and assuming the initial conditions (7) at the point x i , where e −βR(xi) ∈ (10 −9 , 10 −7 ) and our model mimics ΛCDM (for more details see Ref. [34]).Then, we confront the model with the observational data by fitting the free parameters and keeping in mind that H(z) = H * 0 E(z) with the true value for the Hubble parameter today being H 0 = H * 0 E(z = 0) and also the relation (8) for the matter density.
Following the same procedure, a second F (R) model with a power-law of the Ricci scalar is analysed, which is described by the gravitational action [37][38][39]: Note that here the so-called early dark energy (EDE) term F EDE is included where α, ℓ, m, n, R 0 are constants, and the curvature scale R 0 corresponds to the Ricci scalar value for the epoch 1000 ≤ z ≤ 3000 (see Ref. [37]).This term can generate a quasi-stable de Sitter stage at R = R 0 as far as n is an odd integer, and ℓ, m are large enough in absence of matter.However, in the presence of matter, the limitations on the parameters ℓ, m, n are connected with the behavior of F RR in Eq. (12).
By considering the model ( 12) without the EDE term (α = 0), the model does not recover purely the ΛCDM model at the limit R → ∞.However, this power-law model may mimic the ΛCDM asymptotic behavior (7) at large curvature 10 ≤ R ≤ 10 10 , whose solutions are free of divergences and singularities.In this approach we can numerically solve the system of equations ( 3) by fixing the initial conditions (7) at large redshifts, corresponding to 1000 ≤ z ≤ 3000.The system of equations (3) for this case (α = 0) yields: By solving these equations, the corresponding solutions show an undesirable oscillatory behavior at large R (see Refs. [37,39]), especially in the most interesting limit δ ≪ 1.An example of these oscillations is depicted in Fig. 1.Such behavior may be controlled via a choice of the initial conditions (they were optimized in the case shown in the figure) but can not be completely suppressed in the framework of this model.Nevertheless, by including the EDE term, these oscillations can be effectively suppressed, that is with a regular evolution for R and H. Nevertheless, one should note that for n ≥ 1 and (or) ℓ ≥ 1 the second derivative of F EDE (R) changes its sign many times close to R = R 0 , so F RR (R) in the denominator of Eq. ( 3) can lead to singularities if α is not small enough.Due to this reason, if the numbers n, m or ℓ are large, the value for α should be small.In this case, F EDE practically does not influence on describing observational data.Nevertheless, for the case m = 1, ℓ = n = 0, the corresponding EDE term becomes: And the denominator F RR (R) in Eq. ( 3) behaves well for such a case (see the top-left panel in Fig. 1) and we can use this term with rather large α to suppress oscillations during the epoch when R ≃ R 0 and later.Figure 1 also depicts the evolution of the Ricci scalar for this choice.
Thus, non-oscillating and non-diverging solutions of the model ( 12) with the EDE term ( 14) can be obtained and confronted with the observational data.

III. OBSERVATIONAL DATA
Here we are interested to confront the models described in the previous section in order to obtain the best fit for the free parameters and particularly the best fit for H 0 when using different data sources.As mentioned above, these observations include: (a) Type Ia supernovae (SNe Ia) data from Pantheon sample [45]; (b) estimates of the Hubble parameter H(z) from cosmic chronometers and line-of-sight BAO and (c) CMB data from Planck 2018 [3,46].
For the SNe Ia we use the Pantheon sample database [45] with N SN = 1048 points and compare the corresponding SNe Ia distance moduli µ obs i at redshift z i from the catalog with their theoretical values by minimizing the χ 2 function: Here C SN is the covariance matrix [45] and p j are the free model parameters, whereas the distance moduli is given by: In the expression (15) the Hubble constant H 0 is considered as a nuisance parameter [34,35,[41][42][43] for SNe Ia data, so can be marginalized and estimations can not be obtained for H 0 from χ 2 SN .However, this provides very important information when fitting the other model parameters.
On the other hand, the Hubble parameter data H(z) are obtained by two different ways of estimation [34,35,[41][42][43][44].The first one is the cosmic chronometers (CC), i.e. estimations of H(z) by using galaxies of different ages ∆t located closely in terms of the redshift ∆z, Here we consider 31 CC H(z) data points given in Ref. [47].In the second method H(z) values are estimated from the baryon acoustic oscillation (BAO) data along the line-of-sight direction.In this paper we use 36 H BAO (z) data points from Refs.[48,49] that can be found in [50].For a particular cosmological model with free parameters p 1 , p 2 , . . ., we calculate the χ 2 function by using the CC H(z) data and the full set CC + H BAO separately, as follows: function Note that H BAO data points are correlated with BAO angular distances, such that are not considered in other analysis (see Refs. [34,35,41]).Nevertheless, here we do not use BAO angular distances, such that we avoid any correlation.
The last source used along the paper is the data from the cosmic microwave background radiation (CMB) that are given by the following observational parameters [46] x where z * = 1089.80± 0.21 is the photon-decoupling redshift [3], while h = H 0 /[100 kms −1 Mpc −1 ], the radiationmatter ratio X r = Ω 0 r /Ω 0 m is given in( 6), and we consider the current baryon fraction Ω 0 b as the nuisance parameter to marginalize over.The corresponding χ 2 function is: where [46] x P l = R P l , ℓ P l A , ω P l b = 1.7428 ± 0.0053, 301.406 ± 0.090, 0.02259 ± 0.00017 , with free amplitude for the lensing power spectrum, from Planck collaboration 2018 data [3].The covariance matrix C CMB = Cij σ i σ j , the expression r s (z * ) and other details are well described in Refs.[35,41] and [46].

IV. RESULTS AND DISCUSSION
Let us now fit the corresponding models parameters through the χ 2 functions as given in ( 15),( 16) and ( 18) for each F (R) model.We consider separately the SNe Ia and H(z) CC (or CC + H BAO ) data, and the same SNe Ia and H(z) data with the CMB contribution ( 18) We follow this procedure as the CMB data (18) with narrow priors (19) produce the most tight limitations on the model parameters, particularly on the density matter parameter due to the factor Ω 0 m in Eq. ( 17) (see Fig. 2).Thus, for the two F (R) models we analyze four different sets of data: Following the way of maximizing the likelihood, we obtain the corresponding distributions and contour plots for the free parameters for both F (R) models.
The exponential model (10) owns four free parameters: However, Ω Λ may be considered as a conditionally free parameter because the χ 2 functions ( 20) and ( 21) have sharp minimums along the line Ω 0 m + Ω Λ ≃ ξ(β), where ξ(β) → 1 in the ΛCDM limit β → ∞. Figure 2 depicts the contours at 1σ (68.27%), 2σ (95.45%) and 3σ (99.73%) for the four data sets given in (22) when considering the exponential model (10).The planes H 0 − Ω 0 m and H 0 − β are obtained by maximizing the likelihood (minimizing the χ 2 ) over the other parameters, while the absolute maximums are described by circles, stars etc.The right panels depict the same contour plots but with additional details, in particular, the second one is re-scaled along the Ω 0 m axis.The corresponding one-parameter distributions shown in the top left panel corresponds to the likelihood function for H 0 after maximizing over the other parameters: Following the same procedure, the fits of the free parameters for the power-law model (12) with the EDE term ( 14) are obtained by the data sets (22).This model has the free parameters: H 0 , Ω 0 m , Ω Λ , γ, δ, α and R 0 .Although the EDE factor α is denoted by α * = log α in Table II.We fix R 0 /(2Λ) = 1.2 • 10 −7 corresponding to the epoch before or near the recombination and work with the remaining 6 parameters.Here Ω Λ can be considered as a conditionally free parameter, because the functions χ 2 j (Ω Λ , . . . ) behave like in the previous exponential model.The 1σ and 2σ contour plots are shown in Fig. 3.In these contour plots and in the likelihood functions (23) we also minimize χ 2 j over the other parameters.
Table II summarizes the results for both F (R) models together with ΛCDM model, where the minimum χ 2 , the best fits of the free parameters and their errors are shown for the different data sets considered here (22).We can evaluate the three models in Table II from the point of view of information criteria depending on the number N p of the free model parameters.In this sense, the Akaike information criterion [51] AIC = min χ 2 tot + 2N p provides a way to compare the goodness of the fits.Then, ΛCDM model with N p = 2 gives a better estimation in comparison with the exponential F (R) model with N p = 4 and the power-law model with N p = 6.On the other hand, the minimum χ 2 j for the exponential model (10) shows a better fit than ΛCDM model for the four sets (22), with the largest difference for the SNe Ia + CC + H BAO data, where additionally the 1σ region for the β parameter does not include ΛCDM model (recall that this is recovered for β → ∞).The other three sets show similar fits, including ΛCDM model within the errors for β.The value of H 0 for the best fit depends on the data set, with its maximum given by H 0 = 70.28 ± 1.6 km s −1 Mpc −1 (SNe Ia + CC + CMB) and its minimum by H 0 = 68.15(SNe Ia + CC + H BAO ).
One can see that for the power-law model (12), which owns 6 free parameters, the absolute minimum min χ 2 j is similar in comparison with the exponential F (R) and the ΛCDM models for all data sets (22), with a slightly better fit for the case with for SNe Ia + CC + H BAO data.This may be explained as follows: only in this case the best fit value for the parameter δ = 0.17 +0.11  −0.095 is large enough for an effective role of the dark energy F (R) term in this model (12).While the other data sets (22), especially for SNe Ia + CC + H BAO data, the best fit for δ is close to zero, but as far as δ → 0 the power-law model (12) tends to the ΛCDM model and the EDE term should be strongly limited.
Concerning the H 0 tension problem, Table II also gathers the estimations for the Hubble constant for both F (R) models and ΛCDM model when confronting the models with the data sets given in (22).By a first look, the best fit values for H 0 are very similar for both the exponential model (10) as for the power-law model (12), such that they are depicted together in the box plots of Fig. 4.Moreover, these estimations are very close to the ΛCDM model predictions for three combinations of the observational data (22), whereas for the case SNe Ia + CC + H BAO , error bars are larger for the F (R) models and the best fit value is shifted to smaller values of H 0 .As shown in Table II, the case SNe Ia + CC + H BAO gives the most disparate results when compared to ΛCDM model, as the errors on the free parameters does not include ΛCDM model within the error bars in both F (R) models.The same applies to the Ω 0 m parameter, which is almost the same value for three of the data sets, and slightly different for the SNe Ia + CC + H BAO set.In addition, the inclusion of CMB data implies much smaller errors on the matter density parameter, as natural due to the factor appearing in (17).
Hence, the results draw an scenario that despite the H 0 tension is not alleviated in the F (R) models, it gives an interesting result when considering the SNe Ia + CC + H BAO data set, as excludes ΛCDM from the 1σ region.

V. CONCLUSIONS
In this paper we have explored two F (R) gravity models, where the cosmological evolution is obtained by solving a dynamical system of equations and then compare to observational data.The models explored along this paper consist of an exponential correction to the Hilbert-Einstein action (10) and the power-law model (12) with the EDE term (14).Keeping in mind their capabilities in alleviating the H 0 tension between the Planck [3] and SH0ES [5] estimations of H 0 , these models were confronted with observational data including SNe Ia, 2 types of H(z) estimations and CMB data, and by combining this data in four different sets in order to analyze the differences and the possible biased introduced by some of the sets on the parameter estimations.
The results are summarized in Table II, which shows the best fits for the Hubble constant H 0 .As discussed along the paper, the values obtained for each model are very similar, with no particular differences among the F (R) models.In the box plot depicted in Fig. 4, the 1σ error for both F (R) gravities is the same, and in The best fit values for H0 (in km s −1 Mpc −1 ) and other parameters for the cosmological scenarios: the exponential F (R) model (10) and the power-law model (12) in comparison with the ΛCDM model.comparison with the one from the ΛCDM model, are also very close for the all data sets (22).The H 0 estimations for the two F (R) models are close to the ΛCDM model.Moreover, the best fits for three of the data sets are achieved within the range where the F (R) models tend asymptotically to ΛCDM, i.e. β → ∞ for the exponential model (12) and δ → 0 for the power-law model (12).Nevertheless, for the SNe Ia + CC + H BAO data set, both F (R) models do not include ΛCDM model within the 1σ region, and the best fits in terms of the minimum value of χ 2 shows a goodness of the fits much better than ΛCDM model, and better than the other data sets.
In addition, as shown in Fig. 4, where the vertical bands refer to H 0 estimations from Planck 2018 release [3], SH0ES (HST) group [5] and the intermediate recent value by CCHP group with the red giant branch (TRGB) method [9], the H 0 tension does not show a better behavior within these F (R) models but the tension problem remains.In particular, while the data sets that excludes CMB data fit well the estimations for H 0 from Planck 2018 and TRGB, it does not when CMB rulers is included.The same applies to the ΛCDM model, as shown in Fig. 4.
Hence, we may conclude that the F (R) models considered in the paper, described by the gravitational actions given in (10) and (12), can not exhaustively explain the tension between the Planck [3] and SH0ES [5] H 0 estimations, but they can alleviate this tension to some extent.The most interesting result lies on the analysis of both F (R) models with SNe Ia + CC + H BAO data set, as excludes the ΛCDM limit from the best fit region, a possible signal of the deviations from ΛCDM and/or of the issue of the data sets.

FIG. 4 :
FIG.4: Box plots for the exponential and power-law F (R) models,together with the ΛCDM model in comparison with Planck18, TRGB20 and SH0ES (R19) H0 estimations.

TABLE I :
Recent estimations of the Hubble constant H0