The Megamaser Cosmology Project. XIII. Combined Hubble Constant Constraints

Ninety years after Hubble's seminal work (Hubble 1929), observational cosmology remains focused on obtaining a precise and accurate value of the Hubble constant, H₀. Today, measurements of the cosmic microwave background (CMB) at high redshift (z ≈ 1100) determine the angular-size distance to the surface of last scattering, and within the context of the standard ΛCDM cosmological model, predict a precise value for H₀ of 67.4 ± 0.5 km s⁻¹ Mpc⁻¹ (Planck Collaboration et al. 2018). Because this "early-universe" prediction is model-dependent, complementary "late-universe" measurements of H₀ provide an important test of the assumed cosmological model (e.g., Hu 2005).

Currently, measurements of H₀ at low redshifts (z ≪ 10) are in statistical tension with the early-universe prediction. For example, distance ladder measurements using Cepheid variables to calibrate the absolute luminosities of Type Ia supernovae find 74.03 ± 1.42 km s⁻¹ Mpc⁻¹ (Riess et al. 2019), and time-delay strong lensing measurements from multiply imaged quasars currently yield ${73.3}_{-1.8}^{+1.7}$ km s⁻¹ Mpc⁻¹ (Wong et al. 2019). Though this discrepancy between the early- and late-universe H₀ measurements is being taken increasingly seriously by the cosmological community (see Verde et al. 2019), the tantalizing prospects that it holds for heralding physics beyond ΛCDM (e.g., Raveri et al. 2017; Poulin et al. 2018) require that we subject it to a correspondingly strict evidence threshold. Independent avenues for constraining H₀ are thus necessary to provide cross-checks against unrecognized systematics in the measurements.

Water megamasers residing in the accretion disks around supermassive black holes (SMBHs) in active galactic nuclei (AGNs) provide a unique way to bypass the distance ladder and make one-step, geometric distance measurements to their host galaxies. The archetypal AGN accretion disk megamaser system is located in the nearby (7.6 Mpc; Humphreys et al. 2013; Reid et al. 2019) Seyfert 2 galaxy NGC 4258 (Claussen et al. 1984; Nakai et al. 1993; Herrnstein et al. 1999). Very long baseline interferometric (VLBI) observations reveal that the masers trace the accretion disk on sub-parsec scales, where the SMBH dominates the gravitational potential. The masing gas parcels act as test particles in this potential and exhibit the ordered, Keplerian motion expected for orbits about a point mass (Greenhill et al. 1995; Miyoshi et al. 1995). By combining the VLBI position and velocity information with centripetal accelerations measured from spectral monitoring observations (e.g., Argon et al. 2007), the typical degeneracy between mass and distance is broken and precise constraints can be placed on both quantities (e.g., Herrnstein et al. 1999; Humphreys et al. 2013; Reid et al. 2019).

The Megamaser Cosmology Project (MCP) is a multi-year campaign to find, monitor, and map AGN accretion disk megamaser systems (Braatz et al. 2007, 2008). The primary goal of the MCP is to constrain H₀ to a precision of several percent by making geometric distance measurements to megamaser galaxies in the Hubble flow (Kuo et al. 2013, 2015; Reid et al. 2013; Gao et al. 2017). Distance measurements made using the megamaser technique do not rely on distance ladders¹⁴ or the CMB, and they have different systematics than lensing-based techniques. The megamaser measurements thus provide an independent handle on H₀.

In this paper we present a revised analysis that improves the distance measurements for several megamaser systems that have been previously published by the MCP. We then combine all distance measurements made using the revised analysis into a single megamaser-based constraint on H₀. The paper is organized as follows. In Section 2 we discuss the new disk modeling and compare our revised distance estimates with the previously published results. In Section 3 we describe how we combine the distance measurements into one H₀ constraint, detailing several different peculiar velocity treatments, and we present our resulting H₀ measurement. We summarize and conclude in Section 4.

In this section we present updated distance measurements for the megamaser-hosting galaxies UGC 3789, NGC 6264, NGC 6323, and NGC 5765b. Each of these galaxies has had a maser-derived distance measurement published previously by the MCP (see Kuo et al. 2013, 2015; Reid et al. 2013; Gao et al. 2017).

The improvements we present here stem primarily from an update to the fitting procedure that incorporates the "error floor" systematic uncertainties as model parameters, thereby enabling marginalization over a previous source of systematic uncertainty. These error floors get added in quadrature with the data errors, such that an error floor value of zero indicates that the data errors already capture the true measurement uncertainties well. This updated model has already been applied to the galaxies NGC 4258 (Reid et al. 2019) and CGCG 074-064 (Pesce et al. 2020), so for these galaxies we simply use the corresponding published results.

We use a Hamiltonian Monte Carlo (HMC) sampler as implemented within the PyMC3 package (Salvatier et al. 2016) to perform all of the disk fitting described in this section; for a comprehensive overview of the disk model and fitting procedure, see Pesce et al. (2020). The HMC code yields results that match well with those from the Metropolis–Hastings disk-fitting code employed in many previous MCP publications (see, e.g., Reid et al. 2013), but has an improved convergence efficiency. We use broad truncated Gaussian priors for all error floor parameters, with a mean of 10 μas and a standard deviation of 5 μas for both the x- and y-position error floor priors, a mean of 2 km s⁻¹ and a standard deviation of 1 km s⁻¹ for both the systemic and high-velocity¹⁵ error floor priors, and a mean of 0.3 km s⁻¹ yr⁻¹ and a standard deviation of 0.15 km s⁻¹ yr⁻¹ for the acceleration error floor priors. In all cases, the prior distribution is truncated at zero so that the error floors remain strictly positive.

2.1. UGC 3789

The disk megamaser in the galaxy UGC 3789 was discovered by Braatz & Gugliucci (2008), and the first VLBI map of the system was presented in Reid et al. (2009). This map was combined with acceleration measurements by Braatz et al. (2010) to produce an angular-size distance estimate of 49.9 ± 7.0 Mpc to the system. Reid et al. (2013) obtained additional data and improved the distance measurement to 49.6 ± 5.1 Mpc, and we use the data from that latest paper to produce the updated measurements presented in Table 1 and the Appendix.

Table 1.  Results from Maser Disk Modeling

Galaxy	Distance (Mpc)	Velocity (km s−1)	Reference
UGC 3789	${51.5}_{-4.0}^{+4.5}$	3319.9 ± 0.8	this work
NGC 6264	${132.1}_{-17}^{+21}$	10192.6 ± 0.8	this work
NGC 6323	${109.4}_{-23}^{+34}$	7801.5 ±  1.5	this work
NGC 5765b	${112.2}_{-5.1}^{+5.4}$	8525.7 ±  0.7	this work
CGCG 074-064	${87.6}_{-7.2}^{+7.9}$	7172.2 ±  1.9	Pesce et al. (2020)
NGC 4258	7.58 ± 0.11	679.3 ± 0.4	Reid et al. (2019)

Note. Maser galaxy distances and velocities as measured from modeling the maser disks; for each value we quote the posterior median and 1σ confidence interval (i.e., 16th to 84th percentile). For NGC 4258, where the systematic uncertainty in the distance measurement is comparable to its statistical uncertainty, we have added the two in quadrature. All velocities are quoted in the CMB reference frame using the optical convention. The values for the other parameters measured from the disk model are given in the Appendix.

Download table as:  ASCIITypeset image

Our updated fit indicates that the error floors used by Reid et al. (2013) for the position measurements were too conservative; we find a best-fit x-position error floor of 5 ± 1 μas and a best-fit y-position error floor of 6 ± 1 μas, compared to the previous (fixed) values of 10 μas for both x- and y-positions. We find that the data do not provide strong constraints on the velocity measurement error floors for either the systemic or the high-velocity masers, and that the posteriors for both largely follow the Gaussian-distributed prior range of 2.0 ± 1.0 km s⁻¹. Our best-fit acceleration error floor is 0.34 ± 0.06 km s⁻¹ yr⁻¹, again indicating that the Reid et al. (2013) value of 0.57 km s⁻¹ yr⁻¹ was too conservative.

The net result of the updated modeling is a distance measurement of $D={51.5}_{-4.0}^{+4.5}$ Mpc, which represents a 3.8% ≈ 0.4σ increase over the previous value and an improvement in the measurement precision from ±10% to $(+8.7 \% ,-7.8 \% )$.

2.2. NGC 6264

2.3. NGC 6323

2.4. NGC 5765b

Our disk model returns an angular-size distance measurement ${\hat{D}}_{i}$ and a redshift measurement ${\hat{z}}_{i}$ for the SMBH in each megamaser-hosting galaxy. From considerations of kinetic energy equipartition between the SMBH and surrounding stars, we expect the relative velocity of a ∼10⁷ M_⊙ black hole to be ≪1 km s⁻¹ with respect to the system barycenter (Merritt et al. 2007). Measured upper limits on the magnitude of this relative motion are on the order of several km s⁻¹ for the SMBH in the center of the Milky Way (Reid & Brunthaler 2004, 2020) and no more than a few tens of km s⁻¹ for the sources considered in this work (Pesce et al. 2018). We thus proceed under the assumption that each galaxy effectively shares a distance and redshift with its SMBH, and we seek to determine what values of H₀ are compatible with these measurements.

For each galaxy, its expected angular-size distance D_i is related to its expected cosmological recession redshift z_i and H₀ by

$$\begin{eqnarray}\begin{array}{rcl}{D}_{i} & = & \displaystyle \frac{c}{{H}_{0}\left(1+{z}_{i}\right)}{\displaystyle \int }_{0}^{{z}_{i}}\displaystyle \frac{{dz}}{\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+\left(1-{{\rm{\Omega }}}_{m}\right)}}\\ & \approx & \displaystyle \frac{{{cz}}_{i}}{{H}_{0}\left(1+{z}_{i}\right)}\left(1-\displaystyle \frac{3{{\rm{\Omega }}}_{m}{z}_{i}}{4}+\displaystyle \frac{{{\rm{\Omega }}}_{m}(9{{\rm{\Omega }}}_{m}-4){z}_{i}^{2}}{8}\right),\end{array}\end{eqnarray} \tag{ 1 }$$

where by "expected" here we refer to the values that the distance and redshift would take if the galaxy were perfectly following the Hubble flow. Equation (1) assumes a flat ΛCDM cosmology, and the series expansion is accurate to one part in ∼10⁵ for the range of redshifts covered by our observations. We set Ω_m = 0.315, from Planck Collaboration et al. (2018), though we note that any choice in the range 0 ≤ Ω_m ≤ 0.5 would yield a distance that differs by ≲1% for the galaxies in our sample. For each of the H₀-fitting approaches described in this section, the contribution to the likelihood from the distance constraints is given by the product of the posterior distributions ${ \mathcal P }({\hat{D}}_{i}| {D}_{i})$ from the independent disk fits,

$$\begin{eqnarray}&&{{ \mathcal L }}_{D}=\prod _{i}{ \mathcal P }({\hat{D}}_{i}| {D}_{i}),\end{eqnarray} \tag{ 2 }$$

where ${\hat{D}}_{i}$ is the distance measured from disk modeling.

The expected cosmological recession velocity v_i = cz_i differs from the measured galaxy velocity ${\hat{v}}_{i}$ both because of statistical uncertainty in the measurement and because of a systematic uncertainty in the form of peculiar motion. For our measurements, in which the statistical uncertainties in velocity are quite small (typically ∼1–2 km s⁻¹), peculiar motions dominate the recession velocity uncertainty. Because the galaxies in our sample all reside at low redshifts (z ≪ 1), we proceed under the assumption that peculiar velocities are independent of redshift.

When fitting for H₀ we have several options for treating these peculiar velocities, and in this section we describe the approaches we have taken to construct the velocity contribution to the likelihood, ${{ \mathcal L }}_{v}$. In all cases, our combined likelihood ${ \mathcal L }$ is ultimately given by the product of the velocity and distance likelihoods,

$$\begin{eqnarray}&&{ \mathcal L }={{ \mathcal L }}_{D}{{ \mathcal L }}_{v},\end{eqnarray} \tag{ 3 }$$

and the posterior distribution is given by the product of ${ \mathcal L }$ with the prior via Bayes's theorem. We assume flat priors for all model parameters, and we explore the posterior space using the dynesty nested sampling package (Speagle 2020). The right panel of Figure 1 shows the posterior distributions for all model fits.

Zoom In Zoom Out Reset image size

3.1. Treating Peculiar Velocities as Inflated Measurement Uncertainties

The simplest way to take peculiar velocities into account is to incorporate them into the velocity measurement uncertainties. Typical values for galaxy peculiar velocities lie in the range ∼150–250 km s⁻¹ (e.g., Davis et al. 1997; Zaroubi et al. 2001; Masters et al. 2006; Hoffman et al. 2015), so we conservatively take the upper end of this range and add σ_pec = 250 km s⁻¹ in quadrature to our velocity uncertainties. The velocity contribution to the likelihood is then given by a Gaussian distribution,

$$\begin{eqnarray}&&{{ \mathcal L }}_{v}=\prod _{i}\displaystyle \frac{1}{\sqrt{2\pi \left({\sigma }_{v,i}^{2}+{\sigma }_{\mathrm{pec}}^{2}\right)}}\exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\left({v}_{i}-{\hat{v}}_{i}\right)}^{2}}{{\sigma }_{v,i}^{2}+{\sigma }_{\mathrm{pec}}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

where ${\sigma }_{v,i}$ is the statistical uncertainty in velocity measurement ${\hat{v}}_{i}$ and the true velocities v_i are treated as nuisance parameters in the model.

The result from fitting this model to all six maser galaxies simultaneously is H₀ = 73.9 ± 3.0 km s⁻¹ Mpc⁻¹, and Table 2 lists the values obtained from leave-one-out jackknife tests. We assess the goodness-of-fit using a chi-squared statistic,

$$\begin{eqnarray}&&{\chi }_{\nu }^{2}=\displaystyle \frac{1}{\nu }\sum _{i}\left[\displaystyle \frac{{\left({v}_{i}-{\hat{v}}_{i}\right)}^{2}}{{\sigma }_{v,i}^{2}+{\sigma }_{\mathrm{pec}}^{2}}+\displaystyle \frac{{\left({D}_{i}-{\hat{D}}_{i}\right)}^{2}}{{\sigma }_{D,i}^{2}}\right],\end{eqnarray} \tag{ 5 }$$

where ${\sigma }_{D,i}$ is the standard deviation of the distance measurement posterior and ${\chi }_{\nu }^{2}$ is the chi-squared per degree of freedom, ν. For this fit, ν = 5 and ${\chi }_{\nu }^{2}=0.6$, which is consistent with unity within the expected standard deviation of a chi-squared distribution with five degrees of freedom (see also Table 3).

Table 2.  Hubble Constant Constraints and Jackknife Tests

Peculiar Velocity Treatment	Galaxies Excluded from the Fit	H0 (km s−1 Mpc−1)
(1) Assign a fixed velocity uncertainty of 250 km s−1	UGC 3789	${75.8}_{-3.3}^{+3.4}$
	NGC 6264	${73.8}_{-3.2}^{+3.2}$
	NGC 6323	${73.8}_{-3.0}^{+3.1}$
	NGC 5765b	${74.1}_{-4.4}^{+4.5}$
	CGCG 074-064	${72.5}_{-3.2}^{+3.4}$
	NGC 4258	${73.6}_{-3.0}^{+3.1}$
	Fit using all galaxies:	${73.9}_{-3.0}^{+3.0}$
(2) Fit for σpec using the maser data and assuming an outlier-robust form for the peculiar velocity distribution	UGC 3789	${76.4}_{-3.8}^{+4.2}$
	NGC 6264	${74.4}_{-3.8}^{+4.4}$
	NGC 6323	${74.5}_{-3.6}^{+4.0}$
	NGC 5765b	${75.8}_{-5.6}^{+6.6}$
	CGCG 074-064	${73.1}_{-3.9}^{+4.3}$
	NGC 4258	${74.2}_{-3.7}^{+4.5}$
	Fit using all galaxies:	${74.4}_{-3.4}^{+3.9}$
(3) Use galaxy group recession velocities	UGC 3789	${75.0}_{-3.0}^{+3.1}$
	NGC 6264	${73.1}_{-2.7}^{+2.8}$
	NGC 6323	${73.2}_{-2.7}^{+2.8}$
	NGC 5765b	${72.2}_{-4.1}^{+4.2}$
	CGCG 074-064	${73.3}_{-3.0}^{+3.1}$
	NGC 4258	${72.8}_{-2.7}^{+2.8}$
	Fit using all galaxies:	${73.3}_{-2.7}^{+2.8}$
(4) Use 2M++ (Carrick et al. 2015) recession velocities	UGC 3789	${73.3}_{-3.0}^{+3.0}$
	NGC 6264	${71.8}_{-2.8}^{+2.8}$
	NGC 6323	${71.9}_{-2.7}^{+2.8}$
	NGC 5765b	${71.1}_{-3.9}^{+4.0}$
	CGCG 074-064	${70.9}_{-2.9}^{+3.0}$
	NGC 4258	${72.1}_{-2.7}^{+2.7}$
	Fit using all galaxies:	${71.8}_{-2.7}^{+2.7}$
(5) Use CF3 (Graziani et al. 2019) recession velocities	UGC 3789	${73.6}_{-2.9}^{+3.1}$
	NGC 6264	${71.5}_{-2.7}^{+2.8}$
	NGC 6323	${71.7}_{-2.6}^{+2.8}$
	NGC 5765b	${71.5}_{-4.0}^{+4.1}$
	CGCG 074-064	${70.5}_{-2.9}^{+3.0}$
	NGC 4258	${72.0}_{-2.7}^{+2.7}$
	Fit using all galaxies:	${71.8}_{-2.6}^{+2.7}$
(6) Use M2000 (Mould et al. 2000) recession velocities	UGC 3789	${79.3}_{-3.1}^{+3.3}$
	NGC 6264	${76.8}_{-2.9}^{+2.9}$
	NGC 6323	${76.9}_{-2.9}^{+2.9}$
	NGC 5765b	${76.2}_{-4.1}^{+4.3}$
	CGCG 074-064	${75.5}_{-3.0}^{+3.2}$
	NGC 4258	${76.8}_{-2.9}^{+2.9}$
	Fit using all galaxies:	${76.9}_{-2.9}^{+2.9}$

Note. Hubble constant measurements made using various subsets of the megamaser distances and different treatments for the peculiar velocities, as described in Section 3. For each peculiar velocity treatment, we list seven H₀ values: six of these values correspond to "leave-one-out" jackknife tests, in which we fit the data (under the given peculiar velocity prescription) after removing the galaxy specified in the second column; the seventh value corresponds to that obtained from fitting all galaxies simultaneously. For each fit we quote the posterior median and 1σ confidence interval (i.e., 16th to 84th percentile).

Download table as:  ASCIITypeset image

Table 3.  H₀ Goodnesses-of-fit and Comparisons with Other Measurements

Peculiar Velocity Treatment	${\chi }_{\nu }^{2}$	$P({H}_{0}\leqslant {H}_{0,\mathrm{Planck}})$	$P({H}_{0}\geqslant {H}_{0,\mathrm{SH}0\mathrm{ES}})$
(1)	0.60	0.02	0.48
(2)	1.52	0.03	0.55
(3)	0.62	0.01	0.41
(4)	0.55	0.05	0.24
(5)	0.75	0.05	0.23
(6)	0.75	<0.01	0.82

Note. Statistics for the Hubble constant fits described in Section 3, with the different peculiar velocity treatments numbered as in Table 2. The second column lists the chi-squared per degree of freedom for each fit, computed using Equation (5). For treatments (1), (3), (4), (5), and (6), the number of degrees of freedom ν = 5 and the expected standard deviation of the ${\chi }_{\nu }^{2}$-distribution is $\sqrt{2/5}\approx 0.63$, while for treatment (2), ν = 4 and the expected standard deviation in ${\chi }_{\nu }^{2}$ is $\sqrt{2/4}\approx 0.71$. The third column lists one-sided comparison statistics computed using Equation (10), which give the probability that our H₀ measurement is at least as low as the Planck measurement (Planck Collaboration et al. 2018). The fourth column is analogous to the third, and lists the probability that our H₀ measurement is at least as high as the SH0ES measurement (Riess et al. 2019; the statistic is computed using Equation (11)).

Download table as:  ASCIITypeset image

3.2. Modeling Peculiar Velocities as Being Drawn from a Global Distribution

Rather than assuming a typical dispersion for the peculiar velocity distribution of σ_pec, we can instead fit for it as part of the model given some assumption about the form of the underlying distribution from which peculiar velocities are drawn. Because we expect that some unknown fraction of galaxies may have particularly large peculiar velocities (e.g., if the galaxy lives in a cluster), we test the (Sivia & Skilling 2006, Section 8.3.1) "conservative formulation" for the velocity uncertainties as an alternative to a Gaussian distribution. Under this formalism, σ_pec is interpreted as a lower bound on the velocity error σ associated with peculiar velocities, with a distribution for this error given by

$$\begin{eqnarray}{ \mathcal P }(\sigma | {\sigma }_{\mathrm{pec}})=\left\{\begin{array}{ll}\displaystyle \frac{{\sigma }_{\mathrm{pec}}}{{\sigma }^{2}}, & \sigma \gt {\sigma }_{\mathrm{pec}}\\ 0, & \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ 6 }$$

The marginal likelihood contribution from the velocity constraints after integrating out σ is then

$$\begin{eqnarray}&&{{ \mathcal L }}_{v}=\prod _{i}\displaystyle \frac{1}{\sqrt{2\pi \left({\sigma }_{v,i}^{2}+{\sigma }_{\mathrm{pec}}^{2}\right)}}\left(\displaystyle \frac{1-{e}^{-{R}_{i}^{2}/2}}{{R}_{i}^{2}}\right),\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{R}_{i}=\displaystyle \frac{{v}_{\mathrm{pec},i}^{2}}{{\sigma }_{v,i}^{2}+{\sigma }_{\mathrm{pec}}^{2}},\end{eqnarray} \tag{ 8 }$$

as we have assumed that the peculiar velocity distribution has zero mean.

Because we are now modeling peculiar velocity rather than recession velocity, we have to relate the two. We relate the expected cosmological recession redshift z_i to the measured redshift ${\hat{z}}_{i}$ by adding the peculiar velocity contribution via (see, e.g., Davis & Scrimgeour 2014),

$$\begin{eqnarray}&&1+{\hat{z}}_{i}=\left(1+{z}_{i}\right)\left(1+\displaystyle \frac{{v}_{\mathrm{pec},i}}{c}\right).\end{eqnarray} \tag{ 9 }$$

The z_i values are then plugged into Equation (1) to compute the distances D_i.

The result from fitting this model to all six maser galaxies simultaneously is ${H}_{0}={74.4}_{-3.4}^{+3.9}$ km s⁻¹ Mpc⁻¹, and Table 2 lists the values obtained from leave-one-out jackknife tests. We find ${\sigma }_{\mathrm{pec}}={141}_{-80}^{+185}$ km s⁻¹ and use it via Equation (5) to compute ${\chi }_{\nu }^{2}=1.5$ (see Table 3), which is consistent with unity within the expected standard deviation for a chi-squared distribution with ν = 4 degrees of freedom. We note also that the explicitly non-Gaussian form of the likelihood in Equation (7)—in particular the $1/{R}_{i}^{2}$ tails of this distribution—are expected to drive the chi-squared above unity.

3.3. Using Galaxy Group Velocities in Place of Individual Galaxy Velocities

3.4. Using Galaxy Flow Models to Correct for Peculiar Velocities

We have applied an improved approach for fitting maser data to obtain more precise distance estimates to four previously published MCP galaxies: UGC 3789, NGC 6264, NGC 6323, and NGC 5765b. We find that previous maser disk modeling efforts have typically overestimated the systematic measurement uncertainties associated with maser positions from VLBI maps. By incorporating these error floors into our disk model as model parameters, we are able to fit for and then marginalize over them, eliminating a source of systematic uncertainty and generically improving the measurement precision. We have combined the revised distance estimates for UGC 3789, NGC 6264, NGC 6323, and NGC 5765b with the recently published distances to CGCG 074-064 and NGC 4258—both of which included error floors as parameters in the model—to derive constraints on H₀. Assuming a global velocity uncertainty of 250 km s⁻¹ associated with peculiar motions, we find H₀ = 73.9 ± 3.0 km s⁻¹ Mpc⁻¹. A Hubble diagram is shown in the left panel of Figure 1.

Our fiducial H₀ measurement is determined exclusively using megamaser-based distance and velocity measurements, and it thus represents an independent cosmological probe from standard candles, gravitational lenses, and the CMB. The primary source of systematic uncertainty in this measurement comes from the unknown peculiar motions of the maser galaxies, and we have considered three different treatments¹⁸ for determining how these peculiar velocities could modify the H₀ value:

• 1.  
  We permit the global velocity uncertainty to be a free parameter that we fit alongside H₀.
• 2.  
  We replace each galaxy's recession velocity with the velocity of the group that galaxy is a member of.
• 3.  
  We replace each galaxy's recession velocity with the velocity predicted by a galaxy flow model evaluated at that galaxy's location and redshift. We use the 2M++, CF3, and M2000 peculiar velocity models.

The first two of the above treatments modify the best-fit H₀ value by less than 1 km s⁻¹ Mpc⁻¹, though the measurement precision suffers when permitting the global velocity uncertainty to be a free parameter because of the reduced degrees of freedom. Using recession velocities from either the 2M++ or CF3 galaxy flow models systematically reduces the best-fit H₀ value by 2.1 km s⁻¹ Mpc⁻¹, and the measurement precision improves because of the smaller uncertainty (150 km s⁻¹) associated with the catalog velocities. The recession velocities from the M2000 model, on the other hand, result in a substantially increased best-fit H₀ value, which is larger by 3.0 km s⁻¹ Mpc⁻¹ than the fiducial measurement and by 5.1 km s⁻¹ Mpc⁻¹ than the measurement made from correcting for peculiar velocities using the other two flow models.

We have also performed a series of leave-one-out jackknife tests for each of the different peculiar velocity treatments. We find that the removal of any single galaxy from the sample never modifies the best-fit H₀ value by more than 1σ, indicating that our measurement is not being unduly influenced by a single outlying value.

We test the prior empirical claim that the local value of H₀ exceeds the early-universe value (e.g., Riess et al. 2019; Verde et al. 2019; Wong et al. 2019) by calculating the probability of rejecting the null hypothesis that our measurement does not exceed Planck's, i.e.,

$$\begin{eqnarray}&&P({H}_{0}\leqslant {H}_{0,\mathrm{Planck}})={\int }_{-\infty }^{\infty }{ \mathcal P }(H)\left[{\int }_{-\infty }^{H}{ \mathcal P }({H}_{0}^{{\prime} }){{dH}}_{0}^{{\prime} }\right]{dH},\end{eqnarray} \tag{ 10 }$$

where ${ \mathcal P }({H}_{0})$ is our measured posterior distribution for H₀ and we treat the Planck measurement probability distribution ${ \mathcal P }(H)$ as a Gaussian with mean and standard deviation given by the published measurement of ${H}_{0,\mathrm{Planck}}=67.4\pm 0.5$ km s⁻¹ Mpc⁻¹ (Planck Collaboration et al. 2018). The value of $P({H}_{0}\leqslant {H}_{0,\mathrm{Planck}})$ is listed for all four peculiar velocity treatments in Table 3. The first treatment, our default, gives a 2% chance that our value is lower than Planck's, corroborating the sense of the present tension in H₀ at 98% confidence. The other peculiar velocity treatments give confidences of 95%–99%. Performing an analogous comparison with a late-universe measurement from SH0ES of ${H}_{0,\mathrm{SH}0\mathrm{ES}}=74.03\pm 1.42$ km s⁻¹ Mpc⁻¹ (Riess et al. 2019),

$$\begin{eqnarray}&&P({H}_{0}\geqslant {H}_{0,\mathrm{SH}0\mathrm{ES}})={\int }_{-\infty }^{\infty }{ \mathcal P }(H)\left[{\int }_{H}^{\infty }{ \mathcal P }({H}_{0}^{{\prime} }){{dH}}_{0}^{{\prime} }\right]{dH},\end{eqnarray} \tag{ 11 }$$

we find that our result is consistent with the SH0ES measurement, with little preference for a higher or lower value (see Table 3).

The ∼4% H₀ constraint presented in this paper comes from consideration of only six megamaser-hosting galaxies, and the precision is ultimately limited by the quality and quantity of the available distance measurements. Future H₀ measurements from the MCP will improve on this precision by incorporating distance measurements from additional megamaser-hosting galaxies.

We thank R. B. Tully and R. Graziani for help with the Cosmicflows-3 peculiar velocities, M. J. Hudson for advice regarding the 2M++ peculiar velocities, L. Blackburn for modeling discussions, and Erik Peterson for help with understanding redshift uncertainties. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work made use of the Swinburne University of Technology software correlator, developed as part of the Australian Major National Research Facilities Programme and operated under license. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work was supported in part by the black hole Initiative at Harvard University, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation to Harvard University.

Facilities: GBT - Green Bank Telescope, VLA - , VLBA - , Effelsberg.

Software: AIPS, CASA, GBTIDL, PyMC3 (Salvatier et al. 2016), dynesty (Speagle 2020).

In Table 5 we list the values for all model parameters from the updated disk fits to UGC 3789, NGC 6264, NGC 6323, and NGC 5765b. A comprehensive description of the model is given in Pesce et al. (2020).

Table 5.  Updated Disk-fitting Results

Parameter	Units	Galaxy
		UGC 3789	NGC 6264	NGC 6323	NGC 5765b
D	Mpc	${51.5}_{-4.0}^{+4.5}$	${132.1}_{-17.3}^{+21.2}$	${109.4}_{-23.4}^{+34.2}$	${112.2}_{-5.1}^{+5.4}$
${M}_{\mathrm{BH}}$	107 M⊙	${1.19}_{-0.09}^{+0.10}$	${2.76}_{-0.36}^{+0.45}$	${1.02}_{-0.22}^{+0.32}$	${4.15}_{-0.19}^{+0.20}$
v	km s−1	3319.9 ± 0.8	10192.6 ± 0.8	7801.5 ± 1.5	8525.7 ± 0.7
${x}_{0}$	mas	−0.4014 ±  0.0010	0.0050 ± 0.0012	0.0161 ±  0.0010	−0.0440 ±  0.0014
${y}_{0}$	mas	−0.4615 ±  0.0011	0.0076 ± 0.0016	0.0073 ±  0.0024	−0.0995 ±  0.0019
${i}_{0}$	degree	84.9 ± 0.6	91.3 ± 2.3	91.5 ± 0.3	72.4 ± 0.5
$\displaystyle \frac{{di}}{{dr}}$	degree mas−1	7.7 ± 1.0	−1.5 ± 4.0	...	12.5 ± 0.5
${{\rm{\Omega }}}_{0}$	degree	222.4 ± 0.4	84.7 ± 1.3	184.4 ± 0.6	149.7 ± 0.3
$\displaystyle \frac{d{\rm{\Omega }}}{{dr}}$	degree mas−1	−1.6 ± 0.6	16.7 ± 2.2	12.9 ± 1.2	−3.2 ± 0.2
${\sigma }_{x}$	mas	0.0045 ±  0.0011	${0.0012}_{-0.0008}^{+0.0011}$	0.0035 ± 0.0010	${0.0028}_{-0.0011}^{+0.0012}$
${\sigma }_{y}$	mas	0.0063 ±  0.0013	0.0056 ± 0.0020	${0.0037}_{-0.0022}^{+0.0025}$	0.0034 ± 0.0009
${\sigma }_{v,\mathrm{sys}}$	km s−1	${1.7}_{-0.8}^{+0.9}$	${1.6}_{-0.8}^{+0.9}$	${2.1}_{-0.9}^{+1.0}$	1.1 ± 0.5
${\sigma }_{v,\mathrm{hv}}$	km s−1	${1.8}_{-0.7}^{+0.9}$	${1.3}_{-0.6}^{+0.7}$	${1.9}_{-0.7}^{+0.9}$	${1.5}_{-0.5}^{+0.6}$
${\sigma }_{a}$	km s−1 yr−1	${0.34}_{-0.05}^{+0.06}$	${0.08}_{-0.05}^{+0.07}$	0.21 ± 0.09	0.041 ± 0.014

Note. Top: fitting results for the global parameters describing the maser disk, marginalized over all other parameters; for each value we quote the posterior median and 1σ confidence interval (i.e., 16th to 84th percentile). Here, D is the angular-size distance to the galaxy, M_BH is the mass of the SMBH, v is the line-of-sight CMB-frame velocity of the SMBH, (x₀, y₀) is the coordinate location of the SMBH, i₀ is the inclination angle of the disk at r = 0, $\tfrac{{di}}{{dr}}$ is the first-order inclination angle warping parameter, Ω₀ is the position angle of the disk at r = 0, and $\tfrac{d{\rm{\Omega }}}{{dr}}$ is the first-order position angle warping parameter. For the uncertainties we quote 1σ (i.e., 16% and 84%) confidence intervals from the posteriors. The SMBH coordinate locations are referenced to the coordinate zero-point used in the respective data paper: for UGC 3789, see Reid et al. (2009); for NGC 6264, see Kuo et al. (2013); for NGC 6323, see Kuo et al. (2015); and for NGC 5765b, see Gao et al. (2017). Bottom: fitting results for the error floor parameters; σ_x is the x-position error floor, σ_y is the y-position error floor, σ_v,sys is the error floor for the systemic feature velocities, σ_v,hv is the error floor for the high-velocity feature velocities, and σ_a is the acceleration error floor.

Download table as:  ASCIITypeset image