Spectral Chebyshev Approximation of Cosmic Expansion in $f(R)$ Gravity

TL;DR

We address the stiff, nonlinear background evolution in $f(R)$ gravity and its confrontation with late-time data. We introduce a global spectral Chebyshev collocation method that writes $E(z)$ as $E(z)\approx \sum_{k=0}^N a_k T_k(x(z))$ with $x(z)=2z/z_{\max}-1$ and enforces the dimensionless master equation at Chebyshev–Gauss–Lobatto nodes, converting to a nonlinear algebraic system for the coefficients solved by Powell–hybr. Applied to Hu–Sawicki and Starobinsky models and constrained with cosmic chronometer $H(z)$ data plus Union 3.0 SN, the reconstructed expansion histories agree with observations within $2\sigma$ for $0<z<2$, with best-fit parameters around $\Omega_{m0}\simeq 0.29$, $H_0\simeq 68$ km s$^{-1}$ Mpc$^{-1}$, and $\Lambda_{\mathrm{eff}}/H_0^2\simeq 1.2$–$2.5$. The method provides a precise, computationally efficient framework for testing modified gravity in the precision-data era and offers a foundation for incorporating perturbations and Boltzmann-code integration in future work.

Abstract

We present a numerical framework to study the cosmological background evolution in $f(R)$ gravity by employing a \textit{spectral Chebyshev collocation approach}. Unlike standard integration methods such as Runge--Kutta that often encounter stiffness and accuracy issues, this formulation expands the normalized Hubble function $E(z) = H(z)/H_0$ as a finite Chebyshev series. The modified Friedmann equation is then enforced at selected Chebyshev--Gauss--Lobatto points, converting the original nonlinear differential equation into a system of algebraic relations for the series coefficients. This transformation yields exponentially convergent and numerically stable solutions over the entire redshift domain, $0<z<z_{max}$, eliminating the need for adaptive step-size control. We apply the method to two widely studied $f(R)$ models, Hu--Sawicki and Starobinsky, and perform a combined analysis using cosmic chronometer $H(z)$ data and the Union~3.0 supernova compilation. The reconstructed expansion histories match observations to within $2σ$ over $0 < z < 2$, producing best-fit parameters of approximately $(Ω_{m0}, H_0, Λ_{\mathrm{eff}}) \simeq (0.29, 68, 1.2\text{--}2.5\,H_0^2)$. These results indicates that both models reproduce the observed late-time acceleration while permitting small geometric corrections to $Λ$CDM. Overall, the spectral Chebyshev method provides a precise and computationally efficient framework for probing modified-gravity cosmologies in the precision-data era.

Figures (4)

• Figure 1: Equation Flow for $f(R)$ Cosmology (Dimensionless Form): Flowchart summarizing the dimensionless master equation of $f(R)$ cosmology in redshift space and its related quantities ($\mathcal{R}$, $\mathcal{R}'$, $\tilde{f}$, $\tilde{f}_{\mathcal{R}}$, $\tilde{f}_{\mathcal{RR}}$) for the Hu--Sawicki and Starobinsky models. For a given matter density $\Omega_{m0}$ and model parameters $\theta$, this nonlinear equation can be solved by determining $E(z)$ together with its derivatives $E'(z)$ and $E"(z)$. Hence, solving Eq. \ref{['eq:frfriedmann_dimless']} is the primary aim of background cosmological analysis in $f(R)$ gravity.
• Figure 2: Flowchart of the Chebyshev collocation method for $f(R)$ cosmology, showing the master equation $F(E,E',E";z,\Omega_{m0},\theta)=0$, Ricci scalar definitions, algebraic system setup, solver, and inputs/models required for coding. Refer to Appendix \ref{['app:cheb_solver']} for the stepwise process.
• Figure 3: Corner plots showing the posterior distributions and covariances between cosmological parameters for two $f(R)$ gravity models obtained using Cosmic chronometers $H(z) +$ Unity 3.0 SNe Ia observational datasets. The left panel corresponds to the Hu–Sawicki model having parameters $\Omega_m$, $b$, $\Lambda$, and $H_0$ , while the right panel shows the Starobinsky model having parameters $\Omega_m$, $\mathcal{R}_c$, $\Lambda$, and $H_0$ . Contours represent the 68% and 95% confidence levels for each parameter pair.
• Figure 4: Reconstructed expansion histories $H(z)$ for two $f(R)$ gravity models compared with the standard $\Lambda$CDM scenario. Panel (a) shows the Hu-Sawicki model, while panel (b) corresponds to the Starobinsky model. The solid blue curve represents the reconstructed median evolution of $H(z)$ with the associated $68\%$ credible region shaded in blue, obtained from the MCMC posterior samples. The dashed black line denotes the best-fit $\Lambda$CDM prediction, and the data points with error bars correspond to cosmic chronometer measurements of $H(z)$ at various redshifts.