The MBHBM⋆ Project. I. Measurement of the Central Black Hole Mass in Spiral Galaxy NGC 3504 Using Molecular Gas Kinematics

Dieu D. Nguyen; Mark den Brok; Anil C. Seth; Timothy A. Davis; Jenny E. Greene; Michelle Cappellari; Joseph B. Jensen; Sabine Thater; Satoru Iguchi; Masatoshi Imanishi; Takuma Izumi; Kristina Nyland; Nadine Neumayer; Kouichiro Nakanishi; Phuong M. Nguyen; Takafumi Tsukui; Martin Bureau; Kyoko Onishi; Quang L. Nguyen; Ngan M. Le

doi:10.3847/1538-4357/ab77aa

1. Introduction

Supermassive black holes (SMBHs, ${M}_{\mathrm{BH}}\gt {10}^{6}{M}_{\odot }$ ) have been found to reside at the centers of massive galaxies, and their masses correlate with macroscopic properties of the host (e.g., bulge velocity dispersion (σ), bulge mass (M_bulge), and bulge luminosity (L_bulge). These remarkable discoveries are based on large observational efforts (e.g., Kormendy & Richstone 1995; Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Graham et al. 2001; Marconi & Hunt 2003; Häring & Rix 2004; Gültekin et al. 2009; Beifiori et al. 2012; Kormendy & Ho 2013; McConnell et al. 2013; Saglia et al. 2016; van den Bosch et al. 2016) and theoretical studies of self-regulated mechanisms of active galactic nuclei (AGN) feedback onto the outer gas reservoirs (Silk & Rees 1998; Di Matteo et al. 2008; Fabian 2012; Silk & Mamon 2012; Barai et al. 2014; Netzer 2015; Naab & Ostriker 2017). The effect of this feedback is recorded in the M_BH–galaxy scaling relations (e.g., Meza et al. 2003; Vogelsberger et al. 2014). Feedback can also shape the galaxy stellar-mass function (e.g., Croton 2006; Schaye et al. 2015) and the star formation history (SFH; Martín-Navarro et al. 2018) and metallicity (e.g., Choi et al. 2017) of individual galaxies. These findings suggest SMBHs may play a pivotal role in the growth and evolution of galaxies (e.g., Schawinski et al. 2007; Kormendy & Ho 2013; McConnell et al. 2013; Saglia et al. 2016; van den Bosch et al. 2016).

However, our understanding of the scaling relations in lower-mass galaxies is lacking. This includes galaxies near the break of the galaxy stellar-mass function ${M}_{\star }\sim 5\times {10}^{10}$ (e.g., Baldry et al. 2012) like the Milky Way (MW), as well as galaxies at lower masses. Measurements of black hole (BH¹⁸ ) masses are lacking in these galaxies. Despite the MW's precisely measured central BH mass of ${M}_{\mathrm{BH}}$ $\sim \,4\times {10}^{6}\,{M}_{\odot }$ (Boehle et al. 2016; Gillessen et al. 2017), we actually know very little about the demographics of BHs in galaxies of this mass. In particular, there are only 11 published ${M}_{\mathrm{BH}}$ in spiral galaxies within ∼1 dex in stellar mass of the MW, but there are strong indications that the relations between ${M}_{\mathrm{BH}}$ and galaxy properties that hold at higher mass (Kormendy & Ho 2013; McConnell et al. 2013; Saglia et al. 2016; van den Bosch et al. 2016) break down for M_⋆ spiral galaxies (Greene et al. 2010, 2016; Läsker et al. 2016). For instance, Andromeda is only ≲2 times more massive than the MW, but its SMBH is two orders of magnitude higher, ${M}_{\mathrm{BH}}$ ∼1.4 × 10⁸M_⊙ (Bender et al. 2005). At even lower galaxy masses, the changes in the inferred ${M}_{\mathrm{BH}}$ are twofold: not only does the scatter increase enormously at low mass (e.g., Reines et al. 2013), but it also seems that the ${M}_{\mathrm{BH}}$ are lower at fixed stellar mass or σ of the lower-mass systems (e.g., Kormendy & Ho 2013; Scott et al. 2013; Greene et al. 2016; Saglia et al. 2016; Nguyen et al. 2017, 2018, 2019, hereafter N17, N18, N19). Despite the incomplete sample of ${M}_{\mathrm{BH}}$ at lower masses, there are already two intriguing findings from the limited data available: (1) a factor of 2 difference in normalization seen in the ${M}_{\mathrm{BH}}$ –σ correlations for early-type (ETGs) and late-type (LTGs) galaxies (McConnell et al. 2013) and (2) a large scatter up to two orders of magnitude of M_BH seen around the global scaling relations for low-mass systems (Greene et al. 2010; Scott et al. 2013; Graham & Scott 2015; Greene et al. 2016; Läsker et al. 2016; Chilingarian et al. 2018).

Currently, the ${M}_{\mathrm{BH}}$ measurements in both samples of M_⋆ and sub-M_⋆ galaxies are small. Apart from 11 ${M}_{\mathrm{BH}}$ measurements in spiral MW-like galaxies from water masers (Greene et al. 2016; Gao et al. 2017; Zhao et al. 2018), recent progress in the stellar/ionized gas dynamics and the virial ${M}_{\mathrm{BH}}$ estimate techniques have led to the detections of 28 BHs with solar masses below 1 M in the sub-M_⋆ sample, whose ${M}_{\mathrm{BH}}$ are a remarkably tiny fraction of the galaxy stellar masses (Seth et al. 2010; Reines et al. 2013; Nguyen et al. 2014; Baldassare et al. 2015; den Brok et al. 2015; Nguyen 2017; Thater et al. 2017; Chilingarian et al. 2018; Thater et al. 2019; N17, N18, N19), suggesting that these objects fall below the ${M}_{\mathrm{BH}}$ –M_bulge relation extrapolated from higher-mass galaxies in a wide range of 1–3 orders of magnitude of M_BH (N19). Many hypotheses have been proposed to explain this change, which may be due to (1) the formation history of the bulge (Kormendy & Bender 2012; Krajnović et al. 2018), (2) the SFH of the galaxy (Caplar et al. 2015; Terrazas et al. 2017), or (3) the bimodality of accretion efficiency of SMBH seeds (Pacucci et al. 2015, 2017, 2018; Inayoshi & Haiman 2016; Park et al. 2016).

There are also disturbing hints that, at fixed galaxy stellar velocity dispersion or bulge stellar mass, the distribution in ${M}_{\mathrm{BH}}$ may differ depending on the methods used to determine the dynamical ${M}_{\mathrm{BH}}$ , with ${M}_{\mathrm{BH}}$ measurements from water masers seeming to be systematically lower than stellar dynamical measurements (see Figure 1 of Greene et al. 2016). Additionally, only a few galaxies have multiple dynamical estimates: the differing ${M}_{\mathrm{BH}}$ for M87 (Gebhardt et al. 2011; Walsh et al. 2013) inferred from stellar versus ionized gas kinematics, which has only just been settled by the Event Horizon Telescope observations, highlights the importance of comparing multiple measurements for the same BH (Event Horizon Telescope Collaboration et al. 2019a, 2019b, 2019c). The lack of measurements in M_⋆ mass spirals and sub-M_⋆ galaxies in general is in part due to the complexity of the stellar dynamical modeling in these systems (e.g., Thomas et al. 2016); independent dynamical tracers are crucial to understand the systematic scatter and the true distribution in the ${M}_{\mathrm{BH}}$ measurements at the low-mass regime. Thus, we turn to Atacama Large Millimeter/submillimeter Array (ALMA) observations, to exploit the gas-rich nature of these galaxies to determine their ${M}_{\mathrm{BH}}$ .

ALMA observations of molecular gas at millimeter/submillimeter wavelengths offer a promising way to characterize the full spectrum of BH populations across the Hubble sequence from LTGs to ETGs, both active and inactive (Davis et al. 2013, 2017, 2018; Davis 2014; Onishi et al. 2015, 2017; Barth et al. 2016a, 2016b; Yoon 2017; Boizelle et al. 2019; Combes et al. 2019; Nagai et al. 2019; North et al. 2019; Smith et al. 2019; Thater 2019; T. Davis et al. 2020, in preparation; D. Nguyen et al. 2020, in preparation) because of its high angular resolution and sensitivity. Moreover, ALMA can improve on many problems affecting ${M}_{\mathrm{BH}}$ estimates via existing optical/IR instruments including (1) high angular resolution capable of resolving the sphere of influence (SOI) and (2) the ability to obtain dynamical measurements in dusty/obscured nuclei, which are inaccessible at optical wavelengths. The physical idea behind the cold gas-dynamical method is that the ${M}_{\mathrm{BH}}$ is derived by detecting and resolving the Keplerian motion of the gas disk at the galactic center directly (within the SOI). This method is only possible in galaxies that host well-defined and rotating circumnuclear gas disks (CNDs).

This article is the first of a series of the Measuring Black Holes in below Milky Way (M_⋆) Mass Galaxies (MBHBM_⋆) Project. This project aims to gather a large sample of low-mass gas-rich galaxies and then measure their central dark masses, which are likely BHs, using molecular gas tracers (e.g., ¹²CO(2 − 1)) observed with ALMA. In this project, we select targets based on the presence of well-defined CNDs of gas and dust, which serve as morphological evidence for rotating dense gas about galaxy centers based on previous low–spatial resolution surveys. Seven targets for ${M}_{\mathrm{BH}}$ measurements have already been observed by ALMA in cycle 5. Here we start the project with the ${M}_{\mathrm{BH}}$ measurement for NGC 3504, the first galaxy in the sample that has been observed.

The paper is organized into eight sections. We summarize the properties and determine a new distance to the galaxy NGC 3504 in Section 2. In Section 3, we present the Hubble Space Telescope (HST) images and ALMA observations of ${}^{12}\mathrm{CO}(2-1)$ nucleus gas and data reduction. We also report the evidence of a dense gas tracer CS(5 − 4) at the center of the galaxy in this section. The mass modeling of NGC 3504 is discussed in Section 4. We model the ${}^{12}\mathrm{CO}(2-1)$ gas disk and estimate the central ${M}_{\mathrm{BH}}$ and uncertainties via a kinematic molecular simulation (KinMS; Davis et al. 2013) model and a tilted-ring model (e.g., Neumayer et al. 2007; den Brok et al. 2015) in Sections 5 and 6, respectively. We discuss our results in Section 7 and conclude in Section 8. Throughout the paper, unless otherwise indicated, all quoted quantities have been corrected for a foreground extinction A_V = 0.072 (Schlafly & Finkbeiner 2011) using the interstellar extinction law of Cardelli et al. (1989).

2. The Galaxy NGC 3504

2.1. Determination of the Distance

The distance measurements for NGC 3504 in the literature span a wide range from 9 to 26 Mpc from Tully–Fisher measurements (Bottinelli et al. 1984; Tully & Fisher 1988; Russell 2002). When we began the MBHBM_⋆ project (Nguyen 2017) and defined the target sample, we assumed a distance of 14 ± 5 Mpc for NGC 3504 based on Russell (2002). Due to the significant distance uncertainty (and because this has a large impact on our BH mass estimates), we concluded that the systematic uncertainty was 9 Mpc, which is a large fraction of the total distance. We therefore decided to use the existing HST data to attempt to derive a distance estimate using the surface brightness fluctuation (SBF) technique.

The weak constraint on the NGC 3504 distance from existing distance measurements directly affects the uncertainty in the calculated ${M}_{\mathrm{BH}}$ (M_BH ∝ D; see Section 5.3.2, e.g., de Vaucouleurs et al. 1981; Davis 2014; Davis et al. 2017). To reduce this uncertainty, we determined the SBF distance to this galaxy using the methodology outlined in detail by Jensen et al. (2015), updated with current best procedures (Cantiello et al. 2018). The SBF technique takes advantage of the Poisson variations in the number of red giant stars from pixel to pixel (convolved with the instrumental point-spread function [PSF]) to determine the distance-dependent fluctuation magnitude $\overline{m}$ . While galaxies can have similar surface brightnesses at a wide range of distances, distant galaxies will have a greater number of stars per pixel than nearby galaxies, and the resulting variance σ² will be smaller. Distant galaxies look smoother than nearby ones. When properly calibrated (using Cepheid distances to the Virgo and Fornax galaxy clusters, in this case), the fluctuation magnitude can be converted to a distance modulus $(\overline{m}-\overline{M})$ .

Using the HST WFC3/IR F110W and F160W images from program GO-12450 (Section 3), we reprocessed the images following the procedures used for SBF (Cantiello et al. 2018). Near-IR images are used for SBF because the luminosity fluctuations are dominated by red giant stars, which are much most luminous at near-IR wavelengths. The extinction due to clumpy dust in the target galaxy or in the MW is also significantly lower than at optical wavelengths. The SBF technique is typically used to determine distances to smooth, old, elliptical galaxies by measuring the luminosity variations in surface brightness due to the Poisson statistics of the discrete population of red giant stars. NGC 3504 is a poor SBF candidate because of extensive dust lanes and star formation in the spiral arms; usually these features add power to the spatial Fourier power spectrum, from which the SBF signal is measured. Experience shows that the SBF signal from spiral galaxies like NGC 3504 is biased to larger fluctuations and closer distances (Jensen et al. 2003, 2015). Nevertheless, the lack of other reliable distances for this galaxy made it worth the effort, and we were successful at measuring the SBF signal in the regions outside the bulge and between the spiral arms well enough to get a distance good to about 10%. Figure 1 shows the F110W image, the residual image with a smooth model of the galaxy subtracted, and the spatial Fourier power spectrum. The regions used for SBF analysis are indicated in the center panel, where the bumpiness can easily be seen for this relatively nearby galaxy. The fit to the spatial power spectrum from which the fluctuation power and $\overline{m}$ were measured for NGC 3504 is also shown. The SBF signal-to-noise ratio of the F110W fit is 26.

**Figure 1.** Left panel: the full *HST* WFC3/IR F110W image of NGC 3504 with a logarithmic stretch to better show the disk and spiral structure outside the nucleus. Center: residual image of NGC 3504 with a smooth model of the galaxy subtracted. The bumpiness in the regions between the spiral arms is the SBF signal arising from the Poisson noise of the stellar component. Contaminating sources, including star-forming regions, dust lanes, and globular clusters, have been masked. Right panel: the SBF measurement is a fit of the spatial Fourier power spectrum in two dimensions to the sum of a white noise component and a normalized power spectrum of the PSF (red dashed lines show the individual components; the blue line is the sum, and the data are shown with + symbols).
Download figure:
Standard image High-resolution image

The SBF distance calculation also requires an accurate g–z or J–H color measurement to account for stellar population age and metallicity variations in the empirical calibration (Jensen et al. 2015). The g–z color was measured using PanSTARRS images (Magnier et al. 2016), which were combined into a g–z color map and used to determine color values in the regions of the galaxy specifically used for SBF measurement. $J-H$ color maps were constructed from the HST F110W and F160W images. Because this galaxy has active star formation in spiral arms and extensive dust, the uncertainties in the color measurements were the largest contributors to the overall uncertainty in $\overline{m}$ and distance. The colors we measured in the SBF region were (g–z) = 1.31 (AB) from PanSTARRS, scaled to match the ACS (g–z) color calibration, and (J–H) = 0.257 from the WFC3/IR F110W and F160W photometry. These values are within the color calibration range recommended as being reliable by Jensen et al. (2015). The SBF distances derived from the (g–z) and (J–H) colors are highly self-consistent (differences of 0.73 and 0.23 Mpc for the F110W and F160W measurements, respectively), giving confidence that the colors have been measured accurately.

The distance to NGC 3504 using the F110W images is 34.8 ± 3.8 Mpc. The F160W distance is 30.7 ± 3.3 Mpc, which is consistent at the 1σ level and reflects the varying quality of the power spectrum fits and high variable extinction in this spiral galaxy. These values were determined using the PanSTARRS colors and are entirely consistent with the distances derived using J−H instead. The weighted average of the four independent distance measurements (SBF in two filters and using two colors) gives a final distance of 32.4 ± 2.1 Mpc. Given that clumpy dust and young stars bias the SBF measurement, it is wise to consider the SBF distance derived herein as a lower limit on the actual distance.

The SBF distance gives a physical scale of 157 pc arcsec⁻¹ assuming H₀ = 73 km s⁻¹ Mpc⁻¹, Ω_M = 0.27, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.73$ cosmology (Freedman et al. 2001; Freedman & Madore 2010). Although this galaxy is in a region of complex peculiar motions and is relatively nearby, this newly measured distance is roughly consistent with the distance inferred from the Hubble flow. The redshift for this galaxy corrected for Virgo, Great Attractor, and Shapley cluster infall motions is 2033 ± 30 km s⁻¹ (Mould et al. 2000 NASA/IPAC Extragalactic Database, NED¹⁹ ). The velocity we would predict using the SBF distance and assuming H₀ = 73 km s⁻¹ Mpc⁻¹ is 2365 ± 153 km s⁻¹, a 2σ difference, suggesting that the distance to NGC 3504 cannot be much larger than what we have measured and that any bias due to young populations or clumpy dust is insignificant.

2.2. Central Properties

NGC 3504 is a nearby double-barred spiral LTG, Hubble-type (R)SAB(s)ab (de Vaucouleurs et al. 1991) in the Leo Minor Group (de Vaucouleurs 1975).

Laurikainen et al. (2004) decomposed the bulge mass of NGC 3504 using a bulge-disk-bar decomposition model in the H–band image and found the bulge-to-disk ratio, B/D = 0.356, and the effective radius of the bulge, ${R}_{e,{\rm{b}}}\sim 25^{\prime\prime}$ . Salo et al. (2015) recently also used infrared (IR, 3.6 and 4.5 μm) images from the Spitzer Survey of Stellar Structure in Galaxies (S⁴G; Sheth et al. 2010) to decompose the bulge-disk-bar structure of the galaxy and found the total apparent magnitude of the bulge in the 3.6 μm band of ${m}_{{\rm{b}},3.6\mu {\rm{m}}}=11.59$ mag. Taking into account the distribution of the disk and bar within the region of dominant bulge and assuming $M/{L}_{3.6\mu {\rm{m}}}\sim 0.73$ (M_⊙/ ${L}_{\odot }$ ) (the Bell et al. (2003) (g-r)– $M/{L}_{K}$ relation gives M/L_K = 0.88 (M_⊙/ ${L}_{\odot }$ )), the bulge mass is estimated as ${M}_{\mathrm{Bulge}}=2.6\times {10}^{9}$ M_⊙. We also estimate the total stellar mass of NGC 3504 using the apparent magnitude in the B band from the RC3 catalog, m_B = 11.69 ± 0.16 mag (Corwin & Buta 1994) and our distance estimated in Section 2.1, which give M_B = −20.86 ± 0.30 mag and ${M}_{* ,\mathrm{total}}\sim {1.17}_{-0.26}^{+0.35}\times {10}^{11}{M}_{\odot }$ after accounting for the $M/{L}_{B}\approx M/{L}_{g}\sim 3.53$ (M_⊙/ ${L}_{\odot }$ ) estimated from the $(g-r)$ – $M/{L}_{g}$ relation of Bell et al. (2003). These masses are consistent with the masses estimated using the H-band HST image and the dynamical M/L_H found in Sections 4 and 5, respectively. Here, ${M}_{* ,\mathrm{total}}\gt {M}_{\star }$ because our new distance estimate to NGC 3504 is much further than the distance found in the literature (∼14 Mpc; Russell 2002) and the limited volume ( $\lesssim 20$ Mpc)³ of our selected sample for the MBHBM_⋆ Project when we proposed it for ALMA observations.

The galaxy center is determined at (R.A., decl.) = ( ${11}^{{\rm{h}}}{03}^{{\rm{m}}}{11}^{{\rm{s}}}.210$ , $27^\circ 58^{\prime} 21\buildrel{\prime\prime}\over{.} 00$ ) in the J2000-equatorial coordinate system and has a systemic velocity of ${v}_{\mathrm{sys}}=1525.0\pm 2.1$ km s⁻¹ (Epinat et al. 2008). Photometry shows the inclination between the line of sight (LOS) and the polar axis of the galaxy is i = 26 fdg 4²⁰ and oriented with a position angle (PA) = 150° (Paturel et al. 2000).

The Palomar spectroscopic survey measures the stellar velocity dispersion of the bulge as σ = 119.3 ± 10.3 km s⁻¹ (Ho et al. 2009), suggesting a central SMBH with ${M}_{\mathrm{BH}}\sim {3.0}_{-2.5}^{+7.0}\,\times {10}^{7}$ M_⊙ for this galaxy based on the M_BH − σ relation for LTGs (Kormendy & Ho 2013). On the other hand, with a total K-band apparent magnitude of 8.27 mag (Verga 2017), this provides a ${M}_{\mathrm{BH}}\sim {2.84}_{-2.5}^{+20.5}\times {10}^{7}$ M_⊙ for the SMBH (Kormendy & Ho 2013) after accounting for the Laurikainen et al. (2004) B/D ratio. These M_BH estimates are consistent with each other.

Both optical and radio data suggest the galaxy has a composite nucleus with both AGN and emission from star formation (Keel 1984). Hα obtained from a 0.9 m telescope at Kitt Peak National Observatory (Kenney et al. 1993) shows a central starburst localized within 4'' (628 pc) and peaks in a ring of 1''–2'' (157–314 pc). A central starburst is clearly seen in Hα imaging (Kenney et al. 1993), with a star formation rate ∼2.3 ± 0.4 M_⊙ yr⁻¹ within the central 8 farcs 2 (990 pc) of the galaxy (D. Nguyen et al. 2020, in preparation). Moreover, the Palomar optical survey defines the nucleus of NGC 3504 as a transition object (Ho et al. 1993, 1997) and suggests the central starburst component probably powered by hot O-type stars (Ho & Filippenko 1993). NGC 3504 shows a compact nuclear unresolved source observed by very long baseline interferometry with a luminosity of ${L}_{1.4\mathrm{GHz}}=0.6\times {10}^{21}$ W Hz⁻¹ sr⁻¹ at 1.4 GHz (Condon et al. 1998; Deller & Middelberg 2014).

3. Data and Data Reduction

3.1. HST Imaging

We use HST observations in the WFC3 IR F160W band to create a mass model (Section 5), which will be used as an input ingredient for dynamical models in Sections 5 and 6. The HST data were observed on 2012 May 1 (GO-12450, PI: Kochanek) with a total exposure time of 1398 s. We download these images from the Hubble Legacy Archive (HLA²¹ ) directly and use them throughout our analysis. However, to test the sky background level accurately, we also download the flat-fielded (flt) images from the HST/Barbara A. Mikulski Archive for Space Telescopes, combine these images in the same filters without sky subtraction using drizzlepac/Astrodrizzle (Avila & Hack 2012), and then add the differences back into the images. Finally, we compare the combined image with the HLA image. We choose the F160W filter as the default image to model the mass-follows-light map, while the F110W image is an alternative mass model to estimate the uncertainty caused by different wavelengths and extinction (Section 5.3.1).

Figure 2 shows the galaxy NGC 3504 seen in the F160W image with bars connecting the bulge and galactic disk on the left and its zoom-in on the right. There are a few dust lanes that circle the galaxy center, suggesting a CND extends to a radius of at least 4''.

**Figure 2.** Left panel: *HST*/WFC3 IR F160W image of NGC 3504 within the field of view (FOV) of 86'' × 86'' (13.5 kpc × 13.5 kpc). Right panel: the zoom-in for the nucleus region within the FOV of $9\buildrel{\prime\prime}\over{.} 7\times 9\buildrel{\prime\prime}\over{.} 7$ (1.5 kpc × 1.5 kpc) overlaid with black ${}^{12}\mathrm{CO}(2-1)$ integrated intensity contours from our ALMA observation from the low-resolution measurement set.
Download figure:
Standard image High-resolution image

The photocenter of the HLA images is at (R.A., decl.) = ( ${11}^{{\rm{h}}}{03}^{{\rm{m}}}{11}^{{\rm{s}}}.210,+27^\circ {58}^{{\prime} }21\buildrel{\prime\prime}\over{.} 00$ ) in the J2000-equatorial coordinate system as presented in HLA images. This is offset $\sim 0\buildrel{\prime\prime}\over{.} 1$ compared with the peak of the compact ${}^{12}\mathrm{CO}(2-1)$ continuum emission and the optical center of the galaxy derived from optical images (Section 3.2.1). We therefore assume the continuum emission is aligned with the photocenter because the size of this offset is comparable to the astrometry errors of the HLA data.

We use Tiny Tim PSFs (Jedrzejewski 1987; Jedrzejewski et al. 1987) for the WFC3 IR F160W and F110W images to decompose the multi-Gaussian expansion (MGE) model (Emsellem et al. 1994; Cappellari 2002) in Section 4.

3.2. ALMA Observations

The observations of the ¹²CO(2 − 1) line in the nucleus of NGC 3504 are carried out with ALMA as a part of the Weighing Black Hole Masses in Low-mass Galaxies Project (Program 2017.1.00964.S, PI: Nguyen, Dieu). The correlators cover the data in four spectral windows (SPWs) including one 1875 MHz frequency division mode (FDM) SPW that covers over the ${}^{12}\mathrm{CO}(2-1)$ line and three 2 GHz time division mode (TDM) SPWs added simultaneously to detect continuum emission. The observations use 46 of ALMA's 12 m antennas with the C43-9 and C43-6 configurations so that the data have a maximum recoverable scale of 8 farcs 2 in diameter. The raw ALMA data are calibrated by the ALMA regional center staff using the standard ALMA pipeline. Flux and bandpass calibration is conducted using the quasar J1058+0133, while the atmospheric phase offsets of the data are determined by a phase calibrator using J1102+2757. We summarize more details of these observations in Table 1.

Table 1. ALMA Observation Parameters

Phase center	R.A.	Decl.
	${11}^{{\rm{h}}}{03}^{{\rm{m}}}{11}^{{\rm{s}}}.205$	$+27^\circ {58}^{{\prime} }20\buildrel{\prime\prime}\over{.} 80$
	High Res.	Low Res.
Configurations	C43-9	C43-6

Obs. Date	2017 Oct 24	2018 Jan 1
Exposure time	46.6 minutes	25.7 minutes
Beam size	0042 × 0030	0221 × 0164
	6.6 pc × 4.7 pc	34.5 pc × 25.7 pc
Beam PA	−2°.0	32°.3

SPWs	FDM	TDM

Velocity res.	1.5 km s⁻¹	39.5 km s⁻¹
Frequency res.	1.13 MHz	31.25 MHz

Download table as: ASCII Typeset image

Continuum emission is detected at the center of NGC 3504 only (top left panel, Figure 3) and measured over the full line-free bandwidth, fitted by a power-law function (Section 3.2.1) and then subtracted from the data in the uv-plane using the task uvcontsub of the Common Astronomy Software Applications (CASA; McMullin et al. 2007) package version 5.1.1.

**Figure 3.** The top left panel shows the zoom-in map of the 1.3 mm continuum emission at the center of NGC 3504. The other three panels are the moment maps of the detected ${}^{12}\mathrm{CO}(2-1)$ emission created using the masked moment technique described in Section 3.2.2 including the zeroth (integrated intensity, top right panel), the first (intensity-weighted mean velocity field, bottom left panel), and the second moment (intensity-weighted velocity dispersion, bottom right panel) maps. The integrated intensity map shows a CND attenuation hole, which has a diameter of ∼12.6 pc (008) and is centered at the location of the putative SMBH. The ellipse at the bottom left of each panel shows the synthesized beam of the image.
Download figure:
Standard image High-resolution image

We first combine the visibility files of these two measurement sets into a final continuum-free calibrated data set using the CASA task concat with the optimal combination ratios of 2/3 and 1/3 for the high– and low–spatial resolution measurement set during the visweightscale mode, respectively. Second, we create a three-dimensional (R.A., decl., velocity) combined cube from the continuum-free calibrated file using the clean task. To model the CND, estimate the mass of compact central objects, and optimize the sensitivity of diffuse gas and resolution, we image the combined cube using the channel width of 10 km s⁻¹ (Davis 2014), Briggs weighting with a robust parameter of 0.5, and pixel size of 0 farcs 013. We reduce the side lobes on the image by using a mask during the clean interactive mode. We estimate the rms noise, rms = 42 μJy beam⁻¹, in a few signal-free channels, then set 3× rms as the clean threshold in regions of source emission for dirty channels. We also do primary beam correction after cleaning. The final calibrated ${}^{12}\mathrm{CO}(2-1)$ cube has a synthesized beam of 0 farcs 042 × 0 farcs 030 (6.6 pc × 4.7 pc), PA ∼ 9 fdg 7, and rms ∼0.187 mJy beam⁻¹ km s⁻¹. The nuclear molecular gas emission of NGC 3504 is detected from 1340 to 1680 km s⁻¹ with the systemic velocity of 1525 km s⁻¹.

3.2.1. Continuum Emission

The continuum emission is clearly resolved and centrally peaked as seen in the top left panel of Figure 3. We identify the continuum peak as the galaxy/kinematic center. To determine the position, size, and total integrated intensity of this source, we fit this continuum profile with a Gaussian using the CASA task imfit. The emission is centered at (R.A., decl.) = ( ${11}^{{\rm{h}}}{03}^{{\rm{m}}}{11}^{{\rm{s}}}.115\pm {0}^{{\rm{s}}}.040$ , $+27^\circ {58}^{{\prime} }20\buildrel{\prime\prime}\over{.} 80\pm 0\buildrel{\prime\prime}\over{.} 18$ ) with the error bars including both of the imfit fit and ALMA astrometric uncertainties. The SDSS data release 14 (Abolfathi et al. 2018) also reports an optical center at ( ${11}^{{\rm{h}}}{03}^{{\rm{m}}}{11}^{{\rm{s}}}.112$ , $+27^\circ {58}^{{\prime} }20\buildrel{\prime\prime}\over{.} 77$ ), and imfit estimates the size of the emission source of $(0\buildrel{\prime\prime}\over{.} 011\pm 0\buildrel{\prime\prime}\over{.} 004)\times (0\buildrel{\prime\prime}\over{.} 011\pm 0\buildrel{\prime\prime}\over{.} 004)$ or (1.73 ± 0.63) pc × (1.73 ± 0.63) pc, with a PA of ∼30° and the total integrated intensity of 4.70 ± 0.16 (random error) ± 0.24 (systematic from flux calibration) mJy over the emissions free of upper sideband (242–246 GHz) and lower sideband (227.5–231.5 GHz) frequency windows.

3.2.2. The ¹²CO(2–1) and CS(5–4) Emission Lines

We create the integrated intensity, intensity-weighted mean velocity field, and intensity-weighted velocity dispersion maps for the nuclear ${}^{12}\mathrm{CO}(2-1)$ gas disk from the combined cube directly in Figure 3. The emission is significant within the radius of 4''. We enhance the quality of these maps using the moments masking technique (Dame 2011; Davis et al. 2017, 2018; Onishi et al. 2017; Nagai et al. 2019; North et al. 2019; Smith et al. 2019).

The velocity field map reveals the presence of a warped, nuclear rotating disk with a total velocity width of ∼300 km s⁻¹. The velocity dispersion map is quite flat, with constant values of ∼39 and ∼25 km s⁻¹ at $0\buildrel{\prime\prime}\over{.} 1\lt r\lesssim 0\buildrel{\prime\prime}\over{.} 3$ and $0\buildrel{\prime\prime}\over{.} 3\lt r\lt 0\buildrel{\prime\prime}\over{.} 5$ , respectively. However, there is a suddenly increasing peak of ∼65 km s⁻¹ at the galactic center ( $r\lesssim 0\buildrel{\prime\prime}\over{.} 1$ ). We interpret this not as an intrinsic velocity dispersion but rather as a beam smearing caused by the nuclear velocity gradient and LOS integration of an inclined CND (i ≳ 50°).

The intensity map shows a central attenuated hole, which is marginally resolved on parsec scales as seen in the ${}^{12}\mathrm{CO}(2-1)$ integrated intensity map (top right panel, Figure 3). This hole could indicate there is no flux present at the center and the observed flux is due to foreground emission or due to the extent of the synthesized beam. Alternatively, the CO gas could be present near within the hole, with the lack of the ${}^{12}\mathrm{CO}(2-1)$ emission being due to changing excitation conditions.

Figure 4 shows the integrated ${}^{12}\mathrm{CO}(2-1)$ spectrum of NGC 3504 that has a classical double-horn shape of a rotating disk, while in Figure 5, we plot the position–velocity diagram (PVD) extracting from a cut through the major axis of the gas disk (PA ∼ 335°). We show both PVDs of the combined cube (high–spatial resolution) and the low–spatial resolution cube to demonstrate the increased rotation toward the center at different spatial scales. We interpret this rotation so far as a Keplerian curve caused either by an SMBH existing at the heart of NGC 3504 or by noncircular motions of an inner ring/spiral within the innermost regions, which we usually see in barred galaxies.

**Figure 5.** Left panel: the PVD of the ${}^{12}\mathrm{CO}(2-1)$ emission of the combined cube extracted along the major axis with a slit of five pixels in width (0065 or 10.2 pc). We do not show the beams of 0042 × 0030 (6.6 pc × 4.7 pc) and velocity channel width of 10 km s⁻¹ explicitly in this plot due to their small size. Right panel: The same PVD for the low-resolution cube with the beam size of 0221 × 0164 (34.5 pc × 25.7 pc).
Download figure:
Standard image High-resolution image

We also detect a dense gas tracer $\mathrm{CS}(5-4)$ in one of the continuum SPWs, which is centrally peaked and fills in the attenuated hole of the ${}^{12}\mathrm{CO}(2-1)$ map (see the top right panel of Figure 3 and the left panel of Figure 6). In Figure 6, we show the integrated intensity map, spectrum, and radial profile of the $\mathrm{CS}(5-4)$ line in the left, middle, and right panels, respectively. The detection of $\mathrm{CS}(5-4)$ is significant above 20σ (σ ∼ 0.3 mJy beam⁻¹ km s⁻¹) but in a very low-velocity-resolution spectral window with only 128 channels over 2 GHz bandwidth (∼40 km s⁻¹). We estimate its total flux to be 1.76 ± 0.42 Jy km s⁻¹ with 10% of the error budget coming from the flux calibration uncertainty of ALMA data. The line also follows the same kinematic signatures of the ${}^{12}\mathrm{CO}(2-1)$ emission as shown at the moment 1, moment 2, and PVD maps in Figure 7. These features suggest the $\mathrm{CS}(5-4)$ line is an alternative transition that could provide a better constraint on the mass of a central dark object in NGC 3504 than the ${}^{12}\mathrm{CO}(2-1)$ line. This is because the central peak of the $\mathrm{CS}(5-4)$ line would recover the high-velocity upturn in the data at the very central region, much further within the accretion disk and SOI that is somewhat missing in the current ${}^{12}\mathrm{CO}(2-1)$ map due to the centrally attenuated hole. However, the velocity resolution of the $\mathrm{CS}(5-4)$ emission line is too coarse to perform such dynamical models. A higher-velocity-resolution observation for the $\mathrm{CS}(5-4)$ line with at least ∼4× improvement in velocity resolution is required to reduce the uncertainty of the central dark object's mass determination.

**Figure 6.** The dense gas tracer $\mathrm{CS}(5-4)$ detected in one of the continuum spectral windows (left), concentrated in the region of the ${}^{12}\mathrm{CO}(2-1)$ attenuated hole. The integrated spectrum (middle) and radial profile (right) of $\mathrm{CS}(5-4)$ show that it is centrally filled rather than a hole.
Download figure:
Standard image High-resolution image

**Figure 7.** The moment maps of the detected $\mathrm{CS}(5-4)$ emission created using the masked moment technique described in Section 3.2.2 including the first (intensity-weighted mean velocity field, left panel) and the second moment (intensity-weighted velocity dispersion, middle panel). Right panel: the PVD of the $\mathrm{CS}(5-4)$ emission line, which has the same observed velocity and velocity width of the ${}^{12}\mathrm{CO}(2-1)$ emission line, as well as the rotational signature around the galaxy center. These maps have a low-velocity resolution of ∼40 km s⁻¹.
Download figure:
Standard image High-resolution image

The physical radius of the hole observed in the integrated intensity map is estimated to be $\sim 0\buildrel{\prime\prime}\over{.} 04$ (or 4.8 pc) and seems to be larger than would be expected for a torus (e.g., Rieke & Low 1972; Antonucci 1993; Tristram & Schartmann 2011), which is typical of <1 pc (Barvainis 1987). The reason is that the ${}^{12}\mathrm{CO}(2-1)$ emission is sometimes deficient in the nuclear region and therefore may not be the best tracer for the torus (Imanishi et al. 2018; Izumi et al. 2018). To this end, this central deficit emission of ${}^{12}\mathrm{CO}(2-1)$ is also found in other studies by Davis et al. (2017), North et al. (2019), Smith et al. (2019), and T. Davis et al. (2020, in preparation), suggesting that the prevalence of such central holes may be caused by changing excitation conditions rather than the true absence of molecular gas there.

4. Creating a Mass Model

In galaxies, the stellar mass-to-light (M/L) ratio can vary due to the presence of complex nuclear stellar populations (Seth et al. 2010; McConnell et al. 2013; Ahn et al. 2017; Nguyen et al. 2017; Ahn et al. 2018; N18; N19). This variation causes large uncertainties in dynamical ${M}_{\mathrm{BH}}$ estimates (N17; N19). Unlike our work in N17 and N19, we do not have the nuclear stellar spectroscopic information for NGC 3504. However, the examination of the nuclear color map assuming HST/WFC3 F110W $\approx \,J$ and F160W ≈ H finds a constant color of $J-H\sim 0.45$ mag across the FOV where the ${}^{12}\mathrm{CO}(2-1)$ kinematics were measured, although there is some variability in the central dust lanes and twice the color at the large scale of the galaxy where $J-H\sim 0.257$ . Photometry taken with the Jacobus Kapteyn Telescope images also finds a constant $B-I\sim 1.3$ within the FOV of 5'' (Knapen et al. 2002). These uniform colors in the region of interest suggest that it is acceptable to make the simple assumption of using a constant M/L to create the mass model for NGC 3504. We utilize an MGE technique (Emsellem et al. 1994; Cappellari 2002) with the fitting method and software²² mge_fit_sectors IDL version 4.14 by Cappellari (2002) to decompose the photometric surface brightness density of the F160W image and deconvolve the effects of a PSF into individuals Gaussian components.

We first parameterize the PSF (Section 3.1) using Gaussian functions in the first MGE fit and then use them as an input during the second MGE fit to obtain a deconvolved MGE model of the galaxy. This PSF MGE model is tabulated in Table 2. During the second MGE fit, we supply a mask that removes pixels that are contaminated by prominent dust lanes, bright stars, and a strongly asymmetric stellar component seen at 2'' toward the south of the nucleus in both the visible and IR.²³ This asymmetric stellar component and the bar located further out at ∼25'' cause a twist and edge-on orientation in the MGE model if the PA and axis ratio (b/a) are left as free parameters in the fit. Thus, we fix the PA = 335° based on the orientation of the galactic disk for a default axisymmetric mass MGE model and set an optimal axis ratio $b/a=0.65\mbox{--}1$ to optimize data versus the MGE model, which is shown in Figure 8. This model generally fits the data worse than the nonaxisymmetric one at the radius $r\gt 25^{\prime\prime}$ , where it is well outside of the ALMA-detected nuclear ${}^{12}\mathrm{CO}(2-1)$ gas and completely irrelevant for our modeling results. However, it is the best representation of the data in the region of $r\lesssim 10^{\prime\prime}$ , where the MGE model impacts the dynamical results dramatically. A similar axisymmetric model based on the F110W imaging provides a check on systematic errors from our mass models. Note that the incorporation of gas mass is discussed in Section 5.3.1.

**Figure 8.** Panel (a): The comparison between the *HST*/WFC3 F160W photometry of NGC 3504 (open squares) and its corresponding best-fit MGE models (red solid line), which are summed up from multiple Gaussians (color thin lines). We show the best-fit model that is projected along with a 5°-wide sector, which has an angle of 30° to the major axis. Panel (b): The fractional residuals `(Data-Model)/Data` illustrate the agreement between data and model within 15% across 100'' (15.7 kpc). Panels (c) and (d): The comparison between the F160W light distribution and its best-fit MGE model is illustrated in the form of 2D surface brightness density contours at the whole galactic scale and the zoom-in of the $30^{\prime\prime} \times 30^{\prime\prime}$ FOV. Black contours show the data, while the red contours show the model to highlight the agreement between data and model at the same radii and contour levels.
Download figure:
Standard image High-resolution image

$30^{\prime\prime} \times 30^{\prime\prime} $ — **Figure 8.** Panel (a): The comparison between the *HST*/WFC3 F160W photometry of NGC 3504 (open squares) and its corresponding best-fit MGE models (red solid line), which are summed up from multiple Gaussians (color thin lines). We show the best-fit model that is projected along with a 5°-wide sector, which has an angle of 30° to the major axis. Panel (b): The fractional residuals `(Data-Model)/Data` illustrate the agreement between data and model within 15% across 100'' (15.7 kpc). Panels (c) and (d): The comparison between the F160W light distribution and its best-fit MGE model is illustrated in the form of 2D surface brightness density contours at the whole galactic scale and the zoom-in of the $30^{\prime\prime} \times 30^{\prime\prime}$ FOV. Black contours show the data, while the red contours show the model to highlight the agreement between data and model at the same radii and contour levels.
Download figure:
Standard image High-resolution image

Table 2. MGE Parameters of the HST/WFC3 IR F160W PSF

j	Total Count	σ
	of Gaussian_j	(arcsec)
(1)	(2)	(3)
1	0.343	0.049
2	0.538	0.130
3	0.060	0.403
4	0.033	0.897
5	0.033	1.716

Note. Column 1: Gaussian component number. Column 2: The MGE model represents for the total light of each Gaussian. Column 3: The Gaussian dispersion along the major axis.

Download table as: ASCII Typeset image

The observed MGE parameters are listed in Table 3, assuming the HST/WFC3 F160W Solar Vega absolute magnitude system.²⁴ Each Gaussian can be deprojected analytically with a specific axis ratio (or i) to reconstruct a three-dimensional light distribution model.

Table 3. The HST/WFC3 IR 160W MGE Model of NGC 3504

j	(Luminosity Density)	σ	b/a
	( ${L}_{\odot }/{\mathrm{pc}}^{2}$ )	(arcsec)
(1)	(2)	(3)	(4)
1	398848.849	0.081	1.00
2	24716.727	0.379	1.00
3	170.703	0.561	0.69
4	626.707	1.026	0.69
5	2663.322	1.571	0.69
6	13091.046	1.982	0.69
7	15532.957	2.168	0.69
8	643.406	0.277	0.67
9	4102.848	5.773	0.67
10	518.000	9.147	0.65
11	208.449	23.588	0.65

Note. MGE models using KinMS and tilted-ring model fits in Sections 5 and 6. Column 1: Gaussian component number. Column 2: The MGE model represents for the galaxy luminosity model. Column 3: The Gaussian dispersion along the major axis. Column 4: The axial ratios.

Download table as: ASCII Typeset image

5. KinMS Dynamical Modeling

This section describes the KinMS²⁵ dynamical model (Davis 2014) that we employ to measure the central dark mass in NGC 3504 and present our modeling results.

The KinMS model is a millimeter-wave observational simulation tool developed by Davis et al. (2013) that uses the Markov Chain Monte Carlo (MCMC) method to simulate the dynamical motion of molecular/atomic cold gas distributions under the gravity of all material within a galaxy (gas/dust, stars, BH, etc.). First, an initial guess is made, and the gas distribution and kinematics are calculated assuming that the gas rotates on circular orbits governed by the gravitational potentials using the mge_circular_velocity procedure within the IDL Jeans anisotropic modeling (Cappellari 2008, see footnote 5) package. KinMS then creates a simulated cube that is compared to the observed data to compute the likelihood function (Davis et al. 2017, 2018; Onishi et al. 2017; North et al. 2019; Smith et al. 2019). We use Bayesian analysis with a set of walkers through the emcee algorithm (Foreman-Mackey et al. 2013) and affine-invariant ensemble sampler (Goodman & Weare 2010) to explore the parameter space. The relative likelihood for each walker at each step will determine the next move through parameter space. From the posterior distribution of the full pool of model parameters calculated at each step, we obtain the best-fit model parameters. In this work, we use the python code KINMSpy_MCMC²⁶ to find the best set of model parameters.

The KinMS model allows us to input a radial function, Σ(r), to describe the nuclear molecular gas distribution. We assume an axisymmetric morphology for the ${}^{12}\mathrm{CO}(2-1)$ gas in NGC 3504 distributed continuously out to the radius of ∼4 farcs 1 (644 pc).

Although an attenuated central hole is observed in the combined cube, it is filled in with denser gas as seen in Figure 6, suggesting the total distribution of the nuclear gas does not have such a hole. Therefore, we extract the flux from the low–spatial resolution cube (where no hole is seen) in a slit of 7 pixel width and ∼0 farcs 7 radius along the major axis in Figure 9. We model the CND morphology using a single-Gaussian function with only one free parameter of dispersion (G_σ). For the low-resolution cube, we tie the Gaussian center to the kinematic center. However, for the high-resolution cube, due to the central attenuated/deficit hole, we fix the Gaussian center off the kinematic center by ∼0 farcs 11 (17.3 pc). We also incorporate the normalization factor into the total flux parameter (f); this factor scales both the observed gas intensity and the gas mass. We note that the KinMS model matches the observations by fitting a set of free parameters including G_σ, f, i, kinematic center in R.A. (x_c) and decl. (y_c), velocity offset ( ${v}_{\mathrm{off}.}$ ), M/L, M_BH, gas disk thickness (d_t, see Appendix A), and turbulent velocity dispersion of the gas (σ₀), which is assumed to be a spatially constant parameter. A clear kinematic warp is seen in Figure 3. We characterize this warp using a radial PA profile extracted along the major axis with the Kinemetry²⁷ code (Krajnović et al. 2006), which varies by ∼15° with a transition at ∼0 farcs 3. This does not add any free parameters. Our final KinMS model has 10 free parameters; these are listed in Table 4.

**Figure 9.** The nuclear gas morphological distribution of the low–spatial resolution cube is shown with a cut along the major axis through the center of NGC 3504 integrated intensity map $(1\buildrel{\prime\prime}\over{.} 4\times 1\buildrel{\prime\prime}\over{.} 4)$ or (220 pc × 220 pc). The data are plotted in black, while our best-fit KinMS model for our chosen surface brightness profile is plotted in red.
Download figure:
Standard image High-resolution image

**Figure 9.** The nuclear gas morphological distribution of the low–spatial resolution cube is shown with a cut along the major axis through the center of NGC 3504 integrated intensity map $(1\buildrel{\prime\prime}\over{.} 4\times 1\buildrel{\prime\prime}\over{.} 4)$ or (220 pc × 220 pc). The data are plotted in black, while our best-fit KinMS model for our chosen surface brightness profile is plotted in red.
Download figure:
Standard image High-resolution image

Table 4. Best-fit KinMS Model Parameters and Statistical Uncertainties for the Combined ALMA Data Cube (High Resolution) within the FOV of $0\buildrel{\prime\prime}\over{.} 8\times 0\buildrel{\prime\prime}\over{.} 8$

Parameter (Units)	Search Range			Best Fit	1σ Error³	3σ Error^a	3σ Error^b	3σ Error^c
		Uniform			(16%–84%)	(0.14%–99.86%)	(0.14%–99.86%)	(0.14%–99.86%)
(1)		(2)		(3)	(4)	(5)	(6)	(7)
Black Hole:
${\mathrm{log}}_{10}{M}_{\mathrm{BH}}$ (M_⊙)	4.00	$\longrightarrow$	9.00	7.19	−0.05, +0.05	−0.15, +0.15	−0.07, +0.07	−0.12, +0.12
M/L_H (M_⊙/ ${L}_{\odot }$ )	0.01	$\longrightarrow$	3.00	0.44	−0.04, +0.04	−0.12, +0.12	−0.06, +0.06	−0.10, +0.10
f (Jy km s⁻¹)	10.0	$\longrightarrow$	500.0	76.94	−1.33, +1.34	−4.00, +4.02	−1.66, +1.68	−3.21, +3.26
Molecular Gas:
i (°)	45.0	$\longrightarrow$	90.0	53.38	−1.39, +1.22	−4.20, +3.68	−1.65, +1.66	−3.03, +2.95
${\sigma }_{0}$ (km s⁻¹)	1.0	$\longrightarrow$	80.0	38.56	−1.02, +1.03	−3.08, +3.10	−1.54, +1.54	−2.65, +2.68
d_t ('')	0.01	$\longrightarrow$	1.50	0.08	−0.01, +0.01	−0.03, +0.03	−0.00, +0.00	−0.01, +0.01
G ${}_{\sigma }$ ('')	0.01	$\longrightarrow$	1.00	0.14	−0.01, +0.01	−0.03, +0.03	−0.01, +0.01	−0.02, +0.02
Nuisance:
x_c ('')	−0.10	$\longrightarrow$	+0.10	−0.01	−0.00, +0.00	−0.01, +0.01	−0.00, +0.00	−0.01, +0.01
y_c ('')	−0.10	$\longrightarrow$	+0.10	−0.01	−0.00, +0.00	−0.01, +0.01	−0.00, +0.00	−0.01,+0.01
${v}_{\mathrm{off}.}$ (km s⁻¹)	−10.0	$\longrightarrow$	+10.0	5.85	−0.27, +0.31	−0.82, +0.93	−0.31, +0.32	−0.60, +0.64

Note. Column 1: a list of the fitted model parameters. Column 2: a list of the priors of the fitted model parameters and their search ranges. The prior is constructed in the uniform linear space except for the SMBH, which is in logarithmic space. Columns 3–5: the best-fit value of each parameter and their uncertainties at the 1σ and 3σ confidence levels (CLs). The R.A., decl., and velocity offset nuisance parameters are defined relative to the ALMA data phase center position ( ${11}^{{\rm{h}}}{03}^{{\rm{m}}}{11}^{{\rm{s}}}.205$ , $+27^\circ {58}^{{\prime} }20\buildrel{\prime\prime}\over{.} 80$ , ${V}_{\mathrm{sys}.}=1525$ km s⁻¹).^aUncertainty constrained by increasing the rms cube by a factor of ${(2N)}^{0.25}$ as presented in Mitzkus et al. (2017) and Smith et al. (2019); we adopt this approach for our canonical uncertainties estimates in this work. ^bUncertainty constrained from the covariance matrix. ^cUncertainty constrained from the binned data cube where the spatial scale corresponds roughly to 1 pixel per synthesized beam. See detailed discussions of these uncertainty constraints in Appendix B.

Download table as: ASCII Typeset image

5.1. Fitting Process

We run the first KinMS model in an area of 160 pixels × 160 pixels (2'' × 2''), which covers most of the features of the molecular gas disk. From this run, we obtain rough estimates of the model parameters. The second KinMS fit then fits in the central 64 pixels × 64 pixels ( $0\buildrel{\prime\prime}\over{.} 8\times 0\buildrel{\prime\prime}\over{.} 8$ or 126 pc × 126 pc) area and 40 velocity channels (10 km s⁻¹ for each channel). This fit starts with flat priors over a reasonable range (see column 2 of Table 4) to ensure our kinematic fitting process converges; these priors are determined to be the best fit from the first fit. We note that the prior on the ${M}_{\mathrm{BH}}$ is flat in log-space, while the inclination of the gas disk varies over the full physical range allowed by the MGE model. The best-fit models are always found well within the range of the priors.

We next include the mass of the interstellar material (gas and dust) within the fitting region into our mass model. The total flux in this region is 79 ± 2 (random error) ±8 (systemic from flux calibration) Jy km s⁻¹, giving ${M}_{{{\rm{H}}}_{2}}=(6.5\pm 2.6)\times {10}^{8}$ M_⊙ by assuming the line ratio ${}^{12}\mathrm{CO}(2-1){/}^{12}\mathrm{CO}(1-0)\,=0.8$ (Bigiel et al. 2008) and H₂-to-CO conversion factor for starburst galaxies: $N({{\rm{H}}}_{2})/{I}_{(1-0)}={X}_{\mathrm{CO}}=0.5\times {X}_{\mathrm{CO}}^{\mathrm{Milky}\ \mathrm{Way}}\,=(1.0\pm 0.3)\times {10}^{20}$ cm⁻² (K km s⁻¹)⁻¹ (Kuno et al. 2000, 2007; Bolatto et al. 2013). The dust mass of (8.2 ± 3.9) × 10⁶M_⊙ is calculated by D. Nguyen et al. (2020, in preparation) using the ALMA continuum (Section 3.2.1) and far-IR Herschel archival data. In Figure 10, we plot the 1D stellar-mass and gas + dust mass profiles as radial functions simultaneously within the radius of 1''. The gas + dust mass distributions are comparable to the stellar mass, suggesting the gravitational effects of these components play an important role in determining the ${M}_{\mathrm{BH}}$ in NGC 3504 accurately. We add these masses in the KinMS gasGrav mechanism, assuming the dust and gas are cospatially distributed.

**Figure 10.** One-dimensional cumulative mass profile of stellar mass (red) and gas + dust (green) components plotted within the radius of 10 from the center.
Download figure:
Standard image High-resolution image

farcs — **Figure 10.** One-dimensional cumulative mass profile of stellar mass (red) and gas + dust (green) components plotted within the radius of 10 from the center.
Download figure:
Standard image High-resolution image

Finally, we run the final model with an iteration of 3 × 10⁶, and the first 25% of iterations are considered as the burn-in phases to produce our final posterior probability distributions of 10 free parameters.

5.2. Model Results

We detect a compact dark mass that causes the Keplerian rise toward the center. The best-fit parameters and their uncertainties are identified directly from the Bayesian analysis, relying on the likelihood probability distribution functions (PDFs) generated via MCMC (see Appendix B). We choose the best fit of each parameter to be the median of its posterior PDF. This value is not much different (<7%) from the best fit identified from the ${\chi }_{\min }^{2}$ formalism.

In particular, the probability is marginalized over to produce the best-fit value as the median of the marginalized posterior for each parameter. The 1σ, 2σ, and 3σ CLs are estimated from all models within 16%–84%, 2.3%–97.7%, and 0.14%–99.86% of the PDFs, respectively. Although we accept the standard practice in almost astronomical literature to quote systematics at 1σ or 2σ errors, the ALMA noise covariance discussed in Appendix B tends to make our parameter uncertainties unphysically small. We, therefore, use the 3σ error to demonstrate the "actual" systemic errors. At 3σ CL, the KinMS model gives ${M}_{\mathrm{BH}}={1.6}_{-0.4}^{+0.6}\times {10}^{7}$ M_⊙, M/L ${}_{{\rm{H}}}\,=0.44\pm 0.12$ (M_⊙/ ${L}_{\odot }$ ), and $i=53.38^\circ {}_{-4.20}^{+3.68}$ . The best-fit model has ${\chi }_{\min }^{2}=193,659$ or ${\chi }_{\mathrm{red}}^{2}=1.182$ for $(64\times 64\times 40)\,-10=163,830$ degree of freedom (dof), indicating a quite good fit (Bevington 1969). The nonunity reduced ${\chi }_{\mathrm{red}}^{2}$ could be due to data model mismatch, for instance, the model being unable to replicate some observed structures in the CND or the assumption of constant turbulent gas velocity being a poor one. The best-fit parameters and their likelihoods are listed in Table 4.

Figure 11 shows the 1D and 2D marginalization of the physical parameters included in the fit. Covariances are visible in several panels. First, the covariance of ${M}_{\mathrm{BH}}$ versus $M/{L}_{{\rm{H}}}$ , which describes the degeneracy between the potentials of the SMBH and galaxy itself, happens when the observational scales (e.g., HST and ALMA) are not resolved deeply within the SMBH SOI. Second, the covariance of G_σ versus i implies the positive dependence of the CND's inclination angle on the assumed Gaussian of the nuclear molecular gas distribution ( ${\rm{\Sigma }}(r)$ ). Finally, the covariance of ${\sigma }_{0}$ versus f happens if f tends to be normalized to the surface brightness profile (Smith et al. 2019). These covariances are usually seen in a simultaneous constraint.

We show the observed PVD overlaid with the PVDs of the model without a BH, with the best-fit BH, and with an overmassive BH in Figure 12 to illustrate how well the model fits the data. The best-fit model with an appropriate ${M}_{\mathrm{BH}}$ clearly reproduces the kinematics of the molecular gas better than the other two. We show the observed and the best-fit model velocity map and the velocity-residual map (Data-Model) in the plots in the bottom row of Figure 12. The residuals are ∼10 km s⁻¹ except for some noncircular motions that resulted in high residuals on the south/southwest side (∼24 km s⁻¹) along the minor axis of the nucleus; these could be signatures of outflows, inflows, or other noncircular motions induced by a bar. In Appendix C, we present the channel maps from the best-fit model overlaid on the channel maps of the data. We also demonstrate the same results for the low–spatial resolution data. The full list of the best-fit parameters and their likelihoods is also presented in Appendix C.

The best-fit ${M}_{\mathrm{BH}}={1.6}_{-0.4}^{+0.6}\times {10}^{7}$ M_⊙ with the stellar bulge velocity dispersion $\sigma =119.3\pm 10.3$ km s⁻¹ (Ho et al. 2009) and the BH in NGC 3504 has an intrinsic SOI radius ${R}_{\mathrm{SOI}}={{GM}}_{\mathrm{BH}}/{\sigma }^{2}\sim 5$ pc ( $\sim 0\buildrel{\prime\prime}\over{.} 032;$ Merritt 2013); we thus resolved the SOI of the SMBH by a factor of ∼1.5 larger than the resolvable scale of our ALMA observations ( $\sim 0\buildrel{\prime\prime}\over{.} 042$ ). For a discussion of the reliability of this measurement in the context of our observational scales, see Section 7.2.

The ${M}_{\mathrm{BH}}$ and M/L_H estimated from the combined and low–spatial resolution cubes are both consistent with each other within 22% and 2%, respectively (see Tables 4 and 6) even at 5× the difference of the observational scales.

5.3. Uncertainty Sources on the ${M}_{\mathrm{BH}}$ estimate

5.3.1. Stellar-Mass Models

The uncertainties of our analysis so far are based on the ALMA kinematic errors only. Here we examine the mass model uncertainty by analyzing independent models derived under the various assumptions:

(1) Mass model from different photometric filter: We test the dependence on the stellar-mass model by using photometry from a different band to build the MGE mass model. We use the HST/WFC3 F110W image (approximate to the J band), which gives the best-fit ${M}_{\mathrm{BH}}$ of ${1.3}_{-0.6}^{+0.8}\times {10}^{7}$ M_⊙ and $M/{L}_{J}={0.47}_{-0.10}^{+0.10}$ (M_⊙/ ${L}_{\odot }$ ).

(2) Mass models relative to large galaxy structure: The way we fit the light MGE model from HST images also causes some errors on the mass model. Both the large and smaller scale structures of NGC 3504 host outer Lindblad resonances (OLRs; Buta & Crocker 1992; Buta & Combes 1996), especially for the central elongated bars at a distance of ≲25'' from the center (Kuno et al. 2000), suggesting nonaxisymmetric structures. We, therefore, model two axisymmetric MGE fits as follows: (i) there is no constraint on the axis ratio q to get an MGE with elongated barred Gaussian dominant, and (ii) set $q=0.9-1$ to get an MGE with OLR Gaussian dominant. We have double checked these MGEs within the $30^{\prime\prime} \times 30^{\prime\prime}$ central region, and they are not significantly different from each other. The KinMS model with the former MGE gives the best fit ( ${M}_{\mathrm{BH}}$ , ${\text{}}M/{L}_{H}$ , i) = ( ${2.3}_{-0.8}^{+0.7}\times {10}^{7}$ M_⊙, ${0.36}_{-0.12}^{+0.11}$ (M_⊙/ ${L}_{\odot }$ ), ${{82.54}_{-3.32}^{^\circ }}^{+3.54}$ ), while the KinMS model of the latter MGE derives to ( ${1.2}_{-0.8}^{+0.8}\times {10}^{7}$ M_⊙, ${0.48}_{-0.11}^{+0.10}$ (M_⊙/ ${L}_{\odot }$ ), ${{44.05}_{-4.61}^{^\circ }}^{+4.46}$ ). The former stellar-mass model tends to constrain an edge-on morphology for the molecular gas and then increases the ${M}_{\mathrm{BH}}$ due to adding high eccentric Gaussian components that match the bar.

(3) Mass model explicitly includes the stellar asymmetric component at 2'' south of the nucleus: So far, our stellar-mass model excludes this strongly asymmetric stellar distribution via the masking algorithm (Section 4). In this test, we include these pixels in the fit and adapt the twist-MGE model (asymmetric) that allows the PA of each Gaussian to change freely. Other setups are kept similar to the default axisymmetric MGE model. Although the twist-MGE model generates a somewhat better representation of the data in this asymmetric stellar distributed region, some discrepancies are seen in the opposite north (see Figure 13). The dynamical model of this twist-MGE model gives ${M}_{\mathrm{BH}}={1.8}_{-0.7}^{+0.9}\times {10}^{7}$ M_⊙ and ${\text{}}M/{L}_{H}={0.42}_{-0.12}^{+0.10}$ (M_⊙/ ${L}_{\odot }$ ), suggesting the stellar mass of this asymmetric component south of the nucleus is not significantly affected by the inferred ${M}_{\mathrm{BH}}$ . This proves our simplified assumption of axisymmetry in both stellar-mass distribution and dynamical models.

**Figure 13.** The comparison between the F160W light distribution and its best-fit twist-MGE model includes the asymmetric stellar distribution at 2'' south of the nucleus in the form of 2D surface brightness density contours within the central region of $30^{\prime\prime} \times 30^{\prime\prime}$ . Other illustrations are the same as panel (d) in Figure 8.
Download figure:
Standard image High-resolution image

We thus conclude that our ${M}_{\mathrm{BH}}$ and ${\text{}}M/{L}_{H}$ estimates are robust to the systematic errors of our mass models and that the nonaxisymmetric mass component contributed by the bar and the stellar asymmetric component at 2'' south of the nucleus, which are seen in both optical and IR, is not dynamically significant.

5.3.2. Distances

The ${M}_{\mathrm{BH}}$ estimate is systematically affected by the assumed distance to the galaxy based on the relation ${M}_{\mathrm{BH}}\propto D$ . We assume $D=32.4\pm 2.1$ Mpc for NGC 3504. Thus, the systemtic error caused by the distance uncertainty on ${M}_{\mathrm{BH}}$ is 6.5%, similar to the random uncertainties.

NGC 3504 has a wide range of existing distance estimates from 9 to 30 Mpc (de Vaucouleurs et al. 1981; Bottinelli et al. 1984; Tully & Fisher 1988; Russell 2002) after correcting for all the expected peculiar velocities mostly due to being behind Virgo cluster, but it seems that the Great Attractor pulls in the same direction as well. This large uncertainty of the distance adds ∼70% systematic to the ${M}_{\mathrm{BH}}$ , the largest fraction of the error budge.

5.3.3. The Galaxy Bar

The orientation of the bar along the LOS could affect the CND's circular motion dramatically (Randriamampandry et al. 2015). For instance, there are under-/overestimates of the circular motions when the bar is parallel/perpendicular to the projected major axis. As a SAB(s)ab galaxy, the rotation of NGC 3504 could be overestimated because the gas streaming motion along the bars compensates for the contribution of gas in the bulge as suggested by Figure 2 of Randriamampandry et al. (2015, top row plots). To calculate the effect of the bar on our dynamical results, we run the KinMS model with the velocity profile extracted along the major axis of the moment 1 map (Figure 3) using the IDL Kinemetry code (Krajnović et al. 2006) instead of the velocity model built from the mass MGE model in Section 5. The best-fit model gives ${M}_{\mathrm{BH}}$ of ${2.4}_{-0.8}^{+0.7}\times {10}^{7}$ M_⊙ and $M/{L}_{H}={0.41}_{-0.10}^{+0.11}$ (M_⊙/ ${L}_{\odot }$ ).

5.3.4. Turbulent Velocity Dispersion of the Gas

For the above analysis, we assumed a constant turbulent velocity dispersion for the gas, but in general, it could vary with radius and azimuth within the CND. Additionally, the increasing velocity dispersion at the galaxy center due to beam smearing could lead to overestimated ${M}_{\mathrm{BH}}$ . To quantify these effects, we allow a variable velocity dispersion as a function of radius. We test the turbulent velocity dispersion profile using the following prescriptions for $\sigma {(r)}_{\mathrm{gas}}$ :

(a) Linear gradient: $\sigma {(r)}_{\mathrm{gas}}=a\times r+b$ , where a and b are free parameters. We find $a\approx 0$ and b = 39.03 km s⁻¹, and the KinMS results are consistent with the default model of constant velocity dispersion.
(b) Exponential: $\sigma {(r)}_{\mathrm{gas}}={\sigma }_{0}\exp (-r/{r}_{0})+{\sigma }_{1}$ , where ${\sigma }_{0}$ , r₀, and ${\sigma }_{1}$ are free parameters. As discussed in Barth et al. (2016a), we set the lower boundary for $\sigma {(r)}_{\mathrm{gas},\min }=1$ km s⁻¹ during the fit to prevent the line-profile widths becoming arbitrarily small. The best-fit KinMS model provides ${M}_{\mathrm{BH}}={1.7}_{-0.6}^{+0.5}\,\times {10}^{7}$ M_⊙ and $M/{L}_{H}={0.42}_{-0.10}^{+0.12}$ (M_⊙/ ${L}_{\odot }$ ) with an exponential dispersion model with ${\sigma }_{0}\,={25.45}_{-3.07}^{+3.05}$ km s⁻¹, ${r}_{0}=0\buildrel{\prime\prime}\over{.} {18}_{-0.04}^{+0.05}$ , and ${\sigma }_{1}={33.73}_{-2.83}^{+2.98}$ km s⁻¹.
(c) Gaussian: $\sigma {(r)}_{\mathrm{gas}}={\sigma }_{0}\exp (-{(r-{r}_{0})}^{2}/2{\mu }^{2})+{\sigma }_{1}$ , where ${\sigma }_{0}$ , r₀, μ, and ${\sigma }_{1}$ are free parameters. We allow the parameter r₀ to vary over positive and negative values because the line width is sometimes offset from the center and also set the lower boundary for $\sigma {(r)}_{\mathrm{gas},\min }=1$ km s⁻¹ during the fit. The best-fit KinMS model provides ${M}_{\mathrm{BH}}={1.5}_{-0.4}^{+0.5}\,\times {10}^{7}$ M_⊙ and $M/{L}_{H}={0.43}_{-0.09}^{+0.09}$ (M_⊙/ ${L}_{\odot }$ ) with a Gaussian dispersion model with ${\sigma }_{0}={24.92}_{-2.80}^{+2.75}$ km s⁻¹, ${r}_{0}=0\buildrel{\prime\prime}\over{.} {04}_{-0.01}^{+0.01}$ , $\mu =0\buildrel{\prime\prime}\over{.} {15}_{-0.06}^{+0.06}$ , and ${\sigma }_{1}={32.89}_{-3.02}^{+2.98}$ km s⁻¹.

The linear, exponential, and Gaussian dispersion profile gives the minimum ${\chi }_{\mathrm{red}}^{2}$ at 1.208, 1.272, and 1.246, respectively. So, the assumption of constant ${\sigma }_{\mathrm{gas}}$ is good to describe the ${}^{12}\mathrm{CO}(2-1)$ CND's kinematics and to dynamically model ${M}_{\mathrm{BH}}$ .

5.3.5. Different Observational Scales and Distributed Assumption of the Gas

Because our ${}^{12}\mathrm{CO}(2-1)$ observations do not resolve the BH's SOI, the molecular mass estimates within the central synthesized beams, that is, $0\buildrel{\prime\prime}\over{.} 042\times 0\buildrel{\prime\prime}\over{.} 030$ for the combined cube and $0\buildrel{\prime\prime}\over{.} 221\times 0\buildrel{\prime\prime}\over{.} 164$ for the low–spatial resolution cube, are also other systemic mass uncertainties. We convert the total fluxes in the central beams of these cubes into molecular masses of $3.4\times {10}^{7}$ M_⊙ and $1.6\times {10}^{8}$ M_⊙, respectively. These masses' uncertainties are both much larger than the 3σ ${M}_{\mathrm{BH}}$ uncertainty provided by the KinMS model (Tables 4 and 6).

To test the impact of this gas in various assumptions of spatial distribution, we convert the ${}^{12}\mathrm{CO}(2-1)$ intensity map to the molecular gas map and then parameterize it into the MGE form in the same manner of the stellar-mass component, assuming the gas has an axisymmetric distribution (T. Davis et al. 2020, in preparation). In the KinMS fit, we turn off the gasGrav mechanism (which accounts the gas mass) but use the gas mass-MGE alternatively, giving ${M}_{\mathrm{BH}}={1.0}_{-0.8}^{+0.9}\,\times {10}^{7}$ M_⊙ and $M/{L}_{H}={0.47}_{-0.15}^{+0.12}$ (M_⊙/ ${L}_{\odot }$ ).

Our observations do not resolve the SOI that causes a large systemic uncertainty on the ${M}_{\mathrm{BH}}$ ; the ambiguous gas mass within the most central observational beam size is twice the BH mass. However, different assumptions of gas distributions do not significantly change our dynamical measurement values of ${M}_{\mathrm{BH}}$ , suggesting the model results are robust. Since we found the central attenuated hole of ${}^{12}\mathrm{CO}(2-1)$ , a factor of two higher spatial resolution observations of the $\mathrm{CS}(5-4)$ line are highly recommended to remove this large systemic uncertainty cleanly.

6. Thin Disk Tilted-ring Dynamical Model

We also constrain the M/L_H, i, and ${M}_{\mathrm{BH}}$ independently using a different dynamical model, that is, a thin disk tilted-ring model (Begeman 1987; Nicholson et al. 1992; Quillen et al. 1992; Neumayer et al. 2007; den Brok et al. 2015). We use the same kinematic measurements from the nuclear ${}^{12}\mathrm{CO}(2-1)$ gas and stellar-mass model derived in Sections 3.2.2 and 4, respectively. The purpose of this test is to examine the robustness of ${M}_{\mathrm{BH}}$ measurements in various assumptions and dynamical models.

6.1. The Tilted-ring Models

We model the kinematics of ${}^{12}\mathrm{CO}(2-1)$ CND with tilted-ring models using a similar approach as for the H₂1-0 S(1) transition kinematics in Seth et al. (2010), den Brok et al. (2015), and N17. The basic idea of this model is that we assume the emitting gas (e.g., CO or H₂) is rotating in thin rings on concentric circular orbits around the galaxy center with a velocity interpolated between the discrete points on the model grid linearly. Each ring of gas is described by three parameters: radius R, inclination angle i, and azimuthal angle θ (relative to the projected major axis) projected along the LOS allowing the rings to become ellipses. Particularly, we fit model to the data by assuming the ellipses change their geometry smoothly with radius. We determine the radially varying PA for this model with the Kinemetry routine (see footnote 10; Krajnović et al. 2006) but allow i to change linearly with radius.

We model the standard MGE deprojection by assuming a flattened mass distribution. The gravitational potential of each model includes the ${M}_{\mathrm{BH}}$ as a point source, the gas and dust distribution as described above in Section 5.1 parameterized by a new MGE, and the stellar-mass-component distribution, ${M}_{\star }(r)$ , which is modeled using the MGE in Table 3 with M/L_H as a free parameter identical to that used in the KinMS models.

During the fit, our model generates a spectral cube with the same spectral sampling as the ALMA data cube. The model contributes flux over the cube based on its spatial distribution, velocity field, and velocity dispersion for each ring across the coplanar disk and then replicates the observations. We calculate the ${\chi }^{2}$ for both the predicted rotational velocity and velocity dispersion fields but use only the ${\chi }^{2}$ of the rotational velocity field to determine our best-fit model. We note that the velocity errors are determined from the intensity-weighted mean velocity (mom 1) map directly. To compare the results of the tilted-ring model to the KinMS model's results, we use the same fitting area as for our KinMS models ( $0\buildrel{\prime\prime}\over{.} 8\times 0\buildrel{\prime\prime}\over{.} 8$ ). Moreover, we also correct for the strong pixel correlation in the ALMA data in the tilted-ring model by utilizing the third solution of scaling the input measurement uncertainties presented in Appendix B. Finally, we perform 10⁵ calculations and set 30% of these as the burn-in phase to find the best-fit tilted-ring model.

6.2. Results

We summarize the best-fit results of the tilted-ring model in Table 5. All uncertainties of the best-fit parameters quoted in this section are determined within 3σ CL. Here we quote the best-fit results derived at ${M}_{\mathrm{BH}}={1.1}_{-0.2}^{+0.4}\times {10}^{7}$ M_⊙ and $M/{L}_{H}\,={0.46}_{-0.21}^{+0.33}$ (M_⊙/ ${L}_{\odot }$ ), which has ${\chi }_{\mathrm{red}}^{2}=1.58$ for dof $\approx (64\times 64$ ). The ${M}_{\mathrm{BH}}$ is ∼31% lower, while the M/L_H is ∼4.3% higher, than those best-fit values found by the KinMS model. In Figure 14 we plot the best-fit velocity map of the ${}^{12}\mathrm{CO}(2-1)$ gas disk, its data, and the residual $({\mathtt{Data}}-{\mathtt{Model}})$ of the best-fit model. Similar to the KinMS model, the tilted-ring model finds evidence of some noncircular motions on the velocity residual map toward the southwest. We should also note that although the model describes the data at large radii ( $r\gtrsim 0.1^{\prime\prime}$ ) well, it tends to predict a steeper velocity toward the center due to the assumption of thin disks (d_t = 0). Figure 15 shows the 1D histogram and 2D marginalization PDF of a few important physical parameters of the tilted-ring models.

**Figure 14.** Comparison of the velocity fields of our best-fit tilted-ring model (left panel) with the ${}^{12}\mathrm{CO}(2-1)$ velocity field (center panel) and velocity residuals (right panel) with an area of $0\buildrel{\prime\prime}\over{.} 8\times 0\buildrel{\prime\prime}\over{.} 8$ (116 pc × 116 pc). The white contours denote the integrated intensity of ${}^{12}\mathrm{CO}(2-1)$ emission in the fitted area.
Download figure:
Standard image High-resolution image

**Figure 15.** The multidimensional parameter space posterior distributions explored by the tilted-ring dynamical model to the combined cube (high resolution) from the central $0\buildrel{\prime\prime}\over{.} 8\times 0\buildrel{\prime\prime}\over{.} 8$ FOV of NGC 3504. Other demonstrations are similar to Figure 11.
Download figure:
Standard image High-resolution image

**Figure 15.** The multidimensional parameter space posterior distributions explored by the tilted-ring dynamical model to the combined cube (high resolution) from the central $0\buildrel{\prime\prime}\over{.} 8\times 0\buildrel{\prime\prime}\over{.} 8$ FOV of NGC 3504. Other demonstrations are similar to Figure 11.
Download figure:
Standard image High-resolution image

Table 5. Best-fit Tilted-ring Model Parameters and Statistical Uncertainties for the Combined ALMA Data Cube (High Resolution)

	Range	Best Fit	1σ Error	3σ Error
(1)	(2)	(3)	(4)	(5)
${M}_{\mathrm{BH}}$	$6.0\to 8.0$	7.06	−0.10, +0.10	−0.30, +0.30
$M/{L}_{H}$	$0.0\to 5.2$	0.46	−0.07, +0.11	−0.21, +0.33
i	$30\to 90$	39.41	−4.88, +4.01	−14.20, +11.67
PA	$30\to 90$	82.29	−0.96, +0.92	−3.10, +2.90

Note. Column 1: a list of the fitted model parameters. Column 2: a list of the priors of the fitted model parameters and their search ranges. The priors are constructed in the uniform linear space with only SMBH in logarithmic space. Columns 3–5: the best-fit value of each parameter and their uncertainties at the 1σ and 3σ CLs. Units: ${M}_{\mathrm{BH}}$ (M_⊙), $M/{L}_{H}$ (M_⊙/ ${L}_{\odot }$ ), PA (°), and i (°); 1σ (16%–84%) and 3σ (0.14%–99.86%).

Download table as: ASCII Typeset image

7. Discussion

7.1. BH Mass and Scaling Relations

NGC 3504 has no previous ${M}_{\mathrm{BH}}$ measurement to compare with our gas-dynamical results except for the estimates based on the standard ${M}_{\mathrm{BH}}$ –σ relation for the bulge of LTGs (see Section 2.2). Our measurement is fully consistent with that prediction of the ${M}_{\mathrm{BH}}$ –σ scaling relations (see the right plot of Figure 16; Kormendy & Ho 2013; McConnell et al. 2013; Saglia et al. 2016).

**Figure 16.** Our NGC 3504 ${M}_{\mathrm{BH}}$ (red circle) in the context of the ${M}_{\mathrm{BH}}-{M}_{\mathrm{Bulge}}$ (left plot) and ${M}_{\mathrm{BH}}-\sigma$ (right plot) scaling relations. Seventeen measurements using both ALMA (Davis et al. 2013, 2017, 2018; Onishi et al. 2015; Barth et al. 2016a, 2016b; Boizelle et al. 2019; Combes et al. 2019; Smith et al. 2019; Nagai et al. 2019; North et al. 2019, T. Davis et al. 2020 in preparation; D. Nguyen et al. 2020 in preparation) and Combined Array for Research in Millimeter-wave Astronomy (CARMA; Davis et al. 2013; Onishi et al. 2017) observations are plotted in cyan, while ETGs (black dots within open circles) and LTGs (black open circles) are taken from Saglia et al. (2016). The scaling relations of Scott et al. (2013), Kormendy & Ho (2013), McConnell et al. (2013), and Saglia et al. (2016) for ETGs and LTGs are plotted in the dotted, dashed, and long-dashed lines, respectively. We also include the theoretical predictions of a bimodality from Pacucci et al. (2018a, red solid lines) and their 1σ (red dashed line) uncertainties. Bulge masses for objects with molecular gas estimated BH masses are taken from McConnell et al. (2013), Krajnović et al. (2013), Salo et al. (2015), Savorgnan et al. (2016), and Sani et al. (2018).
Download figure:
Standard image High-resolution image

We also examine our ${M}_{\mathrm{BH}}$ measurement of NGC 3504 in the context of the empirical compilations of ${M}_{\mathrm{BH}}-{M}_{\mathrm{Bulge}}$ and ${M}_{\mathrm{BH}}-\sigma$ scaling relations (Kormendy & Ho 2013; McConnell et al. 2013; Scott et al. 2013; Saglia et al. 2016) in Figure 16. The stellar-bulge velocity dispersion of NGC 3504 is determined from Ho et al. (2009), while the bulge mass of ${M}_{\mathrm{Bulge}}=2.5\,\times {10}^{9}$ M_⊙ was derived from Section 2.2. For a sanity check, we estimate the bulge mass of NGC 3504 using our H-band MGE model and adapting the effective radius of the bulge derived from the bulge-disk-bar decomposition model (Laurikainen et al. 2004). Accounting for the dynamical M/L_H (Section 5.2), we obtain ${M}_{\mathrm{Bulge}}=(3.5\pm 0.8)\times {10}^{9}$ M_⊙. Our result shows that the best-fit ${M}_{\mathrm{BH}}$ of NGC 3504 is fully consistent with the empirical ${M}_{\mathrm{BH}}-{M}_{\mathrm{Bulge}}$ and ${M}_{\mathrm{BH}}-\sigma$ relations of Kormendy & Ho (2013), McConnell et al. (2013), and Saglia et al. (2016), but outside $+1\sigma$ uncertainty of the Scott et al. (2013) and Savorgnan et al. (2016) empirical ${M}_{\mathrm{BH}}-{M}_{\mathrm{Bulge}}$ relations for "Sérsic" galaxies (those without central cores). Compared with the theoretical predictions of the bimodality in the BH accretion efficiency model (e.g., Pacucci et al. 2015, 2018), our measurement is consistent with the ${M}_{\mathrm{BH}}-\sigma$ correlation, but a positive outlier is the ${M}_{\mathrm{BH}}-{M}_{\mathrm{Bulge}}$ relation up to 1σ. At the mass of $\sim 1.6\times {10}^{7}$ M_⊙, the SMBH of NGC 3504 lies within the same mass regime as the BHs derived in the Combes et al. (2019) sample and lies between the samples of lower- and higher-mass galaxies previously and currently studied with ALMA (Davis et al. 2013, 2017, 2018; Onishi et al. 2015, 2017; Barth et al. 2016a, 2016b; Boizelle et al. 2019; Nagai et al. 2019; North et al. 2019; Smith et al. 2019, T. Davis et al. 2020, in preparation, D. Nguyen et al. 2020, in preparation), respectively. All of these works prove that the cold gas-dynamical method observed with ALMA at high spatial resolution now can work well in a wide range of BH masses covering six orders of magnitude from 10⁵ to ${10}^{10}{M}_{\odot }$ .

7.2. Angular Resolution Limit for Weighing Nearby Low-mass Galaxies with ALMA

ALMA observations provide a promising path for gas-dynamical ${M}_{\mathrm{BH}}$ measurement and exploration of BH demographics in various types of nearby low-mass galaxies across the Hubble sequence, especially for optical/IR opaque galaxies that host large fractions of gas and dust in their nuclei, making kinematics inaccessible at those wavelengths. As discussed in Barth et al. (2016b), to have accurate measurements of ${M}_{\mathrm{BH}}$ , the cold gas observations should satisfy the following criteria: (1) the CND has a simple disk-like rotation; (2) the gas kinematics are well supported by high SB line emission, and (3) the observational beam size should be at least as small as $2\times {R}_{\mathrm{SOI}}$ (e.g., Davis 2014; Boizelle et al. 2019). Rusli et al. (2013) and Davis (2014) suggest the ratio $\zeta =2{R}_{\mathrm{SOI}}/{\theta }_{\mathrm{FWHM}}$ of the SMBH's diameter of influence to the average beam size projected along the major axis, indicating the relative resolution of ${R}_{\mathrm{SOI}}$ . Observations with $\zeta \gtrsim 2$ are more reliable to produce a precise ${M}_{\mathrm{BH}}$ determination (Davis 2014; Boizelle et al. 2019). While observations with $\zeta \lt 2$ can yield an ${M}_{\mathrm{BH}}$ measurement, they are more susceptible to systematic bias from uncertain stellar masses, resulting in larger uncertainties on the ${M}_{\mathrm{BH}}$ (Kormendy & Ho 2013; Rusli et al. 2013; Barth et al. 2016a, 2016b; Boizelle et al. 2019). Particularly, targets with published CO imaging and ${M}_{\mathrm{BH}}$ measurements using ALMA and CARMA observations so far have relatively low values of $\zeta \lesssim 2$ for both ETGs (Davis et al. 2013, 2017, 2018; Onishi et al. 2017, 2015; Smith et al. 2019; Nagai et al. 2019) and LTGs (Combes et al. 2019) with the three ETG exceptions being the high–spatial resolution observations of NGC 1332 ( $\zeta \sim 10$ ; Barth et al. 2016a), NGC 3258 ( $\zeta \sim 17$ ; Boizelle et al. 2019), and NGC 0383 ( $\zeta \sim 15$ ; North et al. 2019). Our imaging of NGC 3504 is marginally resolved at ${R}_{\mathrm{SOI}}$ with $\zeta \sim 1.5$ for the high–spatial resolution observation. Testing with the low–spatial resolution case ( $\zeta \sim 0.3$ or ${\theta }_{\mathrm{FWHM}}=6.7\times {R}_{\mathrm{SOI}}$ ) gives consistent results on ${M}_{\mathrm{BH}}$ within 15% (Section 5.2). This similarity of our low-resolution data may suggest that reliable measurements on ${M}_{\mathrm{BH}}$ can be obtained even with lower-resolution data, although this result should be verified in additional systems using data at multiple resolutions.

In principle, a high–spatial resolution observation deeply within ${R}_{\mathrm{SOI}}$ is highly required for deriving ${M}_{\mathrm{BH}}$ with high precision and minimal susceptibility to systematic error, enabling a detailed dynamical model accounting for detailed disk structure (e.g., such as nearly edge-on or moderately warped disks; Herrnstein et al. 2005; Humphreys et al. 2013; Gao et al. 2016), isolating the locally irregular kinematics, and being less sensitive to the host galaxy mass profile. However, high–spatial resolution observations often reveal holes in the disk (Davis et al. 2017; North et al. 2019; Smith et al. 2019) and thus do not always provide significantly better constraints on the ${M}_{\mathrm{BH}}$ . As discussed in Section 3.2.2, dense gas tracers may be good alternative transitions for ${M}_{\mathrm{BH}}$ measurements in case of ${}^{12}\mathrm{CO}(2-1)$ deficits, especially for low-mass galaxies, which require very high–spatial resolution observations because of their small ${R}_{\mathrm{SOI}}$ . This problem can be solved by an efficient observational strategy of including dense gas tracers such as $\mathrm{CS}(5-4)$ or CH₃OH, which can be observed simultaneously with the ${}^{12}\mathrm{CO}(2-1)$ line in a single exposure. On the other hand, beam smearing makes the kinematics along the major and minor axes degenerate, as well as between rotation and dispersion at the very central region of a CND (Barth et al. 2016a, 2016b; Boizelle et al. 2019), which creates difficulties for data versus model optimization.

In the near future, observations of galaxies with dominant rotation within the SOI will increase rapidly. Measuring ${M}_{\mathrm{BH}}$ in such galaxies is important to anchor the local BH demographics and BH–host galaxy correlations (Maoz et al. 1998) and constrain the BH seed formations in the early universe (e.g., Volonteri et al. 2008, 2015; Greene 2012; Bonoli et al. 2014; Fiacconi & Rossi 2016, 2017; Ricarte & Natarajan 2018a, 2018b; Bellovary et al. 2019). Also, the growing number of high-resolution observations of ALMA for nearby targets will provide an excellent opportunity to examine molecular gas kinematics in galaxy nuclei in far greater detail than any previous surveys.

8. Conclusions

We present a dynamical mass measurement for the SMBH in NGC 3504 using ${}^{12}\mathrm{CO}(2-1)$ emission observed by ALMA and HST imaging. We highlight our main results below:

1.
We measured a distance to the galaxy NGC 3504 of 32.4 ± 2.1 Mpc, further than the existing distance estimates reported on NED.
2.
NGC 3504 hosts a CND cospatially distributed with obscuring dust lanes visible in the HST imaging within $5^{\prime\prime}$ .
3.
Both KinMS and tilted-ring dynamical models suggest the detection of a central SMBH whose mass and M/L are constrained simultaneously using the HST/H-band image at ${M}_{\mathrm{BH}}\sim 1.6\times {10}^{7}$ M_⊙ and $M/{L}_{H}\sim 0.45$ (M_⊙/ ${L}_{\odot }$ ). The agreement between the two models suggests they are both powerful tools to probe the BH–galaxy scaling relations and BH demographics, especially at the low-mass regimes of both SMBHs and their hosts. Moreover, such a low M/L_H in the nucleus of NGC 3504 is typically seen in normal LTGs (Martinsson et al. 2013; McGaugh & Schombert 2014), suggesting there are young stars present in the nucleus.
4.
The BH mass is consistent with the empirical ${M}_{\mathrm{BH}}\,-{M}_{\mathrm{Bulge}}$ and ${M}_{\mathrm{BH}}-\sigma$ correlations. Compared with the theoretical predictions of the bimodality in BH accretion efficiency model, our measurement is consistent with the ${M}_{\mathrm{BH}}-\sigma$ correlation, but ∼+1σ consistent with the ${M}_{\mathrm{BH}}-{M}_{\mathrm{Bulge}}$ relation.
5.
Both combined (high spatial resolution) and low–spatial resolution observations give consistent constraints on the ${M}_{\mathrm{BH}}$ and M/L_H with systemic <22% and ∼2%, respectively, although their observational scales are ∼5× different from each other. This result confirms our observational strategy for measuring ${M}_{\mathrm{BH}}$ in nearby low-mass galaxies to explore their BH demographics and anchor the BH–galaxy scaling relations at the low-mass ends in a large sample, which cannot be achieved dynamically with stars or ionized/warm gas.
6.
The gas in the CND of NGC 3504 has a relatively high velocity dispersion with ∼39 km s⁻¹ in the region $0\buildrel{\prime\prime}\over{.} 1\,\lt r\lt 0\buildrel{\prime\prime}\over{.} 3$ and ∼25 km s⁻¹ at larger radii ( $0\buildrel{\prime\prime}\over{.} 3\lt r\lt 0\buildrel{\prime\prime}\over{.} 5$ ). The centrally high value (>60 km s⁻¹, $r\lt 0\buildrel{\prime\prime}\over{.} 1$ ) of dispersion is not intrinsic but is a beam smearing effect along with the LOS integration through the nearly edge-on orientation of the CND. This low-velocity dispersion of the nuclear gas in an LTG like NGC 3504 is similar to what is found in ETGs (e.g., Barth et al. 2016a; Davis et al. 2017; Boizelle et al. 2019), suggesting an insignificant role of turbulent velocity dispersions on our gas-dynamical models and proving that high–spatial resolution observations of ALMA can be applied for lower-mass LTGs to measure their SMBH masses accurately.
7.
We find a central attenuated hole that has a radius of 6.3 pc in the ${}^{12}\mathrm{CO}(2-1)$ integrated intensity map. However, the dense gas tracer $\mathrm{CS}(5-4)$ found in one continuum SPW is centrally peaked, fills this hole, and has similar kinematic features to the ${}^{12}\mathrm{CO}(2-1)$ emission line within the CND. The origin of this ${}^{12}\mathrm{CO}(2-1)$ hole could be the changing of excitation conditions and not a true absence of molecular gas.

The authors would like to thank the anonymous referee for his/her careful reading and useful comments, which helped to improve the paper greatly. We also would like to thank NAOJ and NINS for supporting this work. D.D.N. expresses his gratitude for the support from the Willard L. and Ruth P. Eccles Foundation for their Eccles Fellowship during the 2017–2018 academic year at the Department of Physics and Astronomy of the University of Utah and the TAIZAI Visiting Fellowship during the Spring 2018 at NAOJ. A.C.S. acknowledges financial support from the National Science Foundation (NSF) grant AST-1350389, T.A.D. acknowledges support from a Science and Technology Facilities Council Ernest Rutherford Fellowship, and M.C. acknowledges support from a Royal Society University Research Fellowship (RSURF). S.T. acknowledges the TAIZAI Visiting Fellowship during the Spring 2018 at NAOJ. T.I. sincerely thanks the Japan Society for the Promotion of Science (JSPS) KAKENHI, grant No. 17K14247. M.I. is supported by JSPS KAKENHI grant No. 15K05030. In addition, the authors greatly acknowledge publication support from the ALMA Japan Research Grant of the NAOJ Chile Observatory, NAOJ-ALMA-3313.

The authors would like to thank Doctors Naoteru Gouda, Makoto Miyoshi, Yano Taihei, Taiki Kawamuro, and Daisuke Iono at NAOJ, as well as Professor Aaron J. Barth of the University of California Irvine for enlightening discussions and Dr. Benjamin Boizelle of the George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas ${\rm{A}}\&{\rm{M}}$ University, for sharing code and discussions. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2017.1.00964.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We thank the ALMA operators and staff and the ALMA help desk as well for diligent feedback and invaluable assistance in processing these data.

Facility: - ALMA and HST WFC3.

Software: IDL, CASA, astropy, and emcee.

Appendix A: The Impact of Disk Structure

The assumption of molecular gas moving on circular orbits is valid only for the thin disk model where $\sqrt{{\sigma }^{2}+{V}_{\mathrm{rot}}^{2}}\approx {V}_{\mathrm{circ}}$ (Barth et al. 2016a) so that the dynamical pressure distributed effectively by the turbulent velocity dispersion of the gas that supports the disk against gravity is negligible. The nuclear gas disk is then treated as dynamically flat (disk thickness ${d}_{t}\approx 0$ ) and cold ( $\sigma \approx 0$ , ${V}_{\mathrm{rot}}={V}_{\mathrm{circ}}$ ). For NGC 3504, we interpret the central turbulent velocity dispersion of the ${}^{12}\mathrm{CO}(2-1)$ gas as a beam smearing effect within the radius of the central hole ( $r\lesssim 0\buildrel{\prime\prime}\over{.} 04$ ); then we can assume the velocity dispersion in this region is identical to the dispersion of the CND. However, as discussed in Section 3.2.2 the CND turbulent velocity dispersion is relatively high at ∼39 km s⁻¹, suggesting the thin disk structure assumption may not correct for the nuclear gas in NGC 3504. We therefore examine the $\sigma /{V}_{\mathrm{rot}}$ ratio within the CND extracted along the CND's major axis; the disk's physical thickness depends on the profile reaching the mean value of 0.34 at $0\buildrel{\prime\prime}\over{.} 4$ (48 pc). Therefore, we include a spatial constant thickness parameter (d_t) in the KinMS model to describe the disk dominated by gas velocity dispersion.

**Figure 17.** The channel maps of the detected ${}^{12}\mathrm{CO}(2-1)$ emission in our combined (high–spatial resolution) cube overlaid with the best-fit KinMS model. The colored regions with black contours show the areas detected with more than 2.5σ significance, while the thick red contours are plotted with the same contour levels from our best-fitting KinMS model.
Download figure:
Standard image High-resolution image

Appendix B: Accounting for Spatial Correlation in the ALMA Data

Each pixel of ALMA data is strongly correlative within the synthesized beam due to the Nyquist sampling of the spatial dimensions (Barth et al. 2016a; Davis et al. 2017, 2018; Onishi et al. 2017; North et al. 2019; Smith et al. 2019). In other words, the noise covariance mostly corresponds to a local correlation among neighboring pixels, on the scale of the synthesized beam FWHM. The pixels closer than one beam FWHM will be strongly correlated, while pixels more distant than that from each other will be mostly independent. This issue makes our model parameter uncertainties unphysically small; however, we solve this by using the following approaches:

(1)
Method 1: using a response function that characterizes the ALMA beam of the observed cube, we convolve this function to a "true" spaxel uncorrelated image to mimic the observational image. Davis et al. (2017) develop a functional form to relate the response function and the rms noise of the data cube analytically by using a full covariance matrix to calculate the ${\chi }^{2}$ and likelihood distribution (Equation (18) in Cappellari 2017). In practice, because the condition number of the covariance matrix is too large, its inverse-likelihood calculation is impossible. To keep our numerical calculation accurate when calculating the inversion, we apply an alternative approach of using a modified Cholesky factorization step (Gill et al. 1981) and restrict the fits within an optimally small area of 64 pixels × 64 pixels in the MCMC framework.
(2)
Method 2: accounting for the synthesized beam correlation by rebinning the data spatially so that the spatial scale of the obtained cube corresponds roughly to 1 pixel per synthesized beam. The spatial noise in the binned cube is therefore much closer to Gaussian random noise and there is much less covariance between pixels at this rebinned larger pixel scale. We then measure the rms noise from the rebinned cube and use it to calculate the ${\chi }^{2}$ for each model fit.
(3)
Method 3: utilizing an approximate approach introduced in van den Bosch & van de Ven (Section 3.2 of 2009) in which one tries to crudely account for systematic errors. This method assumes that the systematic uncertainties give a similar contribution to the ${\chi }^{2}$ as the statistical uncertainties. Under this assumption, one increases the ${\rm{\Delta }}{\chi }^{2}$ CL for an acceptable fit by the variance of the ${\chi }^{2}$ , which is $2(N-P)\sim 2N$ (for large $N\gt {10}^{5}$ ), where N is the number of constraints and P are the number of model parameters. This implies that a CL becomes defined by the contour level with ${\rm{\Delta }}{\chi }^{2}=\sqrt{2N}$ . The same idea can be applied in our Bayesian framework by multiplying the likelihood by a factor $\sqrt{2N}$ or equivalently scaling the input measurement uncertainties by a factor of ${(2N)}^{0.25}$ (Mitzkus et al. 2017).

We test these ideas for NGC 3504 and record their inferred parameters' uncertainties in Tables 4 and 6. Of these three methods, the first provides the smallest errors on the estimated parameters, while the last provides the larger errors. To be conservative, we therefore use the errors from the final method in our estimates below. This approach is consistent with previous works (Mitzkus et al. 2017; Nagai et al. 2019; North et al. 2019; Smith et al. 2019).

Table 6. Best-fit KinMS Model Parameters and Statistical Uncertainties for the Low-resolution ALMA Cube within the FOV $0\buildrel{\prime\prime}\over{.} 8\times 0\buildrel{\prime\prime}\over{.} 8$

Parameters (Units)	Search Range			Best Fit	1σ Error³	3σ Error³	3σ Error²	3σ Error¹
		Uniform			(16%–84%)	(0.14%–99.86%)	(0.14%–99.86%)	(0.14%–99.86%)
(1)		(2)		(3)	(4)	(5)	(6)	(7)
Black Hole:
${\mathrm{log}}_{10}{M}_{\mathrm{BH}}$ (M_⊙)	4.00	$\longrightarrow$	9.00	7.30	−0.06, +0.06	−0.20, +0.20	−0.08, +0.08	−0.14, +0.14
M/L_H (M_⊙/ ${L}_{\odot }$ )	0.01	$\longrightarrow$	3.00	0.43	−0.04, +0.04	−0.12, +0.12	−0.06, +0.06	−0.10, +0.10
f (Jy km s⁻¹)	10.0	$\longrightarrow$	500.0	149.52	−1.36, +1.38	−4.03, +4.21	−1.67, +1.67	−3.18, +3.17
Molecular Gas:
$i{(}^{^\circ })$	45.0	$\longrightarrow$	90.0	59.18	−1.02, +1.06	−3.06, +3.18	−1.36, +1.35	−2.28, +2.26
${\sigma }_{0}$ (km s⁻¹)	1.0	$\longrightarrow$	80.0	34.61	−1.00, +0.99	−3.00, +3.01	−1.52, +1.53	−2.68, +2.67
d_t ('')	0.01	$\longrightarrow$	1.50	0.05	−0.01, +0.01	−0.03, +0.03	−0.00, +0.00	−0.01, +0.01
G ${}_{\sigma }$ ('')	0.01	$\longrightarrow$	1.00	0.23	−0.01, +0.01	−0.03, +0.03	−0.01, +0.01	−0.02, +0.02
Nuisance:
x_c ('')	−0.10	$\longrightarrow$	+0.10	−0.01	−0.00, +0.00	−0.01, +0.01	−0.00, +0.00	−0.01, +0.01
y_c ('')	−0.10	$\longrightarrow$	+0.10	−0.03	−0.00,+0.00	−0.01, +0.01	−0.00, +0.00	−0.01, +0.01
${v}_{\mathrm{off}.}$ (km s⁻¹)	−10.0	$\longrightarrow$	+10.0	7.39	−0.30, +0.31	−0.92, +0.92	−0.35, +0.35	−0.63, +0.65

Note. All notations and parameters in this table are the same those in Table 4.

Download table as: ASCII Typeset image

Appendix C: Supplementary Figures and Tables

Figure 17 shows the channel maps from the best-fit model overlaid on the channel maps of the data. Figure 18 demonstrates the same best-fit KinMS results for the low-spatial resolution data. The full list of the best-fit parameters of the KinMS model and their likelihood for the low-spatial resolution data are presented in Table 6.

The MBHBM_⋆ Project. I. Measurement of the Central Black Hole Mass in Spiral Galaxy NGC 3504 Using Molecular Gas Kinematics

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. The Galaxy NGC 3504

2.1. Determination of the Distance

2.2. Central Properties

3. Data and Data Reduction

3.1. HST Imaging

3.2. ALMA Observations

3.2.1. Continuum Emission

3.2.2. The ¹²CO(2–1) and CS(5–4) Emission Lines

4. Creating a Mass Model

5. KinMS Dynamical Modeling

5.1. Fitting Process

5.2. Model Results

5.3. Uncertainty Sources on the ${M}_{\mathrm{BH}}$ estimate

5.3.1. Stellar-Mass Models

5.3.2. Distances

5.3.3. The Galaxy Bar

5.3.4. Turbulent Velocity Dispersion of the Gas

5.3.5. Different Observational Scales and Distributed Assumption of the Gas

6. Thin Disk Tilted-ring Dynamical Model

6.1. The Tilted-ring Models

6.2. Results

7. Discussion

7.1. BH Mass and Scaling Relations

7.2. Angular Resolution Limit for Weighing Nearby Low-mass Galaxies with ALMA

8. Conclusions

Appendix A: The Impact of Disk Structure

Appendix B: Accounting for Spatial Correlation in the ALMA Data

Appendix C: Supplementary Figures and Tables

Footnotes

The MBHBM⋆ Project. I. Measurement of the Central Black Hole Mass in Spiral Galaxy NGC 3504 Using Molecular Gas Kinematics

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. The Galaxy NGC 3504

2.1. Determination of the Distance

2.2. Central Properties

3. Data and Data Reduction

3.1. HST Imaging

3.2. ALMA Observations

3.2.1. Continuum Emission

3.2.2. The 12CO(2–1) and CS(5–4) Emission Lines

4. Creating a Mass Model

5. KinMS Dynamical Modeling

5.1. Fitting Process

5.2. Model Results

5.3. Uncertainty Sources on the {M}_{\mathrm{BH}} estimate

5.3.1. Stellar-Mass Models

5.3.2. Distances

5.3.3. The Galaxy Bar

5.3.4. Turbulent Velocity Dispersion of the Gas

5.3.5. Different Observational Scales and Distributed Assumption of the Gas

6. Thin Disk Tilted-ring Dynamical Model

6.1. The Tilted-ring Models

6.2. Results

7. Discussion

7.1. BH Mass and Scaling Relations

7.2. Angular Resolution Limit for Weighing Nearby Low-mass Galaxies with ALMA

8. Conclusions

Appendix A: The Impact of Disk Structure

Appendix B: Accounting for Spatial Correlation in the ALMA Data

Appendix C: Supplementary Figures and Tables

Footnotes

The MBHBM_⋆ Project. I. Measurement of the Central Black Hole Mass in Spiral Galaxy NGC 3504 Using Molecular Gas Kinematics

3.2.2. The ¹²CO(2–1) and CS(5–4) Emission Lines

5.3. Uncertainty Sources on the ${M}_{\mathrm{BH}}$ estimate