Cosmological constraints of Palatini $f(\mathcal{R})$ gravity

Palatini gravity belongs to the category of metric-affine models of gravity, representing the simplest example of this framework. In this approach, the metric tensor $g$ and the connection $\hat{\Gamma}$ are considered as independent entities. Consequently, to derive the field equations governing this gravitational theory, one must perform variations with respect to both of these variables. The action functional for this theory is given as follows:

$$S=\frac{1}{2\kappa^{2}}\int d^{4}x\sqrt{-g}f(\mathcal{R})+S_{m}[g_{\mu\nu},\psi_{m}].$$

(2.1)

Here, $\sqrt{-g}$ represents the determinant of the metric tensor, $f(\mathcal{R})$ is a function of the curvature scalar $\mathcal{R}$, and $S_{m}$ represents the action for matter fields $\psi_{m}$, which exclusively depends on the metric tensor $g_{\mu\nu}$. Importantly, the matter action $S_{m}$ remains entirely independent of the connection $\hat{\Gamma}$. Although it is conceivable to extend the matter action to encompass a dependence on $\hat{\Gamma}$, such an augmentation is deemed pertinent primarily in the context of fermionic interactions as opposed to classical fluid dynamics on a large scale, and thus, we shall omit its consideration in this context.

Furthermore, the curvature scalar $\mathcal{R}$ is constructed utilizing both the metric tensor and the Palatini-Ricci curvature tensor $\mathcal{R}_{\mu\nu}(\Gamma)$. Specifically, it can be expressed as $\mathcal{R}=g^{\mu\nu}\mathcal{R}_{\mu\nu}(\Gamma)$. It’s noteworthy that $\mathcal{R}_{\mu\nu}$ must maintain symmetry to prevent instabilities [87, 88, 89].

Should the functional form of $f$ exhibit linearity with respect to $\mathcal{R}$, we recover GR. In scenarios involving vacuum or spacetime dominated by pure radiation, we effectively deal with GR featuring a cosmological constant. However, when considering $f$ as an arbitrary function of $\mathcal{R}$, we encounter a distinct spacetime structure [90, 91, 54]. Significantly, within our context, these modifications carry implications for the description of stellar phenomena [92] and early and late cosmology [90, 91, 93, 94]. Let us delve deeper into this subject.

Variation of Eq. (2.1) solely with respect to the metric tensor $g_{\mu\nu}$ results in the following field equations:

$$f^{\prime}(\mathcal{R})\mathcal{R}_{\mu\nu}-\frac{1}{2}f(\mathcal{R})g_{\mu\nu}=\kappa^{2}T_{\mu\nu},$$

(2.2)

where $T_{\mu\nu}$ represents the energy-momentum tensor for the matter field and can be derived in the conventional manner:

$$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta S_{m}}{\delta g_{\mu\nu}}.$$

(2.3)

In our scenario, we assume the matter content to be in the form of a perfect fluid. In the aforementioned equations, the primes denote derivatives with respect to the argument of the function, denoted as $f^{\prime}(\mathcal{R})=df(\mathcal{R})/d\mathcal{R}$. Evidently, for linear $f(\mathcal{R})$, these equations reduce to the familiar equations of GR.

Conversely, the variation with respect to the independent connection $\hat{\Gamma}$ can be expressed as [93, 95]:

$$\hat{\nabla}_{\beta}(\sqrt{-g}f^{\prime}(\mathcal{R})g^{\mu\nu})=0.$$

(2.4)

Here, $\hat{\nabla}_{\beta}$ denotes the covariant derivative calculated with respect to $\hat{\Gamma}$. Notably, by introducing a new metric tensor defined as:

$$\hat{g}_{\mu\nu}=f^{\prime}(\mathcal{R})g_{\mu\nu},$$

(2.5)

we can reformulate Eq. (2.4) as follows:

$$\nabla_{\beta}(\sqrt{-\hat{g}}\hat{g}^{\mu\nu})=0,$$

(2.6)

provided that the independent connection $\hat{\Gamma}$ corresponds to the Levi-Civita connection associated with the metric $\hat{g}_{\mu\nu}$. It is also evident that, for linear $f$, the connection $\hat{\Gamma}$ reduces to the Levi-Civita connection $\Gamma$ associated with the metric $g_{\mu\nu}$.

Furthermore, we can gain deeper insights into the characteristics of this gravitational theory by scrutinizing the $g$-trace of Eq. (2.2):

$$f^{\prime}(\mathcal{R})\mathcal{R}-2f(\mathcal{R})=\kappa^{2}T,$$

(2.7)

where $T$ represents the trace of the energy-momentum tensor. This expression unveils the possibility of algebraically expressing $\mathcal{R}$ as a function of matter fields, contingent on a given $f(\mathcal{R})$. In instances featuring vacuum and/or pure radiation-dominated spacetime, where $T=0$, we can solve Eq. (2.7), resulting in the theory converging to GR with a cosmological constant. Additionally, it can be established that this theory does not introduce any supplementary degrees of freedom beyond those already inherent in GR. This contrasts with the case of metric $f(R)$ gravity, where the relationship between $R$ and $T$ arises through a differential equation, thereby introducing an extra degree of freedom, often interpreted as a scalar field, subject to appropriate transformations.

In our subsequent analysis, we will consider the simplest extension of GR, often referred to as the Starobinsky or quadratic model, which is expressed as:

$$f(\mathcal{R})=\mathcal{R}+\alpha\mathcal{R}^{2},$$

(2.8)

where $\alpha$ serves as a parameter. In this scenario, the solution to Eq. (2.7) adopts the form of GR, i.e., $\mathcal{R}=-\kappa^{2}T$. This signifies that any modifications manifesting in the field equations are contingent on matter fields and are parametrized by a solitary parameter, $\alpha$.

As was already mentioned, Palatini $f(\mathcal{R})$ gravity does not introduce any additional degrees of freedom and reduces to GR with a cosmological constant in the vacuum. Due to this fact, Palatini theories find applications in modeling the late-time acceleration of the Universe [54, 96], where curvature functions of type $f(\mathcal{R})=\mathcal{R}+\alpha_{n}\mathcal{R}^{n}$ are of particular interest. It was also found that the choice $f(\mathcal{R})=\mathcal{R}^{n(t)}$, where the exponent depends on the temperature and varies from $n=1$ in the low-energy limit, and $n=2$ at the early Universe, yields a scale-invariant description of gravity in high-energy limit, reduces to GR at later times, and evades Ostroradsky’s instability [97].

In what follows, we will assume the quadratic model in the Palatini approach, i.e. $f(\mathcal{R})=\mathcal{R}+\alpha\mathcal{R}^{2}$, and the standard Friedmann-Lemaître-Robertson-Walker metric:

$$g_{\mu\nu}=\text{diag}\left(-1,\frac{a^{2}}{1-kr^{2}},a^{2}r^{2},a^{2}r^{2}\sin^{2}\theta\right)$$

(2.9)

and the stress-energy tensor for the perfect fluid:

$$T_{\mu\nu}=(p+\rho)u_{\mu}u_{\nu}+pg_{\mu\nu},\quad T=3p-\rho,$$

(2.10)

satisfying the continuity equation $\nabla^{(g)}_{\mu}T^{\mu\nu}=0$, where $\nabla^{(g)}$ is the covariant derivative calculated using the metric tensor $g_{\mu\nu}$. The continuity equation allows us to write:

$$\mathcal{R}=3H_{0}^{2}(\Omega_{m,0}a^{-3}+4\Omega_{\Lambda,0})\,.$$

(2.11)

It must be noted that the two Ricci tensors, Palatini $\mathcal{R}_{\mu\nu}$ and metric $R_{\mu\nu}$, are related to each other via a conformal transformation since $\mathcal{R}_{\mu\nu}=\mathcal{R}_{\mu\nu}(f^{\prime}(T)g)$, so that¹¹1Let us notice that the curvature scalar $\mathcal{R}$ is not the curvature scalar built from the conformal metric only, $\hat{R}(f^{\prime}(T)g)$, but rather it is a contraction of $g^{\mu\nu}$ and the Palatini Ricci tensor $\mathcal{R}_{\mu\nu}(f^{\prime}(T)g)$. [98]

$$\begin{split}\mathcal{R}_{\mu\nu}&=R_{\mu\nu}+\frac{3}{2}\frac{1}{(f^{\prime}(T))^{2}}\partial_{\mu}f^{\prime}(T)\partial_{\nu}f^{\prime}(T)\\ &-\frac{1}{f^{\prime}(T)}\left(\nabla^{(g)}_{\mu}\nabla^{(g)}_{\nu}-\frac{1}{2}g_{\mu\nu}\Box^{(g)}\right)f^{\prime}(T)\,.\end{split}$$

(2.12)

To get the Friedmann equation, one needs to plug relations (2.11) and (2.12) to the field equations (2.2) and get rid of $a^{\prime\prime}$ from one of the formulas. Upon doing so, one arrives at the following equation for the quadratic model [94, 99, 56, 100]:

	$\displaystyle\frac{H^{2}}{H_{0}^{2}}$	$\displaystyle=\frac{b^{2}}{\left(b+\frac{d}{2}\right)^{2}}\left[\Omega_{\alpha}(\Omega_{\text{m},0}a^{-3}+\Omega_{\Lambda,0})^{2}\times\frac{(K-3)(K+1)}{2b}\right.$
		$\displaystyle\left.+(\Omega_{\text{m,0}}a^{-3}+\Omega_{\Lambda,0})+\frac{\Omega_{\text{r},0}a^{-4}}{b}+\Omega_{k}\right],$		(2.13)

where

	$\displaystyle\Omega_{k}$	$\displaystyle=-\frac{k}{H_{0}^{2}a^{2}},$		(2.14)
	$\displaystyle\Omega_{\text{r},0}$	$\displaystyle=\frac{\kappa^{2}\rho_{\text{r,0}}}{3H_{0}^{2}},$		(2.15)
	$\displaystyle\Omega_{\text{m},0}$	$\displaystyle=\frac{\kappa^{2}\rho_{\text{m,0}}}{3H_{0}^{2}},$		(2.16)
	$\displaystyle\Omega_{\Lambda,0}$	$\displaystyle=\frac{\kappa^{2}\rho_{\Lambda,0}}{3H_{0}^{2}},$		(2.17)
	$\displaystyle K$	$\displaystyle=\frac{3\Omega_{\Lambda,0}}{(\Omega_{\text{m,0}}a^{-3}+\Omega_{\Lambda,0})},$		(2.18)
	$\displaystyle\Omega_{\alpha}$	$\displaystyle=3\alpha H_{0}^{2},$		(2.19)
	$\displaystyle b$	$\displaystyle=f^{\prime}(\hat{R})=1+2\Omega_{\alpha}(\Omega_{\text{m,0}}a^{-3}+4\Omega_{\Lambda,0}),$		(2.20)
	$\displaystyle d$	$\displaystyle=\frac{1}{H}\frac{db}{dt}=-2\Omega_{\alpha}(\Omega_{\text{m,0}}a^{-3}+\Omega_{\Lambda,0})(3-K)\,.$		(2.21)

From the above one can check that the standard model (2.22) can be reconstructed in the limit $\alpha\mapsto 0$:

$$\frac{H^{2}}{H_{0}^{2}}=\Omega_{\text{r},0}a^{-4}+\Omega_{\text{m},0}a^{-3}+\Omega_{\Lambda,0}+\Omega_{k}.$$

(2.22)

In what follows, we will skip the spatial curvature term $\Omega_{k}$.

We devote this part of the manuscript to the description of the data sets being used in our analysis. Our baseline data set consists of Hubble expansion data which we then complement with supernova Type Ia (SNIa) and Baryonic acoustic oscillation (BAO) data. We also include the recent release of the SH0ES value for the Hubble constant [101], $H_{0}^{\rm S22}=73.04\pm 1.04\,{\rm km\,s}^{-1}{\rm Mpc}^{-1}$, which has an associated reported absolute magnitude of $M^{\rm S22}=-19.253\pm 0.027\,{\rm mag}$. The other data sets are

• •

  Cosmic Chronometers (CC) – We use the public list of $31$ data points [102, 103, 104, 105, 106, 107, 108] that utilize the CC method. In this approach, a spectroscopic dating technique is used for passively evolving galaxies. In this way, direct measurements are obtainable for the Hubble parameter up to $z\lesssim 2$. CC also has the advantage that it is independent of any cosmological model as well as the the Cepheid distance scale. On the other hand, some assumptions are made on the modeling of stellar ages. For two passively aging galaxies, observations can constrain $\Delta z/\Delta t$ which can inturn be used to estimate the Hubble parameter at that redshift through $H(z)=-(1+z)^{-1}\Delta z/\Delta t$. The corresponding $\chi^{2}_{H}$ estimator is represented by

  $$\chi^{2}_{H}=\sum^{31}_{i=1}\frac{\left(H(z_{i},\Theta)-H_{\mathrm{obs}}(z_{i})\right)^{2}}{\sigma^{2}_{H}(z_{i})}\,,$$

  (3.1)

  where $H(z_{i},\Theta)$ are the theoretical model Hubble parameter values at redshift $z_{i}$ with model parameters $\Theta$ while $H_{\mathrm{obs}}(z_{i})$ are the corresponding Hubble data values at $z_{i}$ with an associated uncertainty of $\sigma_{H}(z_{i})$.

• •

  Pantheon+ Compilation – This compilation is based on SNIa observations that are calibrated using a Cepheid distance scale. These events occur in binary stellar systems with uniform intrinsic brightness making them ideal standard candles in order to measure expansion profilers for distant galaxies. The Pantheon+ (SN) compilation [109] consists of 1701 SNIa samples which reports distance modulus $\mu$ of an object together with the redshift $z$ as measured in the frame of the cosmic microwave background radiation (CMB). The distance modulus depends on the observed apparent magnitude of an object, $m$, together with its absolute magnitude, $M$ which is a measure of intrinsic brightness. For an object at redshift $z_{i}$, the distance modulus is given as

  $$\mu(z_{i},\Theta)=m-M=5\log_{10}[D_{L}(z_{i},\Theta)]+25\,,$$

  (3.2)

  where $D_{L}(z_{i},\Theta)$ is the luminosity distance defined as

  $$D_{L}(z_{i},\Theta)=c(1+z_{i})\int_{0}^{z_{i}}\frac{dz^{\prime}}{H(z^{\prime},\Theta)}\,.$$

  (3.3)

  The SN data set needs to be calibrated for the absolute magnitude $M$, which is treated as a nuisance parameter in the MCMC and is marginalized over. Finally, the cosmological parameters are constrained by minimizing a $\chi^{2}$ likelihood specified by [110]

  $$\chi^{2}_{\mathrm{SN}}=(\Delta\mu(z_{i}),\Theta))^{T}C^{-1}(\Delta\mu(z_{i}),\Theta))\,,$$

  (3.4)

  where $(\Delta\mu(z_{i}),\Theta))=(\mu(z_{i}),\Theta)-\mu(z_{i})_{\mathrm{obs}}$ and $C$ is the corresponding covariance matrix which accounts for the statistical and systematic uncertainties. The SN compilation used in this analysis incorporates the SH0ES Cepheid host distance anchors ($H_{0}^{\rm S22}$) [101] in the likelihood which helps to break the degeneracy between the parameters $M$ and $H_{0}$ when analyzing SNIa alone. The redshift range of the SN compilation is $0.01<z<2.5$ which includes much more information at lower redshift as compared with previous SNIa compilations.

• •

  Baryon Acoustic Oscillations - Our analysis is nuanced by a joint BAO data set consisting of separately independent data points. This includes measurements from the SDSS Main Galaxy Sample at $z_{\mathrm{eff}}=0.15$ [111], the six-degree Field Galaxy Survey at $z_{\mathrm{eff}}=0.106$ [112], and the BOSS DR11 quasar Lyman-alpha measurement at $z_{\mathrm{eff}}=2.36$ together with the BOSS DR10 measurements at $z_{\mathrm{eff}}=\{0.32,0.57\}$ [113]. We also incorporate the angular diameter distances and $H(z)$ measurements of the SDSS-IV eBOSS DR14 quasar survey at $z_{\mathrm{eff}}=\{0.98,1.23,1.52,1.94\}$ [114], along with the SDSS-III BOSS DR12 consensus BAO measurements of the Hubble parameter and the corresponding comoving angular diameter distances at $z_{\mathrm{eff}}=\{0.38,0.51,0.61\}$ [115]. These two BAO data sets include covariance matrix information which is included in our MCMC analysis. The BAO data sets involve various measures of cosmic expansion include the Hubble distance $D_{H}(z)$, comoving angular diameter distance $D_{M}(z)$, and volume-average distance $D_{V}(z)$, which are respectively defined through

  	$\displaystyle D_{H}(z)$	$\displaystyle=\frac{c}{H(z)}\,,$
  	$\displaystyle D_{M}(z)$	$\displaystyle=(1+z)D_{A}(z)\,,$
  	$\displaystyle D_{V}(z)$	$\displaystyle=\left[(1+z)^{2}D_{A}(z)^{2}\frac{cz}{H(z)}\right]^{1/3}\,,$		(3.5)

  and where $D_{A}(z)=(1+z)^{-2}D_{L}(z)$ is the angular diameter distance. In accordance with the BAO results, we calculate the combination of parameters $\mathcal{G}(z_{i})=D_{V}(z_{i})/r_{s}(z_{d})$, $r_{s}(z_{d})/D_{V}(z_{i}),$ $D_{H}(z_{i}),$ $D_{M}(z_{i})(r_{s,\mathrm{fid}}(z_{d})/r_{s}(z_{d})),$ $H(z_{i})(r_{s}(z_{d})/r_{s,\mathrm{fid}}(z_{d})),$
  $D_{A}(z_{i})(r_{s,\mathrm{fid}}(z_{d})/r_{s}(z_{d}))$ for which the comoving sound horizon at the end of the baryon drag epoch, at redshift $z_{d}\approx 1059.94$, has a fiducial value of $z_{d}\approx 1059.94$ [116]. This lets us define the corresponding $\chi^{2}$ for the BAO data as

  $$\chi^{2}_{\text{BAO}}(\Theta)=\Delta G(z_{i},\Theta)^{T}C_{\text{BAO}}^{-1}\Delta G(z_{i},\Theta)$$

  (3.6)

  where $\Delta G(z_{i},\Theta)=G(z_{i},\Theta)-G_{\text{obs}}(z_{i})$and $C_{\text{BAO}}$ is the covariance matrix of all the considered BAO observations.

Finally, we include statistical criterions through which to assess the performance of the model against $\Lambda$CDM for the various data sets as well as the inclusion of the $H_{0}^{\rm S22}$ prior. To do this, we use the Akaike Information Criteria (AIC) [117] which includes the best performing minimum $\chi^{2}_{\mathrm{min}}$ values and a measure of the model complexity, through the number of parameters $n$. This goodness-of-fit measure is then defined as

$$\mathrm{AIC}=\chi^{2}_{\mathrm{min}}+2n\,.$$

(3.7)

AIC penalizes high numbers of parameters and so lower values of AIC are indicative of better performing models.

We also examine the Bayesian Information Criterion (BIC) which nuances the AIC by putting more emphasise on the level of complexity that a model exhibits. BIC is defined as

$$\mathrm{BIC}=\chi^{2}_{\mathrm{min}}+n\ln m\,,$$

(3.8)

where $m$ is the sample size of the observational data combination. The goals of both parameters is the same but BIC helps penalize higher parameters models in a stronger way. Absolute values of these statistical measures are difficult to relate to the real performance of individual models. For this reason, we consider the difference in each criterion for each individual data set for both the model under consideration together with the MCMC analysis of $\Lambda$CDM.

In this section, we consider how the quadratic Palatini model fairs with observational data through an MCMC constraint analysis. We do this for a combination of all the data sets described previously with CC being our baseline for all analyses. We also perform these analyses again with the $H_{0}^{\rm S22}$ prior on the Hubble constant, as well as including statistical criterions in order to compare and contrast the different performances.

Considering the Friedmann equation in Eq. (2.13), we can fit for the set of free parameters that make up this model. We assume a flat FLRW geometry to avoid additional degeneracies. By evaluating this relation at current times, we find that

	$\displaystyle\Omega_{\Lambda}=\Big{[}-1-4\Omega_{\alpha}(-2+\Omega_{m,0})\pm$
	$\displaystyle\sqrt{1+\Omega_{\alpha}^{2}\left(64-96\Omega_{m,0}\right)-8\Omega_{\alpha}\left(-2+3\Omega_{m,0}+4\Omega_{r,0}\right)}\Big{]}/16\Omega_{\alpha}\,,$		(3.9)

which has a healthy $\Lambda$CDM as $\Omega_{\alpha}\rightarrow 0$. The analytical nature of the solution in Eq. (2.13) makes it very fast, and more able to be incorporated into precision calculations. We show the constraints for the full set of parameters, namely $\{H_{0},\Omega_{m,0},M,\Omega_{\alpha}\}$ in Fig. 1. The high level of consistency in the posterior of the $\Omega_{\alpha}$ shows promising results for the low-energy Palatini model. On the other hand, the parameter seems to be degenerate with the other parameters to a certain degree, except for the CC+SN+BAO case where the contours are rather small. For the other parameters, similar relationships are observed as in the standard cosmological model case. For instance, the nuisance parameter is pushed to lower values when BAO data is incorporated, while an anti-correlation is observed between the Hubble constant and the matter density parameter. The matter density parameter is also very consistent across the different data sets, while the posteriors of the Hubble constant show a mild tension between the joint data sets.

The outputs of the MCMC constraint analysis are shown in Table 1, while the corresponding $\Lambda$CDM fits are contained in Appendix A. One notes that higher Hubble constants are obtained in the quadratic Palatini cases except when BAO data is included. On the other hand, smaller matter density parameters are obtained except for the baseline case. As for the $\Omega_{\alpha}$ parameter, the observational data appears to prefer values that are in the $1-2\sigma$ range as compared with the $\Lambda$CDM value where this parameter is zero. Indeed, for the case where the statistical distance is largest from $\Lambda$CDM, namely for the CC+BAO data set, the statistical criterions give a negative value indicating that the model is slightly preferred, whereas the other cases show a marginal preference for the standard cosmological model.

We consider the case of adding the $H_{0}^{\rm S22}$ prior to all the data sets under consideration in Fig. 2 to assess the consistency of the analyses with this value of the Hubble constant reported in the literature. The introduction of the prior makes the absolute magnitude and Hubble constant anti-correlation more pronuanced while the consistency of the matter density parameter is slightly perturbed. On the other hand, the posterior values of the $\Omega_{\alpha}$ parameter are generically more consistent with $\Lambda$CDM in this scenario.

The outputs of the MCMC analyses of the Palatini model for the $H_{0}^{\rm S22}$ prior case are shown in Table 2. Generically, we find higher values of the fits for the Hubble constant, as one would expect, but we also find lower values of the matter density parameters with reasonable uncertainties, which is consistent with how one would expect the cosmology to evolve. On the other hand, the precise fits for the $\Omega_{\alpha}$ parameter are consistent with $\Lambda$CDM but there is a preference for a nontrivial $\Omega_{\alpha}$. By and large the statistical performance of the model against the different data sets remains largely the same. The superior performance of the model for the instance where we consider CC+BAO data set is considered remains the case. The various fits for the model parameters and the different data set combinations are shown in Fig. 3 which contains a whisker plot of the final constraints together with an indicative region for the $H_{0}^{\rm S22}$ prior.