On the evolution of the Hubble constant with the SNe Ia Pantheon Sample and Baryon Acoustic Oscillations: a feasibility study for GRB-cosmology in 2030

The difference from 4 to 6 $\sigma$ in the Hubble constant ($H_0$) between the values observed with the local (Cepheids and Supernovae Ia, SNe Ia) and the high-z probes (CMB obtained by the Planck data) still challenges the astrophysics and cosmology community. Previous analysis has shown that there is an evolution in the Hubble constant that scales as $f(z)=\mathcal{H}_{0}/(1+z)^{\eta}$, where $\mathcal{H}_0$ is $H_{0}(z=0)$ and $\eta$ is the evolutionary parameter. Here, we investigate if this evolution still holds by using the SNe Ia gathered in the Pantheon sample and the BAOs. We assume $H_{0}=70 \textrm{km s}^{-1}\textrm{Mpc}^{-1}$ as the local value and divide the Pantheon into 3 bins ordered in increasing values of redshift. Similar to our previous analysis but varying two cosmological parameters contemporaneously ($H_0$, $\Omega_{0m}$ in the $\Lambda$CDM model and $H_0$, $w_a$ in the $w_{0}w_{a}$CDM model), for each bin we implement a MCMC analysis obtaining the value of $H_0$ $[...]$. Subsequently, the values of $H_0$ are fitted with the model $f(z)$. Our results show that a decreasing trend with $\eta \sim10^{-2}$ is still visible in this sample. The $\eta$ coefficient reaches zero in 2.0 $\sigma$ for the $\Lambda$CDM model up to 5.8 $\sigma$ for $w_{0}w_{a}$CDM model. This trend, if not due to statistical fluctuations, could be explained through a hidden astrophysical bias, such as the effect of stretch evolution, or it requires new theoretical models, a possible proposition is the modified gravity theories, $f(R)$ $[...]$. This work is also a preparatory to understand how the combined probes still show an evolution of the $H_0$ by redshift and what is the current status of simulations on GRB cosmology to obtain the uncertainties on the $\Omega_{0m}$ comparable with the ones achieved through SNe Ia.


Introduction
The ΛCDM model is one of the most accredited models, which implies an accelerated expansion phase [1,2]. Although it represents the favored paradigm, it is affected by great challenges: the fine-tuning, the coincidence [3,4], and the Dark Energy nature's problems.
More importantly, the H 0 tension represents a big challenge for modern cosmology. Indeed, the 4.4 up to 6.2 σ discrepancy, depending on the sample used [5][6][7], between the local value of H 0 obtained with Cepheids observations and SNe Ia, H 0 = 74.03 ± 1.42 km s −1 Mpc −1 [8,9], and the Planck data of Cosmic Microwave background radiation (CMB), H 0 = 67.4 ± 0.5 km s −1 Mpc −1 from the Planck Collaboration [10] requires further investigation. From now on, H 0 will be in the units km s −1 Mpc −1 .
We stress that other probes report values of H 0 ≈ 72 ± 2, similar to the value obtained with the SNe Ia. Surely, to solve the Hubble tension it is necessary to use probes that are standard candles. SNe Ia, considered one of the best standard candles, are observed only up to a low redshift range: the farthest so far discovered is at z = 2.26 [11].
It is important for studying the evolution of the cosmological parameters to investigate probes at high redshift. One of the best candidates in this regard is represented by the Gammaray Bursts (GRBs).
Indeed, it is important to stress that once we have established if the Hubble constant undergoes redshift evolution, the Pantheon sample can safely be combined with other probes. Surely the advantage of the use of the SNe Ia is that their emission mechanism is pretty clear, namely they originate from the thermonuclear explosion of carbon-oxygen white dwarfs (C/O WDs).
For GRBs, more investigation about their progenitor mechanism is needed. We here stress that this work can be also preparatory to the work of future application of GRBs as cosmological tools together with SNe Ia and Baryon Acoustic Oscillations (BAOs) through well-established correlations among the prompt variables, such as: the Amati relation [19], which connects the peak in the νF ν spectrum to the isotropic energy in the prompt emission (E iso ), the Yonetoku relation [20,21] between E peak and the peak luminosity of the prompt emission, L peak , the Liang and Zhang relation [22] between E iso , the rest-frame break time of the GRB t b and the peak energy spectrum in the rest frame E p , the Ghirlanda relation (E peak − E jet = E iso × (1 − cosθ)) [23], and the prompt-afterglow relations for the GRBs with the plateau emission investigated in [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38], which have as common emission mechanism most likely the magnetar model, where a neutron star with an intense magnetic field undergoes a fast-spinning down [39][40][41][42][43].
A feasibility study shows that GRBs can give relevant constraints on the cosmological parameters [17,44]. We here give a list of examples of other probes used for measuring the Hubble constant tension. One of them is the use of data from time-delay measurements and Ref. [196] test seven cosmological models through the constraints of SNe Ia, BAO, CMB, Planck lensing, and Cosmic Chronometers with the outcome that in the ΛCDM scenario a flat universe is favored.
Ref. [197] discuss how the new physical scenarios before the recombination epoch imply the shift of cosmological parameters and how these shifts are related to the discrepancy between the local and non-local values of H 0 .
Ref. [198] proposes that the H 0 tension may be solved if the speed of light is treated as a function of the scale factor (as in [199]), and applies this scenario to SNe Ia data.
Ref. [200,201] discuss the implementation of the alternative Phenomenologically Emergent Dark Energy model (PEDE), which can be also extended to a Generalized Emergent Dark Energy model (GEDE) with the addition of an extra free parameter. This shows the possibility of obtaining the PEDE or the ΛCDM cosmology as sub-cases of the GEDE scenario.
Ref. [202] consider a scenario of modified gravity predicting the increase of the expansion rate in the late-universe, thus proving that in this scenario the Hubble tension reduces significantly.
Ref. [203] study the ΛCDM model constrained, at the early-time universe, by the presence of the early Integrated Sachs-Wolfe (eISW) effect, proving that the early-time models aimed at attenuating the Hubble tension should be able to reproduce the same eISW effect just like the ΛCDM does. The observations of a locally higher value for H 0 led to the discussion of local measurements, constraints, and modeling. In this regards, the assumption of a local void [204][205][206] may produce locally an increased value for H 0 .
The Universe appears locally inhomogeneous below a scale of roughly 100 Mpc. The question some cosmologists are attempting to solve is whether local inhomogeneities have impact on cosmological measurements and the Hubble diagram. Many observables are related to photons paths, which may be directly affected by the matter distribution. Many theoretical attempts were made during the last few decades to develop the necessary average prescription to evaluate the photon propagation on the observer's past light cone based on covariant and gauge-invariant observables [207][208][209][210]. Local inhomogeneities and cosmic structure cause scattering and bias effects in the Hubble diagram, which are due to peculiar velocities, selection effects, and gravitational lensing, but also to non-linear relativistic corrections [210][211][212]. This question was addressed in [213] utilising the N-body simulation of cosmic structure formation through the numerical code gevolution. This non-perturbative approach pointed out discrepancies in the luminosity distance between a homogeneous and inhomogeneous scenario, showing, in particular, the presence of non-Gaussian effects at higher redshifts. These studies related to distance indicators will become even more significant considering the large number of the forthcoming surveys designed to the observations on the Large Scale Structure of the Universe in the next decade (for instance, the Euclid survey [214,215] and the Vera C. Rubin Observatory's LSST [216]). The effect of local structures in an inhomogeneous universe should be considered in the locally measured value of H 0 [217,218]. The local under-density interpretation was also studied in Milgromian dynamics [131,219,220], but in [221] it is shown how this interpretation does not solve the tension. Ref. [219] study the KBC local void which is in contrast with the ΛCDM, thus proposing the Milgromian dynamics as an alternative to standard cosmology. Milgromian dynamics are studied also in [222] where, through the galactic structures and clusters, it is shown how this model can be consistent at different scales and alleviate the Hubble tension.
Ref. [223] describe the late time approaches and their effect on the Hubble parameter. The bulk viscosity of the universe is also considered the link between the early and late universe values of H 0 [224].
Ref. [225] explain how the local measurements over-constrain the cosmological models and propose the graphical analysis of the impact that these constraints have on the H 0 estimation through ad hoc triangular plots.
Ref. [220,226,227] describes the effects of inhomogeneities at small scales in the baryon density. Ref. [228] find out that the late time modifications can solve the tension between the H 0 SH0ES and CMB values through a parametrization of the comoving distance.
Ref. [229] propose to alleviate the Hubble tension considering an abrupt modification of the effective gravitational constant at redshift z ≈ 0.01. Other proposals are focused on the existence of different approaches.
Ref. [230,231] show how H 0 evolves with redshift at local scales.
Ref. [232] discuss how the breakdown of Friedmann-Lemaitre-Robertson-Walker (FLRW [233]) may be a plausible assumption to alleviate the Hubble tension. Ref. [234] investigate the binary neutron stars mergers and, with the analysis of simulated catalogs, show their potential to help to alleviate the H 0 tension.
Ref. [235] explain how Gaussian process (GP) and locally weighted scatter plot smoothing are used in conjunction with simulation and extrapolation (LOESS-Simex) methods can reproduce different sets of data with a high level of precision, thus giving new perspectives on the Hubble tension through the simulation of Cosmic Chronometers, SNe Ia, and BAOs data sets.
Ref. [236] focus on the GP and state the necessity of lower and upper bounds on the hyperparameters to obtain a reliable estimation of H 0 .
On the other hand, Ref. [237] suggested a novel approach to measure H 0 based on the distance duality relation, namely a method that connects the luminosity distance of a source to its angular diameter. In this case, data do not require a calibration phase and the relative constraints are not dependent on the underlying cosmological model.
Ref. [238] showed how the tension can be solved with a modified weak-field General Relativity theory, thus defining a local H 0 and a global H 0 value.
Ref. [239] investigated how a specific Dark Energy model in the generalized Proca theory can alleviate the tension.
In [240], the Horndeski model can describe with significantly good precision the late expansion of the universe thanks to the Hubble parameter data. The same model is considered promising for the solution of the H 0 tension in [241].
Ref. [242] described how the transition observed in Tully-Fisher data could imply an evolving gravitational strength and explain the tension.
Ref. [197] explain how the physical models of the pre-recombination era could cause the observed H 0 values discrepancy and suggest that if the local H 0 measurements are consistent then a scale-invariant Harrison-Zeldovich spectrum should be considered to solve the H 0 issue. The Dynamical Dark Energy (DDE) models are the object of study in [243,244]: in the former, the DDE is proposed as an alternative to ΛCDM, while in the latter it is shown how the Chevallier-Polarski-Linder (CPL) parametrization [245,246] is insensitive to Dark Energy at low redshift scales.
Ref. [247,248] propose Dark Radiation as a new surrogate of the Standard Model.
In [249], the scalar field cosmological model is used, together with the parametrization of the equation of state, to obtain H 0 and investigate the nature of Dark Energy. The possibility of a scalar field non minimally coupled to gravity as a probable solution to the H 0 tension is investigated in [250].
Ref. [251] highlight the advantage of the braneworld models to predict the local higher values of H 0 and, contemporaneously, respect the CMB constraints. Another approach is to solve the H 0 tension by allowing variations in the fundamental constants [252].
Ref. [253] propose a non-singular Einstein-Cartan cosmological model with a simple parametrization of spacetime torsion to alleviate the tension, while [254] propose a model where the Dark Matter is annihilated to produce Dark Radiation.
Ref. [255] introduce a hidden sector of atomic Dark Matter in a realistic model that avoids the fine-tuning problem. The observed weak effect of primordial magnetic fields can create clustering at small scales for baryons and this could explain the H 0 tension [256].
Ref. [257] test the General Relativity at galactic scales through Strong Gravitational Lensing. The Strong Lensing is a promising probe for obtaining new constraints on H 0 , thanks to the next generation DECIGO and B-DECIGO space interferometers [258].
In [259], the cosmological constant Λ is considered a dynamical quantity in the context of the running vacuum models and this assumption could tackle the H 0 tension. [260] show the singlet Majoron model to explain the acceleration of the expansion at later times and prove that this is consistent with large-scale data: this model has been subsequently discussed in other works [261]. The vacuum energy density value is affected by the Hubble tension as well and its measurement may cast more light on this topic [262].
Ref. [263] discuss the outcomes of the Oscillatory Tracker Model with an H 0 value that agrees with the CMB measurements. In [264], it is explained how the Generalized Uncertainty Principle and the Extended Uncertainty Principle can modify the Hubble parameter. [265] explore the implication of the Mirror Twin Higgs model and the need for future measurements to alleviate the tension.
The artificial neural networks can be applied to reconstruct the behavior of large scale structure cosmological parameters [266].
Another alternative is given by the gravitational transitions at low redshift which can solve the H 0 tension better than the late-time H(z) smooth deformations [84,229].
Another comparison between the late-time gravitational transition models and other models which predict a smooth deformation of the Hubble parameter can be found in [267].
Ref. [268], as modifications to the ΛCDM model, consider as plausible scenarios or a Dark Matter component with negative pressure or the decay of Dark Energy into Dark Matter.
Ref. [269] does not observe the H 0 tension through the Effective Field Theory of Large Scale Structure and the Baryon Oscillation Spectroscopic Survey (BOSS) Correlation Function.
Considering the Dark Matter particles with two new charges, Ref. [270] reproduce a repulsive force which has similar effects to the Λ cosmological constant. Furthermore, the models where interaction between Dark Matter and Dark Energy is present are promising for a solution of the Hubble constant tensions, see [271].
In [272], it is shown how two independent sets of cosmological parameters, the background (geometrical) and the matter density (growth) component parameters, respectively, give consistent results and how the preference for high values of H 0 is less significant in their analysis.
Ref. [273] introduce a global parametrization based on the cosmic age which rules out the early-time and the late-time frameworks.
Ref. [274] point out, through the use of non-parametric methods, how the cosmological models may induce biases in the cosmological parameters. In the same way, the statistical analysis of galaxies' redshift value and distance estimations may be affected by biases which could, in turn, affect the estimation of H 0 Ref. [275]. This consideration holds also for the quadruply lensed quasars which are another method to measure H 0 [276].
Ref. [277] use the machine learning techniques to measure time delays in lensed SNe Ia, these being an independent method to measure H 0 .
Additionally, in [231] it is explained how an evolution of H 0 with the redshift is to be expected. If a statistical approach on the different H 0 values is used instead, together with the assumption of an alternative cosmology, another solution to the tension could be naturally implied [278].
Ref. [279] use data to reconstruct the f (T) gravity function without assuming any cosmological model: this f (T) could in turn represent a solution to the H 0 tension.
Ref. [124] discuss how the addition of scalar fields with particle physics motivation to the cosmological model which predicts Dark Matter can retrieve the observed abundances of the Big Bang Nucleosynthesis.
In [280], a Dark Matter production mechanism is proposed to alleviate the H 0 tension. A general review of the perspectives and proposals concerning the H 0 tension can be found in [281] and [282,283].
SNe Ia represents a very good example of standard candles. Here we consider also the contribution of geometrical probes, the so-called standard rulers: while standard candles show a constant intrinsic luminosity (or obey an intrinsic relation between their luminosity and other physical parameters independent of luminosity), standard rulers are characterized by a typical scale dimension. This property allows estimating their distance according to the apparent angular size. Among the possible standard rulers, the BAOs assume great importance for cosmological purposes.
We here investigate the H 0 tension in the Pantheon sample (hereafter PS) from [284] and we add the contribution of BAOs to the cosmological computations to check if the trend of H 0 found in [36] is present also with the addition of other probes. We here point out that the current analysis is not meant to constrain Ω 0m or any other cosmological parameters, but it is focused to study the reliability of the trend of H 0 as a function of the redshift.
We here point out that this analysis is not meant to constrain Ω 0m or any other cosmological parameters, but it is focused to study the reliability of the trend of H 0 as a function of the redshift. The range of redshift in the PS goes from z = 0.01 to z = 2.26. We tackle the problem with a redshift binning approach of H 0 , the same used in [51], but here we adopt a starting value of H 0 = 70 instead of 73.5: if a trend with redshift exists, it should be independent on the initial value for H 0 . The systematic contributions for the PS are calibrated through a reference cosmological model, where H 0 is 70.0 [284]. In the current paper, the aforementioned systematic uncertainties are considered for the analysis. Our approach has a two-fold advantage: on the one hand, it is relatively simple and on the other hand, it avoids the re-estimation of the SNe Ia uncertainties and may be able to highlight a residual dependence on the SNe Ia parameters with redshift.
While a slow varying Einstein constant with the redshift, as it emerges in a modified f (R) gravity, appears as the most natural explanation for a trend H 0 (z), the analysis of Section 7 seems to indicate that such effect is not necessarily related with the Dark Energy contribution of the late universe. Since the Hu-Sawicki gravity lacks of reproducing the correct profile H 0 (z) shows that a Dark Energy model in the late Universe may not be enough to explain the observed effect since the scalar mode dynamics can not easily conciliate the Dark Energy contribution with the decreasing trend of H 0 (z). Thus, it may be necessary a modified gravity scenario more general than a Dark Energy model in the late Universe.
The current paper is composed as expressed in the following: in Section 2 the ΛCDM and w 0 w a CDM models are briefly introduced together with SNe Ia properties; Section 3 describes the use of BAOs as cosmological rulers; Section 4 contains our binned analysis results, after slicing the PS in 3 redshift bins for the aforementioned models, and assuming locally H 0 = 70; in Section 5, we investigate, through simulated events, how the GRBs will be contributing to cosmological investigations by 2030; in Section 6 we discuss the results; in Section 7 we test the Hu-Sawicki model through a binning approach; in Section 8 we report an overview on the requirements that a suitable f (R) model should have to properly describe the observed trend of H 0 and in Section 9 our conclusions are reported.

SNe Ia Cosmology
SNe Ia are characterized by an intrinsic luminosity that is almost uniform. Because of this, SNe Ia are considered reliable standard candles. We compare the theoretical distance moduli µ th with the observed distance moduli µ obs of SNe Ia belonging to the PS. The theoretical distance moduli are defined through the luminosity distance d L (z) which we need to define based on the cosmological model of interest. We here show the CPL parametrization which describes the w parameter as a function of redshift (w(z) = w 0 + w a × z/(1 + z)) in the w 0 w a CDM model. In the usual assumptions w 0 ∼ −1 and w a ∼ 0, and d L (z) is defined as the following [285]: where Ω 0Λ is the Dark Energy component, c is the speed of light, and z is the redshift. We stress that in this context the relativistic components are ignored. Moreover, since in the present universe the radiation density parameter Ω 0r ≈ 10 −5 , this contribution can be neglected.
If we substitute w a = 0, w 0 = −1 in Equation (1) the luminosity distance expression for ΛCDM model is automatically retrieved. According to the distance luminosity expression, the theoretical distance modulus can be written in the following form: which is usually expressed in Megaparsec (Mpc). The observed distance modulus, µ obs = m B − M, taken from PS contains the apparent magnitude in the B-band corrected for statistical and systematic effects (m B ) and the absolute in the B-band for a fiducial SN Ia with a null value of stretch and color corrections (M). Considering the color and stretch population models for SNe Ia, in our approach we average the distance moduli given by the [286] (G2010) and [287] (C2011) models. We here remind the reader that H 0 and M are degenerate parameters: in the PS release, M = −19.35 such that H 0 = 70.0. Ref. [51] obtain information on H 0 by comparing µ obs in [284] 1 with µ th for each SN. Moreover, they fix Ω 0m to a fiducial value to better constrain the H 0 parameter. Furthermore, according to [288], we consider the correction of the luminosity distance keeping into account the peculiar velocities of the host galaxies which contain the SNe Ia. To perform our analysis, we define the χ 2 for SNe: Here ∆µ = µ obs − µ th , and C denotes the 1048 × 1048 covariance matrix, given by [284]. As for the µ obs values of G2010 and C2011, the systematic uncertainty matrices of the two models have been averaged. After building the C total matrix from Equation (16) in [51], we slice the PS in redshift bins, and then we divide C into submatrices considering the order in redshift. More in detail, starting from the 1048 SNe Ia redshift-ordering, we divide the SNe Ia into 3 equally populated bins made up of ≈349 SNe Ia. Concerning only D stat , it is trivial to build its submatrices considering that the statistical matrix is diagonal. Hence, a single matrix element is related to a given SN of the PS. On the other hand, if the non-diagonal matrix C sys is included, a customized code will be used 2 to build the submatrices. Our code was developed to select only the total covariance matrix elements related to SNe Ia having redshift within the considered bin.
The choice of three bins is justified by the high number of SNe Ia (around hundreds of SNe per bin) that can still constitute statistically illustrative subsamples of the PS and that can properly consider the contribution of systematic uncertainties. Subsequently to the bins division, we focus on the optimal values of H 0 to minimize the χ 2 in Equation (3). H 0 is regarded as a nuisance parameter, which is free to vary, to better analyze a possible redshift function of H 0 . We follow the assumptions on the fiducial value of M = −19.35: while in [51] M was estimated assuming a local (z = 0) value of H 0 = 73.5, we here consider the conventional H 0 value of the PS release, namely H 0 = 70.0 for three bins. Our choice of a starting value of H 0 = 70 is dictated by the presence in the current literature of more than 50 papers that are using the PS in combination with other probes to estimate the value of H 0 , see [172,198,204,205,[288][289][290][291][292][293][294][295][296][297][298][299][300][301][302][303][304][305][306][307]. Thus, if an evolutionary effect is present, it is necessary to investigate to which extent this can affect current and future results largely based on the PS sample. Conversely, we fix Ω 0m = 0.298 ± 0.022 according to [284] for a standard flat ΛCDM model. More specifically, after the minimization of χ 2 , we extract the H 0 value in each redshift bin, via the Cobaya code [340]. To this end, we execute an MCMC using the D'Agostini method to obtain the confidence intervals for H 0 at the 68% and 95% levels, in three bins.

The Contribution of BAOs
The environment of relativistic plasma in the early universe was crossed by the sound waves that were generated by cosmological perturbations. At redshift z d ∼ 1059.3, which marks the ending of the drag period [341], the recombination of electrons and protons into a neutral gas interrupted the propagation of the sound waves while the photons were able to propagate further [342]. In the period between the formation of the perturbations and the recombination, the different modes produced a sequence of peaks and minima in the anisotropy power spectrum. Given the huge fraction of baryons in the universe, it is expected by cosmological models that the oscillations may affect also the distribution of baryons in the late universe. As a consequence, the BAOs manifest as a local maximum in the correlation function of the galaxies distribution in correspondence of the comoving sound horizon scale at the given redshift z d , namely r s (z d ): this is associated with the stopping of the propagation of the acoustic waves.
To use the BAOs data for cosmology, we first need to define the following variables: The value of the redshift z d , which corresponds to the drag era ending and marks the decoupling of the photons, allows estimating the sound horizon scale: where we use the adimensional ratio h = H 0 /100(km s −1 Mpc −1 ). To estimate r s , the following approximated formula [343] can be applied: where ω i = Ω i · h 2 , and i = m, ν, b represent matter, neutrino and baryons. We here assume ω ν = 0.00064 [344] and ω b = 0.02237 [10]. Given these quantities, we define the χ 2 for BAOs as follows: where ∆d = d obs z (z i ) − d theo z (z i ) and M is the covariance matrix for the BAO d obs z (z i ) values. In this binned analysis, a subset of the 26 BAO observations set available in [341] will be employed.

Multidimensional Binned Analysis with SNe Ia and BAOs
To investigate the H 0 tension through the SNe Ia and BAOs data, we combine the χ 2 Equations (3) and (7) to obtain the total χ 2 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 w a 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 0 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 0 and Ω 0m vary simultaneously, while in the w 0 w a CDM model the varying parameters are H 0 and w a . We decided to leave w a free to vary since, according to the CPL parametrization, w a gives direct information about the evolution of the w(z) while w 0 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 Ω 0m = 0.298 ± 0.022 and H 0 = 70.393 ± 1.079 for Ω 0m and H 0 , respectively, and with 1 σ = 2 * 0.022 and 1 σ = 2 * 1.079 for Ω 0m 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 w a CDM model, we fix w 0 = −0.905 and we consider the priors on w a 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 χ 2 minimization for each bin, we perform a MCMC simulation to draw the mean value of H 0 and its uncertainty. Once H 0 is obtained for each bin, we perform a fit of H 0 using a simple function largely employed to characterize the evolution of many astrophysical objects, such as GRBs and quasars [17,29,31,35,[345][346][347][348][349]. More specifically, the fitting of H 0 is given by in which H 0 and η are the fitting parameters. The former H 0 ≡ H 0 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 0 = 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 η/σ η ratio for three bins down to 1.2. Concerning the lower half of Table 1 with the w 0 w a 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 η/σ η is observed when BAOs are added in the case of w 0 w a CDM model in three bins. The results can be visualized in Figure 1. Comparing the η/σ η 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 w a CDM model are compatible in 3 σ with the α values in the same reference paper.

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][20][21][22][23] we here choose to apply the fundamental plane (or Dainotti relation) [30,[352][353][354], namely the three-dimensional relation between the end of the plateau emission's luminosity, L a , its time in the rest-frame, T * a , and the peak luminosity of the GRB, L peak : 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: where a, b are the parameters of the plane and c is the normalization constant. It is important to stress that here the variables L a , T * a , and L peak 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: (log 10 L a,th − log 10 L a ) 2 where L a,th is the theoretical luminosity computed through the fundamental plane in Equation (10), σ is the intrinsic scatter of the plane and δ log 10 T * a , δ log 10 L peak , and δ log 10 L a 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, σ, Ω 0m , 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 Ω 0m = 0.310 with a symmetrized uncertainty of σ Ω 0m = 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 Ω 0m = 0.300, σ Ω 0m = 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 11100 steps, Ω 0m = 0.300, σ Ω 0m = 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 Ω 0m = 0.310, σ Ω 0m = 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 Ω 0m 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.

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

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) = α 0 + (α 1 × z), β(z) = β 0 + (β 1 × z). According to [284], there is no clear dependence on the redshift for α(z) and β(z), thus α 1 and β 1 are set to zero. Only the selection effect for color is noteworthy and [284] consider the uncertainty on β 1 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 0 (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 w a CDM, 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][366][367][368][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.

Theoretical Interpretations
We now investigate possible theoretical explanations for our results, focusing particular attention on modified gravity models. We first discuss a general scalar-tensor formulation and, then, we concentrate our attention on the so-called metric f (R) gravity.

The Scalar Tensor Theory of Gravity
The action of the scalar tensor theories (STTs) of gravity is given by S = S JF + S m [370][371][372][373][374] with the Jordan Frame (JF) action whereR is the Ricci scalar obtained with the physical metricg µν , while the matter fields Ψ m couple to the metric tensorg µν and not to Φ, i.e., S m = S m [Ψ m ,g µν ].
In this Section we adopt natural units such that c = 1 and G = 1. Different STTs follow with the appropriate choice of the two functions ω(Φ) andṼ(Φ): e.g., the Brans-Dicke (BD) theory [375][376][377][378] can be obtained for ω(Φ) = ω (const.) andṼ(Φ) = 0, while the metric f (R) gravity, discussed in the next subsection, would correspond to ω ≡ 0. The action S JF can be rewritten in the Einstein Frame (EF), where one definesg µν ≡ A 2 (ϕ)g µν , Matter is coupled to ϕ only through a purely metric coupling, S m = S m [Ψ m , A 2 (ϕ)g µν ] and M * is the Planck mass. The physical quantities in the Jordan and Einstein frame are related by dτ = A(ϕ)dτ,ã = A(ϕ)a,ρ = A(ϕ) −4 ρ,p = A(ϕ) −4 p, where τ is the synchronous time variable. Defining N ≡ log a a 0 , λ ≡ V(ϕ) ρ , w ≡ p ρ , and ϕ = dϕ dN = a dϕ da , the combination of cosmological equations allows to write the equation for ϕ in the form (for a flat Friedmann-Robertson-Walker geometry) [379] 2 3 Moreover, the Jordan-and Einstein-frame Hubble parameters,H ≡ d logã/dτ and H ≡ d log a/dτ, respectively, are related asH For our purpose, we consider A(ϕ) = A 0 e c 1 ϕ+c 2 ϕ 2 /2 , which implies γ(ϕ) = c 1 + c 2 ϕ, where c 1,2 are constants. Under the following conditions ϕ /ϕ 1, ϕ 2 /ϕ 2 1, and the solution of Equation (14) is ϕ(z) = C(1 + z) K − c 1 c 2 , where K = 1−3w 1+w √ 2 c 2 , and C is an integration constant. We are looking for solutions such that H = f (ϕ)H 0 , so thatH = H 0 (1+z) η , whereH 0 is constant. These relations and (15) allow to derive f (ϕ) (the expression of f (ϕ) is quite involved, and in the case in which c 1,2 1, it is a polynomial in ϕ). The scalar field Φ in the (physical) JF can be cast in the form Φ(z) , and z < 1 has been used (note:K is positive for c 1 or c 2 negative). The scalar field Φ reduces to φ for Φ 0 → 1 andK → 2η. From the Friedmann Equation [379] ȧ a with ρ given by matter (ρ = ρ 0m /a 3 = ρ 0m (1 + z) 3 ), and c 1,2 1, one infers the effective potentialṼ where we recall that Ω 0m = ρ 0m /ρ cr , ρ cr = 3M 2 * H 2 0 , f 0 = f (ϕ = 0), and m 2 = Ω 0mH 2 0 . For redshift 0 z < 0.3, to which we are interested, the scalar field varies slowly with z, Φ ∼ Φ 0 , so that the effective potential behaves like a cosmological constant. We see how the proposed scalar-tensor formulation has the right degrees of freedom to reproduce, in the JF, the required behavior of the (physical) trend of H 0 (z). In the next subsection, we analyze a sub-case of the general paradigm discussed above, which leads to the well-known f (R) gravity, which is among the most popular modified gravity formulations.

Metric f(R) Gravity in the Jordan Frame
The observed decaying behavior of the Hubble constant H 0 with the redshift draws significant attention for an explanation and, if it is not due to selection effects or systematics in the sample data, we need to interpret our results from a physical point of view. As already argued in [380] and [51], the simplest way to account for this unexpected behavior of H 0 (z) is that the Einstein constant χ = 8πG (where G denotes the gravitational constant), mediating the gravity-matter interaction, is subjected itself to a slow decaying profile with the redshift. In this Section, we consider c = 1 for the speed of light. More specifically, since the critical energy density ρ c0 = 3 H 2 0 /χ today must be a constant, we need an evolution for χ ∼ (1 + z) −2 η , considering the function H 0 (z) given by Equation (9). The evolution of χ(z) is not expected within the cosmological Einsteinian gravity, therefore we are led to think of it as a pure dynamical effect, associated with a modified Lagrangian for the gravitational field beyond the ΛCDM cosmological model. Ref. [143] obtained cosmological constraints within the Brans-Dicke theory considering how the evolution of the gravitational constant G, contained in χ, affects the SNe Ia peak luminosity. The most natural extended framework is the f (R)-gravity proposal [162,163,167,381] which contains only an additional scalar degree of freedom. For instance, Ref. [323] try to alleviate the H 0 tension considering exponential and power-law f (R) models.
The formulation of the f (R) theories in an equivalent scalar-tensor paradigm turns out to be particularly intriguing for our purposes: the function f (R) is restated as a real scalar field φ, which is non-minimally coupled to the metric in the JF. The information about the function f turns into the expression of the scalar field potential V(φ). The relevance of modified gravity models relies on the possibility that this revised scenario for the gravitational field can account for the physics of the so-called "dark universe"component without the need for a cosmological constant. Indeed, the observed cosmic acceleration in the late universe via the SNe Ia data is a pure dynamical effect, i.e., associated with a modification of the Einsteinian gravity at very large scales (in the order of the present Hubble length).
According to the standard literature on this field (which includes a large number of proposals), three specific f (R) models, i.e., the Hu-Sawicki [382], the Starobinsky [383], and Tsujikawa models [384,385], successfully describe the Dark Energy component (say an effective parameter for the Dark Energy w = w(z) < −1/3) and overcome all local constraints. The difference in the form of the Lagrangian densities associated with f (R) models is reflected in the morphology of the potential term governing the dynamics of the scalar field. For instance, the scalar field potential related to the Hu-Sawicki f (R) proposal, with the power index n = 1, in the JF is given by where we have two free parameters c 1 and c 2 , while m 2 = χ ρ 0m /3. The scalar-tensor dynamics in the JF for a flat FLRW metric with a matter component is summarized by which are the generalized Friedmann equation, the generalized cosmic acceleration equation and the scalar field equation, respectively [167]. We recall that φ = φ(t) is a function of the time (or the redshift z) only for an isotropic universe. Considering the first term on the right-hand side of Equation (19), it is possible to recognize that φ mediates the gravity-matter coupling, and therefore it mimics a space-time varying Einstein constant. Hence, to account for our observed decay of H 0 (z), we have to require that the scalar field assumes a specific behavior with the redshift, i.e.
Moreover, the remaining terms contained in the gravitational field equations must be negligible. This situation is naturally reached when the potential term is sufficiently slow-varying in a given time interval. We see that the hypothesis of a near-frozen scalar field evolution is a possible assumption, as far as the potential term should provide a dynamical impact, sufficiently close to a cosmological constant term. These simple considerations lead us to claim that this scenario is worth to be investigated for the behavior of H 0 (z) here observed. The specific cosmological models affect the expression of the luminosity distance and this should be the starting point of a careful test of a f (R) theory versus the comprehension of the H 0 tension. A new binned analysis of the PS, using the corrected luminosity distance obtained through a reliable f (R), may in principle shed new light on the observed decaying trend of H 0 (z), testing also new physics. This analysis is performed in the next Section.
As a preliminary approach, we try to understand which profile we could expect for the scalar field potential, inferred from the behavior of H 0 (z). This is quite different from a standard analysis of f (R) models. Generally, a specific f (R) function is defined a priori, and then the dynamical equations are studied to obtain constraints on the free parameters. Here, instead, starting from the observed decreasing trend of H 0 (z) and assuming φ(z) from Equation (22), we wonder what the scalar field potential would be in a scalar-tensor dynamics. Eventually, we should have a scalar field in near-frozen dynamics, i.e., a slow-roll of the scalar field potential, mimicking a cosmological constant term (φ → 1). To this end, we rewrite the generalized Friedmann Equation (19) and calculate V(φ): where we have used the standard definition of redshift and the relation (22) for φ(z). Moreover, we recall that for a matter component ρ ∼ (1 + z) 3 . As a final step, we need to calculate the term dz dt . Starting again from the redshift definition, it is well known that In principle, we would need to compute the Hubble parameter H(z) from the field equations, and then replace H(z) in the term dz dt . However, this procedure is not viable, since we need to fix a well-defined V(φ) to solve the field equations. Moreover, H(z) appears also in the right-hand-side of Equation (19), because of the non-minimal coupling with the scalar field. Therefore, we can not calculate exactly dz dt to get V(φ) in the JF. Then, to obtain V(φ) inferred from the trend of H 0 (z), we require that the Hubble function provides the same physical mechanism suggested from our binned analysis in Section 4, i.e., simply replacing H 0 with H 0 (z) given by Equation (9) in the standard Friedmann equation in the ΛCDM model. With this new definition of H 0 , we write the following condition on the Hubble function: In doing so, using Equations (23), (24) and (25), we determine the form of the scalar field potential inferred from the decreasing trend of H 0 (z). In other words, the potential Equ. 26 might provide an effective Hubble constant that evolves with redshift. In the computation, we have used the expression Ω 0m = m 2 /H 2 0 . In Figure 3, we plotted this potential profile, observing that, as expected, a flat region consistently appears, validating our guess on the feasibility of f (R)-gravity in the JF to account for the observed behavior of H 0 (z). We set η = 0.009 in Figure 3, according to our binned analysis results for three bins (see Table 1). We stress that the flatness of the potential does not emerge throughout the Pantheon sample redshift range, 0 < z < 2.3, but it appears only in a narrow region for 0 < z z * , where z * = 0.3 is the redshift at the Dark Energy and Matter components equivalence of the universe. This form of V(φ) is reasonable since the Dark Energy contribution, provided by the scalar field in the JF gravity, dominates the matter component only for 0 < z z * . It is the weak dependence of H 0 on z that ensures the existence of a flat region of the potential, according to the theoretical scenario argued above.
Finally, we can calculate the form of the f (R) function associated with the potential profile. Recalling the following general relations in the JF [167]: we can obtain: Note that the formula above provides a generalization of the Einstein theory of gravity, as it should be in the context of a f (R) model. Indeed, if η = 0, then f (R) ≡ R reproduces exactly the Einstein-Hilbert Lagrangian density in GR with a cosmological constant Λ, as soon as you recognize that Λ = 3m 2 (1 − Ω 0m ) / Ω 0m for a flat geometry, using m 2 = H 2 0 Ω 0m . In particular, expanding the function (29) for η ∼ 0, we can see explicitly the deviation from the Einstein-Hilbert term: The first term at the zero-th order in η is exactly the Einstein-Hilbert Lagrangian density, while the linear term in η provides the correction to GR. Therefore, η, in addition to being the  Figure 3. Profile of the scalar field potential V(φ) in the JF equivalent scalar-tensor formalism of the f(R) modified gravity. The form of V(φ)) is inferred from the behavior of H 0 (z) (Equation (9)). Note that V(φ)/m 2 is a dimensionless quantity. A flat profile of V(φ) occurs only at low redshifts, for 0 < z 0.3 or equivalently φ 1.005. Note, also, the non-linearity of the scale for the redshift axis on top, considering the relation (22) between φ and z. In this plot, η = 0.009.
physical parameter that describes the evolution of H 0 (z), also denotes the deviation from GR and the standard cosmological model. It is worthwhile to remark that the expression above may not be the final form of the underlying modified theory of gravity, associated with the global universe dynamics, but only its asymptotic form in the late Universe, i.e., as the scalar of curvature approaches the value corresponding to the cosmological constant in the ΛCDM model. In all these computations we do not consider relativistic or radiation components at very high redshifts, but it may be interesting to test this model with other local probes in the late Universe.
In this discussion, we infer that the dependence of H 0 on the redshift points out the necessity of new physics in the description of the universe dynamics and that such a new framework may be identified in the modified gravity, related to metric theories.

The Binned Analysis with Modified f(R) Gravity
To try to explain the observed trend of H 0 (z), we focus on f (R) theories of gravity, and then we perform the same binned analysis, using the correction for the distance luminosity according to the modified gravity. We start from the gravitational field action [167]: where f (R), as a function of the Ricci scalar R, is an extra degree of freedom compared to General Relativity. We rewrite f (R) = R + F(R) to highlight the deviation from the standard gravity. Varying the total action with respect to the metric, we obtain the flat FLRW metric field equations: where F R ≡ dF(R) dR . The Ricci scalar R can be cast in the form where the Hubble parameter H is expressed as a function of γ ≡ ln(a), and the prime indicates the derivative with respect to γ. Now, we introduce two dimensionless variables [382] which denote the deviation of H 2 and R with respect to the matter contribution when compared to the ΛCDM model. We rewrite the modified Friedmann Equation (32) and the Ricci scalar relation (33) in terms of y H and y R . Then, we have a set of coupled ordinary differential equations: The solution of this coupled first-order differential equations system above can not be obtained analytically, but can be numerically calculated. We need initial conditions such that this scenario mimics the ΛCDM model in the matter dominated universe at initial redshift z i z * . Hence, we impose the following conditions for y H and y R at the redshift z i : The standard ΛCDM model is reached for z = z i or asymptotically, and we consider a flat geometry, such that Ω 0Λ = 1 − Ω 0m . Finally, the luminosity distance can be written as including the solution y H (z) from Equation (34) [386].

Hu-Sawicki Model
We focus on the Hu-Sawicki model with n = 1, considering a late-time gravity modification, described by the following function [382]: corresponding to the potential V(φ) in Equation (18). The parameters c 1 and c 2 are fixed by the following conditions [382] where F R0 is the value of the field F R ≡ dF/dR at the present time, and F(R) is the deviation from the Einstein-Hilbert Lagrangian density. Cosmological constraints provide |F R0 | ≤ 10 −7 from gravitational lensing and |F R0 | ≤ 10 −3 from Solar system [387,388]. We explore several choices of F R0 .
To simplify the numerical integration of the modified luminosity distance (39), we approximate the numerical solution y H , obtained from the system (35), (36), by a polynomial of order 8. This function is an accurate representation of y H when we restrict the solution to the range of PS (see Figure 4).
As a consequence, we obtain constraints on c 1 and c 2 , according to Equations (41) and (42). Then, we perform the same binned analysis of Section 4 using the Hu-Sawicki model and the modified luminosity distance (39). We here run the analysis both for the case of Ω 0m fixed to a fiducial value of 0.298 and for several values of F R0 = −10 −7 , −10 −6 , −10 −5 , −10 −4 (see Table 2 and Figure 5) or we let Ω 0m vary with the two values of F R0 = −10 −7 , −10 −4 (see Table 3 and Figure 6) for the SNe alone and with SNe +BAOs. Note also that the η parameters are all consistent for the several values of F R0 in 1 σ, as you can see in Table 2, for both SNe Ia and SNe Ia + BAOs. Moreover, the values of η are consistent in 1 σ with the ones obtained from the analysis of the ΛCDM model (see also Table 1). We consider the cases F R0 = −10 −7 , F R0 = −10 −4 and, to study how these results may vary according to the different values of Ω 0m chosen, we tested the model with four values of Ω 0m = (0.301, 0.303, 0.305) taken from the 1 σ from a Gaussian distribution centred around the most probable value of 0.298, see [368].
We show in Figure 6 the comparison between the different applications of the Hu-Sawicki model: in the left panels (upper and lower), we consider SNe Ia only, while in the right panels (upper and lower) we combine SNe Ia+BAOs. We here remind that the assumed values for |F R0 | of 10 −4 and 10 −7 are well constrained by the f (T) theories. [167,382]. Table 2. Fitting parameters of H 0 (z) for three bins within the Hu-Sawicki model, with SNe only and SNe + BAOs with a fixed value of Ω 0m = 0.298 and with several values of F R0 : −10 −4 , −10 −5 , −10 −6 , −10 −7 . The columns contains: (1) the number of bins; (2) H 0 , (3) is η, according to Equation (9); (4) how many σs η is compatible with zero (namely, the ratio η/σ η ); (5) F R0 values; (6) the sample used. Thus, the existence of this trend is, once again, confirmed, and it remains unexplained also in the modified gravity scenario. Indeed, a suitable modified gravity model which would be able to predict the observed trend of H 0 , would allow observing a flat profile of H 0 (z) after a binned analysis. Further analysis must be carried out with other Dark Energy models or other modified gravity theories to investigate this issue in the future, for instance focusing on the proposed model in Section 6.2.2. Table 3. Fitting parameters of H 0 (z) for three bins within the Hu-Sawicki model, with SNe and SNe + BAOs by fixing several values of Ω 0m = 0.298, 0.303, 0.301, 0.305 and values of F R0 = −10 −4 and F R0 = −10 −7 . The columns are as follows: (1) the Ω 0m value; (2) H 0 , (3) η, according to Equation (9); (4) how many σs the evolutionary parameter η is compatible with zero (namely, η/σ η ); (5) F R 0; (6) the sample used.    The same as the lower left, but with the contribution of BAOs. The orange color refers to Ω 0m = 0.301, the red to Ω 0m = 0.303, the magenta to Ω 0m = 0.305, and the blue to Ω 0m = 0.298.

Requirements for a Suitable f(R) Model
Since the Hu-Sawicki model seems to be inadequate to account for the observed phenomenon of the decaying H 0 (z), in what follows, we provide some general properties that an f (R) model in the JF must possess to induce the necessary scenario of a slowly varying Einstein constant. Now, we consider again the dynamical impact of the scalar field φ, related to the f (R) function. Let us observe that the following relation holds in the following way: In order to get the desired behavior φ (1 + z) 2η , we must deal with a dynamical regime where the following request is satisfied:φ We consider a slow-rolling evolution of the scalar field φ in the late universe, near enough to φ 1. Then, we consider in Equation (19) ρ ∼ 0, because we are in the Dark Energy dominated phase, and we consider H 0φ small with respect to the potential term V(φ 1). We neglect, also, the termφ. Under these conditions, Equations (19) and (21) become andφ respectively. Referring to Equation (45) at z ∼ 0, we make the identification H 2 0 ≡ V(φ 1)/6. Hence, in order to reproduce Equation (44), we must require that for φ → 1, the following relation holds: The analysis above states the general features that a f (R) model in the JF has to exhibit to provide a viable candidate to reproduce the observed decay behavior of H 0 (z) (Equ. 36). We conclude by observing that the picture depicted above relies on the concept of a slow-rolling phase of the scalar field, when it approaches the value φ 1 and, in this respect, the potential term should have for such value a limiting dynamics, which remains there confined for a sufficiently long phase. It is just in such a limit that we are reproducing a ΛCDM model, but with the additional feature of a slowly varying Einstein constant. As far as the value of z increases, the deviation of the considered model from General Relativity becomes more important, but this effect is observed mainly in the gravity-matter coupling. In other words, the motion of the photon, as observed in the gravitational lensing, is not directly affected by the considered deviation, since the geodesic trajectories in the space-time do not directly feel the Einstein constant value. This consideration could allow for a large deviation of φ from the unity that is expected in studies of the photons' propagation.

An Example for Low Redshifts
As a viable example for the Dark Energy dominated Universe (slightly different from the traced above), we consider a potential term (and the associated slow-rolling phase) similar to the one adopted in the so-called chaotic inflation [389,390], i.e.,: where δ is a positive constant, such that δ 6H 2 0 . From Equation (45), we immediately get where, we recall that we are considering the slow-rolling phase near φ → 1. Analogously, from Equation (47), we immediately get: The negative value of η is coherent with the behavior H 2 ∝ φ. Hence, we can reproduce the requested behavior of φ(z) by properly fixing the value of δ to get η as it comes out from the data analysis of Section 4. Specifically, we get δ ∼ 10 −3 H 2 0 to have η ∼ 10 −2 .
Furthermore, it is easy to check that, for φ → 1, Equation (44) and Equation (49), we find the relationφ ∼| η | H 0φ 3H 0φ , (51) which ensures that we are dealing with a real slow-rolling phase. Finally, we compute the f (R) function corresponding to the potential in Equation (48), recalling the relation (28): We conclude by observing that this specific model is reliable only as far the universe matter component is negligible, z < 0.3. The Einstein constant in front of the matter-energy density ρ would run as (1 + z) 2η . The example above confirms that the f (R) gravity in the JF is a possible candidate to account for the observed effect of H 0 (z), but the accomplishment of a satisfactory model for the whole ΛCDM phase requires a significant effort in further investigation, especially accounting for the constraints that observations in the local universe provided for modified gravity.

Discussion
Let us now try to summarize the physical insight that we can get from the analysis above, about the possible theoretical nature of the observed H 0 (z) behavior. We can keep as a reliably good starting point the idea that the origin of a modified scaling of the function H(z) with respect to the standard ΛCDM model can be identified in a slowly varying Einstein constant with the redshift. Furthermore, it is a comparably good assumption to search, in the framework of a scalar-tensor formulation of gravity, the natural explanation for such a varying Einstein constant. As shown in Section 6.2, a scalar-tensor formulation can reproduce the required scaling of the function H(z), which we observe as an H 0 (z) behavior in the standard ΛCDM model. Hence, we naturally explored one of the most interesting and well-motivated formulations of a scalar-tensor theory, namely the f (R) gravity in the JF. In this respect, in Section 6.2.2, we first evaluated the form of the scalar field potential inferred from the observed decreasing trend of H 0 (z), and our data analysis suggested a model described in Equations (26) and (29). Then, we investigated if, one of the most reliable models for reproducing the Dark Energy effect with modified gravity, i.e., the Hu-Sawicki proposal, was able to induce the requested luminosity distance to somehow remove the observed effect, thus accounting for its physical nature. The non-positive result of this investigation leads us to explore theoretically the question of reproducing simultaneously the Dark Energy contribution and the observed H 0 (z) effect, by a single f (R) model of gravity in the JF. In Section 8, it has been addressed this theoretical question, by establishing the conditions that a modified gravity model has to satisfy to reach the simultaneous aims mentioned above. Finally, we considered a specific model for the late universe, based on a slow-rolling picture for the scalar field near its today value φ 1. This model was successful in explaining the Dark Energy contribution and the necessary variation of the Einstein constant, but it seems hard to be reconciled with the earlier Universe behavior, when the role of the matter contribution becomes relevant. Thus, based on this systematic analysis, we can conclude that the explanation for H 0 (z) is probably to be attributed to modified gravity dynamics, but it appears more natural to separate its effect from the existence of a Dark Energy contribution. In other words, we are led to believe that what we discovered about the SNe Ia+BAOs binned analysis must be regarded as a modified gravity physics of the scalar-tensor type, but leaving on a standard Universe, well represented by a ΛCDM model a priori.

Conclusions
We analyzed the PS together with the BAOs in three bins in both the ΛCDM and w 0 w a CDM models to investigate if an evolutionary trend of H 0 persists also with the contribution of BAOs and by varying two parameters contemporaneously with H 0 (Ω 0m and w a for the ΛCDM and w 0 w a CDM, respectively). The persistence of the trend of H 0 as a function of redshift is also shown in the case of the Hu-Sawicki model. We here stress that the main goal of the current analysis is to highlight the reliability of the trend of H 0 (z) and not to further constrain Ω 0m or any other cosmological parameters. With the subsequent fitting of H 0 values through the model g(z) = H 0 /(1 + z) η , we obtain η ∼ 10 −2 , as in the previous work [51]: those are compatible with zero from 1.2 to 5.8 σ (see Table 1). The multidimensional results could reveal a dependence on the redshift of H 0 , assuming that it is observable at any redshift scale. If this evolution is not caused by statistical effects and other selection biases or hidden evolution of SNe Ia parameters [364], we show how H 0 (z) could modify the luminosity distance definition within the modified theory of gravity. If we consider a theoretical interpretation for the observed trend, new cosmological scenarios may explain an evolving Hubble constant with the redshift. For instance, we test in Section 6.2, and Section 7 a simple class of modified gravity theories given by the f (R) models in the equivalent scalar-tensor formalism. In principle, this could be due to an effective varying Einstein constant governed by a slow evolution of a scalar field which mediates the gravity-matter interaction. However, the slow decreasing trend of H 0 has proven to be independent of the Hu-Sawicki model application. Indeed, if this theory had worked we would have observed the trend of the η parameter to be flattened out and be compatible with 0 in 1 σ at any redshift bin. This is not the case, thus new scenarios must be explored within the modified theories of gravity or slightly alternative approaches (see Sec. 8.2). We can state that this evolving trend of H 0 is independent of the starting values of the fitting for H 0 (we here have considered H 0 = 70) and, thus, on the fiducial M and on the redshift bins and even when we consider two cosmological parameters changing contemporaneously (Ω 0m and w a in ΛCDM and w 0 w a CDM models, respectively). Thus, we need to further investigate the nature of this trend. In addition, the implementation of GRBs as cosmological probes together with SNe Ia and BAOs has proven to be not only possible in a near future but also necessary since the redshift range that GRBs cover is much larger than the one typical of SNe Ia. This last characteristic will surely allow GRBs to give further information on the nature of the early universe and pose new constraints in the future measurements of H 0 .

Conflicts of Interest:
The authors declare no conflict of interest. Notes 1 https://github.com/dscolnic/Pantheon (accessed on 21 December 2020). 2 The code is available upon request.