DE models with combined $H_0 \cdot r_d $ from BAO and CMB dataset and friends

It has been theorized that Dynamical Dark Energy (DDE) could be a possible solution to the Hubble tension. To avoid the degeneracy between the Hubble parameter $H_0$ and the sound horizon scale $r_d$, in this article we use their multiplication as one parameter $c/\left(H_0 r_d\right)$ and we use it to infer cosmological parameters for 6 different models - $\Lambda$CDM and 5 DDE parametrizations -- the Chevallier-Polarski-Linder (CPL), the Barboza-Alcaniz (BA), the Low correlation (LC), the Jassal-Bagla-Padmanabhan (JBP) and the Feng-Shen-Li-Li model. We choose a dataset that treats this combination as one parameter, that includes the Baryon Acoustic Oscillation (BAO) data $0.11 \le z \le 2.40$ and additional points from the Cosmic Microwave Background (CMB) Peaks ($z \approx 1090$). To them, we add the marginalized Pantehon dataset and GRB dataset. We see that the tension is moved from $H_0$ and $r_d$ to $c/\left(H_0 r_d\right)$ and $\Omega_m$. There is only one model that satisfies the Planck 2018 constraints on both parameters and this is LC with a huge error. The rest cannot fit into both constraints. $\Lambda$CDM is preferred with respect to the statistical measures.


Introduction
The quest for understanding cosmological tensions has driven research for years now. It seems that the tension between the direct measurements of the Hubble constant (H 0 ) from the late universe ( [1][2][3][4][5]) and that from the early universe (i.e., from measuring the temperature and polarization anisotropies in the cosmic microwave background [6][7][8][9]) is only aggravated with the increase in the precision and knowledge of the systematics of the data and has reached 5 σ. This discrepancy has spurred many works, trying to resolve whether the dark energy is a constant energy density or with a dynamical behavior, and if so, of what origin, leading to many different theories and possible explanations [10][11][12][13][14][15][16][17][18][19][20][21].
There are many dark energy (DE) parametrizations [22][23][24][25] that can be used in the search for deviations from the cosmological constant, Λ. Some of them fall in the group of early dark energy models [26][27][28][29][30], which modify physics of the early universe. Others modify the latetime universe physics, such as in the phantom dark energy [31,32] models, and emergent dark energy [33,34], add interaction in the DE sector, as in the interacting dark energy [35][36][37], or add exotic species or scalar fields [38][39][40][41][42]. For a review on the taxonomy of DE models, see Refs. [43][44][45]. Finally, there is the generalized emergent dark energy (GEDE) [45] found to be able to compete with ΛCDM for some BAO datasets [46]. In this work, we take the so-called dynamical dark energy parametrizations, which allow for a non-constant DE contribution, regardless of the origin behind it. We take a number of models, namely the Chevallier-Polarski-Linder (CPL), the Barboza-Alcaniz (BA), the low correlation (LC), the Jassal-Bagla-Padmanabhan (JBP) and the Feng-Shen-Li-Li (FSLLI) model, and we use statistical measures to judge their performance in fitting the data.
An important part of the tensions debate revolves around the role of the sound horizon at drag epoch r d . At recombination, after the onset of CMB at z * 1090, the baryons escape the drag of photons at the drag epoch, z d 1059 (Planck 2018 [7]). This sets the standard ruler for the baryon acoustic oscillations (BAO)-the distance (r d ) at which the baryon-photon plasma waves oscillating in the hot universe froze at z = z d . The sound horizon at drag epoch is given by where c s ≈ c(3 + 9ρ b /(4ρ γ )) −0.5 is the speed of sound in the baryon-photon fluid with the baryon ρ b (z) and the photon ρ γ (z) densities, respectively [47,48]. Many papers discuss the relation between the H 0 and the sound horizon scale r d for different models [49][50][51]. Any DE model claiming to resolve the H 0 tension should also be able to resolve the r d tension since they are strongly connected [51][52][53]. In other words, setting a prior on r d has a very strong effect on H 0 and vice versa. In this paper, to avoid this problem, we combine the H 0 and r d into one parameter. We choose measurements that combine the H 0 and r d from the BAO and the prior distance from the CMB peaks [54][55][56][57][58][59][60][61][62] and we use them to infer the cosmological parameters for ΛCDM and 5 DDE models. To the BAO+CMB dataset, we add the gamma-ray bursts (GRB) dataset and the Pantheon dataset with similarly marginalized dependence on H 0 (and M B ). We do this to expand the redshift considered by the models. In a previous work [46], we used a similar approach in which we integrated H 0 r d in the χ 2 of the model, while here, we use them as one single quantity without modifying the χ 2 . In the marginalized version, we saw an interesting possibility for some DE model to fit the data better than ΛCDM. We continue this investigation with new models and a new approach in this paper.
Historically, the approach of using the combination H 0 r d is not new. It has been used in [63] with BAO and SN data to find consistency with the Planck 2015 best-fit ΛCDM cosmology; Ref. [64] used the BAO data to fit the growth measurement, again finding consistency with the Planck 2015; Ref. [65] used the Cepheids and the Tip of the Red Branch measurements to calibrate BAO and SN measurements and find significant tension in both H 0 and r d , despite testing the ΛCDM and DE models (EDE, wCDM, pEDE). The implication is that modifications of the physics after recombination fail to solve both tensions. The overall conclusion is that the H 0 tension should not be considered separately from the r d measurement implied by it [66]. In the current work, we choose a different approach. We repeat the analysis on H 0 r d used in earlier works, but we also take the ratio r * /r d as an independent parameter. This means that we do not use the known analytical formulas for them, but instead we use MCMC to infer them. This avoids using explicit prior knowledge on the baryon load of the universe. This way, we avoid both the degeneracy on H 0 r d from the BAO data, but also we do not use as a hidden prior the Planck measurements.
The plan of the work is as follows: Section 2 formulates the relevant theory. Section 3 describes the method. Section 5 shows the results, and Section 6 summarizes the results.

Theory
A Friedmann-Lemaître-Robertson-Walker metric with the scale parameter a = 1/(1 + z) is considered, where z is the redshift. The evolution of the universe for it is governed by the Friedmann equation, which connects the equation of the state for the ΛCDM background: where in standard ΛCDM, Ω DE (z) → Ω Λ , with the expansion of the universe E(z) = H(z)/H 0 , where H(z) :=ȧ/a is the Hubble parameter at redshift z, and H 0 is the Hubble parameter today. Ω r , Ω m , Ω DE and Ω k are the fractional densities of radiation, matter, dark energy and the spatial curvature at redshift z = 0. We take into account the radiation energy density as Ω r = 1 − Ω m − Ω Λ − Ω k . The spatial curvature is expected to be zero for a flat universe, Ω k = 0, and we set it to zero because we focus on DE models. We will consider a number of different DE models, all of which will feature a dark energy component depending on z. This can be done with a generalization of the Chevallier-Polarski-Linder (CPL) parametrization [67][68][69]: which allows for three possible models from which we will consider only the CPL: and ΛCDM is recovered for w 0 = −1, w a = 0.
To this parametrization, we add another model [43,70], which is the Barboza-Alcaniz (BA) model with This model is good for describing the whole universe history because it does not diverge for z → −1. It gives Next, we use the low correlation model (LC) [43,71] with where w 0 = w(0) and w c = w(z c ) where z c is the redshift at which w 0 and w z are uncorrelated. The effective entry into the EOS is where, here, are replaced w c with w a for consistency with the other models. The Jassal-Bagla-Padmanabhan (JBP) parametrization [44,72] which gives with w 0 = w(z = 0) and w 1 = (dw/dz) |(z=0) .
Finally, we will also test the Feng-Shen-Li-Li parametrization [44,73] which is divergencefree for the entire history of the universe. It has two cases: with the final contribution to the EOS of each of them being, accordingly, In this work, the plus case (i.e., Ω + DE ) is denoted as FSLLI, and the minus case (i.e., Ω − DE ) is denoted as FSLLII.
The distance priors provide effective information of the CMB power spectrum in two aspects: the acoustic scale l A characterizes the CMB temperature power spectrum in the transverse direction, leading to the variation of the peak spacing, and the "shift parameter" R influences the CMB temperature spectrum along the line-of-sight direction, affecting the heights of the peaks. The popular definitions of the distance priors are [74] where z * is the redshift at the photon decoupling epoch with z * 1089 according to the Planck 2018 results [7]. r * is the co-moving sound horizon at z = z * . Ref. [57] derives the distance priors in several different models using Planck 2018 TT,TE,EE + lowE which is the latest CMB data from the final full-mission Planck measurement [7]. We use the correlation matrices given in Table 1 in [57] to obtain the covariance matrices for l A and R corresponding to each model. The angular diameter distance, D A , needed for both the distance priors and the BAO points, is given by where sinn(x) ≡ sin(x), x, sinh(x) for Ω k < 0, Ω k = 0, Ω k > 0, respectively. We see that for the measured D A /r d , one can isolate the variable b = c/(H 0 r d ). Below, we set Ω k = 0, so this formula simplifies to Finally, for the SN and GRB datasets, we define the distance modulus µ(z), which is related to the luminosity distance ( where d L is measured in units of Mpc, and M B is the absolute magnitude.

Methods
In this paper, we use three datasets, which we treat differently. For the BAO dataset, the definition of χ 2 , which we minimize, is the standard one since we do not use the covariance matrix for it.
where v obs is a vector of the observed points (i.e., the values of D A /r d at each z in Table A1), v model is the theoretical prediction of the model calculated with Equation (16) and σ is the error of each measurement. Additionally, we use the SN and the GRB datasets to further constrain the models. For them, we use the following marginalized over H 0 and M B formula, taken from [46] so that we avoid setting priors on H 0 and M B .
Following the approach used in ( [75][76][77][78]), the integrated χ 2 is for where µ i is the observed luminosity, σ i is its error, d L (z) is the luminosity distance, ∆µ = µ i − 5 log 10 [d L (z i )), E is the unit matrix, and C −1 cov is the inverse covariance matrix of the dataset. For the GRB dataset, C −1 → 1/σ 2 i since there is no known covariance matrix for it. For the Pantheon dataset, the total covariance is defined as C cov = D stat + C sys , where D stat = σ 2 i comes from the measurement and C sys is provided separately [79]. Note, in the so-defined marginalized χ 2 , the values of M and H 0 do not change the marginalizedχ 2 SN . The final χ 2 is
To estimate the possible correlations in the BAO dataset, we use the methodology in [90,91]. This method avoids the use of N-body mocks to find the covariance matrices due to systematic errors and replaces it with an evaluation of the effect of possible small correlation on the final result. We add to the covariance matrix for uncorrelated points C ii = σ 2 i symmetrically a number of randomly selected nondiagonal elements C ij . Their magnitudes are set to C ij = 0.5σ i σ j , where σ i σ j are the published 1σ errors of the data points i, j. We introduce positive correlations in up to 6 pairs of randomly selected data points (more than 25% of the data). Figure A1 in Appendix A shows the corner plots with different randomized points for all the models we employ in this article. From the plots, one can see that the effect from adding the correlations is below 10% on average. This indicates that we can consider the chosen set of BAO points for being effectively uncorrelated.
To run the inference, we use a Monte Carlo Markov Chain (MCMC) nested sampler to find the best fit. We use the open-source package Polychord [92] with the GetDist package [93] to present the results.
The prior is a uniform distribution for all the quantities: Ω m ∈ [0, 1.], 35], w 0 ∈ [−1.5, −0.5] and w a ∈ [−0.5, 0.5]. Since the distance prior is defined at the decoupling epoch (z * ) and the BAO-at drag epoch (z d ), we parametrize the difference between r s (z * ) and r s (z d ) as rat = r * /r d , where the prior for the ratio is rat ∈ [0.9, 1.1].  Table A3 in the Appendix A, where also the corner plots can be found. We see that the models differ seriously in their estimations for the physical quantities c/(H 0 r d ), Ω m and r d /r s , probably due to the very wide prior imposed on Ω m . Since we avoid the degeneracy between r d and H 0 by considering the combined quantity c/(H 0 r d ) this leads to an explicit correlation with Ω m for some models and rather strict bounds on the error. The values of Ω m closest to the ones published by Planck 2018 [7] Ω m = 0.315 ± 0.007 are for the BA, JPB and FSLLI models for the BAO dataset and ΛCDM 1 and BA for the BAO+SN+GRB. The rest significantly overestimate Ω m . For the ratio r * /r d Planck 2018 gives 0.98, the closest models are BA, JPB and FSLLI/FSLLII models for the BAO dataset and (flat) ΛCDM, JPB and FSLLI/FSLLII for the BAO+SN+GRB. For c/(H 0 r d ), the Planck 2018 values is 30.26 ± 0.06. Here, the models closest to this value are ΛCDM, CPL and LC for the BAO dataset and ΛCDM, CPL and LC for the BAO+SN+GRB.

Results
The DE parameters seem to be constrained to different level for the different models. As a whole, the trend to better constrain w 0 than w a , which we observed in [46] (and the referenced inside other works), is confirmed in this case as well. Notable exceptions are the BA and LC models, where the error of w a is much smaller. For them, however, the other parameters seem to be outside of the expected boundaries. ΛCDM performs as expected under both datasets.
To compare the different models, we use well-known statistical measures. The results can be seen in Table 1. In it, we publish four selection criteria: Akaike information criterion (AIC), Bayesian information criterion (BIC), deviance information criterion (DIC) and the Bayes factor (BF). Since for small datasets, both AIC and BIC are dominated by the number of parameters in the model (which are 3 for ΛCDM, and 5 for the DE models), we emphasize here on the DIC and the BF which rely on the numerically evaluated likelihood and evidence, making them more unbiased. The DIC criterion, just like the AIC, selects the best model to be the one with the minimal value of the DIC measurement. The reference table we use for DIC is ∆DIC > 10 shows strong support for the model with lower DIC , ∆DIC = 5-10 shows substantial support for the model with lower DIC, and ∆DIC < 5 gives ambiguous support for the model with lower DIC. Here, we use the logarithmic scale for the BF, for which ln(BF) > 1 shows support for the base model (ΛCDM), while ln(BF) < −1 for the other hypothesis. |ln(BF)| < 1 shows an inconclusive result. From Table 1, we see that the AIC and BIC for all models show a preference for ΛCDM. For the DIC criterion, we see a slight possibility for a preference for other models in the case of the CPL and BA models for both tested datasets. For the BF, we see that there is some possible preference for BA, JPN and FSLLI/FSLLII for the BAO+CMB case and for CPL and BA, JPN and FSLLI/FSLLII in the BAO+CMB+SN+GRB case. The results of the LC model show that it is underfitting the data (from the χ 2 /do f ∼2) and the statistics for it is not reliable. This demonstrates another benefit of performing the statistical analysis.
The preference for the BA and LC models which we observe was also observed in the results of [43], where the authors studied a dataset consisting of SN, cosmic chronometers and gravitational waves.
The BAO dataset we use combines the H 0 and the r d into one quantity. Therefore, we estimate the new variable c/(H 0 r d )∼30. Figure 2 shows the values of the c/(H 0 r d ) for different models vs. the result from Planck 2018: 30.24 ± 0.08. For comparison, the most recent local measurement by SH0ES is 30.19 ± 0.53, corresponding to H 0 = 73.01 ± 0.99 km s −1 Mpc −1 [5]. We do not put it on the plot, because the r d used to obtain it is the indirect result from inference on the H0LiCOW+SN+BAO+SH0ES dataset [65]. It is, however, clearly very close to the Planck value, as expected.
On Figure 2, we superimpose the BAO+CMB-only result with the BAO+CMB+SN+GRB one. This figure enables us to visually track the tension between the Planck 2018 results and the datasets we use, which are mostly local universe ones (except for the 2 CMB points). We see that the tension is now between c/(H 0 r d ) and Ω m . The models whose bounds cross with the Planck 2018 one for c/(H 0 r d ) are ΛCDM, CPL and LC for BAO+CMB and only LC for the BAO+CMB+SN+GRB dataset. For Ω m , the models that enter the interval are all but ΛCDM and CPL for the BAO+CMB dataset and ΛCDM, BA, LC for the BAO+CMB+SN+GRB dataset. We see that the inclusion of the new datasets decreases the number of models satisfying the constraints. The only model that is not in tension is LC because of its huge error. Notably, in this approach, ΛCDM, while satisfying the bounds for c/(H 0 r d ), does not satisfy them for Ω m . From the plot, one can see that in general adding the new datasets decrease the errors, but they do not move the mean values in the same direction and the overall effect is not very big. This may be due to unknown errors in the SN+GRB dataset or to the fact that this dataset is not sensitive toward the combined variable c/(H 0 r d ) since we have marginalized over H 0 so that we do not have to impose a prior on r d . Because of this, the only effect the SN+GRB dataset has on the combined variable is indirect, through Ω m and the other parameters. It could also point to some inconsistency in the µ(M B ) relation such as the ones considered in [95][96][97][98][99][100]) questioning the assumption that M B = const.

Discussion
This paper uses the combination H 0 · r d to avoid the degeneracy between H 0 and r d which has plagued the use of BAO measurements and could be part of the resolution of cosmological tensions. The use of a combined parameter avoids imposing separate priors on H 0 and r d and thus it avoids additional assumptions on them. We use points from the late universe, the BAO dataset (z < 2.4)), few points from the early universe (the CMB distant priors, (z 1089)), to which we add SN data and GRB datasets, properly marginalized, to make a statistical comparison between different DE models.
The results show that the tension is now between the new parameter c/(H 0 r d ) and Ω mthe only model that fits in the constraints set by Planck 2018 is LC, which comes with the biggest error. For the rest of the models, one of the two parameters do not fit the constraints, even if some of them somewhat reduce the tension. Statistically, there is a preference for the ΛCDM model over the DE models in most cases. It is worth noting that there is strong evidence in support of ΛCDM compared to all other models only when using AIC and BIC, while from DIC and BF, the support is not substantial, and it even slightly favors other models. This result raises the question of the use of different statistical measures when comparing DE models, and also it opens the possibility that a better DE model may eventually help in reducing both the H 0 tension and the r d tension.
Another interesting point is that for some models, the known impossibility to constrain w a is eliminated and w a has very tight bounds. These models, LC and BA and somewhat FSLLI, show interesting new possibilities for DE models. Furthermore, the choice of datasets and models make explicit the degeneracy between H 0 · r d and Ω m , emphasizing the need to find a way to disentangle the three quantities-H 0 , r d and Ω m -if we are to understand the cosmological tensions. The results show that adding the SN and GRB datasets decrease the errors on the constrained parameters, but they do not move them in the same direction for each model. We see that combining different datasets and different marginalization techniques, along with the use of statistical measures, is a promising tool to study new cosmological models.
Funding: Bulgarian National Science Fund research grants KP-06-N58/5/19.11.2021. Data Availability Statement: All the data we used in this paper were taken from the corresponding citations and available to use.
Acknowledgments: D.S. thanks David Benisty for the useful comments and discussions. D.S. is thankful to Bulgarian National Science Fund for support via research grants KP-06-N58/5.

Conflicts of Interest:
The authors declare no conflict of interest.   In Sections 5 and 6 we discuss only the flat ΛCDM model. The effect of the spatial curvature on DE models has been considered recently in [94].