**4. Multidimensional Binned Analysis with SNe Ia and BAOs**

To investigate the *H*<sup>0</sup> tension through the SNe Ia and BAOs data, we combine the *χ* 2 Equations (3) and (7) to obtain the total *χ* 2

$$
\chi^2 = \frac{1}{2}\chi\_{\text{SN}}^2 + \frac{1}{2}\chi\_{\text{BAO}}^2 \tag{8}
$$

In our work, we combine each SNe bin with only 1 BAO data point which has a redshift value within the SNe bin: this approach of using one BAO comes from [18]. In this way we do not have the problem of a different number of BAOs in different bins. Through Equation (8), we investigate if a redshift evolution of *H*0(*z*) is present, obtaining it from the binning of SNe Ia+BAOs considering three bins with the ΛCDM and *w*0*wa*CDM models. A feasibility study done in [51] performed with different bins selections has highlighted how the maximum number of bins in which the PS should be divided is 3, otherwise the statistical fluctuations would dominate on a multi-dimensional analysis, leading to relatively large uncertainties which would mask any evolving trend, if present. Furthermore, for the same reason, it is not advisable to leave free to vary more than two parameters at the same time, thus in the current section, we will analyze the behavior of *H*<sup>0</sup> in three bins when it is varied together with a second cosmological parameter. The same considerations make necessary the choice of more tight priors since we are basing the current analysis on the prior knowledge, avoiding the degeneracies among the parameter space, and letting the priors have more weight in the process of posteriors estimation. Differently from [51], for the ΛCDM model, we will let the parameters *H*<sup>0</sup> and Ω0*<sup>m</sup>* vary simultaneously, while in the *w*0*wa*CDM model the varying parameters are *H*<sup>0</sup> and *wa*. We decided to leave *w<sup>a</sup>* free to vary since, according to the CPL parametrization, *w<sup>a</sup>* gives direct information about the evolution of the *w*(*z*) while *w*<sup>0</sup> is considered a constant in the same model. Concerning the fiducial values and the priors assignment for the MCMC computations, we apply Gaussian priors with mean equal to the central values of Ω0*<sup>m</sup>* = 0.298 ± 0.022 and *H*<sup>0</sup> = 70.393 ± 1.079 for Ω0*<sup>m</sup>* and *H*0, respectively, and with 1 *σ* = 2 ∗ 0.022 and 1 *σ* = 2 ∗ 1.079 for Ω0*<sup>m</sup>* and *H*0,

respectively. In summary, to draw the Gaussian priors, we consider the mean value of the parameters as the expected one of the Gaussian distribution and we double the *σ* value which is then considered the new standard deviation for the distribution. Concerning the *w*0*wa*CDM model, we fix *w*<sup>0</sup> = −0.905 and we consider the priors on *w<sup>a</sup>* with the mean = −0.129 taken from Table 13 of [284], while 1 *σ* = is the 20% of its central value. Such an assumption with small prior is needed since we need to assume that *w*(*z*) > −1.168 as the value tabulated in [284]. Besides, since we are here dealing with standard cosmologies, with this constraint we are avoiding some of the phantom Dark Energy models.

After the *χ* <sup>2</sup> minimization for each bin, we perform a MCMC simulation to draw the mean value of *H*<sup>0</sup> and its uncertainty. Once *H*<sup>0</sup> is obtained for each bin, we perform a fit of *H*<sup>0</sup> using a simple function largely employed to characterize the evolution of many astrophysical objects, such as GRBs and quasars [17,29,31,35,345–349]. More specifically, the fitting of *H*<sup>0</sup> is given by

$$f(z) = H\_0(z) = \frac{\mathcal{H}\_0}{(1+z)^{\eta}},\tag{9}$$

in which H<sup>0</sup> and *η* are the fitting parameters. The former H<sup>0</sup> ≡ *H*<sup>0</sup> at *z* = 0, while the latter *η* coefficient describes a possible evolutionary trend of *H*0. We consider the 68% confidence interval at, namely 1 *σ* uncertainty.

In the current treatment, we consider the calibration of the PS with *H*<sup>0</sup> = 70 as provided by [284]. Results are presented in the panels of Table 1. We here stress that the fiducial magnitude value is assumed to be *M* = −19.35 for each SNe bin, thus it will not be mentioned in the same Table. All the uncertainties in the tables in this paper are in 1 *σ*. As reported in the upper half of Table 1, namely with the ΛCDM model, if we do not include the BAOs then the *η* coefficient is compatible with 0 in 2.0 *σ* for the three bins case. When we introduce the BAOs within the ΛCDM model, we observe again a reduction of the *η*/*σ<sup>η</sup>* ratio for three bins down to 1.2. Concerning the lower half of Table 1 with the *w*0*wa*CDM model, when BAOs are not included we have *η* non compatible with 0 in 5.7 *σ* and, including the BAOs, the compatibility with 0 is given in 5.8 *σ*. The increasing of the ratio *η*/*σ<sup>η</sup>* is observed when BAOs are added in the case of *w*0*wa*CDM model in three bins. The results can be visualized in Figure 1. Comparing the *η*/*σ<sup>η</sup>* ratios with the ones reported in [51] (Table 1) we have that for the ΛCDM model the current *η* values are compatible in 1 *σ* with the *α* reported in [51], while the *η* estimated in the *w*0*wa*CDM model are compatible in 3 *σ* with the *α* values in the same reference paper.

**Table 1. Upper half.** Fit parameters of *H*0(*z*) for three bins (flat ΛCDM model, varying *H*<sup>0</sup> and Ω0*m*) in the cases with SNe only and with the SNe + BAOs contribution. The columns are: (1) the number of bins; (2) H0, (3) *η*; (4) how many *σ*s the evolutionary parameter *η* is compatible with zero (namely, *η*/*ση*). **Lower half.** Similarly to the upper half, the lower half shows the fit parameters of *H*0(*z*) (flat *w*0*wa*CDM model, varying *H*<sup>0</sup> and *wa*) without and with the BAOs.



**Figure 1. Left panel** The *H*0(*z*) vs. z wit varying also Ω0*M*. The red color indicates the case with only SNe Ia as probes, while the blue refers to the case of SNe + 1 BAO per bin. This color-coded will be applied also in the right panel. **Right.** The same plot for the *w*0*wa*CDM model, considering the local fiducial value *H*<sup>0</sup> = 70, where both *H*<sup>0</sup> and *w<sup>a</sup>* are left free to vary.

#### **5. Perspective of the Future Contribution of GRB-Cosmology in 2030**

The discussion of GRBs as possible cosmological tools has been going on for more than two decades [350,351]. The best bet is yet to come since we need first to identify the tightest correlation possible with a solid physical grounding. Among the many correlations proposed [19–23] we here choose to apply the *fundamental plane* (or Dainotti relation) [30,352–354], namely the three-dimensional relation between the end of the plateau emission's luminosity, *La*, its time in the rest-frame, *T* ∗ *a* , and the peak luminosity of the GRB, *Lpeak*: it is possible to estimate how many GRBs are needed to obtain constraints for the cosmological parameters that are comparable with the ones obtained from the other probes, such as SNe Ia and BAOs. After a selection of the best fundamental plane sample through the trimming of GRBs, a simulation of a sample of 1500 and 2000 GRBs according to the properties of the fundamental plane relation has been performed. The fundamental plane relation can be expressed as the following:

$$
\log\_{10} L\_a = a \times \log\_{10} T\_a^\* + b \times \log\_{10} L\_{peak} + c,\tag{10}
$$

where *a*, *b* are the parameters of the plane and *c* is the normalization constant. It is important to stress that here the variables *La*, *T* ∗ *a* , and *Lpeak* have been corrected for evolutionary effects with redshift applying the Efron and Petrosian method [355]. Based on Equation (10), we perform the maximization of the following log-likelihood for the simulated sample of GRBs:

$$\begin{split} \ln \mathcal{L}\_{GRB} &= -\frac{1}{2} \Big( \ln(\sigma^2 + (a \ast \delta\_{\log 10} T\_a^\*)^2 + (b \ast \delta\_{\log 10} L\_{\text{pank}})^2 + \delta\_{\log\_{10} L\_{\text{u}}}^2) \Big) \\ &- \frac{1}{2} \Big( \frac{(\log 10 L\_{a, \text{th}} - \log 10 L\_{\text{a}})^2}{\sigma^2 + (a \ast \delta\_{\log 10} T\_a^\*)^2 + (b \ast \delta\_{\log 10} L\_{\text{pank}})^2 + \delta\_{\log\_{10} L\_{\text{u}}}^2} \Big) . \end{split} \tag{11}$$

where *La*,*th* is the theoretical luminosity computed through the fundamental plane in Equation (10), *σ* is the intrinsic scatter of the plane and *δlog*10*<sup>T</sup>* ∗ *a* , *δlog*10*Lpeak* , and *δlog*10*L<sup>a</sup>* are the errors on the rest-frame time at the end of the plateau emission, the peak luminosity and the luminosity at the end of the plateau, respectively.

After performing an MCMC analysis using the D'Agostini method [356] and letting vary the parameters *a*, *b*, *c*, *σ*, Ω0*m*, the results are shown in Figure 2. Through the simulations of 1000 GRBs, with 9500 steps and keeping the same errors (errors undivided) as the ones observed in the fundamental plane (*n* = 1, see the upper left panel of Figure 2) we obtain a value of Ω0*<sup>m</sup>* = 0.310 with a symmetrized uncertainty of *σ*Ω0*<sup>m</sup>* = 0.078. In the case of 2000 GRBs with 13,000 steps and *n* = 1 (see the upper right panel of Figure 2) instead, we have Ω0*<sup>m</sup>* = 0.300, *σ*Ω0*<sup>m</sup>* = 0.052. If we consider the division of the errors on the variables of the fundamental plane by a factor 2 (halved errors, *n* = 2) we obtain, in the case of 1500 GRBs with 11,100 steps, Ω0*<sup>m</sup>* = 0.300, *σ*Ω0*<sup>m</sup>* = 0.037 (see the lower-left panel of Figure 2), while through 2000 simulated GRBs in 14,600 steps (still with *n* = 2, see the lower right panel of Figure 2) we have Ω0*<sup>m</sup>* = 0.310, *σ*Ω0*<sup>m</sup>* = 0.034. The idea of considering halved errors comes from the prospects for improvement in the fitting procedures of GRB light curves. Through this approach, the GRBs have provided constraints on the value of Ω0*<sup>m</sup>* that are compatible with the ones of previous samples of SNe Ia: in the *n* = 1 cases, the values of the uncertainties are comparable with the ones from [357], while for the *n* = 2 cases the values are close to the ones found in [358] with 2000 GRBs. Furthermore, the GRBs have proven to be promising standardizable candles and, given the bigger redshift span they can cover if compared with SNe Ia, GRBs will provide more complete information about the structure and the evolution of the early universe after the Big Bang, together with quasars [359,360]. After discussing the potentiality of GRBs as future standard candles, we estimate the frequency of GRBs with a plateau emission over the total number of GRBs observed to date. We can expect that by 2030 we will have reached several GRB observations such that these—as standalone probes that respect the properties of the GRB platinum sample [35]—will give constraints as precise as the ones from [357] in the case of not halved errors. In case of halved errors, we can reach the level of precision of [357] even now. In addition, if we consider a machine learning analysis [361,362] for which we can double the size of the sample we are able to reach the precision of [357] now with the case of *n* = 1. If we consider the case of reconstructing the light curves and thus we have a sample which has the 47% of cases with halved errors we can reach the limit of [357] in 2022 if *n* = 1 and now if *n* = 2.

**Figure 2. Upper left.** An example of 1000 simulated GRBs with the posterior distribution of the fundamental plane parameters *a*, *b*, *c*, and its intrinsic scatter *σ* together with the total matter density parameter Ω0*m*. In this case, the steps of the simulation are 9500 and the errors on the variables of the fundamental plane have not been divided by any factor (*n* = 1). **Upper right.** The same case of the upper left panel, but considering 2000 GRBs and a number of steps of 13,000. **Lower left.** The results of 1500 simulated GRBs, dividing by two the errors on the fundamental plane variables (halved errors, *n* = 2): here the steps are 11,100. **Lower right.** The same result of the lower-left panel, considering 2000 GRBs and 14,600 steps instead.

m

## **6. Discussions on the Results**

m

Our results can be interpreted because of astrophysical selection biases or theoretical models alternatives to the standard cosmological models.

## *6.1. Astrophysical Effects*

The main effect that has a stake in the SNe Ia luminosity variation is the presence of metallicity and the difference in stellar ages. Indeed, ref. [284] correct the PS with a mass-step contribution (∆*M*). Despite this term improving the results, other effects need to be accounted for. Considering the stretch and the color, ref. [363] claim that the Hubble residuals, after being properly corrected according to the stretch and color observations, for SNe Ia in low mass and high mass host galaxies show a difference of 0.077 ± 0.014 mag, compatibly with the result of [284]. SNe Ia age metallicity and age are believed to be responsible for the observed behavior: those can replicate the Hubble residual trends consistent with the ones of [363]. In the PS, to account for the evolutions of stretch (*α*) and color (*β*), the parametrization utilized is the following: *α*(*z*) = *α*<sup>0</sup> + (*α*<sup>1</sup> × *z*), *β*(*z*) = *β*<sup>0</sup> + (*β*<sup>1</sup> × *z*). According to [284], there is no clear dependence on the redshift for *α*(*z*) and *β*(*z*), thus *α*<sup>1</sup> and *β*<sup>1</sup> are set to zero. Only the selection effect for color is noteworthy and [284] consider the uncertainty on *β*<sup>1</sup> as a statistical contribution. Concerning the stretch evolution in the PS calibration, it appears to be negligible and is not included at any level.

Conversely, ref. [364] recently studied the SALT2.4 lightcurve stretch and showed that the SN stretch parameter is redshift-dependent. According to their analysis, the asymmetric Gaussian model assumed by [284] for describing the populations of SNe Ia does not take into consideration the redshift drift of the PS, thus leaving a residual evolutionary trend that manifests at higher redshifts. Indeed, the simulations performed by [284] for studying the systematics calibration reach redshifts up to *z* = 0.7: this threshold is present in the third bin of our analysis. The effect from 0.7 ≤ *z* ≤ 2.26 needs additional investigations. It is worth noting that this decreasing trend of *H*<sup>0</sup> (with a given value of *η*) found in [51] is consistent in 1 *σ* for the ΛCDM both in the cases of SNe Ia only and SNe+BAOs. When we consider the *w*0*waCDM*, the *η* values are compatible in 3 *σ* with the ones with SNe Ia only and SNe + BAOs. We here have two cosmological parameters varying at the same time, differently from [51]. Therefore, one of the possible astrophysical reasons behind the observed trend is the residual stretch evolution with redshift. If so, in our work the effect is simply switched from stretch to *H*0. The forthcoming release of the Pantheon+ data [107,365–369] will give the chance to test if these evolutionary effects may be still visible, but this analysis goes far beyond the scope of the current paper. The astrophysical interpretation seems to be favored, but also many theoretical explanations may be possible to describe the outcome of these results.
