Evidence of a decreasing trend for the Hubble constant

The value of Hubble constant $H_{0,z}$ is measured by extrapolating the Hubble parameter, $H(z),$ or the luminosity distance, $d_{L}(z)$, from observational data at higher $z$ to $z=0$, by choosing a particular cosmological model. There is no prior knowledge of the redshift evolution of $H_{0,z}$. Similarly to the principal-component approach used to study the redshift evolution of the EoS of dark energy (Huterer & Cooray, 2005), in order to avoid adding some priors on the nature of $H_{0}(z)$, we do not assume that it follows some particular functions. We just allow the value of $H_{0,z}$ in each redshift bin to remain a constant.

Under the assumption of a piece-wise function, $H_{0}(z)$ can be expressed as:

$$H_{0}(z)=\begin{cases}H_{0,z_{1}}&\text{ if }0\leq z<z_{1},\\ H_{0,z_{2}}&\text{ if }z_{1}\leq z<z_{2},\\ \cdots&\cdots,\\ H_{0,z_{i}}&\text{ if }z_{i-1}\leq z<z_{i},\\ \cdots&\cdots,\\ H_{0,z_{N}}&\text{ if }z_{N-1}\leq z<z_{N}.\end{cases}$$

(5)

The parameter $i$ means the $i$th redshift bin, and $N$ is the number of total redshift bins. As described in the previous definition of $H_{0,z}$, we use $H_{0,z_{i}}$ to represent the value of $H_{0}(z)$ between $z_{i-1}$ to $z_{i}$.

A simple way to model the possible evolution of $H_{0,z}$ is through a modification of the standard cosmological model. In the $\Lambda$CDM model, the Hubble parameter is given by

$$H(z)=H_{0}\sqrt{\Omega_{m0}(1+z)^{3}+\Omega_{k0}(1+z)^{2}+\Omega_{\Lambda 0}}.$$

(6)

The parameter $\Omega_{\mathrm{m0}}$ is the cosmic matter density, $\Omega_{\mathrm{k0}}$ is the spatial curvature and $\Omega_{\Lambda 0}$ is the energy density parameter of the cosmological constant. According to the result from Planck CMB measurements (Planck Collaboration, 2020), the space is extremely flat, which gives $\Omega_{\mathrm{k0}}=0$.

The integral form of Equation (6) is

$$\begin{split}H(z)=&H_{0}\sqrt{\Omega_{m0}(1+z)^{3}+\Omega_{\Lambda 0}}\\ =&H_{0}\left(\int_{0}^{z}\frac{3\Omega_{m0}(1+z^{\prime})^{2}}{2\sqrt{\Omega_{m0}(1+z^{\prime})^{3}+\Omega_{\Lambda 0}}}dz^{\prime}+1\right)\\ =&\int_{0}^{z}\frac{H_{0}3\Omega_{m0}(1+z^{\prime})^{2}}{2\sqrt{\Omega_{m0}(1+z^{\prime})^{3}+\Omega_{\Lambda 0}}}dz^{\prime}+H_{0}.\end{split}$$

(7)

The constant 1 is determined by the equation $\Omega_{m0}+\Omega_{\Lambda 0}=1$ when $z=0$. The function $H(z)$ in Equation (7) is identical to that in Equation (6), and it is also more convenient to segment redshift intervals.

Replacing $H_{0}$ in Equation (7) with the piece-wise function in Equation (5), the Hubble parameter is expressed as

$$\begin{split}H(z_{i})=&H_{0,z_{1}}\int_{0}^{z_{1}}\frac{3\Omega_{m0}(1+z)^{2}}{2\sqrt{\Omega_{m0}(1+z)^{3}+\Omega_{\Lambda 0}}}\\ &+H_{0,z_{2}}\int_{z_{1}}^{z_{2}}\frac{3\Omega_{m0}(1+z)^{2}}{2\sqrt{\Omega_{m0}(1+z)^{3}+\Omega_{\Lambda 0}}}\\ &+\cdots\\ &+H_{0,z_{i}}\int_{z_{i-1}}^{z_{i}}\frac{3\Omega_{m0}(1+z)^{2}}{2\sqrt{\Omega_{m0}(1+z)^{3}+\Omega_{\Lambda 0}}}+H_{0,z_{i}}.\end{split}$$

(8)

The last term $H_{0}$ in Equation (7) is replaced by $H_{0}(z)$, so its redshift should be $z_{i}$. Meanwhile, the value of $H_{0}(z)$ at low redshifts is the result of the evolution of high-redshift $H_{0}(z)$. If $H_{0}(z)$ does not exhibit any evolutionary trend, the result will revert to $H_{0}$. We show the value of $H(z)$ calculated by Equations (6) and (8) in Figure 1, respectively. Here, $H_{0}=H_{0,z_{1}}=H_{0,z_{2}}=\cdots=H_{0,z_{i}}=70$ km s^-1 Mpc^-1, $\Omega_{k0}=0,$ and $\Omega_{m0}=0.3$ are assumed. It is clear to see that the value of $H(z)$ calculated via Equation (6) is the same as the one as in Equation (8). Therefore, the expansion of Hubble parameter in Equation (8) is correct.

By assuming the value $H_{0,z}$ is a constant in each redshift bin, previous works have directly used Equation (6) to derive $H_{0,z}$ (Krishnan et al., 2020; Dainotti et al., 2021). First, this approach violates the assumption that the value of $H_{0,z}$ is constant in each bin. Equation (6) actually specifies that the value of $H_{0,z}$ is constant from $z=0$ to $z=z_{i}$ – and not from $z=z_{i-1}$ to $z=z_{i}$. Therefore, the assumption leads $H_{0,z}$ as a constant from $z=0$ to each $z=z_{i}$. During the process of fitting the data in the $i$th redshift bin, the derived value of $H_{0}$ is in the redshift range $0-z_{i}$, not the value of $H_{0,z_{i}}$ in $z_{i-1}-z_{i}$. Their measured $H_{0,z_{i}}$ is not the value in the $i$th redshift bin. Thus, it is problematic to use the result to represent the value of $H_{0,z}$ in the $i$th redshift bin. Second, the correlation among each $H_{0,z_{i}}$ was not considered in previous works. Due to the fact that the Hubble parameter and luminosity distance depend on the summation and the integration over redshifts, the value of $H_{0,z_{i}}$ in low-redshift bins will affect the value in high-redshift bins. In our method, Equation (8) reveals the redshift-evolution of $H_{0}(z)$. The most obvious difference between these two methods is that Equation (6) demonstrates the value of $H_{0,z}$ from $z=0$ to $z=z_{i}$, but Equation (8) reveals the value of $H_{0,z}$ at the $i$th redshift bin. The correlations can be removed by the principal component analysis.

When estimating cosmological parameters, we used the $\chi^{2}$ statistic for a particular model with the parameter set of $\theta$ ($H_{0,z_{i}}$):

$$\chi_{\theta}^{2}=\chi_{H(z)}^{2}+\chi_{BAO}^{2}+\chi_{SNe}^{2},$$

(9)

with

$$\chi_{H(z)}^{2}=\sum_{i=1}^{N}\frac{\left[H_{\mathrm{obs}}\left(z_{i}\right)-H_{\mathrm{th}}\left(z_{i}\right)\right]^{2}}{\sigma^{2}_{i}},$$

(10)

where $H_{obs}(z_{i})$ and $\sigma_{i}$ are the observed Hubble parameter and the corresponding 1$\sigma$ error. Then, $H_{th}(z_{i})$ is the value calculated from Equation (8).

The value of $\chi^{2}_{BAO}$ is

$$\chi_{BAO}^{2}=\left[\nu_{\mathrm{obs}}\left(z_{i}\right)-\nu_{\mathrm{th}}\left(z_{i}\right)\right]\mathbf{C_{BAO}}^{-1}\left[\nu_{\mathrm{obs}}\left(z_{i}\right)-\nu_{\mathrm{th}}\left(z_{i}\right)\right]^{T},$$

(11)

where $\nu_{obs}$ is the vector of the BAO measurements at each redshift $z$ (i.e. $D_{V}\left(r_{s,{\rm fid}}/r_{s}\right),D_{M}/r_{s},D_{H}/r_{s},D_{A}/r_{s}$). The comoving distance $D_{M}(z)$ is:

$$D_{M}(z)=(1+z)D_{A}(z),$$

(12)

where $D_{A}(z)$ is the angular size distance. The radial BAO projection $D_{H}(z)$ is:

$$D_{H}(z)=\frac{c}{H(z)}.$$

(13)

The angle-averaged distance $D_{V}(z)$ is

$$D_{V}(z)=\left[czD^{2}_{M}(z)/H(z)\right]^{1/3}.$$

(14)

The sound horizon of the fiducial model is $r_{s,fid}=147.5$ Mpc.

The full covariance matrix of the Pantheon+ SNe Ia sample is considered in the process of calculating $\chi_{SNe}^{2}$ (Scolnic et al., 2022). This covariance 1701$\times$1701 matrix $\mathbf{C_{SNe}}$ is defined as

$$\mathbf{C_{SNe}}=\mathbf{D_{stat}}+\mathbf{C_{sys}},$$

(15)

where $\mathbf{D_{stat}}$ is the statistical matrix and $\mathbf{C_{sys}}$ is the systematic covariance matrix . The value of $\chi^{2}_{SNe}$ is:

$$\chi_{SNe}^{2}=\left[\mu_{\mathrm{obs}}\left(z_{i}\right)-\mu_{\mathrm{th}}\left(z_{i}\right)\right]\mathbf{C_{SNe}}^{-1}\left[\mu_{\mathrm{obs}}\left(z_{i}\right)-\mu_{\mathrm{th}}\left(z_{i}\right)\right]^{T}.$$

(16)

The parameter $\mu_{\mathrm{obs}}(z_{i})$ is the distance module from the Pantheon+ sample. The theoretical distance modulus, $\mu_{\textrm{th}}$, is defined as

$$\mu_{th}=5\log_{10}d_{\mathrm{L}}+25,$$

(17)

Taking $H_{th}(z)$ from Equation (8), the luminosity distance, $d_{L}$, is expressed as

$$d_{L}(z)=c(1+z)\int_{0}^{z}\frac{\mathrm{d}z^{\prime}}{H_{th}(z^{\prime})}.$$

(18)

The constraints are derived by the Markov Chain Monte Carlo (MCMC) code emcee (Foreman-Mackey et al., 2013). The adopted prior is $H_{0}\in$ [50,80] $\textrm{km}\leavevmode\nobreak\ \textrm{s}^{-1}\textrm{Mpc}^{-1}$ for all $H_{0,z_{i}}$. Observations of SNe Ia and CMB set tight constraints on the cosmic matter density $\Omega_{m0}$, which are both around $0.3$ (Brout et al., 2022; Planck Collaboration, 2020). Thus, a fiducial value $\Omega_{m0}=0.3$ for cosmic matter density was used during the fitting process. We repeated our analysis with the Gaussian prior $\Omega_{m0}=0.315\pm 0.007$ from Planck CMB measurements (Planck Collaboration, 2020) and found that it gives very similar results. Thus, the results are largely insensitive to the choice of either prior.

Although each redshift bin is treated as independent during the fitting process, the results of $H_{0,z_{i}}$ are still correlated. This is expected, as the integration and summation over low-redshift bins in Equation (8) definitely affects the model fit in the middle and higher redshift bins. In order to remove these correlations, we calculate the transformation matrix (Huterer & Cooray, 2005). The principal component analysis presents a compressed form of the result with all information about the constraint from observations.

The fitting results impose constraints on the parameters $H_{0,z_{i}}(i=1...N)$. The covariance matrix can be generated by taking the average over the chain, namely,

$$\mathbf{C}=\langle\mathbf{H}\mathbf{H}^{\rm T}\rangle-\langle\mathbf{H}\rangle\langle\mathbf{H}^{\rm T}\rangle,$$

(19)

where $\mathbf{H}$ is a vector with components $H_{0,z_{i}}$ and $\mathbf{H}^{T}$ is the transpose. We transformed the covariance matrix to decorrelate the $H_{0,z_{i}}$ estimations (Huterer & Cooray, 2005). This was achieved by finding the uncorrelated basis by diagonalizing the inverse covariance matrix. The Fisher matrix is defined as

$$\mathbf{F}\equiv\mathbf{C}^{-1}\equiv\mathbf{O}^{\mathrm{T}}{\Lambda}\mathbf{O},$$

(20)

where the matrix $\Lambda$ is the diagonalized covariance for transformed bins. The transformation matrix $\mathbf{T}$ is used as the square root of the Fisher matrix, which is defined as

$$\mathbf{T}=\mathbf{O}^{\mathrm{T}}{\Lambda}^{\frac{1}{2}}\mathbf{O}.$$

(21)

One advantage of this method is that the rows of $\mathbf{T}$ are almost positive across all bands. After being normalized, the rows of $\mathbf{T}$ which means the weights for $H_{0,z_{i}}$ sum to unity. Finally, it is clear that the new parameters,

$$\tilde{\mathbf{H}}=\mathbf{TH,}$$

(22)

are uncorrelated, because they have the covariance matrix ${\Lambda}^{-1}$.

After the correlations among $H_{0,z}$ are removed, the corner plots of $H_{0,z}$ are shown in supplementary Figures 5 and 6 for the equal-number binning case and equal-width binning case, respectively. The constraints on $H_{0,z}$ are tight and the probability distributions are almost Gaussian.

Before the final result, we estimated the possible evolutionary trend only with the Pantheon+ sample. The whole sample contains 1701 data, most of which are located at low redshift. Therefore, we chose six bins with the upper boundaries of $z_{i}$ = 0.1, 0.2, 0.3, 0.6, 1.0, 2.4. On account of the relatively few data points in high-redshift bins, some intervals have to be loose. The value of $H_{0,z}$ as a function of redshift (panel a) and the normalized probability distributions (panel b) are shown in Figure 2. The decreasing trend is clear at $z>0.3$, but the scarce data in high-redshift bins give poor constraints. The 1 $\sigma$ uncertainties for high-redshift bins are so large that we decided to add the observed Hubble parameter data and BAO data (Yu et al., 2018; Cao & Ratra, 2022). The joint dataset gives a strict constraint and can be used to research the evolution by different binning methods.

Two redshift-binning methods were adopted. The first one is equal-number method, namely, the number of data points in each bin is almost equal (Dainotti et al., 2021). In contrast, the second one is equal-width method, namely, the bins are equally spaced in redshift (Huterer & Cooray, 2005). For the equal-number method, we chose nine bins with the upper boundaries of $z_{i}$ = 0.0122, 0.025, 0.037, 0.108, 0.199, 0.267, 0.350, 0.530, and 2.40. Each bin contains about 194 data points (listed in Table 3). More binnings were also performed and we find that the likelihood distributions of $H_{0,z}$ for some bins deviate from Gaussian, due to the scarcity of data points in some high-redshift bins. In the case of non-Gaussian $H_{0,z}$ distributions, the decorrelation process may introduce some biases. Therefore, nine bins were adopted.

Altogether, there are nine free parameters in the Markov chain Monte Carlo (MCMC) approach, namely, $H_{0,z_{1}},...,H_{0,z_{9}}$. The Pantheon+ SNe Ia sample (Scolnic et al., 2022), BAO data and observed Hubble parameter data were used to estimate the value of $H_{0,z}$ (see Sect. 3). Given this joint dataset, we used the MCMC code emcee (Foreman-Mackey et al., 2013) to sample the parameter set $\theta$ ($H_{0,z_{i}}$). The prior of $H_{0}\in$ [50,80] $\textrm{km}\leavevmode\nobreak\ \textrm{s}^{-1}\textrm{Mpc}^{-1}$ for all $H_{0,z_{i}}$ was adopted. After removing the correlation among $H_{0,z}$ by diagonalizing the covariance matrix (see Methods), the best-fit results are shown in Table 3. The value of $H_{0,z}$ as a function of redshift (panel a) and the normalized probability distributions (panel b) are shown in Figure 3. Due to the large number of data in each bin and a fixed $\Omega_{m0}$, the uncertainty of $H_{0,z}$ is less than 1.0, which is consistent with previous work (Dainotti et al., 2021; Wang, 2022b). In the first seven bins, the fitting results are nearly constant and the value is consistent with the one derived by the local distance ladder (Riess et al., 2022), however, it drops rapidly thereafter. It is worth noting that the result in the last bin is consistent with the value from Planck measurements at a 2$\sigma$ confidence level (Planck Collaboration, 2020).

There is an apparent decreasing trend in $H_{0}(z)$ as a function of redshift. To quantify the significance of this trend, we used the null hypothesis method (Wong et al., 2020; Millon et al., 2020). The hypothesis in this situation posits that the values of $H_{0}(z)$ in each redshift bin are consistent with each other. We first fit a linear regression through each redshift bin. Next, we generated sets of nine mock $H_{0}$ values with their own uncertainty probability distribution centered around the value of $H_{0}=73.04\pm 1.04$ km s^-1 Mpc^-1 (Riess et al., 2022). The weight of each mock value was calculated as follows (Wong et al., 2020). First, the uncertainties’ probability distributions are rescaled so that their maximal values are equal to 1. Then, the area under each distribution is calculated and the areas are rescaled by their median. Last, the inverse square of the rescaled areas is taken as the weight for each mock value. We also fit a linear regression through the mock value. The slope of the data falls 3.8$\sigma$ away from the mock slope distribution. In other words, the decreasing trend in $H_{0}(z)$ with increasing redshift has a significance of 3.8$\sigma$.

In the second case, the bins are equally spaced with a redshift-width 0.1 in the redshift range [0, 0.4]. Because there are so few data points from $z=0.4$ to $z=2.4$, wider intervals should be adopted in this range. We also tried to equally bin $z$ with a width of $0.1$ in the redshift range [0, 1.0]. The poor constraints in some intervals lead to biases in the decorrelation process. Finally, ten bins are adopted with the upper boundaries $z_{i}$ = 0.10, 0.20, 0.30, 0.40, 0.6, 0.8, 1.1, 1.5, 2.0, and 2.4. The number of data points and best-fit results are given in Table 4. The uncertainties of the constraints increase gradually as the number of data decreases. The value of $H_{0,z}$ as a function of redshift (panel a) and the normalized probability distributions of $H_{0,z}$ (panel b) are given in Figure 4. The fitting results of the first three bins are consistent with the value from the local distance ladder within 1$\sigma$ confidence level (Riess et al., 2022), and the last two bins are consistent with the value from Planck CMB measurements within 1$\sigma$ confidence level (Planck Collaboration, 2020). There is a gradually decreasing trend in $H_{0}(z)$ from the third bin to the last bin. Using the same method as above, the significance of the decreasing trend is 5.6$\sigma$. More importantly, the decreasing trend begins at the same redshift $z\sim 0.3$ for the two binning methods.

Finally, we tested whether additional parameters to describe $H_{0}(z)$ are actually needed to improve on a flat $\Lambda$CDM model fit to the data. For the flat $\Lambda$CDM model, the same value of $\Omega_{m0}=0.3$ was adopted and $H_{0}$ was left as a free parameter. Two kinds of standard information criteria were considered: Akaike information criterion (AIC) (Akaike, 1974) and the Bayesian information criterion (BIC) (Schwarz, 1978). Their definitions are: $\mathrm{AIC}=-2\ln L+2k,\leavevmode\nobreak\ \mathrm{and}\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{BIC}=-2\ln L+k\ln n,$ where $L\propto e^{-\chi^{2}/2}$ is the value of the maximum likelihood function, $k$ is the number of model parameters, and $n=1746$ is the total number of data points. From Table 5, there is a significant improvement for both binning methods relative to the flat $\Lambda$CDM model, with $\Delta$AIC of $-44.49$ and $-37.41$, $\Delta$BIC of $-0.77$ and $-11.77$ for equal-number and equal-width binning methods, respectively. Thus, these values are sufficient to favor the $H_{0}(z)$ model over $\Lambda$CDM.