Can $f(R)$ gravity relieve $H_0$ and $\sigma_8$ tensions?

To investigate whether $f(R)$ gravity can relieve current $H_0$ and $\sigma_8$ tensions, we constrain the Hu-Sawicki $f(R)$ gravity with Planck-2018 cosmic microwave background and redshift space distortions observations. We find that this model fails to relieve both $H_0$ and $\sigma_8$ tensions, and that its two typical parameters $\log_{10}f_{R0}$ and $n$ are insensitive to other cosmological parameters. Combining the cosmic microwave background, baryon acoustic oscillations, Type Ia supernovae, cosmic chronometers with redshift space distortions observations, we give our best constraint $\log_{10}f_{R0}<-6.75$ at the $2\sigma$ confidence level.

constrained specific f (R) models with joint cosmological observations in recent years, there is still a lack of a direct test of the ability to alleviate H 0 and σ 8 tensions for f (R) gravity in light of Planck CMB data. Especially, due to three reasons: (i) the data of Planck-2018 full mission is released; (ii) H 0 tension becomes more serious than before; (iii) richer data from large scale galaxy survey to study DM clustering is gradually obtained, this is an urgent issue needed to be addressed. By implementing numerical analysis, we find that the Hu-Sawicki f (R) gravity cannot reduce H 0 and σ 8 tensions.
This work is organized as follows. In the next section, we introduce the basic equations of f (R) gravity and a specific f (R) model to be investigated in this analysis. In Section III, we display the data and analysis method. In Section IV, the numerical results are presented. The discussions and conclusions are exhibited in the final section.

II. f (R) GRAVITY
To construct a modified theory of gravity, one can introduce some terms such as R 2 , R µν R µν , R µναβ R µναβ , or R n R, when quantum corrections are taken into account. In f (R) gravity, different from the above high-order derivative gravity, the modification is just a function of Ricci scalar R. f (R) gravity was firstly introduced by Buchdahl [15] in 1970 and the readers can find more details in recent reviews [16,17]. The action is written as where f (R), L m and g denote a function of R, the standard matter Lagrangian and the trace of the metric, respectively. By varying Eq.(1), one can obtain the modified Einstein field equation where f R ≡ df /dR denotes an extra scalar degree of freedom, i.e., the so-called scalaron and T µν is energy-momentum tensor. In a spatially flat Friedmann-Robertson-Walker (FRW) universe, the equation of background evolution in f (R) gravity is expressed as where f RR ≡ df R /dR, N ≡ ln a, H is Hubble parameter, a is scale factor and ρ m is matter energy density.
We are also of interests to study the perturbations in f (R) gravity and just consider the linear part here. For sub-horizon modes (k aH) in the quasi-static approximation, the linear growth of matter density perturbations is shown as [18] where Ω m denotes the effective matter density ratio at present. The function X has the following form It is noteworthy that the function X in Eq.(4) induces a scale dependence of linear growth factor δ(k, a) in f (R) gravity, when the growth factor is just a function of scale factor a in GR. In general, a viable f (R) model should be responsible for the inflationary behavior in the very early universe, reproduce the late-time cosmic acceleration, pass the local gravity test, and satisfy the stability conditions. To efficiently investigate cosmological tensions in f (R) gravity, we consider the viable Hu-Sawicki f (R) model (hereafter HS model) [19] in this work and it is given by where µ and n are free parameters characterizing this model. By adopting R µ 2 , the approximate f (R) function shall be written as where R 0 is the present-day value of Ricci scalar and f R0 = f R (R 0 ) = −2Λµ 2 /R 2 0 . For the purpose of constraining this model with data, one should first obtain the evolution of background and perturbation by inserting Eq.(7) into Eqs. (3)(4).
To the best of our knowledge, there are three main methods to confront f (R) gravity with cosmological observations. The first is numerically solving the above equations in a direct way [20][21][22][23][24][25][26][27][28][29]. The second is adopting an approximate framework to obtain the analytic solutions of the above equations and this method, to a large extent, can save computational cost [30][31][32][33][34]. The third one is studying the effects of viable f (R) gravity on the large scale structure formation by using N-body and hydrodynamical simulations [35]. Note that the last method always spend more computational cost and storage space than two previous ones. The constrained 2-dimensional parameter spaces (Ωm, σ8) from the "C" (red) and "R" (blue) datasets are shown for ΛCDM, HS f (R) models with n = 1 and free n, respectively.

III. DATA AND METHOD
Since our aim is to study whether HS f (R) gravity can alleviate the H 0 and σ 8 tensions, first of all, we use the following two main datasets.
CMB: Although the mission of the Planck satellite is completed, its meaning for cosmology and astrophysics is extremely important. It has measured many aspects of formation and evolution of the universe such as matter components, topology and large scale structure effects. Here we shall use updated Planck-2018 CMB temperature and polarization data including likelihoods of temperature at 30 2500 and the low-temperature and polarization likelihoods at 2 29, namely TTTEEE+lowE, and Planck-2018 CMB lensing data [5]. We denote this dataset as "C". The constrained 2-dimensional parameter spaces (log 10 fR0, σ8) from the "C" dataset are shown for HS f (R) models with n = 1 (red), 2 (green), 3 (grey), 4 (orange) and free n (blue), respectively. The constrained 2-dimensional parameter spaces (log 10 fR0, σ8) from the data combination "CBSHR" are shown for HS f (R) models with n = 1 (red) and free n (blue), respectively. The magenta dashed line denotes log 10 fR0 = −6.

RSD:
To study the alleviation of σ 8 tension in f (R) gravity, we adopt the redshift space distortions (RSD) as our reference probe which is sensitive to large scale structure formation. Specifically, we use the so-called "Gold-2018" growth-rate dataset [36]. This dataset is denoted as "R".
Furthermore, to break the parameter degeneracy and give tight constraints on on free parameters of HS model, we also employ the following four probes.
BAO: By measuring the position of these oscillations in the matter power spectrum at different redshifts, the BAO, a standard cosmological ruler, can place constraints on the expansion history of the universe after decoupling and break the parameter degeneracy better. It is unaffected by errors in the nonlinear evolution of the matter density field and other systematic uncertainties. Specifically, we take the 6dFGS sample at the effective redshifts z ef f = 0.106 [37], the SDSS-MGS one at z ef f = 0.15 [38] and the BOSS DR12 dataset at three effective redshifts z ef f = 0.38, 0.51 and 0.61 [39]. Specifically, to constrain the HS f (R) gravity, we use the background quantity D A /r d as a function of scale factor a in the numerical analysis, where D A and r d are angular diameter distance and comoving BAO scale, respectively. To calculate the comoving sound horizon r d , we use the fitting formula given by Ref. [40]. This dataset is identified as "B". The marginalized constraints on the HS f (R) models with n = 1, 2, 3, 4 and free n using the "C" dataset are shown, respectively. For the typical parameter log 10 fR0, we quote 2σ (95%) uncertainties or bounds. The symbol "♦" denotes the parameter that cannot be well constrained by observed data. SNe Ia: SNe Ia, the so-called standard candle, is a powerful distance indicator to study the background evolution of the universe, particularly, the Hubble parameter and EoS of DE. In this analysis, we use the largest SNe Ia "Pantheon" sample today, which integrates the SNe Ia data from the Pan-STARRS1, SNLS, SDSS, low-z and HST surveys and encompasses 1048 spectroscopically confirmed points in the redshift range z ∈ [0.01, 2.3] [41]. In our numerical analysis, we use the full Pantheon sample and marginalize over the absolute magnitude parameter M . We refer to this dataset as "S".
Cosmic Chronometers: As a complementary probe to investigate the late-time evolution of the universe, we also include the cosmic chronometers in our numerical analysis. Specifically, we employ 30 chronometers to constrain the HS model [42]. Hereafter we denote this dataset as "H".
It is worth noting that we take the first method (see Section II), namely numerically solving the background and perturbation equations, to implement constraints on HS f (R) model. In order to obtain the posterior probability density distributions of model parameters, we incorporate the modified equations governing the evolution of background and perturbation of HS f (R) model into the public online packages CAMB and CosmoMC [43,44]. Specifically, we roughly calculate the Hubble expansion rate H(a) and the linear growth factor δ(k, a) at each step of a, and use a interpolating scheme to obtain the solutions H(a) and δ(k, a) with varying a. As a consequence, we can numerically obtain the corresponding cosmological observables to be confronted with data. The latter package can be used for implementing a standard Bayesian analysis via the Markov Chain Monte Carlo (MCMC) method to infer the posterior probability density distributions of parameters. We use the Gelman-Rubin statistic R − 1 = 0.1 as the convergence criterion of MCMC analysis. Meanwhile, to analyze the MCMC chains, we take the public package GetDist [45]. For To investigate comprehensively both H 0 and σ 8 tensions in HS model, we carry out the following numerical analysis. For H 0 tension, respectively, we constrain five models, i.e., n = 1, 2, 3, 4 and free n with the "C" dataset when keeping the typical parameter log 10 f R,0 free. For σ 8 tension, we just present the constraining results of the representative case n = 1 and the general one free n from "C" and "R" datasets, respectively. We also display the comprehensive constraints on two models (n = 1 and free n) by using the data combination "CBSHR". The corresponding χ 2 expressions for all the datasets can be found in Ref. [5].

IV. NUMERICAL RESULTS
For the purpose of studying the alleviation of two important cosmological tensions in the framework of HS f (R) models, our main numerical results are displayed in Fig.1 and Fig.2 and marginalized constraining results are presented The marginalized constraints on the HS f (R) models with n = 1 and free n are shown by using the "R" and "CBSHR" datasets, respectively. Similarly, for the typical parameter log 10 fR0, we quote 2σ (95%) uncertainties. The symbol "♦" denotes the parameter that cannot be well constrained by observed data. To some extent, one can predict its H 0 behavior via Eq.(7). In Fig.1, we have exhibited the constrained 2-dimensional parameter spaces (Ω m , σ 8 ) for five HS models, it is easy to see the large H 0 gap between CMB and HST observations. Only using the CMB data, we conclude that H 0 is insensitive to typical model parameter log 10 f R0 in all five models (see the right panel of Fig.1). To investigate the σ 8 tension, in Fig.2, first of all, we display the constrained Ω m -σ 8 plane for ΛCDM as a reference. Then, we present constrained Ω m -σ 8 planes for the commonly used HS model with n = 1 and for the complete HS model with free n, respectively. We find that the relatively small σ 8 discrepancy in both considered HS scenarios with a little larger parameter spaces (Ω m , σ 8 ) cannot be resolved and is still over 1σ level. This implies that the HS f (R) gravity cannot reduce current H 0 and σ 8 tensions, which is the key result of this work. It is also interesting to study the parameter degeneracy between log 10 f R0 and σ 8 . When only using CMB data, for five HS models, one may find that log 10 f R0 is positively correlated with σ 8 , which indicates that stronger deviations in HS f(R) gravity from GR lead to larger effects of matter clustering (see Fig.3). However, when using the combined dataset CBSHR, this positive correlation disappears and σ 8 seems to be insensitive to log 10 f R0 (see Fig.4). Meanwhile, we are also of interests to study the degeneracy between the additional parameter n and σ 8 , and find that the amplitude of matter clustering is insensitive to n regardless of the use of C or CBSHR datasets (see Fig.5). Furthermore, to study the degeneracies between parameters better, we exhibit the marginalized constraints on HS f (R) models with n = 1 and free n in Fig.6 and Fig.7, and obtain the following conclusions: (i) to a large extent, the parameter spaces are compressed when combining C with BSHR datasets; (ii) in all cases, two typical parameters log 10 f R0 and n are insensitive to other cosmological parameters, which is clarified for the first time in the literature.
We also find that when using only CMB data, the case of free n has the smallest χ 2 = 12958.0 but close to other ones, when using only RSD data, the case of free has a relatively better fitting than HS model with n = 1, and that when using the combined datasets CBSHR, these two cases present almost same χ 2 value. Therefore, we can not easily distinguish these HS f (R) variants from currently statistical analysis.
Moreover, in Tab.I, we can find that the best constraint log 10 f R0 < −4.02 at the 2σ confidence level originates from the case of n = 1 by only using CMB data, while two typical parameters log 10 f R0 and n in the free n case cannot be well constrained (see also Fig.6 and Fig.7). Subsequently, in Tab.II, we find that, when using RSD data alone, constraints on typical parameters of HS models are poor and smaller σ 8 values are obtained, which indicates that this RSD dataset gives a smaller effect of matter clustering at late times than the CMB observation. Interestingly, although the mean value of the constraint H 0 = 75 ± 10 km s −1 Mpc −1 from RSD data is consistent with the HST result, it has a much larger uncertainty. Finally, at the 2σ confidence level, we give our best constraint on the typical parameter log 10 f R0 < −6.75 in the case of n = 1, while log 10 f R0 < −6.60 in the free n case. It is worth noting that we still cannot provide good constraint on n even using the joint dataset CBSHR.
It is noteworthy that there are two interesting and tight constraints from large scale structure observations. In Ref. [21], the authors uses the galaxy clustering ratio, a sensitive probe of the nature of gravity in the cosmological regime, gives f R0 < 4.6 × 10 −5 at the 2σ level. Recently, in Ref. [29], the authors place constraints on chameleon-f (R) gravity from galaxy rotation curves and find that f (R) models within the range −7.5 < log 10 f R0 < −6.5 seem to be favored with respect to ΛCDM. Interestingly, our best constraint just lies in this range and this may give a clue of the correct living range for the HS f (R) gravity.

V. DISCUSSIONS AND CONCLUSIONS
Recently, the H 0 and σ 8 tensions under the standard cosmological paradigm have re-activated a wide variety of alternative cosmological models. However, all the time, there is a lack of direct tests of f (R) gravity in resolving both tensions. To address this urgent issue, we confront the popular HS f (R) gravity with current observations. By testing five specific HS f (R) models with observational datasets, we obtain two main conclusions: (i) HS f (R) gravity cannot resolve both H 0 and σ 8 tensions; (i) the typical parameters log 10 f R0 and n are insensitive to other cosmological parameters. Meanwhile, in the HS f (R) model with n = 1, we give our best constraint log 10 f R0 < −6.75 at the 2σ confidence level.
It is noteworthy that a coupling between matter and geometry in the framework of f (R) gravity may help resolve these tensions, and that other f (R) gravity models may relieve both discrepancies much better than the considered HS f (R) one. We expect that future high-precision CMB and SNe Ia observations and independent probes such as gravitational waves could help reduce or even solve these intractable cosmological tensions.