On the Hubble Constant Tension in the SNe Ia Pantheon Sample

The standard cosmology is based on the well-known ΛCDM (cold dark matter) model, which relies on the existence of a cosmological constant Λ (Carroll 2001) with an equation-of-state parameter w = −1 and a CDM component. This model is the most widely accepted paradigm to explain the structure and evolution of the late universe. The discovery of the accelerating expansion phase (Riess et al. 1998; Perlmutter et al. 1999) has suggested the presence of a cosmological constant as the most viable scenario to account for the observations of Type Ia supernovae (SNe Ia; Betoule et al. 2014; Scolnic et al. 2018), Cepheids (Riess et al. 2019), cosmic chronometer probes for the expansion rate of the parameter H(z) (Stern et al. 2010; Moresco et al. 2012; Gomez-Valent & Amendola 2018), cosmic microwave background (CMB) fluctuations (Bennett et al. 2003; Hinshaw et al. 2013; Ade et al. 2016; Aghanim et al. 2020), baryon acoustic oscillations (BAO; Aubourg et al. 2015; Alam et al. 2017), large-scale matter perturbations observed through redshift space distortions (Basilakos & Nesseris 2016; Aghanim et al. 2020; Quelle & Maroto 2020), and weak lensing (WL; Baxter et al. 2016, Efstathiou & Lemos 2018). Among these probes, Cepheids and SNe Ia are considered the most appealing standard candles: astrophysical objects whose luminosity is known or can be derived from well-known intrinsic relations between distance-independent and distance-dependent observables. The luminosity depends on the luminosity distance, a quantity expressed by a given cosmological model (see Equation (8)).

Despite the outstanding results and predictions, the ΛCDM model must deal with open problems of a theoretical and observational nature. One of the biggest challenges in modern astrophysics is the so-called Hubble constant (H₀) tension: the 4.4σ discrepancy between the local value of H₀ based on Cepheids in the Large Magellanic Cloud (LMC), ${H}_{0}=74.03\,\pm 1.42\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ (Riess 2020), and the Planck data of the CMB radiation, ${H}_{0}=67.4\,\pm 0.5\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ (Aghanim et al. 2020). We here note that the discrepancy ranges from 4.4σ to more than 6σ in Riess et al. (2019), Wong et al. (2020), and Camarena & Marra (2020), depending on the combination of the local data used. A value similar to the one of the SNe Ia, ${H}_{0}\approx 72\pm 2\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, is reported by strong lens systems and time-delay measurements (Liao et al. 2019, 2020; Keeley et al. 2021). On the other hand, independent measurements of cosmic chronometers (based on models of evolving galaxy star luminosity) report the best-fit value of ${H}_{0}=67.06\,\pm 1.68\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ in Gomez-Valent & Amendola (2018), favoring the CMB and BAO measurements. Moreover, estimates of H₀ based on a combination of cosmological data, including calibration of the tip of the red giant branch on SNe Ia (Freedman et al. 2019), quasars (Risaliti & Lusso 2019), time-delay measurements, cosmic chronometers, and gamma-ray bursts (GRBs; Cardone et al. 2009, 2010; Dainotti et al. 2013b; Dutta et al. 2019; Krishnan et al. 2020a; Yang et al. 2020), report a value of H₀ that is between the CMB, BAO, and local measurements; for a review, see Di Valentino et al. (2021). Besides these results, there are also discrepancies greater than 2σ between the values of ${{\rm{\Omega }}}_{0m}$ from the ΛCDM model in Amati et al. (2019), when GRBs are adopted, and the ones obtained by the Pantheon sample, which is the largest compilation of spectroscopically confirmed SNe Ia so far. However, the discrepancy is only visible for some combination of a sample of GRBs and SNe Ia. For example, in the works of Cardone et al. (2009, 2010), Dainotti et al. (2013b), and Postnikov et al. (2014), where a combination of GRBs and SNe is adopted, the results agree with the SNe ones within 1σ. The status of all these discrepancies is summarized in Figure 1, where the values of H₀ between different probes are shown. Very recently, some works (Sharov & Sinyakov 2020) have investigated the Hubble constant tension by using different probes (cosmic chronometers) and different samples (the Pantheon sample, H(z), and the CMB data). Moreover, Renzi & Silvestri (2020) suggested a new method to measure H₀, which does not depend on the underlying cosmological model. To date, the reason for this discrepancy requires further investigation. Among several attempts, this discrepancy has been discussed in terms of a variation of the equation of the state for dark energy, dependent on redshift, $w=w(z)$, with a wide range of theories (Armendariz-Picon et al. 2000; Zimdahl et al. 2001; Moffat 2005; Grande et al. 2006; Caldera-Cabral et al. 2009; Yoo & Watanabe 2012; Benisty & Guendelman 2018; Mörtsell & Dhawan 2018; Anagnostopoulos et al. 2019). In particular, these theories include early (Karwal & Kamionkowski 2016; Agrawal et al. 2019; Hazra et al. 2019; Lambiase et al. 2019; Poulin et al. 2019; Keeley et al. 2020) and late dark energy models (Di Valentino et al. 2017; Li & Shafieloo 2019, 2020; Yang et al. 2019a, 2019b; Alestas et al. 2020; Vagnozzi 2020). Modified gravity theories (Ballardini et al. 2016; Lin et al. 2019; Rossi et al. 2019; Braglia et al. 2020; Escamilla-Rivera & Said 2020; Kazantzidis & Perivolaropoulos 2020) and alternative theories in which the speed of light is variable (Nguyen 2020) also attempt to explain this discrepancy.

Zoom In Zoom Out Reset image size

The ΛCDM model suffers from two main serious problems (Weinberg 1989; Peebles & Ratra 2003). The first one deals with the fact that the observed vacuum energy compared to the predicted theoretical expectation from quantum physics (fine-tuning problem) is extremely small: the ratio between the two values is 10¹²⁰. The second problem is that it is still unclear why the constant-energy and dark-matter densities are of the same order of magnitude today, while in the past their difference was 10⁹ in the CMB epoch (this is the so-called coincidence problem). Given the Friedmann acceleration equation in Weinberg (2008), the condition for the late-time acceleration is provided by the equation-of-state parameter of dark energy $w\lt -1/3$. However, the nature of dark energy is still unknown. On the other hand, different from the past, it is now possible to precisely constrain w and its evolution to more than 10%, due to improved determinations of cosmological distances. For instance, a useful parameterization is given by $w(z)\,={w}_{0}+{w}_{a}\times z/(1+z)$, according to the Chevallier–Polarski–Linder model (CPL; Chevallier & Polarski 2001; Linder 2003), where w₀ and w_a are parameters. To clarify the H₀ tension, it is crucial to analyze SNe Ia with a binning in redshifts so that the evolutionary patterns of either cosmological parameters or astrophysical observables may be revealed, if any.

This tension could be explained by internal inconsistencies in Planck data or SNe Ia systematics in the local determination of H₀ or with a new physics that lies beyond the standard cosmological model. There may be many ways to explain theoretically this difference between the various values of H₀ coming from many probes and methodologies. A probable scenario entails the use of teleparallel equivalent general relativity (Nunes et al. 2016; Nunes 2018; Benetti et al. 2021), where the torsion tensor assumes a key role, while another possible explanation relies on modified theories of gravity, thus requiring the use of a function of the Ricci scalar f(R), instead of R, which is the usual gravitational Lagrangian density (Sotiriou 2006; Nojiri & Odintsov 2007; Sotiriou & Faraoni 2010; Capozziello & De Laurentis 2011). Very recently, but after the submission of the current paper, Odintsov et al. (2021) proposed that a possible explanation of the H₀ tension can rely on some specific f(R) models with a functional form of the f(R) exponential or as a power law. This proposal, although similar to our suggestion, takes into account the full Pantheon sample and not a binned analysis, as we have here performed. A day after the submission of our paper, a similar conclusion was reached on the evolution of H₀ by Garcia-Quintero et al. (2020) considering a redshift binning on several probes including the SNe Ia. Another suggestion is to make changes in the early-time universe without modifying the late-time cosmology (Bernal et al. 2016). Other explanations are presented by Kazantzidis & Perivolaropoulos (2020) and Kazantzidis et al. (2021), who focus their analysis on the Pantheon sample. They claim that if their results are not due to statistical fluctuations, caused by the sample size inside each bin, the redshift evolution of the absolute magnitude, M, can be due to a local underdensity of matter that causes a higher value of H₀, or as a time variation of Newton's constant that implies an evolving Chandrasekhar mass and, thus, an evolving absolute magnitude of low-z SNe Ia. Indeed, the effects of the inhomogeneities on the local measurements of H₀ have been studied by Ben-Dayan et al. (2014) and Fanizza et al. (2021) and, in particular, the underdensities in the framework of Milgromian dynamics by Haslbauer et al. (2020) and Asencio et al. (2021). Additional analysis with the full Pantheon SNe Ia sample and other probes like the joint light-curve analysis (JLA; Betoule et al. 2014) have been performed to test the dependence of H₀ on redshift, called the redshift evolution, but their results have been obtained without the binning analysis and do not constrain the redshift evolution of H₀ (Di Valentino et al. 2020). To better clarify the existence of a possible evolutionary trend for H₀, we investigate this issue in a binned analysis of H₀ (in 3, 4, 20, and 40 bins in redshifts) with the SNe Ia Pantheon sample. This compilation includes surveys that can extend the Hubble diagram out to z = 2.26, from a dark-energy- to a dark-matter-dominated universe. Our approach is similar to what has been done by Kazantzidis & Perivolaropoulos (2020), but we investigate the evolution of H₀ instead of M. Moreover, our method is different from the one in Di Valentino et al. (2020), in which H₀ varies together with the free parameters of the evolutionary functions, and a full Pantheon sample analysis is performed, instead of a binning analysis. The advantage of our approach is that we can fit the value of H₀ in each redshift bin, and we show that there is a weak dependence of H₀ on redshift, which is not consistent with zero within more than 1σ.

This paper is built as follows. In Section 2 we introduce the ΛCDM and ${w}_{0}{w}_{a}$CDM models. In Section 3 we detail our method. In Section 4 we present the results for H₀ by dividing the sample into 3, 4, 20, and 40 redshift bins for the flat ΛCDM and ${w}_{0}{w}_{a}$CDM models starting from a local value of ${H}_{0}\,=73.5\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$. We choose this value because it is among the highest value of H₀ found in the SNe Ia probes at low redshift and other local probes, as shown in Figure 1. In Section 5 we discuss our results through astrophysical and theoretical interpretations. In Section 6 we present a summary and conclusions.

A homogeneous and isotropic universe is described by the Friedmann equations (Weinberg 2008; Montani et al. 2011), which are given by the evolution of the scalar factor $a\left(t\right)$:

$$\begin{eqnarray}&&{H}^{2}\left(t\right)={\left(\displaystyle \frac{\dot{a}\left(t\right)}{a\left(t\right)}\right)}^{2}=\displaystyle \frac{{c}^{2}\,\chi \,\rho \left(t\right)}{3}-\displaystyle \frac{{c}^{2}\,k}{{a}^{2}\left(t\right)}\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&\frac{\ddot{a}\left(t\right)}{a\left(t\right)}=-\frac{{c}^{2}\,\chi }{6}\left(\rho +3p\right)\end{eqnarray} \tag{ 1b }$$

where c is the speed of light, $\chi =8\pi G\,/{c}^{4}$ is the Einstein constant, k is the curvature parameter (k = 0 for a flat cosmology), and $\rho \left(t\right)$ and $p\left(t\right)$ are the total energy density and pressure, respectively.

Furthermore, the continuity equation

$$\begin{eqnarray}&&\dot{\rho }+3\,H\left(\rho +\displaystyle p\right)=0\end{eqnarray} \tag{ 2 }$$

allows us to describe the evolution of the matter–energy sources, and it can be obtained by combining Friedmann equations.

Within the framework of the ΛCDM model, the total energy density is given by

$$\begin{eqnarray}&&\rho ={\rho }_{m}+{\rho }_{r}+{\rho }_{{\rm{\Lambda }}},\end{eqnarray} \tag{ 3 }$$

where we denote different components with m = matter, r = radiation, and Λ = cosmological constant, while the pressure $p\left(\rho \right)=w\,\rho$ is referred to corresponding barotropic fluids (w = 0 for matter, $w=1/3$ for radiation, and $w=-1$ for the cosmological constant). Moreover, the total energy density (see Equation (3)) can also be expressed in terms of dimensionless cosmological density parameters ${{\rm{\Omega }}}_{i}\left(t\right)={\rho }_{i}\left(t\right)/{\rho }_{c}\left(t\right)$, where i denotes different components of ρ and ${\rho }_{c}\left(t\right)\,=3{H}^{2}\left(t\right)/{c}^{2}\,\chi$ denotes the critical energy density of the universe. Therefore, we can rewrite the dimensionless first Friedmann Equation (1a) as the following:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{m}\left(t\right)+{{\rm{\Omega }}}_{r}\left(t\right)+{{\rm{\Omega }}}_{{\rm{\Lambda }}}\left(t\right)+{{\rm{\Omega }}}_{k}\left(t\right)=1,\end{eqnarray} \tag{ 4 }$$

where ${{\rm{\Omega }}}_{k}\left(t\right)=-k\,{c}^{2}/{a}^{2}\left(t\right){H}^{2}\left(t\right)$. Equation (1a) can be recast to express the Hubble function H, as in Peebles (1971), in terms of the redshift z, which within the framework of the ΛCDM model is given by

$$\begin{eqnarray}&&H\left(z\right)={H}_{0}\,\sqrt{{{\rm{\Omega }}}_{0m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{0r}{\left(1+z\right)}^{4}+{{\rm{\Omega }}}_{0{\rm{\Lambda }}}+{{\rm{\Omega }}}_{0k}{\left(1+z\right)}^{2}}.\end{eqnarray} \tag{ 5 }$$

We here stress how the scale factor is related to the definition of the cosmological redshift z:

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{0}}{a\left(t\right)}=1+z,\end{eqnarray} \tag{ 6 }$$

where the subscript 0 denotes the present time (z = 0), and it should be emphasized that H₀ and ${{\rm{\Omega }}}_{0i}$ are constants. Since the relativistic components are subdominant in the present universe, ${{\rm{\Omega }}}_{0r}\approx {10}^{-5}$ is usually neglected at late times. Once $H\left(z\right)$ is known, the general expression of the luminosity distance as in Weinberg (2008) for a flat cosmology is given by

$$\begin{eqnarray}&&{d}_{L}(z)=c(1+z){\int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{H(z^{\prime} )}.\end{eqnarray} \tag{ 7 }$$

In particular, d_L (z) in the ΛCDM model becomes

$$\begin{eqnarray}&&{d}_{L}(z)=\displaystyle \frac{c(1+z)}{{H}_{0}}{\int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{\sqrt{{{\rm{\Omega }}}_{0m}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{0{\rm{\Lambda }}}}}.\end{eqnarray} \tag{ 8 }$$

In this framework, we ignore the relativistic components.

Regarding instead the wCDM models in which w evolves with redshift, that is, $w=w(z)$, the continuity equation associated with a dark energy component is written as

$$\begin{eqnarray}&&{\dot{\rho }}_{{DE}}+3\,H\,{\rho }_{{DE}}\left[1+w\left(z\right)\right]=0.\end{eqnarray} \tag{ 9 }$$

Solving this differential equation for ${\rho }_{{DE}}$ in terms of $w\left(z\right)$, we can obtain the corresponding ${{\rm{\Omega }}}_{{DE}}\left(z\right)$ and finally the Hubble function,

$$\begin{eqnarray}&&\,H\left(z\right)={H}_{0}\,\sqrt{{{\rm{\Omega }}}_{0m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{0{DE}}\,\exp \left[3\,{\int }_{0}^{z}\left[1+w\left({z}^{{\prime} }\right)\right]\displaystyle \frac{{{dz}}^{{\prime} }}{1+{z}^{{\prime} }}\right]},\,\end{eqnarray} \tag{ 10 }$$

for a flat cosmology in the late universe, neglecting the relativistic components. Assuming the condition $w=-1$, a cosmological constant is reproduced and the Hubble function for the ΛCDM model is recovered (see Equation (5)). For a simple linear model $w\left(z\right)={w}_{0}+{w}_{1}\left(1+z\right)$, the exponential term in Equation (10) would grow increasingly unsuitable at redshift $z\gg 1$. Thus, for high z, although many models have been proposed, we will focus on the CPL parameterization or ${w}_{0}{w}_{a}\mathrm{CDM}$ model: $w(z)={w}_{0}+{w}_{a}\times z/(1+z)$. To have a slight deviation from the cosmological constant and a slow evolution with redshift, the values for the parameters are usually ${w}_{0}\sim -1$ and ${w}_{a}\sim 0$. The Hubble function $H\left(z\right)$ in Equation (10) using the CPL parameterization for $w\left(z\right)$ becomes

$$\begin{eqnarray}&&H\left(z\right)={H}_{0}\,\sqrt{{{\rm{\Omega }}}_{0m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{0{DE}}{\left(1+z\right)}^{3\left(1+{w}_{0}+{w}_{a}\right)}\,{e}^{-3{w}_{a}\tfrac{z}{1+z}}}.\end{eqnarray} \tag{ 11 }$$

Once $H\left(z\right)$ is known, one can easily write the luminosity distance ${d}_{L}\left(z\right)$ using Equation (7).

The peculiarity of SNe Ia is their nearly uniform intrinsic luminosity with an absolute magnitude around M ∼ −19.5 (Carroll 2001), and this allows us to promote SNe Ia to a well-established class of standard candles. To evaluate the best cosmological model underlying our universe, we make use of the distance modulus, μ, derived from the observations of SNe Ia, and we compare it with the theoretical ${\mu }_{\mathrm{th}}$, defined as follows:

$$\begin{eqnarray}&&{\mu }_{\mathrm{th}}=m-M=5\hspace{0.5ex}{\mathrm{log}}_{10}\ {d}_{L}(z,{{\rm{\Omega }}}_{0m},{H}_{0},{w}_{0},{w}_{a})+25,\end{eqnarray} \tag{ 12 }$$

where m is the apparent magnitude of the source, M is the absolute magnitude, and d_L is the luminosity distance expressed in Mpc, defined in Equation (7).

The difficulty in the determination of the cosmological parameters lies in the identification of M, due to different sources of systematics and statistical errors, like Milky Way extinction, microlensing effects, and selection biases, as detailed in Scolnic et al. (2018). The Pantheon sample is a compilation of 1048 spectroscopically confirmed SNe Ia that gathers different surveys. In Scolnic et al. (2018), the ${\chi }^{2}$ approach needs the definition of the distance moduli, ${\mu }_{\mathrm{obs}}$, obtained by the observations:

$$\begin{eqnarray}&&{\mu }_{\mathrm{obs}}={m}_{B}-M+\alpha {x}_{1}-\beta c+{\rm{\Delta }}M+{\rm{\Delta }}B,\end{eqnarray} \tag{ 13 }$$

where x₁ is the stretch parameter, c is the color, m_B is the B-band apparent magnitude, M is the absolute magnitude in the B band of a reference SN with ${x}_{1}=0$ and c = 0, ${\rm{\Delta }}M$ is a distance correction based on the host-galaxy mass of the SN, and ${\rm{\Delta }}B$ is a bias correction based on previous simulations. The coefficients α, β, and ${\rm{\Delta }}M$ are allowed to freely vary to be optimized in the approach presented by Scolnic et al. (2018). This study on the Pantheon sample requires the implementation of the beams with bias correction (BBC) method (Scolnic & Kessler 2016) to create a Hubble diagram corrected for selection biases. As explained in Scolnic et al. (2018) and Tripp (1998), there is a degeneracy between H₀ and M. It should be emphasized that in the Pantheon release, the absolute magnitude is fixed to $M=-19.35$ such that ${H}_{0}=70.0\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$. The value of $M=-19.35$ can be derived from Scolnic et al. (2018), computing M from Equation (13).

In this paper, H₀ is not derived through the BBC method, but it is obtained by fixing the value of ${{\rm{\Omega }}}_{0m}$ to a fiducial value found in Scolnic et al. (2018) and comparing directly the quantity ${\mu }_{\mathrm{obs}}$ tabulated in Scolnic et al. (2018) with the ${\mu }_{\mathrm{th}}$ for each SN.

We now introduce a slight modification for the computation of luminosity distance (see Equation (8)), which in the case of SNe Ia is more precise according to Kenworthy et al. (2019):

$$\begin{eqnarray}&&{d}_{L}({z}_{{hel}},{z}_{{HD}})=\displaystyle \frac{c(1+{z}_{{hel}})}{{H}_{0}}{\int }_{0}^{{z}_{{HD}}}\displaystyle \frac{{dz}^{\prime} }{\sqrt{{{\rm{\Omega }}}_{0m}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{0{\rm{\Lambda }}}}},\end{eqnarray} \tag{ 14 }$$

where z_hel is the heliocentric redshift, and z_HD is the corrected CMB redshift, or "Hubble diagram" redshift, which takes into account the peculiar velocity corrections.

Also, different models for the stretch and color of the SNe population can be applied to a given sample: C11 (Chotard et al. 2011) and G10 (Guy et al. 2010) are the most appropriate models. C11 is composed of 75% chromatic variation and 25% achromatic variation, and the result of C11 suggests that the dispute in interpreting SN Ia colors and their compatibility with a classical extinction law can be solved with the dispersion in colors and by the variability of features present in SN Ia spectra. The G10 model is composed of 30% chromatic and 70% achromatic variation. G10 concludes that there is no clear evidence for a possible redshift evolution of the slope β of the color–luminosity relation. Scolnic et al. (2018) point out that for the G10 scatter model the relative bias of m_B with redshift is small compared to the relative bias of color with redshift, while the opposite is true for the C11 model. The distance biases in Scolnic et al. (2018) agree to 1% when applied with the G10 and the C11 models, except that at high redshift the divergence of distance bias between the two models starts to weigh more. This is caused by the selection criteria used to gather SNe according to their color and magnitude. Since, in principle, there are no reasons to prefer one model over the other, the average of G10 and C11 bias corrections is taken, and this constitutes the systematic part of the covariance matrix, denoted as C_sys (Scolnic et al. 2018). The average of the two models is here considered when we apply the binned approach. In our analysis for consistency, we also average the two values of ${\mu }_{\mathrm{obs}}$ according to G10 and C11 presented in Equation (13). In this work, we base our analysis on Scolnic et al. (2018), and the idea of binning is similar to the one used by Kazantzidis & Perivolaropoulos (2020). The main part of our investigation is different from Kazantzidis & Perivolaropoulos (2020) because we use more bins, we average the systematic uncertainties of G10 and C11, and we investigate H₀ rather than M.

In this section, we present our analysis. We define the distance residuals ${\rm{\Delta }}\mu ={\mu }_{\mathrm{obs}}-{\mu }_{\mathrm{th}}({H}_{0},\ldots )$ where ${\mu }_{\mathrm{obs}}$ is given by Equation (13), and it is taken from the repository in Scolnic et al. (2018). ¹¹ So, the ${\chi }^{2}$ is defined as

$$\begin{eqnarray}&&{\chi }^{2}={\rm{\Delta }}{\mu }^{T}{{ \mathcal C }}^{-1}{\rm{\Delta }}\mu .\end{eqnarray} \tag{ 15 }$$

To reduce the uncertainty on the derivation of H₀ in each redshift bin, we directly use the values of ${\mu }_{\mathrm{obs}}$ detailed in Equation (13), while the values of ${\mu }_{\mathrm{th}}$ are given by Equation (12). In Equation (15), ${\rm{\Delta }}\mu$ is the distance residual vector containing 1048 SNe from the Pantheon sample, and ${ \mathcal C }$ is the 1048 × 1048 full covariance matrix defined as follows:

$$\begin{eqnarray}&&{ \mathcal C }={C}_{\mathrm{sys}}+{D}_{\mathrm{stat}},\end{eqnarray} \tag{ 16 }$$

where C_sys is a matrix that contains the systematic sources of errors, and D_stat is a diagonal matrix that includes the total distance errors associated with every SN. The latter takes into account the contributions from photometric error, mass step correction, bias, peculiar velocity and redshift in quadrature, stochastic gravitational lensing, and intrinsic scatter (Scolnic et al. 2018). We have reconstructed the full covariance matrix as expressed in Equation (16), and we have divided it into submatrices of three and four different bins in redshift. More specifically, we have taken 1048 SNe, ordered them by redshift, and then divided the ordered 1048 SNe into three and four bins with equally populated subsamples composed of $\approx 349$ SNe and 262 SNe, respectively. The ranges of redshift are $0.0101\lt z\lt 0.1769$, $0.1771\lt z\lt 0.3374$, and $0.3375\lt z\lt 2.2600$ for the three bins, while we have $0.0101\lt z\lt 0.1299$, $0.1323\lt z\lt 0.2485$, $0.2486\lt z\,\,\lt 0.4224$, and $0.4235\lt z\lt 2.2600$ for the four bins. In Equation (15), it is clear that the full covariance submatrices for a bin with N SNe should have dimensions N × N to compute the ${\chi }^{2}$.

If we focus only on the statistical contribution given in D_stat, it is straightforward to build the submatrices, since D_stat is diagonal, so we can easily associate a D_stat element with a single SN. However, the presence of the C_sys matrix, which is not diagonal, led us to write a customized code ¹² that extracts the submatrices, including also systematic errors. This procedure has been developed by extracting only the full covariance matrix elements associated with the SNe with redshift inside the considered bin. This is an improvement in the precision of the results compared to the case of only statistical uncertainties (Kazantzidis & Perivolaropoulos 2020). This analysis is similar to Kazantzidis et al. (2021), in which also a full matrix of systematics has been taken into account.

The choice of three bins is dictated by the fact that these bins have a relatively high number of SNe to still provide statistical representative subsamples, while the four-bin analysis is performed to compare our results with the ones from Kazantzidis & Perivolaropoulos (2020). We also stress here that we use the submatrices that contain both statistical and systematic uncertainties, and, according to Scolnic et al. (2018), many of the evolutionary effects associated with systematic uncertainties are on the 1% level. Given that SNe Ia ${\mu }_{\mathrm{obs}}$ are usually measured with a precision of ≈15%, it is challenging to adequately account for these effects without hundreds of SNe in a given sample or bin. This means that if we have a redshift evolution of the parameters, such as the stretch and color, and we have a sample with less than hundreds of SNe, this evolution would be underestimated; in fact, it would be partially, or even totally, masked out.

After we have performed the bin selection, we use the Cobaya package available in Python (Torrado & Lewis 2020) to minimize the ${\chi }^{2}$ in Equation (15). We find the best values for the cosmological parameter H₀, which is left to vary as a nuisance parameter, while we fix in each bin ${{\rm{\Omega }}}_{0m}=0.298\pm 0.022$ to the fiducial value taken from Scolnic et al. (2018) for a flat ΛCDM cosmology. The choice of performing a one-dimensional analysis by fixing the value of ${{\rm{\Omega }}}_{0m}$ is necessary to constrain the parameters of a given redshift evolutionary function, g(z), for H₀. This function will be discussed in detail in Section 4.1.

Then, we perform a Markov Chain Monte Carlo (MCMC) analysis using the D'Agostini method (D'Agostini 1994) to sample a posterior distribution and obtain the confidence intervals of the H₀ parameter at the 68% and 95% levels. This analysis is repeated for three and four bins. However, one might argue that this approach does not always guarantee that the value of ${{\rm{\Omega }}}_{0m}$, which we fix as a fiducial, remains consistent within 1σ or 2σ with the values obtained by applying the same analysis to the total Pantheon sample. To this end, we perform the binned analysis with three and four bins by varying ${{\rm{\Omega }}}_{0m}$ and H₀ contemporaneously for a flat ΛCDM model. The results obtained for ${{\rm{\Omega }}}_{0m}$ are consistent in 1σ for three bins and in 2σ for four bins with the values obtained when we employ the total Pantheon sample. Our derived value of ${{\rm{\Omega }}}_{0m}=0.298\pm 0.016$ from the Pantheon sample is consistent within 1σ with the value obtained in Scolnic et al. (2018) from the Pantheon sample itself, thus validating our approach. Then, this method guarantees us that the precompiled covariance matrix ${ \mathcal C }$, given by Scolnic et al. (2018), can be reliably used.

We here clarify the selection criteria for the choice of the binning division: the contours for the parameters H₀ and ${{\rm{\Omega }}}_{0m}$ should constrain these parameters so that ${{\rm{\Omega }}}_{0m}$ for every single bin in the groups of three and four bins must be compatible at least within 2σ with the value of ${{\rm{\Omega }}}_{0m}$ relative to the full Pantheon sample. The results are summarized in Figure 2. Looking at three bins (see the left panel of Figure 2), we have the most favored framework: all three bins show closed contours in these intervals: $0\,\lt \,{{\rm{\Omega }}}_{0m}\lt 1$, 60 km s⁻¹ Mpc⁻¹ $\lt {H}_{0}\,\lt 80$ ${{\rm{km}}{\rm{s}}}^{-1}\,{{\rm{Mpc}}}^{-1}$. Furthermore, in this case the value of ${{\rm{\Omega }}}_{0m}$ is consistent within 1σ with the value of ${{\rm{\Omega }}}_{0m}$ related to the full sample, shown in red in the left panel of Figure 2. The four-bin scenario (see the right panel of Figure 2) has been introduced to mimic the bin division of Kazantzidis & Perivolaropoulos (2020), and the results are quite satisfactory, but not all bins are compatible within 1σ with the reference value of the total Pantheon, rather reaching compatibility only in 2σ, as shown in the red contours of the right panel of Figure 2. We performed an additional analysis that has demonstrated that the division of the Pantheon sample into more than four bins leads to incompatible values of ${{\rm{\Omega }}}_{0m}$ with the fiducial value of the total Pantheon sample in 2σ. Thus, this analysis leads us to conclude that in the case of standard cosmological parameters (${{\rm{\Omega }}}_{0m}=0.298,{H}_{0}\,=73.5\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$), the optimal number of bins in which we can divide the Pantheon sample is at most four, thus strongly disfavoring the possibility of subsequent divisions in bins. Nevertheless, in Section 4.1 we show what the results of our analysis would look like in the cases of 20 and 40 bins, here considered very extreme cases of binning.

Zoom In Zoom Out Reset image size

In the Pantheon release and the subsequent analysis, the absolute magnitude of SNe Ia is set to $M=-19.35$, related to a value of ${H}_{0}=70.0\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, but here we assume a different reference value for M according to the measured value of H₀ in the local probes.

4.1. Evolution of H₀ with Redshift in the Flat ΛCDM and w₀ w _a CDM Models

We here investigate if there is a redshift evolution of ${H}_{0}\left(z\right)$, obtained from the redshift binned analysis of SNe Ia in three and four bins, both in the ΛCDM and in the wCDM models. For the latter, we adopt the CPL parameterization, namely the ${w}_{0}{w}_{a}$CDM model. Analogous to the approach for the ΛCDM model, we fix the parameters in the ${w}_{0}{w}_{a}$CDM model to achieve constraints and study the evolution of ${H}_{0}\left(z\right)$. Specifically, the fiducial values of ${w}_{0}=-1.009$, ${w}_{a}\,=-0.129$, and ${{\rm{\Omega }}}_{0m}=0.308$ are chosen according to the results in Scolnic et al. (2018) for the w₀ w_aCDM model combining the SNe and CMB analysis. We compare our procedure with the one obtained by Di Valentino et al. (2017), in which a possible redshift dependence of the intrinsic SNe Ia luminosity, M, is investigated. Our computational approach is different from Di Valentino et al. (2017), where M is varied simultaneously with other parameters. We note that by reducing the degrees of freedom of the parameter space in our analysis, we can constrain a value for the parameter of the evolution: different from Di Valentino et al. (2017), we focus on a single parameter. Some alternative solutions have been proposed in Zhang et al. (2017), where more probes are included such as Cepheids, rather than only SNe Ia, and the fit is performed over the difference between H₀ and M.

For both the ΛCDM and the ${w}_{0}{w}_{a}$CDM models, we set the following priors for MCMC: 60 ${\mathrm{kms}}^{-1}\,{\mathrm{Mpc}}^{-1}$ $\lt {H}_{0}\,\lt 80$ ${\mathrm{kms}}^{-1}\,{\mathrm{Mpc}}^{-1}$. Once we have obtained the values of H₀ for our bins, we perform a nonlinear fit of H₀ with the following functional form:

$$\begin{eqnarray}&&g(z)={H}_{0}(z)=\displaystyle \frac{{\tilde{H}}_{0}}{{\left(1+z\right)}^{\alpha }},\end{eqnarray} \tag{ 17 }$$

where ${\tilde{H}}_{0}$ and α are the fitting parameters, and the α coefficient indicates the evolutionary trend. We here remark that the choice of this function g(z) is standard for characterizing the evolution of many astrophysical sources, and it is widely used for GRBs and quasars (Lloyd & Petrosian 2000; Singal et al. 2011, 2013; Dainotti et al. 2013b, 2015, 2016, 2017b, 2020; Petrosian et al. 2015; Dainotti & Del Vecchio 2017 ). According to the functional form of the fit for z = 0, ${\tilde{H}}_{0}={H}_{0}$, we estimate the confidence interval at 68%, thus obtaining a 1σ error. In Table 1 we show the fit parameters, ${\tilde{H}}_{0}$ and α, as a function of three and four redshift bins in the case of ΛCDM and the w₀ w_aCDM models. In the case of the ΛCDM model, we observe a decreasing trend of H₀ in three bins for which the α parameter is compatible with zero only in 2.0σ, while for four bins it is consistent with zero only at 1.5σ, as shown in the left panel of Figure 5. Similarly, for the ${w}_{0}{w}_{a}$CDM model, a decreasing trend of ${H}_{0}(z)$ with redshift can be seen; see the right panel of Figure 5. It should be emphasized that also for the ${w}_{0}{w}_{a}$CDM model, α is consistent with zero only at the level of 1.9σ for three bins and only in 1.2σ for four bins, thus showing an evolution of H₀ with redshift. Thus, we have investigated this evolutionary trend also within the framework of a wCDM model to verify the hypothesis of whether this trend disappears when considering a modified equation-of-state parameter w for the dark energy. In other words, we are wondering if the observed evolution of ${H}_{0}\left(z\right)$ is due to the $w\left(z\right)$. Our results point out that H₀ exhibits an evolution in both the ΛCDM and ${w}_{0}{w}_{a}$CDM models. Thus, a varying equation-of-state parameter $w\left(z\right)={w}_{0}+{w}_{a}\ast z/(1+z)$ could not imply our results for ${H}_{0}\left(z\right)$.

Table 1. Fit Parameters for ${H}_{0}(z)$ for 3, 4, 20, and 40 Bins, Assuming a Flat ΛCDM Model (Upper Part) with Fixed ${{\rm{\Omega }}}_{0m}=0.298$ and a Flat ${w}_{0}{w}_{a}$CDM Model (Lower Part) with Fixed Parameters ${w}_{0}=-1.009$, ${w}_{a}=-0.129$, and ${{\rm{\Omega }}}_{0m}=0.308$

Flat ΛCDM Model, Fixed ${{\rm{\Omega }}}_{0m}$, with Full Covariance Submatrices ${ \mathcal C }$
Bins	${\tilde{H}}_{0}$	α	$\tfrac{\alpha }{{\sigma }_{\alpha }}$	M	${H}_{0}\left(z=11.09\right)$	${H}_{0}\left(z=1100\right)$	$\%$ Tension
	$\left({{\rm{km}}{\rm{s}}}^{-1}\,{{\rm{Mpc}}}^{-1}\right)$				$\left({{\rm{km}}{\rm{s}}}^{-1}\,{{\rm{Mpc}}}^{-1}\right)$	$\left({{\rm{km}}{\rm{s}}}^{-1}\,{{\rm{Mpc}}}^{-1}\right)$	Reduction
3	73.577 ± 0.106	0.009 ± 0.004	2.0	−19.245 ± 0.006	72.000 ± 0.805	69.219 ± 2.159	54%
4	73.493 ± 0.144	0.008 ± 0.006	1.5	−19.246 ± 0.008	71.962 ± 1.049	69.271 ± 2.815	66%
20	73.222 ± 0.262	0.014 ± 0.010	1.3	−19.262 ± 0.014	70.712 ± 1.851	66.386 ± 4.843	68%
40	73.669 ± 0.223	0.016 ± 0.009	1.8	−19.250 ± 0.021	70.778 ± 1.609	65.830 ± 4.170	57%
Flat ${w}_{0}{w}_{a}$CDM Model, Fixed ${{\rm{\Omega }}}_{0m}$, with Full Covariance Submatrices ${ \mathcal C }$
Bins	${\tilde{H}}_{0}$	α	$\tfrac{\alpha }{{\sigma }_{\alpha }}$	M	${H}_{0}\left(z=11.09\right)$	${H}_{0}\left(z=1100\right)$	$\%$ Tension
	$\left({{\rm{km}}{\rm{s}}}^{-1}\,{{\rm{Mpc}}}^{-1}\right)$				$\left({{\rm{km}}{\rm{s}}}^{-1}\,{{\rm{Mpc}}}^{-1}\right)$	$\left({{\rm{km}}{\rm{s}}}^{-1}\,{{\rm{Mpc}}}^{-1}\right)$	Reduction
3	73.576 ± 0.105	0.008 ± 0.004	1.9	−19.244 ± 0.005	72.104 ± 0.766	69.516 ± 2.060	55%
4	73.513 ± 0.142	0.008 ± 0.006	1.2	−19.246 ± 0.004	71.975 ± 1.020	69.272 ± 2.737	65%
20	73.192 ± 0.265	0.013 ± 0.011	1.9	−19.262 ± 0.018	70.852 ± 1.937	66.804 ± 5.093	72%
40	73.678 ± 0.223	0.015 ± 0.009	1.7	−19.250 ± 0.022	70.887 ± 1.595	66.103 ± 4.148	59%

Notes. The first column indicates the number of bins, and the second and third columns denote the fit parameters, ${\tilde{H}}_{0}$ and α, contained in the g(z) function, according to Equation (17). The fourth column denotes the consistency of the evolutionary parameter α with zero in terms of 1σ, which is represented by the ratio $\alpha /{\sigma }_{\alpha }$. The fifth column indicates the new fiducial absolute magnitude M, such that ${H}_{0}=73.5\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$. In the sixth and seventh columns, we show the extrapolated values of ${H}_{0}\left(z\right)$ and the corresponding errors at the redshift of the most distant galaxies, z = 11.09, and the last scattering surface, z = 1100. The last column denotes the percentage Hubble constant tension reduction: $\% \mathrm{Diff}=1-{x}_{f}/{x}_{i}$. The definition of x_i and x_f follows Equations (19) and (18). All uncertainties are given in 1σ.

Download table as:  ASCIITypeset image

To investigate the extent of this trend of H₀, we extrapolate the function ${H}_{0}(z)$ up to the redshift of the most distant galaxies, z = 11.09 (Oesch et al. 2016), assuming that this trend could be observed also in other high-redshift probes. We find that the extrapolated value at z = 11.09 is consistent within 1σ with the value of H₀ obtained with the Planck measurement for both the ΛCDM and ${w}_{0}{w}_{a}$CDM models. The values of this extrapolation are reported in Table 1. We here stress that the error bars on the extrapolated H₀ values are large essentially because we propagate the errors on both ${\tilde{H}}_{0}$ and α. The choice of using the most distant galaxies is adopted because there are objects such as GRBs that can be hosted in the most distant galaxies and can be observed in principle up to z = 20 (Lamb 2007), thus allowing us to add them to the SNe Ia to further discuss the H₀ tension. Indeed, very recently, there was work (Jiang et al. 2021) claiming that the most distant galaxy at z = 11.09 hosts a GRB. Although this is far beyond the scope of this paper, it is important to show in the future the contribution of the GRBs in the last redshift bins of the SNe Ia, where the redshift of the sample is higher.

Very interestingly, if we extrapolate ${H}_{0}(z)$ to the last scattering surface at z = 1100, we obtain again values of H₀ compatible with the ones from Planck measurements, in 1σ for three and four bins for the ΛCDM and ${w}_{0}{w}_{a}$CDM models (see Table 1). To evaluate how much our results from fitting may reduce the H₀ tension, we compute a percentage difference (%Diff), taking into account the local and high-redshift values of ${H}_{0}(z)$, in the following way: ${\rm{ \% }}{\rm{Diff}}=1-{x}_{f}/{x}_{i}$. We denote with x_f our fit tension between the local value ${\tilde{H}}_{0}$, at z = 0, and the extrapolated value of H₀ at z = 1100, obtained from our fitting procedure. More precisely, x_f is given in terms of σ by

$$\begin{eqnarray}&&{x}_{f}=\displaystyle \frac{{\tilde{H}}_{0}(z=0)-{H}_{0}(z=1100)}{\sqrt{{\sigma }_{{\tilde{H}}_{0}(z=0)}^{2}+{\sigma }_{{H}_{0}(z=1100)}^{2}}}.\end{eqnarray} \tag{ 18 }$$

We denote with x_i the following term:

$$\begin{eqnarray}&&{x}_{i}=\displaystyle \frac{{H}_{0}^{(\mathrm{Cepheids})}(z\sim 0)-{H}_{0}^{(\mathrm{CMB})}(z\sim 1100)}{\sqrt{{\sigma }_{{H}_{0}^{(\mathrm{Cepheids})}(z\,\sim \,0)}^{2}+{\sigma }_{{H}_{0}^{(\mathrm{CMB})}(z\,\sim \,1100)}^{2}}},\end{eqnarray} \tag{ 19 }$$

where ${H}_{0}^{(\mathrm{Cepheids})}(z\sim 0)=74.03\pm 1.42\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ is the local measurement provided by Cepheids in the LMC, and ${H}_{0}^{(\mathrm{CMB})}(z\sim 1100)=67.4\pm 0.5\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ is the Planck data of the CMB radiation.

In the last column of Table 1 we summarize these percentage variations, indicating a sensitive reduction of the H₀ tension. Furthermore, it should be emphasized that the decreasing trend ${H}_{0}(z)$ not only reduces the tension but also provides a new way to approach the problem.

Besides, assuming that the trend we found for H₀ is intrinsic and not caused by hidden evolution or selection biases of other parameters at play in the SNe Ia sample, we would need to account for this new definition in the luminosity distance for a flat ΛCDM model in the following way:

$$\begin{eqnarray}&&{d}_{L}(z)=\displaystyle \frac{c(1+z)}{{\tilde{H}}_{0}}\,{\int }_{0}^{z}\displaystyle \frac{{\left(1+{z}^{{\prime} }\right)}^{\alpha }\,{dz}^{\prime} }{\sqrt{{{\rm{\Omega }}}_{0m}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{0{\rm{\Lambda }}}}}.\end{eqnarray} \tag{ 20 }$$

In Figure 3 we show how the corrected luminosity distance (blue line for three bins, green line for four bins), which takes into account the dependence of ${H}_{0}\left(z\right)$, deviates from the standard luminosity distance (red line) defined in Equation (8) in the ΛCDM model. It is visible that at high z (z = 11.09), there is an overestimation of $2.2 \%$ and 2% of the corrected luminosity distance computed in a ΛCDM model from Equation (20) compared to the standard luminosity distance in Equation (8) for three and four bins, respectively.

Zoom In Zoom Out Reset image size

To show that our results are reliable and independent of the number of bins, we have also added more bins, 20 and 40, although we have previously clarified that the optimal number of bins should not exceed three or four in order to not reduce considerably the number of SNe Ia in each bin. It should be remembered that we are using equally populated redshift bins. The last bin of 40 would leave us indeed with only 26 SNe Ia in each of the redshift ranges and 34 in the last bin. Note that this sample in the 40 bins is an extreme choice because it is even smaller than the sample adopted by Perlmutter et al. (1999), who used 42 SNe Ia. Nevertheless, this analysis confirms the reliability of our results. Since there is a degeneracy between the absolute magnitude M of SNe Ia and H₀, as already mentioned in Section 3, we set a reference value for ${H}_{0}=73.5\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ and we use the same systematic submatrices, despite those being calibrated with a standard cosmology of the previous analysis where ${H}_{0}=70\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$. To properly find the best value for M, a careful analysis is performed with a fixed value of ${H}_{0}=73.5\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ for the first redshift bin fixed coherently according to each binning division. We choose the first bin of the nearest SNe Ia among other possible bins in our analysis. The redshift ranges for the first bins are $0.0101\lt z\lt 0.1769$, $0.0101\lt z\lt 0.1299$, $0.0101\lt z\lt 0.0209$, and $0.0101\,\lt z\lt 0.0157$ for 3, 4, 20, and 40 bins, respectively. This work has been done for both ΛCDM and ${w}_{0}{w}_{a}$CDM models, using a ${\chi }^{2}$ minimization according to Equation (15) in each first bin. We obtain the absolute magnitude M for each binning division. The results for M are summarized in Table 1, and the values are all compatible with each other in 1σ. These analyses are summarized in Figure 4.

Zoom In Zoom Out Reset image size

Once the value of M is determined for the first bin, then this value of M is fixed for the analysis in the other bins, and new values of H₀ are computed in each bin. The results are shown in Figure 5 and in Figure 6. Interestingly, we have highlighted again a decreasing trend with redshift for H₀, according to the same functional form shown in Equation (17). The fit parameter results are listed in Table 1. It should be noted that the evolutionary parameters α are consistent with zero only at the 1.3σ level for the ΛCDM model and only at 1.7σ for the ${w}_{0}{w}_{a}$CDM model. The α parameters for the case of 20 and 40 bins are steeper because more bins in the fit procedure highlight small deviations of H₀ that are not visible in fewer bin divisions. Indeed, when we consider bins with hundreds of SNe, these small deviations are averaged, thus resulting in a flatter trend. Nevertheless, the α coefficients in the four cases (3, 4, 20, and 40) are all compatible with each other in 1σ, thus highlighting a persistent decreasing trend in the data. We here point out that we measure a change in H₀ of 0.4 over ${{\rm{\Delta }}}_{z}=1.0$ for three and four bins.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We obtain the extrapolated values of ${H}_{0}(z)$ at the redshift of the most distant galaxies (z = 11.09) and the redshift of the last scattering surface (z = 1100) with the respective H₀ tension percentage variations. The results are given in Table 1. Note that all of the extrapolated values of H₀ at z = 1100 are consistent within 1σ with the Planck measurement. However, these values are also consistent with local probes.

This double compatibility with the Planck results and the SNe Ia data could be derived from a not-fully-suited calibration of the systematic uncertainties with a different cosmology, and it requires further investigation with a new fitting of the SNe Ia light curves. We here stress that the systematic uncertainties are significantly affected by the particular cosmological model, as pointed out in Scolnic et al. (2018), and consequently by the value of the Hubble constant. In principle, the new systematic contribution should be computed rigorously through a reanalysis of the Pantheon sample data with the lightcurve fitter SALT2 and the BBC method according to a new reference cosmological model. However, this analysis goes far beyond the scope of the current paper.

In the following subsections, we discuss the role of selection biases in the Pantheon sample and the several possible theoretical interpretations related to our results.

5.1. Astrophysical Selection Biases

We here consider several effects that can play a role in our results. One of them is the presence of metallicity in SNe Ia. The average stellar ages and metallicities evolve with redshift, so it may happen that the average corrected SN Ia brightness at higher redshift will be fainter than the one at lower redshift if the observed bias is caused by the progenitor age or metallicity (Childress et al. 2013). This bias could affect in a nonnegligible way the estimation of cosmological parameters. Many authors have suggested different methods to encompass this problem. Sullivan et al. (2010) have suggested using host-galaxy mass as a third SN Ia brightness-correction parameter (after stretch and color), and this is done in Scolnic et al. (2018): many of the associated systematic uncertainties of these effects are on the 1% level. This tactic might improve any effects of luminosity caused by the SNe progenitor, while the possible intrinsic color discrepancy between SNe Ia in hosts of different metallicities may explain the observed bias of corrected SN Ia luminosity with host mass (and metallicity). Even if a strong correction for color–luminosity factor is applied (β ∼ 3), a considerable Hubble residual step, ${\rm{\Delta }}M$, is still observed (see Equation (13)) between high- and low-metallicity hosted SNe. This could mean that other effects are contributing to biases (Childress et al. 2013).

According to the stretch and color evolution, Childress et al. (2013) conclude that the stretch-corrected and color-corrected SN Ia Hubble residuals of SNe Ia in high-mass and low-mass host galaxies differ by 0.077 ± 0.014 mag, a result compatible with the one determined by Scolnic et al. (2018). The physical reason for that behavior may be linked to the interstellar dust, age, and metallicity of SNe: while the first cannot contribute alone to the observed bias, the last two can produce Hubble residual trends compatible with the ones observed in Childress et al. (2013). The most important concept that must be considered regarding the evolution of color and stretch factor concerns the changes of SNe Ia in the Pantheon sample. In Scolnic et al. (2018), to account for the stretch and color parameter evolution, both α and β are parameterized as a function of redshift in the following ways: $\alpha (z)={\alpha }_{0}+{\alpha }_{1}\ast z$ and $\beta (z)={\beta }_{0}+{\beta }_{1}\ast z$. This parameterization is not included in the ${\mu }_{\mathrm{obs}}$ from Equation (13) because there is no explicit evolutionary trend for $\alpha (z)$ and $\beta (z)$, so ${\alpha }_{1}=0$ and ${\beta }_{1}=0$. The only exception is that $\beta (z)$ shows different trends for high- and low-z SNe subsamples, and this is due to selection effects on SNe. So, Scolnic et al. (2018) included this uncertainty in the value of ${\beta }_{1}$ as an additional statistical uncertainty.

Regarding the properties of the host galaxies and the selection effects associated with these, several tests include probing the relations between luminosity and the properties of the host galaxies of the SNe (Kelly et al. 2010; Lampeitl et al. 2010; Sullivan et al. 2010) and analysis of the light-curve fit parameters of SNe and how these parameters relate to luminosity (Kessler & Scolnic 2017). Thus, in our analysis following the treatment of Scolnic et al. (2018), we can state that the color evolution is included, while the stretch evolution is not included because it turns out to be negligible. In previous studies (Tutusaus et al. 2017), the effect of considering or neglecting the evolution of SN parameters has been investigated, showing that the impact on cosmological parameters and models is not negligible. Recent studies on the SALT2.4 light-curve stretch in Nicolas et al. (2020) show that the basic SN stretch distribution evolves with redshift. In this work, the authors extract a particular sample of SNe starting from the Pantheon. The previously cited BBC method adopted by Scolnic et al. (2018) requires the implementation of an asymmetric Gaussian model based on the surveys as an SN stretch distribution, but, according to Nicolas et al. (2020), this modeling does not account for the redshift drift of the sample, and this effect is still present with the surveys reaching higher redshift values. Different from Scolnic et al. (2018), Nicolas et al. (2020) suggest another functional form for the stretch population:

$$\begin{eqnarray}\begin{array}{rcl}{x}_{1}(z) & = & \delta (z)\ast { \mathcal N }({\mu }_{1},{\sigma }_{1}^{2})+(1-\delta (z))\ast [a\ast { \mathcal N }({\mu }_{1},{\sigma }_{1}^{2})\\ & & +(1-a)\ast { \mathcal N }({\mu }_{2},{\sigma }_{2}^{2})]\end{array}\end{eqnarray} \tag{ 21 }$$

where a = 0.51, ${\mu }_{1}=0.37$, ${\mu }_{2}=-1.22$, ${\sigma }_{1}=0.61$, ${\sigma }_{2}\,=0.56$, and $\delta {(z)=({K}^{-1}\ast {(1+z)}^{-2.8}+1)}^{-1}$ with K = 0.87. Nevertheless, Nicolas et al. (2020) state that for cosmological purposes a great number of degrees of freedom for a given model can still be accepted in a sample that contains a large number of SNe Ia, like the Pantheon, but this choice does not allow us to successfully extract the SNe property distribution from a Malmquist-biased (Malmquist 1920) sample. The Malmquist bias effect is indeed an astrophysical effect that deals with all sources that are at a cosmological redshift. According to this effect, we cannot detect faint sources at high redshift, so there is a larger population of brighter sources at lower redshift. We point out that this effect could enter our results especially when we consider the highest redshift bin because the distances derived from SN are affected by selection effects, and at higher redshift, this selection weighs more. The simulations performed by Scolnic et al. (2018) for investigating the redshift evolution are limited up to z = 0.7. This redshift corresponds to the 18th bin out of 20 and the 36th bin out of 40, respectively. Thus, the remaining effect from $0.7\leqslant z\,\leqslant 2.26$ still needs to be investigated. In conclusion, many factors may cause the observed evolutionary trend for H₀ in our binning analysis: the drift of the stretch parameter with redshift is an excellent candidate given that the function we use here for the evolution, g(z), is a general function that can be used for any astrophysical source. The work we have performed is a way indeed to switch the evolution from the stretch to H₀, although surely the function $\delta (z)$ is different from g(z). Thus, with the current binned analysis, we have investigated the impact of the redshift evolution of the SN Ia population's intrinsic properties on cosmology. In particular, we have focused on H₀, thus adding more information and insights to the open discussion initiated by Nicolas et al. (2020). This scenario seems to be the most favorable, but it is not the only possible one. Thus, we would like to allow additional theoretical explanations that will be discussed in the next section.

5.2. Theoretical Interpretation

We here discuss some theoretical interpretations regarding the evolution of H₀, for instance by assuming possible hidden functions of redshift that are somehow related to H₀. In this section, we mainly explore two possible scenarios: the effect of the local inhomogeneity of the universe and the modified gravity theories.

From a theoretical point of view, the existing tension between the Planck observations and the SNe Ia data could suggest an effect associated with a local inhomogeneity of the universe. Since the observer measurements of cosmological observables are potentially sensitive to the local spacetime around the observer, the presence of inhomogeneity in the local universe could affect the cosmological parameters. For instance, the so-called void models place the observer inside a local underdensity of radius about $z\lt 0.15$, and this fact implies a locally measured Hubble constant that is larger than the global expansion rate: a perturbation in density causes a perturbation in the expansion rate, as shown in Kolb & Turner (1990), Marra et al. (2013), and Colgáin (2019). We could infer that the local universe may be underdense on spatial scales of several hundreds of megaparsecs if we consider that the matter underdensity is defined as underdense compared to the universe average density. In Keenan et al. (2013), Shanks et al. (2019a, 2019b), Boehringer et al. (2020), and Luković et al. (2020), there are shreds of evidence for a local matter underdensity on scales of roughly ∼300 Mpc with $\delta {\rho }_{0}/{\rho }_{0}\,\sim -0.2$. Particularly, using low-redshift distance estimators such as SNe Ia, the existence of a local matter underdensity would certainly reduce the H₀ tension, but it would not remove it (Kenworthy et al. 2019).

A similar approach consists of adopting the spherically symmetric Lemaître–Tolmann–Bondi inhomogeneous model (Lemaître 1933; Tolman 1934; Bondi 1947) to outline a possible dependence of the Hubble function also on the radial coordinate, in principle able to theoretically account for the presented behavior of ${H}_{0}\left(z\right)$. See, for instance, the analysis in Cosmai et al. (2019), where inhomogeneities are inferred to explain the whole universe acceleration phenomenon (Bonnor 1974; Ellis 2008; Buchert 2000). Such a scenario, when applied to the H₀ measurement tension, is expected to be associated with a typical spatial scale on which the weak underdensity must manifest itself. Moreover, we have to recover the ΛCDM model and the homogeneous value of H₀ at greater spatial scales, according to the Planck observations. We try another approach, different from that in Equation (17), to build empirically a reasonable fit function that can mimic the presence of a spatial scale associated with a local underdensity: ${H}_{0}(z)={g}_{2}(z)$, where ${g}_{2}(z)={K}_{1}\times (1+\exp (-{K}_{2}(1+z))$ and K₁ and K₂ are free parameters determined by the fit performed in the same way as in Section 4.1 and Section 4. If we have z = 0, this function, different from the previous g(z), does not lead to H₀. Here, ${K}_{2}^{-1}$ indicates a cutoff redshift scale such that at $z\gg {K}_{2}^{-1}$ the exponential term vanishes. The only possibility to obtain a constraint on the scale of the underdensity matter is to fix ${K}_{1}={H}_{0}$, obtained from Planck as an asymptotic value.

Moreover, the local underdensity requires the assumption of our unlikely special reference system inside the underdense region. On the contrary, we observe a natural decay of the H₀ value with redshift according to Equation (17), up to reconciling with the Planck measurement at the redshift of the last scattering surface $z\sim 1100$. In other words, we are observing a phenomenon whose origin does not seem to be associated with an intrinsic weak universe inhomogeneity, but it could be due to a different physical reason. If the H₀ tension will still persist, a new cosmology beyond the ΛCDM model may be necessary, since not even the void models can reconcile such a difference.

In this respect, since a proportionality in the Friedmann equation exists between the Hubble function $H\left(z\right)$ and the ΛCDM model density sources, mediated by the Einstein constant $\chi \equiv 8\,\pi \,G/{c}^{4}$ (G being the Newton constant), a possible dependence of H₀ on redshift can be naturally restated in terms of a dependence of χ on this same variable, that is, ${H}_{0}\left(z\right)\propto \sqrt{\chi \left(z\right){\rho }_{0}}$. Namely, we see that the value of the constant H₀ is decaying when measured by astrophysical sources at increasing z values, because the Einstein constant might decay with increasing redshift, and this effect is not accounted for in the ΛCDM model. The possible evolution of H₀ could point to an evolution of the Einstein constant χ. A similar scenario has also been inferred in Kazantzidis & Perivolaropoulos (2020), but by analyzing the SN Ia data with the scope of obtaining the evolution on ${ \mathcal M }$ rather than H₀.

Given this possible theoretical explanation, we investigate which modified gravity theory can mimic such a behavior of the Einstein constant with redshift. For instance, the evolution of G could be due to a Brans-Dicke theory (Li et al. 2015) or assuming a rapid transition at low redshift (Marra & Perivolaropoulos 2021).

According to the performed fitting, we note that we are predicting a very slow decay, as shown in Table 1, since we would have $\chi \propto {\left(1+z\right)}^{-2\alpha }$ with $\alpha \sim {10}^{-2}$.

A theoretical framework able to justify an effective dependence of the Einstein constant on time and, hence, on redshift, is provided by the so-called $f\left(R\right)$ model in the Jordan frame (Sotiriou 2006; Nojiri & Odintsov 2007; Sotiriou & Faraoni 2010; Capozziello & De Laurentis 2011). This theory restates a generalized Einstein–Hilbert action in terms of a generic function of the Ricci scalar R within the framework of a scalar–tensor formalism. In such a scheme, the degree of freedom associated with the form of $f\left(R\right)$ is translated into a scalar field ϕ nonminimally coupled to standard gravity with the resultant observed effect of ${H}_{0}(z)$. The action that describes this theory assumes the form

$$\begin{eqnarray}\begin{array}{rcl}S & = & -\displaystyle \frac{1}{2\,c\,\chi }\,\displaystyle \int {d}^{4}x\,\sqrt{-g}\left[\phi \,R-V\left(\phi \right)\right]\\ & & +\displaystyle \frac{1}{c}\,\displaystyle \int {d}^{4}x\,\sqrt{-g}\,{{ \mathcal L }}_{m},\end{array}\end{eqnarray} \tag{ 22 }$$

where $V\left(\phi \right)$ is the scalar field potential, related to $f\left(R\right)$, and a Lagrangian density ${{ \mathcal L }}_{m}$ for the matter source has been added to the action. Note that in the Jordan frame the quantity χ affects only the gravitational part of the action in Equation (22). It should be emphasized that the standard gravity interaction with a matter source is mediated via the scalar field $\phi \left(x\right)$, which implies an Einstein constant as a function of the coordinates, that is, $\chi \to \chi \,/\,\phi$. However, in addition to the rescaling of the Einstein constant, the nonminimally coupled scalar field in this scheme has its intrinsic evolution characterized by second-order derivatives entering the Einstein field equations.

Among all possible choices for the form of the $f\left(R\right)$ function, three proposals stand out for their ability to account for the universe acceleration: the so-called Hu-Sawicki (Hu & Sawicki 2007; Song et al. 2007), Starobinsky (Starobinsky 2007), and Tsujikawa (Amendola et al. 2007; Tsujikawa 2008) models.

Even more controversial is the question concerning the contribution carried by the kinetic term of the scalar field to the universe's total energy density. We suggest that this scenario can interpret our results, because the possibility of dealing with a significant universe acceleration requires a slow dynamics of the field, allowing that its potential term mimics a cosmological constant. In such a context, the nonminimally coupled scalar field has two main effects. The first relates to the rescaling of the Einstein constant, and the latter is able to provide a quasi-constant energy density responsible for the late acceleration of the universe. Moreover, the scalar field kinetic term, which contains second-order derivatives, is negligible and thus remains close to the ΛCDM model.

Comparing this proposal with our results, we see that $\phi \left(z\right)\sim {\left(1+z\right)}^{2\alpha }$ with $\alpha \sim {10}^{-2}$. This trend exactly suggests a scalar-field near-frozen dynamics, ensuring a very slow kinetic contribution to the universe energy density.

An additional proposal is provided by an equivalent description of the gravitational interaction, considering that the torsion accounts for gravitation. The f(R) theory of gravity and analog models are essentially metric theories. However, an open debate is related to identifying the correct variables describing the gravitational field. If tetrads, instead of metric, describe the gravitational field, then dynamics are given by torsion instead of curvature (notice that, in this context, the equivalence principle is not the foundation of the gravitational field, and affinities assume a fundamental role). These considerations led to the teleparallel formulation of general relativity, called teleparallel equivalent general relativity (TEGR). While curvature is used to geometrize the spacetime according to general relativity (GR), the TEGR attributes gravitation to torsion. In this perspective, the straightforward extension of the curvature gravity f(R) is now f(T), which extends TEGR (where T is the torsion scalar). The action for this theory is given in the following form:

$$\begin{eqnarray}&&{S}_{\mathrm{TEGR}}=\displaystyle \frac{1}{16\pi G}\int {d}^{4}{xe}\left[T+f(T)\right]+\displaystyle \frac{1}{c}\int {d}^{4}x\sqrt{-g}{{ \mathcal L }}_{m},\end{eqnarray} \tag{ 23 }$$

where e takes the place of $\sqrt{-g}$ and stands for the determinant of tetrad fields. In a Friedmann universe, the TEGR field equations lead to the following equation (Nunes et al. 2016; Nunes 2018; Benetti et al. 2021):

$$\begin{eqnarray}&&\displaystyle \frac{H{\left(a\right)}^{2}}{{H}_{0}}\equiv E{\left(a\right)}^{2}=\left[{{\rm{\Omega }}}_{0m}{a}^{-3}+{{\rm{\Omega }}}_{0r}{a}^{-4}+\displaystyle \frac{1}{{T}_{0}}[f-2{Tf}^{\prime} ]\right].\end{eqnarray} \tag{ 24 }$$

The above background evolution recovers the standard model for $\tfrac{1}{{T}_{0}}[f-2{Tf}^{\prime} ]\to {{\rm{\Omega }}}_{{\rm{\Lambda }}}$. These interpretations surely deserve further investigation, which goes beyond the scope of the current analysis.

We have estimated the values of H₀ in the Pantheon sample divided into 3, 4, 20, and 40 redshift bins, where each bin is populated with an equal number of SNe Ia. We have then fitted the values obtained in these bins with a function $g{(z)={H}_{0}(z)={\tilde{H}}_{0}/(1+z)}^{\alpha }$, where the parameter α denotes the evolution of the observables in many astrophysical endeavors (mainly used for studies of GRBs and active galactic nuclei), and ${\tilde{H}}_{0}$ coincides with H₀ at z = 0. We find that there is a slow evolution of H₀ with redshift resulting from the fitting of 3, 4, 20, and 40 bins. A key feature of our data analysis is provided by the evolutionary parameter α, which is of the order of $\sim {10}^{-2}$, and it is compatible with zero within only at least the level of 1.2σ (see Table 1). Although we considered a different number of bins, we obtain the same results for a decreasing ${H}_{0}(z)$ and with the evolutionary parameters α consistent with those cases of 3, 4, 20, and 40 bins. These results make our analysis significantly reliable because it holds in different bins, and thus with a smaller number of SNe Ia in each bin. Furthermore, this evolutionary trend of H₀ was pointed out also in Krishnan et al. (2020a, 2020b) using different probes (Masers, Pantheon, BAO, and CC), a different bin division (six bins), and an alternative fitting function for H₀(z).

Interestingly, if the evolutionary pattern of ${H}_{0}(z)$ is extrapolated at the redshift of the most distant galaxy, z = 11.09, and of the last scattering surface, z = 1100, we obtain a value of ${H}_{0}(z)$ that is compatible within 1σ with the H₀ found by Planck in both the ΛCDM and the ${w}_{0}{w}_{a}$CDM models, thus reducing the H₀ tension, albeit with larger errors.

These results indicate that in the case of a fiducial absolute magnitude value $M=-19.25$ (ΛCDM) or $M=-19.24$ (${w}_{0}{w}_{a}$CDM) such that ${H}_{0}=73.5\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ at z = 0, considering three bins we have reduced the H₀ tension by 54% and 55% for the ΛCDM and ${w}_{0}{w}_{a}$CDM models, respectively. Similar results are obtained with a different number of bins, and they are listed in Table 1; they range in a reduction of the tension of 54% to 72% when considering both the ΛCDM and ${w}_{0}{w}_{a}$CDM models.

Our results not only reduce the so-called H₀ tension, but specifically could highlight an intrinsic evolutionary behavior of ${H}_{0}(z)$: it is no longer a discrepancy between SNe Ia and Planck data, but an effect, in principle, observable at any redshift. If this evolution is not due to the statistical fluctuations of the division in redshift bins and other selection biases or evolution of the color and the stretch that have not been fully accounted for (Nicolas et al. 2020), we show how ${H}_{0}(z)$ could affect the definition of the luminosity itself, especially at high z. More specifically, the new definition of luminosity overestimates the standard distance luminosity definition at z = 11.09 in the ΛCDM model by ≈2%.

Regarding viable theoretical interpretations of this result, the possibility of a slow but monotonic decaying of the H₀ parameter opens interesting theoretical scenarios for its physical or dynamical interpretation. If this result is due to an intrinsic effect associated with the universe expansion, we are led to consider alternative physical and dynamical paradigms for its explanation, outlined in Section 5.2.

The result of our data analysis may indicate that the ΛCDM model needs a possible modification. One viable interpretation consists of regarding ${H}_{0}(z)$ as a consequence of an evolution of the Einstein constant with redshift. This could be an exotic proposal for the late universe dynamics, but a varying Einstein constant is naturally predicted in the so-called f(R) modified theories of gravity, as a result of the presence of an additional scalar field nonminimally coupled to ordinary gravity, as discussed in Section 5.2. This scalar field has frozen dynamics, whose consequence is observable as an effective decreasing of the Einstein (gravity–matter) coupling constant with increasing z. Namely, the observed decaying profile of H₀ with z measures the discrepancy between a real ΛCDM model for the late universe and its representation in terms of an f(R) modified theory of gravity. Our analysis outlines how the dependence of the H₀ value on redshift could arise from hidden astrophysical effects or alternative cosmological models. We have paved the way to a new perspective on the interpretation of the H₀ tension, thus encouraging further data analysis on the physical properties of SNe Ia as standard candles or suggesting the investigation of a new fundamental physics.

In conclusion, regardless of the astrophysical or theoretical interpretation we can provide, we show that our data analysis might lead to a possible solution of the H₀ tension.

This work made use of data supplied by Scolnic et al. (2018) in the GitHub repository https://github.com/dscolnic/Pantheon. We are particularly grateful to D. Scolnic for supporting us with useful suggestions about the estimation of statistical uncertainties. We are grateful to S. Nagataki and T. Hatsuda for discussions and for the suggestions that made possible the beginning of this work. We are particularly grateful to I. Rolf Seitenzahl and F. Roepke for the discussion on the metallicity. We are thankful to A. Lenart, G. Sarracino, and S. Savastano for the support on the cosmological computations.