Tests of Standard Cosmology in Hořava Gravity, Bayesian Evidence for a Closed Universe, and the Hubble Tension

Standard cosmology, which is usually formulated as the Lambda Cold Dark Matter (LCDM) model, is based on general relativity (GR) with a positive cosmological constant Frie:1922 ; Lema:1927 ; Ries:1998 and has been quite successful in describing the observational data WMAP:2010 . However, with the increased accuracy of data, the significant deviations from LCDM are becoming clearer Ries:2018 ; Ries:2019 ; Asga:2019 . In particular, regarding the recent discrepancies of the Hubble constant from the Cosmic Microwave Background (CMB) data and the direct (local) measurements at the lower redshift, which corresponds to the mismatches between the early and late universes, there have been various proposals (for recent reviews, see DiVa:2021 ; Peri:2021 ; Scho:2021 ) but it is still a challenging problem to find a resolution at the fundamental level.

On the other hand, if our universe was created from the Big Bang, we need a quantum gravity to describe the early universe or later space-time. But, it has been well known that a renormalizable quantum gravity can not be realized in GR or its (relativistic) higher-derivative generalizations, due to the ghost problem Stelle:1976 . Recently, Hořava has proposed a renormalizable, higher-derivative, Lorentz-violating quantum gravity model without the ghost problem, due to the high-energy (UV) Lorentz-symmetry violations from the different scaling dimensions for space and time à la Lifshitz, DeWitt, and Hořava Lifs:1941 ; DeWi:1967 ; Hora:2009 . In the last 12 years, there have been many works on its various aspects (see Wang:2017 for a review and extensive literatures). Theoretically, there are still several fundamental issues, like the full/complete symmetry structure, full dynamical degrees of freedom, renormalizability, and the very meaning of black holes and Hawking radiation, etc. (see Deve:2021 for related discussions and their current status). However, phenomenologically, Hořava gravity is one of the (viable) modified gravity theories and it can be tested from astrophysical or cosmological observations Olmo:2019 . In particular, from the recent detections of gravitational waves from black holes/neutron stars, the importance of quantum gravity and its non-GR behaviors is increasing LIGO:2021 . There are some interesting results on testable Hořava gravity effects, such as the increased maximum mass of neutron stars Kim:2018 , reduced light deflection Liu:2010 and black hole shadow Li:2021 , etc. There are also some constraints on its low-energy limit or Einstein-Aether theory from astrophysical data Emir:2017 ; Gong:2018 ; Khod:2020 ; Gupt:2021 or cosmological data with the “assumed” spatially-flat universe in standard cosmology Frus:2015 ; Frus:2020 . But still, for a renormalizable Hořava gravity with the desired higher-derivative terms, there are no systematic and significant constraints from the observational data.

In this paper, we test the spatially non-flat universe in standard cosmology in the context of Hořava gravity. A peculiar property of the cosmology based on Hořava gravity is that the spatially-curved universe may be more “natural” due to contributions from higher spatial derivatives. For the spatially-flat universe, on the other hand, the usual FLRW background cosmology in Hořava gravity Frie:1922 ; Lema:1927 is the same as in GR and hence there are no observable differences in the data analysis which means the same (background) tensions in LCDM model. In other words, Hořava gravity is a “natural laboratory” for the tests of the spatially non-flat universe in standard cosmology.

Recently, it has been found DiVa:2019 ; Hand:2019 ; DiVa:2020 ; eBOS:2020 that the tensions get worse in LCDM as $-\Omega_{k}$ increases, i.e, a rather lower value of Hubble constant for a closed universe ($\Omega_{k}<0$), which is preferred in the recent Planck CMB data Planck:2018 , without being combined with lensing and Baryon Acoustic Oscillations (BAO). In this paper, we show that the situation in Hořava gravity is the opposite, i.e., tensions get better with an increasing $-\Omega_{k}$, due to some non-linear corrections from (spatial) higher-derivative terms. From an MCMC analysis with the observational data, including BAO, we obtain (a) preference of a closed universe and (b) reduction of the Hubble tension for our two Hořava-gravity based cosmological models A and B, depending on the choice of parametrization, with “very strong” and “strong” Bayesian evidences, respectively, against the standard, spatially-flat, LCDM cosmological model.

To this ends, we consider the ADM (Arnowitt-Deser-Misner Arno:1962 ) metric

and the Hořava gravity action with $z=3$, à la Lifshitz, DeWitt, and Hořava Lifs:1941 ; DeWi:1967 ; Hora:2009 , given by (up to boundary terms)

which is viable ³³3At the perturbation level, even for the spatially-flat background, Hořava gravity produces notable differences from GR. In order to obtain a (nearly) scale-invariant CMB power spectrum for the spatially-flat universe in Hořava gravity, where the “inflation without inflation era” is possible Muko:2009 ; Gesh:2011 , we need a proper form of the six-derivative (UV) terms which break the detailed balance condition in UV as in the action (2) Shin:2017 . It would be interesting to see the curvature-induced effect on the scale-invariant power spectrum. However, in this paper, we consider the action (2) without any UV conditions for the generality of our approach. Shin:2017 and power-counting renormalizable Viss:2009 , with the extrinsic curvature

(the dot $(\dot{~{}})$ denotes a time derivative), the Ricci tensor $R_{ij}$ of the (Euclidean) three-geometry, their corresponding traces $K=g_{ij}K^{ij},R=g_{ij}R^{ij}$, coupling constants ⁴⁴4One might consider extension terms which depend on $a_{i}=\partial_{i}N/N$ and $\nabla_{j}a_{i}$ generally Blas:2009 . However, since these terms considerably affect the IR physics compared to those in the standard action (2) Deve:2018 ; ONea:2020 , we do not consider those terms in this paper by assuming that the IR physics is well described by GR. This implies that the gravity probe $E_{G}$ approaches the GR prediction in the current epoch, i.e., low $z$, though currently it would not be tested due to large statistical errors Zhan:2020 . $\kappa,\lambda,\xi,\alpha_{i}$, and a cosmological constant parameter ${\sigma}$.

In order to study standard cosmology for the Hořava gravity action (2), we consider the homogeneous and isotropic Friedmann-Lemaitre-Robertson-Walker (FLRW) metric ansatz

with the (spatial) curvature parameter $k=+1,0,-1$ for a closed, flat, open universe, respectively, and the curvature radius $R_{0}$ in the current epoch $a(t_{0})\equiv 1$. Assuming the perfect fluid form of matter contributions with energy density $\rho$ and pressure $p$, we obtain the Friedmann equations as

where we have used the conventional parametrization of the coupling constants $\sigma,\xi,{\alpha}_{1},{\alpha}_{2}$ for the lower-derivative terms Hora:2009 ; Park:2009

with an IR-modification parameter ${\omega}$ Park:2009 , $\mu^{2}>0~{}(<0)$ for a positive (negative) cosmological constant $(\sim{\Lambda}_{W})$, and the Hubble parameter $H(t)\equiv\dot{a}/a$. However, we take the coupling constants ${\alpha}_{3},\cdots,{\alpha}_{6}$ for higher-derivative (UV) terms to be arbitrary so that a (nearly) scale-invariant CMB spectrum with respect to the background universe, as well as power-counting renormalizability, can be obtained Shin:2017 ⁵⁵5For the cosmological perturbation around the spatially-flat FLRW background, the scale-invariant spectrum for the tensor modes depend only on the coupling ${\alpha}_{4}$, whereas the scalar mode depends on the combination $\widetilde{{\alpha}}_{4}\equiv{\alpha}_{4}+2{\alpha}_{5}/3+8{\alpha}_{6}/3$.. However, it is important to note that there are no contributions in the above Friedmann equations from the fifth and sixth-derivative UV terms in the Hořava gravity action (2) due to $R_{ij}={2k}g_{ij}/{R_{0}^{2}a^{2}(t)},~{}K_{ij}=H(t)g_{ij}$, but only from the fourth-derivative terms, which leads to $k^{2}/a^{4}$ terms ⁶⁶6If we include non-derivative higher-curvature terms, like $R^{3},R^{2}R_{ij}R^{ij}$, etc., we have $a^{-6}$ terms as well, which correspond to stiff matter Dutt:2009 ; Son:2010 . However, in this paper we do not consider those terms for simplicity. in the Friedmann equations (6) and (7). On the other hand, for the spatially-flat universe with $k=0$, all the contributions from the higher-derivative terms disappear and we recover the same background cosmology as in GR, which means a return to the LCDM model with the Hubble constant tension. In this sense, Hořava gravity is a “natural laboratory” for the tests of the spatially non-flat universe in standard cosmology.

Introducing dust matter (non-relativistic baryonic matter and (non-baryonic) cold dark matter with $p_{m}=0$) and radiation (ultra-relativistic matter with $p_{r}=\rho_{r}/3$), which satisfy the continuity equations $\dot{\rho}_{i}+3H(\dot{\rho}_{i}+p_{i})=0~{}(i=m,r)$, we conveniently define the canonical density parameters at the current epoch $a_{0}=1$ as ⁷⁷7We adopt the common convention ${\Omega}_{i}$ for the current values and ${\Omega}_{i}(a)$ for the fully time-dependent values.

where ${\beta}\equiv\kappa^{2}/2(3\lambda-1)$, ${\gamma}\equiv\kappa^{4}\mu^{2}{\Lambda}_{W}/8(3\lambda-1)^{2}$ are positive constant parameters. Then, we can write the (first) Friedmann equation (6) as

where we have introduced the (dynamical) dark-energy (DE) component as

which is defined as all the extra contributions to the first-three GR terms in (10) Park:2009 . Here, we note that this dynamical dark-energy component includes the dark radiation (DR) and dark curvature (DC) components as

which play the roles of the extra radiation and curvature terms, respectively, as well as the usual cosmological constant component $\Omega_{\Lambda}$, so that $\Omega_{\rm DE}(a)=\Omega_{\rm DR}(a)+\Omega_{\rm DC}(a)+\Omega_{\Lambda}$.

So far we have not assumed any specific information about the early universe. The whole physics of early universe in the context of Hořava gravity would be quite different from our known physics and needs to be revisited for a complete analysis. However, as in the standard cosmology, we may introduce the phenomenological parametrization so that all the unknown (early-universe) physics can be taken into account. For example, regarding the Big Bang Nucleosynthesis (BBN) at the decoupling epoch $a_{\rm dec}=(1+z_{\rm dec})^{-1}=(1091)^{-1}$, the early dark radiation can be expressed as the contribution from the hypothetical excess $\Delta N_{\rm eff}$ in the standard model prediction of the effective number of relativistic species $N_{\rm eff}=3.046$ as Stei:2012

where we have used the present radiation density parameter for the standard model particles with negligible masses (photon and three species of neutrinos), ${\Omega}_{r}=[1+\frac{7}{8}(\frac{4}{11})^{4/3}N_{\rm eff}]~{}{\Omega}_{{\gamma}}=\frac{7}{8}(\frac{4}{11})^{4/3}\cdot(0.13424)^{-1}{\Omega}_{{\gamma}}$ with the photon density ${\Omega}_{{\gamma}}=2.4730\times 10^{-5}h^{-2}$ for the present CMB temperature $T_{0}=2.7255K$ Planck:2018 and $h\equiv H_{0}/100{\rm~{}kms^{-1}Mpc^{-1}}$. If we use the relation (13) for the dark radiation formula (12) in Hořava gravity, we can express the cosmological constant component ${\Omega}_{{\Lambda}}$ as

so that the Friedmann equation (10) can be written by $\Delta N_{\rm eff}$, instead of ${\Omega}_{{\Lambda}}$ Dutt:2009 ; Nils:2018 . Here, it is important to note that $\Delta N_{\rm eff}$ needs not to be an integer but can be an arbitrary and positive (negative) real number for a positive (negative) ${\Omega}_{{\Lambda}}$, i.e., asymptotically de Sitter (Anti de Sitter) universe with a cosmological constant $\sim{\Lambda}_{W}$. Moreover, we note that the relation (14) implies an intriguing correlation between ${\Omega}_{k},{\Omega}_{r},\Delta N_{\rm eff}$, and ${\Omega}_{{\Lambda}}$, which are otherwise unrelated.

In this paper, we consider two models, A and B, depending on whether we implement (14) or not, to see the effectiveness of the BBN-like parametrization in terms of $\Delta N_{\rm eff}$. Then, for model A, the Friedmann equation (10) reduces to

which is one of our main equations for comparing with cosmological data, with the assumption of non-vanishing ${\Omega}_{k}$ and $\Delta N_{\rm eff}$ for the well-defined equation (15). Here we note that, considering ${\Omega}_{r}$ is a function of $H_{0}$ as given above, this model is described by five non-linearly coupled parameters $H_{0},{\Omega}_{m},{\Omega}_{k},{\Omega}_{\omega},\Delta N_{\rm eff}$, in contrast to the three linear, decoupled parameters $H_{0},{\Omega}_{m},{\Omega}_{\Lambda}$ for the standard, spatially-flat (background) cosmology, LCDM ⁸⁸8From the current-epoch constraint, i.e., $H(t_{0})\equiv H_{0}$ in (10), the number of independent parameters reduces to 4 and 2 for our two models based Hořava gravity and the standard models based on GR, respectively. In this paper, we conveniently take $H_{0},{\Omega}_{k},{\Omega}_{\omega},\Delta N_{\rm eff}$(or ${\Omega}_{\Lambda}$) and $H_{0},{\Omega}_{k},{\Omega}_{\Lambda}$ for the models A, B, and (k)LCDM, respectively, by considering ${\Omega}_{m}$, as well as ${\Omega}_{r}$, as the derived (dependent) parameter. However, we need to introduce the baryonic parameter ${\Omega}_{b}$ as an independent parameter when we analyze CMB and BAO data later so that only the (non-baryonic) dark matter sector in ${\Omega}_{m}$ is the derived parameter..

On the other hand, for model B, the Friedmann equation (10) is simply given by

without using the relation (14) from the BBN-inspired formula (13) for the dark radiation in Hořava gravity. This model has two (non-linearly) coupled parameters ${\Omega}_{k}$ and ${\Omega}_{\omega}$, in addition to those of standard flat cosmology, $H_{0},{\Omega}_{\Lambda}$, with assuming a non-vanishing ${\Omega}_{\Lambda}$ for the well-defined equation (16). We also consider the spatially flat $(k=0)$ LCDM as well as non-flat $(k\neq 0)$ kLCDM for comparison.

We probe our universe by using a Markov-Chain Monte Carlo (MCMC) method for the four models, A, B, LCDM, kLCDM and determine the independent parameters with the statistical inferences for the models. We use the Metropolis-Hastings algorithm ⁹⁹9We use the Metropolis-Hastings algorithm at the background level since there is no known analysis of the cosmological perturbations for a non-flat universe, contrary to the flat universe Shin:2017 , due to computational complications in Hořava gravity. However, the algorithm has the limitation in our case, though not a general feature, that does not show the proper $\chi^{2}$ values for the separate data sets, while it shows the convergent results for the whole data sets. Robe:2015 for the posterior parameter distributions, the statistical method in Dunk:2004 for convergence tests of the MCMC chains, and the public code GetDist Lewi:2019 for the visualization of the results. The cosmological data sets we consider (for the details, see Appendix A) are CMB (Planck 2018 Zhai:2018 ), BAO (SDSS-BOSS BOSS:2016 , SDSS-IV Ata:2017 , Lyman-$\alpha$ forest deSa:2019 , and WiggleZ Blak:2012 ), SNe Ia (Pantheon Scol:2017 ), GRBs (Mayflower Liu:2014 ), Lensed Quasars (H0liCOW Wong:2019 ), and Cosmic Chronometers (CC) More:2015 . Here, it is important that we need to include CMB lensing data in combination with BAO or SNe Ia in order to test the curvature of the universe Efst:2020 ; Gonz:2021 .

Our main results are shown in Table 1 and their essentials are plotted in Fig. 1 (for the full plots, see Appendix B). Some noticeable results are as follows:

1. For our two main models A and B, constraints with all the data sets, including BAO, show some reduced tensions in the Hubble constant $H_{0}$, from the direct local measurements using Cepheids Ries:2018 (see Fig. 2 for a comparison with other previous measurements Park:2018 ; Nune:2020 ; Beni:2020 ; Vagnozzi:2020 ; Vagnozzi:2020dfn ). In particular, for model A, the tension is reduced to within $1{\sigma}$ even if BAO data is included.

On the other hand, for ${\Omega}_{m}$ and ${\Omega}_{\Lambda}$, there are only some slight differences from the standard LCDM or kLCDM, and the conventional constraint of ${\Omega}_{m}\approx 0.3$ and ${\Omega}_{\Lambda}\approx 0.7$ Ries:1998 is still robust even in Hořava cosmology. This result is important because it may indicate a resolution of the discordance problem in the LCDM paradigm, which constrains ${\Omega}_{m}\sim 0.5,{\Omega}_{\Lambda}\sim 0.5$ when a closed universe is considered as preferred by Planck DiVa:2019 ; Hand:2019 , but incompatible with other local measurements. However, since our result may not be its complete confirmation because we do not have the results for separate data sets but only for their combinations, as noted in the footnote No. 7.

The constraint of $\Omega_{r}$, for model A, is not overlapping within $2\sigma$ CL of the Planck 2018 data, $({9.179^{+0.441}_{-0.431}})\times 10^{-5}$ (CMB+lensing+BAO) CMB Wiki . However, it is not surprising that $\Omega_{r}$, which is a derived quantity from the standard formula below Eq. (13) and reduced by a factor $H_{0}^{-2}$ from the increased constraint on $H_{0}$, which is in about $4\sigma$ tension with the Planck 2018 data $H_{0}=67.27\pm 0.60~{}{\rm kms^{-1}Mpc^{-1}}$ Planck:2018 , has a similar tension with the Planck data.

The constraints of ${\Omega}_{\omega}$ are rather different for models A and B. ${\Omega}_{\omega}$ is a newly introduced parameter in our models and it has no a priori known constraint. But our results show that ${\Omega}_{\omega}<0$ is preferred at $1.6{\sigma}$ for model A so that ${\Omega}_{k}{\Omega}_{\omega}>0$, whilst poorly constrained for model B.

2. ${\Omega}_{k}$ is peculiar in two aspects as follows. First, the result for model A shows a higher precision, i.e., narrower MCMC contour (see ${\Omega}_{k}$ vs. ${\Omega}_{m}$ in Fig. 1, for example) and, as a result, a closed universe, i.e., ${\Omega}_{k}<0$, is more strongly preferred than in kLCDM (cf. DiVa:2021 ). (See Fig. 3 for the comparisons with other results based on a minimal extension to LCDM) This is peculiar because the lower precision (i.e., wider distribution) is normally expected with the addition of more parameters, as can be seen in ${\Omega}_{b}h^{2},{\Omega}_{r}$, and $H_{0}$ for A and B, in contrast to kLCDM. Second, the correlation between $H_{0}$ and ${\Omega}_{k}$ is quite different for A or B and kLCDM: As $-{\Omega}_{k}$ increases, i.e., the universe tends more towards a closed universe, $H_{0}$ increases for A and B, i.e., the two parameters are anti-correlated, whereas the situation is the opposite for kLCDM. This property for $H_{0}$ resolves the problem of Hubble constant tension in kLCDM DiVa:2019 ; DiVa:2020 ¹⁰¹⁰10It is interesting that a similar effect of curvature has been noted earlier in Ries:2019 [see Fig. 4] though it would be hardly justified in the data sets with the CMB, based on GR Vagnozzi:2020 ; Vagnozzi:2020dfn (see Fig. 5 for the $H_{0}$ vs. $\Omega_{k}$ tension in the measurements of Figs. 2 and 3). These peculiar properties seem to be due to the non-linear coupling of ${\Omega}_{k}$ in the Friedmann equations (15) and (16).

3. In our analysis, we use the CMB distance priors, or shift parameters, instead of the full CMB data (for the details, see Appendix A). It is possible that this choice may lead to a lower statistical weight generally, though it has been widely used and is a convenient substitute of the full CMB data. However, the comparison of our analysis for the kLCDM or LCDM case and the full CMB analysis in Fig. 2 when combining other data sets, including BAO, shows that it has just a small effect as can be seen. For example, CMB+BAO+SN eBOS:2020 , CMB+BAO Full Vagnozzi:2020 , CMB+lensing+BAO Planck:2018 (with or without $\Omega_{k}$) that use the full CMB data, show a few percent effects ($<2\%$), in comparison with our LCDM, kLCDM results with the CMB priors ¹¹¹¹11There is a more systematic analysis on the test of the CMB priors which supports the use of the CMB priors for the background analysis of dark-energy dynamics Zhai:2019 . It would be an interesting problem to confirm the accuracy of the CMB priors for our Hořava cosmology model also in a more systematic way..
4. The constraint of $\Delta N_{\rm eff}$, for model A (considering its central value), is consistent with SPT-3G 2018 data (alone), $0.65\pm 0.70$ (CMB)¹²¹²12In this case, the constraint of $H_{0}$ is given by $73.5\pm 5.2~{}{\rm kms^{-1}Mpc^{-1}}$. SPT:2021 by $0.31\sigma$, while in a tension with Planck 2018 data, for example, $2.7-2.8\sigma$ of $-0.06^{+0.34}_{-0.33}$ (CMB+lensing+BAO) Planck:2018 . Since, in our MCMC analysis, we have not included SPT-3G 2018 data, which gives better constraints at the higher angular multipoles $l>1000$, the agreement of our constraint of $\Delta N_{\rm eff}$ might give an independent support of our model A. However, due to the lack of data with both $\Delta N_{\rm eff}$ and $\Omega_{k}$, which would be quite (anti-)correlated (see Fig. 4), as the varying parameters in SPT and Planck, the exact quantification of (in)compatibility is not available yet. But we can generally expect some tensions with Planck 2018 in $\Delta N_{\rm eff}$ as in the case of $H_{0}$, due to the correlation of $\Delta N_{\rm eff}$ and $H_{0}$ in Planck 2018 Planck:2018 ; SPT:2021 , as well as in SPT-3G 2018 SPT:2021 and our case (Fig. 4).

It is known that primordial gravitational waves may add to the effective relativistic degrees of freedom Mang:2001 . For example, in a recently proposed scenario linking $\Delta N_{\rm eff}$ with the gravitational waves from primordial black holes Flor:2020 , it is possible to have $\Delta N_{\rm eff(GW)}\sim 0.1-0.2$ which may increase the standard value $N_{\rm eff}=3.046$ and its associated radiation parameter ${\Omega}_{r}$, slightly. In some different scenarios, it has been also found that the primordial gravitational waves do not solve but alleviate the $H_{0}$-tension problem Grae:2018 . However, in all those cases, a spatially-flat universe is assumed implicitly and the effect of a non-flat universe seems to still be an open problem. Moreover, due to some inaccuracies of the (as adopted by us) CMB priors for the primordial power spectrum in a non-flat universe Zhai:2019 , in this paper, we will not quantify its effect in the $H_{0}$-tension problem further which needs another full data analysis as well.

5. From Bayes theorem, the Bayesian evidence $E({\cal D}|{\cal M})$ for a model ${\cal M}$ with the total data sets ${\cal D}$ is given by the integration over the model parameters ${\bf\theta}$

where ${\cal L}({\cal D}|{\bf\theta},{\cal M})$ is the likelihood ${\cal L}({\cal D}|{\bf\theta},{\cal M})\equiv exp[-\chi^{2}({\cal D}|{\bf\theta},{\cal M})/2]$ in which the total $\chi^{2}$ is obtained by summing $\chi^{2}({\cal D}_{i}|{\bf\theta},{\cal M})$ over all the data sets ${\cal D}_{i}$. $\pi({\bf\theta}|{\cal M})$ is the prior probability, which we have assumed to be flat, i.e., no prior information on $\bf\theta$, in order to be as agnostic as possible: our only priors ¹³¹³13In general, if the prior range is too small, the parameter chain can be seen to ‘hit a wall’ at one end of the prior (which serves as a hard cutoff), but we have not observed this behavior in our analysis. For the 2D contours in Fig. 4, in particular for $\Omega_{k}$ and $\Omega_{w}$, we zoom out enough to show the details of model A as clearly as possible but at least without cutting out $2\sigma$ CL, in compatible with $2\sigma$ contours in Fig. 1. However, for the 1D contours of $1\sigma$ CL in Fig. 4, we have the desired vanishing tails of the posterior distributions at the boundary of the priors, and all of our results have passed the convergence tests in Dunk:2004 . are $0<\Omega_{b}<\Omega_{m},~{}0<H_{0}<100$ for LCDM, kLCDM; in addition, $\Omega_{k}\neq 0,~{}\Delta N_{\rm eff}\neq 0$ for model A and $\Omega_{\Lambda}\neq 0$ for model B, in order to avoid the singularity of the corresponding Friedmann equations (15) and (16), respectively ¹⁴¹⁴14The infinitesima widths of the excluded parameter regions, which are basically discrete, around the singularities are not fixed but randomely chosen in the MCMC analysis. But, the smooth contours around the singuraities in Fig. 4 show that our chosen priors work well.. The differences of $\chi^{2}_{\rm min}$ for the models A and B, with respect to LCDM, $\Delta\chi^{2}_{\rm min}={-14.3,-7.8}$ indicate the notable improvements of fitting to the given data sets, in contrast to the smaller difference $\Delta\chi^{2}_{\rm min}={-2.9}$ for kLCDM.

For a more quantitative comparison of the models, we consider the Bayes factor

which quantifies the preference for model ${\cal M}_{i}$ against model ${\cal M}_{j}$, using the (revised) Jeffrey’s scale Jeff:1939 ; Kass:1995 ; Ness:2012 : weak ($0\leq{\rm ln}B_{ij}<1.1$), definite ($1.1\leq{\rm ln}B_{ij}<3$), strong ($3\leq{\rm ln}B_{ij}<5$), very strong ($5\leq{\rm ln}B_{ij}$). Table 1 gives the Bayes factors for model ${\bf A}$ and ${\bf B}$ with respect to LCDM, $\ln{B_{ij}}={+6.2,+3.1}$, i.e., “very strong” and “strong” evidences, respectively, against the flat LCDM, in contrast to the “definite”¹⁵¹⁵15This is in contrast to the “strong” evidence for kLCDM with $\ln{B_{ij}}=+3.3$ for CMB data alone DiVa:2019 . evidence for kLCDM with $\ln{B_{ij}}={+1.1}$ ¹⁶¹⁶16The detailed decisiveness of the results can depend on the adopted scales. For example, in Trotta’s revisited scale Trot:2008 , our results show “strong”, “moderate”, and“weak” eveidence, respectively. .

6. The theory parameters can be written in terms of cosmological parameters as

where $M_{P}$ and $L_{P}$ are the Planck mass and length, respectively, with $M_{P}/L_{P}=c^{2}/8\pi G={\alpha}^{-1}$. Then we obtain their best-fit values¹⁷¹⁷17It is interesting to note that this is the case of ${\omega}<0,~{}{\omega}<2{\Lambda}_{W}$ in which the observer region is located between the inner and outer (black hole) horizons with the cutting-edge (surface-like) singularity of the space-time (or, the end of world) inside the inner horizon Argu:2015 ., $(\bar{{\omega}},~{}\bar{{\Lambda}}_{W},~{}\bar{\mu})=(-275.94,~{}270.99,~{}0.0062),~{}(173.51,~{}390.24,~{}0.0043)$, where $\bar{{\omega}}\equiv{\omega}R_{0},~{}\bar{{\Lambda}}_{W}\equiv{\Lambda}_{W}R_{0},~{}\bar{\mu}\equiv\left(\mu^{2}L_{P}/H_{0}R_{0}M_{P}\right)^{1/2}$, for the models A and B, respectively¹⁸¹⁸18These correspond to the CPL parameters, $({\omega}_{0},~{}{\omega}_{a})=(-1.005,~{}-0.010),~{}(-0.998,~{}{0.004})$ for the expansion of ${\omega}_{\rm DE}={\omega}_{0}+{\omega}_{a}(1-a)+{\omega}_{b}(1-a)^{2}+\cdots$, near the current epoch $a=1$.. There are large errors (see Table 2 and Fig. 6) due to the non-linearity of the relation (19) but their best-fit values are distributed near what can be obtained by plugging in the obtained best-fit values of observational parameters in Table 1, $(\bar{{\omega}},~{}\bar{{\Lambda}}_{W},~{}\bar{\mu})=(-299.20,~{}281.20,~{}0.00596),~{}(167.50,~{}347.50,~{}0.00479)$. These would be the first full (cosmological) determination of the theory parameters (for the earlier determinations ¹⁹¹⁹19There are sign errors for the results of $\bar{{\omega}}$ and $k$ in Park:2009 ., see Park:2009 ).

In conclusion, we have tested the spatially non-flat universe in standard cosmology within Hořava  gravity. We have obtained the “very strong” and “strong” Bayesian evidences against flat LCDM, for our two models A and B, respectively. Moreover, the MCMC analysis shows (a) the preference of a closed universe, (b) a reduction of the Hubble constant tension, and (c) a positive result on the discordance problem, even if BAO data is included. It is remarkable that just the use of (14) for model A, which gives a novel relation between ${\Omega}_{\Lambda},{\Omega}_{k},{\Omega}_{r}$ via $\Delta N_{\rm eff}$ which are otherwise unrelated, produces the difference in the results of the two models.

The reduced Hubble tension may be related to a natural inclusion of the (early) dark radiation within the dynamical dark energy and its associated contribution to $\Delta N_{\rm eff}=0.87$ for model A, and the peculiar contribution of the curvature ${\Omega}_{k}<0$, as have been previously considered on phenomenological ground Ries:2018 ; Knox:2019 . However, the resolution does not seem to be quite complete yet Ries:2020 ; this might be due to the fact that our observational data may not be completely model independent and are based on known physics. For example, as noted earlier, the BBN-like parametrization (13), which has been used for model A, is based on standard particle physics model with the phenomenological parameter $\Delta N_{\rm eff}$. Our MCMC analysis implies that the phenomenological approach of model A is a good approximation of the BBN, constraining ${\Omega}_{k}$ more precisely and reducing the Hubble constant tension, i.e., increasing $H_{0}$, with the almost doubled Bayesian evidence compared to model B. It would be a challenging problem to thoroughly consider the effect of anisotropic scaling beyond the standard model BBN and revisit the Hubble tension problem to see whether a complete resolution can be found or not. Finally, the analysis of the cosmological perturbations for the non-flat universe and with the full use of Boltzmann solvers such as CAMB/Class Lewi:1999 ; Lesg:2011 would also be an important arena for studying standard cosmology, like the cosmic-shear or ${\sigma}_{8}$ tension DiVa:2019 .