An analytical late–Universe approach to the weaving of modern cosmology

Combining cosmological probes has consolidated the standard cosmological model with percent precision, but some tensions have recently emerged when certain parameters are estimated from the local or primordial Universe. Whether this is due to hidden systematics or it points towards new physics is still under debate, however, it is crucial to study as many probes as possible to cross–check the results with independent methods and provide additional pieces of information to the cosmological puzzle. In this work, by combining several late–Universe probes sampling the redshift range 0 < 𝑧 < 10, namely, Type Ia SuperNovae, Baryon Acoustic Oscillations, Cosmic Chronometers and Gamma–Ray Bursts, we aim to derive cosmological constraints independently of local or early–Universe anchors. To test the standard cosmological model and its various extensions, considering an evolving Dark Energy Equation of State and the curvature as a free parameter, we analyse each probe individually and all their possible permutations. Assuming a flat Λ CDM model, the full combination of probes provides 𝐻 0 = 67 . 2 + 3 . 4 − 3 . 2 km s − 1 Mpc − 1 and Ω 𝑚 = 0 . 325 ± 0 . 015 (68% C.L.). Considering a flat 𝑤 CDM model, we measure 𝑤 0 = − 0 . 91 + 0 . 07 − 0 . 08 (68% C.L.), while by relaxing the flatness assumption ( Λ CDM model, 95% C.L.) we obtain Ω 𝑘 = 0 . 125 + 0 . 167 − 0 . 165 . Finally, we analytically characterize the degeneracy directions and the relative orientation of the probes’ contours. By calculating the Figure–of–Merit, we quantify the synergies among independent methods, estimate the constraining power of each probe and identify which provides the best contribution to the inference process. Pending the new cosmological surveys, this study confirms the exigency for new emerging cosmological probes in the landscape of modern cosmology.


INTRODUCTION
At the beginning of the twenties of the third millennium, the discovery of the accelerated expansion of the Universe radically changed our understanding of its origin and evolution.To date, the physical nature of the energy driving this expansion -commonly called Dark Energy -and of most of the matter components in the Universe -referred as Cold Dark Matter (CDM) -still animate the debate within the scientific community.
Despite the fact that we are still unsure of the exact nature of these components, with the advanced technologies at our disposal we are able nowadays to precisely measure their effects on the observable Universe.For example, it has been observed that the total energy budget of the cosmos roughly matches its critical value, so that the flatness of the Universe is inferred with an extremely high level of precision (Planck Collaboration et al. 2021).The properties of the cosmological fluid responsible for the current accelerated expansion can also be measured, although to accurately study its possible temporal evolution we have to wait until more advanced surveys start to observe the sky, such as Euclid (Laureĳs et al. 2011), Vera Rubin ★ E-mail: fabrizio.cogato@inaf.it† E-mail: michele.moresco@unibo.it‡ E-mail: lorenzo.amati@inaf.it§ E-mail: a.cimatti@unibo.itObservatory (LSST Science Collaboration et al. 2009) and Nancy Grace Roman Space Telescope (Spergel et al. 2013).Moreover, another new fact emerged recently in the already complicated weaving of modern cosmology: the present-day expansion velocity (the Hubble constant,  0 ) presented a significant discrepancy when measured with independent probes, leading to the wellknown Hubble tension (Verde et al. 2019).Currently, this tension is mostly driven by the difference between the  0 estimated with Cepheids and Type Ia SuperNovae (SNe,  0 = 73.04 ± 1.04 km s −1 Mpc −1 , Riess et al. 2022) and with the analysis of Cosmic Microwave Background radiation from the ESA mission Planck (CMB,  0 = 67.36±0.54km s −1 Mpc −1 , Planck Collaboration et al. 2021), but some evidence suggest a dichotomy between late-Universe and early-Universe cosmological probes (Abdalla et al. 2022).To solve this tension, several alternative models have been proposed (Di Valentino et al. 2021) that, however, need a deeper and more detailed comparison with observations (Schöneberg et al. 2022a).Of course, before abandoning the standard cosmological model, a detailed assessment of the possible systematic effects is necessary.To address this delicate topic, one of the most studied methods is definitively the statistical combination of cosmological probes.
The combination of probes is hardly a new concept in cosmology.Since an additional cosmological component was measured by Riess et al. (1998) and Perlmutter et al. (1999), this method has been widely used to increase the precision of cosmological parameters inference (Bridle et al. 2001;Linder 2006;Davis et al. 2007;Lampeitl et al. 2010;Suzuki et al. 2012;Huterer & Shafer 2018;Schöneberg et al. 2022b;Brieden et al. 2023).Essentially, by combining the likelihoods obtained from independent probes, it is possible to exploit the different constraint powers of each observable to alleviate the degeneracy among cosmological parameters and thus improve our capability to investigate a wide range of models.Currently, only a few probes have been so extensively studied that they now represent the standard for any cosmological analysis, including Cosmic Microwave Background, Type Ia SuperNovae and Baryon Acoustic Oscillations.While the former samples the early stages of the cosmos ( ∼ 1100), the others are observed at low-redshift up to  ∼ 2.5.
In this paper, Type Ia SuperNovae and Baryon Acoustic Oscillations are analysed in synergy with two of the new emerging methods recently developed in the context of late-Universe cosmological probes (Moresco et al. 2022).On the one hand, Cosmic Chronometers provide independent measurements of the Hubble parameter and, thus, are able to probe the cosmic expansion history up to  ∼ 2. On the other, Gamma-Ray Bursts offer the opportunity to extract the cosmological signal up to  ∼ 10.Hence, the scope of this work is to provide a detailed study of these late-Universe probes, analysing the most recent data collections and combining them to obtain new and precise cosmological constraints that are completely independent of the standard approach in modern cosmology, namely, the Cosmic Microwave Background (Planck Collaboration et al. 2021) or the three-rung distance ladder (Brout et al. 2022).
Our approach essentially consists of a probe-by-probe combination through which we aim to monitor the potential systematic effects and compare the constraining power of each probe.In fact, our general purpose is to consolidate and extend the constraints from these late-Universe probes (especially for the emerging ones), ensure the robustness of the combination technique, and finally obtain reliable estimates of cosmological parameters that are noteworthy in the modern cosmology landscape.
Finally, we seek to define a mathematical framework to properly assess the complementarity between different cosmological probes and effectively exploit the synergies among their constraint powers.In practice, starting from the idea developed by Linder (2006), the peculiar degeneracy affecting each probe is evaluated by looking at the confidence contours' orientation on the parameter planes.From the largest eigenvalue of the covariance matrix associated with these contours, we trivially calculate the relative orientation of different probes and, by means of the Figure-of-Merit defined by Albrecht et al. (2006) and Wang (2008), evaluate the effects of probes combination.From this perspective, we derive some useful insights on the most effective combination to improve the inferential process as a function of the cosmological model.
To summarize, the physics explored by an experiment is strictly defined by the orientation of confidence contours in a particular parameter space.Hence, by fixing the parameter space and the probes to constrain it, the falsifiable models are closely related to the degeneracy directions affecting those probes.This is the key concept behind our work: different experiments are sensitive to different physical phenomena, therefore, it is crucial to correctly combine as many independent cosmological probes as possible in order to explore the whole parameter space and avoid any biased determination of the underlying cosmology.
The paper is organised as follows: in Section 2 we briefly describe the main features of each probe and model, as well as the Bayesian framework within which the analysis was carried out.Section 3 presents the results of the cosmological inference, with special attention to the impact and reliability of the probe-by-probe strategy.Then, in Section 4 an original study of synergies and complementar- ities between probes is presented, while Section 5 summarizes the main outcomes of our work.

Late-Universe data
Combining "standard rulers, candles and clocks" is an old recipe (Heavens et al. 2014) to extract the cosmological signal using exclusively low-redshift data with a model-independent approach (see also Moresco et al. 2016b andBenisty &Staicova 2021).Aiming to extend this pioneering work, we select the most updated and complete late-Universe datasets available today.More specifically, this work is based on the analysis of two primary probes, such as Type Ia SuperNovae and Baryon Acoustic Oscillations, which are combined with two new emerging methodologies, namely Cosmic Chronometers and Gamma-Ray Bursts, that have recently proven their strength and reliability (Moresco et al. 2022).
As shown in Figure 1, our dataset allows us to span a wide redshift range 0 <  < 2 with most of the probes, where Gamma-Ray Bursts observations extend it up to  ∼ 10.Hence, we extract the cosmological signal from a time frame currently not deeply explored in the literature.

Type Ia SuperNovae
Type Ia SuperNovae (SNe) are one of the most powerful and suggestive astrophysical events observed in the cosmos.The most recent and complete SNe collection is the so-called Pantheon+ sample (Brout et al. 2022), in which 1701 observations (0.001 <  < 2.26) of 1550 SNe are collected.Following Kessler & Scolnic (2017), the brightness of these objects is standardized through the definition of the distance modulus  as: where  corr  is the apparent magnitude taking into account several correction terms, such as selection bias, dust extinction, light-curve colour and stretch parameter, while the nuisance parameter  is the absolute magnitude of a fiducial SNe.Since  is related to the luminosity distance   through the following equation: by measuring the apparent magnitudes   of an SNe sample it is possible to constrain the cosmological parameters thanks to the generic expression: where Ω  is the curvature parameter,  () is the Hubble parameter and the function   is defined as: Due to the intrinsic degeneracy between  0 and the nuisance parameter , SNe are not able to constrain the expansion rate of the Universe.At the same time, SNe have been widely used in the past three decades to precisely infer the energy budget of the Universe.
Following the same procedure reported in Brout et al. (2022), we do not consider the nearby Hubble Diagram ( < 0.01) in order to avoid any systematics due to unmodeled peculiar velocities.Therefore, our final cosmological sample includes 1590 observations with their associated covariance matrix.

Baryon Acoustic Oscillations
Density fluctuations in the primordial photo-baryonic fluid affect the formation of large-scale structures in the more recent stages of the Universe.In fact, after the photo-baryonic decoupling ( ∼ 1100), the baryons at the boundary of the gravitational potential wells preferably cluster at a distance from the centre of the CDM haloes fixed by the length of the radius   of the sound horizon calculated at the time of decoupling.As a consequence of this phenomenon, in the latetime matter distribution it is most probably to find two galaxies separated by a distance   , and this particular feature, known as Baryon Acoustic Oscillations (BAOs), represents an incredible and versatile way to infer cosmological parameters.Considering   as a standard ruler, it is possible to measure different types of observables: the Hubble distance   () ≡ / (), the transverse comoving distance   () ≡   ()/(1 + ) or the volume-averaged distance   () ≡ [  () 2  ()] 1/3 , depending on the direction (radial and/or transverse) along which the BAOs signal is observed.
In this framework, we select a sample (0 <  < 2.4) of "BAOonly" measurements -collected in Alam et al. (2021) -carried out from different surveys by Ross et al. (2015), Alam et al. (2017), Gil-Marín et al. (2020), du Mas des Bourboux et al. (2020), de Mattia et al. (2021) and Hou et al. (2021).Here, "BAO-only" means that the measurements are uncalibrated, i.e., the value of   is assumed as a constant parameter to be inferred and is not calibrated through early-time Universe observations.The observables   ,   , and   are indeed scaled by a factor  −1  which, as explained by Brieden et al. (2023), is assumed to be a constant and isotropic length.
Hence, like SNe, BAOs are a very useful tool for measuring the energy budget of the Universe whilst they are insensitive to the Hubble constant  0 .

Cosmic Chronometers
The method of Cosmic Chronometers (CCs), introduced for the first time in Jimenez & Loeb (2002), measures the Hubble parameter as: (5) While the redshift  is a directly observable quantity, the variation of the look-back time d is determined by analysing the spectra or photometry of a particular population of massive and passive galaxies whose majority of stars were formed in the early stages of their evolution.The reason why this method has been increasingly used in cosmological analyses relies on its independence from cosmology.Assuming only the validity of the General Relativity, the Cosmological Principle, and the Weyl postulate, with CCs measurements we are able to directly track the history of the cosmic expansion  ().Thus, CCs are consolidating their role in the modern cosmological framework thanks to the relatively simple needs and the minimal cosmological assumptions they use.
As summarized in the review by Moresco et al. (2022), there are several methods to derive d measurements from the CCs observations: the full-spectrum fitting (see, e.g., Jiao et al. 2023;Tomasetti et al. 2023), the Lick indices analysis (see, e.g., Borghi et al. 2022), and the calibration of specific spectroscopic features (see, e.g., Moresco et al. 2012;Moresco 2015;Moresco et al. 2016a).More recently, it has been demonstrated the possibility of using photometric surveys to retrieve accurate  () measurements by selecting CCs with a Machine Learning approach (Jimenez et al. 2023).
Until now, this method has provided 32 measurements of the Hubble parameter over the redshift range 0.07 ≤  ≤ 1.965 (Moresco et al. 2022).

Gamma-Ray Bursts
Despite their great success in measuring cosmological parameters, the SNe observations do not go beyond  ∼ 2. Thanks to the observational efforts of the last decades, a new interesting kind of distance indicator emerged allowing us to investigate stages of the Universe out of the reach of the standard cosmological probes.This is the case of long Gamma-Ray Bursts (GRBs), the most energetic explosive events observed in the cosmos, produced by the core-collapse of peculiar massive stars (see, e.g., Piran 2004, Kumar & Zhang 2015and Levan et al. 2016).The combination of their origin, huge brightness and redshift distribution extending up to more than  ∼ 9 makes these phenomena very powerful probes for cosmology.GRBs are characterized by a prompt emission, lasting typically from a few seconds up to a few minutes, during which most of the energy is radiated in X/-rays, and a multi-wavelength afterglow emission spanning from -rays to radio and fading on time scales ranging from several hours to several days.
Certainly, GRBs are not standard candles.Also, the emission mechanisms at work, especially during the prompt phase, are not yet fully understood and only some aspects of the progenitors models are known.But, the discovery and deep investigation of empirical correlations between radiated energy (or peak luminosity) and spectral (or temporal) properties is consolidating the role of these events in modern observational cosmology (see, e.g., Amati & Della Valle 2013, Dainotti & Amati 2018, and for a recent review Moresco et al. 2022).In this context, the most investigated correlation for GRBs physics and cosmology is by far the "Amati relation" (Amati et al. 2002;Amati 2006;Amati et al. 2008;Amati & Della Valle 2013) between the photon energy at which the time-integrated F  spectrum of the X/-rays prompt emission peaks, i.e., the cosmological restframe peak energy  p,i =  p (1 + ), and the isotropic-equivalent radiated energy  iso .The  p,i - iso ("Amati") correlation takes the form: log  p,i keV =  +  log  iso 10 52 erg (6) where  and  are constants to be inferred.While  p,i is directly extracted from the measured prompt spectrum, the quantity  iso is related to the luminosity distance through the following relation: where  () is the Band function (Band et al. 1993) and  bolo is the measured bolometric fluence.Therefore, such a relation represents a remarkable tool to infer cosmological parameters since it correlates two observable quantities ( p and  bolo ), one of which ( bolo ) depends on them according to Equation 7.Although it is very highly significant, the Amati relation (Equation 6) is characterized by a scatter of the data around the best-fit power-law which significantly exceeds the level expected from Poissonian fluctuations and underestimated systematics in the measurement of the two quantities.Thus, an additional parameter  int has to be considered to take into account the intrinsic scatter of the correlation.Moreover, as already stated for BAOs and SNe, the intrinsic degeneracy between  0 and the nuisance parameter  does not allow us to infer the Hubble constant value from the "Amati" relation.
In our multi-probes cosmological analysis, we introduce GRBs through the Amati relation based on the updated  p,i - iso dataset by Amati et al. (in prep.), which includes 264 GRBs observations up to redshift values  ∼ 9 and is represented in Figure 2.
The values of  p,i and  iso of these GRBs are based on measurements of fluence and spectral parameters by the main GRBs detectors operated since the first measurements on GRBs redshifts in 1997 up to end of 2022: CGRO/BATSE, BeppoSAX/GRBM, HETE-2, Swift/BAT, Fermi/GBM and Konus-WIND.The sample was built by updating and integrating the early compilation by Amati et al. (2008) with the measurements reported in more recent systematic spectral analysis and catalogues (Amati et al. 2009;Gruber et al. 2011;Atteia et al. 2017;Tsvetkova et al. 2017;Fana Dirirsa et al. 2019;Katsukura et al. 2020;Minaev & Pozanenko 2020;Tsvetkova et al. 2021) and, for the latest GRBs, in GCN Circulars1 .In order to get the most robust possible sample, selection criteria were applied to the considered datasets, including the total duration of the event, the energy band of the detector, the signal-to-noise ratio of the fluence measurement, and the integration time of the GRBs spectrum with respect to the total burst duration.More details will be provided in Amati et al. (in prep.).

Cosmological models
To provide a detailed study on the extensions of the standard cosmological ΛCDM model, we focus our attention on a specific family of cosmological models that, in the CDM framework, investigate the null curvature hypothesis (Ω  ≡ 0) and the assumption of Dark Energy as a Cosmological Constant Λ ( ≡ −1).In the more general case, the DE EoS is described by the CPL model developed by Chevallier & Polarski (2001) and Linder (2003) in the early 2000s.This is a two-parameters model describing the variation with redshift  of the parameter  as: Without any a priori assumption on the curvature of the Universe, the evolution of the cosmic expansion is thus described through the following equation: where Here, even if not indicated, all the components to the energy budget Ω  are referred to their value at  = 0, e.g., Ω  ≡ Ω ,0 .
Therefore, three different DE EoS are derived by simply fixing particular parameters in Equation 8, namely: and, for each of them, the null (or flat) curvature case is investigated as well by imposing Ω Λ ≡ 1 − Ω  in Equation 9.

Bayesian framework
To constrain cosmological parameters we base our strategy on a statistical approach known as Markov Chain Monte Carlo (MCMC).In this framework, once the likelihood functions are defined and the posterior distributions are constructed through the priors definition (reported in Table 1), a random walk into the parameter space determines the best-fit values of the parameters.
In the present paper, we apply the widely used emcee software, i.e., a pure-Python implementation of the Goodman & Weare's MCMC sampler (Foreman-Mackey et al. 2013).We verified the convergence of our chains considering the Gelman-Rubin criterion with  < 0.01.
We consider more conservative priors on nuisance parameters w.r.t.those from the literature.In particular, the prior on  is roughly 5 times wider than the one imposed by Brout et al. (2022).Moreover, to avoid biases, the starting points of the MCMC have been selected with a random extraction within the prior intervals.
The best-fit value θ and the 68% (95%) Confidence Level (C.L.) of each parameter are extracted from the sampled distribution by marginalizing over the other ones and taking the 50 th percentile and 16 th (2.5 th ) and 84 th (97.5 th ) percentiles, respectively.

Single probes
As a first step, we test the relative strength and constraining power of each separate cosmological probe.
Table 1.Priors on the cosmological parameters explored with the MCMC method in the analyses involving BAOs, CCs, SNe and GRBs data.The symbol U indicates a uniform prior between the indicated extremes.

Cosmological Parameters
Nuisance Parameters

Baryon Acoustic Oscillations
As already discussed, our BAOs sample is composed of a diverse set of observations coming from a wide variety of spectroscopic surveys and cosmological quantities, namely   ()/  ,   ()/  and   ()/  .We first analysed the individual BAOs datasets and compared the derived constraints, to verify the consistency of the results.Alam et al. (2021) verified that the correlation between data is negligible and, moreover, provided a collection of covariance matrices for those sub-samples whose systematics were computed, namely, Alam et al. (2017), Gil-Marín et al. ( 2020) and Hou et al. (2021).To demonstrate how the combination of different BAOs probes allows us to significantly improve the cosmological measurements, in Figure 3 we show the constraints obtained from the main (separate) BAOs data, revealing how the synergies between various observables could be fundamental to achieving better precision in the inference process.Then, to analyse the full BAOs datasets we combine the likelihood functions of the different measurements as: ln L BAO = ln L Ross + ln L Alam + ln L Gil-Marín + + ln L du Mas de Bourboux + ln L de Mattia + ln L Hou (10) where the exact form of the likelihood function depends on whether or not the sub-sample is provided with a covariance matrix, that is: where the sum runs on the  measurements in the dataset.  ≡ (,   ) represents the theoretical value -depending on the cosmological parameters  -of the observable   ≡ (  ).Finally,   is the Gaussian-assumed error of the −th measurement and    is the matrix describing the variation of the −th measurement with respect to the −th one.
In Table 2 we report the main cosmological constraints from BAOs data, among which we highlight the significance of the inferred values of Ω  .On the other hand, this BAOs collection does not precisely constrain Ω Λ and the DE EoS parameters  0 and   , because of the strong degeneracy between parameters that quickly grows as the model dimension increases.Our results fully agree with those carried out by Alam et al. (2021) 2 .

Type Ia SuperNovae
As explained above, Brout et al. (2022) applied a cut to the full Pantheon+ sample and to the associated covariance matrix.Following their approach, our analysis is restricted to 1590 observations from which we construct the likelihood function as: (13) where   ℎ = 5 log 10   (, ) + , with  and  respectively representing the cosmological and nuisance parameters.Table 2 shows the main results of the SNe analysis.Through the most updated SNe sample, we derive precise constraints on both the density parameters Ω  and Ω Λ but also on the DE EoS parameter  0 .Obviously, such precision shows a decreasing trend as the complexity of the model increases and the cosmological parameters become more degenerate with each other.As already mentioned and widely demonstrated in the literature, SNe do not constrain the Hubble constant as this parameter is strongly degenerate with .A more detailed analysis of this feature is given in Section 3.3.2.
The agreement between our results and those obtained by the original authors 3 further demonstrates the reliability of our implementation of the Bayesian approach.Moreover, we extend their analysis also to the CDM and  0   CDM models, finding no significant deviations from the standard cosmological model.

Cosmic Chronometers
As proof of the great effort in the last decade to make this probe increasingly reliable and robust, Moresco et al. (2020) compute the full covariance matrix for CCs taking into account several potential systematic errors, such as young stellar population, uncertainty in the stellar population synthesis models, and estimate of the stellar metallicity.Exploiting the framework built into that work, we write the likelihood function as: ) from which it is possible to infer cosmological parameters considering both statistical and systematic uncertainties.The analyses of the ΛCDM models show how this method could provide  0 measurements with a precision up to ∼ 9%, although far from the one obtained by the current main probes -such as local Cepheid variables (SH0ES, Riess et al. 2022;Brout et al. 2022) or CMB (Planck Collaboration et al. 2021).Unfortunately, the current number of observations is not enough to keep such good precision when increasing the number of free parameters of the model.Therefore, when used to probe more complicated extensions, this method becomes less sensitive to the physical properties of the Universe, such as the curvature or the DE EoS.It is worth noting how CCs are more sensitive to the degeneracy between Ω  and Ω Λ than to the one between Ω  and the DE EoS parameters.In fact, the inferred value of Ω  in the flat  0   CDM case shows a significantly higher level of precision compared to the one obtained in the ΛCDM model.The results reported in Table 2 are in agreement with those coming from the literature (e.g., Moresco et al. 2016a,b andMoresco et al. 2022).

Gamma-Ray Bursts
To extract the cosmological signal from this probe, we exploit the Amati relation based on the most updated and robust dataset, as described in Section 2.1.Following the same approach adopted by Amati et al. (2008) and Amati & Della Valle (2013), we construct the likelihood function as: where  p,i =   p,i ln(10) • p,i  iso =   iso ln(10) • iso and   p,i (  iso ) represent the uncertainty on the measurement of  p,i ( iso ).We refer the reader to Reichart et al. (2001) for more details about the definition of the likelihood function.Since the values of  iso shown in Figure 2 are calculated from Equation 7 by assuming a fiducial flat ΛCDM cosmology with  0 = 70 km s −1 Mpc −1 and Ω  = 0.3, we impose: where  are the parameters of the generic model and  fid iso are those measurements shown in Figure 2. In this way, we are thus able to avoid any problem of circularity and to constrain the cosmological parameters in a way as reliable as robust.
The results from the analyses of GRBs are reported in the bottom panel of Table 2. Our work extends the ones from the literature by investigating a wider class of cosmological models, finding results fully consistent with those obtained by Amati et al. (2008), Amati & Della Valle (2013) and Moresco et al. (2022).Considering the flat ΛCDM model we further confirm that also GRBs indicate a value of Ω  around 0.3, even if the level of precision is not comparable with the other probes.We find, as expected, that for more complicated models the constraining power of this probe decreases.At the same time, the prospect of probing the cosmos up to a very high redshift ( ≫ 1) makes GRBs one of the most promising emerging probes.Indeed, while on the one side, they allow us to study a phase of the Universe not yet sampled by standard probes, on the other side, these events reveal completely different degeneracy directions compared to other probes, which can be exploited to break degeneracies among parameters and achieve more accurate constraints.For example, looking at the orientation of the two-dimensional contours in the Ω  − Ω Λ plane (ΛCDM) reported in Figure 5b, it should be evident the different orientations of the confidence contours of GRBs and SNe.Although the cosmological observable -the luminosity distance   () -is the same, such a difference in the orientation of the contour is essentially due to the different redshift distributions (see Figure 1), as well as to the distinct nuisance parameters these two probes adopted to constrain the cosmological ones.

Towards the full probes combination
After analysing the individual probes separately, we apply a probeby-probe combination technique to highlight the different constraining powers attainable from the various combinations of cosmological probes.
In Figure 4 the distribution of the best-fit values (and 68% C.L.) obtained with this approach is shown as a function of the entire set of cosmological parameters.Two things catch immediately the eye.First, we note that there is a small variance and good agreement between the constraints derived from different combinations of probes.Then, it is evident how the error bars decrease as the number of probes involved in the combination increases.An exception can be found in the  0 column where the red dots (triplets) corresponding to the BAO+SN+GRB combination show an uncertainty significantly greater than the yellow dots (doublets).However, this behaviour is expected since, as already highlighted, we are considering probes mostly insensitive to the Hubble constant.
This analysis is fundamental to assess the degeneracies affecting each probe (and their combination) and, also, how these features impact cosmological inference and mutate according to the models.Indeed, the dispersion and the uncertainties of measurements increase as the cosmological model becomes more complicated, as can be seen from the comparison of the flat ΛCDM and  0   CDM distributions of the Ω  measurements.
As can be also deduced from Table 2, we note a marginal trend of our results (at 68% C.L.) of preferring solutions with hyperbolic geometry (Ω  > 0, Ω Λ < 0.7) when the curvature of the Universe is left free, as also seen by Moresco et al. (2016b); Brout et al. (2022); Alam et al. (2017).However, at 95% C.L. the curvature constraints are compatible with a flat geometry.Furthermore, regarding the  0 distribution, it is worth noticing that increasing the number of combined probes the inferred value tends to  0 ∼ −1.Moreover, we note that there are combinations of probes that are less constraining (e.g., CC+GRB) but, also in these cases, the results agree with a Dark Energy component in the form of a Cosmological Constant Λ.
Finally, looking at the   distribution and the large error bars, it should be evident that current cosmological data are not able to accurately constrain this parameter, even if they exclude a significant part of the parameter space, namely,   < −1 and   > 2.  3).Regions with diagonal line patterns indicate those cosmologies not constraining those specific parameters, e.g., the Ω Λ parameter is not constrained by "flat" models.

The full combination of probes
From the probe-by-probe combination analysis, we found a generic consistency between the results obtained from the various methods, so the next step is to simultaneously combine all probes to further increase the precision of the inference.
The likelihood function is constructed through the standard approach: where lnL BAO , lnL SN , lnL CC and lnL GRB are built from Equation 10, 13, 14 and 15, respectively.The cosmological constraints from the BAO+CC+SN+GRB combination (full, in short) are reported in Table 3, while in Figure 5 we show some specific projections of the sampled posterior distributions.The combination of these late-Universe probes shows immediately its potential in the modern cosmological context.Assuming a flat ΛCDM model, we precisely constrain (at 68% C.L.) the standard cosmological parameters as: showing remarkable precision on both parameters, especially if we consider that they are obtained without any restrictive prior.As suggested by Figure 5a, CC is the only considered probe sensitive to the Hubble constant.For this reason, it is maybe more interesting to note how by combining CCs with BAOs, SNe, and GRBs we are able to improve the percentage precision on  0 from 8% to 5%.On the other hand, the noticeable result achieved for Ω  is essentially driven by the high precision of BAOs and SNe constraints.
We note also the tendency of the BAO+CC+SN+GRB combination to prefer a Universe with negative curvature ( = −1 ⇒ Ω  > 0), but this trend is not significant at 95% C.L, and therefore does not represent a statistically significant deviation from the standard paradigm.Finally, considering the most complicated model, i.e.,  0   CDM, the DE EoS seems to be consistent with a Cosmological Constant Λ, even though Figure 5f demonstrates how the one-dimensional marginalized distribution of   appears particularly asymmetric, with the peak of the distribution around   ∼ 1.

Comparison with literature results
It is very interesting to compare our results to the ones in the literature considering other (combinations of) probes.
First of all, our measurement of the Hubble constant in the flat ΛCDM cosmology is fully consistent with the early-Universe (TT,TE,EE+lowE+lensing) inference derived by Planck Collaboration et al. ( 2021   variables and Pantheon+ SNe: Although we do not achieve enough precision to disentangle between the  0 current values of the Hubble tension, we note that our results prefer solutions leaning towards a CMB-like measurement.Moreover, they show remarkable stability (in terms of precision  0 ) as the model dimension increases.Since it is known the strong model dependence of the constraints from CMB power spectrum (TT,TE,EE+lowE), we also compare our  0 measurement to a cosmological independent result, such as the combination between CMB lensing and BAOs (Planck Collaboration et al. 2021).With a prior on the baryon density (Ω  ℎ 2 = 0.0222 ± 0.0005), the standard cosmological model is constrained (68% C.L.) as: with respect to which our late-Universe analysis, without applying any restrictive prior, is extremely competitive for what concerns both  0 and Ω  .
By further looking at the literature results, it clearly emerges how our approach represents a worthy rival to the standard probes usually combined in modern cosmology.In fact, by assuming a flat CDM model and combining CMB (TT,TE,EE+lowE) with some BAOs data -including Ross et al. (2015), Alam et al. (2017) In all these cases, the CMB power spectrum acts as the main driver of the high precision achieved.Especially for high dimensional models, it is hence essential to emphasize once again the extreme relevance of a combination of late-Universe probes, since this method has proven its reliability in providing cosmological results that are competitive with those from standard cosmological literature.Such a relevant outcome is essentially driven by the different cosmological dependencies depicted by each probe selected for this work.As shown in Figure 5, according to which observable and empirical relationship they make use of to constrain cosmological parameters, distinct probes show different oriented contours.So, the combination of many different probes can exploit several constraining powers in order to break (or, at least, to alleviate) some of the parametric degenerations characterizing every cosmological model.The striking examples of this behaviour are observed in Figure 5c and 5e where the pronounced orthogonality between BAOs and SNe (supported also by the contribution of CCs and GRBs constraints) significantly increases the precision on the cosmological parameter estimation.All these considerations are underlying the fundamental role played by probes combination in the modern cosmological context and its innovative potential to address some of the still open questions.As we will study and deeply analyse in Section 4, in a scenario dominated by strong parametric degenerations it is worth implement- ing the combination of probes as a tool to break them and improve the inference of cosmological parameters.

A nuisance parameters' perspective on the Hubble tension
In this work, in addition to the cosmological ones, we also constrain some nuisance parameters, such as the absolute magnitude  of a fiducial SNe, the radius   of the sound horizon evaluated at the decoupling epoch, as well as the slope  and the intercept  of the "Amati" relation.The dependence of the observables on these quantities introduces a degeneracy with the Hubble constant.For this reason, it is possible to constrain  0 only when SNe, BAOs or GRBs are combined with  0 -sensitive probes, such as the CCs (Moresco et al. 2022).Another way to break the degeneracy is to use an external calibrator able to constrain the nuisance parameters, namely, the observation of local Cepheid variables (Brout et al. 2022) for calibrating the absolute magnitude  of SNe, or the CMB power spectrum (Planck Collaboration et al. 2021) to directly measure the length of   .
Starting from this consideration, we derive two important consequences.First of all, including CC data allows us to obtain more precise measurements of nuisance parameters and to calibrate the SNe, BAOs and GRBs methods.The effects of the CC-calibration on other probes can be seen by observing how the precision on nuisance parameters increases passing from the individual probes analysis (Table 2) to the full combination one (Table 3).At the same time, considering the full probes combination, we notice that the inferred precision on  0 is almost constant, contrary to what would be expected as the dimension of the analysed cosmological model increases.It is interesting to notice (looking at Table 2) that similar behaviour is observed also for nuisance parameters from individual probe analyses (i.e.,  for SNe,   for BAOs,  for GRBs).This effect points to the fact that these nuisance parameters do not depend on the cosmological model, but only on observational data.Therefore, a possible explanation would be that in the case of the full combination of probes, the uncertainties on  0 cannot be reduced below a given threshold given by the intrinsic scatter in the SNe, BAO and GRB data, representing a plateau in the  0 error (see Table 3).2022) may also be read as a discrepancy between the inferred values of the nuisance parameters of the BAOs and SNe methods, i.e.,   and , respectively.From this perspective, Figure 6 shows how our combination of late-Universe probes seems to prefer values of   and  consistent with those obtained with the inverse distance ladder, i.e., through a calibration from CMB data (Planck Collaboration et al. 2021), rather than those found from the three-rung distance ladder in the local Universe (Brout et al. 2022).Indeed, considering a flat ΛCDM model, from the full combination we measured these two parameters as: showing a small deviation to results obtained by Alam et al. (2021): and by Brout et al. (2022): from the calibration of the sound horizon   and the absolute magnitude  through the three-rung distance ladder.

SYNERGIES AND COMPLEMENTARITIES
We conclude our analysis with an analytical quantification of the contributions of different datasets to the inference of cosmological parameters.Our intention is to figure out the role of each probe in the improvement of precision of the various parameters, quantify the complementarities between them, and study how to exploit synergies in order to maximize the constraint powers of cosmological data.

Degeneracy and orthogonality
As implicitly assumed, every dataset has a peculiar constraining power reflected in different orientations of the contours that each probe shows in specific planes (see Figure 5).
Here, we define a method to describe the geometric properties of contours within the parameter space.The first step is the determination of the degeneracy direction, i.e., the orientation of the contours on a defined parameter plane.For this purpose, we extract from the MCMC samples the marginalised 2D distribution of the posterior for a given combination of cosmological parameters and then estimate the counter-clockwise angle  of the contours as the inclination of the major axis of the ellipse defined by the eigenvalues and eigenvectors of the 2D contour.We then define the degeneracy direction as the line -passing through the 50−th percentile of the distribution -with slope: It is then straightforward to determine the relative orientation of the contours of different probes by considering the acute angle  between two different degeneracy directions.As a result, geometrically speaking, we quantify the orthogonality factor  between two probes by defining: so that two perfectly orthogonal (parallel) contours are characterized by  = 1 ( = 0).Because of the different ranges spanned by the various cosmological parameters (see Table 1), when calculating the relative orientation , we normalize the sampled distributions to one to avoid scale bias and improve the accuracy on .In practice, while we do nothing for the energy density parameters, we renormalise the range of  0 and  0 to one (equivalently to considering ℎ).Specifically: where  MC is the original MCMC distribution,  min and  max are the extremes of the prior intervals.
As an example, in Figure 7, the relative orientation  of the SNe and CCs contours is shown in the  0 − Ω  plane for the flat ΛCDM model.As reported in Table 4, these two probes have an orthogonality factor  = 0.35 corresponding to a relative orientation  = 31.4• .
Obviously, because the accuracy of this estimate strongly depends on the calculation of the covariance ellipse from the MCMC sampling, the method implemented here is not sensitive to distributions significantly deviating from the multivariate Gaussian.For this reason, in Table 4 we show only those parameter planes in which contours are approximately ellipsoidal and not extremely widespread, i.e., flat ΛCDM, ΛCDM and flat CDM (see Figure 5).

Figure-of-Merit
Since the evidence of a new component emerged in the cosmological energy budget, the scientific community has developed several tools to quantify how future experiments could improve the measurement of Dark Energy properties (Linder 2006;Albrecht et al. 2006;Albrecht & Bernstein 2007;Wang 2008;Sartoris et al. 2016).Considering a cosmological model C described by a set of N parameters  C ≡ {  1 ,  2 , ...,   }, the Figure-of-Merit (FoM) associated  to the constraints provided by a given probe is defined as: where  is the covariance matrix between the measured  C parameters.Generally speaking, if the posterior is a multivariate Gaussian distribution, the volume that it occupies within the entire cosmological parameter space is proportional to the inverse of the FoM.As suggested by Linder (2006), the volume of the probability distribution represents a suitable way to evaluate how different probes constrain the parameter space.But we shall keep in mind that it can be used as a proxy of the precision achieved by a chosen dataset only in the Gaussian approximation, i.e., as long as the posterior distribution is a multivariate Gaussian.
Here, we focus on the main set of cosmological parameters relevant to our analysis.In particular, as already done for the contours shown Table 5. Figure-of-Merit obtained for particular planes within the parameter space by BAOs, CCs, SNe, GRBs and their combination BAO+CC+SN+GRB. in Figure 5, the results obtained for those specific parameters planes are shown in Table 5.

Figure-of-Merit
The most relevant results we obtained are the following: 1. Through this approach, we quantified the gain in terms of FoM that the full combination achieves w.r.t. the individual probes.Indeed, the FoM values retrieved by the BAO+CC+SN+GRB inference are systematically higher than the sum of FoM from the individual datasets, with a maximum gain of about 3 times for the flat CDM and  0   CDM cases.
2. The capability of CCs to measure  0 makes them the most powerful probe (among the one considered) in the  0 − Ω  plane (flat ΛCDM), even if the constraint on Ω  is weaker w.r.t. the ones from BAOs and SNe (see Table 2).Although SNe are not able to constrain  0 , the FoM of SNe and CCs are essentially the same, emphasizing the fundamental contribution that SNe provide for those analyses aiming to measure Ω  precisely.
3. BAOs and SNe are undoubtedly the most powerful probes that one should exploit whatever the cosmological model considered.While in the flat ΛCDM, ΛCDM and flat CDM models, SNe represent by far the most constraining probe; in the other models, the situation is reversed and BAOs achieve higher precision.In terms of contribution to the FoM, this makes the BAO+SN combination dominant on all parameter planes, but flat ΛCDM.For example, in the ΛCDM case, this combination reaches around 90% of the FoM achieved by BAO+CC+SN+GRB.
4. The contribution of GRBs to the FoM in the  0 − Ω  plan (flat ΛCDM) is minimal.This can also be confirmed by looking at the FoM obtained with the BAO+CC+SN combination, which is roughly the same as the full combination.Conversely, their contribution increases for the other models, with a maximum of 25.6 in the Ω  − Ω Λ plane (ΛCDM).In fact, the wider the range of redshift sampled by the probe, the more sensitive the probe is to changes in density and DE EoS parameters, remarking how the extension in redshift of GRBs observations (see Figure 1) is fundamental to better constrain both the energy budget of the Universe and the Dark Energy properties.
However, we emphasize that the Gaussian approximation hypothesis is not consistent with the more complicated models (CDM, flat  0   CDM and  0   CDM) where one-dimensional marginal distributions begin to be strongly asymmetric (see Figure 5).Therefore, for these models, the FoM values reported in Table 5 should be taken as an indication of the relative gain of the full combination, but they are not informative of the volume occupied by the posteriors in the parameter space.

CONCLUSIONS
In this work, we demonstrated how the combination of some late-Universe cosmological probes, such as Baryon Acoustic Oscillations (BAOs), Cosmic Chronometers (CCs), Type Ia Supernovae (SNe), and Gamma-Ray Bursts (GRBs), could be profitably exploited in modern cosmological analysis to obtain precise constraints completely independent of early-Universe or local anchors, such as the Cosmic Microwave Background or the three-rung distance ladder.
The main results can be summarized as follows.
• Firstly, we developed, tested and validated a Bayesian code to derive the posterior distribution of a variety of cosmological parameters (considering different cosmological models of increasing complexity) considering several different probes.Once we demonstrated the reliability of our analysis by reproducing the main results coming from the literature for each one of the individual probes, we applied a probe-by-probe combination procedure aiming at extracting the cosmological signal from any combination achievable with our selection of probes.This approach has been useful to ensure that no remarkable biases and no systematic effects were affecting the Bayesian analyses.
• Surely, BAOs and SNe, as the most robust and widely studied probes, are found to contribute significantly to most of the results.For example, in the ΛCDM case, they are respectively found to be responsible for ∼ 30% and ∼ 50% of the FoM and reach around 90% when combined together.At the same time, they are restricted to a relatively small redshift range ( < 1.5) and, if uncalibrated, they are not sensitive to the Hubble constant  0 .From this perspective, it is not surprising the need to combine additional probes and better capitalise different sensitivities to cosmological parameters in order to break degeneracies and go beyond individual constraint powers.In this study, CCs and GRBs are the new emerging probes strongly contributing to thinning out that fog not allowing to see possible deviations to the standard cosmology.
• CCs offer a cosmology-independent constraint on  0 , Ω  , and curvature, even if with a smaller FoM compared to BAOs and SNe.It is however interesting to underline that their constraints allowed us to break the degeneracies between the Hubble constant and some of the nuisance parameters of other probes, like  for SNe and   for BAOs, and enhance the constraints on all of them.For example, in a flat ΛCDM model, by including CCs in the combination of BAOs and SNe, the uncertainty on  0 decreases from about 5.5 (CC-only) to about 3.5 km/s/Mpc (BAO+CC+SN) while the determination of Ω  improves from 0.330 +0.078  −0.061 to 0.326 +0.015 −0.014 .• The wide redshift range covered by GRBs (0 <  < 10) has been pivotal to improving the constraints on the energy composition of the Universe and the Dark Energy Equation of State, as these quantities are proven to be better measured by having a larger redshift leverage.In particular, adding GRBs to the combination of the other probes enhances the constraints making the posterior distributions more Gaussian.For example, considering a flat CDM model, while BAO+CC+SN measures  0 = −0.906+0.068  −0.084 , the full combination provides  0 = −0.907+0.072  −0.077 .• By analysing the full combination of probes (BAO+CC+SN+GRB) that sample an extremely wide range of cosmic redshift (0 <  < 10), we derived interesting and precise constraints on the properties of the Universe.Without imposing any restrictive priors, some of the main results (at 68 % C.L.) are: • By relaxing the assumption of flatness, we obtained remarkable but no statistically significant deviations from the null curvature paradigm.Indeed, such a combination seems always to prefer (at 68% C.L.) a hyperbolic Universe: Ω  = 0.125 ± 0.082 (ΛCDM), Ω  = 0.122±0.113(CDM) and Ω  = 0.148±0.115( 0   CDM).
• As a matter of fact, current data are not ready (yet, see Moresco et al. 2022) to provide competitive constraints on  0 to solve the well-known tension.At the same time, the great contribution to  0 constraints provided by CCs reveals an important feature of modern cosmological analyses.When combining probes that lack direct constraining power on  0 , careful consideration must be given to the impact of nuisance parameters used to calibrate their cosmological observables.Notably, the full combination analysis reveals that the precision on  0 is inherently constrained by the inherent scatter within data derived from probes insensitive to  0 , like BAOs, SNe, and GRBs.Consequently, probes like CCs play an essential role in calibrating these  0 -insensitive probes, if one wants to explore independent kinds of calibrators.At the same time, precisely and accurately measuring cosmological observables like the Hubble diagram or the Amati relation directly impacts the data scattering, thereby enhancing the precision of  0 measurements obtained through the combination of both sensitive and insensitive probes.
• For the first time in literature, by introducing the orthogonality factor  between contours, we quantified the degeneracy directions and the relative orientation of probes on each parameter plane.Studying the Figure-of-Merit, we provided a clear quantification of the constraint powers of each probe.
In conclusion, the precision achieved by combining a limited sample of late-Universe probes has shown extreme competitiveness compared to the combination of probes sampling from both the earlyand late-Universe, such as CMB, BAOs and SNe.This is extremely important in view of currently ongoing and future large surveys, like Euclid (Laureĳs et al. 2011), Vera Rubin Observatory (LSST Science Collaboration et al. 2009) and Nancy Grace Roman Space Telescope (Spergel et al. 2013), which data exploitation will further increase the effectiveness of the inference process and, thus, facilitate the transition from precision cosmology to accurate cosmology (Peebles 2002;Verde 2014).
We acknowledge here two possible follow-ups of this work for the next future, that are currently beyond the scope of this paper.First, the amplitude of matter perturbation (through the parameter  8 ) could be included amongst the late-Universe analyses by adding the Redshift Space Distortion probe (Benisty 2021;Aubert et al. 2022) to the combination process.This step would be fundamental to studying the degeneration of  8 with the other cosmological parameters and could provide an independent perspective to the current tension on  8 ≡  8 (Ω  /0.3) 0.5 (Abdalla et al. 2022;Poulin et al. 2023).Second, studying other emergent cosmological probes, especially the ones able to measure the Hubble constant independently -such as Gravitational Waves (Abbott et al. 2017) or Kilonovae (Sneppen et al. 2023), would certainly be an excellent way to verify the systematics involved in the inference process, and contemporaneously enhance the precision of the estimation of the current expansion rate of the Universe.

Figure 1 .
Figure 1.Redshift distribution of data collected in this work.

Figure 4 .
Figure 4. Distribution of the best-fit value (and 68% C.L.) of cosmological parameters inferred from the progressive combination of probes.Each column shows the constraints obtained for the different cosmological parameters as a function of the cosmological model analysed.Dots with different colours indicate a different number of probes involved in the analysis: individual probes (1x) in dark blue, a couple of probes (2x) in orange, and a triplet of probes (3x) in red.In particular, from top to bottom, each panel (i.e., each model) reports the measurements obtained by GRB, SN, CC, BAO, SN+GRB, CC+GRB, BAO+GRB, CC+SN, BAO+SN, BAO+CC, BAO+CC+SN, CC+SN+GRB, BAO+SN+GRB and BAO+CC+GRB.Vertical shaded regions report the value of cosmological parameters inferred from BAO+CC+SN+GRB in the most general  0   CDM model (see Table3).Regions with diagonal line patterns indicate those cosmologies not constraining those specific parameters, e.g., the Ω Λ parameter is not constrained by "flat" models.
Secondly, as already established by several works, such as Efstathiou (2021), Alam et al. (2021) or Favale et al. (2023), the difference between values of the current expansion rate of the Universe found by Planck Collaboration et al. (2021) and Riess et al. (

Figure 7 .
Figure 7. Example of the calculation of degeneracy directions and relative orientation .The contours obtained from the CCs and SNe analyses of the flat ΛCDM model are shown in the  0 − Ω  plane.

Table 2 .
from the analyses of BAOs, SNe, CCs, and GRBs data.Note that   and  0 values are respectively in units of Mpc and km s −1 Mpc −1 .Empty values correspond to parameters not constrained by the corresponding cosmological model.The values of Ω  (≡ 1 − Ω  − Ω Λ ) are directly extracted by combining the marginalized distributions of Ω  and Ω Λ .

Table 3 .
Constraints on cosmological and nuisance parameters (best-fit and 68% C.L. values) from the combination of BAOs, CCs, SNe, and GRBs data.Note that   and  0 values are respectively in units of Mpc and km s −1 Mpc −1 .Empty values correspond to parameters not constrained by the related cosmological model.The values of Ω  (≡ 1 − Ω  − Ω Λ ) are directly extracted by combining the marginalized distributions of Ω  and Ω Λ .

Table 4 .
Relative orientation  and orthogonality  between BAOs, CCs, SN and GRBs on particular parameter planes within the parameter space.