{'fraction_non_alphanumeric': 0.03838099998649032, 'fraction_numerical': 0.029221437159724944, 'mean_word_length': 4.129018847006652, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Refining the M_BH-V_c scaling relation with HI rotation curves of water megamaser galaxies

Black hole - galaxy scaling relations provide information about the coevolution of supermassive black holes and their host galaxies. We compare the black hole mass - circular velocity (MBH - Vc) relation with the black hole mass - bulge stellar velocity dispersion (MBH - sigma) relation, to see whether the scaling relations can passively emerge from a large number of mergers, or require a physical mechanism, such as feedback from an active nucleus. We present VLA H I observations of five galaxies, including three water megamaser galaxies, to measure the circular velocity. Using twenty-two galaxies with dynamical MBH measurements and Vc measurements extending to large radius, our best-fit MBH - Vc relation, log MBH = alpha + beta log(Vc /200 km s^-1), yields alpha = 7.43+/-0.13, beta = 3.68+1.23/-1.20, and intrinsic scatter epsilon_int = 0.51+0.11/-0.09. The intrinsic scatter may well be higher than 0.51, as we take great care to ascribe conservatively large observational errors. We find comparable scatter in the MBH - sigma relations, epsilon_int = 0.48+0.10/-0.08, while pure merging scenarios would likely result in a tighter scaling with the dark halo (as traced by Vc) than baryonic (sigma) properties. Instead, feedback from the active nucleus may act on bulge scales to tighten the MBH - sigma relation with respect to the MBH - Vc relation, as observed.


INTRODUCTION
The observed scaling relations between supermassive black hole (BH) mass and properties of the host galaxy, intensively studied over the past decade, suggest that black hole growth is tied to the growth of the surrounding host galaxy. These galaxy properties include the bulge/spheroid stellar velocity dispersion σ * (e.g., Ferrarese & Merritt 2000;Tremaine et al. 2002;Gültekin et al. 2009;Beifiori et al. 2012;McConnell & Ma 2013), the mass and luminosity of galaxy bulges (e.g., Marconi & Hunt 2003;Häring & Rix 2004;McConnell & Ma 2013), and the circular velocity V c (e.g., Ferrarese 2002;Kormendy & Bender 2011;Beifiori et al. 2012), which is the rotation velocity measured at large radius to probe the dark matter halos potential. It is intriguing that these power-law relations, especially the M BH − σ * relation, hold over several orders of magnitude in BH mass with small scatter, even though the black hole accounts for only a few thousandths of the mass of the galaxy (e.g., Häring & Rix 2004).
There are a wide array of theories attempting to explain the BH-galaxy scaling relations (e.g., Silk & Rees 1998;Ciotti & Ostriker 2001;Murray et al. 2005;Hopkins et al. 2006;Peng 2007). Two of the most popular models include variants of feedback from active galactic nuclei (AGN) and scenarios in which merging alone can lead to BH-galaxy scaling laws. In AGN feedback models, the central BH accretes mass and grows until it is massive enough to expel gas from the galaxy potential well and quench its own growth. BH growth in this picture is regulated by the depth of the galaxy potential well (Silk & Rees 1998;Fabian 1999;Di Matteo et al. 2005;Robertson et al. 2006). On the other hand, the pure merging scenario suggests that the correlation between linear quantities, for example the BH mass M BH and the halo mass M DM , can emerge from a large number of mergers based on the central limit theorem, even without a physical mechanism linking the two (Peng 2007;Hirschmann et al. 2010a;Jahnke & Macciò 2011).
Although both the feedback and merging phenomena may occur in galaxy evolution, it is unclear whether either of the mechanisms is essential in establishing the scaling relations. The most important physical scale for feedback is also a matter of debate (Booth & Schaye 2010a;Debuhr et al. 2010). Furthermore, we do not know how the AGN output couples to the gas, whether via thermal energy (Silk & Rees 1998) or momentum (Ostriker et al. 2010), nor do we know the average efficiency of the feedback.
Therefore, empirical evidence that distinguishes the relative importance of different physical processes in establishing the scaling relations is key to constructing the coevolution history of black holes and galaxies. In this paper, we investigate the origin of BH-galaxy scaling relations by comparing the M BH − V c relation to the M BH − σ * relation. Circular velocity V c is a good indicator of dark matter halo mass, and velocity dispersion σ * serves as its counterpart on bulge scales. While some AGN feedback scenarios (e.g., Debuhr et al. 2010) suggest that BH mass will be most tightly linked to baryons (as opposed to dark matter), a pure merging scenario suggests that the M BH − M DM (or M BH − V c ) relation should be the cleanest and tightest relation, as it is free from the baryonic physics (e.g., star formation) that occurs during merging. Comparison of the two relations, especially their scatter, can help determine the mechanism that drives BH-galaxy coevolution (Ferrarese 2002;Novak et al. 2006;Kormendy & Bender 2011). Ferrarese (2002) first proposed that BH mass may correlate with dark matter halo mass, based on the M BH −σ * relation and the correlation between σ * and V c . Later, a number of papers (Pizzella et al. 2005;Courteau et al. 2007;Ho 2007) pointed out that the σ * − V c relation depends on surface brightness, light concentration, and morphology, and suggested that the M BH − σ * relation, not the M BH − V c relation, is most fundamental. Kormendy & Bender (2011) compiled a sample of 25 galaxies with both dynamical BH mass measurements and V c from spatially resolved rotation curves. From this direct M BH − V c correlation they concluded that the dark matter halo mass alone cannot determine the BH mass, given that the BH mass can range from < 10 3 − 10 6 M at a circular velocity of 120 km s −1 (for a different view see Volonteri et al. 2011). Beifiori et al. (2012) also found the scatter in the M BH − V c relation to be about twice as large as that in the M BH − σ * relation using a large sample of M BH upper limits from Hubble Space Telescope spectra (Beifiori et al. 2009) and V c from unresolved H I line-width measurements. For a comprehensive review, see Kormendy & Ho (2013).
In this work, we aim to refine the M BH − V c statistics. We start with the most up-to-date galaxy sample with dynamically measured black hole masses, and study all of these that also have spatially resolved circular velocity measurements, predominantly disk galaxies observed with H I. We present five new H I rotation curves measured with the Karl G. Jansky Very Large Array 1 (VLA) during the Expanded Very Large Array (EVLA) construction period. Three of these galaxies have dynamical M BH measured from the kinematics of water megamasers in sub-pc disks (Reid et al. 2009;Greene et al. 2010;Kuo et al. 2011). With BH mass errors smaller than 11 percent, these megamaser measurements are especially useful in constraining the intrinsic scatter of the M BH − V c relations. In total our sample contains thirtythree galaxies. We assign V c to all galaxies in a consistent way and investigate V c reliability as a function of the spatial extension of the rotation curve. Using only reliable V c measurements at large radius, we quantify the M BH − V c correlation and compare it with the M BH − σ * relation. We investigate whether M BH is correlated more tightly with V c or σ * , to discriminate AGN feedback scenarios from those in which merging alone establishes the scaling relations.

H I OBSERVATIONS
We observed 5 spiral galaxies in H I with the VLA.
The observations were taken in the L-band and the C configuration under project 10B−220 between October and December 2010 (for details see Table 1). The observations used dual circular polarizations and a single spectral window with 256 channels across the 4MHz (852 km s −1 ) bandwidth, giving a channel width of 3.3 km s −1 . Twenty-one edge channels on each side were flagged due to a higher noise level. Removing these channels did not affect our results since none of the sources had emission in these parts of the band. The full width at half power of the primary beam is 32 , and the synthesized beam size ranges from 23 to 52 depending on the source declination and number of antennas used ( Table  2).
The duration of each observation was 5 hours and the on-source time was ∼ 200 minutes per source. During a track, the gain calibrator was observed for 3 − 6 minutes every 20 − 27 minutes, while the flux calibrator was observed once, for 15 − 23 minutes. The observations were carried out during the EVLA upgrade phase and some of the L-band receivers were not yet installed. For the observations of NGC 1194, NGC 2748, and UGC 3789, four antennas were not in the array, three were missing during the NGC 7582 observations, and two during the NGC 2960 track. One additional antenna was flagged out of the data set for NGC 1194, NGC 2748, NGC 7582, and UGC 3789 for high noise levels. The total number of antennas used was between 22 and 25, as listed in Table 1. Radio Frequency Interference (RFI) was found in three observations, and the contaminated data were excluded, which accounts for only a few percent of the data.
Data calibration and image processing were carried out with the Common Astronomical Software Applications (CASA) package. Time averaging of 10 seconds and online flagging were applied before calibration. The amplitude loss due to time averaging is less than 1% according to the VLA Observational Status Summary. The antenna position correction, antenna-based delays, bandpass, phase gain, amplitude gain and flux scale calibrations were carried out consecutively. The minimum acceptable signal-to-noise ratio (SNR) for the bandpass and gain calibration tables were all set to 3. The bandpass solutions were solved for each channel, which gave stable solutions with amplitude variations 1%. The phase gain and amplitude gain were solved for each calibrator scan. Gain elevation curves were applied to correct for elevation-dependent amplitude gain due to antenna deformation. Data weighting was not used throughout the calibration process, since the VLA does not monitor system temperature. The two polarizations were calibrated separately and combined in the imaging process.
For the H I emission line analysis, the continuum was estimated from line-free channels with more than 20 channels on each side and subtracted in the uv plane with the task uvcontsub. A constant continuum level was used for all galaxies but NGC 1194 and NGC 2960. In these two galaxies, an adapted linear spectrum was used in order to subtract a nearby bright source with an inclined spectrum. The nearby bright continuum sources were well-subtracted at this stage, and the contamination from their side-lobes was negligible.
For imaging we used the CASA task clean in velocity mode. This mode corrects for the doppler shift due to rotation of the Earth. In this process the channels were also regridded linearly into wider spectral bins -the image planes. The width of the image planes were chosen to be 10 km s −1 for higher signal-to-noise (SNR) observations (NGC 2748 and NGC 7582), and 20 km s −1 for lower SNR observations (NGC 1194, NGC 2960, andUGC 3789). The SNR ratios for each data cube, defined as the peak intensity per channel divided by the RMS noise, are listed in Table 2. For NGC 2748 and NGC 7582 we applied Briggs weighting with the parameter robust set to 0.5, and natural weighting was applied for NGC 1194, NGC 2960, and UGC 3789. All images were cleaned to the 3 σ level. All velocities in this paper use the optical convention and the Local Standard of Rest Kinematic (LSRK) reference frame. Position angles (P.A.) in this paper are measured east of north.

ANALYSIS
In this section we describe the procedures that we used to extract velocity fields from the H I data cubes and measure rotation curves and V c . The H I properties of individual galaxies are discussed in Appendix A.
3.1. Non-parametric moment-0 and moment-1 maps Non-parametric integrated intensity moment-0 maps of the five galaxies are shown in the left columns of Figures 1 to 5. In order to improve the SNR, we used the CASA tool image.moments, which masked pixels with no signal in producing the moment maps. The masks were produced as follows: the data cube was first smoothed spatially over one beam and spectrally over two image planes (20 km s −1 for NGC 2748 and NGC 7582; 40 km s −1 for NGC 1194, NGC 2960, and UGC 3789). Then, pixels with smoothed intensity below a certain threshold were masked. The threshold was set to 3 σ for NGC 2748, NGC 7582, and NGC 1194, and 2.5 σ for NGC 2960 and UGC 3789, chosen to optimize the signal seen in the moment-0 maps.

Parametric Fitted Velocity Fields
To extract rotation curves from the H I data cube, it is conventional to fit a velocity-field model. The twodimensional velocity-field model is parametrized by geometrical parameters and the rotation curve. Then, the model is fit to the two-dimensional velocity-field data derived from the 3-D data cube. To construct the velocity field, one assigns a representative velocity to each spectrum in the data cube. In the ideal case of a simple rotating disk with small random motions, the projected rotational velocity will correspond to the peak velocity in the spectrum. In practice there are different ways to measure the velocity at each position, including the peak velocity, the intensity-weighted mean velocity, and parametrized fitting (e.g., Gaussian fitting, deBlok et al. 2008). We seek the velocity measure that can best recover the true projected rotation velocity in the case of finite resolution, sensitivity, disk warping and overlapping velocities along the line of sight.
Measuring the peak velocity is the most straightforward approach, but the measurement is sensitive to the noise. It may pick up a noise spike instead of the real signal when the SNR is low. Moreover, it is discretized to the velocity channel width, which is 10 or 20 km s −1 in our case. The intensity-weighted mean (moment-1) velocity is less sensitive to the noise, but will be biased if the spectrum is not symmetric. For example, beam smearing will produce a wing towards the systemic velocity, causing the intensity-weighted mean velocity to be skewed towards lower values (see the left panel of Fig.  6).
Another alternative is parametrized fitting, such as with a Gaussian line profile. The Gaussian fit picks up the intensity peak if the line profile is symmetric, and it is not as sensitive to the noise as is the peak velocity, so Gaussian fitting is our preferred approach. However, at high SNR if the line shape is asymmetric, then the mean velocity of the best-fit Gaussian profile will be biased toward the wing (see the right panel of Figure 6). In this case, we can use a Gauss-Hermite expansion to capture the non-Guassian line shape (van der Marel & Franx 1993). For example, the skewness asymmetry is taken into account using the third order Hermite polynomial H 3 (x) Here the skewness parameter h 3 quantifies the asymmetry of the line shape. The intensity-weighted mean velocity v m,3 of this line shape f (v) is calculated by van der Marel & Franx (1993) to be v m,3 =v + √ 3σh 3 to first order. Compared to a Gaussian, the mean velocity v m,3 measured by Gauss-Hermite fitting is less biased towards the wing and closer to the peak, while compared to the peak velocity the Gauss-Hermite fit is also less sensitive to the noise. Differences in the velocity fields derived from these different methods and their effects on the resulting rotation curves are discussed in Appendix B. We find that for high SNR data, the rotation curve derived from the moment-1 velocity field is systematically lower than from the Gaussian/Gauss-Hermite velocity field at the ∼15% level, and that the peak-velocity produces noisy rotation curves for low SNR data.
Considering the advantages of Gaussian and Gauss-Hermite fitting, we utilized these two methods to construct the velocity field. In order to minimize contamination from noise, we applied the noise masking described in §3.1. We then performed a Markov Chain Monte Carlo (MCMC) fitting of a Gaussian line shape to the five H I data cubes. We assume the likelihood function of intensity is Gaussian, and the sigma of this Gaussian is a constant throughout the data cube. The sigma was measured from the RMS intensity of line-free channels in the cube.
For the three highest-SNR galaxies, i.e., NGC 2748, NGC 7582, andNGC 1194 (SNR > 14), there were regions along the major axis with high reduced-χ 2 values, i.e., where the cumulative distribution function of χ 2 is larger than 95 percent, indicating deviations of the line shape from a Gaussian, mostly because of the line asymmetry described above. To better capture the skewness of the line, we fitted these high χ 2 spectra with a thirdorder Gauss-Hermite function. The final velocity fields of NGC 2748, NGC 7582, and NGC 1194 were constructed using a mixture of Gaussian and Gauss-Hermite fitting. The lower SNR cubes of NGC 2960 and UGC 3789 (SNR ∼ 6) do not have enough SNR to show deviations from single Gaussians.
Higher-order deviations from the Gauss-Hermite line shape, such as the fourth-order kurtosis (h 4 ), still exist in smaller particular regions after using the third-order polynomial; along the major axis of the edge-on galaxy NGC 1194 the line-of-sight overlap is severe, and along the spiral arms in NGC 7582 the gas kinematics are complex. We did not apply higher-order fits as the inclusion of the H 4 polynomial would not change the mean velocity to the first order (van der Marel & Franx 1993).
One advantage of the MCMC fitting is that the uncertainty in each fitted parameter can be estimated from the standard deviation of the probability distribution probed by the chain. The uncertainties in the velocities are especially useful. A large velocity uncertainty may indicate that the spectrum contains no emission line, only noise. We therefore performed a further masking to exclude pixels with velocity uncertainty larger than 80 km s −1 . The final velocity fields are shown in the right columns of Figures 1 to 5.

Rotation Curves
Rotation curves are derived by fitting the velocity fields ( §3.2) with a rotating disk model. This model consists of coplanar concentric rings, each with its own rotation velocity. The width of each ring is chosen to be 14 , corresponding to the typical beam size. We do not use tilted-ring modeling, but we show in Appendix B that the rotation curves derived with a tilted ring analysis are consistent with those presented here.
The rotation-curve fitting was performed with the same MCMC method. We assume that the likelihood function of velocity is Gaussian, and use the velocity uncertainty estimated in §3.2, such that noisier pixels with larger velocity uncertainties are down-weighted naturally. We fit five geometrical parameters, the center x 0 , y 0 , the systemic velocity v sys , the inclination i, and the position angle P.A., together with the rotation velocity, at each ring. These parameters describe the velocityfield model. For the megamaser galaxies (NGC 1194, NGC 2960, andUGC 3789) we fix the central positions (x 0 , y 0 ) to the VLBI maser positions, and we fit the kinematic center of NGC 2748 and NGC 7582 as free parameters. The v sys , i, P.A., and rotation velocities are fitted as free parameters with the exceptions of the inclination in the case of NGC 2960 and UGC 3789, and the P.A. of UGC 3789, where these values cannot be constrained by the data. Therefore, we fix these three values to the HyperLeda values (Paturel et al. 2003), derived from the optical images. The best-fit geometrical parameters are listed in Table 3.
The best-fit rotation curves are plotted in Figure 7. Various error sources in the rotation curves are considered, including the fitting errors estimated by the standard deviation in the Markov chains, the RMS errors estimated by the RMS variation of the residuals in a ring, and the errors induced by uncertainties in inclination and P.A. It is commonly found that the formal fitting errors in the rotation curve fits are unrealistically small (e.g., deBlok et al. 2008), likely because the rotating disk model fails to capture features in the observed velocity field such that the residuals are non-Gaussian, causing the fitting errors to appear artificially small (5-20 times smaller than the RMS error). As shown in Appendix B, the RMS error is larger than or comparable to the differences between various velocity assignment methods and different parameterized rotating disk models, and therefore provides a conservative estimate of the potential systematic uncertainties. For NGC 2960 and UGC 3789, we also include the uncertainties in the literature inclination and/or P.A. values. The uncertainties contributed by these values are estimated by recalculating the rotation velocities assuming an inclination and/or P.A. within ±5 • of the HyperLeda value. The final error is the quadratic addition of the considered error sources.

The Sample
Our primary sample is an updated list of galaxies with both dynamical M BH and spatially resolved V c measurements, as tabulated in Table 4 and 5. There have been a number of compilations of M BH and V c (Kormendy & Bender 2011;Beifiori et al. 2012); we introduce a few significant improvements over these. First, we update five maser and three stellar/gas dynamical M BH measurements. In addition to our five VLA H I rotation curve V c measurements, there are a few other V c measurements included in this sample from the literature that are not compiled in previous M BH − V c studies. We also reexamine the literature V c values with a careful reliability analysis, to ensure that our measurements and the literature values are derived in a consistent manner. Thus, our V c values may differ from previous work. In addition to this primary sample, we also consider a secondary sample composed of galaxies with three spatially unresolved V c (single dish H I measurements) and fifteen black hole mass upper limits. This secondary sample is described in Appendix C and plotted in Fig. 8 for comparison, but is not used in our scaling relation fitting analysis. Table 4 lists the V c values, sources, and rotation-curve properties of our primary sample. We examined each rotation curve from the original literature to assign a V c value. Because V c is used as a tracer of the dark matter halo potential in the outer parts of the galaxy, we assign V c as the rotation velocity at the largest observed radius R o , which is consistent with the definition in Ferrarese (2002). If the inclination correction is applied in the original literature, we use this inclination-corrected rotational velocity, otherwise we apply a simple 1/ sin(i) inclination correction using the optically determined inclination from HyperLeda (Paturel et al. 2003). The uncertainty in V c that results from our inclination correction treatment is estimated to be less than ±10% for a typical inclination error of ±5 • in the optically derived inclination (e.g., Dressler & Sandage 1983) at an average inclination of ∼ 60 • . This is comparable to or smaller than the V c measurement errors described below.

Circular Velocities Vc
The observational error in V c , listed first, is taken from the original literature. If it is not presented then we assume a 10% observational error, which is typical among our sample. The second error stands for the uncertainties due to rotation curve variation. In the ideal case where the rotation curve is flat at large radius, there is one uniquely defined V c . However, if the rotation curve is not flat but keeps rising or starts to fall, then V c depends on where the outermost observed radius R o is located. To take this uncertainty into account, we assign a second error equivalent to the amplitude of variation in the rotation curve. While these error assignments may overestimate the real error between the observed V c and the true halo potential, we hope to avoid overestimating the intrinsic scatter int in the M BH − V c relation. The final error used in the M BH − V c relation fit is taken as the larger of the observational and rotation curve variation error, as listed in Table 5.
There are two galaxies, NGC 3998 and NGC 5128 (Centaurus A) in our primary sample that have no rotation curve measurements but have spatially resolved H I data. For NGC 3998 we adopt the V c value and observational error estimated by the original literature . For NGC 5128, we determine V c using the published P-V diagram, and apply an inclination correction using the inclination from HyperLeda (Paturel et al. 2003). To reflect our lack of knowledge about the rotation curve trends, we assign conservative rotation curve variation errors of 20% in these two cases.
How faithfully V c reflects the potential of the dark matter halo depends on whether the rotation curve is probing the halo-dominated part of the galaxy. Therefore we compare the outermost observed radius of the rotation curve R o with the optically determined R 25 , which is the B = 25 mag arcsec −2 isophote from RC2 (de Vaucouleurs et al. 1995) shown in Column (6) of Table 4. It is worth noting that if R o /R 25 is small, as in the case of many optical rotation curves, then the rotation curve may be dominated by baryons rather than the dark matter halo, in turn biasing our measurement of intrinsic scatter in the M BH − V c relation. Thus, we divide galaxies into four groups according to the spatial extent and shape of their rotation curves. The first group has the largest spatial extent (R o > R 25 ), followed by the second group (R 25 > R o > R 25 /2). The third and fourth groups both have short rotation curves (R o < R 25 /2), but in the third group the rotation curves flatten while those ones in the fourth group are still rising until the last bin. We discuss the reliability of the V c measurements for these groups in §4.3 and will conclude that only the first two groups (R o > R 25 /2) have reliable V c for our M BH − V c analysis.  Table 5 lists V c , M BH , σ * and other galaxy quantities. The black hole mass M BH , stellar velocity dispersion σ * , morphology, and distance are adopted from McConnell & Ma (2013), except for NGC 4526 (Davis et al. 2013). The method of M BH measurement is also listed. We note potential caveats in the M BH measurements of three galaxies as pointed out by McConnell & Ma (2013). Lodato & Bertin (2003) measured the black hole mass in NGC 1068 and found a non-Keplerian maser disk. For NGC 2748, Atkinson et al. (2005) cautions that the determination of the disk center may be affected by heavy extinction. Finally, the black hole mass in NGC 7457 as measured by Gebhardt et al. (2003) could be overestimated if the central bright source excluded in their dynamical modeling is a nuclear cluster instead of an AGN. While we keep these three galaxies in our sample, we note these points of caution. The σ * that we adopt from McConnell & Ma (2013) is the light-weighted stellar velocity dispersion within one effective radius. However, for some galaxies the σ * is lower than in previous studies (e.g., Gültekin et al. 2009;Beifiori et al. 2012), as the BH sphere of influence is excluded in the velocity integration to avoid contamination. NGC 3998 and NGC 4594 in our primary sample, as well as NGC 1399 and NGC 4486 in our secondary sample, have σ * updated by McConnell & Ma (2013) with this correction. We also note that there is a new distance update from megamaser measurements for UGC 3789 (Reid et al. 2013), but we adopt the distance in McConnell & Ma (2013), which is consistent with Reid et al. (2013).

Fitting Method
To quantitatively analyze the M BH − V c and M BH − σ * relations, we fit each scaling relation with a power-law functional form Here α is the intercept and β is the logarithmic slope of the relation. We use √ 2σ * to compare with V c , motivated by the widely used singular isothermal sphere model where V c = √ 2σ * (Binney & Tremaine 2008). Observationally Ferrarese (2002) also found an almost linear relation between V c and σ with a ratio close to √ 2. We use a maximum-likelihood fitting method similar to Gültekin et al. (2009) that takes intrinsic (or cosmic) scatter int into account, with modifications regarding the V c error treatment. For simplicity we assume that the probability distribution of both the observational errors and intrinsic scatter are normal in log space, and all the observational errors are symmetrized in log space before fitting. Here we denote and the normalized Gaussian error distribution to be G (x) with mean zero and standard deviation .
To take the error in the independent variable (V c and σ * in this case) into account, we have an error term σ ν in the likelihood function. We write the likelihood for observing one galaxy (µ i , ν i ) given an intrinsic scaling relation µ = α + βν with intrinsic scatter int to be This treatment of the errors in the independent variable is less computationally expensive than the Monte Carlo method used in Gültekin et al. (2009). We then maximize the total likelihood to get the best-fit scaling relation parameters, the intercept α, the slope β, and the intrinsic scatter int . This maximum-likelihood method is similar to the minimum χ-square method described in Tremaine et al. (2002), except for the 1/ √ 2π normalization factor of the Gaussian function. We adopt the ∆χ 2 = 1 confidence limit as an error estimate for the fitted parameters α, β, and int .

Fitting Results
We find that M BH is correlated with both V c and σ * , with Spearman rank correlations of ρ = 0.62 (p-value 2 × 10 −3 ) and ρ = 0.67 (6 × 10 −4 ) respectively, for the reliable sample of 22 R o > R 25 /2 galaxies discussed below, meaning that a correlation exists for both of the relations. Using this sample, the best-fit M BH −V c relation is α = 7.43 +0.13 −0.13 , β = 3.68 +1.23 −1.20 , and int = 0.51 +0.11 −0.09 , and is plotted as blue lines in Figure 8. There is no significant correlation between the residuals in the M BH − V c relation and σ * (Spearman ρ = 0.33, p-value= 0.13). The parameters of different M BH − V c samples as well as the comparison with the M BH − σ * relation are shown in Table 6.
To determine a reliable sample and investigate how the spatial extension of the rotation curves affect the reliability of the V c measurements, we derive the scaling relations using sample with different R o criteria. The first group of fourteen galaxies (R o > R 25 ) should be the most reliable sample with rotation curves that probe the outer, halo-dominated, region of the galaxy. The second sample of eight galaxies (R 25 > R o > R 25 /2) also gives the same intrinsic scatter as the first group, showing that the V c measurements with R 25 > R o > R 25 /2 are as reliable as the first group for our purpose. When we include groups three and four, the eleven shorter rotation curves with R o < R 25 /2, which are predominantly measured in the optical, the intrinsic scatter increases significantly. It is a sign that for these short rotation curves the halo potential is not represented by the observed V c , and the intrinsic scatter will be inflated artificially if they are included in the analysis. We therefore rely only on the R o > R 25 /2 sample of twenty-two galaxies from the first two groups for our main scientific discussion.
Using this R o > R 25 /2 sample, the M BH − V c relation has an intercept consistent with the M BH − σ * relation, if V c corresponds to √ 2σ * . While the slopes of the M BH − V c and M BH − σ * relations are consistent within our sample, our M BH −V c relation slope falls on the lower end of reported slopes from studies of the M BH − σ * relation using larger samples, (e.g., Gültekin et al. 2009;Mc-Connell & Ma 2013), possibly due to our limited dynamic range. The intrinsic scatter in the M BH − V c relation of int = 0.51 +0.11 −0.09 is similar to that in the M BH − σ * relation of int = 0.48 +0.10 −0.08 when using matched R o > R 25 /2 samples, but is larger than the M BH − σ * relation scatter int = 0.38 for the entire sample in McConnell & Ma (2013), which includes both early and late-type galaxies. The implications of our fits are discussed in §6.
Some elliptical galaxies also have V c measured from dynamical modeling. Although these measurements may involve different systematics than the rotation-curvederived V c for disk galaxies, we plot six such ellipticals from Beifiori et al. (2012) in Fig. 8 as an attempt to probe the M BH − V c trend at the high mass end. NGC 1399, NGC 3379, NGC 4374, NGC 4472, and NGC 4486 (M87) have V c measured by Kronawitter et al. (2000) using stellar kinematics, while NGC 3608 is measured by Coccato et al. (2009) (2013) found a similar steepening in the M BH −σ * relation when considering both early and late-type galaxies. On the other hand, the intrinsic scatter decreases slightly by about 0.5 σ to int = 0.46 +0.10 −0.09 and int = 0.45 +0.08 −0.07 for the M BH − V c and M BH − σ * relations respectively. These results will not be used in our scientific discussion, as their interpretation awaits better understanding of the correspondence between the V c from the dynamical modeling and from H I rotation curves.  (2011), which have dynamical M BH measurements and V c from spatially resolved H I or stellar kinematics, using the methods of Gültekin et al. (2009). They constrain the M BH − V c relation to have α = 7.39 ± 0.14, β = 4.33 ± 0.93, and int = 0.53 ± 0.10, consistent with our results. Beifiori et al. (2012) study the M BH − V c relation using 28 galaxies with M BH compiled by Gültekin et al. (2009), and V c from resolved H I kinematics, unresolved H I line width, or dynamical models of stellar kinematics in early type galaxies. Their M BH − V c relation parameters are α = 7.82 ± 0.15, β = 3.29 ± 0.61, and intrinsic scatter int = 0.51 ± 0.09. Their intercept α is different from our result, possibly because more elliptical galaxies are included in their sample. Their fit for the slope has a larger dynamic range because they include more elliptical galaxies, but is still consistent with ours, as is their intrinsic scatter measurement.

POSSIBLE SYSTEMIC UNCERTAINTIES
Intrinsic scatter in the scaling relations provides an important discriminant between different BH scaling relation scenarios. Various observational errors in V c have been taken into account in §4 to avoid contaminating the intrinsic scatter measurement, including the inclination correction, the observational error, rotation curve variations, and the uncertainties due to short rotation curves. Here we discuss whether other systematic uncertainties, such as sample bias, may affect the interpretation of, and the comparison between, the M BH − V c and M BH − σ * relations. Also, in order to interpret the M BH − M DM relation from the observed M BH − V c relation, we discuss the scatter that may be introduced in the V c − M DM conversion. We conclude first that the comparison between the M BH − V c and the M BH − σ * relations should be fair even in the face of selection effects, and second, that the scatter introduced by the V c − M DM conversion or halo triaxiality is not important compared to other sources of uncertainty.

Selection Effects
If our sample is subject to selection effects, it may not represent the true distribution of objects. Selection effects may exist for both BH and circular velocity measurements. First, it has been argued that analysis excluding M BH upper limits might be biased toward more mas-sive BHs at a given velocity dispersion, especially for latetype galaxies, since only massive BHs have a sphere of influence large enough to be spatially resolved . Beifiori et al. (2012) includes a large number of M BH upper limits in their fits to the M BH − V c and M BH − σ * relations, which should be closer to the true distribution. They find a slope β for the M BH − V c relation consistent with ours, but cannot derive a reliable intrinsic scatter. While concerns over M BH selection biases cannot be excluded definitively, we use the same sample for both M BH − σ * and M BH − V c relations, such that if M BH biases exist, they have the same effect on the two relations.
In addition to selecting against low-mass BHs, we also select against galaxies that have no V c from H I rotation curves, which includes massive elliptical (gas-poor) galaxies, and disk galaxies that either have only optical (mostly covering a limited radial distance) or no rotationcurve measurements. With no stellar dynamical V c measurements for ellipticals, our dynamic range is limited compared to the M BH −σ * relation in other studies, (e.g., Gültekin et al. 2009;McConnell & Ma 2013). While selection biases in V c measurements also exist, at least there is no obvious bias on the intrinsic scatter inherent with these selections. To estimate the magnitude of this scatter, we take the Navarro-Frenk-White (NFW; Navarro et al. 1996) model as an example. This model is a mass profile that is parametrized by the inner density ρ s and the characteristic inner radius r s , or equivalently the virial mass M vir and the concentration c vir ≡ R vir /r s . The virial radius R vir is defined as the radius within which the averaged halo density exceeds the background density by a certain factor, the virial mass M vir is the mass enclosed within R vir , and the virial velocity V vir is the rotation velocity at R vir . Because there is a one-to-one correspondence between M vir and V vir , the main uncertainty actually comes from the conversion between V vir at the virial radius (typically a few hundred kpc) and V c at the observed radius (typically 10-20 kpc, Ferrarese 2002). The observed V c is close to the maximum velocity V max occurring at radius r max ∼ 2.16 r s , which is a few to a few tens of kpc. The ratio V max /V vir depends on c vir .
The scatter in concentration is estimated by Bullock et al. (2001) and corrected by Wechsler et al. (2002) to be ∆ log 10 c vir = 0.14 dex, which translates to an error of 0.035 dex in V max or 0.14 dex in M BH for a slope of β = 4 in the M BH − V c relation. There is also a weak dependence of concentration on the virial mass, which affects the slope of the M BH − V c relation in principle. Additionally, if the halo is not spherically symmetric, as assumed by the NFW profile, but is prolate or triaxial, there will be extra uncertainty in the V c − M DM conversion. Franx & de Zeeuw (1992) estimate that the low observed scatter in the Tully-Fisher relation constrains the ellipticity of the halo potential in the disk plane to be ≤ 0.1. Triaxiality can at most contribute a scatter of 0.026 dex in V c or 0.1 dex in M BH in the M BH − V c relation.
Removing these additional sources of uncertainty in quadrature, assuming that they are uncorrelated, changes the final intrinsic scatter very little (0.48 rather than 0.51). Thus we conclude that those sources of scatter are small compare to outstanding measurement uncertainty.

DISCUSSION
We will now consider the implications of our new fits for the evolution of BHs and galaxies. In §6.1, we hope to gain new insight into whether AGN feedback is required to explain the BH-bulge scaling relations (e.g., Silk & Rees 1998;Hopkins et al. 2006), or whether galaxy merging alone may lead to tight BH-bulge correlations in massive galaxies (e.g., Peng 2007; Jahnke & Macciò 2011). In §6.2 we discuss the implications of M BH upper limits as outliers in the M BH − V c relation.

Is It Feedback or Merging That Sets the Scaling
Relations? We address the relative importance of feedback and merging empirically by asking which can account for features in the observed M BH − V c and M BH − σ * relations. We start by reviewing the expected slope, outlier behavior, and intrinsic scatter of the scaling relations in the context of each scenario. We will find that while more theoretical guidance and a larger sample size are needed to make inferences from the outliers and slope measurements, the intrinsic scatter is already an available and useful discriminant between the two scenarios.
In the galaxy merging scenario (Peng 2007;Hirschmann et al. 2010b;Jahnke & Macciò 2011), the correlations between M BH and galactic quantities (e.g., M * , M bulge , M DM ) emerge from a large number of mergers simply by the central limit theorem, even without a physical mechanism linking the two. In the context of this scenario, we expect the following features of the scaling relations. First, as the number of mergers increases and the galaxies grow larger, the scatter between M BH and galaxy properties will decrease and the distribution will converge to a tighter correlation. Therefore, it is expected that the intrinsic scatter int will decrease as the number of mergers increases towards larger V c and σ * . Outliers, with BH masses deviating from a power-law scaling relation, become increasingly unlikely towards higher masses. Second, the logarithmic slope β of the correlation between two linear quantities, e.  (Ferrarese 2002), we should find the slope of the M BH − V c relation, β, to be close to 3 − 3.32, somewhat smaller than, but allowed by, our observations. Lastly, the M BH − M DM relation should be tighter than other scaling relations, e.g., the M BH − M * or M BH − M bulge relations, because M DM simply adds during merging and does not depend on the baryonic physics that M * and M bulge are subject to. For example, during merging, extra factors of star formation and disk to bulge conversion have to be considered for the evolution of M bulge , and this will make the scatter in the M BH − M bulge relation larger than that in the M BH − M DM relation, if the baryonic regulation mechanism is not present. However, we note that diffuse DM accretion in principle may also enhance the scatter in M DM .
Feedback, on the other hand, provides a physical mechanism that couples BH mass directly with the galaxy potential well (e.g., Silk & Rees 1998). If the BH is massive enough, its energy or momentum output can expel gas from its vicinity and quench the further growth of both the BH and the galaxy. This feedback loop sets an upper limit on the BH mass for a given potential well depth, measured by the stellar velocity dispersion σ * of the bulge or V c of the dark matter halo. However, it is unclear whether this self-regulation process takes place on the scale of the BH sphere of influence, the bulge (Debuhr et al. 2010), or the dark matter halo (Di Matteo et al. 2003;Booth & Schaye 2010b). Some predictions of the feedback scenario are as follows. First, the scaling relations form an upper boundary for M BH . The BH cannot grow above the relation via accretion (although over-massive BHs could already be in place from massive seeds Volonteri & Natarajan (2009)). Second, the slope of the M BH − σ * or the M BH − V c relation is predicted to be five for energy feedback (Silk & Rees 1998) and four for momentum feedback (Murray et al. 2005). Third, the correlation between BH mass and the potential indicator (σ * or V c ) on the feedback scale should be tighter than on other scales.
We here compare the statistics of our V c sample with the predictions discussed above regarding the outlier properties, slope, and the intrinsic scatter, respectively. There are no obvious outliers above the M BH − V c relation in our sample, possibly due to limited sample size. However, there are bulgeless galaxies with M BH upperlimits that are low compared to their V c , which will be further discussed in §6.2. Regarding the measured slopes, with such limited dynamic range our data are still inadequate to discriminate the two scenarios. Both scenarios are consistent with the wide range of allowed M BH − V c slopes (β = 2.48 -4.91).
However, the comparison of intrinsic scatter between the two relations can provide useful constraints on the origin of the scaling relations. Our measurements show that the scatter in the M BH − V c relation is similar to that in the M BH − σ * relation, within the errors. It is worth noticing that the intrinsic scatter in M BH − V c , int = 0.51 +0.11 −0.09 dex, is a conservatively low estimate, as we have assigned large errors to the V c measurements to account for various observational uncertainties, including inclination uncertainties and rotation curve variations. Even so, we still find a value that is consistent with our measured M BH − σ * scatter, int = 0.48 +0.10 −0.08 dex, and that is not smaller than the M BH − σ * intrinsic scatter int = 0.46 for late-type galaxies constrained by McConnell & Ma (2013) or int = 0.38 for the entire sample in McConnell & Ma (2013).
We find no evidence that the M BH − V c (or M BH − M DM ) relation is significantly tighter than the M BH − σ * relation, in contradiction with naive expectations from the pure merging scenario. Therefore, our result disfavors merging as the only mechanism that ties together BHs and galaxies. Note that we use σ * to trace the baryonic mass in the bulge, while Jahnke & Macciò (2011) looked at the M BH − M bulge relation instead. Unfortunately, observations of M bulge for lower-mass galaxies with dynamical M BH are scarce, so a direct comparison between M BH − M DM and M BH − M bulge is not yet possible. If the real scatter in the M BH − M bulge relation is even smaller than in M BH −σ * , our conclusion is stronger. Even if the scatter in the M BH −M bulge relation turns out to be larger than in the M BH − σ * relation, the fact that the scatter in the M BH − V c relation is no smaller than that in the M BH − σ * relation still requires (baryonic) mechanisms beyond the pure merging scenario.
On the other hand, if AGN feedback is important in regulating BH growth, and it acts primarily on the bulge scale (Debuhr et al. 2010), the smaller intrinsic scatter in the M BH − σ * relation could be explained. We speculate that on top of the correlations formed by merging, feedback further regulates the BH mass according to the bulge potential σ * , tightens the M BH − σ * relation (not the M BH − V c relation), and decreases the intrinsic scatter in a relative sense. Alternatively, other baryonic mechanisms that connect the BH to the bulge, e.g., via feeding of bulge stars onto the accretion disk (Miralda-Escudé & Kollmeier 2005), or star formation-regulated BH growth (Burkert & Silk 2001), could also contribute to the tighter M BH − σ * relation.

M BH Upper Limits in the Two Scaling Relations
M BH upper limits in bulgeless galaxies have also been used to differentiate the dependence of black hole mass on galaxy halos or bulges. Kormendy & Bender (2011) found that M BH cannot be uniquely determined by V c . Bulgeless galaxies have low (or zero) M BH and low σ * but relatively high V c . Therefore, σ * is a better predictor of M BH than is V c for these galaxies.
We expand the sample of bulgeless galaxies with M BH upper limits by including the measurements from Neumayer & Walcher (2012) in our secondary sample. This sample is described in Appendix D and listed in Tables 8 and 9. As can be seen from Figure 8, there are a few upper limits lying below the M BH − V c relation but not the M BH − σ * relation. With a larger number of galaxies, we confirm Kormendy's statement that the halo mass alone does not determine the black hole mass. However, it is also possible that these bulgeless galaxies do not host a massive black hole at all. The presence or absence of a BH seed may involve different physical mechanisms than those that couple BHs to galaxies (e.g., Volonteri et al. 2011). In any case, the upper limits strengthen our conclusion that bulge scales matter, likely both in seeding and in BH growth with cosmic time.

SUMMARY
We refine the measured M BH − V c relation and compare it to the M BH − σ * relation to gain insight into the mechanisms that drive BH scaling relations. We perform VLA observations to measure the circular velocities V c for five galaxies with dynamical M BH . Together with a thorough literature search, we increase the sample size of galaxies with both M BH and V c measurements to thirtythree. Twenty-two of these have V c evaluated at large enough radius (> R 25 /2) to be reliable for our scientific purpose.
With this sample, we constrain the power-law M BH −V c relation to have an intercept α = 7.43 +0.13 −0.13 , slope β = 3.68 +1.23 −1.20 , and intrinsic scatter int = 0.51 +0.11 −0.09 . The intrinsic scatter in the M BH −V c relation is not significantly smaller than that in the M BH − σ * relation, showing that M BH does not correlate better with dark matter halo mass than with bulge properties. This contradicts naive expectations from pure merging scenarios. Furthermore, we consider a number of M BH upper limits in bulgeless galaxies that lie significantly below the M BH − V c relation, suggesting that BH masses are better predicted by the bulge, not the halo, either via its seeding or accreting mechanism. We thus suggest that pure merging is not likely to be the only mechanism that drives the scaling relations. AGN feedback may also be an essential ingredient to tighten the correlation between BH mass and bulge properties.
We highlight possible future improvements to this work. First, modeling of the dark matter halo mass distribution, to separate the halo mass from the baryons using both rotation curves and light distributions, can improve the halo mass estimation and better constrain the M BH − M DM relation. Second, V c measured by stellar dynamical modeling for elliptical galaxies, if proven to be comparable to those from H I rotation curves, can improve the dynamic range and better constrain the slope of the relation. Third, quantitative theoretical predictions for the scaling relations, especially the intrinsic scatter, from both feedback and pure merging scenarios, will enable more direct comparisons with the empirical relations. Ultimately, a large sample unbiased with respect to M BH and V c , which awaits next-generation telescopes, would be the most ideal data set for this study.
We thank the referee for a very thorough report that helped us improve this manuscript. We thank K. Gültekin, L. C. Ho, N. McConnell, and K. Jahnke for helpful discussions, M. Strauss for constructive suggestions, and J. Kormendy for providing valuable references. We acknowledge the usage of the HyperLeda database (http://leda.univ-lyon1.fr) and NASA/IPAC Extragalactic Database (http://ned.ipac.caltech.edu).

A. H I PROPERTIES OF INDIVIDUAL GALAXIES
In addition to the rotation curves kinematics, there is also rich information about the H I gas properties in our VLA data cubes. We discuss the H I fluxes and masses and the H I properties of each individual galaxy in this section. H I moment maps of these galaxies are shown in Figures 1 to 5 and the integrated spectra of lower SNR galaxies are shown in Figure 9.

A.1. H I Flux and Mass
The H I fluxes and mass are tabulated in Table 7. We measured the H I flux by integrating the data cube over the spatial regions as shown in the moment-0 maps. To avoid contamination, for NGC 1194, we excluded the northwestern cloud, which is detached from the main galaxy as seen from the moment-0 map. We also removed the central absorption region in NGC 7582. As an estimate of the H I flux error, we also measured the flux from the masked moment-0 maps, which exclude the noise-dominated channels and therefore can underestimate the total flux. We take the difference between the two as the systematic flux error. The H I masses are calculated using the equation from Draine (2011). For consistency, we adopt the same distances as those used for the M BH measurements in §4.1.2 and Table 2. Distance uncertainties of 10% are assumed, and the uncertainties in M HI are propagated from the uncertainties in flux and distance. It is common in the literature to compare the H I mass to the total stellar mass, to understand the gas fractions and available fuel for future star formation. Here we use the B-band luminosity (L B ) as a proxy for the stellar mass and calculate the H I mass-to-light ratio (M HI /L B ) as a proxy for the H I-to-stellar mass ratio. Using the B-band luminosity from HyperLeda (Paturel et al. 2003), the M HI /L B ratios of our galaxies range from 3% to 22% (Table 7). This is similar to the range of typical Sa galaxies, which have an average value ∼ 10% with a factor of ∼ 3 dispersion (Roberts & Haynes 1994), as well as AGN hosts which also have M HI /L B ∼ 10% (Fabello et al. 2011;Ho et al. 2008). Therefore, our sample galaxies have roughly similar H I-to-stellar mass ratios as typical disk galaxies and AGN hosts.
A.2. NGC 2748 NGC 2748 is an SAbc galaxy ) at a distance of 21.6 ± 1.4 Mpc. The H I observation of this galaxy has the highest SNR among our five galaxies, with an H I-based inclination of 72.6 • , compared to 68.1 • from optical data in HyperLeda (Paturel et al. 2003). The fitted H I systemic velocity is 1482 km s −1 . The kinematics show the typical signature of a rotating disk (Fig. 1). The asymmetric velocity field suggests that slight disk warping may exist. The angular diameter along the major axis at the 3 σ level, i.e., where H I is detected above the 3 σ level, is 3. 8. The H I flux is measured to be 37 ± 5 Jy km s −1 , which corresponds to an H I mass of 4.2 ± 0.9 × 10 9 M .
A.3. NGC 7582 NGC 7582 is an SB(s)ab galaxy ) at a distance of 20.6 ± 2.4 Mpc. Our observation of this galaxy also has high SNR. The inclination derived from the H I kinematics is 67.9 • , compared with 68.2 • from HyperLeda (Paturel et al. 2003). The systemic velocity is 1588 km s −1 as fit by H I. Its velocity field shows the clear signature of rotation and is asymmetric (Fig. 2), suggesting a warped disk. The velocity discontinuity at the edge of the warp coincides with the location of the spiral arms. Its angular diameter is 3. 7 at the 3 σ level. The H I flux is 20 ± 3 Jy km s −1 , which corresponds to an H I mass of 2.0 ± 0.5 × 10 9 M .
NGC 7582 shows H I absorption features against the central continuum source. The continuum-subtracted spectrum averaged over the central beam is shown in Figure 9. There is an emission peak at a velocity of 1500 km s −1 coinciding in velocity with the emission in nearby regions. Therefore, we suspect that we are observing the superposition of a wide absorption feature with narrow emission. The wide absorption feature has a FWHM of ∼ 400 km s −1 centered at 1580 km s −1 , close to the systemic velocity. Because of contamination from H I emission, we can only estimate a lower limit on the absorbed flux. Using the line width and depth, that flux limit is > 0.94 Jy km s −1 . The lower limit on the H I absorption optical depth is estimated to be > 0.027.
Finally, there are two H I companions of NGC 7582 observed in the same velocity range and sitting 9 and 12 to the north-east respectively. These are the companion galaxies NGC 7590 and NGC 7599. There is faint diffuse H I emission ∼ 13 long that extends from NGC 7590 to the west of NGC 7582 with a closest distance to NGC 7582 of 5. 5. This is possibly tidally stripped gas due to interactions between these galaxies. Larger scale moment-0 and moment-1 maps showing these structures are plotted in Figure 10.
A.4. NGC 1194 NGC 1194 is an SA0 galaxy ) at a distance of 53.2 ± 3.7 Mpc. It is one of our water megamaser galaxies (Kuo et al. 2011). The H I SNR is sufficiently high to determine an inclination of 69.1 • , compared to 71.1 • from HyperLeda (Paturel et al. 2003). The systemic velocity as measured from the H I map is 4075 km s −1 . The H I moment-0 map shows an elongated morphology with angular diameter 2. 8 at the 3 σ level, showing that the H I gas is organized within the galaxy disk. The velocity field shows a clear velocity gradient due to rotation (Fig.  3). The integrated spectrum (Fig. 9) also has a clear double peaked rotation signature with the width consistent with our V c measurement. The H I flux is measured to be 7 ± 1 Jy km s −1 , corresponding to a mass of 5 ± 1 × 10 9 M . There is one H I cloud to the north-west side of the galaxy that is detached from the main galaxy body (see Fig. 3). This cloud has a mass of 5 ± 2 × 10 8 M .
A.5. NGC 2960 NGC 2960, also called Mrk 1419, is an Sa ) galaxy at 72.2 ± 5.1 Mpc and also a water megamaser galaxy from Kuo et al. (2011). The systemic velocity from H I is 4939 km s −1 . Our map of NGC 2960 has lower SNR (SNR = 6). Although H I emission is only seen in discrete patches of the galaxy, an overall velocity gradient is observed (Fig 4). The detected diameter is 2. 5 at the 1 σ level from the moment-0 map. The inclination cannot be constrained by the H I data alone, and we derive the rotation curve by adopting the optical inclination from HyperLeda (Paturel et al. 2003) of 41.5 • . The integrated spectrum (Fig. 9) shows a double peaked rotation signature and a width that is consistent with our V c measurement. The H I flux is measured to be 2.2 ± 0.3 Jy km s −1 , corresponding to a mass of 2.7 ± 0.5 × 10 9 M . Previous D-array VLA observations (Kuo et al. 2008) are consistent with the spatial extent and velocity from our data, but their H I flux of 1.7 ± 0.3 Jy km s −1 is lower than ours. This is understandable as their flux was calculated from the moment-0 map and should be lower than from the data cube.

A.6. Jet in NGC 2960
Our H I data also provides images of the radio continuum at 20 cm. We find that the radio continuum of NGC 2960 is slightly extended, suggesting that there is a radio jet launched from the central black hole (Fig. 11). The continuum image is made from line-free channels on both sides of the H I line with a total velocity range of 300 km s −1 , and is cleaned to the 5 σ level using robust = 0.5 weighting. We measure the size and position angle of the jet using two different methods. First, we fit a central point source (15. 7) and subtract it. We detect residual emission to the south-east. Taking into account the positional uncertainty due to the finite beamsize, this extended emission is 20 ± 3 away from the center at a position angle of 125 • . The flux in the extended emission is 1.6 ± 0.2 mJy. To estimate the size and P.A. error we also try a second method to quantify the elongated structure. We fit the whole continuum image with a two-dimensional Gaussian while allowing the semi-major and semi-minor axis, central position, and P.A. to vary, without deconvolution of the beam. The fitted Gaussian has a major axis of 22 ± 4 , and a minor axis of 18 ± 4 with P.A. 119 ± 7 • , consistent with the first method. Therefore we conclude that the jet size is 20 ± 5 at a P.A. of 125 ± 10 • . As in all megamaser galaxies studied to date, the jet axis is coincident with the rotation axis of the maser disk (Greene et al. 2013).
A.7. UGC 3789 UGC 3789 also hosts a water maser disk, is an SA(r)ab galaxy , and is at a distance of 49.9 ± 7.0 Mpc . The H I systemic velocity is 3229 km s −1 . The H I observation also has low SNR = 6 and shows a ring-like structure with a diameter of 1. 3 at the 1 σ level on the moment-0 map (Fig. 5). There is a clear velocity gradient. Neither inclination nor position angle can be constrained from the H I data alone. The optical inclination i = 43.2 • and position angle P.A. = 164.7 • (Paturel et al. 2003) are adopted to derive the rotation curve. The measured H I flux is 1.3 ± 0.4 Jy km s −1 , corresponding to a mass of 0.7 ± 0.3 × 10 9 M . Companion UGC 3797, which is 5 away to the east, is detected in H I at the same angular distance but we do not detect any H I tidal tails.

B. DETAILS OF ROTATION CURVES ESTIMATION
In this section we investigate two complications that may affect rotation curve estimation. First, we consider differences in rotation curves arising from different methods of assigning velocity. Second, we look at tilted ring modeling, which captures warps in the rotation disk, to see if our circular velocities are changed.
As discussed in §3.2, there are various methods to assign velocity to a spectrum, e.g., the peak velocity, the intensityweighted mean (moment-1), or Gaussian/Gauss-Hermite velocities. The first two methods have some drawbacks. The peak velocity is sensitive to the noise in the data, and is discretized to the channel width of the data cube, and the moment-1 velocity is biased towards the wing if the spectral line is asymmetric. Our preferred method is therefore the Gaussian/Gauss-Hermite fit. The averaged difference between the Gaussian/Gauss-Hermite and the peak velocity (moment-1) velocity field is 17 (18) km s −1 for the high SNR galaxies (NGC 2748, NGC 7582, andNGC 1194) and 37 (33) km s −1 for low SNR galaxies (NGC 2960 and UGC 3789).
We now propagate these different velocity assignments to investigate the differences in our inferred rotation curves. We focus on two representative galaxies, NGC 7582 (high SNR = 23) and NGC 2960 (low SNR = 6). We use the same MCMC procedure as discussed in §3.3, assuming a homogeneous error in the velocity field of 20 km s −1 for NGC 7582 and 40 km s −1 for NGC 2960, and applying the same masking as in the Gaussian/Gauss-Hermite velocity field to mask out noisy pixels. The fitted rotation curves are shown in Fig. 12. The average differences between Gaussian/Gauss-Hermite and the peak velocity (moment-1) rotation curves are 11 (23) km s −1 for NGC 7582 and 67 (39) km s −1 for NGC 2960. Compared to the RMS errors of the Gaussian/Gauss-Hermite rotation curves, which are 24 km s −1 for NGC 7582 and 37 km s −1 for NGC 2960, we found that: First, the peak velocity rotation curve of NGC 7582 is consistent with the Gaussian/Gauss-Hermite rotation curve within the error. However, the peak velocity rotation curve deviates from the Gaussian significantly for the lower SNR galaxy NGC 2960.
These tests confirm our concern that the peak velocity is highly sensitive to noise and is not suitable for estimating the rotation curves for low SNR data. Also, NGC 7582 has a moment-1-derived rotation curve that is systematically lower than the Gaussian/Gauss-Hermite one, with a deviation comparable to the RMS error. This is consistent with our expectation that the moment-1 velocity is biased by the wing of the spectral line. Moment-1 velocities usually underestimate the rotation velocity for high SNR data, due to beam smearing and higher intensity at the inner part of the galaxy. This bias was also demonstrated in the left panel of Fig. 6. NGC 2690 has a moment-1 rotation curve that differs from the Gaussian rotation curve at a level that is comparable to the RMS errors as well. In this case we see no systematic bias, possibly because the SNR is not high enough to manifest an asymmetric line shape. Considering these effects, we find that the Gaussian/Gauss-Hermite method is relatively reliable in the face of asymmetric line shapes and noisy data.
In §3.3 we adopt a coplanar disk model for rotation curve fitting. However, this model does not capture the asymmetry in the velocity field that arises from warping of the disk, as seen in NGC 2748 and NGC 7582. Here we assess whether these asymmetries have an effect on our V c estimates. We take our most asymmetric galaxy NGC 7582 as an example, and use the Kinemetry method (Krajnović et al. 2006) to fit a tilted-ring model. In this model, the position angle and the inclination of each ring are allowed to vary. Furthermore, second and third order harmonic terms, e.g., sin(2ψ) and sin(3ψ), are included to capture higher frequency variations along each ring. The best-fit models are plotted in Fig. 13 with their residuals, and the fitted rotation curves are plotted in Fig. 14. The tilted-ring model better captures warps in the disk and reduces the residuals, but the best-fit rotation velocities stay unchanged within the uncertainties. The difference at the outer-most bin is less than 5% of the RMS error. This demonstrates the robustness of the rotation curve fitting against warps and higher order variations of the velocity field, and shows that the RMS error is a conservative estimate of the potential systematics in the rotation curve fitting.

C. Vc FROM SINGLE DISH MEASUREMENTS
In our secondary sample we have three galaxies (NGC 3368,NGC 3393,NGC 3489) with dynamical M BH but no available rotation curves. However, V c can also be inferred from the line width of integrated H I spectra taken with single dish radio observations (V c,SD ), and the single dish V c for these three galaxies from HyperLeda (Paturel et al. 2003) are listed in Table 8.
Since single-dish measurements are more readily available for large samples, they have been used in previous scalingrelation studies (Pizzella et al. 2005;Courteau et al. 2007;Ho 2007;Beifiori et al. 2012). Roberts (1978) and Ho (2007) find that V c,SD is a robust substitute for V c,RC . On the other hand, without spatial information, single-dish circular velocities may contain large uncertainties due to the distribution of atomic gas, irregular rotation-curve shape, or contamination from companion galaxies. In general, the values may skew towards lower values, since single-dish measurements are biased to the inner part of the rotation curve.
To understand how much V c,SD can deviate from V c,RC , we compare the two numbers for all galaxies in our sample that have both measurements (Figure 15). For most of the galaxies, V c,SD is consistent with V c,RC . However for a few of them the two numbers can deviate by up to a factor of two. These are preferentially S0 galaxies, suggesting that in these cases low gas fractions are skewing the V c,SD values low. On the other hand, there are only three galaxies in our sample that only have single-dish H I measurements. In practice, including or excluding these three galaxies from our fitting does not change the result. We decide not to include them in our primary sample.

D. M BH UPPER LIMITS
For completeness, dynamically constrained M BH upper limits for bulgeless galaxies, mostly from Neumayer & Walcher (2012), are listed in the second section of Table 8 and 9. The V c values are assigned as described in §4.1.1, except for the single-dish V c , which are from HyperLeda (Paturel et al. 2003). These upper limits are also plotted in Figure 8 as grey triangles in both of the scaling relations. Some of them are outliers in the M BH − V c relation (see discussion in §6.2). The NGC 2748 moment-0 map made by the CASA image.moments tool. The masking is described in §3.1. We mask out the pixels below 3 σ in a map that is smoothed spatially and over two velocity channels (20 km s −1 ). We then construct the moment-0 map from the original data using this mask. The ellipse in the bottom left corner represents the beam. Right: NGC 2748 Gaussian/Gauss-Hermite fitted velocity field. Velocity is in optical LSRK. The NGC 7582 moment-0 map. The masking is described in §3.1. We mask out the pixels below 3 σ in a map that is smoothed spatially and over two velocity channels (20 km s −1 ). We then construct the moment-0 map from the original data using this mask. The ellipse in the bottom left corner represents the beam. Right: The NGC 7582 Gaussian/Gauss-Hermite fitted velocity field. The central region is masked because of the absorption feature. The NGC 1194 moment-0 map. The masking is described in §3.1. We mask out the pixels below 2.5 σ in a map that is smoothed spatially and over two velocity channels (40 km s −1 ). We then construct the moment-0 map from the original data using this mask. The ellipse in the bottom left corner represents the beam size. Right: The NGC 1194 Gaussian/Gauss-Hermite fitted velocity field. The NGC 2960 moment-0 map. The masking is described in §3.1. We mask out the pixels below 2.5 σ in a map that is smoothed spatially and over two velocity channels (40 km s −1 ). We then construct the moment-0 map from the original data using this mask. The ellipse in the bottom left corner represents the beam. Right: The NGC 2960 Gaussian fitted velocity field. The UGC 3789 moment-0 map. The masking is described in §3.1. We mask out the pixels below 2.5 σ in a map that is smoothed spatially and over two velocity channels (40 km s −1 ). We then construct the moment-0 map from the original data using this mask. The ellipse in the bottom left corner represents the beam size. Right: The UGC 3789 Gaussian fitted velocity field.  -Rotation curves of the five galaxies. Red error bars are the RMS errors, which dominate the error budget in all cases. Yellow is the MCMC fitting error. Green/blue is the error contribution from the inclination/position angle uncertainty. Black error bars are the total error, calculated by the quadratic sum of all the error sources. The circular velocity Vc as measured from the outer-most radial bin is denoted with the black horizontal lines with the errors marked with the grey shaded region. This Vc error includes the observational error and the rotation curve variation. The R 25 radii are marked as the solid black vertical lines and dashed black horizontal lines are the R 25 /2. All of the five galaxies have rotation curves extending beyond R 25 /2, and NGC 2748, and NGC 1194, and NGC 2960 are beyond R 25 . For NGC 7582 and NGC 1194 the last two bins, linked by the dark gray lines, are noisy bins and thus not used for the Vc measurements.   Table 5. Our primary sample with dynamical M BH and spatially resolved Vc measurements is plotted in blue or green depending on the rotation curve extent. The blue circles (Ro > R 25 ) and blue squares (R 25 > Ro > R 25 /2) have long rotation curves and are used for constraining the two relations. The green circles/squares have short rotation curves (Ro < R 25 ) with flat/rising trends. They are not used to constrain the two relations because of their lower reliability reflected by the larger scatter in the M BH − Vc relation. Also plotted is our secondary sample described in Appendix C and listed in Table 9. The grey diamonds and triangles represent the single-dish Vc measurements and M BH upper limits respectively. The red circles are the elliptical galaxies with Vc measured by dynamical modeling and M BH from McConnell & Ma (2013) as discussed in §4.3. Our five observed H I galaxies are marked by magenta dots labeled as A (NGC 2748), B (NGC 7582), C (NGC 1194), D (NGC 2960), and E (UGC 3789). The fitted M BH − Vc and M BH − σ * scaling relations using the Ro > R 25 /2 primary sample, plotted in blue, is shown by the thick blue line with the intrinsic scatter plotted by the two thin lines. We use √ 2σ * for the abscissa of the M BH − σ * relation as a direct comparison to Vc.   Figure 2. We mask out the pixels below 3 σ in a map that is smoothed spatially and over two velocity channels (20 km s −1 ). We then construct the moment-0 map from the original data using this mask. Two companions, NGC 7590 and NGC 7599, and an elongated tidal stream can be seen. Right: The moment-1 map of the same field. -Comparing the rotation curves derived from different velocity fields -Gaussian/Gauss-Hermite (blue), moment-1 (green), and peak velocity (red). The left is NGC 7582 with high SNR = 23, and the right is NGC 2960 with lower SNR = 6. The thick error bars represent the fitted error and the thin error bars are the RMS error. NGC 7582 has Vc, the rotation velocity at large radius, from the three velocity fields consistent with each other. However, NGC 2960 has peak-velocity Vc significantly differs from the other two methods. Therefore, the peak velocity may incur large errors in the rotation curves for low SNR data. Fig. 13.-Comparing coplanar (top) and tilted-ring (bottom) models for the velocity field of NGC 7582. From left to right we show the data (Gauss-Hermite velocity field), the best-fit model, and the residuals. The coplanar model is described in §3.3 and the tilted-ring model in Appendix B. The color bar represents the color scheme of the residuals. While the tilted-ring model captures the warping feature in the velocity field, we show in Fig. 14 (5): The angular separation between the source (galaxy) and the phase calibrator. Col. (6): Total on-source observation time. Col. (7): Average length of each source scan, which is the separation between two phase calibrator scans. Col. (8): Number of antennas used in the observation. Some antennas were not used because the L-band receiver was not yet installed or because the antenna had unstable or noisy data quality. Col. (9): Whether or not radio frequency interference (RFI) was found in the data. The RFI was visually inspected and flagged. After the flagging, there was negligible or minor contamination from RFI in the data cube. The most severe case was NGC 1194, where some faint elongated stripes parallel to the galaxy can be seen. Note.   (NGC 1194, NGC 2960, andUGC 3789) the center is fixed at the maser position (Kuo et al. 2011). Col. (4): Fitted systemic velocity. The random error from fitting is less than 0.1%, including determining the center for NGC 2748 and NGC 7582, which have no megamaser positions. Col. (5): Fitted inclination angle, except for NGC 2960 and UGC 3789, which are fixed at HyperLeda values (Paturel et al. 2003). Col. (6): Fitted position angle, except for UGC 3789, which is fixed at HyperLeda value (Paturel et al. 2003).     Table 2. A distance error of 10% is assumed. Col. (4): Absolute B-band magnitude corrected for extinction from HyperLeda (Paturel et al. 2003). Col. (5) Table 8). Col. (9): The observational method of Vc. SD stands for spatially unresolved single dish Vc, RC for H I rotation curve Vc, H I V-field for spatially resolved data but unavailable rotation curves, and optical for optical rotation curves. References: (1)