Joint Constraints on the Hubble Constant, Spatial Curvature, and Sound Horizon from the Late-time Universe with Cosmography

In this paper, using the latest Pantheon+ sample of Type Ia supernovae, baryon acoustic oscillation measurements, and observational Hubble data, we carry out a joint constraint on the Hubble constant H 0, the spatial curvature ΩK, and the sound horizon at the end of the drag epoch r d. To be model-independent, four cosmography models—i.e., the Taylor series in terms of redshift y 1 = z/(1 + z), y2=arctan(z) , y3=ln(1+z) , and the Padé approximants—are used without the assumption of a flat Universe. The results show that H 0 is anticorrelated with ΩK and r d, indicating that smaller ΩK or r d would be helpful in alleviating the Hubble tension. The values of H 0 and r d are consistent with the estimate derived from the Planck cosmic microwave background data based on the flat ΛCDM model, but H 0 is in 2.3 ∼ 3.0σ tension with that obtained by Riess et al. in all these cosmographic approaches. Meanwhile, a flat Universe is preferred by the present observations under all approximations except the third order of y 1 and y 2 of the Taylor series. Furthermore, according to the values of the Bayesian evidence, we found that the flat ΛCDM remains to be the most favored model by the joint data sets, and the Padé approximant of order (2,2), the third order of y 3 and y 1 are the top three cosmographic expansions that fit the data sets best, while the Taylor series in terms of y 2 are essentially ruled out.


INTRODUCTION
In modern cosmology, the Hubble tension has become one of the most prominent issues.The Hubble constant H 0 , as a fundamental cosmological parameter, represents the expansion rate of the present Universe, which can be obtained by both global and local measurements, i.e., based on the standard cosmological model(ΛCDM), Planck Collaboration inferred H 0 = 67.4± 0.5 km s −1 Mpc −1 (hereafter P18) (Planck Collaboration et al. 2020) from Planck satellite measurements of Cosmic Microwave Background (CMB) temperature and polarization anisotropies, while, Supernovae and H 0 for the Equation of State of dark energy (SH0ES) collaboration found H 0 = 73.04 ± 1.04 km s −1 Mpc −1 (hereafter R22) (Riess et al. 2022) via local measurements from type Ia supernovae (SNe Ia) calibrated by the distance ladder without any cosmological model.The 5σ tension between these two independent H 0 estimates could be pointing towards new physics beyond the standard model or residual systematics (Freedman 2017;Feeney et al. 2018;Di Valentino et al. 2018;Handley 2021).To alleviate the discrepancy, fundamental physics beyond ΛCDM are investigated, such as time-dependent dark energy equation of state (Huang & Wang 2016;Di Valentino et al. 2016, 2021;Zhao et al. 2017;Miao & Huang 2018;Yang et al. 2019;Poulin et al. 2019;Vagnozzi et al. 2021;Colgáin et al. 2021;Yang et al. 2023), modified gravity (Capozziello et al. 2003;Nunes 2018;Farrugia et al. 2021;Koussour et al. 2022;Sultana et al. 2022), and additional relativistic particles, see Kumar & Nunes (2016); Xu & Huang (2018); Carneiro et al. (2019); Pandey et al. (2020);D'Eramo et al. (2022).Alternately, new approaches to determine the Hubble constant from local direct measurements are proposed.For example, through time-delay cosmography, the H 0 Lenses in COSMOGRAIL's Wellspring (H0LiCOW) Collaboration (Suyu et al. 2017) measured the Hubble constant from strong gravitational lens systems with time delays between the multiple images.Using detected gravitational waves (GW) as a 'standard siren', binary neutron-star (BNS) system detection GW170817 and subsequent observations in the electromagnetic (EM) domain provide another independent method of measuring the Hubble constant (Abbott et al. 2017).Adopting revised measurement, Freedman et al. (2020) found H 0 = 69.9± 0.8(±1.1% stat)±1.7(±2.4% sys)km s −1 Mpc −1 , which provided one of the most accurate means of measuring the distances to nearby galaxies by the red giant branch method.
However, if the systematic uncertainties are not the main drivers of H 0 tension, then the discrepancy coming from the concrete cosmological model assumption can not be ignored.In order to extricate from dependency on a cosmological model and study the expansion of the Universe directly from the observations, various model-independent techniques are used, such as cosmography, the Bézier parametric curve, the Parameterization based on cosmic Age (PAge), the Gaussian process (GP) method and so on (Wojtak & Agnello 2019;Capozziello et al. 2019;Yang et al. 2020;Zhang & Huang 2021;Cai et al. 2022a,b;Hu & Wang 2022;Jalilvand & Mehrabi 2022;Liu et al. 2022).For examples, Zhang & Huang (2021) reconstructed H(z) in cubic expansion and polynomial expansion respectively, and constrained H 0 and r d with the joint data of SNe Ia, Baryon Acoustic Oscillation (BAO) measurements, Observational Hubble Data (OHD), and GW data, and found the H 0 value is in 2.4-2.6σtension with SH0ES 2019 (Riess et al. 2019).Cai et al. (2022a,b) applied the PAge approximation to consistently use the OHD and the late-time matter perturbation growth data at high redshifts, and found the Hubble tension can't be solved by introducing the new physics at the late time beyond the ΛCDM model.Hu & Wang (2022) investigated the redshift-evolution of H 0 with 36 Hubble parameter H(z) data based on the GP method, and found there was a late-time transition of H 0 which effectively alleviated the Hubble crisis by 70%.It should be pointed out that these works are all taken under the assumption of a spatially flat Universe.However, there is still one observational probe that is in agreement with a negative curvature, sparking the debate about the flatness of the Universe.From the combination of CMB temperature and polarization power spectra, the constraint on curvature suggests a closed Universe at more than three standard deviations (Ω K = −0.044+0.018  −0.015 ) (Planck Collaboration et al. 2020).Whereas once the Planck CMB data are combined with the BAO measurements, a flat Universe is preferred with a high precision Ω K = 0.001 ± 0.002 (Planck Collaboration et al. 2020), indicating there is a discrepancy between CMB and BAO data for the constraint on Ω K , namely "curvature tension" (Handley 2021).As allowing for spatial curvature may significantly affect constraints on cosmological parameters (Dossett & Ishak 2012), the curvature parameter also needs to be considered when exploring Hubble tension.For instance, it is found that a negative curvature preferred by P18 will exacerbate H 0 tension (Di Valentino et al. 2020, 2021), while Wang et al. (2021) and Cao et al. (2022) found that a spatially flat Universe is favored using the mock GW data in combination with the Hubble parameter or strong gravitational lensing time delay (SGLTD) data respectively.Therefore, due to the significant influence of spatial curvature on the H 0 measurement and the inconsistency of the spatial curvature measurement, it is very necessary for us to measure H 0 with various observations under the assumption of a non-flat Universe.
Additionally, since that the value of H 0 obtained from the Planck CMB data is directly tied to the sound horizon at last scattering, which is closely related to the sound horizon r d at the baryon decoupling, a number of amendments to the ΛCDM model have been proposed, aiming to release the Hubble tension by reducing r d and increasing H 0 (Karwal & Kamionkowski 2016;Buen-Abad et al. 2018;Anchordoqui & Bergliaffa 2019;Niedermann & Sloth 2020;Berghaus & Karwal 2020;Archidiacono et al. 2020).However, the most recent works demonstrated that any model which only reduces r d can never fully resolve the Hubble tension (Pogosian et al. 2020;Jedamzik et al. 2021).In addition, in order to break the measured degeneracy between r d and H 0 , many works have therefore attempted to measure the H 0 and r d by combining the BAO measurements with other late-time observations of Universe.For instance, Wojtak & Agnello (2019) combined BAO with the SGLTD, Joint Light-Curve Analysis (JLA), and the reconstructed H(z) data via the polynomial expansion, and obtained H 0 = 72.3± 6.9 km s −1 Mpc −1 , r d = 139.2± 13.3 Mpc.Then Pogosian et al. (2020) found r d = 143.7 ± 2.7 Mpc and H 0 = 69.6 ± 1.8 km s −1 Mpc −1 by using the latest BAO data along with a Ω m h 2 prior based on the Planck best-fit ΛCDM model, and similar values were obtained when they combined BAO with the Pantheon supernovae, the Dark Energy Survey Year 1 galaxy weak lensing (Abbott et al. 2018), Planck/SPTPo1 CMB lensing, and OHD.In addition, Cai et al. (2022b) combined the SNe Ia, BAO and OHD, and obtained H 0 = 68.958+1.779  −1.826 km s −1 Mpc −1 and r d = 146.466+3.448  −3.302 Mpc by using the PAge approximation.Although the Hubble tension can be alleviated in these works, it cannot be completely resolved.To reflect the reality of these works taken under the assumption of a spatially flat Universe, it is well worth constraining H 0 and r d with leaving Ω K free by using the latest low-redshift observational data.
In this work, we therefore plan to carry out a joint constraint on H 0 , Ω K and r d with the latest late-time observations of Universe including SNe Ia, BAO, and OHD.Specially, the newest Pantheon+ sample of SNe Ia (Brout et al. 2022;Scolnic et al. 2022) ranging in redshift from z = 0.001 to 2.26, is used.This sample is made of 1701 light curves of 1550 spectroscopically confirmed SNe Ia and significantly enlarges the origin Pantheon sample size from the addition of multiple cross-calibrated photometric systems of SNe Ia.For the BAO measurements, fourteen latest measurements of D V (z)/r d , D M (z)/r d and D H (z)/r d summarized in Ref. (Alam et al. 2021), covering the redshift range 0.15 ≤ z ≤ 2.33, are used.These data are obtained from final observations of clustering using galaxies, quasars, and Lyα forests from the completed SDSS lineage of experiments in a large-scale structure, composing of data from SDSS, SDSS-II, BOSS and eBOSS.Furthermore, to be model-independent, the well-consolidated approach named cosmography (Chiba & Nakamura 1998;Caldwell & Kamionkowski 2004;Visser 2004Visser , 2005Visser , 2015;;Capozziello et al. 2013Capozziello et al. , 2019;;Dunsby & Luongo 2016;Zhang et al. 2017;Yin & Wei 2019;Li et al. 2020) is used, which has attracted lots of attention in the study of the expansion of Universe.The idea of cosmography is to expand the cosmological distances or the Hubble parameter into a Taylor series of redshift z, which performs well at low redshifts but encounters the convergence problems in the high-redshift domain (Cattoën & Visser 2007).To overcome the convergence issues, several improved approaches have been proposed, one of which relies on the use of auxiliary variables (Cattoën & Visser 2007;Aviles et al. 2012;Capozziello et al. 2020), and the other expands observables in terms of rational approximations (Gruber & Luongo 2014;Wei et al. 2014;Shafieloo 2012;Capozziello et al. 2018).Recently, Li et al. (2020) adopted two cosmographic methods, the Taylor series in terms of y = z/(1 + z) and Padé polynomials, to investigate the spatial curvature parameter with the dataset including the Pantheon sample of SNe Ia, BAO, and OHD data, and found the H 0 tension problem can be slightly relaxed by introducing the spatial curvature parameter.In order to get more robust results, we explore the Hubble tension through the extended cosmographic techniques with Taylor series in terms of y 1 = z/(1 + z), y 2 = arctan(z), y 3 = ln(1 + z) and Padé approximations using the latest data, while leaving Ω K and r d free at the same time.Meanwhile, a Bayesian approach based on Bayesian evidence is applied to test which cosmographic method is most favored by the observational data.
This paper is organized as follows: In Section 2, we introduce the cosmographic approaches, dataset, and methodology used in this work.In Section 3, our constraint results and analysis are presented.Finally, the conclusions are drawn in Section 4.

METHODOLOGY AND DATA
The cosmographic approach is an artful combination of kinematic parameters via the Taylor series with the assumption of large-scale homogeneity and isotropy, which can be retained in the Friedmann-Robertson-Walker (FRW) metric, where c is the speed of light, K is the constant curvature of the three-space of the FRW metric, and a(t) is the scale factor with cosmic time t.Considering a photon traveling towards us along a radial path (ds = 0), we can get Integrating the equation, the comoving distance can be obtained by where t 0 and a 0 are the time and scale factor when the photon was observed from us, and t e is the time when it was emitted from the source at r = r e .Substituting 1 + z = a 0 /a, then the comoving distance becomes where E(z) ≡ H(z)/H 0 , and the transverse comoving distance (Hogg 1999) can be expressed by where Ω K = −Kc 2 /(a 0 H 0 ) 2 is the present value of spatial curvature parameter.Now define the cosmographic parameters, i.e., Hubble parameter H, deceleration q, jerk j, snap s and lerk l parameters, as follows, We can expand D M (z) up to the fifth order, where with the subscript "0" denoting the values at the present time.Obviously, the H(z) can be expanded directly with the cosmographic parameters, but the expansion is independent of Ω K .In order to expand H(z) with Ω K included and improve the usefulness of OHD, we obtain the Hubble function in terms of Ω K and D M (z), However, Cattoën & Visser (2007) showed that z-based expansions must break down for z > 1.In order to avoid the convergent problem at high redshifts, improved cosmographic techniques have been proposed, namely the auxiliary y-variables and Padé approximations.A fairly well-known y variable is given by which is proposed by Cattoën & Visser (2007) and well performed all the way back to the big bang with a nice finite range [0, 1) for z ∈ [0, ∞).Moreover, in the work Aviles et al. (2012), another parametrization is adopted, which behaves smoothly with the arctangent function.It has been demonstrated that y 2 can give welldefined limits [0, π 2 ) in the range of z ∈ [0, ∞) (Aviles et al. 2012).Besides, the third parametrization is introduced by Semiz & Kazım C ¸amlıbel (2015).As a natural logarithmic function, y 3 tends to infinity as z → ∞, but it increases slowly with the redshift growing.
In addition to the conventional methodology of cosmography applying Taylor expansions of observables, Gruber & Luongo (2014) employed Padé approximants which have the superior convergence properties.For a given function f (z), the Padé approximant of order (m, n) is given by where m and n are non-negative integers, and a i and b i are constant that satisfy the conditions P mn (0) = f (0), P ′ mn (0) = f ′ (0),..., P m+n mn = f m+n (0).According to the plots of the transverse comoving distance D M (z) and Hubble parameter H(z) for the third-, fourth-and fifth-order Padé approximants in Figure 1 of the work Li et al. (2020), the P 12 , P 22 , P 13 , P 32 expansions can give good approximations to the ΛCDM model over the redshift interval z < 2.4 while other expansions will diverge from ΛCDM model outside a low redshift region.
In this work, we adopt these three parameterizations of Taylor series as well as the Padé approximants to reconstruct the expansion evolution of the Universe and constrain H 0 , Ω K and r d .In order to achieve high accurate performance without introducing too many model parameters, the third-, fourth-and fifth-order approximants, namely y (j) i (i=1,2,3; j=3,4,5) for the y series and P 12 , P 22 , P 13 , P 32 for the Padé series, are used.We also utilize the Bayesian evidence method to study which cosmographic approach performs the best.More details about the explicit cosmography expressions and Bayesian evidence are reported in Appendix A and B.
As a guideline, the Hubble parameter in the ΛCDM+Ω K model is introduced, with corresponding cosmographic parameters being where Ω m represents the current matter density.In particular, when Ω K = 0 corresponds to the flat ΛCDM model, the terms containing Ω K in the above equations will disappear.
In this paper, we use the newest Pantheon+ SNe Ia sample (Brout et al. 2022;Scolnic et al. 2022), the 14 latest BAO measurements (Alam et al. 2021), and the 32 OHD (Wu et al. 2023) obtained from cosmic chronometer method.More details about the data and fitting methods are reported in Appendix C.

CONSTRAINT RESULTS
We use CosmoMC (Lewis & Bridle 2002) together with a nested sampling plug-in, namely PolyChord (Handley et al. 2015a,b) which enables the computation of the Bayesian evidence, to study four cosmography models.With imposing uniform priors on the free parameters as listed in Table 1, we obtain the constraint results by numerical calculations, and summarize the mean and 68.3% confidence limits of cosmological parameters in Table 2. And, the likelihood distributions of y 1 , y 2 , y 3 and Padé approximations from the third order to the fifth order are shown in Figure 1− 3, respectively.Furthermore, the logarithmic values of the Bayes factor lnB ij for the ΛCDM+Ω K and the four cosmography models are given in Table 3 with the flat ΛCDM model as the reference model.First, we focus on the Hubble constant H 0 and spatial curvature parameter Ω K .From Table 2 we can see that the values of H 0 obtained from the four cosmographic approaches are consistent with the P18 value at 1σ confidence level (CL), and lower than the R22 value by 2.3∼3.0σ.Meanwhile, they are in good agreement with that obtained from the flat ΛCDM and ΛCDM+Ω K models.Additionally, it is worth noting that under the ΛCDM+Ω K model and all cosmography models, a flat Universe is preferred by the present observations except for the y (3) 1 and y (3) 2 expansions in which a closed Universe is preferred at more than 2σ confidence level (CL).Meanwhile, from the H 0 − Ω K plane in Figure 1− 3, we can find that there is a visible anti-correlation between H 0 and Ω K , indicating that a larger H 0 can be obtained from a smaller Ω K , and in all the y approximations, the lower the expansion order, the smaller the constraint Ω K .
Second, compared with r d = 147.05± 0.3 Mpc reported by Planck 2018, our constraints on r d using all cosmographic approaches are in agreement with Planck 2018 within 1σ CL.Since the values of r d are essentially the same including the mean value and the uncertainty in four cosmography models, from Table 2 and Figure 1− 3, it can be concluded that the sound horizon r d is not sensitive to the parameterization and truncation of the expansions.In addition, from the H 0 − r d contours, we can see that there is an obvious anti-correlation between r d and H 0 .Therefore, if the expansion evolution of the Universe is determined, a smaller r d would result in a larger H 0 , which can also be seen from Equation (C19).Third, we analyze the constraints on the cosmological parameters q 0 , j 0 , s 0 and l 0 .From Table 2, it can be seen that the observations could provide a good constraint on the lower-order expansions, within which the values of q 0 and j 0 are basically consistent with that derived with the flat ΛCDM model, except the (1,2) order of Padé approximant.However, due to the complexity of high expansion orders and the insufficiency of high-redshift data, there are large uncertainties in the constraints of s 0 and l 0 by all the cosmography techniques.This suggests that new probes into the evolution of the Universe are needed to break the degeneracies of these parameters.Meanwhile, one can see from Table 2 that the uncertainties on the parameters of H 0 , Ω K , and r d remain roughly the same as one increases the expansion order, which is similar to the results obtained by Li et al. (2020).This is mainly because the degeneracy between the parameters of H 0 , Ω K , and r d , and the cosmographic parameters decreases with the expansion order increases as shown in Figures 1−3, and thus the constraint ability of observational data on the parameters of H 0 , Ω K , and r d is not substantially weakened with the increase of the number of model parameters.
Fourth, we apply the Bayesian evidences with the flat ΛCDM as the reference model to determine the most favored cosmography model for the joint datasets.From Table 3, one can see that all the lnB ij values in this table are negative, indicating that the flat ΛCDM is the most preferred model by the observations, and the ΛCDM+Ω K model is also more favored than the four cosmography models.Furthermore, according to the Jeffreys scale list in Table 4, the |lnB ij | values of five models, i.e. all expansions of y 1 , the third order of y 3 and the (2,2) order of Padé approximant, are smaller than 3, indicating that those models could fit the observations well.Especially, y 1 , y (3) 3 and P 22 are the top three models that are most supported by the observations.Meanwhile, we find that all the y 2 expansions have the largest values of |lnB ij | greater than 5, which means these kinds of expansions are strongly disfavored by the observational data, and can therefore be ruled out.In addition, we derive the evolutions of H(z) within the y redshift (z < 1 ∼ 1.5).And although the evolutions of H(z) from the cosmography models have some deviations from the flat ΛCDM model at the higher redshift, it is still in agreement with the flat ΛCDM model within 1σ CL.
Lastly, we compare our results with the previous works on the constraints of H 0 , Ω K and r d by model-independent methods.Adopting Taylor series in terms of y 1 and the Padé approximants along with the SNe Ia+OHD+BAO data, Li et al. (2020) showed that the tension level of H 0 has less than 2σ significance between different approximations and the local distance ladder determination while our constraint results show a tension of more than 2σ using the latest data.Furthermore, their constraint results on Ω K preferred a closed Universe under all the approximations, whereas we find the newest datasets prefer a flat Universe with a higher precision.Therefore, the latest BAO measurements and Pantheon+ sample data can improve the accuracy of constraints on H 0 and Ω K .Moreover, through Bayesian evidence, Li et al. (2020) found that the y (3) 1 has weak evidence against ΛCDM+Ω K model and P 22 has positive evidence, while our constraint results show that y (3) 3 behaves the best in all series of y-variables with the updated BAO data and Pantheon+ SNe Ia sample.In addition, assuming a spatially flat Universe, Zhang & Huang (2021) reconstructed H(z) in cubic expansion and polynomial expansion respectively, and constrained Hubble constant and sound horizon with the joint data of SNe Ia, BAO measurements, OHD, and GW data.And they found H 0 and r d were consistent with the estimate derived from Planck 2018 data, which is similar to our constraint results.Likewise, Cao et al. (2022) combined 55 simulated SGLTD data and 1000 simulated GW data, and obtained high-precision constraint results H 0 = 73.65 ± 0.33 km s −1 Mpc −1 and Ω K = 0.008 ± 0.038.Although those results showed a strong constraint on H 0 comparing to our constraint results, those works use simulated data, whereas we use actual observational data.

CONCLUSION AND DISCUSSIONS
The Hubble tension problem can be relaxed to some extent by introducing Ω K or any cosmological model that assumes a small r d in the early Universe.In this paper, using the Pantheon+ sample of SNe Ia which significantly enlarges the origin Pantheon sample size from 1048 to 1701, the 14 latest BAO measurements, and 32 OHD data, we constrain Hubble constant H 0 while leaving spatial curvature parameter Ω K and sound horizon r d completely free.In order to be model-independent, four cosmography approaches are taken, namely the Taylor expansions in terms of y 1 = z/(1 + z), y 2 = arctan z, y 3 = ln(1 + z) and the Padé approximants.We expand the transverse comoving distance and Hubble parameter with these four variables, and obtain the constraint results on H 0 , Ω K , and r d from the joint datasets.Then we adopt Bayesian evidence to determine which cosmography model gives the best approximation.
From the constraint results, we find that the newest Pantheon+ sample and updated BAO measurements could tighten the constraints on H 0 and Ω K obviously.And the values of H 0 obtained from the four cosmography models are in good agreement with P18 and lower than R22 by 2.3∼3.0σ.Meanwhile, a flat Universe is preferred by the observations under the ΛCDM+Ω K model and all cosmography approaches, except for y (3)   1 and y   (3) 2 .In addition, the constraints on r d are consistent with the estimate derived from the Planck CMB data within 1σ under all expansions, suggesting that r d is not sensitive to the approximations parameterizations.Furthermore, we find that there are anti-correlations between H 0 and Ω K or H 0 and r d , indicating smaller Ω K or r d would be helpful in alleviating the H 0 tension.Moreover, according to the Bayesian evidence, the ΛCDM is still the most favored model by the joint datasets.As for the cosmography models, the Padé approximant of order (2,2), the third order of y 3 and y 1 are the top three models that fit the datasets best, but the Taylor series in terms of y 2 are essentially excluded by cosmological observations.
As the final remarks, the release of H 0 tension in this paper is mainly caused by the increase in its uncertainty.If we take into account the full covariance of OHD (Moresco et al. 2020(Moresco et al. , 2022)), the H 0 tension with R22 will decrease from 2.3∼3.0σ to 1.0∼1.9σdue to the further increased uncertainty.At the same time, the cosmographic parameters can not be constrained effectively when expanding up to high order.This suggests that more high-precision data and additional probes are needed to measure the cosmic expansion due to the degeneracy with Ω K and cosmographic parameters.In the future, with the space-based GW detectors such as DECIGO (Seto et al. 2001;Kawamura et al. 2011), Taiji (Hu & Wu 2017), TianQin (TianQin collaboration 2016) and ET (Einstein Telescope) (Punturo et al. 2010) exploring more GW from BNS systems with EM counterparts, the constraints on H 0 and Ω K will be more accurate and tight.

APPENDIX
A. COSMOGRAPHIC EXPANSIONS In this section, all the cosmographic expansions up to the fifth order which we adopted in this study are reported.(i) the transverse comoving distance within the y 1 -variable model: (ii) the transverse comoving distance within the y 2 -variable model: (iii) the transverse comoving distance within the y 3 -variable model: (iv) P 12 approximation of the transverse comoving distance: (v) P 13 approximation of the transverse comoving distance: (vi) P 22 approximation of the transverse comoving distance: (vii) P 32 approximation of the transverse comoving distance: Bayesian evidence provides a statistical way to evaluate the performances of the cosmography models.For a given model M with parameter space θ and specific observational data d, the Bayesian evidence E is defined as where π(θ|M) is the prior of θ in model M, and p(d|θ, M) is the likelihood.Then, for the two models M i and M j , combing the Bayes theorem, the posterior probability is The ratio between posterior probabilities leads to the definition of the Bayes factor B ij , which is written in a logarithmic scale as Then, one can determine the strength of the preference for one of the competing models over the other by means of the Jeffreys scale listed in Table 4 (Trotta 2008;Kass et al. 1995).

C. DATA AND FITTING METHODS
In this section, we introduce the observational data and fitting methods used in our work.

C.1. SNe Ia
As the standard candle, SNe Ia is widely used to measure the cosmological luminosity distance.In our analysis, the largest SNe Ia sample Pantheon+ is used.Pantheon+ consists of 1701 light curves of 1550 spectroscopically confirmed SNe Ia across 18 different surveys, and covers the redshift range z ∈ [0.00122, 2.26137].The observed distance modulus of each SNe Ia in this sample is defined as Here, m * B is the observed peak magnitude in rest frame B-band, X 1 is the time stretching of the light-curve, C is the SNe Ia color at maximum brightness, M B is the absolute magnitude, α, β are two nuisance parameters, which should be fitted simultaneously with the cosmological parameters.And δ bias is a correction term to account for selection biases, δ host is the luminosity correction for residual correlations between the standardized brightness of an SNe Ia and the host-galaxy mass.In order to avoid the dependence of the nuisance parameters on the cosmological model, Kessler & Scolnic (2017) proposed a new method called BEAMS with bias corrections (BBC) to calibrate the SNe Ia, and the corrected apparent magnitude m * B,corr = m * B + αX 1 − βC − M B − δ bias + δ host for all the SNe Ia is reported in Popovic et al. (2021).Then the observed distance modulus is rewritten as On the other hand, introducing the Hubble-free luminosity d L (z) = H 0 D L (z), the theoretical value of distance modulus can be derived from where µ 0 = 42.38 − 5logh with h = H 0 /100 km s −1 Mpc −1 .Therefore, the χ 2 function for the Pantheon+ SNe Ia data can be written as where ∆µ i = m * B,corr − 5logd L (z) − N , and Cov is the covariance matrix, respectively.Here, N = M B + µ 0 , which can be marginalized over analytically with the method proposed in Conley et al. (2011).Finally, the χ 2 function is rewritten as Here, it should be pointed out that since that the sensitivity of peculiar velocities of SNe Ia is very large at z < 0.01 and thus the Hubble residual bias can not be negligible, the SNe Ia data with z > 0.01 is only used instead of the full sample, in our analysis.We refer the reader to the Ref. (Brout et al. (2022)) for further details on this issue.

C.2. BAO
BAO data have provided another way to probe the expansion rate and the large-scale properties of the Universe, which give information about the imprint in the primordial plasma.The BAO data used in this paper include the measurements from the Sloan Digital Sky Survey Data Release 7 (SDSS DR7) main galaxy sample (MGS) (Ross et al. 2015), the SDSS-III BOSS DR12 galaxy sample (Alam et al. 2017), the SDSS-IV eBOSS DR16 LRG sample (Gil-Marín et al. 2020), the SDSS-IV eBOSS DR16 ELG sample (de Mattia et al. 2021), the SDSS-VI eBOSS DR16 quasars (QSO) sample (Neveux et al. 2020), the eBOSS DR16 auto-correlations and cross-correlations of the Lyα absorption and quasars (du Mas des Bourboux et al. 2020).We summarize these data in Table 5.
We note that the likelihoods of MGS data, DR16 Lyα-Lyα and DR16 Lyα-QSO data used in this work cannot be well approximated by a Gaussian.Thus, their full likelihoods are used.The 4×4 covariance matrix for DR12 LRG data from Gil-Marín et al. ( 2020 In Table 5, the observable D V ≡ [czD M (z)/H(z)] 1/3 is the volume average distance and D H ≡ c/H(z) is the Hubble distance.The parameter r d is comoving sound horizon at the end of radiation drag epoch z d , shortly after recombination, when baryons decouple from the photons, where c s (z) is the sound speed of the photon-baryon fluid.In this paper, the parameter r d is treated as a free parameter.Finally, the BAO data are combined into a χ 2 -statistic Here, the DR12, DR16 LRG and DR16 QSO data χ 2 i are given in the form of χ 2 i = (w i − d i ) T • Cov i −1 • (w i − d i ).And, the vector d i is the observational data of the ith-type data set from Table 5, w i is the prediction for these vectors in a given cosmological model, and Cov i is the covariance matrix of different BAO data set.

C.3. OHD
In this paper, we use the most recent 32 OHD points summarized in Table 6.And the best-fitting parameters are obtained by minimizing this quantity where σ 2 Hi is the error of the i-th measurement.

Figure 1 .
Figure1.One-dimensional and two-dimensional marginalized distributions with 1σ and 2σ contours for the selected parameters under the third order.

Figure 2 .
Figure2.One-dimensional and two-dimensional marginalized distributions with 1σ and 2σ contours for the selected parameters under the fourth order.

Figure 3 .
Figure3.One-dimensional and two-dimensional marginalized distributions with 1σ and 2σ contours for the selected parameters under the fifth order.

Figure 4 .and
Figure 4.The reconstructed evolutions of H(z) with the flat ΛCDM, y (3) 1 , y (3) 3 and P22 models.The black lines show the center values and 1σ errors of H(z) for flat ΛCDM model.The dark green, blue and red lines with light bands show the center values and 1σ errors of H(z) for y (3) 1 , y (3) 3 and P22 models, respectively.

Table 1 .
The priors for the model parameters.

Table 2 .
Constrained cosmographic parameters by the SNe Ia+OHD+BAO dataset under various expansion orders within the 1σ confidence level.Here, '-' denotes there is no constraint result.

Table 3 .
The values of lnBij computed for the selected cosmography model Mi, where the reference scenario is the flat ΛCDM model (Mj).

Table 4 .
Revised Jeffreys scale quantifying the observational viability of any model Mi compared with some reference model Mj.

Table 5 .
BAO measurements used in our work.

Table 6 .
The compilation of OHD (in units of km s −1 Mpc −1 ) and their errors σH at redshift z.