The MBHBM* Project -- I: Measurement of the Central Black Hole Mass in Dwarf Galaxy NGC 3504 Using Molecular Gas Kinematics

TL;DR

This study measures the supermassive black hole mass in galaxy NGC 3504 using ALMA molecular gas observations, revealing a black hole consistent with galaxy scaling relations and identifying a dense gas tracer for improved measurements.

Contribution

First measurement of SMBH mass in NGC 3504 using ALMA CO(2-1) data, introducing CS(5-4) as a potential more accurate tracer for central black hole mass estimation.

Findings

Black hole mass estimated at ~10^7 solar masses.

Detection of a dense gas tracer CS(5-4) at galaxy center.

Black hole mass aligns with galaxy scaling relations.

Abstract

We present the first measurement of the mass of a supermassive black hole (SMBH) in the nearby double-barred spiral galaxy NGC 3504 as part of the Measuring Black Holes Below the Milky Way ($M_{\star}$) mass galaxies (MBHBM$_{\star}$) Project. Our analysis is based on Atacama Large Millimeter/submillimeter Array (ALMA) Cycle-5 observations of the ${\rm ^{12}CO(2-1)}$ emission line. NGC 3504 has a circumnuclear gas disk (CND), which has a relatively high-velocity dispersion of 30 \kms. Our dynamical models of the CND yield a \Mbh~of $M_{\rm BH}=1.02^{+0.18}_{-0.15}\times10^7$\Msun and a mass-to-light ratio in $H$-band of \ml$_{\rm H}=0.66^{+1.44}_{-0.65}$ (\Msun/\Lsun). This black hole (BH) mass is consistent with BH--galaxy scaling relations. We also detect a central deficiency in the ${\rm ^{12}CO(2-1)}$ integrated intensity map with a diameter of 2.7 pc at the putative position of the SMBH. However, this hole is filled by a dense gas tracer ${\rm CS(5-4)}$ that peaks at the galaxy center found in one of the three low-velocity-resolution continuum spectral correlators. The ${\rm CS(5-4)}$ line has the same kinematics with the ${\rm ^{12}CO(2-1)}$ line within the CND, suggesting that it is also an alternative transition for measuring the central \Mbh~in NGC 3504 probably more accurately than the current commonly used of ${\rm ^{12}CO(2-1)}$ due to its centralization.

Figures

Tables

Table 1. Table 1: ALMA observation parameters.

Phase center	R.A.	Decl.
	${11}^h{03}^m{11}^s.205$	$+{27}^{\circ }{58}^{\prime }20\stackrel{\prime \prime }{.}80$
	High Res.	Low Res.
Configurations	C43-9	C43-6
Obs. Date	2017 Oct. 24	2018 Jan. 01
Exposure time	46.6 min.	25.7 min.
Beamsize	$0\stackrel{\prime \prime }{.}042\times 0\stackrel{\prime \prime }{.}030$	$0\stackrel{\prime \prime }{.}221\times 0\stackrel{\prime \prime }{.}164$
(or FWHM)	2.9 pc $\times$ 2.0 pc	15.0 pc $\times$ 11.2 pc
Beam PA	$-$2∘.0	32∘.3
Velocity resolution	1.5 km s-1	39.5 km s-1
Frequency resolution	1.13 MHz	31.25 MHz

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

$j$	Total Count	$\sigma$	$a/b$
	of Gaussianj	(arcsec)
(1)	(2)	(3)	(4)
1	0.343	0.049	0.999
2	0.538	0.130	0.999
3	0.060	0.403	0.999
4	0.033	0.897	0.996
5	0.033	1.716	0.993

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

$j$	$\log$(Luminosity Density)	${\sigma }^{\prime }$	$a/b$
	($L_{\odot }/{\mathrm{pc}}^2$)	(arcsec)
(1)	(2)	(3)	(4)
1⋆	4.641	0.049	0.889
2	4.639	0.210	1.000
3	4.268	0.546	1.000
4	4.217	2.649	0.780
5	3.051	4.325	1.000
6	3.016	9.147	0.650
7	2.620	23.588	0.650
8	1.949	78.853	1.000

Table 4. Table 4: Best-fit KinMS Model Parameters and Statistical Uncertainties for the Combined ALMA Data Cube (High Resolution)

Parameter Names (Notations; Units)	Search Range			Best Fit	$1\sigma$ Error	$3\sigma$ Error
					(68% conf.)	(99.7% conf.)
(1)		(2)		(3)	(4)	(5)
Black Hole:
$M_{\mathrm{BH}}$ in log scale (${\log }_{10}{M}_{\mathrm{BH}}$; $M_{\odot }$)	4.00	$⟶$	9.00	7.01	$-$0.01, +0.01	$-$0.07, +0.07
Mass-to-light ratio in $H$-band (M/LH; $M_{\odot }$/$L_{\odot }$)	0.01	$⟶$	3.00	0.66	$-$0.32, +0.32	$-$0.65, +1.44
Molecular Gas Disc:
Position angle (PA; ∘)	150.0	$⟶$	360.0	334.81	$-$0.50, +0.45	$-$3.04, +2.03
Inclination angle ($i$; ∘)	45.0	$⟶$	90.0	62.13	$-$0.46, +0.49	$-$3.51, +3.11
Gas velocity dispersion ($\sigma$; km s-1)	1.0	$⟶$	80.0	38.30	$-$0.38, +0.37	$-$0.06, +0.06
Disk thickness ($d_t$; arcsec)	0.01	$⟶$	1.50	0.32	$-$0.02, +0.02	$-$0.05, +0.05
Gaussian peak (Gc; arcsec)	0.01	$⟶$	1.00	0.04	$-$0.00, +0.00	$-$0.02, +0.02
Gaussian HWHM (Gh; arcsec)	0.01	$⟶$	1.00	0.11	$-$0.01, +0.01	$-$0.03, +0.03
Exponential disk length scale ($R_e$; arcsec)	0.01	$⟶$	10.0	1.50	$-$0.00, +0.00	$-$0.02, +0.02
Nuisance Parameters:
CO surface brightness scaling factor (Flux)	10.0	$⟶$	500.0	97.26	$-$0.52, +0.51	$-$2.21, +2.75
R.A. offset ($x_c$; arcsec)	$-$0.10	$⟶$	+0.10	$-$0.01	$-$0.00, +0.00	$-$0.01, +0.01
Decl. offset ($y_c$; arcsec)	$-$0.10	$⟶$	+0.10	$-$0.02	$-$0.00, +0.00	$-$0.01, +0.01
Systemic velocity offset ($v_{\mathrm{off}.}$; km s-1)	$-$10.0	$⟶$	+10.0	$-$5.46	$-$0.69, +0.59	$-$2.61, +1.83

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

Parameter Names (Notations; Units)	Search Range			Best Fit	$1\sigma$ Error	$3\sigma$ Error
					(68% conf.)	(99.7% conf.)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
Model 1:
$M_{\mathrm{BH}}$ in log scale ($\log M_{\mathrm{BH}}$; $M_{\odot }$)	6.0	$⟶$	8.0	6.91	$-$0.01, +0.01	$-$0.05, +0.04
Mass-to-light ratio in $H$-band ($M/L_H$; $M_{\odot }$/$L_{\odot }$)	${10}^{-2}$	$⟶$	5.2	0.62	$-$0.01 +0.01	$-$0.04, +0.04
Position angle of the inner disk (PA1; ∘)	65	$⟶$	90	76.20	$-$0.17, +0.13	$-$0.53, +0.30
Position angle of the outer disk (PA1; ∘)	65	$⟶$	90	65.02	$-$0.01, +0.02	$-$0.07, +0.07
Inclination angle of the inner disk ($i_1$; ∘)	25	$⟶$	90	54.75	$-$0.17, +0.18	$-$0.60, +0.52
Inclination angle of the outer disk ($i_2$; ∘)	25	$⟶$	90	61.84	$-$0.20, +0.20	$-$0.60, +0.62
Transition radius between two disks (Rbreak; arcsec)	${10}^{-2}$	$⟶$	1.0	0.29	$-$0.0, +0.0	$-$0.04, +0.04
Model 2:
$M_{\mathrm{BH}}$ in log scale ($\log M_{\mathrm{BH}}$; $M_{\odot }$)	6.0	$⟶$	8.0	6.88	$-$0.02, +0.02	$-$0.06, +0.07
Mass-to-light ratio in $H$-band ($M/L_H$; $M_{\odot }$/$L_{\odot }$)	${10}^{-2}$	$⟶$	5.2	0.63	$-$0.03, +0.02	$-$0.08, +0.07
Inclination angle of the inner ellipse ($i_1$; ∘)	30	$⟶$	90	47.54	$-$1.20, +1.34	$-$3.72, +3.92
Inclination angle of the outer ellipse ($i_2$; ∘)	30	$⟶$	90	60.42	$-$1.06, +0.99	$-$1.21, +1.82
Model 3:
$M_{\mathrm{BH}}$ in log scale ($\log M_{\mathrm{BH}}$; $M_{\odot }$)	6.0	$⟶$	8.0	6.94	$-$0.00, +0.04	$-$0.06, +0.13
Mass-to-light ratio in $H$-band ($M/L_H$; $M_{\odot }$/$L_{\odot }$)	${10}^{-2}$	$⟶$	5.2	0.61	$-$0.00 +0.03	$-$0.03, +0.09
Position angle of the inner disk (PA; ∘)	65	$⟶$	90	65.51	$-$0.83, +0.94	$-$2.35, +2.90
Inclination angle of the inner disk ($i$; ∘)	50	$⟶$	90	57.77	$-$0.96, +1.99	$-$3.23, +5.65
Central hole radius (rinner; arcsec)	${10}^{-2}$	$⟶$	0.1	0.04	$-$0.01, +0.01	$-$0.03, +0.03

Table 6. Table 6: Best-fit KinMS Model Parameters and Statistical Uncertainties for the Low-Resolution ALMA Cube

Parameter Names (Notations; Units)	Search Range			Best Fit	$1\sigma$ Error	$3\sigma$ Error
					(68% conf.)	(99.7% conf.)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
Black Hole:
$M_{\mathrm{BH}}$ in log scale (${\log }_{10}{M}_{\mathrm{BH}}$; $M_{\odot }$)	4.00	$⟶$	9.00	7.18	$-$0.04, +0.04	$-$0.10, +0.13
Mass-to-light ratio in $H$-band (M/LH; $M_{\odot }$/$L_{\odot }$)	0.01	$⟶$	3.00	0.67	$-$0.02, +0.03	$-$0.62, +1.01
Molecular Gas Disc:
Position angle (PA; ∘)	150.0	$⟶$	360.0	330.17	$-$1.42, +0.90	$-$4.24, +3.25
Inclination angle ($i$; ∘)	45.0	$⟶$	90.0	60.31	$-$1.43, +1.34	$-$4.51, +4.01
Gas velocity dispersion ($\sigma$; km s-1)	1.0	$⟶$	80.0	35.50	$-$0.53, +0.45	$-$1.56, +1.42
Disk thickness ($d_t$; arcsec)	0.01	$⟶$	1.50	0.41	$-$0.03, +0.03	$-$0.09, +0.12
Gaussian HWHM (Gh; arcsec)	0.01	$⟶$	1.00	0.35	$-$0.03, +0.03	$-$0.10, +0.13
Nuisance Parameters:
CO surface brightness scaling factor (Flux)	10.0	$⟶$	500.0	149.01	$-$0.02, +0.02	$-$0.07, +0.09
R.A. offset ($x_c$; arcsec)	$-$0.10	$⟶$	+0.10	$-$0.02	$-$0.00, +0.00	$-$0.01, +0.01
Decl. offset ($y_c$; arcsec)	$-$0.10	$⟶$	+0.10	$-$0.03	$-$0.00, +0.00	$-$0.01, +0.01
Systemic velocity offset ($v_{\mathrm{off}.}$; km s-1)	$-$10.0	$⟶$	+10.0	$-$0.02	$-$0.22, +0.02	$-$0.06, +0.07

Equations

$$v_c^2(r)=r\frac{d\Phi (r)}{dr}$$

$$v_c^2(r)=G(M_{\star }(<r)+M_{\mathrm{BH}})r.$$

Full text

galaxy: individual NGC 3504 – galaxy: kinematics and dynamics – galaxy: supermassive black hole – galaxy: interstellar medium – galaxy: spirals – galaxy: active galactic nuclei.

††journal: ApJ

1 Introduction

Supermassive black holes (SMBHs, $M_{\rm BH}\gtrsim 10^{6}M_{\odot}$) are believed to reside at the centers of massive galaxies (M_{\star}\gtrsim 10^{11}$$M_{\odot}) and their masses correlate to macroscopic properties of the host (e.g., bugle velocity dispersion, bulge mass, and bulge luminosity). These remarkable discoveries are based on large efforts of observations (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 works 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; Barai et al., 2014; Netzer, 2015). The underlying physics involve physically small scales but gravitationally large influences of SMBHs on their galactic environments that are recorded in the $M_{\rm BH}$–galaxy scaling relations. These correlations suggest SMBHs may play a pivotal role in the growth and evolution of galaxies (e.g., Schawinski et al., 2007; Kormendy & Ho, 2013; Saglia et al., 2016; van den Bosch et al., 2016).

The observational census on the scaling relations in the regime of low-mass galaxies is very incomplete, as direct measurements of black hole (BH111In this article, we use the acronyms of SMBH and BH exchangeable.) masses are lacking. This includes the increased scatter around the relations for the Milky Way late-type galaxies (LTGs; Greene et al., 2010, 2016; Läsker et al., 2016) and roughly 2–3 orders of magnitude below the bulge mass relation for the lowest mass galaxies (Scott et al., 2013; Graham & Scott, 2015; Nguyen et al., 2017; Chilingarian et al., 2018; Nguyen et al., 2018, 2019, hereafter N18, N19). Many hypotheses have been proposed to explain this change that may be due to (i) the formation history of the bulge (Kormendy & Bender, 2012; Krajnović et al., 2018), (ii) the star formation history (SFH) of the galaxy (Caplar et al., 2015; Terrazas et al., 2017), or (iii) the bimodal population of high-$z$ BH seeds in which the light seeds ($M_{\rm BH}\lesssim 10^{3}M_{\odot}$) accrete inefficiently with largely sub-Eddington rates and duty cycles, while the massive seeds ($M_{\rm BH}\gtrsim 10^{4}M_{\odot}$) grow efficiently (Pacucci et al., 2015; Inayoshi & Haiman, 2016; Park et al., 2016; Pacucci et al., 2017, 2018; Pacucci & Loeb, 2018)

Recently, the number of known intermediate-mass black holes (10^{3}$$M_{\odot}$$<M_{\rm BH}\lesssim 10^{6}$$M_{\odot}), which are so called million/sub-million Solar masses BHs, are increasing dramatically with their masses inferred from a variety of methods including: (1) the velocity widths of their optical broad-line emissions (Barth et al., 2004; Greene & Ho, 2007; Thornton et al., 2008; Dong et al., 2012; Reines et al., 2013; Baldassare et al., 2015; Reines & Volonteri, 2015; Chilingarian et al., 2018), (2) the accretion signatures of the narrow-line emissions (e.g., Moran et al., 2014) and the coronal emission in the mid-infrared (MIR; Satyapal et al., 2009), (3) tidal-disruption events (TDE; e.g., Maksym et al., 2013; Stone et al., 2017), (4) hard X-ray emission (e.g., Gallo et al., 2008; Desroches et al., 2009; Gallo et al., 2010; Miller et al., 2015; She et al., 2017), (5) the dynamics of accretion disks containing megamasers (Miyoshi et al., 1995; Lo, 2005; Kuo et al., 2011; van den Bosch et al., 2016), and (6) the dynamics of nuclear stars and gas to measure $M_{\rm BH}$ in low-mass galaxies (5\times 10^{8}$$M_{\odot}$$<M_{\star}\lesssim 10^{10}$$M_{\odot}; Verolme et al., 2002; Valluri et al., 2005; Neumayer et al., 2007; van den Bosch & de Zeeuw, 2010; Seth et al., 2010; den Brok et al., 2015; Nguyen et al., 2017; Thater et al., 2017, N18, N19) and ultracompact dwarfs (UCDs, 1\times 10^{7}$$M_{\odot}$$<M_{\star}\leq 5\times 10^{8}$$M_{\odot}; Seth et al., 2014; Ahn et al., 2017; Afanasiev et al., 2018; Ahn et al., 2018; Voggel et al., 2018). The importance of the \lesssim$$10^{6}M_{\odot} BH population and their masses are discussed in detail in N18 and N19.

Atacama Large Millimeter/sub-millimeter Array (ALMA) observations of molecular gas at mm/sub-mm wavelengths offer a promising way to characterize the full spectrum of BH populations across the Hubble sequence from LTGs to early-type galaxies (ETGs; Davis et al., 2013; Onishi et al., 2015; Barth et al., 2016a, b; Davis et al., 2017; Onishi et al., 2017; Davis et al., 2018) because of the high angular resolution and sensitivity. Moreover, ALMA can improve on many problems affecting $M_{\rm BH}$ estimates via existing optical/infrared instruments including (1) high angular resolution capable of resolving the sphere of influence (SOI) and (2) ability to obtain dynamical measurements in dusty/obscured nuclei, which are inaccessible at optical wavelengths. The physical idea behind the dynamical method utilizing the observations of ALMA is that the $M_{\rm BH}$ is derived by detecting the Keplerian turnover motion of the cold gas disk at the galactic center directly. However, the cold gas dynamical method at mm/sub-mm works well for gas-rich galaxies those host well-defined and relaxed circumnuclear gas disks (CNDs) only.

This article is the first of a series in the Measuring Black Holes Below the Milky Way ($M_{\star}$) mass galaxies (MBHBM*⋆*) project. This project aims to gather a large sample of low-mass gas-rich galaxies and measure their central dark masses, which are likely BHs, using molecular gas tracer 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_{\rm BH}$ measurements have already been observed by ALMA in Cycle-5. Here, we start the project with the $M_{\rm BH}$ measurement for NGC 3504, the first galaxy in the sample that has been observed.

We also want to test the growing power of ALMA and the capacity of the gas dynamical methods to (i) populate the sub-million Solar masses BH populations, (ii) constrain the scatters and slopes of BH–galaxy scaling relations at the low-mass regime, and (iii) precisely measure the occupation fraction of low-mass galaxies hosting central BHs. The occupation fraction is an important parameter helping to constrain the unknown formation mechanisms of BH seeds in the early Universe, which either form from the direct collapse of primitive gas clouds and produce massive seeds (Lodato & Natarajan, 2006; Bonoli et al., 2014) or from the remnants of the first stars (Population III) and produce lighter seeds (Volonteri et al., 2008; van Wassenhove et al., 2010; Volonteri, 2010; Volonteri & Bellovary, 2012; Volonteri, 2012; Fiacconi & Rossi, 2016, 2017).

The paper is organized into nine Sections. The properties of the galaxy NGC 3504 are presented in Section 2. In Section 3, we present the ALMA observations of ${}^{12}{\rm 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 Section 4. The mass modeling of NGC 3504 is discussed in Section 5. We model the ${}^{12}{\rm CO(2-1)}$ gas disk and estimate the central $M_{\rm BH}$ and uncertainties via Kinematics Molecular Simulation (KinMS; Davis, 2014) model and Tilted-ring model (e.g., Neumayer et al., 2007; den Brok et al., 2015) in Sections 6 and 7, respectively. We discuss our results in Section 8 and conclude in Section 9. Throughout the paper, unless otherwise indicated, all quantities quoted in this work have been corrected for a foreground extinction $A_{V}=0.306$ (Schlafly & Finkbeiner, 2011) using the interstellar extinction law of Cardelli et al. (1989).

2 The Galaxy NGC 3504

NGC 3504 is a nearby double-barred spiral LTG (Hubble type (R)SAB(s)ab) located at the distance of $D=13.6$ Mpc estimated from Tully-Fisher relation (Russell, 2002) in the Leo Minor Group (de Vaucouleurs, 1975), giving a physical scale of 68 pc arcsec*-1*. Using Sloan Digital Sky Survey (SDSS) photometric magnitude: $g=12.10$ mag, $r=11.39$ mag and Roediger & Courteau (2015) $(g-r)$–$M/L_{r}$ relation give $M/L_{r}=2.5$ ($M_{\odot}$/$L_{\odot}$), then the total stellar mass of NGC 3504 is M_{\star}\sim 6\times 10^{9}$$M_{\odot} ($M_{r}=-19.3$ mag). Laurikainen et al. (2004) decomposed the bulge mass of NGC 3504 using bugle-disk-bar decomposition model ($H$–band image of the Ohio State University Bright Galaxies Survey, OSUBGS) and found the bulge/disk ratio is 0.356. Adopting the above total stellar mass gives the bulge mass of M_{\rm Bulge}=2.1\times 10^{9}$$M_{\odot}. Recent work of Salo et al. (2015) also used images at IR wavelengths (3.6 and 4.5 $\mu$m) from the Spitzer Survey of Stellar Structure in Galaxies (S4G) imaging (Sheth et al., 2010) to decompose the bugle-disk-bar structure of NGC 3504. They found the total apparent magnitude of the bulge is $m_{\rm Bulge}=11.59$ mag222[https://www.oulu.fi/astronomy/S4G_PIPELINE4/MAIN/

decomp0087.html#entry0009](https://www.oulu.fi/astronomy/S4G_PIPELINE4/MAIN/%5C%5C%0Adecomp0087.html#entry0009). Taking into account the distribution of the disk and bar within the region of dominant bulge ($\sim$25$\arcsec$) and assuming M/L${}_{3.6\mu{\rm m}}\sim 0.7$ ($M_{\odot}$/$L_{\odot}$) (the Bell et al. (2003) $(g-r)$–$M/L_{K}$ relation gives $M/L_{K}=0.85$ ($M_{\odot}$/$L_{\odot}$)), the bulge mass is estimated as M_{\rm Bulge}=1.7\times 10^{9}$$M_{\odot}. These bulge mass estimates are both close to our bulge mass estimate using our a few first multiple Gaussian expansion (MGE; Emsellem et al., 1994; Cappellari, 2002) components model in Section 8.2 distributed at the effective radius of the bulge (Laurikainen et al., 2004).

The galaxy center is determined at (R.A., Decl.) = ($11^{\rm h}03^{\rm m}11^{\rm s}\!\!.210$, $27\arcdeg 58\arcmin 21\farcs 00$) in the $\alpha$(J2000) Equatorial coordinate system and has a systemic velocity of $1525.0\pm 2.1$ km s*-1* determined from the Fabry–Perot observations in the frame of the Gassendi $H\alpha$ survey of SPirals (GHASP; Epinat et al., 2008). Photometry shows the inclination between the line-of-sight (LOS) and the polar axis of the galaxy is 26.4*∘333http://leda.univ-lyon1.fr/ and oriented with a position angle (PA) of 150∘* (Paturel et al., 2000).

Kinematic study of NGC 3504 finds its velocity field is dominated by circular motions (Peterson, 1982). Palomar spectroscopic survey measures the stellar velocity dispersion of the bulge is $\sigma=119.3\pm 10.3$ km s*-1* (Ho et al., 2009), suggesting a central SMBH with mass of \sim$$0.9^{+3.8}_{-0.4}\times 10^{7}$$M_{\odot} for this galaxy based on the $M_{\rm BH}-\sigma$ relation (Kormendy & Ho, 2013).

Both optical and radio data suggest that the galaxy has a composite nucleus with both AGN and emission from star formation (Keel, 1984). $H\alpha$ obtained from 0.9 m telescope at Kitt Peak National Observatory (Kenney et al., 1993) shows a central starburst localized within 4$\arcsec$ (272 pc) and peaks in a ring of 1–2$\arcsec$ (68–136 pc). The star formation rate (SFR) within a diameter of 11$\arcsec$ of NGC 3504 is estimated as $2.3\pm 0.4$ $M_{\odot}$ yr*-1* (D. Nguyen et al. in preparation), confirming it is a well-known starburst (Kenney et al., 1993; Boselli et al., 2015), while the total SFR of the whole galaxy body is $0.92\pm 0.37$ $M_{\odot}$ yr*-1* (Erroz-Ferrer et al., 2015). Moreover, Palomar optical survey defines the nucleus of NGC 3504 as a transition object (Ho et al., 1993, 1997) and suggests the central starburst component is probably powered by hot O-type stars (Ho & Filippenko, 1993). NGC 3504 shows a compact nuclear unresolved source observed by very-long-baseline interferometry (VLBI) with a luminosity of $L_{\rm 1.4GHz}=0.6\times 10^{21}$ W Hz*-1* sr*-1* at 1.5 GHz (Condon et al., 1998; Deller & Middelberg, 2014).

3 Data and Data Reduction

3.1 Hubble Space Telescope (HST) Imaging

We use HST observations in wide field camera 3 (WFC3) IR band F160W to create a mass model (Section 5), which will be used as an input ingredient for dynamical models in Sections 6 and 7. The HST data was observed in 2012 May 01 (GO-12450, PI: Kochanek) with a total exposure time of 1398 s. We downloaded this image from the Hubble Legacy Archive (HLA) directly and used them throughout our analysis. However, to test the sky background level accurately, we also downloaded the flat-fielded (flt) images from the HST/The Barbara A. Mikulski Archive for Space Telescopes (MAST), combined these images in the same filter using drizzlepac/Astrodrizzle (Avila et al., 2012), and compared the combined image with the HLA image. We choose F160W filter as the default image to model the mass-follows-light map, while the F110W image is used as an alternative mass model to estimate the uncertainty caused by different wavelengths and extinction (Section 6.5.1).

Figure 1 shows the large structure of NGC 3504 in the F160W image with prominent bar connecting the bulge and galactic disk on the left and its zoom-in at the field-of-view (FOV) of $10\arcsec\times 10\arcsec$ on the right. There is a few regular dust lanes circling the galaxy center, suggesting a circumnuclear gas disk (CND), which extends to a radius of at least $5\arcsec$ from the center.

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\farcs 00$) in the $\alpha$(J2000) system as presented in HLA images. This is offset \sim$$0\farcs 21 compare to the peak of the compact ${}^{12}{\rm CO(2-1)}$ continuum emission and the optical center of the galaxy derive from optical images (Section 3.2.1); those are consistent to each other. We therefore align the HST images to the ${}^{12}{\rm CO(2-1)}$ continuum emission map to correct for this astrometric mismatch.

We used Tiny Tim point spread functions (PSFs; Jedrzejewski, 1987; Jedrzejewski et al., 1987) for the WFC3/IF F160W and F110W images to decompose the MGE model (Emsellem et al., 1994; Cappellari, 2002) in Section 5.

3.2 ALMA Observations

The observations of the ${}^{12}{\rm 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 including one 1875 MHz FDM spectral window covers over the ${}^{12}{\rm CO(2-1)}$ line and three 2 GHz TDM spectral windows added simultaneously to detect continuum emission. The observations use 46 ALMA’s 12m antennas with the C43-9 and C43-6 configurations so that the data has a maximum recoverable scale (MRS) of 10$\arcsec$ in diameter. The raw ALMA data are calibrated by the ALMA regional center staff using the standard ALMA pipeline. Flux and bandpass calibration are conducted using the quasar J1058+0133, while the atmospheric phase offsets of the data are determined as a phase calibrator using J1102+2757. More details of these observations are summarized in Table 1.

Continuum emission is detected at the galactic center of NGC 3504 only (panel-a, Figure 2) and measured over the full line-free bandwidth, fitted by a power law function, and then subtracted from the data in the $uv$-plane using the task uvcontsub of the Common Astronomy Software Applications (CASA) package version 5.1.1.

We first combine the visibility files of these two measurement sets into a final continuum-free calibrated data using the CASA task concat with the optimal combination ratios of 1/3 and 2/3 for the high- and low-spatial-resolution measurement set during the visweightscale mode, respectively. Second, we create a three-dimensional (3D) (R.A., Decl., velocity)-combined cube from the continuum-free calibrated file using the clean task. To dynamically model the CND and estimate the $M_{\rm BH}$, we image the combined cube using the channel width of 10 km s*-1* (Davis, 2014), Briggs weighting with a robust parameter of 0.5, and pixel size of $0\farcs$013 in order to optimize the sensitivity of diffuse gas and resolution. The sidelobes on the image are reduced using a mask during interactive mode. We estimate the root-mean square noise, RMS = 42 $\mu$Jy beam*-1*, in a few blank signal channels, then set 3 $\times$ RMS as the clean threshold in regions of source emission for dirty channels. We also do primary beam correction during the imaging cube. The final calibrated ${}^{12}{\rm CO(2-1)}$ cube has a synthesized beam of $0\farcs 044\times 0\farcs 031$ (3.0 pc $\times$2.1 pc), PA $\sim$ 9.7*∘, and RMS $\sim$ 0.187 mJy beam-1* km s*-1*. The nucleus molecular gas emission of NGC 3504 is detected from 1340 to 1680 km s*-1*.

3.2.1 Continuum Emission

The continuum emission is clearly resolved and centrally peaked as seen in the panel-a of Figure 2. 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\farcs 040$, $+27^{\circ}58^{\prime}20\farcs 80\pm 0\farcs 18$) with the error bars include both of the imfit fit and ALMA astrometric uncertainties. The SDSS data release 14 (Abolfathi et al., 2018) also report an optical center at ($11^{\rm h}03^{\rm m}11^{\rm s}.112$, $+27^{\circ}58^{\prime}20\farcs 77$) that is consistent to what we found here. imfit estimates the size of the emission source of $30^{+10}_{-10}\;10^{-3}$ arcsec $\times$ $10^{-3}$ arcsec (or $2.0^{+0.7}_{-0.7}$ pc $\times$ pc), with a PA of $\sim$30*∘* and the total integrated intensity of $4.70\pm 0.16\pm 0.24$ mJy over the emissions free of USB (442–446 GHz) and LSB (227.5–231.5 GHz) frequency windows. Note that the former/latter errors are the systematic error/ALMA flux calibration error, respectively.

3.2.2 12CO(2–1) Line emission

We create the integrated intensity, intensity-weighted mean velocity field, and intensity-weighted velocity dispersion maps for the nuclear ${}^{12}{\rm CO(2-1)}$ gas disk from the combined cube directly in Figure 2. The emission is significant within the radius of 5$\arcsec$. We enhance the quality of these maps using the moments masking technique; details of this technique were described in Davis et al. (2017); Onishi et al. (2017); Davis et al. (2018).

The integrated intensity map reveals the presence of a nuclear rotating disk. However, as seen in the panel-b and c of Figure 2, this disk seems separating into two distinct rotating disks including a dense compact CND distributed within 1$\arcsec$ and a faint extended disk of spiral arms and voids. Kinemetry444http://davor.krajnovic.org/idl/#kinemetry (Krajnović et al., 2006) examination reveals these two disks are clearly misalignment by an angle of $\sim$30*∘*, this is also consistent to what found in Knapen et al. (2002).

The velocity field (panel-c) shows a regular disk-like rotation with a total velocity width of \sim$$300 km s*-1*. The velocity dispersion map (panel-d) are quite flat with constant values of $\sim$30 km s*-1* and $\sim$10 km s*-1* in the inner and outer gas disks, respectively. However, there is a suddenly increasing peak $\sim$65 km s*-1* at the galactic center. This is not an intrinsic velocity dispersion but beam smearing over the velocity gradient at the galaxy center and the LOS integration through a nearly edge-on orientation of the inner gas disk ($i\gtrsim 60^{\circ}$).

Figure 3 shows the integrated ${}^{12}{\rm CO(2-1)}$ spectrum of NGC 3504 that has a classical double-horn shape of a rotating disk. While in Figure 4, we plot the position-velocity diagram (PVD) extracting from a cut through the major axis of the gas disk (PA $\sim 332^{\circ}$). We show both PVDs of the combined cube (high resolution) and the low-resolution cube side by side in this figure to demonstrate the increased rotation towards the center even at different spatial scales. We interpret this motion as a Keplerian curve caused either by an SMBH existing at the heart of NGC 3504 or by non-circular motions of an inner ring/spiral within the innermost regions, which is typically seen in barred galaxies.

The intensity map also shows a central hole, which is marginally resolved on pc scales as seen in the zoom-in ${}^{12}{\rm CO(2-1)}$ intensity map (panel-a, Figure 5). The intensity remains above zero inside the hole is due to either the angular extent of the synthesized beam or some galactic foreground gas distribute further above the torus and align along the LOS. Another possibility of the non-zero intensity in the hole is a genuine existence of ${}^{12}{\rm CO(2-1)}$ in the hole but faint due to high density, making CO gas is excited to higher-J transitions. Based on its physical radius of $\sim 0\farcs 04$ (or 2.7 pc, estimated in Section 6) and centralization, we believe that we are marginally observing a dusty molecular torus surrounding an accreting SMBH, a fuel source of the unified AGN model (Antonucci, 1993). We estimate the size of this dusty molecular torus using MIR flux observations (Rieke & Low, 1972) and the torus size-luminosity relation from Tristram & Schartmann (2011). The radius of the hole is $\sim$2 pc, roughly corresponding to the hole in our ${}^{12}{\rm CO(2-1)}$ observations. This radius is bigger than the inner dust sublimation radius of the torus, which is typically of $<$1 pc (Barvainis, 1987). The reason for this mismatch is the ${}^{12}{\rm CO(2-1)}$ emission is sometimes deficient in the nuclear regions, and therefore may not be the best tracer for the torus (Imanishi et al., 2018; Izumi et al., 2018).

4 Detection of Dense Gas CS(5–4)

We detect a dense gas tracer ${\rm CS(5-4)}$ in one of the continuum spectral windows. This ${\rm CS(5-4)}$ is centrally concentrated, filling in the hole of the ${}^{12}{\rm CO(2-1)}$ map. The integrated intensity map, spectrum, and radial profile of ${\rm CS(5-4)}$ are shown in the panel-b, c, and d of Figure 5. The detection of ${\rm CS(5-4)}$ is significant above 20$\sigma$ ($\sigma\sim 0.3$ mJy beam*-1* km s*-1*) but in a very low-velocity-resolution spectral window with only 128 channels over 2 GHz band width ($\sim$40 km s*-1*). We estimate the total flux is $1.76\pm 0.42$ Jy km s*-1* with 10% of the error budget comes from the flux calibration uncertainty of ALMA data. The ${\rm CS(5-4)}$ line is more centrally concentrated than the ${}^{12}{\rm CO(2-1)}$ emission as shown in the moment 1 and moment 2 maps and the PVD in Figure 6. These features suggest the ${\rm CS(5-4)}$ line is an alternative transition that could provide a better constraint on the $M_{\rm BH}$ than the ${}^{12}{\rm CO(2-1)}$ line. This is because the central concentration of ${\rm CS(5-4)}$ would recover the high-velocity upturn in the data at the very center that is missing in the current ${}^{12}{\rm CO(2-1)}$ map due to the central hole. However, the velocity resolution of ${\rm CS(5-4)}$ is not good enough to perform such dynamical models. A higher-velocity-resolution observation of ${\rm CS(5-4)}$ with $\sim$30 times improvement in velocity resolution is required to reduce the uncertainty on $M_{\rm BH}$ determination.

5 Creating a Mass Model

In galaxies the stellar M/L 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 on dynamical $M_{\rm BH}$ estimates (Nguyen et al., 2017, henceforth N17; N19). Unlike our work in N17 and N19, we are lacking the nuclear stellar spectroscopic information for NGC 3504. However, we examine the nuclear $J-H$ color map of NGC 3504. Here, we assume HST/WFC3 F110W $\approx$ $J$ and F160W $\approx$ $H$. We find a constant color of $J-H\sim 0.85$ mag across the FOV where we measure the ${}^{12}{\rm CO(2-1)}$ kinematics, although there is variability in few central dust lances and at larger scale of the galaxy; those will be masked out. Photometry examination with Jacobus Kapteyn Telescope (JKT) images also found a constant $B-I\sim 1.3$ across the FOV of 5$\arcsec$ (Knapen et al., 2002). The effect of uniform $J-H$ and $B-I$ colors in the region of interest suggests we can use a constant M/L to create the mass model for NGC 3504. We utilize the MGE code555https://www-astro.physics.ox.ac.uk/~mxc/software/#mge (Emsellem et al., 1994; Cappellari, 2002) to decompose the light surface density into individuals Gaussian components that can then be deprojected. We use the mge_fit_sectors IDL version 4.14 (Cappellari, 2002) to fit the F160W image and deconvolve the effects of a PSF.

We first parameterize the PSF using Gaussian functions in the first MGE fit, 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 map that masks out pixels contaminated by prominent dust lanes, a small-jet like feature, and bright stars near the center, as well as the bar that connects the central bulge and the larger structure along the major axis seen in the left panel of Figure 1. This bar causes a twist on the MGE model if we allow the PA changing as a free parameter, which cannot be modeled in an axisymmetric dynamical models. In this work, we fix the PA $=-28^{\circ}$ during the MGE fit based on the orientation of the whole galactic disk for an axisymmetric model.

Each Gaussian can be deprojected analytically with a specific axis ratio (or inclination) to reconstruct a 3D light distribution model. To optimize data vs. MGE model, we set the axis ratio in the range of $q=0.65-1.0$. This MGE parameters are listed in Table 3. Note that we assume HST/WFC3 F160W solar Vega absolute magnitude system666http://mips.as.arizona.edu/~cnaw/sun.html.

We show the radial light surface brightness density and its best-fitting MGE model in the upper-left panel of Figure 7, while the fractional residual indicates the agreement between the best-fit model and the data are shown in the upper right panel. The 2D light surface density is also plotted with their MGE model in the bottom panels of Figure 7 with the whole image on the left and the zoom-in $30\arcsec\times 30\arcsec$ FOV on the right. The figure shows the agreement of the data and its MGE model at the same radii and contours levels. The difference between the best-fit model and the data are $<$15% across the FOV ($>$100$\arcsec$).

6 KinMS Dynamical Modeling

In this section, we describe the KINematic Molecular Simulation dynamical model (KinMS777https://github.com/TimothyADavis/KinMS; Davis, 2014) that we employ to measure the $M_{\rm BH}$ in NGC 3504 and state our results.

The KinMS model is a mm-wave observational simulation tool developed by Davis et al. (2013) using the Markov Chain Monte Carlo (MCMC) method to simulate dynamical motion of molecular/atomic cold gas distributions under the influence of galaxy and central dark massive object’s gravity. In practice, the MCMC technique allows to plug in initial guesses for the true gas distribution and kinematics based on the assumption that the gas is rotating in circular orbits. The model creates a simulated data cube, which can be compared to the observed data via $\chi^{2}$-minimized likelihood function (Davis et al., 2017; Onishi et al., 2017; Davis et al., 2018). The model explores the parameter space using Bayesian analysis technique with a set of walkers using the emcee algorithm (Foreman-Mackey et al., 2013) and affine-invariant ensemble sampler (Goodman & Weare, 2010). The relative likelihood for each walker at each step will determine for their next move through parameter space. The best-fit model parameters are then obtained from the posterior distribution of the full pool of model parameters. In this work, we use the Python code KINMSpy_MCMC888https://github.com/TimothyADavis/KinMSpy_MCMC to find the best set of model parameters.

6.1 Nuclear Gas Morphology

The KinMS model requires surface brightness distribution function of gas, $\Sigma(r)$, that describes the gas morphology radially. We assume an axisymmetric morphology for the nuclear ${}^{12}{\rm CO(2-1)}$ gas in NGC 3504 as a thin disk distributed continuously out to \sim$$5\arcsec.

We examine the morphological distribution of the nuclear ${}^{12}{\rm CO(2-1)}$ gas by extracting the flux of the low-resolution cube in a central box with seven pixels in width along the major axis and plotting the slit within the area of $1\arcsec\times 1\arcsec$ where we actually fit for the KinMS model in Figure 8. Note that we limit our fitting region within $1\arcsec\times 1\arcsec$ to avoid the contamination from the outer rotating disk. We model the inner gas disk morphology using a (1) a central offset single-Gaussian function to describe the ring morphology and central attenuated hole of the ${}^{12}{\rm CO(2-1)}$ integrated intensity map and (2) an exponential disk, which is fitted quite well to the outer disk of the nuclear gas reservoirs. Overall, we describe the nuclear gas distribution by three parameters including the exponential disk length scale ($R_{e}$), the Gaussian FWHM ($G_{h}$), and center peaks ($G_{\rm c}$).

The best-fit profile of the nuclear gas distribution produced by the KinMS model is plotted as a red solid line in Figure 8. We should note that the KinMS model matches the observations by fitting a set of free parameters including the above gas morphology and the total flux (scaling factor), PA, and $i$ of the gas disk, as well as its kinematic center in R.A. ($x_{c}$), Decl. ($y_{c}$), and velocity offset ($v_{\rm off.}$), M/L, $M_{\rm BH}$, gas disk thickness ($d_{\rm t}$), and internal velocity dispersion of the gas ($\sigma$). This dispersion is assumed as a spatial-constant parameter. Due to the axisymmetric model, we constraint $i$ and PA in each fit as a single value throughout the disk to avoid warp. So, our initial KinMS model thus has 13 free parameters listed in Table 4.

6.2 Velocity Model

We assume the gas rotates around the galaxy center in circular orbits and that vary radially. The KinMS model takes in the circular velocity curve in the equatorial plane of the CND to produce a simulated cube. This circular motion of the gas is controlled by the gravitational potential of the galaxy calculated from the stellar mass model (Section 5) and the point source potential of the SMBH described in Cappellari (2002). We calculate this circular velocity profile using the mge_vcirc procedure within the Python Jeans Axisymmetric Modelling (JamPy999http://www-astro.physics.ox.ac.uk/~mxc/software/#jam; Cappellari, 2008) package.

6.3 Fitting Process

We first run the KinMS model with an area of 300 pixels $\times$ 300 pixels ($4\arcsec\times 4\arcsec$), where roughly covers most of important features of the molecular gas disk, and a number of iterations of $10^{5}$. To avoid the possible effect of AGN on our dynamical $M_{\rm BH}$ and M/L estimates simultaneously, we leave out the marginalized resolved Gaussian component during the model fitting; the first MGE component denoted by an upper star in Table 3. From this run, we obtain rough estimates of the model parameters. The second fit then fits in the central 80 pixels $\times$ 80 pixels ($1\arcsec\times 1\arcsec$ or 68 pc $\times$ 68 pc) area only, and starts with reasonable flat priors for all free parameters to ensure our kinematic fitting process converge. These priors are the best fit from the first fit. We note that the prior on the $M_{\rm BH}$ was flat in log-space, while the inclination of the gas disk was allowed to vary over the full physical range allowed by the MGE model. The search ranges of these parameters are shown in column 2 of Table 4, and good fits were always found well within these ranges. To account the covariance between pixels, which describes the neighbouring spaxels strongly correlated by the synthesized beam due to the Nyquist spatial sample of ALMA data (Davis et al., 2017; Onishi et al., 2017; Davis et al., 2018), we increase the RMS by a factor $(2N\times M)^{1/4}$, where $N$ and $M$ are the numbers of spaxels and velocity channels included in the fit (Mitzkus et al., 2017).

We take in to account the distributions of any interstellar material including gas and dust within the fitting region into our mass model. The total flux in this region is $6.9_{\pm 0.7}^{\pm 0.2}$ Jy km s*-1*, where the upper and lower uncertainty are systematic and $\sim$10% flux calibration uncertainty in the measurement and ALMA data, respectively. This total flux provides an estimate on the total mass of H2 gas of M_{\rm H_{2}}=(4.2\pm 0.5)\times 10^{7}$$M_{\odot} by assuming the line ratio ${}^{12}{\rm CO(2-1)}$/${}^{12}{\rm CO(1-0)}=0.8$ (Bigiel et al., 2008) and H2-to-CO conversion factor for starburst galaxies: $N({\rm H_{2}})/I_{(\rm 1-0)}=X_{\rm CO}=0.5\times X_{\rm CO}^{\rm MilkyWay}=1\times 10^{20}$ cm*-1* (K km s*-1*)-1 for converting ${}^{12}{\rm CO(1-0)}$ to H2 mass (Kuno et al., 2000, 2007; Bolatto et al., 2013). The dust mass of (5.1\pm 0.7)\times 10^{5}$$M_{\odot} is calculated in D. Nguyen et al. in preparation. In Figure 9, we plot the 1D stellar mass and gas + dust mass profiles as the radial functions simultaneously within the radius of $6\arcsec$. In our fitting area, the gas and dust mass distributions are dominant over the stellar mass, suggesting the gravitational effect of these components plays an important role in determining the $M_{\rm BH}$ in NGC 3504 accurately. We add these masses in the KinMS gasGrav mechanism, assuming the dust and gas is co-spatial distribution.

We run the model fits with velocity channels are from $-$200 to 200 km s*-1*, w.r.t. the systemic velocity of 1525 km s*-1*. The total number of iterations are set to be $3\times 10^{6}$ and the first 20% of iterations are considered as the burn-in phases to produce our final posterior probability distributions of these 13 model parameters.

6.4 Model Results

We clearly detect an SMBH that causes the increasing Keplerian curve towards the center. The best-fit model parameters are identified directly from our Bayesian analysis, relying on the likelihood probability distribution functions (PDFs) generated via MCMC. We choose the best-fit parameters that are the medians of the parameter posterior PDFs. Particularly, the probability is marginalized over to produce a best-fit value by taking the median of the marginalized posterior samples for each parameter. The 1$\sigma$, 2$\sigma$, and 3$\sigma$ confidence levels (CLs) or uncertainties are estimated from all models within (16% and 84%), (3% and 97%), and (0.2% and 99.8%) of the PDFs, respectively. At $3\sigma$ CL, the central $M_{\rm BH}$ is measured at M_{\rm BH}=1.02^{+0.18}_{-0.15}\times 10^{7}$$M_{\odot}, while the stellar M/L${}_{\rm H}=0.66^{+1.44}_{-0.65}$ ($M_{\odot}$/$L_{\odot}$) and i=61.01^{\circ}$${}^{+3.11}_{-3.51}. The reduced chi-square ($\chi^{2}_{\rm red}$) for the best fit is 1.16, indicating a good fit. We list all the best-fitting parameters and their likelihoods in Table 4.

Figure 10 shows the 1D and 2D marginalization of the physical parameters included in the fit. All correlated parameters are well constrained by the data. We show the observed PVD overlaid with the PVD of the best-fit model in Figure 11. To enhance the illustration of our result and how good the best-fit model fit to the data, we add the models with no SMBH and with an overly large $M_{\rm BH}$. The best-fit model with an appropriate $M_{\rm BH}$ clearly reproduces the kinematics of the molecular gas better than the other two. The velocity map of the data, the velocity map of the KinMS best-fit model, and the velocity-residual map ${\tt(Data-Model)}$ are shown in the bottom-row plots of Figure 11.

We find evidence of (1) some non-gravitational motions resulted in high residual on the southwest side of the nucleus as a signature of either outflows or caused by a warped inner disk/bar and (2) signature of a nuclear small rotating structure at the scale of $r\sim 0\farcs 15$. This central structure is oriented along the minor axis. However, its morphology is very ambiguous. The detailed investigations of these nuclear rotating structure and non-circular motions will be discussed in D. Nguyen et al. in preparation.

In Figure 14 we show our best-fit model overlaid on the channel maps of our ALMA observations. Also, we demonstrate the same results for the low-resolution data cube in Figure 15 of the Appendix. The full list of the best-fit parameters and their likelihoods are presented in Table 6. We note that the central ${}^{12}{\rm CO(1-0)}$ hole does not appear in the low-resolution cube because its size is smaller than the cube’s beamsize. We model the central ${}^{12}{\rm CO(1-0)}$ surface brightness profile within the fitting area by a single Gaussian with the peak tighten to the kinematic center ($G_{c}=0$), leaving the KinMS model with 11 free parameters. From the result of this fit, we do not find any signatures of the inner rotating structure in the velocity residual map as seen in the case of high-resolution cube fit. This is because the spatial scale of this structure is smaller than the spatial angular resolution of this cube.

With the stellar velocity dispersion of the bulge of $\sigma=119.3\pm 10.3$ km s*-1* (Ho et al., 2009), the BH in NGC 3504 has an intrinsic SOI radius $R_{\rm SOI}=GM_{\rm BH}/\sigma^{2}=4.6$ pc ($0\farcs 067$), where $G$ is the gravitational constant, $M_{\rm BH}$ is the $M_{\rm BH}$, and $\sigma$ is the bugle stellar velocity dispersion. However, Davis (2014) argue the spatially angular resolution (in pc scale) that is required to detect a BH with specific $M_{\rm BH}$ and $\alpha$ CL in an object with a circular velocity caused by luminous matter in their equation (9). Using this formula for our NGC 3504 observations and expecting to detect the given $M_{\rm BH}$ at $3\sigma$ CL, we find the essential angular resolution of $\sim$3.1 pc ($0\farcs 045$). This scale corresponds to the FWHM of the synthesized beamsize of the combined cube. All of these scales mean that our ALMA observations are able to constrain $M_{\rm BH}$ and kinematic properties of the CND in NGC 3504.

6.5 Uncertainty Sources on the $M_{\rm BH}$ Estimate

6.5.1 Stellar Mass Models

We test the robustness of the mass model by create a new MGE mass model using the HST/WFC3 F110W image (approximate to $J$-band), which gives the best-fit $M_{\rm BH}$ of 8.05_{-0.21}^{+0.27}\times 10^{6}$$M_{\odot} and M/L${}_{J}=0.76_{-0.61}^{+1.17}$ ($M_{\odot}$/$L_{\odot}$). Moreover, the way we fit the light MGE model from HST images also causes some error on the mass model. Since large structure of NGC 3504 hosts an outer Lindblard resonance (OLR; Buta & Crocker, 1992; Buta & Combes, 1996), while its smaller scale contains central elongated bar (Kuno et al., 2000). We rerun two MGE fits as follows: (1) there is no constrain on $q$ to get a MGE with elongated barred Gaussian dominant, and (2) set $q=0.9-1$ to get a MGE with OLR Gaussian dominant. We have double checked these MGEs within $16\arcsec\times 16\arcsec$ central region do not change significantly. The KinMS model with the former MGE gives the best fit ($M_{\rm BH}$, $\emph{M/L}_{H}$, $i$) = (8.74_{-0.08}^{+0.09}\times 10^{6}$$M_{\odot}, $0.72_{-0.60}^{+1.26}$ ($M_{\odot}$/$L_{\odot}$), 58.34∘${}_{-3.16}^{+3.21}$), while that of the latter MGE gives the best fit ($M_{\rm BH}$, $\emph{M/L}_{H}$, $i$) = (1.16_{-0.21}^{+0.26}\times 10^{7}$$M_{\odot}, $0.61_{-0.58}^{+1.36}$ ($M_{\odot}$/$L_{\odot}$), 60.45∘${}_{-3.76}^{+3.52}$); other nuisance parameters vary within 15% compare to the best fit of the default model. We thus conclude that our $M_{\rm BH}$ estimate is robust to the systematic errors of our mass models.

6.5.2 AGN and Distance

The contribution of AGN at the center of NGC 3504 as seen in optical (Ho & Filippenko, 1993) and radio (Deller & Middelberg, 2014) may attribute to the F160W waveband flux at the very center, which is represented by the first MGE component in Table 3. We rerun the KinMS model included this component to test for the possible effect of AGN distribution, providing $M_{\rm BH}$ =1.30_{-0.05}^{+0.09}\times 10^{7}$$M_{\odot} and M/L**H $=0.58_{-0.51}^{+1.26}$ ($M_{\odot}$/$L_{\odot}$). This systemic uncertainty of the $M_{\rm BH}$ is within our 3$\sigma$ CL, suggesting the impact of AGN on our mass model is insignificant.

In addition, the $M_{\rm BH}$ estimate is systematically affected by the distance to the galaxy based on the relation $M_{\rm BH}\propto D$. In this work, we assumed the Tully-Fisher distance to NGC 3504 is $D=13.6\pm 1.4$ Mpc (Russell, 2002). Therefore, the systematic distance uncertainties on the $M_{\rm BH}$ are of 10% as similar as the random uncertainties. We should note that the six existing distance estimates101010https://ned.ipac.caltech.edu/ to NGC 3504 are in a wide range of 8.7–26.5 Mpc, suggesting the systemic errors due to different distances adopted are much larger than that of the distance estimate which we are using (Russell, 2002) up to $\sim$95%.

6.5.3 Non-circular Motions

In this work, we assume the gas in NGC 3504 is in purely circular motion. However, non-circular motions (e.g., inflow, outflow, streaming) are usually seen in barred galaxies like NGC 3504, where the non-circular motion is primarily caused by gas streaming along the bar and could affect our analysis significantly. The large scale bar ($>$20$\arcsec$) does not strongly affect to our dynamical results because our interested in region is deeply small within the bulge and the nuclear molecular gas disk resides within the radius of $\lesssim$5$\arcsec$.

However, from the residual map of Figure 11 we find the orientation of the nuclear small rotating structure, which is perpendicular to the projected major axis. Randriamampandry et al. (2015) shows the non-circular motions derived from mass profiles of strongly barred galaxies vary dramatically if the bar is orientated at specific angles w.r.t. the LOS. Particularly, there are under/overestimates the circular motions when the bar is parallel/perpendicular to the projected major axis. As a SAB(s)ab type galaxy like NGC 3504, the expected overestimate rotation is minimum as suggested by Figure 2 of Randriamampandry et al. (2015, top-row plots). This is because the contribution of gas from the bulge is compensated for gas streaming motion along the bars.

To calculate the effect of this non-circular motion on our dynamical results, we run the KinMS model with a new velocity model extracted along the major axis of the moment 1 map (Figure 2) using the IDL Kinemetry code (Krajnović et al., 2006) instead of the velocity model built up from the mass MGE model in Section 6.2. The best-fit model gives $M_{\rm BH}$ of 1.12_{-0.08}^{+0.07}\times 10^{7}$$M_{\odot} and M/L${}_{H}=0.64_{-0.56}^{+1.38}$ ($M_{\odot}$/$L_{\odot}$), while all other parameters are variable within 9% compare to the best-fit value listed in Table 4.

6.5.4 Gas Velocity Dispersion

During the above analysis, we assumed a constant gas velocity dispersion. In reality, the velocity dispersion could vary with radius and azimuth within the gas disk. In the central part of the galaxy, where beam smearing is important, an increase of velocity dispersion could lead to overestimated $M_{\rm BH}$. To quantify this effect, 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)_{\rm gas}$:

(a) Linear gradient: $\sigma(r)_{\rm gas}=a\times r+b$, where $a$ and $b$ are free parameters. We find $a\approx 0$ and $b=37.57$ km s*-1*, and KinMS results are consistent with the default model of constant velocity dispersion.

(b) Exponential: $\sigma(r)_{\rm gas}=\sigma_{0}\exp{(-r/r_{0})}+\sigma_{1}$, where $\sigma_{0}$, $r_{0}$, and $\sigma_{1}$ are free parameters. As discussed in Barth et al. (2016a), we set the lower boundary for $\sigma(r)_{\rm gas,\;min}=1$ km s*-1* during the fit to prevent the line-profile widths becoming arbitrarily small. The best-fit KinMS model provides M_{\rm BH}=9.12_{-0.08}^{+0.10}\times 10^{6}$$M_{\odot} and M/L${}_{H}=0.70_{+1.37}^{-0.62}$ ($M_{\odot}$/$L_{\odot}$) with an exponential dispersion model with $\sigma_{0}=36.95_{-2.55}^{+1.40}$ km s*-1*, $r_{0}=1\farcs 41_{-1.13}^{+0.30}$, and $\sigma_{1}=9.89_{-1.27}^{+1.31}$ km s*-1*. Other best-fit parameter values are consistent with the default model.

(c) Gaussian: $\sigma(r)_{\rm gas}=\sigma_{0}\exp{(-(r-r_{0})^{2}/2\mu^{2})}+\sigma_{1}$, where $\sigma_{0}$, $r_{0}$, $\mu$, and $\sigma_{1}$ are free parameters. We allow the parameter $r_{0}$ 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)_{\rm gas,\;min}=1$ km s*-1* during the fit. The best-fit KinMS model provides M_{\rm BH}=8.73_{-0.05}^{+0.09}\times 10^{6}$$M_{\odot} and M/L${}_{H}=0.71_{-0.56}^{+1.28}$ ($M_{\odot}$/$L_{\odot}$) with a Gaussian dispersion model with $\sigma_{0}=36.96_{-3.10}^{+2.11}$ km s*-1*, $r_{0}=0\farcs 1_{-0.15}^{+0.12}$, $\mu=0\farcs 74_{-0.10}^{+0.12}$, and $\sigma_{1}=10.01_{-1.32}^{+1.28}$ km s*-1*. Other best-fit parameter values are similar to the default model.

The minimum $\chi^{2}_{\rm red}$ are determined at 0.91, 0.87, and 0.83 for the linear, exponential, and Gaussian dispersion profile, respectively.

6.5.5 Different Observational Scales $\&$ Synthesis Beamsizes

The dynamical model of $M_{\rm BH}$/M/L**H estimated from the combined cube is higher/lower than that from the low-resolution cube by 32%/16%, respectively. The mass estimate within the synthesis beamsize is also another systemic uncertainty on the $M_{\rm BH}$. We convert the total flux in the beam of the combined cube into molecular mass of 2.5\times 10^{5}$$M_{\odot}. This mass is only 10% of the 3$\sigma$ $M_{\rm BH}$ uncertainty provided by the KinMS model (Table 4), suggesting the mass uncertainty due to the synthesis beamsize of our observation is minimal.

7 Thin Disk Tilted-Ring Dynamical Model

We also constrain the M/L, $i$, and $M_{\rm BH}$ of the central SMBH in NGC 3504 independently using a different dynamical model, which is the so called thin disk tilted-ring model (Begeman, 1987; Quillen et al., 1992; Nicholson et al., 1992; Neumayer et al., 2007; den Brok et al., 2015). We use the same kinematic measurements from the nuclear ${}^{12}{\rm CO(2-1)}$ CND and stellar mass model that are conducted in Sections 3.2.2 and 5, respectively. The purpose of this test is to examine the robustness of the measured $M_{\rm BH}$ in various assumptions and dynamical models.

7.1 The Tilted-Ring Models

We model the kinematics of ${}^{12}{\rm CO(2-1)}$ CND with tilted ring models in a similar approach as for the H21-0 S(1) transition kinematics in Seth et al. (2010); den Brok et al. (2015); N17. The basic idea of these models is that we assume the emitting gas (e.g., CO or H2) is rotating in thin rings on concentric circular orbits around the center of NGC 3504 with a velocity that can be 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 $\theta$ (relative to the projected major axis), which is projected along the LOS allowing the rings become ellipses.

To test for systemics, we fit three different models to the data as follow:

The first model mimics the presence of two rotating disks around the center. Both the inner disk and the outer disk have their own PA and $i$, which are free parameters in the fit, as well as the radius at which the model transitions from the inner to the outer disk. 2. 2.

The second model, we assume that the ellipses change their geometry smoothly with radius. We determine the radially varying PA for this model with the Kinemetry routine (see footnote 6; Krajnović et al., 2006) but allow $i$ to change linearly with radius. 3. 3.

The third model assumes the potential is axially symmetric around the $z$ axis and that the gas is rotating in the $z=0$ plane of the potential.

We assume for the flattening of the mass distribution the standard MGE deprojection. On each ring at a certain radius, the velocity of the gas is determined by the derivative of the potential with respect to radius (Seth et al., 2010; den Brok et al., 2015, and N17) as the following equation:

[TABLE]

which for the spherically symmetric potential is equal to

[TABLE]

The gravitational potential of each model is based on the $M_{\rm BH}$ as a point source and the stellar-mass-component distribution, $M_{\star}(r)$, which is modeled as the deconvolved and deprojected stellar light profile, $L_{\star}(r)$, multiplied with an additional scaling in M/L.

During the fit, our model generates a spectral cube with the same spectral sampling as the ALMA data cube. The model distributes the flux over the cube based on its spatial distribution and the velocity and velocity dispersion for each ring across the coplanar disk to replicate 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. In order to compare to the KinMS model, we optimize the tilted-ring model in the same fitting area. We ignore gas and dust distribution during the fit in this work; however, in the future work, we will develop the code to add in these masses in another form of MGE separately.

7.2 Results

We summarize the best-fit results of these three tilted-ring models in Table 5. All uncertainties of the best-fit parameters we quote in this section are determined within 3$\sigma$ CL. The tilted-ring model derives a best-fit $M_{\rm BH}$ of M_{\rm BH}=7.7_{-1.2}^{+1.5}\times 10^{6}$$M_{\odot}. This mass is 23% lower than the best-fit $M_{\rm BH}$ found by the KinMS model, while the outer disk best-fit inclination is determined at $i\sim 60^{\circ}$. However, the best-fit M/L**H is found at $0.63^{+0.08}_{-0.02}$ ($M_{\odot}$/$L_{\odot}$); less than $\sim$5% of the KinMS’s predictions. In these models, we assume uniform errors for the velocities, the confidence intervals are therefore likely underestimated.

In Figure 12 we plot the best-fit velocity map of of ${}^{12}{\rm CO(2-1)}$ gas disk, its data, and the residual ${\tt(Data-Model)}$ of the Model 2 only. As similar as the KinMS model, the tilted-ring model finds evidence of some non-circular motions on the velocity residual map as well. However, the signature of the nuclear small rotating structure is weak due to the presence of other residuals at the same scale of the structure.

8 Discussions

8.1 $M_{\rm BH}$ Measurements in Nearby Low-mass Galaxies with ALMA

ALMA observations provide a new promising path for gas-dynamical $M_{\rm BH}$ measurements and exploration of BH demographics in a variety types of Hubble sequence in nearby low-mass galaxies; especially for LTGs where host a large fraction of gas and dust in their nuclear regions. As discuss in Barth et al. (2016a), to measure $M_{\rm BH}$ in these galaxies accurately, cold gas observations are required to (i) have a simple disk-like rotation and kinematically dominated, (ii) the beamsize should be at least twice of the BH SOI along the minor axis, and (iii) the gas kinematics have to mapped well supported by high surface brightness molecular line emissions. In addition, accurate modeling of host galaxy luminosity/mass profiles using NIR observations to avoid dust contamination and extinction, with angular resolutions that are as high as the ALMA data are also greatly important.

In the near future, a number of galaxies with detection of gas high-velocity rotation within SOI should increase rapidly. Measuring $M_{\rm BH}$ in such galaxies accurately are important to anchor the local BH demographics and BH–host galaxy correlations (Maoz et al., 1998) and (2) constraint the BH seed formations in the early Universe (e.g., Volonteri et al., 2008; Bonoli et al., 2014; Volonteri et al., 2015; Fiacconi & Rossi, 2016, 2017). So, the future growing number of high-resolution observations of ALMA of nearby targets will make it possible to examine molecular gas kinematics in galaxy nuclei in far greater details than any previous surveys.

8.2 $M_{\rm BH}$ Scaling Relations

We examine our best-fit $M_{\rm BH}$ of NGC 3504 in the context of scaling relations of $M_{\rm BH}-M_{\rm Bulge}$ and $M_{\rm BH}-\sigma$ including the empirical compilations (Kormendy & Ho, 2013; McConnell et al., 2013; Saglia et al., 2016) and theoretical prediction that assumes a bimodality in BH accretion efficiency (Pacucci et al., 2017, 2018; Pacucci & Loeb, 2018) in Figure 13. The stellar velocity dispersion of the bulge of NGC 3504 is determined from Ho et al. (2009). Here, we estimate the bulge mass of NGC 3504 using our $H$-band MGE model (Section 5) and adapt the bulge effective radius from the bugle-disk-bar decomposition model (Laurikainen et al., 2004). After calibrating for the distance to the galaxy and accounting for the dynamical M/L**H (Section 6.4), we obtain M_{\rm Bulge}=(2.3\pm 0.4)\times 10^{9}$$M_{\odot}. Our result shows that the best-fit $M_{\rm BH}$ of NGC 3504 is fully consistent with both empirical and theoretical $M_{\rm BH}-\sigma$ and $M_{\rm BH}-M_{\rm Bulge}$ scaling relations of Kormendy & Ho (2013); Saglia et al. (2016); van den Bosch et al. (2016), but outside $+1\sigma$ uncertainty of the theoretical $M_{\rm BH}-M_{\rm Bulge}$ relation of Pacucci et al. (2018). At the mass of \sim$$10^{7}$$M_{\odot} the central BH of NGC 3504 lies in the middle with one order of magnitude above/below the recent dynamical measurements in low-mass systems (den Brok et al., 2015, N17; N18; N19) and in more massive targets using molecular gas (Davis et al., 2013; Onishi et al., 2015; Barth et al., 2016a, b; Davis et al., 2017; Onishi et al., 2017; Davis et al., 2018)

8.3 Complex Structure of the Molecular Gas

Our both low- and high-resolution (combined) kinematics show a larger rotational disk as seen in the residual maps of velocity fields after subtracting the data to the models when we fit the models to the data with a larger FOV ($1\arcsec\times 1\arcsec$). This larger disk is co-rotating but misaligned by 30∘ with the CND as mentioned in Section 3.2.2. So, both morphology and kinematics suggest a nuclear complex structure of the ${}^{12}{\rm CO(2-1)}$ gas in NGC 3504 with (i) a large faint disk with a size of $\pi\times 5\arcsec^{2}$ consists of spiral arms and voids, which may be associated with bars, (ii) a dominant integrated flux CND within $0\farcs 5$, and (iii) a nuclear rotating structure with a size of $r\sim 0\farcs 15$ that is revealed by the high-resolution kinematic measurement only.

The higher rotational velocities of the outermost disk and CND compare to that of the nuclear rotating structure suggest that this structure is not only oriented along the minor axis but also along the LOS direction, while the larger discs are oriented perpendicular to the LOS direction. Detailed investigations of morphology, kinematics, and role of this nuclear structure on the growth of SMBH will be discussed in D. Nguyen et al. in preparation.

9 Conclusions

We present a dynamical mass measurement for the SMBH in NGC 3504 using ${}^{12}{\rm CO(2-1)}$ emission observed by ALMA and HST imaging. Our main results are highlighted as follow:

NGC 3504 hosts a CND, which is co-spatial distribution with a obscuring dust lanes visible in the HST imaging within $5\arcsec$. 2. 2.

Both KinMS and tilted-ring models suggest the detection of a central SMBH with M_{\rm BH}\sim 1\times 10^{7}$$M_{\odot}, which is consistent with both $M_{\rm BH}-M_{\rm Bulge}$ and $M_{\rm BH}-\sigma$ relations. The agreement between two models suggests they are both powerful tools to probe the BH–galaxy scaling relations. 3. 3.

Our dynamical models also give M/L${}_{H}=0.66$ ($M_{\odot}$/$L_{\odot}$) for the nucleus of NGC 3504, suggesting its nucleus is in an early phase of its evolution. 4. 4.

Both combined (high-) and low-resolution observations give consistent constraints on the $M_{\rm BH}$ and M/L**H. These prove our observational strategy for measuring $M_{\rm BH}$ in nearby low-mass galaxies to explore the BH demographics in a large sample, which cannot achieve with the stellar dynamical method. 5. 5.

The CND of NGC 3504 has a relatively high velocity dispersion with $\sim$30 km s*-1* in the region $<$1$\arcsec$ and $\sim$10 km s*-1* at larger radii. The central high value ($>$60 km s*-1*) of dispersion is not an intrinsic but a beam smearing effect along the LOS integration through the nearly edge-on orientation of the CND. 6. 6.

Our multiple scale observations reveal the complexity of the nuclear molecular gas disk in NGC 3504 comprised of a small nuclear rotating structure, a CND, and a larger disk with different orientations along the LOS and galactic axes. 7. 7.

We find a central hole that has a radius of 2.7 pc in the ${}^{12}{\rm CO(2-1)}$ integrated intensity map. However, this hole is filled by a dense gas tracer ${\rm CS(5-4)}$, which is centrally peaked, coexisting, and has a similar kinematic feature to the ${}^{12}{\rm CO(2-1)}$ line within the CND.

Acknowledgements

The authors would like to thank National Astronomical Observatory of Japan (NAOJ), National Institute of Natural Sciences (NINS) for supporting this work. D.D.N. delivers his gratitude to the supports from 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 NSF grant AST-1350389 and M.C. acknowledges the support from a Royal Society University Research Fellowship. T.I. sincerely thanks supports from the JSPS KAKENHI grant number 17K14247. M.I. is supported by JSPS KAKENHI grant number 15K05030.

The authors would like to thank Drs Naoteru Gouda, Makoto Miyoshi, Yano Taihei, and Daisuke Iono at NAOJ for enlightenment 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.

Facilities: ALMA and HST WFC3.

Software: Astropy, IDL, and emcee

Appendix A Supplementary Figures and Tables