Measurements of the Hubble Constant with a Two-rung Distance Ladder: Two Out of Three Ain't Bad

The current rate of cosmic expansion, identified as the Hubble constant H₀, is a fundamental cosmological parameter which sets the age and scale of the universe. In recent times, the unexplained discrepancy between local measures of H₀ based on local distance indicators (e.g., Riess et al. 2022) and predictions based on measurements of the early universe (e.g., Planck Collaboration et al. 2018) has attracted attention as a "Hubble tension" (Di Valentino et al. 2021). Predictions based on the early universe appear robust, producing the same value at high precision when assuming the Lambda cold dark matter (ΛCDM) model with or without the use of cosmic microwave background (CMB) anisotropy data (Addison et al. 2018). Direct measurements of the Hubble constant at low redshift have been conducted using a variety of distance indicators, including water mega-masers (Pesce et al. 2020), type II supernovae (de Jaeger et al. 2022), surface-brightness fluctuations (SBF; Blakeslee et al. 2021), Tully–Fisher relations (Kourkchi et al. 2020; Schombert et al. 2020), and modeling of gravitationally lensed quasars (Birrer et al. 2020; Wong et al. 2020). However, the most precise measurements in the literature have come from the use of a three-rung distance ladder using Type Ia supernovae (SNe Ia; Riess et al. 2022; Freedman 2021).

Since ∼2005, the SH0ES team (Riess et al. 2005, 2022) has been engaged in a project to refurbish the distance ladder measurement of H₀ using the Wide Field Camera 3 (WFC3) and Advanced Camera for Surveys (ACS) instruments on board the Hubble Space Telescope (HST). These measurements have been based on a three-rung distance ladder formed by (1) geometric measurements of the distance to Cepheid variable stars, (2) Cepheid distance measurements of SN Ia host galaxies, and (3) standardized SNe Ia magnitudes found in the Hubble flow between 0.023 < z < 0.15. This redshift range was chosen to reduce the dependence of the result on peculiar motions at low redshifts and cosmological effects at higher redshifts, while retaining sufficient statistics to constrain the Hubble flow. Hallmarks of the SH0ES project include the use of the same instrument to measure all Cepheids across the ladder to nullify systematic uncertainties, the use of near-infrared (NIR) and three colors to mitigate dust, and accounting for the dependence of Cepheid luminosity on locally measured abundances. This program also extended the range of Cepheid parallax measurements to a few kiloparsecs using spatial scanning of the HST (Riess et al. 2014; Casertano et al. 2016). The SH0ES team has expanded the sample of cosmological SNe Ia calibrated by Cepheids from a few in the Hubble Key Project (Freedman et al. 2001) with high-quality data to 42 in 37 host galaxies in the most recent result (Riess et al. 2022, hereafter R21).

In this work we present a two-rung distance ladder based on geometrically calibrated Cepheid distances at z < 0.01 as a measure of the Hubble flow. In the conventional three-rung distance ladder, extragalactic distances measured by a geometrically calibrated distance indicator (e.g., Cepheid or Mira variables, the tip of the red giant branch) are colocated with one or more SNe Ia to calibrate the SN Ia absolute luminosity. However, the redshifts of these galaxies are not involved in the determination of H₀. Instead, the Hubble flow is constrained by a separate sample of SNe Ia, because their bright absolute magnitude (M_B ∼ −19) allows them to be detected well into the Hubble flow. By using objects at redshifts cz ⪆ 7000 km s⁻¹, the fractional effect of galactic peculiar velocities (∼300 km s⁻¹) will be smaller than SN Ia uncertainties. Wu & Huterer (2017) found that, given the Riess et al. (2016) SN Ia sample, (uncorrected) cosmic variance would affect the SH0ES H₀ measurement by only ∼0.5%. As SNe in calibrator and Hubble flow samples are largely observed by the same telescopes and discovered by similar surveys, many calibration and demographic biases are expected to largely cancel out from the differential measurement, and are further reduced by recalibration using wide-angle surveys (Scolnic et al. 2015; Brownsberger et al. 2021). Semi-independent cross-checks of the third rung of the distance ladder include Burns et al. (2018), Jones et al. (2022), and Dhawan et al. (2022). However the astrophysics of SNe Ia are debated, including the roles of progenitor age (Jones et al. 2018; Rigault et al. 2020; Rose et al. 2021) and host galaxy dust (Brout & Scolnic 2021; González-Gaitán et al. 2021; Thorp et al. 2021; Popovic et al. 2021) in the diversity of observed populations. SN-independent measures are thus a crucial check of possible systematic biases in the SH0ES result.

In this work we remove the third rung, measuring H₀ by directly comparing Cepheid distance measurements to their host galaxy redshifts. Doing so will give us a measurement of H₀ independent of any potential SNe Ia supernova systematics. However, by removing the third rung we are forced to reduce the median redshift of the Hubble flow sample by an order of magnitude from z ∼ 0.04 to z ∼ 0.006 and the maximum redshift from z = 0.15 to z = 0.011. At such low redshifts, peculiar velocities (PVs) are expected to have a magnitude ∼20% of the observed redshifts. Further, the volume probed by these Cepheid hosts is not large compared to cosmological structures on the scale of baryon acoustic oscillations (∼140 Mpc), and as a result the effects of large-scale structure or "cosmic variance" on the measurement must be addressed. The reduction in median redshift by a factor of 6 can be naively expected to increase the effect of PVs on the measurement by at least the same amount, corresponding to an expected uncertainty in H₀ greater than 3%.

Cosmological measurements using SNe Ia as standard candles have made use of PV corrections based on models of the density field in the nearby universe to mitigate systematic PV effects from structure on scales greater than ∼10 Mpc h⁻¹, where h = H₀/100 km s⁻¹ Mpc⁻¹. Redshifts used in the cosmological analysis of Brout et al. (2022) and Riess et al. (2022) also made use of "group" corrections (the replacement of an SN Ia host redshift with a mean redshift of associated galaxies), which attempt to remove noise on small scales from virial motions of galaxies within groups or clusters. Peterson et al. (2021) conclusively demonstrated at high significance that the dispersion of the SN Ia Hubble diagram is reduced with the application of both kinds of PV corrections. In particular two density field reconstructions (Carrick et al. 2015; Lilow & Nusser 2021) performed excellently on these tests, reducing the scatter of the Hubble diagram by ∼10% in a sample of 584 SNe at z < 0.08. ⁸ Both reconstructions are based on galaxy counts and luminosities from redshift surveys of the nearby universe, covering the local universe within ∼200 Mpc h⁻¹ of the Milky Way.

In order to make a precise constraint on H₀ from the Cepheid distances and Pantheon+ redshifts, it will be essential to make use of these reconstructions to extract as much information about the structure of the nearby universe as is possible. To ensure our constraint is accurate we must also quantify any remaining systematic uncertainties in the corrected redshifts, including the covariance between host galaxies. Finally, as the velocity noise relative to the distance noise increases at low redshifts, selection effects become increasingly important; the well-known volumetric or Malmquist bias increases as the square of the relative error. At distances where PV effects make up ∼20% of the signal of the Hubble flow, it will be necessary to evaluate the selection of the SH0ES Cepheid sample as fair markers of the Hubble flow and ensure that any associated uncertainties are included in our error estimates.

In this work, we make use of distances to Cepheid-hosting galaxies measured by SH0ES and Pantheon+ redshifts, in conjunction with two density field reconstructions of the nearby universe to measure the Hubble constant. In Section 2 we discuss the data we make use of, including Cepheid data, redshifts, and PV reconstructions. In Section 3 we present our analysis methodology, which uses a hierarchical Bayesian model to address uncertainties and biases associated with a low-redshift measurement; this section includes estimates of the correlations in predicted redshifts based on modeling and the limitations of the reconstructions, as well as modeling of the SH0ES host selection. We discuss the implementation of this model in Section 4 and provide validations of the methodology on the SH0ES data. We present estimates of H₀ from the two-rung distance ladder in Section 5, and discuss our findings in Section 6.

2.1. Cepheid Distances

2.2. Galaxy Redshifts and Peculiar Velocities

Because uncorrected PVs represent a dominant source of uncertainty, and redshift or distance limits imposed by the sample selection may produce a bias, we proceed with a more extensive effort to model the data than is typically used for a three-rung distance ladder where similar effects are negligible.

We construct a hierarchical Bayesian model for the observed distance moduli of the Cepheid galaxies given the reconstructed density/velocity fields, observed redshifts, and observed positions on the sky. We model our data as determined by the Hubble constant and a set of nuisance parameters, which includes the true comoving distance to each observed galaxy as a latent parameter. This formulation of the problem allows us to easily account for aspects of the analysis such as selection, while budgeting for many of the systematic uncertainties of our analysis. In Section 3.1 we present our account of the observed data in terms of the model. Our treatment of the C15 and LN21 reconstructions are presented in Sections 3.1.1 and 3.1.2. We discuss the priors we use in Section 3.2, and the modeling for systematic effects of selection in Section 3.2.1. We summarize our model in Section 3.3.

3.1. Data Likelihood

We first define the relation between the modeled parameters and the observed data. To track the relationships between distances, redshifts, and PVs, we define the parameter vector of latent comoving distances in Mpc h⁻¹ units for the 35 objects of our Cepheid galaxy sample χ . In a universe described by the Friedman–Lemaître–Robertson–Walker metric, the cosmological redshift for the ith object is

$$\begin{eqnarray}&&\begin{array}{l}{z}_{\mathrm{Cosm}.}({{\boldsymbol{\chi }}}_{i})=100\cdot {{\boldsymbol{\chi }}}_{i}/c\times \left(1+\displaystyle \frac{{q}_{0}+1}{2}\cdot 100\cdot {{\boldsymbol{\chi }}}_{i}/c\right.\\ \left.+\,\displaystyle \frac{{j}_{0}+2{q}_{0}+1}{6}{\left(100\cdot {{\boldsymbol{\chi }}}_{i}/c\right)}^{2}+O\left({{\boldsymbol{\chi }}}_{i}^{3}\right)\right)\end{array}\end{eqnarray} \tag{ 1 }$$

and, concurrently, the distance modulus is

$$\begin{eqnarray}&&\begin{array}{l}\mu ({{\boldsymbol{\chi }}}_{i},{H}_{0})=5{\mathrm{log}}_{10}\left(\displaystyle \frac{{{cz}}_{\mathrm{Cosm}.}({{\boldsymbol{\chi }}}_{i})}{{H}_{0}}\times \left(1+\displaystyle \frac{1-{q}_{0}}{2}{z}_{\mathrm{Cosm}.}({{\boldsymbol{\chi }}}_{i})\right.\right.\\ -\,\left.\left.\displaystyle \frac{1-{q}_{0}-3{q}_{0}^{2}+{j}_{0}}{6}{z}_{\mathrm{Cosm}.}{\left({{\boldsymbol{\chi }}}_{i}\right)}^{2}+O({{\boldsymbol{\chi }}}_{i}^{3})\right)\right)+25,\end{array}\end{eqnarray} \tag{ 2 }$$

with the two kinematic parameters fixed to q₀ = − 0.55 and j₀ = 1, and H₀ a free parameter. Our result is insensitive to changes in the kinematic quantities, as the effect of these corrections to the Hubble law is ∼0.5% for the most distant objects in our sample, and much less for the majority.

The likelihood of the observed distance moduli μ _Ceph with uncertainties Σ_μ is simply defined by the multivariate normal distribution:

$$\begin{eqnarray}&&{{\boldsymbol{\mu }}}_{\mathrm{Ceph}}| {\boldsymbol{\chi }},\mathrm{log}{H}_{0}\sim { \mathcal N }(\mu ({\boldsymbol{\chi }},{H}_{0}),{{\boldsymbol{\Sigma }}}_{\mu }).\end{eqnarray} \tag{ 3 }$$

This term of the model is the component of the likelihood in which H₀ appears.

The observed likelihood of the redshifts is more complicated. The redshift addition formula gives the observed redshift corrected to the CMB frame as z_CMB + 1 = (1 + z_Pec.) × (1 + z_Cosm.). From this the likelihood of the observed vector of CMB-frame, group-corrected redshifts z _CMB are given by the multivariate normal distribution

$$\begin{eqnarray}&&\displaystyle \frac{1+{{\boldsymbol{z}}}_{\mathrm{CMB}}}{1+{z}_{\mathrm{Cosm}.}({\boldsymbol{\chi }})}| {\boldsymbol{\chi }},{\boldsymbol{\theta }}\sim { \mathcal N }(1+{v}_{\mathrm{Pred}.}({\boldsymbol{\chi }},{\boldsymbol{\theta }})/c,{{\boldsymbol{\Sigma }}}_{\mathrm{Pec}.}({\boldsymbol{\chi }};{\boldsymbol{\theta }})/{c}^{2}),\end{eqnarray} \tag{ 4 }$$

where the division on the left-hand side is conducted element-wise. v_Pred.( χ , θ ) is the predicted PV from the reconstruction used as discussed in Section 3.1.1, and Σ_Pec.( χ ; θ ) is our estimated uncertainty in these predictions as presented in Section 3.1.2. The parameters θ that define the velocity predictions and uncertainties will be introduced in the next two sections.

These two likelihood definitions make up the constraints on the model parameters from the data. H₀ is measured as the relative scale between the Cepheid distances and the redshifts.

3.1.1. Peculiar Velocity Predictions

The two PV reconstructions are publicly distributed as data cubes, which give the values of the density contrast $\delta ({\boldsymbol{r}})=(\rho ({\boldsymbol{r}})-\bar{\rho })/\bar{\rho }$ (where ρ is the density of matter) and PVs v _Pec.( r ) at points on a densely sampled grid. These grids are given at a resolution higher than the smoothing scale. These can easily be interpolated to give functions for the density contrast and velocity as a function of position in comoving galactic coordinates δ_Rec.( r ), v _Rec.( r ), respectively.

The two reconstructions also have both made measurements of the parameters that relate the density to observed velocities. We quantify the effect of the stated uncertainties by marginalizing over the values given in the original work. The relevant parameters are as follows:

• 1.  
  The scale parameter β, the ratio of the logarithmic growth rate f(Ω_M ) of density perturbations to the galaxy bias b. This relates predicted PV to galaxy overdensities.
• 2.  
  The external bulk flow V _Ext, a velocity vector added to every point within the reconstructed volume. This models any contribution to PVs sourced by perturbations on scales larger than the volume probed by the data.

The external bulk flow is parameterized by V_Ext, the magnitude in units of kilometers per second, and (l, b) the galactic longitude and latitude of its direction. Note that this dipole is distinct from the ∼600 km s⁻¹ motion of the Milky Way with respect to the CMB, most of which is due to gravitational effects from the local volume (within 200 Mpc h⁻¹). We remove the fiducial values of the quantities {β_Fid., V_Fid.Ext., l_Fid., b_Fid.} from the calculated PVs to define the unscaled local velocities:

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{\mathrm{Rec}.}({\boldsymbol{r}})=\displaystyle \frac{{{\boldsymbol{v}}}_{\mathrm{Rec}.}({\boldsymbol{r}})-{V}_{\mathrm{Fid}.\mathrm{Ext}.}\times \hat{n}({l}_{\mathrm{Fid}.},{b}_{\mathrm{Fid}.})}{{\beta }_{\mathrm{Fid}.}},\end{eqnarray} \tag{ 5 }$$

where $\hat{n}(l,b)$ is the unit position vector as a function of angular coordinate. We can then compute the predicted PV for the ith object with unit position vector ${\hat{n}}_{i}$:

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\mathrm{Pred}.}({{\boldsymbol{\chi }}}_{i},\beta ,{V}_{\mathrm{Ext}.},l,b) & = & {\hat{n}}_{i}\cdot \left[\beta {{\boldsymbol{u}}}_{\mathrm{Rec}.}({{\boldsymbol{\chi }}}_{i}\times {\hat{n}}_{i})\right.\\ & & +\left.{V}_{\mathrm{Ext}.}\times \hat{n}(l,b))\right].\end{array}\end{eqnarray} \tag{ 6 }$$

3.1.2. Covariance of Peculiar Velocity Reconstructions

Previous work, including Riess et al. (2016), Scolnic et al. (2018), Boruah et al. (2020) and others, have analyzed SN Ia distances while accounting for PV uncertainties using an assumption of diagonal scatter in corrected redshifts, corresponding to an assumption that all correlations in the data have been removed. However, the smoothing at the scale of 5 Mpc h⁻¹ of the density fields used here removes any information at smaller distances. Any two galaxies closer together than this scale will have similar predictions for their density/PV, while their true velocities may include unmodelled nonlinear or virial noise. Further, they make use of a finite number of tracers (galaxies) of the underlying density. Lastly, the two reconstructions show some differences in predicted densities and velocities at the scale of ∼60 km s⁻¹, with some differences apparent across the entire volume. It is thus apparent that use of these reconstructions to remove predicted velocities will ameliorate but not totally remove correlations in the data, and a trustworthy estimate of the Hubble constant from data with median distance ∼20 Mpc h⁻¹ will require an estimate of these correlations. Our goal then is to construct an appropriate covariance matrix Σ_Pec. for the predicted PVs. These correlated uncertainties will ultimately contribute ∼2% of the final error budget (see Table 1).

Table 1. Two-rung Simple Estimator and Cumulative Analysis Changes

	H0 − 72.8 km s−1 Mpc−1	σ H0	ΔH0 from last
Simple estimator	−2.5	1.7	⋯
Include Cepheid systematics	−3.1	1.9	−0.7
Include PV corrections	−0.5	2	2.6
Include PV covariance estimate	1.9	2.8	2.4
Include selection effects (distance limited) (see Row 7, Table 2)	0	2.7	−1.9
Simple estimator, cz > 2000 km s−1	−0.6	2.1	⋯

Note. Effect of components of the model, as discussed in Section 4.2, on the H₀ measurement. Except for the last row, we cumulatively add each element of the model to show the impact of each element on the error estimate and central value of H₀. The full model, as shown in the second to last row, is used to produce the full results in Row 7 of Table 2; it uses the C15 reconstruction with assumed distance-limited selection and fixed-PV amplitude parameters. The last row, given for reference, shows the simple estimator when only galaxies with cz > 2000 km s⁻¹ are used.

Download table as:  ASCIITypeset image

To understand the uncertainties in the reconstructed density fields, we first need to understand the scales and amplitudes of the underlying cosmological structures which these attempt to model. If the density field under ΛCDM is approximated as Gaussian, the statistics of the density contrast δ are fully specified by the power spectrum P(k). The power spectrum is defined by the relation

$$\begin{eqnarray}&&\langle \tilde{\delta }({{\boldsymbol{k}}}^{* })\tilde{\delta }({\boldsymbol{k}})\rangle ={\left(2\pi \right)}^{3}P(| {\boldsymbol{k}}| ){\delta }^{3}({\boldsymbol{k}}-{{\boldsymbol{k}}}^{* }),\end{eqnarray} \tag{ 7 }$$

where δ³ is a Dirac delta function and $\tilde{\delta }(k)$ is the Fourier transform of the density contrast. Given the linear approximation that the divergence of the velocity field is proportional to the density and that the velocity is vorticity free, the covariance of the PVs along the line of sight at comoving distances χ_i and χ_j separated by angle θ is given by

$$\begin{eqnarray}\begin{array}{rcl}\langle {v}_{i}{v}_{j}\rangle & = & {\left.\displaystyle \frac{{dD}}{d\tau }\right|}_{{\chi }_{i}}{\left.\displaystyle \frac{{dD}}{d\tau }\right|}_{{\chi }_{j}}\displaystyle \int {dk}/(2{\pi }^{2})\ {P}_{{\rm{\Lambda }}\mathrm{CDM}}(k)\\ & & \times F(k{\chi }_{i},k{\chi }_{j},\cos (\theta )),\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\begin{array}{l}F(u,v,\cos (\theta ))\\ \ \ =\,\left.\displaystyle \frac{{\partial }^{2}\mathrm{sinc}(\sqrt{{x}^{2}-2{xy}\cos (\theta )+{y}^{2}})}{\partial x\partial y}\right|x=u,y=v,\end{array}\end{eqnarray} \tag{ 9 }$$

where τ is conformal time, and D is the linear, scale-independent "growth factor" that describes the amplitude of density perturbations as a function of cosmic time (Davis et al. 2011; Huterer et al. 2015). The term we write as $F(u,v,\cos (\theta ))$ can be evaluated either as a series over spherical Bessell functions or constructed from elementary and trigonometric functions. We evaluate the power spectrum using the CLASS code (Blas et al. 2011) with nonlinear corrections from halofit (Takahashi et al. 2012) using cosmological parameters from Planck (Planck Collaboration et al. 2018). As F(u, u, 1) = 1/3, the predicted diagonal dispersion $\langle {v}_{i}^{2}\rangle$ at z = 0 is proportional to the integral ∫dkP(k). Evaluating $\langle {v}_{i}^{2}\rangle$, we find a predicted ${\left(380\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}$ of PV scatter in the local universe. This quantity represents the uncertainty in uncorrected CMB-frame redshifts as a distance indicator.

To budget for PV uncertainties in redshifts after correction using a reconstruction, we derive two independent sources of remaining noise, which we label ${P}_{\mathrm{Unmodelled}}(k),{P}_{\mathrm{Structure}}(k)$. These quantities are derived and defined in Appendix B. P_Unmodelled(k) represents uncertainties resulting from the Gaussian smoothing applied in the reconstruction procedures, and is predicted to contribute 264 km s⁻¹ of diagonal scatter. This term also accounts for the increased homogeneity of the maps relative to the true density fields (Berlind et al. 2000) by representing the additional variance on scales below the smoothing scale. The size of this prediction is uncertain, as it is dependent on theory/simulation-based accounts of nonlinear gravitational effects. ${P}_{\mathrm{Structure}}(k)$ is empirically estimated from the differences between the results of C15 and LN21, and represents our estimate of the systematic uncertainty in their modeling of large-scale structure. This component contributes 78 km s⁻¹ of diagonal scatter, which is present across larger distances on the sky. The sum of these two components, which we label ${\hat{P}}_{{\rm{\Delta }}}(k)$, is the total power remaining in the density (and thus velocity) field after correction. We compare ${\hat{P}}_{{\rm{\Delta }}}(k)$ to the ΛCDM prediction in Figure 1; the difference between the two represents the uncertainty removed from our measurement by the use of the galaxy catalog information. Accordingly, our fiducial error budget predicts that use of the corrections will reduce the diagonal scatter of the sample from ∼${\left(380\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}$ to ∼${\left(270\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}$.

Zoom In Zoom Out Reset image size

To account for uncertainties in the size of each of the two effects discussed, we introduce two additional parameters, ${\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})},{\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}$, which scale the amplitude of the two components of the budgeted uncertainties:

$$\begin{eqnarray}&&\begin{array}{l}{\hat{P}}_{{\rm{\Delta }}}(k;{\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})},{\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})})={\left({\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})}\right)}^{2}\\ \ \ \times \,{P}_{\mathrm{Unmodelled}}^{(\mathrm{Normalized})}(k)+{\left({\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}\right)}^{2}{P}_{\mathrm{Structure}}^{(\mathrm{Normalized})}(k),\end{array}\end{eqnarray} \tag{ 10 }$$

where each component of the power spectrum has been normalized by the fiducial scatter $\langle {v}_{i}^{2}\rangle$ so the value of each parameter is equal to the scatter contributed by each component. A value of ${\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})}$ greater than the fiducial value of 264 km s⁻¹ would indicate large virial or nonlinear motions on scales too small to be detected by the reconstructions; a value of ${\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}$ larger than the fiducial value of 78 km s⁻¹ would indicate that our estimates of the error in the reconstructed fields are understated.

In two model variants (results shown in the first two rows of Table 2) we attempt to fit the data without the use of external galaxy catalog information. In these cases, we instead use the relation

$$\begin{eqnarray}&&\begin{array}{l}{\hat{P}}_{{\rm{\Delta }}}(k;{\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})},{\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})})={\left({\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})}\right)}^{2}\\ \ \ \times \,{P}_{\mathrm{Unmodelled}}^{(\mathrm{Normalized})}(k)+{\left({\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}\right)}^{2}{P}_{{\rm{\Lambda }}\mathrm{CDM}\,\mathrm{Structure}}^{(\mathrm{Normalized})}(k),\end{array}\end{eqnarray} \tag{ 11 }$$

where ${P}_{{\rm{\Lambda }}\mathrm{CDM}\,\mathrm{Structure}}(k)$ is defined as

$$\begin{eqnarray}&&{P}_{{\rm{\Lambda }}\mathrm{CDM}\ \mathrm{Structure}}(k)={P}_{{\rm{\Lambda }}\mathrm{CDM}}(k)-{P}_{\mathrm{Unmodelled}}(k).\end{eqnarray} \tag{ 12 }$$

This reduces to the full ΛCDM power spectrum in the case when the amplitude parameters ${\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})},{\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}$ are fixed to their fiducial values, but allows the measurement to marginalize over unknown nonlinearities and cosmology.

Table 2. Final Results: Variations on Selection and Peculiar Velocity Assumptions

Pec. Vel. Map	Assumed Selection	H0	dT	${\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})}$	${\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}$	${\widehat{\mathrm{elpd}}}_{\mathrm{loo}}$	p(elpd < Ref.)
		(km s−1 Mpc−1)	(Mpc h−1)	(km s−1)	(km s−1)	log-like
None	Distance limited	${70.8}_{-4.4}^{+4.4}$	${28}_{-2}^{+2}$	Fixed to 264	Fixed to 266	−6.78	0.966
None	Distance limited	${71.6}_{-4.6}^{+4.5}$	${29}_{-2}^{+2}$	${163}_{-35}^{+44}$	${296}_{-69}^{+90}$	−3.73	Ref.
LN21	Distance limited	${71.6}_{-2.8}^{+2.8}$	${29.7}_{-1.4}^{+1.5}$	Fixed to 264	Fixed to 78	−2.46	0.299
LN21	Redshift limited	${75.2}_{-2.5}^{+2.5}$	⋯	Fixed to 265	Fixed to 78	−1.16	0.123
LN21	Distance limited	${71.7}_{-2.0}^{+2.0}$	${30.1}_{-1.2}^{+1.4}$	${136}_{-26}^{+30}$	${73}_{-21}^{+32}$	2.55	0.005
LN21	Redshift limited	${73.2}_{-2.0}^{+2.3}$	⋯	${136}_{-25}^{+29}$	${77}_{-23}^{+35}$	3.10	0.003
C15	Distance limited	${72.7}_{-2.6}^{+2.7}$	${29.8}_{-1.4}^{+1.4}$	Fixed to 264	Fixed to 78	−3.91	0.550
C15	Redshift limited	${76.4}_{-2.4}^{+2.6}$	⋯	Fixed to 265	Fixed to 78	−1.92	0.227
C15	Distance limited	${72.4}_{-1.9}^{+2.0}$	${30.1}_{-1.2}^{+1.3}$	${131}_{-23}^{+27}$	${71}_{-20}^{+30}$	2.11	0.013
C15	Redshift limited	${74.1}_{-2.1}^{+2.4}$	⋯	${130}_{-23}^{+29}$	${79}_{-24}^{+35}$	2.83	0.005
Fiducial result: ${72.9}_{-2.2}^{+2.4}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$

Note. Posterior medians, with upper and lower uncertainties from the 16th and 84th percentiles, of our parameter of interest H₀ and nuisance parameters under permutations of three analysis choices: choice of PV reconstruction, assumed selection of the SH0ES sample, and whether to fix or fit the amplitude of PV covariance. The four bolded values of H₀ are stacked to produce the fiducial result. Final columns show the expected log predictive density (elpd) and associated p-value; large elpd's and small p-values indicate better information on PVs. See Section 4.3 for details.

Download table as:  ASCIITypeset image

Our estimated PV covariance is then

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{\Sigma }}}_{{Pec}}({{\boldsymbol{\chi }}}_{i},{{\boldsymbol{\chi }}}_{j};{\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})},{\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})})={\displaystyle \frac{\partial D}{\partial \tau }}^{2}\\ \ \ \times \,\displaystyle \int {dk}/(2{\pi }^{2})\ {\hat{P}}_{{\rm{\Delta }}}\left(k;{\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})},{\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}\right)\\ \ \ \times \,F\left(k{{\boldsymbol{\chi }}}_{i},k{{\boldsymbol{\chi }}}_{j},{\hat{n}}^{(i)}\cdot {\hat{n}}^{(j)}\right).\end{array}\end{eqnarray} \tag{ 13 }$$

Three Cepheid hosts are identified as members of the Virgo cluster from NED ¹⁰ , but have received different group assignments with distinct redshifts. No other galaxies in our sample appear to be cluster members. To account for the difficulty of the group assignments in the dense environment with the cluster, we add an additional ${\left(550\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}$ of variance to the diagonal elements of the covariance matrix for these three objects, representing the virial noise of the complex (Kashibadze et al. 2020). An alternative approach could assign a single redshift for the cluster center to all three objects, but the width of the cluster is not negligible relative to the distance. Particularly given that identifying a unique estimate for the central redshift of an object with the size and complex structure of Virgo is nontrivial, we conclude our method is more robust.

We calculate the covariance matrix for our sample at fiducial distances using our estimated power spectrum ${\hat{P}}_{{\rm{\Delta }}}(k)$ as well as P_ΛCDM(k) at fixed amplitude. We show the values of the covariance between each pair of objects compared to their physical separation in Figure 2. As can be seen, applying the corrections to remove power at large scales drastically reduces the correlations internal to our sample. As the window function F(u, v, θ) is a complex function of the lines-of-sight of each of the two objects and their angle with respect to to each other, the points fall within an envelope as a function of physical separation. The upper bound is set by a particular geometry of the two objects with respect to the observer with small but nonzero angular separation, while the lower bound will be set by a geometry with the same physical separation but 180° greater angular separation, which will have a similar scale of correlation with opposite sign. The sample of observed Cepheids has few such nearby pairs on opposite sides of the sky at small scales, thus the lower half of the envelope is not seen until larger scales.

Zoom In Zoom Out Reset image size

3.2. Priors and Sample Definition

We use the following priors on the scalar parameters:

• 1.  
  $\mathrm{log}{H}_{0}\sim { \mathcal N }(4,1)$, where H₀ is in km s⁻¹ Mpc⁻¹.
• 2.  
  $\beta \sim { \mathcal N }({\beta }_{\mathrm{Fid}.},{\sigma }_{\beta })$.
• 3.  
  ${V}_{\mathrm{Ext}.}\sim { \mathcal N }({V}_{\mathrm{Fid}.\mathrm{Ext}.},{\sigma }_{V})$.
• 4.  
  $l\sim { \mathcal N }({l}_{\mathrm{Fid}.},{\sigma }_{l})$.
• 5.  
  $b\sim { \mathcal N }({b}_{\mathrm{Fid}.},{\sigma }_{b})$.
• 6.  
  ${\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}\sim \mathrm{InvGamma}(6.52,6.18{\sigma }_{\mathrm{Structure}}^{\mathrm{Fiducial}})$, ¹¹ ${\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}\in (0,5{\sigma }_{\mathrm{Structure}}^{\mathrm{Fiducial}}]$.
• 7.  
  ${\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})}\sim \mathrm{InvGamma}(3.58,3.25{\sigma }_{\mathrm{Unmodelled}}^{\mathrm{Fiducial}})$.

The prior on H₀ was chosen as lognormal because the value is known to be positive, with a width set much larger than the dispersion in modern measurements. The fiducial quantities for the central value and uncertainties on the reconstruction parameters were taken directly from the relevant published work (i.e., C15 or LN21). For the PV amplitude parameters, we chose the inverse gamma distribution as a positive-definite but relatively uninformative distribution. The scale and shape parameters of the distribution were chosen to restrict the model from accessing pathological regions of the posterior at extreme values of these parameters. We solved for the values of these parameters which set the median to the fiducial value and 1% of probability mass at more than 3× the fiducial value for ${\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}$ and 5× the fiducial value for ${\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})}$. ¹² As previously given, the fiducial value of ${\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})}$ is 264 km s⁻¹. The fiducial value of ${\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}$ is dependent on whether we use reconstructions (78 km s⁻¹) or assume the full ΛCDM systematic (264 km s⁻¹). We additionally put an upper bound on ${\sigma }_{\mathrm{Structure}}$ for computational reasons, as in some cases one Markov Chain Monte Carlo (MCMC) chain could become stuck in the pathological region associated with high values of this parameter. In this mode, the PV amplitude parameters become stuck in a local extremum at >10× their fiducial value.

3.2.1. Selection and Distance Priors

For nearby galaxies, the observed redshift as a proxy for distance is significantly biased by volume effects. Since the number density of galaxies per unit distance generally goes as the square of distance, a galaxy observed with a given redshift is likely to have scattered down from the larger volume element with higher true distance than would be naively inferred by the Hubble law. For a PV scatter of 250 km s⁻¹, this bias is approximately given by

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Delta }}}_{z}}{z}\approx \displaystyle \frac{{10}^{-6}}{{z}^{2}},\end{eqnarray} \tag{ 14 }$$

derived from Monte Carlo simulations described below; for our sample, this estimated bias ranges from ∼100% for the lowest redshift objects to ∼1% at the edge of our sample.

In comparison, the Cepheid distances exhibit little statistical noise, resulting in much smaller bias. To quantify and account for the biases of our sample, it is necessary to understand the selection of our sample of Cepheid host galaxies. Specifically, it is necessary to check whether the fact a galaxy was included in our sample provides additional information about the galaxy's distance, independent from the measured values of the redshift and Cepheid distance modulus.

The selection goals of the SH0ES program were ultimately to target every star-forming galaxy at suitable inclination that has hosted a supernova of cosmological quality reported before maximum light and observed digitally in the local volume at d ⪅ 40 Mpc or z ⪅ 0.011. However there is a key ambiguity, in that the indicator used to evaluate this last criterion could be either the redshift of the galaxy or the SN Ia inferred distance. The SH0ES program was conducted over ∼15 years and without any pre-allocation of observing time, so HST proposals were submitted on a cycle-to-cycle basis to target small sets of one to eight hosts of often recently discovered SNe Ia or to take advantage of an extended distance range from a new HST instrument. Near the end of the sample collection, the SH0ES Team targeted the completion of the sample to the above redshift limit, therefore it is difficult for the authors to resolve precisely which limit, distance, or redshift was the most relevant. SH0ES Cepheid program proposals and target lists mix discussion of the redshift or distance limit. Ultimately we conclude that we are unable to determine a priori the selection unambiguously, so we model both and account for this ambiguity as a systematic uncertainty.

To better understand the impact of selection, we ran simple Monte Carlo simulations, in which we sampled galaxies from a three-dimensional volume, using the density values found along SH0ES lines of sight in the Carrick et al. (2015) reconstruction to weight the distances. We simulated redshifts with a 250 km s⁻¹ scatter, and supernova distance moduli with 0.15 mag of scatter. We then compared the observed distribution of redshifts in the SH0ES sample to the simulated redshifts observed when using a strict cut on either redshift or supernova distance in Figure 3 just above the most distant point in the SH0ES sample. We show the predicted bias as a function of simulated CMB redshift from these Monte Carlo simulations in Figure 4 for both scenarios. The bias is very similar for both at z < 0.008 but diverges at higher redshifts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the redshift distribution of the SH0ES sample can be seen to slightly disfavor selection by a redshift cut (a Kolmogorov–Smirnov test gives a p-value = 0.11), neither scenario is inconsistent with the data. As our data cannot rule out either scenario, we produce results assuming two limiting scenarios:

• 1.  
  Assume redshift-limited selection: we will assume that no external distance information of any kind was used in selection, besides the galaxy redshift. This assumption is more conservative as regards final uncertainties, as this reduces the information in our sample.
• 2.  
  Assume distance-limited selection: we assume a strict cut on SN Ia distance modulus was incorporated into the selection. As this assumption entails that the most distant redshifts in our sample are biased high compared to their true distances, this assumption produces a lower value of the Hubble constant. This assumption is more conservative as regards the ongoing Hubble tension; alternative selection models consistent with the data would reduce the dependence of the selection on SN Ia information, reducing the impact.

We account for other aspects of selection for SH0ES in addition to the distance criterion. As SN Ia progenitors are stellar remnants formed in galaxies, our analysis will use the galaxy density distributions found by Carrick et al. (2015) and Lilow & Nusser (2021) to account for the inhomogeneous density of galaxies as a function of distance along each line of sight. Doing so additionally accounts for any bias due to the common practice (until the advent of real-time all-sky surveys) of targeting supernova searches to galaxies selected from catalogs rather than random pointings over the entire sky. ¹³ We do not expect the star-forming galaxies which were then targeted by SH0ES to have galaxy bias parameters significantly different from unity (Swanson et al. 2008), and thus expect the distribution of these galaxies to trace those of the reconstructions.

Lastly, we have checked whether discovery bias could affect our results, i.e., whether that an SN Ia was reported before maximum light provides additional information about the distance of the host galaxy. This would be violated if the magnitude limit for discovery is fainter than the dimmest supernova in our sample. For each supernova we checked who was determined as discoverer by the IAUC. ¹⁴ For each person or collaboration who discovered an SN Ia in our sample, we determined the median magnitude of the supernovae they found at discovery. We deduce that 32 out of 35 SNe were found by discoverers with a median discovery magnitude higher than the peak magnitude of the dimmest SN Ia in our sample. We expect that discovery bias will not significantly affect our measurement, although it may more strongly affect selection at higher distances. If discovery bias were to affect the sample, the effect on the measurement of H₀ would resemble the scenario we call "distance-limited."

Given these considerations, the prior for the latent distances is proportional to the number density of galaxies per unit distance:

$$\begin{eqnarray}&&P(\chi )\propto {\chi }^{2}\times (1+{\delta }_{\mathrm{Rec}.}(\chi \times {\hat{n}}_{i})),\end{eqnarray} \tag{ 15 }$$

using the reconstructed density field δ_Rec.( r ). Tests on simulated data found that the use of this prior tended to bias the recovered distances and H₀ high. We thus modified this prior accordingly to ensure an unbiased measurement. More detail about our simulations is given in Appendix C. The modified prior is

$$\begin{eqnarray}&&f(\chi )={\chi }^{2}\times (1+{\delta }_{\mathrm{Rec}.}(\chi \times {\hat{n}}_{i})),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&P({\boldsymbol{\chi }})=\max {\left({\boldsymbol{\chi }}\right)}^{\alpha }\prod _{i}^{n}\displaystyle \frac{f({{\boldsymbol{\chi }}}_{i})}{{\int }_{0}^{\max ({\boldsymbol{\chi }})}f(x){dx}}\,,\,\,\,{{\boldsymbol{\chi }}}_{i}\in (0,80],\end{eqnarray} \tag{ 17 }$$

where the value α was set to 60 based on the simulation testing. We use the definition of the distance prior given in Equation (17) in our analyses that "assume redshift-limited selection."

If the SH0ES sample was selected using the information from the SNe Ia, then Equation (17) should be modified to include the constraint that the observed SN Ia distance modulus was less than some critical value. As one of the purposes of this work is to test the distance ladder against extreme contaminations in the SN Ia signal, we do not wish to incorporate these values directly into the analysis. We instead adopt a simple model of SN magnitudes as normally distributed with dispersion 0.15 mag and fit a cut-off magnitude/distance as a free parameter. The distance prior becomes

$$\begin{eqnarray}\begin{array}{rcl}f(\chi ,{d}_{T}) & = & {\chi }^{2}\times (1+{\delta }_{\mathrm{Rec}.}(\chi \times {\hat{n}}_{i}))\\ & & \times {\rm{\Phi }}((\mu (\chi )-\mu ({d}_{T}))/.15),\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&P({{\boldsymbol{\chi }}}_{i}| \mathrm{log}{d}_{T})=\displaystyle \frac{f({{\boldsymbol{\chi }}}_{i},{d}_{T})}{{\int }_{0}^{80}f(x,{d}_{T}){dx}}\,,\,\,\,{{\boldsymbol{\chi }}}_{i}\in (0,80\,\mathrm{Mpc}\,{h}^{-1}],\end{eqnarray} \tag{ 19 }$$

where Φ(x) is the cumulative distribution function of the standard normal distribution, and d_T is the assumed distance cutoff in units of Mpc h⁻¹ for inclusion in the SH0ES sample. The upper bound on distances is done purely for performance reasons, and 80 Mpc h⁻¹ is well above the distances of any object in our sample.

In this scenario we also include a prior on the hyperparameter $\mathrm{log}{d}_{T}$ (where d_T is in units of Mpc h⁻¹):

$$\begin{eqnarray}&&\mathrm{log}{d}_{T}\sim { \mathcal N }(\mathrm{log}30,0.4).\end{eqnarray} \tag{ 20 }$$

A lognormal prior was chosen because this parameter is positive-definite, with the width of the prior sufficiently large to allow the data to dominate and sufficiently small to allow efficient convergence of the MCMC. As we are unable to distinguish between these two models of selection a priori, we present results using both priors. In Section 5 we discuss further which model is most appropriate for the data.

3.3. Model Summary

4.1. Computation

We use the code Stan (Carpenter et al. 2017) to sample from the posterior distribution of the parameters as defined by the prior and likelihoods given above. Accounting for the full dependence of the velocity covariance on the latent distances within Stan's Hamiltonian MCMC framework is computationally intractable, requiring evaluation by chain rule of the derivatives of the likelihood with respect to ∼1200 matrix elements for ∼35 parameters. As the effect is small, we use the MCMC to sample from the posterior while keeping the velocity covariance fixed, evaluated at fiducial values for the distances, then check this approximation using importance sampling to reweight each sample appropriately. Given the posterior samples of the parameters θ⁽ⁱ⁾, the relative weight of the ith sample is given by the ratio of the likelihoods of the observed redshifts with the covariance varying/fixed:

$$\begin{eqnarray}&&{w}^{(i)}=\displaystyle \frac{p({{\boldsymbol{z}}}_{\mathrm{CMB}}| {\theta }^{(i)},{\boldsymbol{\Sigma }}({{\boldsymbol{\chi }}}^{(i)}))}{p({{\boldsymbol{z}}}_{\mathrm{CMB}}| {\theta }^{(i)},{{\boldsymbol{\Sigma }}}_{\mathrm{Fid}})}.\end{eqnarray} \tag{ 21 }$$

The overall effect is small, typically ∼0.1σ in H₀.

4.2. Effects of Each Component of Modeling

4.3. Validation

To validate our approach, we wish to ensure that the model is capturing the features found in the real data, as well as compare the use of the C15 and LN21 reconstructions to a more conservative analysis which does not make use of PV information from sources external to SH0ES. To do so we use leave-one-out ("loo") cross-validation based on Vehtari et al. (2017) and Bürkner et al. (2021) to quantify the performance of different models. By refitting the model while withholding one data point from the model, we determine which model most accurately predicts the held-out data. Our metric is the expected log predictive density:

$$\begin{eqnarray}&&{\widehat{\mathrm{elpd}}}_{\mathrm{loo}}^{{ \mathcal M }}=\sum _{i=1}^{N}\mathrm{log}p({{\boldsymbol{\mu }}}_{i}| {{\boldsymbol{\mu }}}_{-i},{ \mathcal M }),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&p({{\boldsymbol{\mu }}}_{i}| {{\boldsymbol{\mu }}}_{-i},{ \mathcal M })=\int p({{\boldsymbol{\mu }}}_{i}| {{\boldsymbol{\mu }}}_{-i},\theta ,{ \mathcal M })p(\theta | {{\boldsymbol{\mu }}}_{-i},{ \mathcal M })d\theta ,\end{eqnarray} \tag{ 23 }$$

where μ _−i is the vector of the distance moduli with the ith entry removed, and ${ \mathcal M }$ denotes the choice of PV reconstructions and other modeling choices. Higher ${\widehat{\mathrm{elpd}}}_{\mathrm{loo}}$ indicates that the model does a better job at predicting the left-out data. As a log-likelihood-based metric, elpd can be interpreted similarly to similar metrics such as Bayesian evidence. Absolute values of this quantity are not informative in themselves, but the differences between models allow discrimination between them. We compare the elpd to a reference, which we choose as the model which uses no SH0ES-external information, assumes distance-limited selection, and fits the amplitude of PV covariance. As the estimated elpd is a summation of N independent terms, we estimate the significance of these differences by the Bayesian bootstrap to evaluate a p-value, which we write as

$$\begin{eqnarray}&&p(\mathrm{elpd}\,\lt \,\mathrm{Ref}.)=p({\mathrm{elpd}}_{\mathrm{True}}^{{ \mathcal M }}\,\lt \,{\mathrm{elpd}}_{\mathrm{True}}^{\mathrm{Ref}.}),\end{eqnarray} \tag{ 24 }$$

which gives the probability that the reference model is more informative than the model under consideration given the data. Small p-values indicate that the reconstructions used are accurate estimates of the underlying PVs.

We can also examine the residuals of the leave-one-out tests by directly comparing the observed value of the distance moduli μ _Obs. to the predicted values when that data point is left out of the model fitting. Since the model is given the correlation structure we expect to be present (through the given covariance matrices and modeled nuisance parameters), if there were any correlations in the corrected residuals, it would show that our model is inadequate to the data. Given that the scatter in the predictions is dominated by Gaussian uncertainties, we expect that the residuals should be approximately normal in distribution. We show a comparison of the residuals relative to the predictions in Figure 5, and a comparison of predicted redshifts to observed distance moduli in Figure 6. We test the corrected residuals for autocorrelation or a linear trend, and find no significant evidence of either at >2σ. The positive tail in the uncorrected residuals seen in the upper panel of Figure 5 at low distances is a distinctive artifact of the volumetric or Malmquist bias due to scatter in redshifts. The same artifact at these redshifts can be seen in the SN Ia Hubble residuals of the Pantheon+ analysis (Brout et al. 2022). These data were not used in the Pantheon+ cosmology fit, which cut the sample at z = 0.01, nor in the SH0ES analysis R21, which cut the sample at z = 0.023. The bias is not expected to affect either result, as it scales as z⁻² and is negligible in both samples.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We first measure the Hubble constant without the use of PV reconstructions. This gives us a measurement independent of PV reconstruction procedures, and allows us to quantify the improvement from the use of this data. We use the full ΛCDM power spectrum as described in Equation (11), set v_pred( χ _i ) and δ( χ _i ) to 0 in Equations (4) and (19), and remove the reconstruction parameters V_Ext., l, b, and β from the model. Without the use of other PV information and assuming distance-limited selection, we measure H₀ to be ${70.8}_{-4.4}^{+4.4}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, a precision of ∼6.3%. Fitting the PV amplitudes leads to a smaller estimate of the small-scale or nonlinear power, yielding ${H}_{0}={71.6}_{-4.6}^{+4.5}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$. This choice is significantly favored by the data. We take the latter model as the reference for a two-rung measurement without the use of PV information. We show the values for H₀ and the other parameters in Table 2. Ignoring the galaxy catalog information in the case of an assumed redshift cut produces pathological behavior in the posterior. This failure mode appeared in tests on simulated data, showing similar properties to the behavior referenced at the end of Section 3.2, wherein the PV amplitude parameters could not be effectively constrained by the data. We thus do not produce results for this scenario.

Next we fit the data using each of the two reconstructions, two choices of prior, and the choice to either fit the amplitude of PV covariances or fix them to fiducial values. We report the appropriate parameter constraints in Table 2. Including information from the reconstructions, the fractional error in H₀ decreases from ∼6% to ∼3%, with the assumption of SN Ia selection decreasing uncertainties further. We show the corner plot of the constraints on the scalar parameters using the C15 reconstruction with fitted PV covariance in Figure 7 under both selection scenarios.

Zoom In Zoom Out Reset image size

As expected, use of either reconstruction is significantly favored by the elpd metrics. Neither reconstruction is significantly favored compared to the other. The choice to fit the PV amplitudes is also strongly favored by the data. The data slightly favors the redshift-limited scenario for selection, but the difference is not significant (p > 0.2 for both LN21 and C15 when PV amplitudes are fitted). We conclude that we are incapable of differentiating the models of selection from the data alone. As the redshift-limited scenario implies a more restricted experimental process, we consider both equally likely given the data.

Inclusion of the reconstructions shifts the central value by ∼1.5 km s⁻¹ Mpc⁻¹, much less than the ∼3 km s⁻¹ Mpc⁻¹ reduction in the uncertainty, and thus within expectation. We note that the two reconstructions produce very similar results regardless of the other analysis choices made, consistent with our error budget. The different approaches to selection shift results by more, increasing the value of H₀ by ∼2–5 km s⁻¹ Mpc⁻¹, depending on whether or not the PV covariance is also fit. No result in Table 2 is more than 1.6σ discrepant with any other value.

The fitted values of the amplitude of PV covariance are consistent with our expectation. As we have used group-corrected redshifts, we expect that much of the diagonal scatter has been removed, lowering the value of ${\sigma }_{\mathrm{Unmodelled}}^{(\mathrm{PV})};$ our fitted values are consistent with the estimated 150 km s⁻¹ for galaxy halos found by comparison to N-body simulations in Carrick et al. (2015), and assumed in Betoule et al. (2014), Boruah et al. (2020), and others. Additionally, this amplitude is roughly consistent with the estimates of the uncertainty shown in Figures 4 and 9 of LN21. However, Brout et al. (2022) finds a PV scatter of 240 km s⁻¹ in the Pantheon+ supernova sample. This number was found in a mixed sample of galaxies with and without group assignments, at larger distances where group assignments are more uncertain. Additionally, PV corrections for the baseline Pantheon+ analysis were not integrated over the line of sight, which tends to increase statistical dispersion in distances (Peterson et al. 2021). The values of ${\sigma }_{\mathrm{Structure}}^{(\mathrm{PV})}$ are in all cases consistent with their fiducial values, showing that our estimates of the magnitude of errors in the reconstructions are reasonable and that the scale of observed local structures is consistent with that expected from ΛCDM.

As the elpd strongly favors the use of PV reconstructions and strongly favors fitting of the PV amplitudes, we use these four models to estimate the final result including the systematic uncertainty. We assume that these four models are each, in the posterior, equally likely to be representative of the true data-generating process. ¹⁵ Under this assumption, and given that the data are common between all models, the posterior distribution of H₀ is found by stacking (adding together) the individual distributions from these four models with equal weight. The resulting median value for H₀ is ${H}_{0}={72.9}_{-2.2}^{+2.4}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ (with 16% and 84% percentile uncertainties). We take this as our fiducial result. We evaluate a p-value for the difference between our value and the Planck measurement of H₀ = 67.4 ± 0.5 km s⁻¹ Mpc⁻¹ (Planck Collaboration et al. 2018). Comparing the MCMC chains, we find p = 0.007, corresponding to a significance of 2.4σ. We compare our result to R21 and Planck Collaboration et al. (2018) in Figure 8.

Zoom In Zoom Out Reset image size

5.1. Comparison to SNe Ia Hubble Diagram

As an additional cross-check, we compare the distances we measure from our models to the SN Ia Hubble diagram. We use the standardized magnitudes of 73 light curves of cosmological SNe from the Pantheon+ sample (Scolnic et al. 2021; Brout et al. 2022) with common host galaxies to our Cepheid observations. To evaluate a goodness-of-fit metric, we compute a χ² by comparing predicted values of μ from the model to SN magnitudes, fitting the absolute magnitude of the SNe Ia, and including the covariance of the predicted μ's and the covariance of the supernovae provided by Brout et al. (2022). χ² values do not vary significantly from model to model, between 49.4 and 51.0, with 72 degrees of freedom. That the values do not change significantly is not surprising, as the Cepheid distances have the least statistical error of any indicator used in this analysis, and the predicted distances are largely driven by the Cepheid measurement. We show a comparison of these two rungs of the distance ladder in Figure 9. Our results are consistent with the results of R21.

Zoom In Zoom Out Reset image size

We find that the Hubble constant can be constrained at extremely low redshifts by taking advantage of measurements of the local density and velocity fields from external sources. Our code and data will be made public as of the acceptance of this paper. Values of the Hubble constant measured when incorporating this information, depending on further modeling assumptions, range from 71.7 km s⁻¹ Mpc⁻¹ to 76.4 km s⁻¹ Mpc⁻¹. Our data supports a credible interval of ${H}_{0}={72.9}_{-2.2}^{+2.4}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, at 2.4σ tension with the measurement of Planck. As a result, we conclude the cause of the Hubble tension is unlikely to be an unknown systematic of SNe Ia.

Our analysis is sensitive to the modeling of the Cepheid variable stars, including reddening, crowding, and metallicity corrections. However any systematic errors or uncertainties in this stage of the analysis will be mitigated by our inclusion of a covariance matrix derived from the analysis variants of R21.

We find that the most significant uncertainty of our analysis is the selection function of the SH0ES sample. Seeking a further reduction of the Hubble constant would require informative distance priors with increased probability mass at low distances. However, as our data show good consistency with the observed distributions of distances, increasing the probability of low distances will require reductions in the information of the prior. Decreasing the sum of the PV amplitude parameters σ_PV reduces the estimated volumetric bias and thus will tend to reduce the value of H₀; however, doing so below 150 km s⁻¹ seems unmotivated by the data or theoretical estimates of the PV noise. We conclude that alternative selection models beyond those discussed are unlikely to significantly reduce our measured value of H₀.

As the R21 analysis is based on comparison of indicators with much lower relative scatter than the redshifts we use here, the analysis is not likely to be strongly affected by selection effects. Distance selection and biases in the Hubble sample is accounted for using the BBC methodology (Kessler & Scolnic 2017). However, any analysis which makes use of high-scatter indicators, including nearby redshifts as a measure of the Hubble flow, must carefully evaluate the selection of their sample. This has been done in measurements of H₀ from "bright sirens" in Abbott et al. (2017) as well as "dark sirens" (Soares-Santos et al. 2019). Boruah et al. (2020) found a nontrivial effect from the use of a volumetric prior in a measurement of H₀ from water mega-masers using data from Reid et al. (2019), and we suggest future works more closely examine the selection of mega-maser, SBF, and Tully–Fisher distance measurements.

We chose here to estimate the sytematic uncertainties in cosmic density/PV reconstructions from the differences between the published data products of Carrick et al. (2015) and Lilow & Nusser (2021). Other methods of estimating these uncertainties may include the use of a fully Bayesian reconstruction procedure such as the BORG method (Jasche & Lavaux 2019) or direct computation of covariances from constrained realizations produced using the CORAS code (Lilow & Nusser 2021). While off-diagonal terms are not provided in either C15 or LN21, both provided an estimate of the diagonal variance. The PV amplitude parameters we derive show good consistency with those shown in Figures 4 and 9 of LN21 for our region of interest and Appendix A of Carrick et al. (2015).

Cosmological uncertainties are unlikely to affect the result. We assume Gaussian behavior in our calculation of the PV uncertainties, while the cosmological density field is known to show nonlinear and non-Gaussian behavior, particularly at scales k > 0.1 Mpc h⁻¹. However, the velocity field is comparatively insensitive to these effects, as the velocity field is a factor of 1/k² less sensitive to small scales as compared to the density. We also neglect any velocity bias, which is absent on large scales in the absence of nongravitational forces as we deal with in this work (see, e.g., Desjacques et al. 2018). We thus expect our result to be insensitive to the effects quantified by higher-order statistics. Our error estimates also assume ΛCDM physics, but are only weakly dependent on alternative physics through changes to the power spectrum at small scales, given that the power spectrum at large scales has been constrained by numerous independent experiments. Our estimates of the uncertainties in PV reconstructions are optimistic, insofar as they assume that the reconstructed velocity fields are no more inaccurate than the differences between the two. At scales larger than the volume of the Cepheid sample, we are unable to effectively constrain these quantities ourselves. At scales smaller than this, we have mitigated these uncertainties by introducing amplitude parameters that are marginalized over in the final result. This data and others (Said et al. 2020; Peterson et al. 2021) strongly support the use of these reconstructions as broadly applicable models of the nearby universe.

This work provides a measurement of H₀ independent of any supernovae sample. At the expense of increased uncertainties, dropping supernova measurements from the analysis has three main advantages. First, this work provides a satisfactory systematic error check of the SH0ES analysis: our result is fully compatible with the findings of R21, potentially ruling out any sizable systematic error contribution from the inclusion of the supernovae. Further, our measurement of H₀ is uncorrelated with the absolute magnitude of higher-redshift supernovae, which circumvents any potential issues when using this value as a prior in analyses which additionally use SNe Ia as a constraint on the evolution of the distance-redshift relation. Lastly, as the sample's maximum redshift is z ∼ 0.01, our result is robust to any late-time effects such as transitions in the dark content of the universe unless they have occurred within the last ∼100 Myr.

Measurements of the Hubble constant have been conducted at low redshift using a variety of distance indicators. The treatment of PVs has in some cases been a matter for discussion in the literature (Howlett & Davis 2020; Boruah et al. 2020; Mukherjee et al. 2021). Our methodology shares several elements with Boruah et al. (2020), with improvements to the treatment of off-diagonal velocity uncertainties, and a more flexible model framework parametrized to reduce cosmology dependence. We anticipate the use of similar methodologies to that presented here in future Cepheid-based measurements of the Hubble constant from the SH0ES team, and suggest these elements should be incorporated into other analyses based on low-redshift distance indicators.

We are thankful to Tamara Davis and Anthony Carr for their contributions to the quality of the Pantheon+ redshifts, and comments on this paper.

Support for this work was provided by NASA through the NASA Hubble Fellowship grant No. HF2-51462.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This research was supported by NASA/HST grant No. 16269 from the Space Telescope Science Institute. D.S. is supported by DOE grant No. DE-SC0010007 and the David and Lucile Packard Foundation. This project/publication was made possible through the support of grant #62314 from the John Templeton Foundation. J.L.B. is supported by the Allan C. and Dorothy H. Davis Fellowship.

Software: AstroPy (Astropy Collaboration et al. 2013, 2018), corner.py (Foreman-Mackey et al. 2021), Matplotlib (Hunter 2007), nbodykit (Hand et al. 2018), NumPy (Harris et al. 2020), pystan (Riddell et al. 2021), SciPy (Virtanen et al. 2020), Stan (Carpenter et al. 2017), tqdm (da Costa-Luis et al. 2021).

In order to quantify the systematics associated with methodological choices made internal to the Cepheid analysis, we include additional systematic uncertainties from the analysis variants of R21. The baseline SH0ES analysis quantified these effects by measuring H₀ for each variant, taking the standard deviation of these as a systematic uncertainty in quadrature. This procedure in general will overestimate the uncertainty, as the analysis is not given the opportunity to mitigate the effect on the measurement by more strongly weighting the data less affected. Since the R21 analysis is generally dependent only on changes in the mean of the calibrator Cepheids relative to the anchor Cepheids, this overestimate would not significantly inflate the final error budget. However, our measurement is more sensitive to changes in the relative values of the Cepheid host galaxies due to the correlations among the host galaxy PVs. By incorporating information about the effect of the analysis variants directly into the model, we minimize the effect on our measurement.

For a specific effect (e.g., treatment of outliers in the Cepheid period–luminosity relation) we define a set of analysis variants which represent reasonable alternative treatments of the baseline SH0ES Cepheid sample. For an effect X represented by n analysis variants, where the ith variant produces distance moduli ${{\boldsymbol{\mu }}}_{{\rm{X}}}^{(i)}$, we define a systematic covariance matrix by the average of the outer products of the residuals relative to the baseline treatment:

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\mu }^{({\rm{X}})}=\displaystyle \frac{1}{n}\sum _{i=0}^{n}({{\boldsymbol{\mu }}}_{{\rm{X}}}^{(i)}-{{\boldsymbol{\mu }}}_{\mathrm{Ceph}})\otimes ({{\boldsymbol{\mu }}}_{{\rm{X}}}^{(i)}-{{\boldsymbol{\mu }}}_{\mathrm{Ceph}}).\end{eqnarray} \tag{ A1 }$$

The individual covariance matrices are added together to evaluate the total uncertainty in the distance moduli:

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\mu }={{\boldsymbol{\Sigma }}}_{\mu }^{(\mathrm{SH}0\mathrm{ES})}+\sum _{X}{{\boldsymbol{\Sigma }}}_{\mu }^{({\rm{X}})}.\end{eqnarray} \tag{ A2 }$$

We briefly discuss the analysis variants we use here; for full details of the treatment of the data under each, refer to Section 6 and Table 5 of R21. We budget for the potential systematic impact of four different steps in the Cepheid analysis methodology:

• 1.  
  Cepheid outlier treatment: the selection of our Cepheid samples from optical data is expected to leave some residual contamination by blended variable stars in the sample. Our baseline analysis therefore includes outlier rejection by individually discarding the largest outlier until no Cepheids deviate from the global period–luminosity fit by >3.3σ. We include variants 2–8 as viable alternative treatments, excluding variant 9, which discards outlier rejection entirely.
• 2.  
  Reddening, extinction, and color relations: our baseline analysis uses dereddend Wesenheit magnitudes based on the Fitzpatrick (1999) reddening law with the absolute-to-selective extinction ratio R_V fixed to 3.3. Alternative treatments include allowing the assumed value of R_V as a free parameter, alternative reddening laws, and individually inferred values of R_V for each Cepheid host galaxy using the empirical dust attenuation framework of Hahn et al. (2022). We include variants 16,17, 19, 20, 21, and 22 as viable alternative treatments. Variant 23 excludes reddening standardization entirely and is omitted from our budget.
• 3.  
  Period standardization: the baseline sample includes all Cepheids with observed periods P > 5 days, and uses a linear P-L relation with a single slope. As some studies have found broken period–luminosity relations for Cepheids at P = 10 days we include a variant (number 24) which allows different slopes above and below P = 10 days. The difference in slope is consistent with 0 within 2σ. Other variants exclude high- or low-period Cepheids from the sample (variants 25 and 26). All three variants are included as viable alternative treatments.
• 4.  
  Metallicity: the baseline variant includes a linear standardization of Cepheid magnitudes for the metallicity of the host galaxy. The free parameter γ, the slope of the metallicity–luminosity relation, is significant at ∼4.6σ. The Cepheid metallicities are inferred based on Teimoorinia et al. (2021). We include a single-variant analysis, which uses instead the metallicity scale of Pettini & Pagel (2004), as discussed further in Section 6.6 and Appendix C of R21. We exclude variant 30, which discards metallization standardization entirely.

In Figure 10, we show both the baseline covariance ${{\boldsymbol{\Sigma }}}_{\mu }^{(\mathrm{SH}0\mathrm{ES})}$ and the full covariance including Σ_μ . We show how including each effect discussed above in the systematic covariance matrix affects the mean and variance of the H₀ measurement in Table 3. The systematic effect with the most impact on the final measurement is the behavior of the Cepheid period–luminosity relation.

Zoom In Zoom Out Reset image size

Table 3. Mean and Standard Deviation of H₀ as Measured with the Baseline Cepheid Covariance ${{\boldsymbol{\Sigma }}}_{\mu }^{(\mathrm{SH}0\mathrm{ES})}$, and with each Component of the Systematic Uncertainty Added Individually

Systematics	Δ 〈 H0 〉 from Baseline	σ H0	Contribution to Error Budget
${{\boldsymbol{\Sigma }}}_{\mu }^{(\mathrm{SH}0\mathrm{ES})}$ only	0.4	2.4	Ref.
Outlier treatment	0.4	2.4	0.2
Reddening and color	0.3	2.4	0.0
Period–luminosity relation	0.8	2.4	0.3
Metallicity scale	0.4	2.4	0.1
Σμ (All Sys. Inc.)	0	2.5	0.7

Note. Error budget contribution defined as $\sqrt{{\sigma }_{\mathrm{Effect}}^{2}-{\sigma }_{\mathrm{Ref}.}^{2}}$. All results used the model variation assuming redshift-limited selection, fixed the PV amplitude, and used the C15 reconstruction, to ensure the assumed errors in the data are kept constant apart from the changes in the Cepheid covariance.

Download table as:  ASCIITypeset image

R21 presented several other analysis variants. Many of these are not treated here as they impact only the SN Ia sample rather than the Cepheid distances. We also do not make use of analysis variants which ignore large amounts of data (the optical Wesenheit variants, which wholly discards all NIR data) or ignore components of the standardization (such as removing metallicity dependence from the Cepheid luminosity relations). As we are interested in making the most accurate measurement possible given all available data, these analysis variants do not represent reasonable treatments. In general our final results will be more sensitive than the R21 result to the exclusion of large amounts of data or the removal of components of the standardization.

We first evaluate the power spectrum of the C15 reconstruction. As the C15 reconstruction has masked some parts of their volume due to a lack of data, we use a 212 Mpc h⁻¹ cube centered on the Milky Way. This volume avoids any masking, and encompasses the whole of our sample. We compare with the predicted ΛCDM power in Figure 11. As can be seen, the reconstruction's power spectrum is well approximated by ΛCDM smoothed with a Gaussian kernel at a scale of 5 Mpc h⁻¹. The dispersion of the line-of-sight PVs within the reconstruction is smaller than the ${\left(380\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}$ expected from ΛCDM: $\langle {v}_{i}^{2}\rangle \sim {\left(270\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}$ evaluated using the corresponding power spectrum. It is clear that the reconstruction cannot correct for the half of the dispersion that has not been modeled. More formally, our quantities of interest are the differences between the estimated and true density/velocity fields, as represented by the density residuals:

$$\begin{eqnarray}&&{\delta }_{{\rm{\Delta }}}({\boldsymbol{r}})={\delta }_{{\rm{C}}15}({\boldsymbol{r}})-{\delta }_{\mathrm{True}}({\boldsymbol{r}}).\end{eqnarray} \tag{ B1 }$$

Our goal is to determine the power spectrum of this quantity, which will enable us to budget for unmodeled PV noise and quantify the improvement in our result from the use of the reconstructions. The power spectrum of the residuals is

$$\begin{eqnarray}&&{P}_{{\rm{\Delta }}}(k)={P}_{\mathrm{True}}(k)+{P}_{\ {\rm{C}}15}(k)-2{P}_{\mathrm{True},{\rm{C}}15}(k),\end{eqnarray} \tag{ B2 }$$

where P_True,C15 is the cross-spectrum term

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{True},{\rm{C}}15}(k) & = & \displaystyle \int {d}^{3}{k}^{* }/{\left(2\pi \right)}^{3}\lt {\tilde{\delta }}_{{\rm{C}}15}^{* }({{\boldsymbol{k}}}^{* }){\tilde{\delta }}_{\mathrm{True}}({\boldsymbol{k}})\\ & & +{\tilde{\delta }}_{{\rm{C}}15}({\boldsymbol{k}}){\tilde{\delta }}_{\mathrm{True}}^{* }({{\boldsymbol{k}}}^{* })\gt /2.\end{array}\end{eqnarray} \tag{ B3 }$$

This term is essentially a measure of the alignment in phase between the true and reconstructed density fields, while the power spectrum is phase independent. In the case where the reconstruction is perfectly correct, this cross-spectrum term is equal to the power spectrum, and, as we expect, the power in the residuals drops to zero. While we can assume an underlying ΛCDM power spectrum for the true density fields and the power spectrum of the reconstruction can be measured, evaluating the cross-spectral term requires us to determine the "correctness" of the reconstructions, which cannot be done from first principles.

Zoom In Zoom Out Reset image size

As an estimate of this term from the data we have available, we have chosen to base our error budget on the differences between the results of C15 and LN21. This is likely to give an underestimate of the true uncertainties, as both are based on the 2MRS redshifts published in Huchra et al. (2012). However, this provides an estimate of the methodological uncertainties that includes off-diagonal terms of the covariance, not otherwise available to us. We test this method on our data in Section 5, as well as parameterizing the scale of covariance in the data and marginalize over these quantities. Ultimately the Cepheid data significantly favor the use of reconstructions with our budgeted uncertainties, and we conclude we have found a reasonable approximation of our true uncertainties.

By assuming ${P}_{\mathrm{True},{\rm{C}}15}(k)={P}_{{\rm{C}}15,\mathrm{LN}21}(k)$, and using the relation

$$\begin{eqnarray}&&2{P}_{{\rm{C}}15,\mathrm{LN}21}(k)={P}_{{\rm{C}}15}(k)-{P}_{{\rm{C}}15-\mathrm{LN}21}(k)+{P}_{\mathrm{LN}21}(k),\end{eqnarray} \tag{ B4 }$$

we derive an estimate of the power spectrum which describes the uncertainty in the reconstructed density fields

$$\begin{eqnarray}&&{\hat{P}}_{{\rm{\Delta }}}(k)={P}_{\mathrm{Unmodelled}}(k)+{P}_{\mathrm{Structure}}(k),\end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{Unmodelled}}(k)={P}_{{\rm{\Lambda }}\mathrm{CDM}}(k)-{P}_{\mathrm{Reconstructed}}(k),\end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{Structure}}(k)={P}_{{\rm{C}}15-\mathrm{LN}21}(k),\end{eqnarray} \tag{ B7 }$$

where ${P}_{{\rm{C}}15\ -\ \mathrm{LN}21}(k)$ is the power spectrum of the residuals between the two reconstructions ${\delta }_{{\rm{C}}15}({\boldsymbol{r}})-{\delta }_{\mathrm{LN}21}({\boldsymbol{r}})$. Rather than directly using the power spectrum evaluated from the data cubes, for both C15 and LN21 we use the smoothed ΛCDM power spectrum,

$$\begin{eqnarray}&&{P}_{\mathrm{Reconstructed}}(k)={P}_{{\rm{\Lambda }}\mathrm{CDM}}(k)\exp \left(-{k}^{2}{\sigma }_{\mathrm{Smooth}}^{2}\right),\end{eqnarray} \tag{ B8 }$$

with a smoothing scale σ_Smooth = 5 Mpc h⁻¹. Our estimate of the uncertainties then has a fairly simple and intuitive form: subtract the modeled physical information from the theoretical prediction and add the power describing the differences between our two models. The relative contribution of the two can be seen in Figure 1; the power from the differences adds ${\left(78\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}$ of variance, while the nonlinear term contributes ${\left(264\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}$ at much smaller scales. As we have no strong evidence to favor one reconstruction over the other, our error estimate is the same for both, and we produce results using both.

To test whether our methodology produced unbiased estimates of H₀, we simulated data using both selection scenarios, comparing the simulated "truth" values to those recovered using our model. We use the C15 reconstruction as our baseline for the simulations. The scalar parameters of the model are sampled according to the priors given in Section 3.2. We then simulate latent distances χ and data for both selection scenarios according to the following scheme:

• 1.  
  Distance-limited selection: 35 galaxies are randomly drawn from the C15 data cube, weighted by the density in each voxel. Each voxel is also weighted by a factor Φ((μ(χ) −μ(d_T ))/.15) as a function of the distance from the Milky Way galaxy. Peculiar velocities and distance moduli are then randomly generated as given in Equations (4) and (3) to produce simulated redshift and distance vectors.
• 2.  
  Redshift-limited selection: for this scenario a somewhat more complex iterative method is necessary to reproduce the correlation structure of the observed redshifts. Proposed galaxies are randomly drawn from the C15 data cube, weighted by the density in each voxel. For each proposed galaxy, a PV is drawn given the PVs of the galaxies accepted to that point, with a covariance matrix as specified by Equation (13). The proposed galaxy is trialed against a cut on the observed redshift z < 0.011; if the galaxy is accepted, it is added to the sample. This algorithm is iterated until the sample reaches 35. The simulated distance moduli are then generated as in Equation (3).

In either case, we analyze the simulated data using the model which fits the PV amplitude parameters and uses the C15 reconstruction to measure H₀. We repeat this procedure 100 times to allow us to evaluate any biases in the recovered values. We inspect the recovered latent distance parameters χ and H₀ as compared to the true values with which the simulated data were generated. In Figure 12 we show that the recovered distances match the distribution of distances simulated. In Figure 13 we show the z-scores of the recovered H₀ values compared to the simulated value are consistent with a mean of 0 in both scenarios. We conclude that our method is unbiased to acceptable precision.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size