Alleviating the $H_0$ tension in Rastall gravity

TL;DR

This paper explores the Rastall $\\Lambda$CDM model as a way to reduce the $H_0$ tension in cosmology, finding it alleviates the discrepancy with local measurements but does not resolve the $\sigma_8$ tension.

Contribution

It introduces and constrains the Rastall gravity framework within cosmology, demonstrating its potential to mitigate the $H_0$ tension compared to standard models.

Findings

The Rastall $\\Lambda$CDM model reduces the $H_0$ tension to about 2 sigma.

Closed geometry in the model further alleviates the $H_0$ discrepancy.

The $\sigma_8$ tension remains unresolved in this framework.

Abstract

The tension between measurements of the Hubble constant $H_0$ locally and the value inferred from Planck satellite have provided a strong motivation to explore theoretical frameworks beyond the standard cosmological model ($\Lambda$CDM). To this aim, we investigate the $H_0$ tension in the Rastall $\Lambda$CDM model (R-$\Lambda$CDM), formulated within the Rastall gravity framework as a phenomenological extension of General Relativity featuring a non-conserved energy-momentum tensor. To constrain the cosmological parameters, we employ a comprehensive combined dataset including Big Bang Nucleosynthesis (BBN), Hubble parameter $H(z)$, Type Ia supernovae (SNe Ia), Baryon Acoustic Oscillations (BAO), Cosmic Microwave Background (CMB) anisotropies, and the growth rate of cosmic structures $f(z)\sigma_8(z)$. The estimation of model parameters is performed via Markov Chain Monte Carlo (MCMC) sampling within a Bayesian statistical framework. Our results show that the R-$\Lambda$CDM model, especially in closed geometry, reduces the $H_0$ tension, bringing the discrepancy with the local SH0ES measurements down to $1.97\sigma$ to $2.18\sigma$. However, the tension related to the $\sigma_8$ parameter remains unresolved.

Tables

Table 1. Table 4: The 68% confidence limits of flat Λ \Lambda CDM model cosmological parameters (CPs) from D 1 , D 2 , D 3 , and D 4 datasets. H 0 has units of km s -1  Mpc -1 . Here Ω m 0 = Ω d ​ m 0 + Ω b 0 \Omega_{m_{0}}=\Omega_{dm_{0}}+\Omega_{b_{0}} .

CPs	D1	D2	D3	D4
${\Omega }_{b_0}$	${0.0440}_{-0.0193}^{+0.0195}$	${0.0588}_{-0.0014}^{+0.0013}$	$0.0494\pm 0.0003$	${0.0488}_{-0.0007}^{+0.0007}$
${\Omega }_{d{m}_0}$	${0.2447}_{-0.0225}^{+0.0229}$	${0.2363}_{-0.0075}^{+0.0078}$	$0.2719\pm 0.0036$	${0.2675}_{-0.0076}^{+0.0079}$
${\Omega }_{{\Lambda }_0}$	${0.7110}_{-0.0128}^{+0.0126}$	${0.7048}_{-0.0072}^{+0.0070}$	${0.6786}_{-0.0039}^{+0.0038}$	${0.6835}_{-0.0086}^{+0.0083}$
${\Omega }_{r_0}$	$\left({8.2172}_{-0.1842}^{+0.1855}\right)\times {10}^{-5}$	$\left({8.2776}_{-0.1651}^{+0.1727}\right)\times {10}^{-5}$	$\left({9.2088}_{-0.0754}^{+0.0772}\right)\times {10}^{-5}$	$\left({9.1236}_{-0.1607}^{+0.1660}\right)\times {10}^{-5}$
$H_0$	${71.30}_{-0.79}^{+0.81}$	${71.04}_{-0.73}^{+0.72}$	$67.35\pm 0.28$	${67.66}_{-0.61}^{+0.60}$
${\sigma }_8$	${0.8177}_{-0.0237}^{+0.0243}$	${0.8150}_{-0.0229}^{+0.0227}$	${0.8030}_{-0.0224}^{+0.0222}$	$-$
$\mathrm{AIC}$	$1581.53132$	$1593.33913$	$1642.92754$	$6.00022$
$\mathrm{BIC}$	$1608.91651$	$1620.76378$	$1670.35780$	$3.29606$

Table 2. Table 5: The 68% confidence limits of open Λ \Lambda CDM model cosmological parameters (CPs) from D 1 , D 2 , and D 3 datasets. H 0 has units of km s -1 Mpc -1 . Here Ω m 0 = Ω d ​ m 0 + Ω b 0 \Omega_{m_{0}}=\Omega_{dm_{0}}+\Omega_{b_{0}} .

CPs	D1	D2	D3
${\Omega }_{b_0}$	${0.0441}_{-0.0192}^{+0.0195}$	${0.0625}_{-0.0027}^{+0.0035}$	$0.0492\pm 0.0003$
${\Omega }_{d{m}_0}$	${0.2116}_{-0.0285}^{+0.0284}$	${0.2244}_{-0.0122}^{+0.0106}$	${0.2702}_{-0.0037}^{+0.0036}$
${\Omega }_{{\Lambda }_0}$	${0.6455}_{-0.0413}^{+0.0400}$	${0.6789}_{-0.0231}^{+0.0178}$	${0.6802}_{-0.0039}^{+0.0040}$
${\Omega }_{r_0}$	$\left({8.3036}_{-0.1938}^{+0.1974}\right)\times {10}^{-5}$	$\left({8.3456}_{-0.1749}^{+0.1788}\right)\times {10}^{-5}$	$\left({9.1536}_{-0.0793}^{+0.0775}\right)\times {10}^{-5}$
${\Omega }_{k_0}$	${0.0983}_{-0.0578}^{+0.0591}$	${0.0333}_{-0.0218}^{+0.0289}$	${0.0002}_{-0.0002}^{+0.0004}$
$H_0$	${70.92}_{-0.83}^{+0.84}$	${70.75}_{-0.75}^{+0.75}$	${67.55}_{-0.28}^{+0.29}$
${\sigma }_8$	${0.8655}_{-0.0361}^{+0.0395}$	${0.8288}_{-0.0250}^{+0.0259}$	$0.8048\pm 0.0224$
$\mathrm{AIC}$	$1582.21109$	$1595.00379$	$1647.55658$
$\mathrm{BIC}$	$1615.07332$	$1627.91337$	$1680.47290$

Table 3. Table 6: The 68% confidence limits of closed Λ \Lambda CDM model cosmological parameters (CPs) from D 1 , D 2 , and D 3 datasets. H 0 has units of km s -1 Mpc -1 . Here Ω m 0 = Ω d ​ m 0 + Ω b 0 \Omega_{m_{0}}=\Omega_{dm_{0}}+\Omega_{b_{0}} .

CPs	D1	D2	D3
${\Omega }_{b_0}$	${0.0436}_{-0.0190}^{+0.0194}$	${0.0569}_{-0.0023}^{+0.0019}$	$0.0505\pm 0.0005$
${\Omega }_{d{m}_0}$	${0.2579}_{-0.0247}^{+0.0256}$	${0.2426}_{-0.0089}^{+0.0099}$	${0.2634}_{-0.0039}^{+0.0040}$
${\Omega }_{{\Lambda }_0}$	${0.7332}_{-0.0202}^{+0.0266}$	${0.7179}_{-0.0112}^{+0.0156}$	$0.6907\pm 0.0046$
${\Omega }_{r_0}$	$\left({8.1853}_{-0.1886}^{+0.1941}\right)\times {10}^{-5}$	$\left({8.2444}_{-0.1677}^{+0.1729}\right)\times {10}^{-5}$	$\left(9.2550\pm 0.0814\right)\times {10}^{-5}$
${\Omega }_{k_0}$	$-{0.0306}_{-0.0393}^{+0.0221}$	$-{0.0156}_{-0.0199}^{+0.0113}$	$-0.0047\pm 0.0013$
$H_0$	${71.43}_{-0.83}^{+0.84}$	$71.18\pm 0.74$	$67.18\pm 0.30$
${\sigma }_8$	${0.8032}_{-0.0268}^{+0.0257}$	${0.8084}_{-0.0236}^{+0.0232}$	${0.8045}_{-0.0223}^{+0.0222}$
$\mathrm{AIC}$	$1583.71522$	$1595.45415$	$1634.50770$
$\mathrm{BIC}$	$1616.57745$	$1628.36373$	$1667.42402$

Table 4. Table 7: Cosmological Parameter Differences (CPD) expressed in units of σ \sigma for the flat, open, and closed Λ \Lambda CDM models. The deviations are computed between D 1 and D 2 , and between D 1 and D 3 , using symmetrized uncertainties. The Planck-only dataset D 4 is treated separately as an early-Universe reference baseline. Values are rounded to one decimal place.

CPD	$\Lambda$CDM (Flat)		$\Lambda$CDM (Open)		$\Lambda$CDM (Closed)
	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>2</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>3</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>2</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>3</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>2</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>3</sub>}\right|$
${\Omega }_{d{m}_0}$	0.4$\sigma$	1.2$\sigma$	0.4$\sigma$	2.0$\sigma$	0.5$\sigma$	0.2$\sigma$
${\Omega }_{b_0}$	0.8$\sigma$	0.3$\sigma$	0.9$\sigma$	0.3$\sigma$	0.7$\sigma$	0.4$\sigma$
${\Omega }_{k_0}$	–	–	1.0$\sigma$	1.7$\sigma$	0.4$\sigma$	0.8$\sigma$
${\Omega }_{{\Lambda }_0}$	0.4$\sigma$	2.4$\sigma$	0.7$\sigma$	0.9$\sigma$	0.6$\sigma$	1.7$\sigma$
${\Omega }_{r_0}$	0.2$\sigma$	5.0$\sigma$	0.2$\sigma$	4.0$\sigma$	0.2$\sigma$	4.8$\sigma$
$H_0$	0.2$\sigma$	4.6$\sigma$	0.2$\sigma$	3.8$\sigma$	0.2$\sigma$	4.8$\sigma$
${\sigma }_8$	0.1$\sigma$	0.5$\sigma$	0.7$\sigma$	1.4$\sigma$	0.2$\sigma$	0.0$\sigma$

Table 5. Table 8: The 68% confidence limits of flat R- Λ \Lambda CDM cosmological parameters from D 1 , D 2 , D 3 , and D 4 datasets. H 0 H_{0} is in km s -1  Mpc -1 . In Rastall gravity we report the effective density parameters Ω i eff ≡ Ω ~ i \Omega_{i}^{\rm eff}\equiv\tilde{\Omega}_{i} that govern the background expansion.

CPs	D1	D2	D3	D4
${\Omega }_{b0}^{\mathrm{eff}}$	${0.03708}_{-0.01600}^{+0.01676}$	${0.05494}_{-0.01278}^{+0.01326}$	$0.04856\pm 0.00039$	${0.05018}_{-0.00429}^{+0.00191}$
${\Omega }_{dm0}^{\mathrm{eff}}$	${0.2749}_{-0.0232}^{+0.0216}$	${0.2378}_{-0.0092}^{+0.0094}$	${0.2749}_{-0.0038}^{+0.0037}$	${0.29934}_{-0.07910}^{+0.08302}$
${\Omega }_{\Lambda 0}^{\mathrm{eff}}$	${0.6880}_{-0.0150}^{+0.0153}$	${0.7070}_{-0.0097}^{+0.0099}$	${0.6765}_{-0.0039}^{+0.0040}$	${0.65038}_{-0.08507}^{+0.08340}$
${\Omega }_{r0}^{\mathrm{eff}}$	$\left(7.16\pm 0.43\right)\times {10}^{-5}$	$\left(8.30\pm 0.18\right)\times {10}^{-5}$	$\left(8.94\pm 0.10\right)\times {10}^{-5}$	$\left({9.33}_{-0.74}^{+0.35}\right)\times {10}^{-5}$
$ϵ$	$-{0.06437}_{-0.02447}^{+0.03287}$	${0.00189}_{-0.00662}^{+0.00615}$	${0.00203}_{-0.00049}^{+0.00048}$	${0.002265}_{-0.006460}^{+0.004823}$
$H_0$	${71.32}_{-0.81}^{+0.84}$	$71.06\pm 0.75$	${68.50}_{-0.41}^{+0.40}$	${67.130}_{-1.041}^{+2.282}$
${\sigma }_8$	$0.8232\pm 0.0240$	${0.8147}_{-0.0228}^{+0.0232}$	${0.8002}_{-0.0225}^{+0.0226}$	$-$
$\mathrm{AIC}$	$1580.01154$	$1595.12210$	$1628.04364$	$8.00019$
$\mathrm{BIC}$	$1612.87377$	$1628.03168$	$1660.95996$	$4.39464$

Table 6. Table 9: Comparison of AIC and BIC values between flat-universe R- Λ \Lambda CDM and Λ \Lambda CDM models for datasets D 1 , D 2 , and D 3 . The differences are defined as Δ ​ AIC = AIC R − AIC Λ \Delta\mathrm{AIC}=\mathrm{AIC}_{\rm R}-\mathrm{AIC}_{\Lambda} and similarly for Δ ​ BIC \Delta\mathrm{BIC} . Negative values indicate a preference for the R- Λ \Lambda CDM model. The comparison is performed within each dataset independently, as AIC and BIC are not directly comparable across datasets with different numbers of data points.

Dataset	AIC			BIC
	$\Lambda$CDM	R-$\Lambda$CDM	$\Delta$AIC	$\Lambda$CDM	R-$\Lambda$CDM	$\Delta$BIC
D1	1581.53	1580.01	$-1.52$	1608.92	1612.87	$+3.96$
D2	1593.34	1595.12	$+1.78$	1620.76	1628.03	$+7.27$
D3	1642.93	1628.04	$-14.88$	1670.36	1660.96	$-9.40$

Table 7. Table 10: The 68% confidence limits of open R- Λ \Lambda CDM cosmological parameters from D 1 , D 2 , and D 3 datasets. H 0 H_{0} is in km s -1  Mpc -1 . In Rastall gravity we report effective density parameters Ω i eff ≡ Ω ~ i \Omega_{i}^{\rm eff}\equiv\tilde{\Omega}_{i} that govern the background expansion; Ω k ​ 0 \Omega_{k0} is not rescaled.

CPs	D1	D2	D3
${\Omega }_{b0}^{\mathrm{eff}}$	${0.03805}_{-0.01656}^{+0.01746}$	${0.05243}_{-0.01240}^{+0.01400}$	${0.04772}_{-0.00056}^{+0.00050}$
${\Omega }_{dm0}^{\mathrm{eff}}$	${0.23881}_{-0.03577}^{+0.03134}$	${0.22733}_{-0.01236}^{+0.01153}$	${0.27424}_{-0.00398}^{+0.00413}$
${\Omega }_{\Lambda 0}^{\mathrm{eff}}$	${0.63181}_{-0.04785}^{+0.03810}$	${0.67908}_{-0.02389}^{+0.01962}$	${0.67675}_{-0.00469}^{+0.00438}$
${\Omega }_{r0}^{\mathrm{eff}}$	$\left({7.29}_{-0.47}^{+0.63}\right)\times {10}^{-5}$	$\left({8.42}_{-0.19}^{+0.20}\right)\times {10}^{-5}$	$\left(8.70\pm 0.12\right)\times {10}^{-5}$
${\Omega }_{k0}$	${0.08911}_{-0.05727}^{+0.07347}$	${0.03931}_{-0.02515}^{+0.03211}$	${0.00108}_{-0.00076}^{+0.00120}$
$H_0$	$70.96\pm 0.84$	$70.79\pm 0.76$	${69.50}_{-0.49}^{+0.52}$
$ϵ$	$-{0.06077}_{-0.02718}^{+0.03640}$	${0.00526}_{-0.00697}^{+0.00604}$	${0.00296}_{-0.00054}^{+0.00053}$
${\sigma }_8$	${0.8685}_{-0.0374}^{+0.0473}$	${0.8307}_{-0.0254}^{+0.0263}$	${0.8009}_{-0.0226}^{+0.0224}$
$\mathrm{AIC}$	$1581.62349$	$1596.38812$	$1621.95742$
$\mathrm{BIC}$	$1619.96276$	$1634.78263$	$1660.35979$

Table 8. Table 11: Comparison of AIC and BIC values between open-universe R- Λ \Lambda CDM and Λ \Lambda CDM models for datasets D 1 , D 2 , and D 3 . The differences are defined as Δ \Delta AIC = AIC R − AIC Λ =\mathrm{AIC}_{\rm R}-\mathrm{AIC}_{\Lambda} and similarly for Δ \Delta BIC.

Dataset	AIC			BIC
	$\Lambda$CDM	R-$\Lambda$CDM	$\Delta$AIC	$\Lambda$CDM	R-$\Lambda$CDM	$\Delta$BIC
D1	1582.21	1581.62	$-0.59$	1615.07	1619.96	$+4.89$
D2	1595.00	1596.39	$+1.38$	1627.91	1634.78	$+6.87$
D3	1647.56	1621.96	$-25.60$	1680.47	1660.36	$-20.11$

Table 9. Table 12: The 68% confidence limits of closed R- Λ \Lambda CDM cosmological parameters from D 1 , D 2 , and D 3 datasets. H 0 H_{0} is in km s -1  Mpc -1 . In Rastall gravity we report effective density parameters Ω i eff ≡ Ω ~ i \Omega_{i}^{\rm eff}\equiv\tilde{\Omega}_{i} that govern the background expansion; Ω k ​ 0 \Omega_{k0} is not rescaled.

CPs	D1	D2	D3
${\Omega }_{b0}^{\mathrm{eff}}$	${0.03646}_{-0.01589}^{+0.01664}$	${0.05684}_{-0.01328}^{+0.01231}$	${0.04870}_{-0.00054}^{+0.00068}$
${\Omega }_{dm0}^{\mathrm{eff}}$	${0.29318}_{-0.02649}^{+0.02764}$	${0.24318}_{-0.01010}^{+0.01044}$	${0.26989}_{-0.00437}^{+0.00407}$
${\Omega }_{\Lambda 0}^{\mathrm{eff}}$	${0.71445}_{-0.02386}^{+0.03240}$	${0.71924}_{-0.01324}^{+0.01589}$	${0.68274}_{-0.00452}^{+0.00498}$
${\Omega }_{r0}^{\mathrm{eff}}$	$\left({7.09}_{-0.39}^{+0.54}\right)\times {10}^{-5}$	$\left(8.25\pm 0.18\right)\times {10}^{-5}$	$\left(8.85\pm 0.13\right)\times {10}^{-5}$
${\Omega }_{k0}$	$-{0.03897}_{-0.04994}^{+0.02832}$	$-{0.01614}_{-0.02057}^{+0.01182}$	$-{0.00129}_{-0.00148}^{+0.00091}$
$H_0$	${71.48}_{-0.83}^{+0.84}$	${71.15}_{-0.74}^{+0.76}$	${68.86}_{-0.54}^{+0.53}$
$ϵ$	$-{0.06681}_{-0.02279}^{+0.03127}$	$-{0.00009}_{-0.00645}^{+0.00660}$	${0.00230}_{-0.00064}^{+0.00057}$
${\sigma }_8$	${0.8030}_{-0.0290}^{+0.0276}$	${0.8078}_{-0.0237}^{+0.0231}$	${0.8025}_{-0.0230}^{+0.0231}$
$\mathrm{AIC}$	$1582.06061$	$1597.27448$	$1621.92194$
$\mathrm{BIC}$	$1620.39988$	$1635.66899$	$1660.32431$

Table 10. Table 13: Comparison of AIC and BIC values between closed-universe R- Λ \Lambda CDM and Λ \Lambda CDM models for datasets D 1 , D 2 , and D 3 . The differences are defined as Δ \Delta AIC = AIC R − AIC Λ =\mathrm{AIC}_{\rm R}-\mathrm{AIC}_{\Lambda} and similarly for Δ \Delta BIC.

Dataset	AIC			BIC
	$\Lambda$CDM	R-$\Lambda$CDM	$\Delta$AIC	$\Lambda$CDM	R-$\Lambda$CDM	$\Delta$BIC
D1	1583.72	1582.06	$-1.65$	1616.58	1620.40	$+3.82$
D2	1595.45	1597.27	$+1.82$	1628.36	1635.67	$+7.31$
D3	1634.51	1621.92	$-12.59$	1667.42	1660.32	$-7.10$

Table 11. Table 14: Updated cosmological parameter differences (CPD) expressed in units of σ \sigma for the R- Λ \Lambda CDM model in flat and non-flat geometries. The CPDs are computed using symmetrized 1 ​ σ 1\sigma uncertainties, σ = ( σ + + σ − ) / 2 \sigma=(\sigma_{+}+\sigma_{-})/2 , and σ comb = σ 1 2 + σ 2 2 \sigma_{\rm comb}=\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}} . In Rastall gravity, all density parameters are the effective ones Ω i eff \Omega_{i}^{\rm eff} governing the background expansion, while Ω k ​ 0 \Omega_{k0} is not rescaled.

CPD	Flat		Open		Closed
	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>2</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>3</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>2</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>3</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>2</sub>}\right|$	$\left|\text{D}\text{<sub>1</sub>}-\text{D}\text{<sub>3</sub>}\right|$
${\Omega }_{dm0}^{\mathrm{eff}}$	$1.53\sigma$		$0.32\sigma$	$1.05\sigma$	$1.73\sigma$	$0.85\sigma$
${\Omega }_{b0}^{\mathrm{eff}}$	$0.85\sigma$	$0.70\sigma$	$0.67\sigma$	$0.57\sigma$	$0.98\sigma$	$0.75\sigma$
${\Omega }_{k0}$	–	–	$0.70\sigma$	$1.35\sigma$	$0.54\sigma$	$0.96\sigma$
${\Omega }_{\Lambda 0}^{\mathrm{eff}}$	$1.05\sigma$	$0.73\sigma$	$0.98\sigma$	$1.04\sigma$	$0.15\sigma$	$1.11\sigma$
${\Omega }_{r0}^{\mathrm{eff}}$	$2.45\sigma$	$4.03\sigma$	$1.94\sigma$	$2.50\sigma$	$2.33\sigma$	$3.65\sigma$
$ϵ$	$2.26\sigma$	$2.32\sigma$	$2.03\sigma$	$2.00\sigma$	$2.40\sigma$	$2.56\sigma$
$H_0$	$0.23\sigma$	$3.07\sigma$	$0.15\sigma$	$1.49\sigma$	$0.29\sigma$	$2.64\sigma$
${\sigma }_8$	$0.26\sigma$	$0.70\sigma$	$0.76\sigma$	$1.41\sigma$	$0.13\sigma$	$0.01\sigma$

Table 12. Table 15: Summary of the inferred Hubble constant H 0 H_{0} (km s -1  Mpc -1 ) for all investigated cosmological models and dataset combinations. Late-time datasets D 1 and D 2 generally favor higher values of H 0 H_{0} , while the inclusion of CMB information in D 3 and D 4 shifts the constraints toward the Planck-preferred region.

Model	D1	D2	D3	D4
Flat $\Lambda$CDM	${71.30}_{-0.79}^{+0.81}$	${71.04}_{-0.73}^{+0.72}$	$67.35\pm 0.28$	${67.66}_{-0.61}^{+0.60}$
Open $\Lambda$CDM	${70.92}_{-0.83}^{+0.84}$	${70.75}_{-0.75}^{+0.75}$	${67.55}_{-0.28}^{+0.29}$	–
Closed $\Lambda$CDM	${71.43}_{-0.83}^{+0.84}$	$71.18\pm 0.74$	$67.18\pm 0.30$	–
Flat R-$\Lambda$CDM	${71.32}_{-0.81}^{+0.84}$	$71.06\pm 0.75$	${68.50}_{-0.41}^{+0.40}$	${67.13}_{-1.04}^{+2.28}$
Open R-$\Lambda$CDM	$70.96\pm 0.84$	$70.79\pm 0.76$	${69.50}_{-0.49}^{+0.52}$	–
Closed R-$\Lambda$CDM	${71.48}_{-0.83}^{+0.84}$	${71.15}_{-0.74}^{+0.76}$	${68.86}_{-0.54}^{+0.53}$	–

Table 13. Table 16: Comparison of the inferred H 0 H_{0} values in Λ \Lambda CDM and R- Λ \Lambda CDM for the CMB-inclusive D 3 dataset.

Geometry	$H_0$ in $\Lambda$CDM	$H_0$ in R-$\Lambda$CDM	$ϵ$ in R-$\Lambda$CDM
Flat	$67.35\pm 0.28$	${68.50}_{-0.41}^{+0.40}$	${0.00203}_{-0.00049}^{+0.00048}$
Open	${67.55}_{-0.28}^{+0.29}$	${69.50}_{-0.49}^{+0.52}$	${0.00296}_{-0.00054}^{+0.00053}$
Closed	$67.18\pm 0.30$	${68.86}_{-0.54}^{+0.53}$	${0.00230}_{-0.00064}^{+0.00057}$

Table 14. Table 17: Summary of H 0 H_{0} tension (in units of σ \sigma ) relative to SH0ES ( 73.2 ± 1.3 73.2\pm 1.3 ) for D 1 , D 2 and D 3 and relative to the Planck 2018 ( 67.4 ± 0.5 67.4\pm 0.5 ) for D 3 .

Model	D1 (SH0ES)	D2 (SH0ES)	D3 (Planck)	D3 (SH0ES)
Flat $\Lambda$CDM	$1.23\sigma$	$1.44\sigma$	$0.09\sigma$	$4.40\sigma$
Open $\Lambda$CDM	$1.48\sigma$	$1.63\sigma$	$0.26\sigma$	$4.25\sigma$
Closed $\Lambda$CDM	$1.15\sigma$	$1.35\sigma$	$0.38\sigma$	$4.51\sigma$
Flat R-$\Lambda$CDM	$1.22\sigma$	$1.43\sigma$	$1.70\sigma$	$3.45\sigma$
Open R-$\Lambda$CDM	$1.45\sigma$	$1.60\sigma$	$2.95\sigma$	$2.65\sigma$
Closed R-$\Lambda$CDM	$1.11\sigma$	$1.37\sigma$	$1.99\sigma$	$3.09\sigma$

Equations

$$R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }=8\pi G{T}_{\mu \nu },$$

$$T_{\mu \nu }=\sum _i{T}_{\mu \nu }^i,$$

$$T_i^{\mu \nu }=({\rho }_i+p_i)u^{\mu }{u}^{\nu }+p_i{g}^{\mu \nu },$$

$$H^2+\frac{k}{a^2}=\frac{1}{3}\sum _i{\rho }_i,$$

$$\frac{\ddot{a}}{a}=-\frac{1}{6}\sum _i({\rho }_i+3p_i).$$

$$\dot{{\rho }_i}=-3H(p_i+{\rho }_i)=-3H(1+{\omega }_i){\rho }_i.$$

$${\rho }_i(a)={\rho }_{i0}{\left(\frac{a_0}{a}\right)}^{3(1+{\omega }_i)},$$

$${\Omega }_{m_0}\equiv \frac{{\rho }_{m_0}}{3H_0^2},{\Omega }_{{\Lambda }_0}\equiv \frac{{\rho }_{{\Lambda }_0}}{3H_0^2},{\Omega }_{r_0}\equiv \frac{{\rho }_{r_0}}{3H_0^2},{\Omega }_{k_0}\equiv -\frac{k}{a_0^2{H}_0^2}.$$

$$E(z;\mathbf{p})\equiv \frac{H(z;\mathbf{p})}{H_0}=[{\Omega }_{m_0}(1+z{)}^3+{\Omega }_{r_0}(1+z{)}^4+{\Omega }_{k_0}(1+z{)}^2+{\Omega }_{{\Lambda }_0}{]}^{1/2}$$

$$G_{\mu \nu }=\kappa \left(T_{\mu \nu }+{\hat{T}}_{\mu \nu }\right),$$

$${\hat{T}}_{\mu \nu }=-ϵR{g}_{\mu \nu }.$$

$$\kappa =8\pi G\left(\frac{1-4ϵ\kappa }{1-6ϵ\kappa }\right).$$

$$6ϵ{\kappa }^2-(1+4ϵ)\kappa +1=0.$$

$$ϵ\le \frac{2-\sqrt{3}}{4}\simeq 0.067\text{or}ϵ\ge \frac{2+\sqrt{3}}{4}.$$

$$|ϵ\kappa |\ll 1,$$

$$ϵ\in (-0.10,0.07).$$

$$\kappa \simeq 1+2ϵ.$$

$$G_{\mu \nu }=T_{\mu \nu }+ϵ\left[2T_{\mu \nu }-(1+2ϵ)R{g}_{\mu \nu }\right].$$

$$R=-\frac{\kappa }{1-4ϵ\kappa }T.$$

$$R=-(1+2ϵ)\left(1+4ϵ+24{ϵ}^2+\cdots \right)T,$$

$$G_{\nu }^{\mu }=T_{\nu }^{\mu }+ϵ\left[2T_{\nu }^{\mu }+(1+8ϵ)g_{\nu }^{\mu }T\right].$$

$${T_{(i)}}_{\nu }^{\mu }=\mathrm{diag}(-{\rho }_i,p_i,p_i,p_i),$$

$$H^2+\frac{k}{a^2}=\frac{1}{3}\sum _i{\rho }_i\left[1+3ϵ(1-{\omega }_i)+8{ϵ}^2(1-3{\omega }_i)\right].$$

$$\frac{\ddot{a}}{a}=\frac{\left[1-6ϵ(1+2ϵ)\right]\frac{1}{3}{\sum }_i{\rho }_i\left[1+3ϵ(1-{\omega }_i)+8{ϵ}^2(1-3{\omega }_i)\right]+(1+2ϵ){\sum }_i{\omega }_i{\rho }_i}{6ϵ(1+2ϵ)-2}.$$

$$\frac{\ddot{a}}{a}=-\frac{1}{6}\sum _i{\rho }_i\left[1+3{\omega }_i+12ϵ{\omega }_i+16{ϵ}^2(3{\omega }_i-1)\right].$$

$$\frac{\ddot{a}}{a}=-\frac{1}{6}\sum _i({\rho }_i+3p_i).$$

$${\dot{\rho }}_i=-3\frac{\dot{a}}{a}\left[\frac{1+{\omega }_i}{1+ϵ(1+6ϵ)(3{\omega }_i-1)}\right]{\rho }_i,$$

$${\rho }_i={\rho }_{i0}{\left(\frac{a_0}{a}\right)}^{\frac{3(1+{\omega }_i)}{1+ϵ(1+6ϵ)(3{\omega }_i-1)}}.$$

$${\tilde{\Omega }}_i={\Omega }_i\left[1+3ϵ(1-{\omega }_i)+8{ϵ}^2(1-3{\omega }_i)\right],$$

$$\frac{H(z;\mathbf{p})}{H_0}=[{\tilde{\Omega }}_{m_0}(1+z{)}^{\frac{3}{1-ϵ(1+6ϵ)}}+{\tilde{\Omega }}_{r_0}(1+z{)}^4+{\Omega }_{k_0}(1+z{)}^2+{\tilde{\Omega }}_{{\Lambda }_0}{]}^{1/2}.$$

Full text

We analyze three primary dataset combinations: D1 (late-time probes: SN + $H(z)$ + $f\sigma_{8}$), D2 (late-time probes combined with DESI DR2 BAO and BBN), and D3 (late-time probes combined with BAO and Planck 2018 CMB distance priors). Parameters are constrained via Markov Chain Monte Carlo sampling, and tensions with SH0ES ($73.2\pm 1.3~{\rm km\,s^{-1}\,Mpc^{-1}}$) and Planck are expressed in units of the combined uncertainty. In addition, we include a Planck-only configuration (D4) as a reference baseline to isolate early-Universe constraints on $H_{0}$. Within $\Lambda$CDM, D1 and D2 yield $H_{0}\simeq 70.75$-$71.43~{\rm km\,s^{-1}\,Mpc^{-1}}$, reducing the SH0ES discrepancy to $1.11\sigma$–$1.63\sigma$ while maintaining a $3.62\sigma$-$4.23\sigma$ tension with Planck. Including CMB distance priors (D3) shifts the result to $H_{0}\simeq 67.18$–$67.55~{\rm km\,s^{-1}\,Mpc^{-1}}$, consistent with Planck at $0.09\sigma$-$0.38\sigma$ but increasing the SH0ES discrepancy to $4.25\sigma$-$4.51\sigma$.

A similar dataset-driven behavior is found in R-$\Lambda$CDM. For D1 and D2, we find $H_{0}\simeq 70.79$–$71.48~{\rm km\,s^{-1}\,Mpc^{-1}}$, with SH0ES tensions of $1.11\sigma$-$1.60\sigma$ and Planck discrepancies of $3.64\sigma$-$4.19\sigma$. For D3, the inferred value becomes $H_{0}\simeq 68.50$–$69.50~{\rm km\,s^{-1}\,Mpc^{-1}}$, reducing the Planck tension to $1.70\sigma$-$2.95\sigma$ and yielding a SH0ES discrepancy of $2.65\sigma$-$3.45\sigma$. In this configuration, $H_{0}$ lies between the Planck and SH0ES determinations, partially alleviating their discrepancy. This behavior can be interpreted as a consequence of the modified matter evolution in Rastall cosmology, which induces a degeneracy between $\epsilon$ and $H_{0}$ at the background level. Small positive values $\epsilon\sim\mathcal{O}(10^{-3})$, favored by CMB-inclusive datasets, lead to a mild but systematic increase in $H_{0}$ relative to $\Lambda$CDM, while the effect vanishes as $\epsilon\to 0$.

Model comparison using AIC and BIC shows statistical equivalence between $\Lambda$CDM and R-$\Lambda$CDM for D1 and D2, while in D3 the R-$\Lambda$CDM extension is preferred in flat and closed geometries and decisively favored in the open case. Overall, the results highlight that any apparent alleviation of the Hubble tension remains strongly dependent on dataset composition and does not provide a universal resolution within the considered framework.

keyword: Hubble tension, Rastall gravity

I Introduction

Modern cosmology is currently confronted with several significant challenges, among which the Hubble constant $H_{0}$ tension and the nature of dark energy stand out as particularly prominent. Understanding these issues may point toward new physics and has therefore motivated extensive theoretical and observational investigations. Various experiments have been conducted to explore the properties of dark energy and the accelerated expansion of the Universe. Cosmologists have been striving to find an appropriate model to explain the late-time positive cosmic acceleration, leading to extensive research on gravitational theories riess1998observational ; perlmutter1999measurements ; spergel2003first ; peebles2003cosmological ; Pavlov:2013nra ; Adil_2023 ; Escamilla_2024 ; Vincenzi_2024 .

Einstein’s General Relativity (GR), supplemented by a cosmological constant, remains the preferred theory of gravity for explaining the accelerated expansion of the Universe. However, GR encounters conceptual and phenomenological challenges on cosmological scales, including singularities and parameter degeneracies. These limitations have stimulated extensive efforts toward extending or modifying the gravitational framework. Inflationary cosmology is one of the most successful theoretical developments in modern cosmology. Although originally proposed to address the classical problems of the Big Bang scenario, its predictions, such as the generation of primordial quantum fluctuations, and its remarkable agreement with observations have established it as a cornerstone of contemporary cosmological theory Starobinsky:1980te ; Guth:1980zm ; Albrecht:1982wi ; Linde:1981mu ; Linde:1983gd ; Spalinski:2007dv ; Bessada:2009pe ; Nazavari:2016yaa ; Amani:2018ueu ; Mohammadi:2018zkf ; berera1995warm ; berera2000warm ; BasteroGil:2004tg ; Rosa:2018iff ; Bastero-Gil:2019gao ; Sayar:2017pam ; Akhtari:2017mxc ; Sheikhahmadi:2019gzs ; Bastero-Gil:2017wwl ; Motohashi:2014ppa ; Odintsov:2017hbk ; Oikonomou:2017bjx ; Mohammadi:2019qeu ; Mohammadi:2020ftb ; Mohammadi:2022vru ; Yogesh:2025wak ; Mohammadi:2025gbu . Although inflation pertains to the earliest moments of cosmic history, its consequences extend to late-time observables, including implications for the current Hubble rate.

To explain the observed cosmic acceleration, the cosmological constant and dark energy remain popular candidates. Dark energy models constitute one of the most well-known approaches to describe the late-time acceleration of the Universe. On the other hand, cosmologists have explored various modified gravity theories, such as $f(R)$ gravity DeFelice:2010aj ; Sotiriou:2008rp ; nojiri2011unified , $f(T)$ gravity Bengochea:2008gz ; Wu:2010mn ; Myrzakulov:2010vz , $f(R,T)$ gravity harko2011f ; Carvalho:2017pgk ; Mohammadi:2023kzd ; Ossoulian:2023moq ; Taghavi:2023ptn , modified Gauss-Bonnet gravity $f(G)$ Li:2007jm ; Davis:2007ta ; Yogesh:2025hll , scalar-tensor gravity Brans:1961sx ; Carloni:2007eu ; Mohammadi:2022fiv , brane gravity Maartens:2010ar ; Borzou:2009gn ; Mohammadi:2020ake , among others. Modified gravitational theories generally extend the Einstein–Hilbert action of General Relativity through modifications of the curvature-based description of gravity. In this context, Rastall gravity has emerged as a phenomenologically interesting alternative, based on a generalized conservation law for the energy–momentum tensor AFSHAR2023101357 ; SALTI2020100630 ; SALEEM2021100808 ; Nazavari:2023khn .

Rastall gravity was originally proposed as a modification of the standard conservation law of the energy–momentum tensor Rastall1972 ; Rastall1976 . In contrast to General Relativity, where the covariant divergence of the energy–momentum tensor vanishes, Rastall theory allows a non-zero divergence proportional to the gradient of the Ricci scalar, leading to an effective non-minimal coupling between matter and geometry.

The theoretical status of Rastall gravity has been the subject of ongoing discussion. In particular, it has been argued that the theory can be formally recast as standard General Relativity with a redefined (non-conserved) effective energy–momentum tensor Visser2017 . On the other hand, several studies have pointed out that, depending on the physical interpretation and the matter sector, the theory may still lead to phenomenological deviations from GR in specific astrophysical and cosmological contexts Darabi2018 . In this work, we adopt the commonly used phenomenological viewpoint and treat Rastall gravity as an effective framework to test possible deviations from the standard conservation law at cosmological scales.

From an observational perspective, constraints on the Rastall parameter $\epsilon$ have been obtained across a wide range of physical systems. Analyses of compact stellar structures typically constrain $\epsilon$ at the level $\mathcal{O}(10^{-2})$ Oliveira2015 ; ElHanafy2022 ; ElHanafy2023 , while galaxy–galaxy strong gravitational lensing studies have reported larger values of the order $\mathcal{O}(10^{-1})$ Li2019 . At the cosmological level, recent analyses combining background expansion and linear perturbations in the flat case indicate that $\epsilon$ is strongly constrained to be very close to zero, with deviations at most of order $\mathcal{O}(10^{-4})$ Akarsu2020 .

Taken together, these results indicate that the allowed magnitude of the Rastall parameter is strongly dependent on the observational scale and the underlying datasets. This scale dependence motivates further tests of Rastall cosmology using different observational probes and dataset combinations.

Among the outstanding issues in modern cosmology, the tension in measurements of the Hubble constant $H_{0}$ remains one of the most persistent challenges. In recent years, this discrepancy has become increasingly prominent due to the divergence between two independent measurement approaches. The Planck collaboration, based on observations of the cosmic microwave background (CMB), has reported a value of $H_{0}=67.4\pm 0.5$ km s*-1* Mpc*-1* aghanim2020planck ; Planck:2013pxb ; Planck:2015fie . In contrast, the SH0ES (Supernovae $H_{0}$ for the Equation of State) team, utilizing distance ladder measurements with the Hubble Space Telescope (HST), has obtained a significantly higher value of $H_{0}=73.2\pm 1.3$ km s*-1* Mpc*-1* riess2021large . This discrepancy, often referred to as the ”Hubble tension”, has profound implications for our understanding of the expansion history of the Universe and may indicate new physics beyond the standard cosmological model Riess:2019cxk ; riess2021large ; CAPOZZIELLO2023101201 .

The HST method involves measurements of the local distance ladder by combining photometry of Cepheid variables (via their period–luminosity relation) with other local distance anchors, such as Milky Way parallaxes and calibration distances to Cepheids in nearby galaxies hosting Type Ia supernovae. Specifically, for estimating the Hubble constant, observations of Cepheids in the Large Magellanic Cloud are employed by the HST group riess2021large . Currently, the discrepancy between $H_{0}$ estimations by the Hubble Space Telescope and those inferred from Planck exceeds $4\sigma$. This tension has persisted over the years, as summarized in Table 1. The problem can be interpreted as a discrepancy between early-Universe and late-Universe cosmological observations. The HST group relies on late-time data, while the Planck collaboration combines observations over a wide redshift range ($0<z<1100$), adopting the $\Lambda$CDM model as the fiducial cosmological framework.

Different methods have been proposed to address the $H_{0}$ tension. Several studies have examined $H_{0}$ estimations using approaches independent of both the Cepheid distance scale and CMB anisotropies verde2019tensions ; di2021snowmass2021 . Among these, the following methods have provided noteworthy results:

• •

Strong gravitational lensing time-delay measurements (e.g., the H0LiCOW project) wong2020h0licow ; denzel2021hubble .

• •

CMB lensing analyses wong2020h0licow .

• •

The tip of the red giant branch (TRGB) method used by the Carnegie–Chicago Hubble Program (CCHP) freedman2019carnegie ; freedman2020calibration .

• •

Megamaser-hosting galaxies baxter2021determining .

• •

Oxygen-rich variable stars (Miras) Pesce_2020 .

Some studies have suggested that either Planck or HST measurements may be affected by systematic uncertainties; however, no consensus explanation has yet resolved the discrepancy.

As shown in Table 1, the reported values of the Hubble constant vary between the Planck 2018 result and the SH0ES measurement. Recent late-Universe measurements tend to favor values closer to the SH0ES determination, while early-Universe probes based on the CMB remain consistent with the lower Planck value.

In light of this persistent $H_{0}$ tension, various extensions to the standard cosmological model have been proposed. These include modified gravity theories ballesteros2020h0 ; zumalacarregui2020gravity ; Odintsov:2020qzd , alternative dark energy scenarios huang2016dark , frameworks involving additional relativistic species d2018hot , and models with interactions between dark energy and dark matter di2017can .

In this paper, we conduct a detailed investigation of the Hubble tension within the context of Rastall gravity. This framework introduces a generalized conservation law for the energy-momentum tensor, allowing for an effective interaction between matter fields and spacetime geometry. By applying the Rastall hypothesis in a cosmological setting, we analyze how such modifications affect the expansion history of the Universe and whether the resulting dynamics can alleviate the $H_{0}$ tension. This approach enables us to confront modified gravity parameters directly with current cosmological observations.

This paper is organized as follows. In Section II, we describe the cosmological models and parameters considered in this study. The observational datasets are presented in Section III. Section IV outlines the statistical methodology used for model comparison and parameter estimation. The results, including constraints on cosmological parameters and implications for the $H_{0}$ tension, are discussed in Section V. Finally, Section VI summarizes the main findings and outlines possible directions for future research. Appendix A provides the statistical estimator used to quantify the tension between datasets.

II Cosmological Models

To determine the Hubble parameter $H(z)$ in different cosmological models, we begin with Einstein’s field equations in general relativity

[TABLE]

where $g_{\mu\nu}$ is the metric tensor, $G$ is Newton’s gravitational constant, and throughout this section we adopt natural units with $8\pi G=1$. The tensors $R_{\mu\nu}$ and $R$ denote the Ricci tensor and Ricci scalar, respectively. The energy-momentum tensor $T_{\mu\nu}$ describes the total energy content of the universe, including matter, radiation, and dark energy, and is defined as

[TABLE]

where the index $i$ labels the individual components: baryonic matter and cold dark matter (collectively denoted as $i=m$), radiation ($i=r$), and dark energy ($i=\Lambda$). Each component is modeled as a perfect fluid, characterized by the energy-momentum tensor

[TABLE]

where $\rho_{i}$ and $p_{i}$ are the energy density and pressure of the $i$-th fluid, and $u^{\mu}$ is the four-velocity. The conservation of energy-momentum, $T^{\mu}_{\nu;\mu}=0$, is compatible with the Bianchi identities, $G^{\mu}_{\nu;\mu}=0$ Petrov:2020wgy . Assuming spatial homogeneity and isotropy, the Einstein equations in a non-flat geometry reduce to the Friedmann equations

[TABLE]

[TABLE]

The Hubble parameter is defined as $H\equiv\dot{a}/a$, where $a(t)$ is the cosmological scale factor and the overdot denotes differentiation with respect to cosmic time. The curvature contribution can be formally rewritten as an effective energy density, $\rho_{k}\equiv-3k/a^{2}$, where $k=0$, $+1$, and $-1$ correspond to flat, closed, and open universes, respectively. The evolution of each fluid component is determined by its equation of state, given by $p_{i}=\omega_{i}\rho_{i}$, where $\omega_{i}$ is the equation of state (EoS) parameter. For pressureless matter and cold dark matter, $\omega_{m}=0$; for radiation, $\omega_{r}=1/3$; for a cosmological constant, $\omega_{\Lambda}=-1$; and for a generic dark energy component $X$, one has $\omega_{X}<-1/3$. In addition, the energy conservation equation is given by

[TABLE]

Assuming a constant EoS parameter $\omega_{i}$, the solution to Eq. (5) is

[TABLE]

where $\rho_{i0}$ and $a_{0}$ denote the present values of the energy density of the $i$-th component and the scale factor, respectively. In what follows, we adopt the normalization $a_{0}=1$. The dimensionless energy density parameters at the present time are defined as

[TABLE]

Using these definitions, the Friedmann equation (3) for the $\Lambda$CDM model can be written as

[TABLE]

where $z=(a_{0}/a)-1$ is the redshift, $H_{0}$ is the present Hubble parameter, and $\bm{p}$ denotes the set of model parameters. In the flat case, $\Omega_{k_{0}}=0$, while for open (closed) geometries $\Omega_{k_{0}}>0$ ($\Omega_{k_{0}}<0$). The density parameters satisfy the closure condition $\Omega_{\Lambda_{0}}+\Omega_{m_{0}}+\Omega_{r_{0}}+\Omega_{k_{0}}=1$.

II.1 Rastall Model

A notable approach to modifying General Relativity consists in relaxing the standard conservation law of the energy-momentum tensor. In this framework, the usual condition $T^{\mu}_{\nu;\mu}=0$ is modified, leading to deviations from standard geodesic motion and the appearance of an additional force Lobo:2022aop . Rastall originally proposed such a modification by allowing a non-vanishing covariant divergence of the matter energy–momentum tensor Rastall1972 ; Rastall1976 .

It is worth noting that the theoretical status of Rastall gravity has been the subject of ongoing discussion. In particular, it has been argued that the theory can be formally recast as standard General Relativity with a redefined non-conserved effective energy–momentum tensor Visser2017 . On the other hand, other studies have emphasized that, depending on the physical interpretation and the matter sector, Rastall gravity may still lead to phenomenological deviations from GR in specific astrophysical and cosmological contexts Darabi2018 .

In the present work, we adopt this phenomenological viewpoint and treat Rastall gravity as an effective framework for probing deviations from the standard conservation law at cosmological scales. Within this approach, the parameter $\epsilon$ can be directly constrained by cosmological observations, as discussed in the following sections.

The aim here is to place constraints on cosmological parameters in both the standard $\Lambda$CDM model and its generalized counterpart, referred to as the R–$\Lambda$CDM model, in order to investigate their effectiveness in alleviating the $H_{0}$ tension. In Rastall gravity, the modified Einstein field equations are written as

[TABLE]

where

[TABLE]

The effective gravitational coupling $\kappa$ satisfies

[TABLE]

Adopting $8\pi G=1$, Eq. (10) reduces to

[TABLE]

The requirement that $\kappa$ be real implies

[TABLE]

We discard the second branch as it corresponds to large deviations from general relativity. In what follows, we focus on the phenomenologically relevant regime

[TABLE]

which motivates the restriction

[TABLE]

Under the perturbative branch continuously connected to general relativity, expanding Eq. (11) around $\kappa=1$ yields

[TABLE]

The field equations then take the approximate form

[TABLE]

Taking the trace of the Rastall field equations yields the exact relation

[TABLE]

Using the perturbative solution given in Eq. (15), the Ricci scalar becomes

[TABLE]

where we consistently retain terms up to $\mathcal{O}(\epsilon^{2})$ throughout the cosmological analysis. Substituting back, we obtain

[TABLE]

For a perfect fluid,

[TABLE]

with equation-of-state parameter $\omega_{i}=p_{i}/\rho_{i}$. In a non-flat FLRW background, the modified Friedmann equation becomes

[TABLE]

The corresponding modified acceleration equation can be written as

[TABLE]

Expanding Eq. (21) up to $\mathcal{O}(\epsilon^{2})$ yields

[TABLE]

In the limit $\epsilon\to 0$, Eq. (22) reduces to the standard GR acceleration equation,

[TABLE]

The modified continuity equation reads

[TABLE]

whose solution is

[TABLE]

Equation (25) shows that the Rastall parameter $\epsilon$ modifies the effective dilution rate of non-relativistic matter relative to the standard $\Lambda$CDM scenario. In particular, for $\epsilon>0$ the matter density evolves faster than the standard $(1+z)^{3}$ scaling, leading to a mild suppression of late-time matter clustering, whereas $\epsilon<0$ produces a slower matter dilution and therefore tends to enhance clustering relative to GR. In contrast, radiation and vacuum energy preserve their standard evolution laws, namely $\rho_{r}\propto(1+z)^{4}$ and $\rho_{\Lambda}=\mathrm{const}$. Consequently, the main cosmological effect of Rastall gravity at the background level arises through the modified matter evolution. We adopt $a_{0}=1$ and consistently retain terms up to $\mathcal{O}(\epsilon^{2})$. For convenience, we define the rescaled density parameters

[TABLE]

which govern the background expansion. The Hubble expansion rate is therefore

[TABLE]

with

[TABLE]

and the closure condition

[TABLE]

We refer to these quantities as effective density parameters, $\Omega_{i}^{\rm eff}\equiv\tilde{\Omega}_{i}$, which reduce to the standard $\Lambda$CDM values in the limit $\epsilon\to 0$.

III Observational data analysis

In the following, we derive constraints on the parameters of both the standard and Rastall versions of the $\Lambda$CDM model. To this end, we apply the Markov Chain Monte Carlo (MCMC) method to explore the parameter space and maximize the likelihood function $L(\bm{p})$, or equivalently minimize $\chi^{2}(\bm{p})$ under the assumption of Gaussian likelihoods farooq2017hubble , with respect to the parameter vector $\bm{p}$, in order to determine the best-fit values $\bm{p}_{0}$. Confidence regions at the $1\sigma$, $2\sigma$, and $3\sigma$ levels are presented in terms of two-dimensional parameter contours.

For the calculation of $\chi^{2}(\bm{p})$, we use the following observational datasets:

• •

Big Bang Nucleosynthesis (BBN) constraints based on the updated PRIMAT framework Pitrou2018 ; PRIMAT ,

• •

$H(z)$ measurements, incorporating 40 data points as listed in Table 2,

• •

Type Ia Supernovae data from the Pantheon+ compilation, including the SH0ES calibration and the full STAT+SYS covariance matrix Brout_2022 ,

• •

CMB distance priors from Planck 2018, using the appropriate flat or curved priors and the associated covariance matrix aghanim2020planck ; Chen:2018dbv ,

• •

BAO measurements from the DESI DR2 Gaussian BAO dataset (GCcomb), employing the full mean vector and covariance matrix DESI2025_DR2_II ,

• •

Growth-rate measurements $f(z)\sigma_{8}(z)$, consisting of 26 data points as listed in Table 3.

Assuming statistical independence among the different datasets, the total likelihood for a given data combination is given by the product of the corresponding likelihood terms,

[TABLE]

and therefore the combined chi-square can be written as

[TABLE]

where each contribution is described in detail in the following subsections.

III.1 Big Bang Nucleosynthesis (BBN data)

Big Bang Nucleosynthesis (BBN) provides an independent constraint on the physical baryon density through the primordial light-element abundances. In this work, we employ the publicly available PRIMAT framework Pitrou2018 , using the most recent grid release (September 2024 update)555PRIMAT public grid: https://www2.iap.fr/users/pitrou/primat.htm. The PRIMAT calculation includes QED corrections and incomplete neutrino decoupling effects Froustey2020 , such that $\Delta N=0$ corresponds to an effective number of relativistic species $N_{\rm eff}\simeq 3.044$. We fix $N_{\rm eff}=3.046$ throughout the analysis, which is fully consistent with the PRIMAT prediction for $\Delta N=0$ and with the standard value adopted in Planck 2018 CMB analyses.

The theoretical prediction for the primordial helium-4 mass fraction, $Y_{p}^{\rm th}$, is obtained by interpolating the PRIMAT grid as a function of the physical baryon density $\omega_{b}\equiv\Omega_{b0}h^{2}$ and $\Delta N$, and evaluated at $\Delta N=0$. We compare this prediction with the observational determination $Y_{p}^{\rm obs}=0.2448\pm 0.0033$ Aver2022 , assuming a Gaussian uncertainty, and define the BBN contribution to the total chi-square as

[TABLE]

where $\sigma_{Y_{p}}=0.0033$. This BBN likelihood provides an observationally driven and independent constraint on the baryon density, complementary to late-time cosmological probes.

III.2 Hubble parameter measurements ($H(z)$ data)

We use 40 independent measurements of the Hubble parameter obtained from the cosmic chronometer (CC) method, as summarized in Table 2, spanning the redshift range $0.07\leq z\leq 2.3$. These measurements are based on the differential age evolution of passively evolving galaxies and provide direct determinations of the expansion rate. No direct $H_{0}$ prior at $z=0$ is included in this compilation.

We emphasize that none of the $H(z)$ measurements used in this work are derived from baryon acoustic oscillation (BAO) techniques. Therefore, the $H(z)$ dataset is statistically independent of the DESI BAO measurements included in our analysis, and no overlap or double counting is present.

Assuming no correlations among the measurements, the contribution of the $H(z)$ data to the total chi-square is defined as

[TABLE]

where $H_{\rm th}(z_{i};\bm{p})=H_{0}\,E(z_{i};\bm{p})$ denotes the theoretical prediction of the Hubble parameter at redshift $z_{i}$.

This dataset provides a direct probe of the late-time expansion history and plays a central role in constraining the Hubble constant within both the $\Lambda$CDM and R-$\Lambda$CDM frameworks.

III.3 Type Ia supernovae: Pantheon+SH0ES sample (SN data)

Type Ia supernovae (SNe Ia) provide precise measurements of the late-time expansion history through the luminosity distance–redshift relation. In this work, we employ the Pantheon+SH0ES compilation Brout_2022 , which consists of 1701 light curves of 1555 spectroscopically confirmed SNe Ia in the redshift range $0.001<z<2.3$. This sample includes Cepheid-calibrated distances from the SH0ES program and provides observed distance moduli together with the full statistical and systematic covariance matrix. The theoretical distance modulus is defined as

[TABLE]

where the luminosity distance is

[TABLE]

with

[TABLE]

where $c$ denotes the speed of light. The transverse comoving distance function is given by

[TABLE]

Distances are evaluated at the CMB-frame redshifts provided in the Pantheon + SH0ES catalogue. All observational uncertainties, including peculiar velocity corrections and calibration systematics, are incorporated through the full covariance matrix ${\sf C}$ (STAT+SYS). The absolute magnitude parameter $M_{B}$ is treated as a nuisance parameter and analytically marginalized over, such that the supernova data do not impose an explicit prior on $H_{0}$. The supernova chi-square is defined as

[TABLE]

Following the standard procedure Conley_2011 , we analytically marginalize over the nuisance parameter associated with the absolute magnitude $M_{B}$. Defining

[TABLE]

the marginalized chi-square can be written as

[TABLE]

where

[TABLE]

and $\mathbf{1}$ denotes a vector of ones with dimension equal to the number of supernovae. For numerical stability and efficiency, the inverse-covariance operations are implemented using a Cholesky decomposition of the covariance matrix.

III.4 Cosmic Microwave Background (CMB data)

We employ the CMB distance priors derived from Planck 2018 data, as provided in Chen:2018dbv . The compressed CMB information is encoded in the shift parameter $R$, the acoustic scale $l_{A}$, and the physical baryon density $\Omega_{b0}h^{2}$. The shift parameter is defined as

[TABLE]

where $D_{M}(z)$ denotes the transverse comoving distance evaluated consistently with the assumed spatial curvature, and $c$ is the speed of light. The acoustic scale is given by

[TABLE]

where $r_{s}(z_{*})$ is the comoving sound horizon at the redshift of recombination $z_{*}$. The sound horizon is computed as

[TABLE]

with the baryon sound speed

[TABLE]

Throughout this analysis, we assume standard radiation content with $N_{\rm eff}=3.046$. The recombination redshift $z_{*}$ is evaluated using the fitting formula of hu1996small . The CMB chi-square is defined as

[TABLE]

where $X_{\rm CMB}$ depends on the assumed spatial geometry. For spatially flat models, we adopt

[TABLE]

with inverse covariance matrix

[TABLE]

For non-flat (open or closed) models, we use

[TABLE]

with inverse covariance matrix

[TABLE]

Since the allowed values of the Rastall parameter correspond to small deviations from general relativity and affect only the background expansion at the perturbative level considered here, we consistently adopt the Planck 2018 distance priors derived within the standard cosmological framework.

III.5 Baryon Acoustic Oscillations (BAO data)

We use the Baryon Acoustic Oscillation (BAO) measurements from the Dark Energy Spectroscopic Instrument (DESI) Data Release 2 (DR2), based on three years of observations DESI2025_DR2_II ; DESI2025_DR2_I ; Lodha2025 . These data include post-reconstruction clustering measurements from over 14 million galaxies and quasars, providing 13 distance constraints over the redshift range $0.295<z<2.33$. The DESI DR2 sample includes: (i) isotropic measurements of $D_{V}(z)/r_{s}(z_{d})$ at low redshift, and (ii) anisotropic measurements of $D_{M}(z)/r_{s}(z_{d})$ and $D_{H}(z)/r_{s}(z_{d})$ at higher redshift, where $r_{s}(z_{d})$ is the comoving sound horizon at the baryon drag epoch. The transverse comoving distance is defined as

[TABLE]

where $E(z)=H(z)/H_{0}$ and $c$ is the speed of light. The curvature-dependent function $S_{k}(\chi)$ is given by

[TABLE]

The Hubble distance is

[TABLE]

The volume-averaged distance used for isotropic BAO measurements is

[TABLE]

which is equivalent to the standard definition expressed in terms of the angular diameter distance. The sound horizon at the drag epoch, $r_{s}(z_{d})$, is computed consistently with the background cosmology, assuming $N_{\rm eff}=3.046$, and using the same radiation content as in the CMB analysis. The BAO likelihood is assumed to be Gaussian in the measured distance ratios relative to $r_{s}(z_{d})$. The corresponding chi-square is defined as

[TABLE]

where $\mathbf{d}_{\rm obs}$ is the vector of the 13 DESI DR2 measurements ($D_{V}/r_{s}$, $D_{M}/r_{s}$, or $D_{H}/r_{s}$ depending on the redshift bin), $\mathbf{d}_{\rm th}$ is the corresponding theoretical prediction, and $\mathbf{C}$ is the full $13\times 13$ covariance matrix provided by the DESI collaboration for the combined galaxy+quasar sample (excluding Ly$\alpha$ forest).

Since the CMB distance priors and DESI BAO measurements probe different physical observables and are derived from statistically independent datasets, their likelihoods are treated as independent in the combined analysis.

III.6 The growth rate of cosmic structures ($f\sigma_{8}(z)$ data)

While the geometrical probes discussed in the previous sections primarily constrain the background expansion history, the growth of cosmic structures provides an independent and complementary test of cosmological models at the perturbation level. Measurements of the quantity $f(z)\sigma_{8}(z)$ directly probe the evolution of matter density fluctuations and are therefore sensitive to both the expansion rate and the underlying theory of gravity.

In this work, we use 26 independent measurements of $f(z)\sigma_{8}(z)$ compiled in Table 3. In the absence of a publicly available full covariance matrix for the combined growth sample, the data points are treated as statistically independent and a diagonal covariance matrix is assumed. The corresponding chi-square is therefore written as

[TABLE]

In the linear regime ($\delta\ll 1$), the evolution of matter perturbations in a homogeneous and isotropic background is governed by relativistic perturbation theory. For sub-horizon modes and in the absence of dark-energy clustering, the matter density contrast $\delta_{m}$ satisfies

[TABLE]

where a prime denotes differentiation with respect to the scale factor $a$, and

[TABLE]

The matter density parameter evolves as

[TABLE]

In the present analysis, we assume homogeneous dark energy with effective sound speed $C_{\rm eff}^{2}=1$ and vanishing dark-energy perturbations, $\delta_{\Lambda}=0$, which is well justified on sub-horizon scales where dark-energy fluctuations are strongly suppressed. For the Rastall model, we allow the modification to enter only through the background expansion rate $H(a)$, while the linear perturbation equation retains its standard form. In this approach, the effective density parameters $\tilde{\Omega}_{i}$ are used only in the background expansion history entering $E(a)$, whereas the perturbation source term is written in terms of the standard matter density parameter $\Omega_{m0}$. Therefore, Rastall corrections affect the growth evolution indirectly through the modified Hubble expansion rate. This approximation is consistent with the perturbative regime $|\epsilon|\ll 1$ considered in this work. The growth rate function and the matter fluctuation amplitude are defined as

[TABLE]

[TABLE]

where $\sigma_{8}(0)$ is treated as a free cosmological parameter in the MCMC analysis. The product $f(z)\sigma_{8}(z)$ constitutes a directly measurable quantity and provides a powerful probe for testing cosmological models at the perturbative level.

IV Constraints on Models

In this section, we employ the Markov Chain Monte Carlo (MCMC) technique to derive the best-fit values and confidence intervals of the cosmological parameters. Both the standard $\Lambda$CDM and the R-$\Lambda$CDM models are investigated under the assumptions of spatially flat and non-flat universes.

To assess the robustness of the constraints and to quantify the impact of different cosmological probes, we consider four complementary combinations of observational data:

• •

$\mathbf{D_{1}}$ = SN + $H(z)$ + $f\sigma_{8}(z)$,

• •

$\mathbf{D_{2}}$ = SN + $H(z)$ + $f\sigma_{8}(z)$ + BAO + BBN,

• •

$\mathbf{D_{3}}$ = SN + $H(z)$ + $f\sigma_{8}(z)$ + BAO + CMB,

• •

$\mathbf{D_{4}}$ = CMB (Planck 2018 distance priors only).

Here, SN denotes the Pantheon+SH0ES Type Ia supernova sample, $H(z)$ corresponds to the compilation of cosmic chronometer and BAO-derived Hubble measurements, $f\sigma_{8}(z)$ represents the growth-rate data, BAO refers to the DESI DR2 clustering measurements, BBN is implemented via the PRIMAT helium abundance constraint, and CMB corresponds to the Planck 2018 distance priors.

The dataset $\mathbf{D_{1}}$ contains exclusively late-time probes and therefore constrains the expansion history and growth of structure without imposing early-Universe information. The combination $\mathbf{D_{2}}$ supplements late-time data with primordial nucleosynthesis constraints, while $\mathbf{D_{3}}$ incorporates compressed CMB information from recombination.

In addition, the dataset $\mathbf{D_{4}}$ is introduced to isolate the contribution of early-Universe information by considering the Planck 2018 CMB distance priors alone. This allows for a direct determination of the Hubble constant from CMB observations, independently of late-time probes, and provides a reference baseline for interpreting the Hubble tension.

This hierarchical dataset construction enables us to explicitly evaluate how different observational regimes affect the inferred cosmological parameters, particularly the Hubble constant $H_{0}$ and the amplitude of matter fluctuations $\sigma_{8}$.

The results of the parameter estimation are summarized in Tables 4–14, while the corresponding two-dimensional confidence contours are illustrated in Figures 1–8.

Model selection criteria

To compare the statistical performance of the models, we employ the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC), which are widely used tools in cosmological model selection capozziello2011comprehensive ; kass1995bayes ; spiegelhalter2002bayesian ; whitehead2007selection ; liddle2007information ; rezaei2021comparison . The AIC and BIC are defined as

[TABLE]

[TABLE]

where $d$ is the number of free cosmological parameters in the model and $N$ is the total number of data points included in a given dataset combination (accounting for all measurements entering the likelihood). For model comparison within each dataset combination, we consider the differences

[TABLE]

where the reference value corresponds to the model with the minimum information criterion for that specific dataset. According to the $\Delta\mathrm{AIC}$ criterion, $\Delta\mathrm{AIC}\leq 2$ indicates substantial support for a model, $4\leq\Delta\mathrm{AIC}\leq 7$ suggests considerably less support, and $\Delta\mathrm{AIC}\geq 10$ implies that the model is strongly disfavored. Similarly, in the context of BIC, $\Delta\mathrm{BIC}\leq 2$ is regarded as weak evidence, $2\leq\Delta\mathrm{BIC}\leq 6$ indicates positive evidence, $6\leq\Delta\mathrm{BIC}\leq 10$ reflects strong evidence, and values exceeding 10 are interpreted as very strong evidence against the model with the higher BIC kass1995bayes .

V Numerical results

V.1 Constraining the flat and non-flat models of $\Lambda$CDM

In this subsection, we present the parameter constraints obtained for both the spatially flat and non-flat $\Lambda$CDM models using the dataset combinations defined in Section IV. For each geometry, we derive the marginalized best-fit values and $1\sigma$ confidence intervals of the cosmological parameters, together with the minimum chi-square and the corresponding information criteria.

For the spatially flat model, the sampled parameter set consists of ($\Omega_{dm0}$, $\Omega_{b0}$, $H_{0}$, $\sigma_{8}$, $M_{B}$), while the radiation density parameter is derived assuming fixed $N_{\rm eff}=3.046$. The dark-energy density parameter $\Omega_{\Lambda 0}$ is obtained from the Friedmann closure condition, $\Omega_{\Lambda 0}=1-\Omega_{m0}-\Omega_{r0}$, with $\Omega_{m0}=\Omega_{dm0}+\Omega_{b0}$. In the non-flat case, the curvature parameter $\Omega_{k0}$ is included as an additional free parameter, and the closure relation becomes $\Omega_{\Lambda 0}=1-\Omega_{m0}-\Omega_{r0}-\Omega_{k0}$.

We first examine the constraints obtained from the late-time dataset $\mathbf{D_{1}}$, which provides a baseline determination of $H_{0}$ and $\sigma_{8}$ without imposing early-Universe priors. We then investigate how the inclusion of primordial information through BBN ($\mathbf{D_{2}}$) or CMB distance priors ($\mathbf{D_{3}}$) modifies the posterior distributions. Particular attention is devoted to the behavior of the Hubble constant, the curvature parameter, and the stability of $\sigma_{8}$ across datasets.

In addition, we include the dataset $\mathbf{D_{4}}$, corresponding to the Planck 2018 CMB distance priors alone, as a reference baseline. This configuration allows us to isolate the constraints driven purely by early-Universe information and to directly compare them with those obtained from late-time and combined datasets. The D4 dataset is therefore used as an anchor to assess how the inclusion of different cosmological probes shifts the inferred value of $H_{0}$ relative to the Planck-preferred region.

To quantify the level of agreement or tension between different dataset combinations, we employ the statistical estimator described in Appendix A. This allows us to assess whether shifts in parameter means between $\mathbf{D_{1}}$, $\mathbf{D_{2}}$, $\mathbf{D_{3}}$, and the Planck baseline $\mathbf{D_{4}}$ are statistically significant or remain within the expected combined uncertainties.

The complete numerical results are summarized in Tables 4–7, while the corresponding two-dimensional marginalized confidence contours are shown in Figures 1–4. The posterior distributions for the Planck-only D4 dataset are presented separately in Fig. 2, illustrating the tight constraints provided by early-Universe observations. These results provide a comprehensive assessment of the internal consistency of the $\Lambda$CDM framework under both flat and non-flat assumptions.

V.1.1 Flat $\Lambda$CDM model

We perform a comprehensive analysis of the dataset combinations D1, D2, and D3 within the framework of the spatially flat $\Lambda$CDM cosmology. For reference, the Planck baseline (denoted here as D4) is included only for external comparison. The resulting 68% confidence level constraints on the cosmological parameters are summarized in Table 4. In this analysis, the sampled parameters are $\Omega_{dm_{0}}$, $\Omega_{b_{0}}$, $H_{0}$, and $\sigma_{8}$, while $\Omega_{\Lambda_{0}}$ and $\Omega_{r_{0}}$ are derived from the flatness condition and the fixed radiation sector ($N_{\rm eff}=3.046$). The corresponding one- and two-dimensional marginalized constraints for D1–D3 are shown in Fig. 1.

The inferred Hubble constant spans the range $67.35\leq H_{0}\leq 71.30$ km s*-1* Mpc*-1*. In particular, dataset D2 yields $H_{0}=71.04^{+0.72}_{-0.73}$ km s*-1* Mpc*-1*. This value differs from the Planck 2018 estimate ($67.4\pm 0.5$) by $4.13\sigma$ and from the SH0ES measurement ($73.2\pm 1.3$) by $1.45\sigma$. For comparison, the direct Planck–SH0ES discrepancy amounts to approximately $4.16\sigma$. Thus, late-time datasets shift the inferred expansion rate toward the local value, though without reconciling both measurements simultaneously.

Comparison of D1 with D2 and D3

Between D1 and D2, all parameters remain consistent within less than $1\sigma$. The deviations for $\Omega_{dm_{0}}$, $\Omega_{b_{0}}$, $H_{0}$, and $\sigma_{8}$ are $0.35\sigma$, $0.76\sigma$, $0.24\sigma$, and $0.08\sigma$, respectively. The derived parameters shift by $0.43\sigma$ ($\Omega_{\Lambda_{0}}$) and $0.24\sigma$ ($\Omega_{r_{0}}$).

In contrast, D3 (including CMB priors) shows larger deviations relative to D1: $1.18\sigma$ for $\Omega_{dm_{0}}$, $0.28\sigma$ for $\Omega_{b_{0}}$, $4.66\sigma$ for $H_{0}$, and $0.45\sigma$ for $\sigma_{8}$, with derived shifts of $2.44\sigma$ in $\Omega_{\Lambda_{0}}$ and $4.96\sigma$ in $\Omega_{r_{0}}$. These shifts primarily reflect the strong constraining power of CMB information on the background expansion, in particular on $H_{0}$, with corresponding shifts in derived quantities.

Comparison with the D4 (Planck) baseline

The dataset D4, corresponding to the Planck 2018 CMB distance priors alone, provides a direct determination of cosmological parameters from early-Universe information and serves as a reference baseline for interpreting the Hubble tension within a consistent framework.

Using the constraints summarized in Table 4, we find that the late-time dataset D1 exhibits a significant deviation from the Planck baseline, with a discrepancy of approximately $3.63\sigma$ in $H_{0}$ and $1.80\sigma$ in $\Omega_{\Lambda_{0}}$. A similar behavior is observed for D2, with differences of $3.58\sigma$ in $H_{0}$ and $1.93\sigma$ in $\Omega_{\Lambda_{0}}$.

In contrast, the CMB-inclusive dataset D3 shows a much closer agreement with the Planck baseline, with deviations reduced to $0.47\sigma$ in $H_{0}$ and $0.53\sigma$ in $\Omega_{\Lambda_{0}}$. This behavior is also clearly reflected in the posterior distributions shown in Figures 1 and 2, where the inclusion of CMB information leads to a significant tightening of the parameter constraints and shifts the inferred value of $H_{0}$ toward the Planck-preferred region.

These results demonstrate that the dataset D4 provides a stable early-Universe anchor for the cosmological parameters, while late-time probes alone (D1 and D2) systematically favor higher values of $H_{0}$. The inclusion of CMB information in D3 effectively bridges these regimes, highlighting that the observed variation in $H_{0}$ is primarily driven by the dataset composition rather than by intrinsic modifications of the background cosmological model.

AIC and BIC analysis for flat $\Lambda$CDM model

For D1, we obtain AIC = 1581.53 and BIC = 1608.92; for D2, AIC = 1593.34 and BIC = 1620.76; and for D3, AIC = 1642.93 and BIC = 1670.36. Since the number of free parameters is identical across datasets, these differences quantify the relative goodness-of-fit under different observational combinations rather than model preference.

Overall, the flat $\Lambda$CDM analysis reveals a pronounced dependence of the inferred $H_{0}$ value on dataset selection, with late-time combinations systematically preferring higher expansion rates while CMB-inclusive constraints restore consistency with Planck. No dataset combination simultaneously reconciles the Planck and SH0ES determinations. Hence, within the minimal spatially flat framework, the Hubble tension remains statistically significant.

Figure 2 shows the posterior distributions of the cosmological parameters for the flat $\Lambda$CDM model using the Planck 2018 CMB distance priors (D4 dataset). The contours are well localized and exhibit nearly Gaussian behavior, indicating that the compressed CMB likelihood provides tight and stable constraints on the standard cosmological parameters.

In particular, the Hubble constant is tightly constrained around $H_{0}\sim 67$–$68~{\rm km\,s^{-1}\,Mpc^{-1}}$, in excellent agreement with the Planck baseline result. The strong degeneracy directions among parameters such as $\Omega_{m0}$ and $\Omega_{\Lambda 0}$ follow the expected behavior of the standard $\Lambda$CDM framework, confirming the internal consistency of the dataset and the robustness of the inference.

V.1.2 Open $\Lambda$CDM model

The open $\Lambda$CDM scenario extends the standard framework by introducing spatial curvature as an additional free parameter. The corresponding 68% confidence level constraints obtained from D1, D2, and D3 are summarized in Table 5. The marginalized one- and two-dimensional likelihood contours for the three dataset combinations are shown in Fig. 3. The inferred Hubble constant spans the range $67.55\leq H_{0}\leq 70.92$ km s*-1* Mpc*-1*. Relative to the Planck 2018 determination ($67.4\pm 0.5$ km s*-1* Mpc*-1*), the deviations amount to $3.62\sigma$ (D1), $3.72\sigma$ (D2), and $0.26\sigma$ (D3). With respect to the SH0ES local measurement ($73.2\pm 1.3$ km s*-1* Mpc*-1*), the corresponding discrepancies are $1.48\sigma$, $1.63\sigma$, and $4.25\sigma$, respectively. Thus, as in the flat case, late-time dataset combinations favor systematically higher values of $H_{0}$, while the inclusion of CMB information restores agreement with the Planck baseline.

The curvature parameter $\Omega_{k_{0}}$ is positive for D1 and D2. In both cases the deviation from spatial flatness remains below the $2\sigma$ level. In contrast, D3 yields $\Omega_{k_{0}}=0.0002^{+0.0004}_{-0.0002}$, which is fully consistent with a spatially flat Universe. This demonstrates the strong constraining power of CMB distance information on spatial curvature.

AIC and BIC analysis for open $\Lambda$CDM model

The AIC and BIC values are (1582.21, 1615.07) for D1, (1595.00, 1627.91) for D2, and (1647.56, 1680.47) for D3, as reported in Table 5.

Since the number of free parameters remains identical across the three dataset combinations, these criteria do not constitute model selection in the usual sense. Instead, they quantify the relative statistical performance of the open $\Lambda$CDM model under different observational inputs.

Allowing for positive spatial curvature does not qualitatively alter the tension pattern observed in the flat case. The inferred curvature parameter remains statistically consistent with spatial flatness, and the shift in $H_{0}$ continues to depend primarily on the inclusion of CMB information. Consequently, introducing positive curvature does not provide a statistically meaningful resolution of the Hubble tension.

V.1.3 Closed $\Lambda$CDM model

Cosmological parameter constraints for the seven-parameter closed $\Lambda$CDM model, derived from the dataset combinations D1, D2, and D3, are summarized in Table 6. The corresponding marginalized likelihood contours for the three dataset combinations are shown in Fig. 4.

The Hubble constant spans the range $67.18\leq H_{0}\leq 71.43$ km s*-1* Mpc*-1*. Relative to the Planck 2018 determination ($67.4\pm 0.5$ km s*-1* Mpc*-1*), the deviations amount to $4.14\sigma$, $4.23\sigma$, and $0.38\sigma$ for D1, D2, and D3, respectively. Thus, only D3 remains statistically consistent with Planck.

With respect to the SH0ES local measurement ($73.2\pm 1.3$ km s*-1* Mpc*-1*), the discrepancies are $1.15\sigma$ (D1), $1.35\sigma$ (D2), and $4.51\sigma$ (D3). As in the flat and open cases, the late-time combinations D1 and D2 prefer higher values of $H_{0}$ closer to the local determination, whereas the inclusion of CMB information in D3 restores agreement with Planck while maintaining a substantial discrepancy with SH0ES. Therefore, the overall Planck–SH0ES tension is not resolved within the closed $\Lambda$CDM framework.

The curvature parameter $\Omega_{k_{0}}$ is negative in all datasets, corresponding to a closed spatial geometry. For D1 and D2, the inferred values ($-0.0306^{+0.0221}_{-0.0393}$ and $-0.0156^{+0.0113}_{-0.0199}$) remain statistically consistent with spatial flatness at approximately the $1\sigma$ level. In contrast, D3 yields $\Omega_{k_{0}}=-0.0047\pm 0.0013$, corresponding to a statistically significant $\sim 3.6\sigma$ deviation from zero. This indicates a statistically significant preference for negative curvature when CMB priors are included.

The dark matter density parameter $\Omega_{dm_{0}}$ varies within the interval $0.243\lesssim\Omega_{dm_{0}}\lesssim 0.263$ across the three dataset combinations, while the baryon density parameter $\Omega_{b_{0}}$ lies in the range $0.044$–$0.057$. The structure growth parameter $\sigma_{8}$ remains close to the $\sim 0.80$ level for the late-time combinations D1 and D2, whereas the inclusion of CMB priors in D3 shifts the preferred value to $\sigma_{8}\simeq 0.87$. This increase reflects the strong correlation between spatial curvature and the growth amplitude in the closed geometry, highlighting the impact of CMB distance information on both background and perturbation parameters.

AIC and BIC analysis for closed $\Lambda$CDM model

The AIC and BIC values for the closed $\Lambda$CDM model are reported in Table 6. For D1, we obtain AIC = 1583.72 and BIC = 1616.58; for D2, AIC = 1595.45 and BIC = 1628.36; and for D3, AIC = 1634.51 and BIC = 1667.42.

Since the number of free parameters remains identical across D1, D2, and D3, these criteria quantify the relative goodness-of-fit of the same cosmological model under different observational combinations rather than indicating intrinsic model preference.

In the closed $\Lambda$CDM scenario, CMB priors indicate a statistically significant ($\sim 3.6\sigma$) preference for negative curvature; however, this geometric extension does not enable a simultaneous agreement between early- and late-universe measurements of $H_{0}$. As in the flat and open geometries, late-time datasets favor higher expansion rates, whereas CMB-inclusive constraints align closely with Planck. Therefore, the persistence of the Hubble tension appears to be driven primarily by dataset dependence rather than by spatial curvature.

V.1.4 Assessment of geometrical models based on parameter consistency: The $\Lambda$CDM model

A comparative examination of the cosmological parameter differences between datasets D1, D2, and D3, summarized in Table 7, enables us to assess the internal consistency of the flat, open, and closed $\Lambda$CDM geometries. The deviations are expressed in units of the combined statistical uncertainty ($\sigma$) between D1 and the other dataset combinations, computed using symmetrized errors. This analysis is independent of external reference values (such as Planck or SH0ES) and instead provides a direct measure of the internal consistency among the dataset combinations themselves.

A clear pattern emerges across all three geometrical scenarios. For the D1–D2 comparison, the majority of parameters exhibit deviations at or below the $1\sigma$ level in flat, open, and closed models. This indicates strong mutual consistency among late-time datasets (SN, $H(z)$, $f\sigma_{8}$, BAO, and BBN), largely independent of the assumed spatial curvature.

In contrast, the comparison between D1 and D3 reveals substantially larger shifts in several parameters, most notably in $H_{0}$ and $\Omega_{r_{0}}$. In the flat model, the deviation in $H_{0}$ reaches $4.6\sigma$, while in the open and closed models it amounts to $3.8\sigma$ and $4.8\sigma$, respectively. A similar pattern is observed for $\Omega_{r_{0}}$, where deviations approach the $4$–$5\sigma$ level across geometries. It is important to note that $\Omega_{r_{0}}$ is not independently sampled in our analysis but is derived assuming a fixed effective number of relativistic species, $N_{\rm eff}=3.046$. Therefore, the apparent large shift in $\Omega_{r_{0}}$ primarily reflects its dependence on $H_{0}$ rather than any modification of the underlying radiation physics. These large shifts arise from the inclusion of CMB distance priors in D3, which strongly constrain the background geometry and expansion history.

The curvature parameter $\Omega_{k_{0}}$ itself shows only moderate deviations (below $2\sigma$) when comparing D1 with D3 in the non-flat cases. This suggests that spatial curvature is not the primary driver of parameter inconsistencies between dataset combinations. Rather, the dominant effect originates from the addition of CMB information, which pulls all geometrical scenarios toward values consistent with Planck.

The structure growth parameter $\sigma_{8}$ remains relatively stable across dataset combinations, particularly in the closed case where deviations are well below $1\sigma$. This further supports the conclusion that late-time probes are mutually consistent and that the principal shifts arise when early-universe information is incorporated.

In the closed $\Lambda$CDM scenario, CMB priors indicate a mild preference for negative curvature; however, this geometric extension does not enable a simultaneous agreement between early- and late-universe measurements of $H_{0}$. As in the flat and open geometries, late-time datasets favor higher expansion rates, whereas CMB-inclusive constraints align closely with Planck. Therefore, the persistence of the Hubble tension appears to be driven by dataset dependence rather than by spatial curvature.

V.2 Constraining the flat and non-flat models of R-$\Lambda$CDM

In this subsection, we investigate the R-$\Lambda$CDM model as a minimal extension of the standard cosmological scenario, characterized by the inclusion of the Rastall parameter $\epsilon$. We consider spatially flat, open, and closed geometries, and derive observational constraints using the four dataset combinations defined in Section IV.

For the spatially flat case, the sampled parameter set is given by ($\Omega_{dm0}$, $\Omega_{b0}$, $H_{0}$, $\sigma_{8}$, $\epsilon$, $M_{B}$), while the radiation density parameter is derived by fixing $N_{\rm eff}=3.046$. The effective density components entering the modified Friedmann equation are consistently constructed up to $\mathcal{O}(\epsilon^{2})$, and the dark-energy density parameter is obtained from the generalized closure relation.

In the non-flat scenario, the curvature parameter $\Omega_{k0}$ is included as an additional free parameter. This allows us to explore potential degeneracies between spatial curvature and the Rastall parameter, thereby providing a direct assessment of the interplay between geometric effects and modified-gravity contributions in determining the expansion history.

We begin by analyzing the constraints obtained from the late-time dataset $\mathbf{D_{1}}$, which isolates the impact of modified background dynamics in the absence of early-Universe priors. We then examine the effect of incorporating primordial information through Big Bang Nucleosynthesis ($\mathbf{D_{2}}$), as well as the combined constraints from baryon acoustic oscillations and CMB distance priors ($\mathbf{D_{3}}$). Finally, we consider the CMB-only dataset ($\mathbf{D_{4}}$), which serves as a direct Planck baseline. This dataset plays a crucial role in assessing the behavior of the Rastall parameter $\epsilon$ and its degeneracy with spatial curvature, independently of late-time observational inputs.

To quantify the level of agreement or tension among different dataset combinations, we employ the statistical estimator described in Appendix A. In particular, the inclusion of $\mathbf{D_{4}}$ enables a direct comparison with the Planck baseline, allowing us to assess the consistency in the determination of $H_{0}$ between early- and late-Universe datasets.

The numerical results are summarized in Tables 8–14, while the corresponding marginalized one- and two-dimensional confidence contours are presented in Figures 5–8. These results provide a systematic comparison between flat and curved geometries within the Rastall framework and explicitly highlight the role of the $\mathbf{D_{4}}$ dataset in establishing the Planck reference constraints, allowing us to assess the extent to which the inclusion of $\epsilon$ alleviates existing cosmological tensions.

V.2.1 R-$\Lambda$CDM model in a flat universe

In the flat R-$\Lambda$CDM framework, the deviation parameter $\epsilon$, commonly referred to as the Rastall parameter, introduces a modification to the standard $\Lambda$CDM dynamics. This modification originates from the non-conservation of the matter energy-momentum tensor due to a non-minimal coupling between matter and geometry (for further details, see subsection II.1). The results derived from this model, utilizing three distinct datasets D1, D2, and D3, are summarized in Table 8 and depicted in Figure 5.

The Hubble constant, $H_{0}$, is estimated to be approximately $71$ km s*-1* Mpc*-1* for datasets D1 and D2, while D3 yields a lower value close to $68.5$ km s*-1* Mpc*-1*. The discrepancies between these estimates and those obtained from Planck and the local SH0ES measurements show a clear dataset dependence. Relative to the Planck value, D1 and D2 remain in tension at roughly the $4\sigma$ level, whereas D3 reduces this discrepancy to about $1.7\sigma$. Conversely, with respect to the local SH0ES determination, D1 and D2 are consistent within approximately $1.3\sigma$, while D3 exhibits a larger deviation of about $3.5\sigma$. These results indicate that the flat R-$\Lambda$CDM scenario shifts the preferred value of $H_{0}$ depending on the adopted data combination, improving agreement with Planck when CMB information is included, but remaining closer to the local determination when only late-time data are considered.

The effective dark matter density parameter, $\Omega_{dm0}^{\rm eff}$, takes nearly identical central values in datasets D1 and D3, while D3 substantially tightens the associated uncertainty. In contrast, D2 prefers a lower value, indicating a noticeable shift in the inferred effective matter sector when the data combination is modified. Similarly, the effective baryonic matter density parameter, $\Omega_{b0}^{\rm eff}$, shows a strong dependence on the inclusion of CMB information. While D1 and D2 allow comparatively broad ranges, D3 yields a very tightly constrained value. The effective dark energy density parameter, $\Omega_{\Lambda 0}^{\rm eff}$, remains within a narrow interval across all datasets, exhibiting moderate shifts between D2 and the other combinations. The effective radiation density parameter, $\Omega_{r0}^{\rm eff}$, remains small in all cases, at the expected order of $10^{-5}$, although its precision improves once additional high-redshift information is incorporated.

The Rastall deviation parameter, $\epsilon$, displays a clear dataset dependence. Dataset D1 favors a negative value, whereas D2 and D3 yield values fully consistent with zero, with D3 providing significantly tighter constraints due to the inclusion of CMB information. This trend suggests that the departure from the standard $\Lambda$CDM limit becomes progressively suppressed as tighter observational constraints, particularly from CMB data, are included.

The parameter $\sigma_{8}$, which characterizes the amplitude of matter fluctuations, is constrained within the interval $0.80$–$0.82$ across the three datasets, with D3 yielding slightly lower central values and tighter uncertainties. These values remain broadly compatible with CMB-based estimates within the quoted uncertainties.

Finally, information criteria, namely AIC and BIC, indicate that the dataset combination D1 yields the lowest values among the three cases considered, while D2 and D3 result in comparatively larger values for their respective data combinations.

Comparison of dataset D1 with D2 and D3

The comparison of constrained parameters derived from the D1 and D2 datasets in the flat R-$\Lambda$CDM model reveals noticeable differences across several effective parameters. For example, the effective dark matter density parameter, $\Omega_{dm0}^{\rm eff}$, shows a shift of approximately $1.53\sigma$, indicating a lower value in D2. Similarly, the effective radiation density $\Omega_{r0}^{\rm eff}$ and the effective cosmological constant $\Omega_{\Lambda 0}^{\rm eff}$ exhibit deviations of approximately $2.45\sigma$ and $1.06\sigma$, respectively. The parameter $\epsilon$ shows a deviation of about $2.26\sigma$, reflecting a significant dataset dependence. In contrast, the Hubble constant $H_{0}$ remains highly consistent across the two datasets, differing only at the $\sim 0.23\sigma$ level. These discrepancies, illustrated in Table 14, reflect the sensitivity of the effective parameter constraints to the adopted data combination.

The cosmological parameters derived from the D1 and D3 datasets show partial consistency but also reveal notable shifts in specific sectors. For instance, $H_{0}$ differs by approximately $3.07\sigma$. The effective radiation density $\Omega_{r0}^{\rm eff}$ and the Rastall parameter $\epsilon$ exhibit deviations of about $4\sigma$ and $2.32\sigma$, respectively. Meanwhile, $\Omega_{dm0}^{\rm eff}$ remains consistent at the central-value level, though D3 significantly tightens its uncertainty. As summarized in Table 14, these results demonstrate that certain effective parameters are strongly affected by the inclusion of additional datasets.

Comparison with the D4 (Planck) baseline

The dataset D4, based on the Planck 2018 CMB distance priors alone, provides a direct early-Universe reference for the flat R-$\Lambda$CDM model. As shown in Table 8, this dataset yields $H_{0}=67.13^{+2.28}_{-1.04}~{\rm km\,s^{-1}\,Mpc^{-1}}$ and $\epsilon=0.0023^{+0.0048}_{-0.0065}$, indicating that the Rastall parameter is consistent with the general relativity limit $\epsilon\simeq 0$, albeit with relatively large uncertainties due to parameter degeneracies in the compressed CMB likelihood.

In comparison, the late-time datasets D1 and D2 prefer higher values of the Hubble constant, $H_{0}\sim 71~{\rm km\,s^{-1}\,Mpc^{-1}}$, with D1 additionally showing a mild preference for negative values of $\epsilon$, while D2 remains fully consistent with $\epsilon=0$. These results reflect the increased flexibility of the model when only late-time information is considered.

A qualitatively different behavior emerges for the CMB-inclusive dataset D3, which yields $H_{0}=68.50^{+0.40}_{-0.41}~{\rm km\,s^{-1}\,Mpc^{-1}}$ and tightly constrains the Rastall parameter to $\epsilon=0.00203^{+0.00048}_{-0.00049}$. In this case, the inclusion of CMB information not only reduces the uncertainty in $H_{0}$, but also selects a small positive value of $\epsilon$ at high significance.

This comparison highlights the central role of dataset composition in the R-$\Lambda$CDM framework. While the Planck-only dataset D4 anchors the model close to the $\Lambda$CDM limit with $\epsilon\approx 0$, it is the combined dataset D3 that breaks parameter degeneracies and reveals a preference for small positive values of $\epsilon$. This, in turn, induces a mild but systematic upward shift in $H_{0}$ relative to the standard $\Lambda$CDM case, providing a direct observational manifestation of the $\epsilon$–$H_{0}$ degeneracy at the level of the background expansion.

Overall, these results demonstrate that the effect of the Rastall parameter on the inferred value of the Hubble constant is not intrinsic to the model alone, but emerges from the interplay between early- and late-Universe datasets. In particular, the enhancement of $H_{0}$ becomes visible only when $\epsilon$ is constrained to small positive values by CMB-inclusive data, while it remains suppressed when $\epsilon$ is weakly constrained or consistent with zero.

Figure 6 presents the corresponding posterior distributions for the flat R-$\Lambda$CDM model. Compared to the standard case, the inclusion of the Rastall parameter $\epsilon$ leads to a noticeable broadening of the parameter space and introduces additional degeneracies, most prominently between $\epsilon$ and $H_{0}$.

The contour structure clearly indicates that small positive values of $\epsilon$ are associated with slightly higher values of the Hubble constant, providing a direct visual manifestation of the $\epsilon$–$H_{0}$ degeneracy discussed in this work. At the same time, the Planck-only dataset strongly constrains $\epsilon$ to remain close to zero, limiting the magnitude of the deviation from the standard $\Lambda$CDM scenario.

This behavior highlights that, although the Rastall framework allows for additional flexibility in the background dynamics, the high precision of CMB observations restricts the allowed deviations, resulting in only a mild shift in the inferred value of $H_{0}$.

Model Selection via AIC and BIC

To assess the relative statistical performance of the flat R-$\Lambda$CDM model, we compare the AIC and BIC values obtained for dataset combinations D1, D2, and D3 (Table 9). The differences are defined as $\Delta$AIC $=\mathrm{AIC}_{\rm R}-\mathrm{AIC}_{\Lambda}$ and similarly for $\Delta$BIC.

For D1, the R-$\Lambda$CDM model yields $\Delta$AIC$=-1.52$, indicating a marginal preference over the standard $\Lambda$CDM model according to the usual $|\Delta\mathrm{AIC}|<2$ criterion. However, the corresponding $\Delta$BIC$=+3.96$ suggests that the additional parameter $\epsilon$ is moderately penalized once model complexity is taken into account.

For D2, both information criteria favor the standard $\Lambda$CDM model. The AIC difference is $\Delta$AIC$=+1.78$, indicating only weak statistical evidence against the R-$\Lambda$CDM model, while the BIC difference $\Delta$BIC$=+7.27$ indicates strong evidence against the more complex R-$\Lambda$CDM scenario within this dataset configuration.

A markedly different behavior is observed for D3. In this case, the R-$\Lambda$CDM model is strongly favored over the standard scenario, with $\Delta$AIC$=-14.88$ and $\Delta$BIC$=-9.40$. According to conventional model selection thresholds, these values correspond to decisive statistical support for the R-$\Lambda$CDM model when BAO and CMB data are included simultaneously.

This preference can be interpreted as a direct consequence of the additional flexibility introduced by the Rastall parameter $\epsilon$, which allows the model to better accommodate the combined constraints from late- and early-Universe observables, particularly through its impact on the background expansion history.

It is important to emphasize that the absolute AIC and BIC values are meaningful only within a fixed dataset configuration. Therefore, comparisons should be interpreted separately for each dataset combination rather than across D1, D2, and D3.

We also note that the D4 (Planck-only) dataset is not included in this comparison, since information criteria are not directly comparable across datasets with substantially different effective numbers of data points, and such a comparison would therefore not be statistically meaningful.

From the perspective of the Hubble tension, these results further reinforce the dataset-dependent nature of the inferred cosmological parameters. In particular, the strong statistical preference for the R-$\Lambda$CDM model in D3 correlates with the regime where small positive values of $\epsilon$ are favored and the inferred value of $H_{0}$ is shifted toward intermediate values between the Planck and SH0ES determinations. This behavior highlights that the apparent improvement in model performance is closely tied to the interplay between early- and late-Universe datasets, rather than representing a universal resolution of the tension.

Overall, the information criteria analysis indicates that the statistical viability of the R-$\Lambda$CDM model strongly depends on dataset composition. While it is moderately disfavored in D2, it receives substantial support in D3, highlighting the key role of early-Universe constraints in assessing deviations from the standard $\Lambda$CDM dynamics.

V.2.2 R-$\Lambda$CDM model in an open universe

The cosmological parameter estimates for the open R-$\Lambda$CDM model, obtained using datasets D1, D2, and D3, are summarized in Table 10 and illustrated in Figure 7. In this framework, effective density parameters $\Omega_{i}^{\rm eff}$ govern the background expansion, while the curvature parameter $\Omega_{k0}$ remains unrescaled.

The Hubble constant $H_{0}$ is estimated to be $70.96\pm 0.84$, $70.79\pm 0.76$, and $69.50\pm 0.51$ km s*-1* Mpc*-1* for datasets D1, D2, and D3, respectively. These values lie between the Planck 2018 reference ($67.4\pm 0.5$ km s*-1* Mpc*-1*) and the local SH0ES measurement ($73.2\pm 1.3$ km s*-1* Mpc*-1*).

In the standard flat $\Lambda$CDM model, the tension between Planck and SH0ES is approximately $4.17\sigma$. In the open R-$\Lambda$CDM model, the tension with Planck is reduced to about $3.63\sigma$, $3.73\sigma$, and $2.92\sigma$ for D1, D2, and D3, respectively. Similarly, the tension with SH0ES is approximately $1.45\sigma$, $1.61\sigma$, and $2.64\sigma$. Thus, the open R-$\Lambda$CDM framework shifts $H_{0}$ toward intermediate values, with D3 moving toward the early-Universe calibration when CMB information is included.

Regarding other parameters, the effective dark matter density $\Omega_{dm0}^{\rm eff}$ is around $0.24$ for D1 and D2, while D3 favors a higher value near $0.274$. The effective baryonic density $\Omega_{b0}^{\rm eff}$ remains broadly consistent across datasets, with the tightest constraint obtained for D3.

The spatial curvature parameter $\Omega_{k0}$ is positive in all cases but exhibits substantial dataset dependence. D1 allows the largest curvature, $\Omega_{k0}\sim 0.09$, D2 prefers a smaller value around $0.04$, and D3 is consistent with a nearly flat Universe, with $\Omega_{k0}$ close to zero and tightly constrained.

The effective dark energy density $\Omega_{\Lambda 0}^{\rm eff}$ lies in the range $0.63$–$0.68$, with D1 yielding the lowest value and D2 and D3 showing mutually consistent results. The effective radiation density $\Omega_{r0}^{\rm eff}$ remains at the expected order of $10^{-5}$ in all cases, with tighter constraints obtained when BAO and CMB information are included.

The Rastall deviation parameter $\epsilon$ exhibits a distinct pattern. Dataset D1 prefers a negative value, $\epsilon\approx-0.061$, corresponding to roughly a $2\sigma$ deviation from zero, whereas D2 and D3 yield values consistent with $\epsilon=0$ within $1\sigma$. This indicates that the inclusion of BAO and especially CMB data drives the model closer to the standard $\Lambda$CDM limit.

The structure growth parameter $\sigma_{8}$ decreases from approximately $0.87$ in D1 to $0.80$ in D3, with D2 taking an intermediate value.

Finally, the information criteria show that D1 yields the lowest AIC and BIC values within the open R-$\Lambda$CDM scenario, while D3 produces the largest values. Since the dataset compositions differ, these criteria should be interpreted separately within each configuration.

Comparison of dataset D1 with D2 and D3

Comparing D1 and D2 reveals moderate deviations in $\Omega_{k0}$ ($0.7\sigma$), $\epsilon$ ($2.0\sigma$), and $\Omega_{r0}^{\rm eff}$ ($2.0\sigma$), while $\Omega_{dm0}^{\rm eff}$, $\Omega_{\Lambda 0}^{\rm eff}$, $H_{0}$, and $\sigma_{8}$ remain statistically consistent within $1\sigma$. Overall, most parameters are compatible within uncertainties.

In contrast, comparing D1 and D3 shows stronger differences in curvature and dark matter density. The curvature parameter $\Omega_{k0}$ differs at approximately the $1.4\sigma$ level, $\Omega_{dm0}^{\rm eff}$ shifts by about $1.1\sigma$, and $\sigma_{8}$ differs by roughly $1.6\sigma$. The Hubble constant differs by approximately $1.5\sigma$. Despite these variations, the parameter sets remain statistically compatible within the open R-$\Lambda$CDM framework.

Model Selection via AIC and BIC

The relative statistical performance of the open R-$\Lambda$CDM model was evaluated using the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) for dataset combinations D1, D2, and D3 (Table 11). The differences are defined as $\Delta$AIC$=\mathrm{AIC}_{\rm R}-\mathrm{AIC}_{\Lambda}$ and similarly for $\Delta$BIC.

For D1, the R-$\Lambda$CDM model yields $\Delta$AIC$=-0.59$, indicating statistical equivalence with the standard $\Lambda$CDM model according to the conventional $|\Delta\mathrm{AIC}|<2$ criterion. However, the corresponding $\Delta$BIC$=+4.89$ shows a moderate penalty for the additional parameter when model complexity is more strongly weighted.

For D2, the AIC difference ($\Delta$AIC$=+1.38$) indicates statistical equivalence between the two models, while the BIC difference ($\Delta$BIC$=+6.87$) provides strong evidence in favor of the standard $\Lambda$CDM model due to the stronger penalty on model complexity.

A qualitatively different behavior is observed for D3. In this case, the open R-$\Lambda$CDM model significantly outperforms the standard $\Lambda$CDM model, with $\Delta$AIC$=-25.60$ and $\Delta$BIC$=-20.11$. According to conventional model selection thresholds, such large negative differences correspond to strong statistical evidence in favor of the R-$\Lambda$CDM scenario when BAO and CMB data are simultaneously included.

As emphasized above, the absolute values of AIC and BIC are meaningful only within a fixed dataset configuration.

Overall, the information criteria analysis indicates that the statistical viability of the open R-$\Lambda$CDM model depends sensitively on dataset composition. While it remains statistically equivalent to $\Lambda$CDM in D1 and statistically equivalent under AIC but disfavored by BIC in D2, it receives decisive support in D3, highlighting the crucial role of CMB constraints in assessing deviations from the standard cosmological model.

V.2.3 R-$\Lambda$CDM model in a closed universe

Table 12 summarizes the cosmological parameter estimates of the closed R-$\Lambda$CDM model obtained from the three independent observational datasets D1, D2, and D3. The corresponding likelihood contours are shown in Figure 8. In this framework, the effective density parameters $\Omega_{i}^{\rm eff}$ govern the background expansion, while the curvature parameter $\Omega_{k0}$ remains unrescaled.

The Hubble constant is constrained to be $H_{0}=71.48\pm 0.84$, $71.15\pm 0.75$, and $68.86\pm 0.54$ km s*-1* Mpc*-1* for D1, D2, and D3, respectively (symmetrized $1\sigma$ uncertainties). Relative to the Planck 2018 determination ($67.4\pm 0.5$ km s*-1* Mpc*-1*), the corresponding tensions are approximately $4.19\sigma$, $4.16\sigma$, and $1.99\sigma$. Compared to the local SH0ES measurement ($73.2\pm 1.3$ km s*-1* Mpc*-1*), the tensions are $1.11\sigma$, $1.37\sigma$, and $3.09\sigma$, respectively. Thus, in the absence of CMB information (D1 and D2), the closed R-$\Lambda$CDM model favors values of $H_{0}$ closer to the local distance-ladder measurement, whereas the inclusion of CMB data in D3 shifts the result toward the Planck value and reduces the Planck tension to the $\sim 2\sigma$ level.

The effective dark matter density parameter $\Omega_{dm0}^{\rm eff}$ ranges from $0.243$ (D2) to $0.293$ (D1), with D3 yielding an intermediate value of $0.270$. The effective baryon density $\Omega_{b0}^{\rm eff}$ lies between $0.036$ and $0.057$, with the tightest constraint obtained for D3.

The spatial curvature parameter $\Omega_{k0}$ is negative in all cases, confirming a closed spatial geometry. Its magnitude decreases as additional early-Universe information is incorporated: $\Omega_{k0}\approx-0.039$ (D1), $-0.016$ (D2), and $-0.0013$ (D3). The D3 result indicates only a mild ($\sim 1\sigma$) preference for negative curvature and remains close to spatial flatness, highlighting the strong constraining power of CMB data on geometry.

The effective dark energy density $\Omega_{\Lambda 0}^{\rm eff}$ varies between $0.68$ and $0.72$, while the effective radiation density $\Omega_{r0}^{\rm eff}$ remains of order $10^{-5}$ in all cases, with progressively tighter uncertainties from D1 to D3.

The Rastall deviation parameter $\epsilon$ shows a clear dataset dependence. D1 favors a negative value, $\epsilon=-0.067\pm 0.027$, corresponding to roughly a $2\sigma$ deviation from zero. In contrast, D2 yields a value fully consistent with $\epsilon=0$ within $1\sigma$, whereas D3 prefers a small but statistically significant positive value, $\epsilon\simeq 0.0023$ ($\sim 3.8\sigma$). This indicates that the inclusion of BAO and especially CMB information drives the model closer to (or tightly constrains deviations from) the standard $\Lambda$CDM limit.

The structure growth parameter $\sigma_{8}$ remains stable across datasets, with values around $0.80$ and differences well within $1\sigma$.

The information criteria indicate that D1 provides the lowest AIC and BIC values within the closed R-$\Lambda$CDM scenario, followed by D2, while D3 yields the largest values. Since the dataset compositions differ, these criteria should be interpreted separately within each configuration.

Comparison of dataset D1 with D2 and D3

The comparison between D1 and D2 reveals a significant shift in the effective dark matter density $\Omega_{dm0}^{\rm eff}$ at the $1.73\sigma$ level. The effective baryon density differs by approximately $0.98\sigma$, while the curvature parameter $\Omega_{k0}$ shifts by about $0.54\sigma$. The Hubble constant remains statistically consistent ($0.29\sigma$), and $\sigma_{8}$ differs by only $0.13\sigma$. The radiation density parameter shows a deviation of approximately $2.33\sigma$, and the Rastall parameter $\epsilon$ differs by about $2.40\sigma$. Overall, most parameters remain compatible within $2\sigma$.

A comparison between D1 and D3 indicates moderate deviations in several parameters. The dark matter density differs by approximately $0.85\sigma$, the baryon density by about $0.75\sigma$, and the curvature parameter by roughly $0.96\sigma$. The Hubble constant shifts by approximately $2.64\sigma$, while $\sigma_{8}$ remains essentially unchanged ($0.01\sigma$). The radiation density parameter shows a deviation close to $3.65\sigma$, reflecting its sensitivity to early-Universe constraints. The Rastall parameter differs by about $2.56\sigma$, driven mainly by the negative value preferred in D1.

These results indicate that, although several parameters exhibit moderate shifts when CMB information is included, the overall parameter sets remain statistically compatible within the closed R-$\Lambda$CDM framework.

Model Selection via AIC and BIC

Table 13 presents a comparative assessment of the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) for the closed R-$\Lambda$CDM and standard closed $\Lambda$CDM models. The differences are defined as $\Delta$AIC$=\mathrm{AIC}_{\rm R}-\mathrm{AIC}_{\Lambda}$ and similarly for $\Delta$BIC.

For dataset D1, the R-$\Lambda$CDM model yields $\Delta$AIC$=-1.65$, indicating statistical equivalence with the standard $\Lambda$CDM model according to the conventional $|\Delta\mathrm{AIC}|<2$ criterion. However, the corresponding $\Delta$BIC$=+3.82$ reflects the penalty imposed by the additional Rastall parameter when model complexity is more strongly weighted.

For D2, the AIC difference ($\Delta$AIC$=+1.82$) indicates statistical equivalence between the two models, while the BIC difference ($\Delta$BIC$=+7.31$) provides strong evidence in favor of the standard $\Lambda$CDM model.

A markedly different behavior is observed for D3. In this case, the closed R-$\Lambda$CDM model significantly outperforms the standard $\Lambda$CDM model, with $\Delta$AIC$=-12.59$ and $\Delta$BIC$=-7.10$. According to standard model selection thresholds, such large negative differences correspond to strong statistical evidence in favor of the R-$\Lambda$CDM extension when BAO and CMB data are simultaneously included.

As emphasized above, the absolute values of AIC and BIC are meaningful only within a fixed dataset configuration.

Overall, the information criteria analysis indicates that the statistical viability of the closed R-$\Lambda$CDM model depends sensitively on dataset composition. While it remains statistically equivalent to $\Lambda$CDM in D1 and statistically equivalent under AIC but disfavored by BIC in D2, it receives strong support in D3, highlighting the crucial role of CMB constraints in assessing deviations from the standard cosmological model.

V.2.4 Assessment of geometrical models based on parameter consistency: The R-$\Lambda$CDM case

Table 14 quantifies the internal consistency of the R-$\Lambda$CDM model across different dataset combinations and spatial geometries. The cosmological parameter differences (CPDs), expressed in units of the combined statistical uncertainty, provide a direct measure of the stability of the inferred parameters when progressively incorporating additional observational information.

A common feature across all geometries is that the largest dataset-induced shifts occur in the effective radiation density, $\Omega_{r0}^{\rm eff}$, and in the Rastall deviation parameter, $\epsilon$. In the flat case, $\Omega_{r0}^{\rm eff}$ exhibits a substantial shift of $4.03\sigma$ between D1 and D3, while $\epsilon$ varies at the $\sim 2.3\sigma$ level. A similarly pronounced CMB-driven effect is observed in the Hubble constant, which changes by $3.07\sigma$ between D1 and D3. It is important to note that the shift in $\Omega_{r0}^{\rm eff}$ primarily reflects its implicit dependence on $H_{0}$, rather than a modification of radiation physics itself. These results indicate that, in the flat R-$\Lambda$CDM scenario, the inclusion of early-Universe information significantly affects both the radiation sector and the inferred expansion rate.

In contrast, the open R-$\Lambda$CDM geometry exhibits the most homogeneous internal behavior. No parameter exceeds the $3\sigma$ level for any dataset comparison. The largest deviations are found in $\Omega_{r0}^{\rm eff}$ ($2.50\sigma$ for $|D_{1}-D_{3}|$) and in $\epsilon$ (approximately $2.0\sigma$), while the Hubble constant remains comparatively stable, with a shift of $1.49\sigma$ between D1 and D3. The curvature parameter $\Omega_{k0}$ itself varies at the $\lesssim 1.35\sigma$ level. Overall, the open geometry provides the most internally uniform parameter reconstruction within the R-$\Lambda$CDM framework.

The closed R-$\Lambda$CDM case exhibits intermediate behavior. While the growth amplitude $\sigma_{8}$ is remarkably stable ($0.01\sigma$ between D1 and D3), the effective radiation density still shows a significant $3.65\sigma$ shift, and the Hubble constant varies by $2.64\sigma$ when CMB information is included. The Rastall parameter also displays a persistent $\sim 2.5\sigma$ dataset dependence. Therefore, although the closed geometry maintains strong consistency in the late-time growth sector, it remains sensitive to early-Universe constraints in the radiation and expansion-rate sectors.

Taken together, the CPD analysis indicates that the dominant source of dataset sensitivity in R-$\Lambda$CDM is associated with $\Omega_{r0}^{\rm eff}$ and, through its coupling, the inferred value of $H_{0}$. Among the three geometrical scenarios, the open R-$\Lambda$CDM model exhibits the most homogeneous internal parameter behavior, as it avoids deviations exceeding the $3\sigma$ level while maintaining stable estimates across matter, curvature, and expansion parameters.

V.3 Dependence of $H_{0}$ constraints on dataset combinations

A key aspect in the interpretation of the Hubble tension is the dependence of the inferred value of $H_{0}$ on the choice of observational datasets. In order to systematically investigate this issue, we consider four different dataset combinations, namely D1, D2, D3, and D4, where D4 corresponds to the Planck 2018 CMB distance priors alone.

The results for all cosmological models are summarized in Table 15. A clear and robust trend emerges across all geometries and models. Late-time datasets (D1 and D2), which include supernovae, cosmic chronometers, and large-scale structure data, systematically favor higher values of the Hubble constant, typically in the range $H_{0}\sim 70$–$71\ {\rm km\,s^{-1}\,Mpc^{-1}}$. In contrast, when CMB information is included (D3), or when CMB data are used alone (D4), the inferred value of $H_{0}$ shifts toward the Planck-preferred region, $H_{0}\sim 67$–$68\ {\rm km\,s^{-1}\,Mpc^{-1}}$.

This behavior is consistently observed in both $\Lambda$CDM and R–$\Lambda$CDM cosmologies, as well as for flat, open, and closed geometries. In particular, the D4 dataset provides a direct estimate of $H_{0}$ from early-Universe information alone. For the flat $\Lambda$CDM model, we obtain $H_{0}=67.66^{+0.60}_{-0.61}\ {\rm km\,s^{-1}\,Mpc^{-1}}$, in full agreement with the Planck 2018 baseline result. Similarly, in the flat R–$\Lambda$CDM scenario we find $H_{0}=67.13^{+2.28}_{-1.04}\ {\rm km\,s^{-1}\,Mpc^{-1}}$, indicating that the inclusion of the Rastall parameter enlarges the allowed parameter space but does not significantly shift the central value when only CMB information is considered.

These results highlight that the Hubble tension can be interpreted as a discrepancy between early- and late-Universe probes, rather than a failure of a specific cosmological model. Within this context, modified gravity effects, such as those introduced by Rastall theory, can alter the inferred value of $H_{0}$ depending on the dataset combination.

To further clarify this point, we note that the Rastall parameter $\epsilon$ modifies the evolution of the matter density according to Eq. (25), leading to a deviation from the standard $(1+z)^{3}$ scaling. For $\epsilon>0$, the matter density dilutes faster with redshift, reducing its contribution to the expansion rate at intermediate redshifts. In order to maintain consistency with distance-based observables, this effect is compensated by a higher inferred value of $H_{0}$, thereby generating a positive correlation between $\epsilon$ and $H_{0}$.

This behavior is explicitly confirmed by our results for the CMB-inclusive dataset D3. For example, in the flat case $H_{0}$ increases from $67.35\pm 0.28$ in $\Lambda$CDM to $68.50^{+0.40}_{-0.41}$ in R–$\Lambda$CDM. In the open geometry, the shift is more pronounced, from $67.55^{+0.29}_{-0.28}$ to $69.50^{+0.52}_{-0.49}$, while in the closed case $H_{0}$ increases from $67.18\pm 0.30$ to $68.86^{+0.53}_{-0.54}$. These shifts are summarized in Table 16.

On the other hand, for $\epsilon<0$, the matter dilution becomes slower than in the standard scenario, leading to a relatively enhanced matter contribution at late times. In this case, the same compensating mechanism does not produce a significant increase in $H_{0}$. This is consistent with our D1 results, where mildly negative values of $\epsilon$ are allowed, but the inferred values of $H_{0}$ remain very close to those of the $\Lambda$CDM model.

Overall, these findings demonstrate that the impact of the Rastall parameter on the Hubble constant is not universal, but depends sensitively on both its sign and the dataset combination. In particular, a mild enhancement of $H_{0}$ occurs primarily for positive $\epsilon$ values in the presence of CMB constraints, while this effect is significantly reduced or absent for negative $\epsilon$ or when only late-time data are considered.

VI Conclusions and Outlook

In this work, we have presented a systematic and internally consistent assessment of the Hubble tension within the standard $\Lambda$CDM framework and its Rastall extension (R-$\Lambda$CDM), considering flat, open, and closed spatial geometries and three complementary dataset combinations: D1 (late-time probes), D2 (late-time probes combined with BAO and BBN), and D3 (late-time probes combined with BAO and Planck 2018 CMB distance priors). Our primary objective was to determine whether spatial curvature or matter–geometry coupling can provide a statistically meaningful and dataset-independent alleviation of the discrepancy between early- and late-Universe determinations of $H_{0}$.

Within $\Lambda$CDM, a clear dataset-driven structure emerges across all geometries. Using only late-time probes (D1 and D2), the inferred expansion rate lies in the range $H_{0}\simeq 70.75$–$71.43~{\rm km\,s^{-1}\,Mpc^{-1}}$, corresponding to a reduced SH0ES discrepancy of $1.11\sigma$–$1.63\sigma$, while remaining in $3.62\sigma$–$4.23\sigma$ tension with the Planck 2018 reference value. When CMB distance priors are included (D3), the posterior shifts to $H_{0}\simeq 67.18$–$67.55~{\rm km\,s^{-1}\,Mpc^{-1}}$, restoring agreement with Planck at the $0.09\sigma$–$0.38\sigma$ level but increasing the SH0ES discrepancy to $4.25\sigma$–$4.51\sigma$. Allowing for spatial curvature does not qualitatively modify this redistribution of the tension, and $\Omega_{k0}$ becomes statistically consistent with spatial flatness once early-Universe information is incorporated. These results indicate that curvature alone does not provide a structural resolution of the Hubble tension.

In the R-$\Lambda$CDM framework, the inclusion of the Rastall parameter $\epsilon$ alters the quantitative behavior while preserving the same global trend. For D1 and D2, the inferred Hubble constant lies in the range $H_{0}\simeq 70.79$–$71.48~{\rm km\,s^{-1}\,Mpc^{-1}}$, yielding a SH0ES discrepancy of $1.11\sigma$–$1.60\sigma$ and a Planck tension of $3.64\sigma$–$4.19\sigma$. For D3, the expansion rate decreases to $H_{0}\simeq 68.50$–$69.50~{\rm km\,s^{-1}\,Mpc^{-1}}$, lowering the Planck discrepancy to $1.70\sigma$–$2.95\sigma$ and producing a SH0ES tension of $2.65\sigma$–$3.45\sigma$. In these configurations, the posterior value of $H_{0}$ lies between the Planck and SH0ES determinations, partially reducing the discrepancy with both datasets, but not eliminating it. The Rastall parameter exhibits moderate dataset dependence, favoring a negative value at approximately the $2\sigma$ level in D1, while becoming tightly constrained and statistically compatible with $\epsilon=0$ once BAO and CMB information are included. The corresponding tension levels for all geometries and dataset combinations are summarized in Table 17, highlighting the dataset-driven redistribution of the discrepancy.

We find that the observational constraints, particularly for the CMB-inclusive dataset combinations D3, strongly favor values of the Rastall parameter close to the GR limit $\epsilon=0$. Consequently, although Rastall gravity modifies the effective matter evolution relative to standard $\Lambda$CDM and can in principle affect the late-time clustering behavior, the observationally allowed departures from the standard cosmological scenario remain relatively mild.

The CPD analysis further clarifies the origin of these shifts. We find that the dominant dataset-induced deviations occur primarily in $H_{0}$, whereas most other cosmological parameters remain statistically stable within $\lesssim 1\sigma$ across dataset combinations. This behavior indicates that the Hubble tension is largely driven by calibration differences between early- and late-Universe probes rather than by a global restructuring of the cosmological background dynamics.

Model comparison based on AIC and BIC shows statistical equivalence between $\Lambda$CDM and R-$\Lambda$CDM for D1 and D2. In contrast, within the D3 configuration the R-$\Lambda$CDM extension becomes statistically preferred in flat and closed geometries and decisively favored in the open configuration. This preference arises from the combined inclusion of BAO and CMB distance priors, highlighting the sensitivity of extended gravity frameworks to early-Universe information.

Overall, neither spatial curvature nor simple matter–geometry coupling provides a dataset-independent reconciliation of early- and late-Universe measurements of $H_{0}$. The alleviation of the discrepancy remains strongly dependent on the adopted dataset combination, reinforcing the interpretation that the Hubble tension is fundamentally calibration-sensitive and possibly structural in nature.

Future progress will require a fully consistent joint analysis employing complete CMB likelihoods rather than distance priors alone, combined with next-generation large-scale structure, BAO, and distance-ladder measurements. In parallel, exploring extended gravity scenarios beyond minimal matter–geometry couplings, as well as model-independent consistency tests and hierarchical Bayesian calibration frameworks, will be essential to determine whether the Hubble tension reflects residual systematics or genuinely signals physics beyond the standard cosmological paradigm.

To further clarify both the dataset dependence of the inferred Hubble constant and the physical origin of its shift in the Rastall framework, we have incorporated two complementary summary tables in the revised manuscript. Table 15 provides a global overview of the $H_{0}$ constraints across all cosmological models and dataset combinations, clearly showing that late-time probes (D1 and D2) systematically favor higher values of $H_{0}$, while the inclusion of CMB information in D3 and D4 shifts the results toward the Planck-preferred region.

In addition, Table 16 isolates the CMB-inclusive case and demonstrates that, for small positive values of the Rastall parameter $\epsilon\sim\mathcal{O}(10^{-3})$, the R-$\Lambda$CDM model consistently yields slightly higher values of $H_{0}$ compared to the standard $\Lambda$CDM scenario across all spatial geometries. This behavior can be directly traced back to the modified matter evolution in Rastall cosmology, which induces a degeneracy between $\epsilon$ and $H_{0}$ at the background level.

Taken together, these results show that the apparent enhancement of $H_{0}$ is not a generic feature of a non-zero Rastall parameter, but rather a controlled and dataset-dependent effect that becomes visible primarily when $\epsilon$ is positive and constrained by CMB data, while remaining suppressed when $\epsilon$ is driven toward zero by high-precision observations.

Appendix A Statistical Estimator for Quantifying Tension Between Datasets

A practical approach for evaluating potential statistical inconsistencies between two independent datasets, denoted as Di and Dj, is to compare the marginalized one-dimensional posterior distributions of a given cosmological parameter $\alpha$. To quantify such a discrepancy, a statistical estimator was proposed by LEIZEROVICH2024138844 , which computes the tension in terms of the number of standard deviations separating the parameter means. The estimator is defined as

[TABLE]

where $\mu_{D_{i}}$ and $\mu_{D_{j}}$ represent the mean values of the parameter estimates obtained from the posterior distributions corresponding to datasets Di and Dj, respectively, and $\sigma_{D_{i}}$ and $\sigma_{D_{j}}$ are the associated standard deviations. This method is widely used due to its simplicity, and we have also employed it to calculate and assess the discrepancies in cosmological parameters between different datasets.