THE MEGAMASER COSMOLOGY PROJECT. IV. A DIRECT MEASUREMENT OF THE HUBBLE CONSTANT FROM UGC 3789

Measurements of the expansion history of the universe, H(z), play a fundamental role in our understanding of cosmological evolution and its far-reaching physical implications, including addressing the nature of dark energy, the curvature of space, the masses of neutrinos, and the number of families of relativistic particles. While detailed measurements of the cosmic microwave background (CMB; e.g., Komatsu et al. 2011) allow us to measure linear scales at redshift z ≈ 1100, this mostly provides information at a single epoch and when the influence of dark energy was still negligible. Complementary information from later cosmic times is therefore essential. This information can come from Type Ia supernovae, galaxy clustering, gravitational lensing, and baryon acoustic oscillations (e.g., Riess et al. 2011; Bonamente et al. 2006; Suyu et al. 2010; Reid et al. 2010).

All of these data provide distances and linear scales at significant redshifts. However, it is the local universe where dark energy is dominant. Thus the local Hubble constant, H₀, provides the largest "lever arm" with respect to the CMB for constraining the time dependence of the equation of state of dark energy. The use of Cepheids (e.g., Freedman et al. 2001; Sandage et al. 2006) as well as the combined use of Cepheids and Type Ia supernovae (e.g., Riess et al. 2011) has traditionally dominated determinations of H₀. What is missing, however, are direct geometric distance estimates that do not require a complex and uncertain ladder of calibration of "standard candles."

Direct geometric distance measurements to water masers in nuclear regions of galaxies that are well into the Hubble flow (roughly >30 Mpc distant) provide a promising new and independent method for refining the value of H₀. Observations of water masers in accretion disks within ∼0.1 pc of a galaxy's central supermassive black hole have been used to measure angular-diameter distances to galaxies, independently of other techniques that often rely on standard candles. Very Long Baseline Array (VLBA) observations of the H₂O masers in the nearby Seyfert 2 galaxy NGC 4258 established the technique and provided an accurate, angular-diameter distance of 7.2 ± 0.5 Mpc to the galaxy (Herrnstein et al. 1999). This galaxy is too close to permit a direct measurement of H₀ (i.e., by dividing its recessional speed by its distance) since its uncertain peculiar velocity could be a large fraction of its recessional speed. However, NGC 4258 has proven extremely valuable as a solid anchor for the extragalactic distance scale (Freedman et al. 2001; Riess et al. 2011).

Recently, the Megamaser Cosmology Project (MCP) reported a distance to UGC 3789, another galaxy with water masers in a nuclear accretion disk, of 49.9 ± 7.0 Mpc (Reid et al. 2009; Braatz et al. 2010; hereafter Papers I and II). This galaxy has a recessional velocity of ≈3481 km s⁻¹ (relativistically corrected and referenced to the CMB), which includes its peculiar velocity of 151 ± 163 km s⁻¹ based on galaxy flow models of Masters et al. (2006) and Springob et al. (2007). Combining the recessional velocity and distance yielded H₀ = 69 ± 11 km s⁻¹ Mpc⁻¹.

Since UGC 3789 is well into the Hubble flow it can provide a direct estimate of H₀ with a potential uncertainty as small as ±5%, limited by the uncertainty in its peculiar motion. Therefore, since reporting our initial results, we have conducted additional observations of UGC 3789, in order to reduce the uncertainty of the H₀ estimate. In total, we have now analyzed nine very long baseline interferometric (VLBI) observations, using the NRAO⁵ 10 antenna VLBA, the 100 m Green Bank Telescope (GBT), and the 100 m Effelsberg⁶ telescope. The VLBI observations are reported in Section 2. We have also extended our monitoring of changes in the water maser spectrum with monthly observations (except during the humid summer) with the GBT for a period spanning 5.5 years. These spectra are documented in Section 3 and used to determine the accelerations of individual maser features. In Section 4, we describe a Bayesian approach for fitting a model of the accretion disk to these data. The model allows for a warped disk with eccentric gas orbits and includes parameters for the central (black hole) mass and position, as well as for H₀.

VLBI observations were conducted at nine epochs (NRAO program codes BB227A, BB227B, BB242A, BB242I, BB242K, BB242L, BB242Q, BB261G, and BB261S) between 2007 December and 2010 April. Six of the nine observations yielded maps with signal-to-noise ratios degraded by factors of two or greater owing to the loss of one of the 100 m telescopes (usually to poor weather) or the weakening of the peak maser emission to levels where self-calibration (phase-referencing) using the maser emission was poor. The three epochs with excellent weather, antenna performance, and strong maser emission included BB227A (2006 December 10; reported in Paper I), BB242L (2008 December 12), and BB261G (2009 April 11). Only results from these three epochs are reported here.

We observed with 16 MHz bands covering five frequencies, three in dual polarization and two in single polarization. In our early observations, these bands were centered at (optical definition) local standard of rest (LSR) velocities (and polarizations) of 3880.0 (LCP and RCP), 3710.0 (LCP), 3265.0 (LCP and RCP), 2670.0 (LCP and RCP), and 2500.0 (LCP) km s⁻¹ for BB227A; for the BB242L and BB261G, we shifted the center velocities of the fourth and fifth bands to 2717.6 and 2513.6 km s⁻¹, in order to map new maser features not covered in the original setup. The data for each polarization of each band were cross-correlated with 128 spectral channels, yielding channels separated by 1.7 km s⁻¹.

Generally, the data were analyzed as described in Paper I. The final calibration step involved selecting a maser feature as the interferometer phase reference, and details of this procedure varied among epochs depending on maser strength and interferometer coherence times. The strongest maser feature in the spectrum usually peaked at ≈0.07 Jy and was fairly broad. For BB227A, we averaged five spectral channels spanning an LSR velocity range of 2685–2692 km s⁻¹ (i.e., channels 52–56 from the blueshifted high-velocity band centered at V_LSR = 2670 km s⁻¹), adding together the data from both polarizations, and fitting fringes over a one-minute period. For BB242L, we averaged 13 spectral channels spanning an LSR velocity range of 2682–2702 km s⁻¹ (i.e., channels 74–86 from the blueshifted high-velocity band centered at V_LSR = 2717.6 km s⁻¹) and fitted fringes over a two-minute period. For BB261G, we averaged seven spectral channels spanning an LSR velocity range of 2685–2696 km s⁻¹ (i.e., channels 78–84 from the blueshifted high-velocity band centered at V_LSR = 2717.6 km s⁻¹), again with two-minute averaging.

After calibration, we Fourier-transformed the gridded (u,v)-data to make images of the maser emission in all spectral channels for each of the five IF bands. The images were deconvolved with the point-source response using the CLEAN algorithm and restored with a circular Gaussian beam with a 0.30 mas FWHM (approximately matching the geometric mean of the dirty beam). The rms noise levels in channel maps were ≈1 mJy. All images appeared to contain single, point-like maser spots. We then fitted each spectral-channel image with an elliptical Gaussian brightness distribution in order to obtain positions and flux densities.

Spectral images for the three epochs are shown in Figure 1 for maser spots stronger than 5 mJy. (Note that in Paper I, we displayed only maser spots stronger than 10 mJy.) These images show nearly identical spatial-velocity patterns of maser emission. This is expected as material in near-circular orbit in the accretion disk at the average masing radius has an orbital period of ∼1000 yr and thus rotates by less than 1° of disk azimuth between the first and the last observation. Most of the spatial scatter in the images is from measurement uncertainty, primarily from signal-to-noise limitations, of ≈25 μas for weaker features of ≈10 mJy.

Zoom In Zoom Out Reset image size

Paper II describes the procedures used to monitor the H₂O maser spectrum of UGC 3789 from 2006 January to 2009 March with the GBT. We continued these observations through 2011 June. Typically, observations were conducted at monthly intervals for nine months per year, avoiding the summer when atmospheric water vapor precluded sensitive observations. These spectra were divided into six yearly blocks and analyzed to determine velocity drifts (accelerations) of individual maser features.

Measurement of the acceleration of systemic-velocity masers is key to determining the angular-diameter distance of the galaxy. In principle, accelerations are straightforwardly obtained from the linear drift in peak velocity with time. In practice, however, blending of adjacent spectral features, coupled with modest signal-to-noise ratio observations, makes this difficult. We applied two methods described below to measure accelerations: method 1 to obtain the best acceleration estimates and method 2 to check for possible fitting biases owing to initial parameter values.

3.1. Method 1

3.2. Method 2

3.3. Acceleration Fitting Results

For systemic-velocity features, we used three velocity windows for acceleration fitting: 3230–3280, 3286–3316, and 3350–3370 km s⁻¹, since emission outside of these windows was generally absent or very weak (<5 mJy). Systemic-velocity maser features in UGC 3789 typically have lifetimes (with flux densities >5 mJy) of ∼12 months. Therefore, we fitted accelerations separately to groups of ≈9 consecutive monthly spectra (from fall through spring) covering 5.5 years and centered at ≈2008.8. These observations covered well the three high-quality VLBI imaging observations discussed in this paper, whose mean epoch was ≈2008.4.

Yearly acceleration measurements generally were consistent with a single value, although for some velocity ranges there was considerable scatter (up to about ±30% for masers near 3290 km s⁻¹). Because of this, and the expectation of small changes in acceleration over our observing period (as masing clouds move by less than 1° of disk azimuth), we velocity-binned and -averaged the fitted accelerations from the six "yearly" groups of spectra. Variance-weighted averages of acceleration as a function of velocity for the systemic features are plotted Figure 2 for both fitting methods outlined above.

Zoom In Zoom Out Reset image size

The results for method 1 are plotted with filled (red) circles in Figure 2. There is moderate scatter among the accelerations within each velocity window, and this scatter greatly exceeds that expected from the formal acceleration uncertainties. The greatest variations occur in the velocity range 3279–3330 km s⁻¹; we can bound this variation by assuming a constant acceleration over this velocity range and calculating a standard deviation of 1.6 km s⁻¹ yr⁻¹. We suspect that some of the variation is real and intrinsic to the source, but also that some of the variations may originate in the fitting process, owing to multiple velocity-blended components and using data with only moderate signal-to-noise ratios.

In order to provide an independent check on the results from method 1, we re-fitted the spectra using method 2. For each velocity window and each "yearly" time group, we selected the four trial fits with the lowest χ² values. We then binned all results in 2 km s⁻¹ bins and calculated a weighted average acceleration for each velocity bin. These results are plotted in Figure 2 with open (blue) squares. There is very good agreement between the two methods; the differences in average accelerations in each velocity window are small; unweighted means for methods 1 and 2 are 1.7 and 1.7 (±0.2) km s⁻¹ yr⁻¹ for the velocity range 3230–3278 km s⁻¹, 5.3 and 5.9 (±0.4) km s⁻¹ yr⁻¹ for the velocity range 3279–3330 km s⁻¹, and 2.7 and 3.0 (±0.3) km s⁻¹ yr⁻¹ for the velocity range 3350–3370 km s⁻¹. The small differences between the slopes of A versus V_LSR for the two methods are within statistical uncertainties after accounting for the high correlation among accelerations within a window for method 2. We conclude that the accelerations given by the solid (red) points in Figure 2 are near optimum values and are not significantly biased by the choice of initial parameter values in the first step of method 1.

Acceleration measurements for high-velocity features are considerably less complicated as these features are not highly blended and have small changes in velocity over time. Therefore, we only used method 1 (described in Section 3.1) for these features. Fitted accelerations are given in Table 1.

Table 1. UGC 3789 H₂O Maser Data

VLSR	Θx	$\sigma _{\Theta _x}$	Θy	$\sigma _{\Theta _y}$	A	σA
(km s−1)	(mas)	(mas)	(mas)	(mas)	(km s−1 yr−1)	(km s−1 yr−1)
2472.8	−0.183	0.009	−0.190	0.009	...	...
2477.9	−0.179	0.008	−0.197	0.021	...	...
2562.9	−0.116	0.005	−0.176	0.028	...	...
2564.6	−0.119	0.004	−0.153	0.005	...	...
2566.3	−0.122	0.003	−0.149	0.003	...	...
2568.0	−0.123	0.003	−0.143	0.009	0.5	0.2
2569.7	−0.130	0.009	−0.137	0.007	...	...
2571.4	−0.093	0.020	−0.161	0.008	...	...
2586.7	−0.104	0.012	−0.137	0.012	...	...
2588.4	−0.089	0.004	−0.118	0.003	0.0	0.2
2590.1	−0.084	0.009	−0.102	0.005	...	...
2591.8	−0.087	0.008	−0.103	0.006	...	...
2593.5	−0.096	0.004	−0.107	0.004	0.0	0.2
2595.2	−0.098	0.004	−0.097	0.007	...	...
2596.9	−0.105	0.007	−0.107	0.015	...	...
2602.0	−0.066	0.010	−0.106	0.012	...	...
2603.7	−0.084	0.003	−0.102	0.009	...	...
2605.4	−0.076	0.007	−0.102	0.004	...	...
2607.1	−0.065	0.005	−0.102	0.004	...	...
2608.8	−0.067	0.015	−0.096	0.004	...	...
2610.5	−0.073	0.005	−0.089	0.007	0.0	0.2
2612.2	−0.077	0.003	−0.087	0.007	0.0	0.2
2613.9	−0.081	0.006	−0.087	0.007	...	...
2615.6	−0.068	0.007	−0.084	0.015	0.1	0.2
2617.3	−0.072	0.007	−0.085	0.007	...	...
2629.2	−0.045	0.012	−0.079	0.013	...	...
2630.9	−0.042	0.009	−0.054	0.010	...	...
2632.6	−0.049	0.004	−0.059	0.004	...	...
2634.3	−0.050	0.004	−0.077	0.004	...	...
2636.0	−0.064	0.008	−0.063	0.005	...	...
2668.3	−0.009	0.006	−0.034	0.007	...	...
2670.0	−0.020	0.006	−0.032	0.005	...	...
2678.5	−0.006	0.004	−0.024	0.004	...	...
2680.2	−0.016	0.003	−0.011	0.003	−0.1	0.2
2681.9	−0.008	0.002	−0.008	0.002	...	...
2683.6	−0.008	0.002	−0.008	0.002	0.0	0.2
2685.3	−0.005	0.002	−0.007	0.002	0.1	0.2
2687.0	0.000	0.002	0.000	0.002	...	...
2688.7	0.000	0.002	0.001	0.002	−0.1	0.2
2690.4	0.005	0.002	0.001	0.002	...	...
2692.1	0.000	0.003	0.005	0.002	−0.1	0.2
2693.8	0.009	0.002	0.006	0.002	0.0	0.2
2695.5	0.009	0.002	0.010	0.003	...	...
2697.2	0.010	0.003	0.007	0.004	−0.1	0.2
2698.9	0.014	0.002	0.010	0.002	0.1	0.2
2700.6	0.012	0.002	0.011	0.003	0.1	0.2
2702.3	0.011	0.002	0.011	0.003	−0.1	0.2
2704.0	0.013	0.003	0.011	0.004	−0.2	0.2
2705.7	0.019	0.003	0.026	0.003	...	...
2707.4	0.022	0.004	0.028	0.010	−0.3	0.2
2709.1	0.026	0.004	0.016	0.014	...	...
2710.8	0.036	0.007	0.033	0.007	−0.4	0.2
2712.5	0.052	0.004	0.029	0.007	...	...
2714.2	0.051	0.004	0.048	0.005	...	...
2715.9	0.036	0.014	0.046	0.005	...	...
2717.6	0.044	0.011	0.048	0.005	...	...
2719.3	0.062	0.006	0.054	0.009	...	...
2721.0	0.077	0.005	0.060	0.006	0.0	0.2
2734.6	0.046	0.003	0.045	0.005	...	...
2736.3	0.040	0.006	0.042	0.012	0.2	0.2
2738.0	0.045	0.008	0.054	0.006	...	...
2739.7	0.044	0.008	0.042	0.003	−0.1	0.2
2741.4	0.050	0.004	0.051	0.004	−0.1	0.2
2743.1	0.049	0.005	0.053	0.004	−0.1	0.2
2744.8	0.053	0.003	0.048	0.006	−0.2	0.2
2746.5	0.048	0.009	0.047	0.003	...	...
2748.2	0.058	0.004	0.050	0.005	−0.4	0.2
2749.9	0.058	0.004	0.052	0.014	...	...
3234.4	−0.403	0.010	−0.406	0.010	1.6	0.2
3239.5	−0.401	0.012	−0.429	0.005	1.3	0.2
3241.2	−0.404	0.005	−0.418	0.005	1.8	0.2
3242.9	−0.393	0.004	−0.433	0.004	2.5	0.2
3244.6	−0.389	0.003	−0.430	0.003	1.1	0.2
3246.3	−0.397	0.003	−0.433	0.004	1.9	0.2
3248.0	−0.397	0.004	−0.431	0.009	1.8	0.2
3249.7	−0.405	0.004	−0.438	0.008	1.6	0.2
3251.4	−0.405	0.003	−0.439	0.005	1.2	0.2
3253.1	−0.401	0.004	−0.437	0.005	1.0	0.2
3254.8	−0.413	0.007	−0.414	0.009	1.6	0.2
3256.5	−0.422	0.006	−0.442	0.011	0.9	0.2
3258.2	−0.411	0.004	−0.443	0.004	1.1	0.2
3259.9	−0.407	0.003	−0.442	0.008	1.5	0.2
3261.6	−0.406	0.004	−0.445	0.005	1.8	0.2
3263.3	−0.403	0.003	−0.463	0.015	1.9	0.2
3266.7	−0.401	0.009	−0.454	0.009	1.7	0.2
3268.4	−0.412	0.011	−0.484	0.007	1.9	0.2
3270.1	−0.402	0.006	−0.470	0.015	2.1	0.2
3271.8	−0.403	0.006	−0.478	0.006	1.9	0.2
3273.5	−0.404	0.006	−0.477	0.006	4.2	0.2
3287.1	−0.410	0.017	−0.482	0.008	2.5	0.8
3288.8	−0.406	0.004	−0.474	0.004	6.6	0.8
3290.5	−0.404	0.005	−0.506	0.005	6.1	0.8
3292.2	−0.409	0.003	−0.485	0.009	7.3	0.8
3293.9	−0.404	0.003	−0.487	0.003	8.4	0.8
3295.6	−0.403	0.003	−0.490	0.004	7.4	0.8
3297.3	−0.403	0.002	−0.492	0.006	7.2	0.8
3299.0	−0.405	0.002	−0.489	0.003	5.9	0.8
3300.7	−0.406	0.002	−0.489	0.004	5.7	0.8
3302.4	−0.409	0.002	−0.490	0.004	4.8	0.8
3304.1	−0.406	0.002	−0.493	0.002	4.7	0.8
3305.8	−0.404	0.002	−0.495	0.005	5.0	0.8
3307.5	−0.406	0.003	−0.494	0.005	5.4	0.8
3309.2	−0.404	0.002	−0.491	0.004	4.9	0.8
3310.9	−0.411	0.003	−0.493	0.003	5.9	0.8
3312.6	−0.410	0.003	−0.499	0.004	6.0	0.8
3314.3	−0.399	0.005	−0.503	0.005	5.3	0.8
3316.0	−0.412	0.010	−0.502	0.004	5.3	0.8
3317.7	−0.405	0.019	−0.499	0.005	5.5	0.8
3319.4	−0.418	0.015	−0.518	0.028	5.9	0.8
3358.5	−0.437	0.007	−0.511	0.006	2.4	0.2
3660.7	−1.198	0.010	−1.517	0.012	0.2	0.2
3735.5	−0.977	0.005	−1.139	0.008	...	...
3737.2	−0.980	0.006	−1.133	0.004	0.1	0.2
3738.9	−0.992	0.012	−1.119	0.011	...	...
3740.6	−0.971	0.005	−1.157	0.005	0.2	0.2
3745.7	−0.965	0.014	−1.143	0.020	0.2	0.2
3747.4	−0.973	0.006	−1.123	0.009	...	...
3749.1	−0.981	0.009	−1.109	0.005	−0.1	0.2
3750.8	−0.988	0.007	−1.095	0.007	...	...
3752.5	−0.965	0.011	−1.111	0.005	...	...
3754.2	−0.968	0.015	−1.094	0.008	−0.2	0.2
3755.9	−0.940	0.008	−1.092	0.004	...	...
3757.6	−0.953	0.006	−1.092	0.005	0.0	0.2
3759.3	−0.945	0.003	−1.104	0.004	...	...
3761.0	−0.944	0.004	−1.101	0.004	0.1	0.2
3764.4	−0.944	0.014	−1.082	0.015	...	...
3808.6	−0.782	0.006	−0.833	0.006	...	...
3837.5	−0.830	0.010	−0.917	0.028	...	...
3839.2	−0.811	0.009	−0.905	0.011	0.2	0.2
3840.9	−0.814	0.010	−0.914	0.005	...	...
3842.6	−0.832	0.006	−0.916	0.007	...	...
3846.0	−0.766	0.012	−0.874	0.014	...	...
3851.1	−0.791	0.007	−0.858	0.029	...	...
3859.6	−0.760	0.004	−0.848	0.005	...	...
3861.3	−0.764	0.004	−0.866	0.004	...	...
3863.0	−0.764	0.004	−0.882	0.004	...	...
3864.7	−0.760	0.007	−0.874	0.005	...	...
3866.4	−0.761	0.009	−0.868	0.007	0.3	0.2
3868.1	−0.765	0.010	−0.870	0.004	...	...
3871.5	−0.751	0.030	−0.874	0.021	...	...
3878.3	−0.761	0.016	−0.845	0.023	...	...
3880.0	−0.758	0.004	−0.871	0.004	...	...
3881.7	−0.750	0.004	−0.863	0.003	0.5	0.2
3883.4	−0.753	0.003	−0.854	0.009	...	...
3885.1	−0.746	0.003	−0.854	0.013	0.0	0.2
3886.8	−0.750	0.007	−0.865	0.005	...	...
3890.2	−0.758	0.024	−0.849	0.017	...	...
3891.9	−0.744	0.008	−0.857	0.007	...	...
3893.6	−0.715	0.006	−0.856	0.007	...	...
3903.8	−0.738	0.007	−0.830	0.005	...	...
3905.5	−0.727	0.005	−0.838	0.011	...	...
3907.2	−0.733	0.004	−0.835	0.008	...	...
3908.9	−0.727	0.004	−0.833	0.009	...	...
3910.6	−0.720	0.003	−0.827	0.004	...	...
3912.3	−0.730	0.007	−0.824	0.007	...	...
3914.0	−0.733	0.006	−0.809	0.014	...	...

Download table as:  ASCIITypeset images: 1 2 3

The position–velocity measurements from the three VLBI maps were combined by binning the velocities in 2.0 km s⁻¹ wide bins (comparable to maser line widths) and calculating variance-weighted positions and velocities. We associated these positions and velocities with maser feature accelerations, from the GBT time monitoring of spectra, by choosing the VLBI velocity closest to that of an acceleration fit. As there are fewer acceleration fits than VLBI position–velocity measurements, not all VLBI measurements have corresponding accelerations. For features lacking acceleration measurements we use only the position–velocity data when modeling. The velocity, position, and acceleration values (when available) are given in Table 1 and used to model the accretion disk.

Rather than using formal fitting uncertainties, which tend to be optimistic for high signal-to-noise data, we adopted more realistic "error floors" for the data uncertainties of ±0.01 mas for positions, ±1.0 and ±0.3 km s⁻¹ for the velocities of systemic and high-velocity maser features, respectively, and ±0.57 km s⁻¹ yr⁻¹ for accelerations. Error floors were added in quadrature to formal fitting uncertainties. For example, systematic errors of ≈0.01 mas between systemic and high-velocity features can be caused by an error in the absolute position of the reference maser spot of a few mas (Argon et al. 2007). Using error floors allows not only for systematic uncertainty not captured by formal estimates but also for slight incompleteness in the model of the accretion disk such as could come from unmodeled spiral structure. For example, there are indications of slight departures from perfect Keplerian gas orbits about a dominant central mass in the nearby, well-studied disk of NGC 4258 at levels of ≈0.5 km s⁻¹ yr⁻¹ for systemic-feature accelerations (Humphreys et al. 2008). In Section 4.1, we explore the sensitivity of the fitted parameters to changes in the magnitudes of these error floors.

Conceptually, were the masers in a perfectly thin and flat accretion disk orbiting circularly about a point mass and viewed edge-on, the position–velocity data for the high-velocity features would trace a Keplerian profile (i.e., $V=\sqrt{GM/r}$), where G is the gravitational constant, M is the mass of the central black hole, and r is the distance of a masing cloud from the black hole. The symmetry of the approaching and receding features can be used to precisely locate the central black hole both in position (x₀, y₀) and velocity (V₀). Thus, for high-velocity features that are in the plane of the sky, the VLBI image directly gives a feature's angular radius, θ = r/D, where D is the distance to the galaxy. Systemic-velocity masers in front of the black hole, moving transversely on the sky, display a change in velocity over time (A) as they orbit the central mass. Features in a thin annulus at radius r will be observed to accelerate at A = V²/r. Thus, from position–velocity and acceleration measurements, one can estimate  D = V²/Aθ and then H₀ ≈ V_r/D, where V_r is the recessional velocity. Note that proper motions of systemic-velocity masers can also be used to estimate distance, and Herrnstein et al. (1999) showed that the proper motion distance agreed with the acceleration distance estimate for the nearby galaxy NGC 4258. However, proper motions decrease in magnitude with source distance and are generally far less accurate than radial acceleration measurements based on changing Doppler shifts (which are distance independent).

In practice, accretion disks are somewhat more complicated than the idealized case just described. The disks are not precisely edge-on (although strong maser amplification prefers disk spin axes to be within about 5° of the plane of the sky (e.g., Watson & Wallin 1994), disks are often somewhat warped, and gas may be on slightly eccentric orbits). Also, galaxies typically have peculiar motions of hundreds of km s⁻¹ with respect to a pure Hubble flow (e.g., Masters et al. 2006). These complications are best addressed by constructing a model for maser orbits in the disk and adjusting the model parameters to best match the observations.

We modeled the system with up to 13 global parameters. The central black hole has a mass, M, is located at (x₀, y₀) on the sky (relative to the reference maser feature), and has a line-of-sight velocity, V₀, in the CMB rest frame. For UGC 3789, V₀ = V_LSR + 60 km s⁻¹. Following Herrnstein et al. (1999), we model the disk warp by a change in inclination with radius, i(r) = i₀ + (∂i/∂r)δr, and position angle (defined east of north) with radius, p(r) = p₀ + (∂p/∂r)δr, where δr = r − r_ref and r_ref is a reference radius assigned to the middle of the maser distribution. (Given the very small warping evident in the VLBI maps, we did not use second-order warping terms.) Our model can allow the masers to have eccentric orbits, with eccentricity e and pericenter rotated in angle ω = ω₀ + (∂ω/∂r)δr with respect to our line of sight. Finally, the galaxy (angular diameter) distance is calculated from H₀, assuming that its observed radial motion deviates from a pure Hubble flow by a peculiar velocity V_p relative to the CMB reference frame. The angular-diameter distance is calculated rigorously from formulae of Hogg (1999), assuming cosmological parameters for matter density Ω_m = 0.27 and dark energy $\Omega _\Lambda =0.73$.

We chose to estimate H₀ directly, rather than fit first for D, followed by a second step to estimate H₀. There are several advantages to our approach. First, unlike NGC 4258, which is not sufficiently distant to be in the Hubble flow and hence cannot be used to accurately estimate H₀, for UGC 3789 our primary interest is H₀ and not D. Second, when solving for D, systemic-feature accelerations, which are the least accurately measured type of data, appear in the denominator (D = V²/Aθ) and lead to an asymmetric posteriori probability density function (PDF) and complicate the error estimation for H₀. Third, we can directly incorporate uncertainty in the galaxy peculiar velocity (V_p) into the uncertainty in H₀ through the prior on V_p.

In addition to the global parameters, each maser feature requires two parameters, its radius, r, and azimuth, ϕ, to specify its location in the disk. Initial values for r and ϕ for high-velocity features were estimated from the position data, assuming ϕ = 90° for redshifted features and ϕ = −90° for blueshifted features. Systemic-velocity features were assigned initial r values based on the acceleration data and ϕ values based on the velocity data. All (r, ϕ) parameters were assigned flat priors and were adjusted in the fitting process.

We evaluated parameter posteriori PDFs with Markov chain Monte Carlo (McMC) trials that were accepted or rejected according to the Metropolis–Hastings algorithm. All parameters except for one were given flat priors and hence could vary freely, constrained only by the difference between data and model. The only parameter with a constraining prior was the galaxy's peculiar velocity, V_p; this comes from large-scale gravitational perturbations (Masters et al. 2006; Springob et al. 2007) which are expected to cause observed velocities for UGC 3789 to be lower by 151 ± 163 km s⁻¹ with respect to the Hubble flow (i.e., cz = V₀ + V_p). All velocities quoted here are non-relativistic, optical definition, in the CMB frame; the model and the data velocities (shifted from the observed LSR frame to the CMB frame by adding 60 km s⁻¹) were converted internally in the fitting program to relativistically correct values. Initial runs indicated very little warping (as is evident in the VLBI images) and maser orbital eccentricities near zero. Hence, for our "basic model," we used only the first-order warping parameters and assumed circular gas orbits.

We ran 10 "burn-in" stages, each with 10⁶ trials, in order to arrive at near-optimum parameter values and parameter step sizes. Parameter step sizes were iteratively adjusted after each burn-in stage so as to come from Gaussian distributions with similar widths as the anticipated posteriori PDFs (estimated from the burn-in McMC trials), multiplied by a global step-size factor. This factor (≈0.02) was also adjusted in the burn-in stages to scale parameter steps so that an optimal Metropolis–Hastings acceptance rate near 23% was achieved. The widths of the anticipated posteriori PDFs of the parameters and the global step-size factor, which together determine parameter step sizes, were taken from the last burn-in stage and then held constant for the final McMC trials.

After discarding the burn-in stage trials, we evaluated 10⁷ McMC trials to obtain final posteriori PDFs. The Pearson product–moment correlation coefficients from these trials are given in Table 2. The only correlation coefficient (r) with magnitude greater than 0.5 is between the Hubble constant and central black hole mass (r$_{{H_0} ,M}\approx -0.85$) (see discussion in Kuo et al. 2013). In order to more optimally sample the PDFs for these parameters, we modified the McMC values for these parameters to be totally anti-correlated for half of the trials and uncorrelated for the remainder (Gregory 2011).

From Bayes' theorem, the probability density of the parameters (ρ), given a model (m), data (d), and priors (I), is given by

$$\begin{equation*} P(\rho |m,d,I) \propto P(d|,\rho ,m,I) \times P(\rho |m,I). \end{equation*}$$

Since we assume Gaussianly distributed uncertainties for the data and priors, maximizing P(d|ρ, m, I) × P(ρ|m, I) is equivalent to minimizing $\chi ^2_d + \chi ^2_\rho$, where the subscripts d and ρ refer to the data and model parameters (given the priors), respectively. The best-fitting trial gave $\chi ^2_d = 1.50$ for 227 degrees of freedom. In Figure 3, we show projections of the data with the best-fitting model superposed. The modeling also yields the (r, ϕ) coordinates of each maser feature in the disk plane, which are displayed in Figure 4 (neglecting the slight warping of the disk).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Optimum values of the model parameters were estimated from the posteriori PDFs marginalized over all other parameters. They were generated from a total of 10⁷ McMC trials, obtained from 10 independent program runs starting with slightly different parameter values and new random number generator seeds. Each run produced 10⁷ trials, but stored only every 10th trial (i.e., "thinned" by a factor of 10). Parameter values given in Table 3 were produced from binned histograms for each parameter and finding the sample median and ±34% (≈ ± 1σ) range. Of some interest is the very accurate estimate of the mass of the central black hole of (1.16 ± 0.12) × 10⁷ M_☉.

The binned posteriori PDF for the Hubble constant, marginalized over all other parameters, is displayed in Figure 5, and the marginalized distributions for the other nine global parameters are shown in Figure 6. The H₀ distribution can be well approximated by a Gaussian with a mean of 68.9 km s⁻¹ Mpc⁻¹ and a ±34% confidence range of σ = ±5.8 km s⁻¹ Mpc⁻¹. Since, the reduced $\chi ^2_\nu = 1.50$, we conservatively inflate this uncertainty by a factor of $\sqrt{1.50}$ to ±7.1 km s⁻¹ Mpc⁻¹. This Hubble constant, coupled with the recessional velocity of cz = V₀ + V_p = 3466 km s⁻¹(optical definition relative to the CMB, and corrected for a peculiar velocity of 151 km s⁻¹), and assuming cosmological parameters for matter density Ω_m = 0.27 and dark energy $\Omega _\Lambda =0.73$, places UGC 3789 at an angular-diameter distance of 49.6 ± 5.1 Mpc. This corresponds to a luminosity distance of 50.8 ± 5.2 Mpc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Sensitivity of H₀ to Data Weighting

4.2. Secondary χ² Minima

4.3. Eccentric Gas Orbits

Direct measurements of H₀, such as from UGC 3789, are especially important for constraining the equation of state of dark energy, w, since the effects of dark energy are greatest at the present epoch. Our estimate of H₀ can be combined with results from the Wilkinson Microwave Anisotropy Probe (WMAP) mission to tighten constraints on the dark energy equation of state, w, independent of other methods such as using SN Ia or baryon acoustic oscillations (Eisenstein et al. 2005).

Figure 7 shows two-dimensional PDFs for w and H₀ with 95% and 68% confidence contours. The gray-scale contours were generated by binning the parameter values from Markov chains (wmap_wcdm_sz_lens_wmap7.2_chains_norm_v4p1.tar.gz) from the WMAP seven-year data (processed with WMAP version 4.1, RECFAST version 1.5, and modeled with a constant-w, ΛCDM model that incorporates the effects of the SZ effect and gravitational lensing). Fitting a Gaussian to the marginalized one-dimensional PDF for w yields w = −1.09 ± 0.37 (±68% confidence).

Zoom In Zoom Out Reset image size

Our new constraint that H₀ = 68.9 ± 7.1 km s⁻¹ Mpc⁻¹, when added to the WMAP PDF, is shown in Figure 7 with dashed red contours. This information for H₀, which is independent of w, improves the w-constraint by nearly a factor of two: w = −0.98 ± 0.20. This is probably about as accurate a result as one can obtain from measurements of the water maser emission from UGC 3789 with current equipment and a reasonable amount of observing time. Significantly improved measurements of this galaxy will likely await the completion of the Square Kilometre Array's high-frequency component.

Progress in the near future will come from measurements of other megamaser galaxies. The goal of the MCP is to determine H₀ with ±3% accuracy via ≈10 megamaser galaxy measurements. For example, were such measurements to yield H₀ = 68.9 ± 2.1 km s⁻¹ Mpc⁻¹, this would further tighten the constraint on the equation of state of dark energy (blue dotted contours in Figure 7) to w = −0.97 ± 0.10, independent of other methods.

We are grateful to the NRAO VLBA and GBT staff for their many contributions to the Megamaser Cosmology Project. We thank the referee for many valuable comments on the original version of this paper.

Facilities: VLBA - Very Long Baseline Array, GBT - Green Bank Telescope, Effelsberg - Effelsberg Radio Telescope