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

Dieu D. Nguyen National Astronomical Observatory of Japan (NAOJ), National
Institute of Natural Sciences (NINS), 2-21-1 Osawa, Mitaka, Tokyo 181-8588,
Japan  Department of Physics and Astronomy, University of Utah, 115 South 1400
East, Salt Lake City, UT 84112, USA  Mark den Brok  Leibniz-Institut für
Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany  Anil
C. Seth  Department of Physics and Astronomy, University of Utah, 115 South
1400 East, Salt Lake City, UT 84112, USA  Satoru Iguchi  National Astronomical
Observatory of Japan (NAOJ), National Institute of Natural Sciences (NINS),
2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Department of Astronomical
Science, Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  Jenny E. Greene  Department of Astrophysics,
Princeton University, Princeton, NJ 08540, USA  Timothy Davis  School of
Physics and Astronomy, Cardiff University, Queens Buildings, The Parade,
Cardiff, CF24 3AA, UK  Masatoshi Imanishi National Astronomical Observatory of
Japan (NAOJ), National Institute of Natural Sciences (NINS), 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  Department of Astronomical Science, Graduate
University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Takuma Izumi NAOJ fellow  National Astronomical Observatory
of Japan (NAOJ), National Institute of Natural Sciences (NINS), 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  Michelle Cappellari Sub-department of
Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, UK  Nadine Neumayer Max Planck Institut
für Astronomie (MPIA), Königstuhl 17, D-69121 Heidelberg, Germany  Kristina
Nyland National Research Council, resident at the Naval Research Laboratory,
Washington, DC 20375, USA  Takafumi Tsukui  National Astronomical Observatory
of Japan (NAOJ), National Institute of Natural Sciences (NINS), 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  Department of Astronomical Science, Graduate
University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Kouichiro Nakanishi National Astronomical Observatory of
Japan (NAOJ), National Institute of Natural Sciences (NINS), 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  Department of Astronomical Science, Graduate
University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Phuong M. Nguyen  Department of Physics, Quy Nhon University,
170 An Duong Vuong, Quy Nhon, Vietnam  Quang L. Nguyen  AA$\&$AI lab, IBM
Canada, 120 Bloor Street East, Toronto, ON, M4Y 1B7, Canada  Graduate School
of Natural Sciences, Nagoya City University, 1 Yamanohata, Mizuho-cho, Nagoya,
Aichi 467-8501, Japan  Sabine Thater Leibniz-Institut für Astrophysik Potsdam
(AIP), An der Sternwarte 16, 14482 Potsdam, Germany  Martin Bureau  Sub-
department of Astrophysics, Department of Physics, University of Oxford, Denys
Wilkinson Building, Keble Road, Oxford OX1 3RH, UK  Kyoko Onishi  Research
Center for Space and Cosmic Evolution, Ehime University, 2-5 Bunkyo-cho,
Matsuyama, Ehime 790-8577, Japan  Karina T. Voggel Department of Physics and
Astronomy, University of Utah, 115 South 1400 East, Salt Lake City, UT 84112,
USA  Ngan M. Le  Department of Space and Aeronautics, University of Science
and Technology of Hanoi, 18 Hoang Quoc Viet, Hanoi, Vietnam  Trung V. Dinh
Institute of Physics, Vietnamese Academy of Science and Technology, 10, Dao
Tan, Ba Dinh, Hanoi, Vietnam

(Received February 11 2018; Revised XXX XX 2018; Accepted XXX XX 2019)

###### 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⋆)
Project. Our analysis is based on Atacama Large Millimeter/sub-millimeter
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 km s-1. Our dynamical models of the CND yield a
$M_{\rm BH}$ of $M_{\rm BH}=1.02^{+0.18}_{-0.15}\times 10^{7}$$M_{\odot}$and a
mass-to-light ratio in $H$-band of _M/L_ ${}_{\rm H}=0.66^{+1.44}_{-0.65}$
($M_{\odot}$/$L_{\odot}$). 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 $M_{\rm BH}$ in NGC 3504 probably more accurately than
the current commonly used of ${\rm{}^{12}CO(2-1)}$ due to its centralization.

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$ mag222https://www.oulu.fi/astronomy/S4G_PIPELINE4/MAIN/
decomp0087.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).

Figure 1: Left panel: _HST_ /WFC3 IR F160W image of NGC 3504 within the field-
of-view (FOV) of $140\arcsec\times 140\arcsec$ (10 kpc $\times$ 10 kpc). Right
panel: the zoom-in for the nucleus region within the FOV of $10\farcs 0\times
10\farcs 0$ (0.68 kpc $\times$ 0.68 kpc) overlaid with black ${}^{12}{\rm
CO(2-1)}$ integrated intensity contours from our ALMA observation from the
low-resolution measurement set (MS). Figure 2: The panel-a shows the zoom-in
map of the 1.3 mm continuum emission at the center of NGC 3504. The
synthesized beam of $0\farcs 044\times 0\farcs 031$ (2.9 pc $\times$ 2.1 pc)
is shown as an ellipse at the bottom left of the panel. The three rest panels
are the moments of the detected ${}^{12}{\rm CO(2-1)}$ emission created using
the masked moment technique described in Section 3.2.2 including the zeroth
moment map (panel-b), the first (panel-c) and the second moment (panel-d)
maps. The synthesized beam cannot be seen in these large field-of-view (FOV)
plots. White plus indicates the kinematic center and also the galaxy center.

##  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.

Figure 3: The integrated ${}^{12}{\rm CO(2-1)}$ spectrum extracted within the
nuclear region of $10\arcsec\times 10\arcsec$ (680 pc $\times$ 680 pc), where
includes all the detected emission. We observe the classic symmetric double
horn shape of a rotating disk.

Figure 4: Left panel: The position-velocity diagram (PVD) of the ${}^{12}{\rm
CO(2-1)}$ emission in the same field-of-view (FOV) as presented in moment
maps, extracted along the kinematic major axis with a slit of five pixels in
width ($0\farcs 065$ or 4.4 pc). We do not show the beamsize of $0\farcs
044\times 0\farcs 031$ (3.0 pc $\times$ 2.1 pc) and velocity channel width 10
km s-1 explicitly in this plot, as they are very small compared to the plotted
ranges. Clear channelization is present in the PVD, suggesting the velocity
dispersion in the gas is as small as compared to our channel width. Right
panel: The same PVD for the low-resolution cube with the beamsize of $0\farcs
221\times 0\farcs 164$ (15.0 pc $\times$ 11.2 pc).

###  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.

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\farcs 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\farcs 042\times 0\farcs 030$ | $0\farcs 221\times 0\farcs 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

††footnotemark:

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).

Figure 5: (a): The nuclear integrated intensity map of the ${}^{12}{\rm
CO(2-1)}$ gas emission zoom-in to the field-of-view (FOV) of $1\arcsec\times
1\arcsec$. The map shows a CND attenuation hole, which has a diameter of
$\sim$2.7 pc and is centered at the location of the putative SMBH (white
crosses). (b): The same FOV for our detected dense gas tracer ${\rm CS(5-4)}$
emission line in one of the continuum spectral windows. The ${\rm CS(5-4)}$
line is clearly centrally concentrated and peaks at the position of the
central SMBH where we observe the ${}^{12}{\rm CO(2-1)}$ attenuated hole. The
integrated spectrum (c) and radial profile (d) of ${\rm CS(5-4)}$, both show
the line is centrally filled rather than a hole. White pluses indicate the
kinematic center.

Figure 6: The moment maps of the detected ${\rm CS(5-4)}$ emission created
using the masked moment technique described in Section 3.2.2 including the
first (rotational velocity–left panel) and the second moment (velocity
dispersion–middle). Right panel: the position-velocity diagram of ${\rm
CS(5-4)}$, which has the same recession velocity and velocity width of
${}^{12}{\rm CO(2-1)}$, as well as the rotational signature around the
galactic center. We note that these maps have low-velocity resolution $\sim$40
km s-1.

##  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.

Figure 7: Top plots: the comparison between the _HST_ /WFC3 F160W photometry
of NGC 3504 (open squares) and its corresponding best-fit light MGE models
(red solid line), which are summed up from multiple Gaussians (color thin
lines), is shown in the left panel. We show the best-fit model that is
projected along 5∘-wide sector and has an angle of 30∘ between the major and
the minor axis. The fractional residuals (Data-Model)/Data is shown in the
right panel with the agreement between data and model within 15% across
100$\arcsec$. Bottom plots: the comparison between the F160W light
distribution and its best-fit MGE model in the form of 2D light surface
density contours for NGC 3504 at the whole galactic scale (left panel) and the
central zoom-in of $30\arcsec\times 30\arcsec$ field-of-view (FOV; right
panel). 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. The mismatch between data and model on the south side of the nucleus
is caused by the small jet due to AGN activity visually in Figure 1, which has
been masked out during the MGE fit.

##  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.

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

††footnotemark:

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$).

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

$j$ |  $\log$(Luminosity Density) | $\sigma^{\prime}$ | $a/b$
---|---|---|---
| ($L_{\odot}/{\rm 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

††footnotemark:

##  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$.

Figure 8: The nuclear gas morphological distribution of the combined cube
(high spatial resolution) shows with a cut along major axis through the center
of NGC 3504 integrated intensity map $(1\arcsec\times 1\arcsec)$. The data are
plotted in black diamond, while and our best-fit KinMS model for our chosen
surface brightness profile is plotted in red solid line. Figure 9: One
dimensional (1D) cumulative mass profile of stellar mass (black) and gas +
dust (red) components plotted within the radius of 6$\arcsec$ from the center.
Figure 10: The multidimensional parameter space posterior distributions that
are explored by the KinMS dynamical model to the combined cube (high
resolution) from the central $1\arcsec\times 1\arcsec$ field-of-view (FOV) of
NGC 3504. The top panel of each column is an one-dimensional (1D) histogram
shows the marginalized posterior distribution of that parameter, with 68%
(1$\sigma$) confidence interval, which is corresponding to the innermost
contour showed in the two-dimensional (2D) marginalization of those fitted
parameters in the panels below. We also show the distribution of the parameter
space within 2$\sigma$ (97.0%) and 3$\sigma$ (99.7%) confidence intervals
within the second and the outermost contours, respectively. See Table 4 for a
quantitative description of the likelihoods of all fitting parameters. Note
that $M_{\rm BH}$ is flat in log scale, others are in linear scales. Figure
11: Top panels: comparisons between the combined cube (high resolution) with a
few specific KinMS dynamical model including the best fit. The position-
velocity diagram (PVD) of the ${}^{12}{\rm CO(2-1)}$ emission in NGC 3504
extracted along the major-axis (orange scale and grey contours) in the same
manner we did in Figure 4. The model PVDs are extracted in an identical
fashion from models that are different by the central $M_{\rm BH}$ (red
contours) only. The left panel shows the case with no SMBH, the central panel
shows our best-fitting $M_{\rm BH}$, and the right panel has an overly massive
SMBH. The models with no/overly large SMBHs are clearly not the good fits to
the data in the central parts. Bottom panels: The central $1\arcsec\times
1\arcsec$ two-dimensional (2D) velocity map of the ${}^{12}{\rm CO(2-1)}$
emission in NGC 3504 extracted from the bottom-left panel of Figure 2 (left),
the velocity map derived from the KinMS model with the best-fit $M_{\rm BH}$
(middle), and the residual velocity map between the data and the model, (Data
- Model) (right). White plus indicates the kinematic center and also the
galaxy center. The synthesized beam of the observation is shown as an ellipse
at the bottom left of the panel.

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.

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_{\rm BH}$ in log scale ($\log_{10}M_{\rm BH}$; $M_{\odot}$) | 4.00 | $\longrightarrow$ | 9.00 | 7.01 |  $-$0.01, +0.01 |  $-$0.07, +0.07
Mass-to-light ratio in $H$-band (_M/L_ H; $M_{\odot}$/$L_{\odot}$) | 0.01 | $\longrightarrow$ | 3.00 | 0.66 |  $-$0.32, +0.32 |  $-$0.65, +1.44
Molecular Gas Disc: |  |  |  |  |  |
Position angle (PA; ∘) | 150.0 | $\longrightarrow$ | 360.0 | 334.81 |  $-$0.50, +0.45 |  $-$3.04, +2.03
Inclination angle ($i$; ∘) | 45.0 | $\longrightarrow$ | 90.0 | 62.13 |  $-$0.46, +0.49 |  $-$3.51, +3.11
Gas velocity dispersion ($\sigma$; km s-1) | 1.0 | $\longrightarrow$ | 80.0 | 38.30 |  $-$0.38, +0.37 |  $-$0.06, +0.06
Disk thickness ($d_{t}$; arcsec) | 0.01 | $\longrightarrow$ | 1.50 | 0.32 |  $-$0.02, +0.02 |  $-$0.05, +0.05
Gaussian peak (Gc; arcsec) | 0.01 | $\longrightarrow$ | 1.00 | 0.04 |  $-$0.00, +0.00 |  $-$0.02, +0.02
Gaussian HWHM (Gh; arcsec) | 0.01 | $\longrightarrow$ | 1.00 | 0.11 |  $-$0.01, +0.01 |  $-$0.03, +0.03
Exponential disk length scale ($R_{e}$; arcsec) | 0.01 | $\longrightarrow$ | 10.0 | 1.50 |  $-$0.00, +0.00 |  $-$0.02, +0.02
Nuisance Parameters: |  |  |  |  |  |
CO surface brightness scaling factor (Flux) | 10.0 | $\longrightarrow$ | 500.0 | 97.26 |  $-$0.52, +0.51 |  $-$2.21, +2.75
R.A. offset ($x_{c}$; arcsec) |  $-$0.10 | $\longrightarrow$ | +0.10 |  $-$0.01 |  $-$0.00, +0.00 |  $-$0.01, +0.01
Decl. offset ($y_{c}$; arcsec) |  $-$0.10 | $\longrightarrow$ | +0.10 |  $-$0.02 |  $-$0.00, +0.00 |  $-$0.01, +0.01
Systemic velocity offset ($v_{\rm off.}$; km s-1) |  $-$10.0 | $\longrightarrow$ | +10.0 |  $-$5.46 |  $-$0.69, +0.59 |  $-$2.61, +1.83

††footnotemark:

###  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.

Figure 12: Comparison of the velocity fields of our best-fit title-ring model
(left) with the ${}^{12}{\rm CO(2-1)}$ velocity field (center) and velocity
residuals (right) with the area of $0\farcs 8\times 0\farcs 8$ (55 pc $\times$
55 pc). The while contours denote the integrated intensity of ${}^{12}{\rm
CO(2-1)}$ emission in the the fitted area. Figure 13: Our NGC 3504 $M_{\rm
BH}$ (red encircle) in the context of the $M_{\rm BH}-M_{\rm Bulge}$ (left)
and $M_{\rm BH}-\sigma$ (right) scaling relations. Six measurements using ALMA
(Onishi et al., 2015; Barth et al., 2016a, b; Davis et al., 2017, 2018) and
Combined Array for Research in Millimeter-wave Astronomy (CARMA; Davis et al.,
2013; Onishi et al., 2017) observations are plotted in blue, while early-type
galaxies (ETGs, black dots within open circles) and late-type galaxies (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); Saglia
et al. (2016) for ETGs and LTGs are plotted in the dotted, dashed, long-dashed
lines, respectively. We also include the theoretical prediction of a
bimodality from Pacucci et al. (2018, red solid lines) and its 1$\sigma$ (red
dashed line) uncertainty.

##  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:

  1. 1.

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:

$v^{2}_{c}(r)=r\frac{d\Phi(r)}{dr}$ (1)

which for the spherically symmetric potential is equal to

$v^{2}_{c}(r)=G(M_{\star}(<r)+M_{\rm BH})r.$ (2)

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.

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_{\rm BH}$ in log scale ($\log M_{\rm BH}$; $M_{\odot}$) | 6.0 | $\longrightarrow$ | 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}$ | $\longrightarrow$ | 5.2 | 0.62 |  $-$0.01 +0.01 |  $-$0.04, +0.04
Position angle of the inner disk (PA1; ∘) | 65 | $\longrightarrow$ | 90 | 76.20 |  $-$0.17, +0.13 |  $-$0.53, +0.30
Position angle of the outer disk (PA1; ∘) | 65 | $\longrightarrow$ | 90 | 65.02 |  $-$0.01, +0.02 |  $-$0.07, +0.07
Inclination angle of the inner disk ($i_{1}$; ∘) | 25 | $\longrightarrow$ | 90 | 54.75 |  $-$0.17, +0.18 |  $-$0.60, +0.52
Inclination angle of the outer disk ($i_{2}$; ∘) | 25 | $\longrightarrow$ | 90 | 61.84 |  $-$0.20, +0.20 |  $-$0.60, +0.62
Transition radius between two disks (Rbreak; arcsec) | $10^{-2}$ | $\longrightarrow$ | 1.0 | 0.29 |  $-$0.0, +0.0 |  $-$0.04, +0.04
Model 2: |  |  |  |  |  |
$M_{\rm BH}$ in log scale ($\log M_{\rm BH}$; $M_{\odot}$) | 6.0 | $\longrightarrow$ | 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}$ | $\longrightarrow$ | 5.2 | 0.63 |  $-$0.03, +0.02 |  $-$0.08, +0.07
Inclination angle of the inner ellipse ($i_{1}$; ∘) | 30 | $\longrightarrow$ | 90 | 47.54 |  $-$1.20, +1.34 |  $-$3.72, +3.92
Inclination angle of the outer ellipse ($i_{2}$; ∘) | 30 | $\longrightarrow$ | 90 | 60.42 |  $-$1.06, +0.99 |  $-$1.21, +1.82
Model 3: |  |  |  |  |  |
$M_{\rm BH}$ in log scale ($\log M_{\rm BH}$; $M_{\odot}$) | 6.0 | $\longrightarrow$ | 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}$ | $\longrightarrow$ | 5.2 | 0.61 |  $-$0.00 +0.03 |  $-$0.03, +0.09
Position angle of the inner disk (PA; ∘) | 65 | $\longrightarrow$ | 90 | 65.51 |  $-$0.83, +0.94 |  $-$2.35, +2.90
Inclination angle of the inner disk ($i$; ∘) | 50 | $\longrightarrow$ | 90 | 57.77 |  $-$0.96, +1.99 |  $-$3.23, +5.65
Central hole radius (rinner; arcsec) | $10^{-2}$ | $\longrightarrow$ | 0.1 | 0.04 |  $-$0.01, +0.01 |  $-$0.03, +0.03

††footnotemark:

###  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:

  1. 1.

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

Figure 14: The inner fitted region channel maps of our combined (high
resolution) ALMA ${}^{12}{\rm CO(2-1)}$ data cube in the velocity range where
emission is detected. The colored regions with thin black contours show the
areas detected with more than 2.5$\sigma$ significance. The thick red contours
are plotted with the same contour levels from our best-fitting KinMS model.
Figure 15: Comparisons between the low-resolution cube with a few specific
KinMS dynamical model including the best fit. All notations and illustrations
in this figure are similar to Figure 11. 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_{\rm BH}$ in log scale ($\log_{10}M_{\rm BH}$; $M_{\odot}$) | 4.00 | $\longrightarrow$ | 9.00 | 7.18 |  $-$0.04, +0.04 |  $-$0.10, +0.13
Mass-to-light ratio in $H$-band (_M/L_ H; $M_{\odot}$/$L_{\odot}$) | 0.01 | $\longrightarrow$ | 3.00 | 0.67 |  $-$0.02, +0.03 |  $-$0.62, +1.01
Molecular Gas Disc: |  |  |  |  |  |
Position angle (PA; ∘) | 150.0 | $\longrightarrow$ | 360.0 | 330.17 |  $-$1.42, +0.90 |  $-$4.24, +3.25
Inclination angle ($i$; ∘) | 45.0 | $\longrightarrow$ | 90.0 | 60.31 |  $-$1.43, +1.34 |  $-$4.51, +4.01
Gas velocity dispersion ($\sigma$; km s-1) | 1.0 | $\longrightarrow$ | 80.0 | 35.50 |  $-$0.53, +0.45 |  $-$1.56, +1.42
Disk thickness ($d_{t}$; arcsec) | 0.01 | $\longrightarrow$ | 1.50 | 0.41 |  $-$0.03, +0.03 |  $-$0.09, +0.12
Gaussian HWHM (Gh; arcsec) | 0.01 | $\longrightarrow$ | 1.00 | 0.35 |  $-$0.03, +0.03 |  $-$0.10, +0.13
Nuisance Parameters: |  |  |  |  |  |
CO surface brightness scaling factor (Flux) | 10.0 | $\longrightarrow$ | 500.0 | 149.01 |  $-$0.02, +0.02 |  $-$0.07, +0.09
R.A. offset ($x_{c}$; arcsec) |  $-$0.10 | $\longrightarrow$ | +0.10 |  $-$0.02 |  $-$0.00, +0.00 |  $-$0.01, +0.01
Decl. offset ($y_{c}$; arcsec) |  $-$0.10 | $\longrightarrow$ | +0.10 |  $-$0.03 |  $-$0.00, +0.00 |  $-$0.01, +0.01
Systemic velocity offset ($v_{\rm off.}$; km s-1) |  $-$10.0 | $\longrightarrow$ | +10.0 |  $-$0.02 |  $-$0.22, +0.02 |  $-$0.06, +0.07

††footnotemark:

## References

  * (1)
  * Abolfathi et al. (2018) Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42
  * Afanasiev et al. (2018) Afanasiev, A. V., Chilingarian, I. V., Mieske, S., et al. 2018, MNRAS, 477, 4856
  * Ahn et al. (2017) Ahn, C. P., Seth, A. C., den Brok, M., et al. 2017, ApJ, 839, 72
  * Ahn et al. (2018) Ahn, C. P., Seth, A. C., Cappellari, M., et al. 2018, ApJ, 858, 102
  * Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473
  * Avila et al. (2012) Avila, R. J., Hack, W. J., & STScI AstroDrizzle Team. 2012, in American Astronomical Society Meeting Abstracts, Vol. 220, American Astronomical Society Meeting Abstracts 220, #135.13
  * Baldassare et al. (2015) Baldassare, V. F., Reines, A. E., Gallo, E., & Greene, J. E. 2015, ApJ, 809, L14
  * Barai et al. (2014) Barai, P., Viel, M., Murante, G., Gaspari, M., & Borgani, S. 2014, MNRAS, 437, 1456
  * Barth et al. (2016a) Barth, A. J., Boizelle, B. D., Darling, J., et al. 2016a, ApJ, 822, L28
  * Barth et al. (2016b) Barth, A. J., Darling, J., Baker, A. J., et al. 2016b, ApJ, 823, 51
  * Barth et al. (2004) Barth, A. J., Ho, L. C., Rutledge, R. E., & Sargent, W. L. W. 2004, ApJ, 607, 90
  * Barvainis (1987) Barvainis, R. 1987, ApJ, 320, 537
  * Begeman (1987) Begeman, K. G. 1987, PhD thesis, , Kapteyn Institute, (1987)
  * Beifiori et al. (2012) Beifiori, A., Courteau, S., Corsini, E. M., & Zhu, Y. 2012, MNRAS, 419, 2497
  * Bell et al. (2003) Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289
  * Bigiel et al. (2008) Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846
  * Bolatto et al. (2013) Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Nature, 499, 450
  * Bonoli et al. (2014) Bonoli, S., Mayer, L., & Callegari, S. 2014, MNRAS, 437, 1576
  * Boselli et al. (2015) Boselli, A., Fossati, M., Gavazzi, G., et al. 2015, A&A, 579, A102
  * Buta & Combes (1996) Buta, R., & Combes, F. 1996, Fund. Cosmic Phys., 17, 95
  * Buta & Crocker (1992) Buta, R., & Crocker, D. A. 1992, AJ, 103, 1804
  * Caplar et al. (2015) Caplar, N., Lilly, S. J., & Trakhtenbrot, B. 2015, ApJ, 811, 148
  * Cappellari (2002) Cappellari, M. 2002, MNRAS, 333, 400
  * Cappellari (2008) Cappellari. 2008, MNRAS, 390, 71
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Chilingarian et al. (2018) Chilingarian, I. V., Katkov, I. Y., Zolotukhin, I. Y., et al. 2018, ApJ, 863, 1
  * Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693
  * Davis (2014) Davis, T. A. 2014, MNRAS, 443, 911
  * Davis et al. (2013) Davis, T. A., Bureau, M., Cappellari, M., Sarzi, M., & Blitz, L. 2013, Nature, 494, 328
  * Davis et al. (2017) Davis, T. A., Bureau, M., Onishi, K., et al. 2017, MNRAS, 468, 4675
  * Davis et al. (2018) Davis et al.. 2018, MNRAS, 473, 3818
  * de Vaucouleurs (1975) de Vaucouleurs, G. 1975, ApJ, 202, 319
  * Deller & Middelberg (2014) Deller, A. T., & Middelberg, E. 2014, AJ, 147, 14
  * den Brok et al. (2015) den Brok, M., Seth, A. C., Barth, A. J., et al. 2015, ApJ, 809, 101
  * Desroches et al. (2009) Desroches, L.-B., Greene, J. E., & Ho, L. C. 2009, ApJ, 698, 1515
  * Di Matteo et al. (2008) Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33
  * Dong et al. (2012) Dong, X.-B., Ho, L. C., Yuan, W., et al. 2012, ApJ, 755, 167
  * Emsellem et al. (1994) Emsellem, E., Monnet, G., & Bacon, R. 1994, A&A, 285, 723
  * Epinat et al. (2008) Epinat, B., Amram, P., Marcelin, M., et al. 2008, MNRAS, 388, 500
  * Erroz-Ferrer et al. (2015) Erroz-Ferrer, S., Knapen, J. H., Leaman, R., et al. 2015, MNRAS, 451, 1004
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455
  * Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9
  * Fiacconi & Rossi (2016) Fiacconi, D., & Rossi, E. M. 2016, MNRAS, 455, 2
  * Fiacconi & Rossi (2017) Fiacconi & Rossi. 2017, MNRAS, 464, 2259
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gallo et al. (2008) Gallo, E., Treu, T., Jacob, J., et al. 2008, ApJ, 680, 154
  * Gallo et al. (2010) Gallo, E., Treu, T., Marshall, P. J., et al. 2010, ApJ, 714, 25
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, Vol. 5, No. 1, p. 65-80, 2010, 5, 65
  * Graham et al. (2001) Graham, A. W., Erwin, P., Caon, N., & Trujillo, I. 2001, ApJ, 563, L11
  * Graham & Scott (2015) Graham, A. W., & Scott, N. 2015, ApJ, 798, 54
  * Greene & Ho (2007) Greene, J. E., & Ho, L. C. 2007, ApJ, 670, 92
  * Greene et al. (2010) Greene, J. E., Peng, C. Y., Kim, M., et al. 2010, ApJ, 721, 26
  * Greene et al. (2016) Greene, J. E., Seth, A., Kim, M., et al. 2016, ApJ, 826, L32
  * Gültekin et al. (2009) Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198
  * Häring & Rix (2004) Häring, N., & Rix, H.-W. 2004, ApJ, 604, L89
  * Ho & Filippenko (1993) Ho, L. C., & Filippenko, A. V. 1993, Ap&SS, 205, 19
  * Ho et al. (1997) Ho, L. C., Filippenko, A. V., & Sargent, W. L. W. 1997, ApJS, 112, 315
  * Ho et al. (2009) Ho, L. C., Greene, J. E., Filippenko, A. V., & Sargent, W. L. W. 2009, ApJS, 183, 1
  * Ho et al. (1993) Ho, L. C., Shields, J. C., & Filippenko, A. V. 1993, ApJ, 410, 567
  * Imanishi et al. (2018) Imanishi, M., Nakanishi, K., Izumi, T., & Wada, K. 2018, ApJ, 853, L25
  * Inayoshi & Haiman (2016) Inayoshi, K., & Haiman, Z. 2016, ApJ, 828, 110
  * Izumi et al. (2018) Izumi, T., Wada, K., Fukushige, R., Hamamura, S., & Kohno, K. 2018, ApJ, 867, 48
  * Jedrzejewski (1987) Jedrzejewski, R. I. 1987, MNRAS, 226, 747
  * Jedrzejewski et al. (1987) Jedrzejewski, R. I., Davies, R. L., & Illingworth, G. D. 1987, AJ, 94, 1508
  * Keel (1984) Keel, W. C. 1984, ApJ, 282, 75
  * Kenney et al. (1993) Kenney, J. D. P., Carlstrom, J. E., & Young, J. S. 1993, ApJ, 418, 687
  * Knapen et al. (2002) Knapen, J. H., Pérez-Ramírez, D., & Laine, S. 2002, MNRAS, 337, 808
  * Kormendy & Bender (2012) Kormendy, J., & Bender, R. 2012, ApJS, 198, 2
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
  * Kormendy & Richstone (1995) Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581
  * Krajnović et al. (2006) Krajnović, D., Cappellari, M., de Zeeuw, P. T., & Copin, Y. 2006, MNRAS, 366, 787
  * Krajnović et al. (2018) Krajnović, D., Cappellari, M., McDermid, R. M., et al. 2018, MNRAS, 477, 3030
  * Kuno et al. (2000) Kuno, N., Nishiyama, K., Nakai, N., et al. 2000, PASJ, 52, 775
  * Kuno et al. (2007) Kuno, N., Sato, N., Nakanishi, H., et al. 2007, PASJ, 59, 117
  * Kuo et al. (2011) Kuo, C. Y., Braatz, J. A., Condon, J. J., et al. 2011, ApJ, 727, 20
  * Läsker et al. (2016) Läsker, R., Greene, J. E., Seth, A., et al. 2016, ApJ, 825, 3
  * Laurikainen et al. (2004) Laurikainen, E., Salo, H., Buta, R., & Vasylyev, S. 2004, MNRAS, 355, 1251
  * Lo (2005) Lo, K. Y. 2005, ARA&A, 43, 625
  * Lodato & Natarajan (2006) Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285
  * Maksym et al. (2013) Maksym, W. P., Ulmer, M. P., Eracleous, M. C., Guennou, L., & Ho, L. C. 2013, MNRAS, 435, 1904
  * Maoz et al. (1998) Maoz, D., Koratkar, A., Shields, J. C., et al. 1998, AJ, 116, 55
  * Marconi & Hunt (2003) Marconi, A., & Hunt, L. K. 2003, ApJ, 589, L21
  * McConnell et al. (2013) McConnell, N. J., Chen, S.-F. S., Ma, C.-P., et al. 2013, ApJ, 768, L21
  * Miller et al. (2015) Miller, B. P., Gallo, E., Greene, J. E., et al. 2015, ApJ, 799, 98
  * Mitzkus et al. (2017) Mitzkus, M., Cappellari, M., & Walcher, C. J. 2017, MNRAS, 464, 4789
  * Miyoshi et al. (1995) Miyoshi, M., Moran, J., Herrnstein, J., et al. 1995, Nature, 373, 127
  * Moran et al. (2014) Moran, E. C., Shahinyan, K., Sugarman, H. R., Vélez, D. O., & Eracleous, M. 2014, AJ, 148, 136
  * Netzer (2015) Netzer, H. 2015, ARA&A, 53, 365
  * Neumayer et al. (2007) Neumayer, N., Cappellari, M., Reunanen, J., et al. 2007, ApJ, 671, 1329
  * Nguyen et al. (2017) Nguyen, D. D., Seth, A. C., den Brok, M., et al. 2017, ApJ, 836, 237
  * Nguyen et al. (2018) Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2018, ApJ, 858, 118
  * Nguyen et al. (2019) Nguyen et al.. 2019, arXiv e-prints, arXiv:1901.05496
  * Nicholson et al. (1992) Nicholson, R. A., Bland-Hawthorn, J., & Taylor, K. 1992, ApJ, 387, 503
  * Onishi et al. (2017) Onishi, K., Iguchi, S., Davis, T. A., et al. 2017, ArXiv e-prints, arXiv:1703.05247
  * Onishi et al. (2015) Onishi, K., Iguchi, S., Sheth, K., & Kohno, K. 2015, ApJ, 806, 39
  * Pacucci & Loeb (2018) Pacucci, F., & Loeb, A. 2018, ArXiv e-prints, arXiv:1810.12302
  * Pacucci et al. (2018) Pacucci, F., Loeb, A., Mezcua, M., & Martín-Navarro, I. 2018, ApJ, 864, L6
  * Pacucci et al. (2017) Pacucci, F., Natarajan, P., Volonteri, M., Cappelluti, N., & Urry, C. M. 2017, ApJ, 850, L42
  * Pacucci et al. (2015) Pacucci, F., Volonteri, M., & Ferrara, A. 2015, MNRAS, 452, 1922
  * Park et al. (2016) Park, K., Ricotti, M., Natarajan, P., Bogdanović, T., & Wise, J. H. 2016, ApJ, 818, 184
  * Paturel et al. (2000) Paturel, G., Fang, Y., Petit, C., Garnier, R., & Rousseau, J. 2000, A&AS, 146, 19
  * Peterson (1982) Peterson, C. J. 1982, PASP, 94, 409
  * Quillen et al. (1992) Quillen, A. C., de Zeeuw, P. T., Phinney, E. S., & Phillips, T. G. 1992, ApJ, 391, 121
  * Randriamampandry et al. (2015) Randriamampandry, T. H., Combes, F., Carignan, C., & Deg, N. 2015, MNRAS, 454, 3743
  * Reines et al. (2013) Reines, A. E., Greene, J. E., & Geha, M. 2013, ApJ, 775, 116
  * Reines & Volonteri (2015) Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82
  * Rieke & Low (1972) Rieke, G. H., & Low, F. J. 1972, ApJ, 176, L95
  * Roediger & Courteau (2015) Roediger, J. C., & Courteau, S. 2015, MNRAS, 452, 3209
  * Russell (2002) Russell, D. G. 2002, ApJ, 565, 681
  * Saglia et al. (2016) Saglia, R. P., Opitsch, M., Erwin, P., et al. 2016, ApJ, 818, 47
  * Salo et al. (2015) Salo, H., Laurikainen, E., Laine, J., et al. 2015, ApJS, 219, 4
  * Satyapal et al. (2009) Satyapal, S., Böker, T., Mcalpine, W., et al. 2009, ApJ, 704, 439
  * Schawinski et al. (2007) Schawinski, K., Thomas, D., Sarzi, M., et al. 2007, MNRAS, 382, 1415
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Scott et al. (2013) Scott, N., Graham, A. W., & Schombert, J. 2013, ApJ, 768, 76
  * Seth et al. (2010) Seth, A. C., Cappellari, M., Neumayer, N., et al. 2010, ApJ, 714, 713
  * Seth et al. (2014) Seth, A. C., van den Bosch, R., Mieske, S., et al. 2014, Nature, 513, 398
  * She et al. (2017) She, R., Ho, L. C., & Feng, H. 2017, ApJ, 842, 131
  * Sheth et al. (2010) Sheth, K., Regan, M., Hinz, J. L., et al. 2010, PASP, 122, 1397
  * Silk & Rees (1998) Silk, J., & Rees, M. J. 1998, A&A, 331, L1
  * Stone et al. (2017) Stone, N. C., Küpper, A. H. W., & Ostriker, J. P. 2017, MNRAS, 467, 4180
  * Terrazas et al. (2017) Terrazas, B. A., Bell, E. F., Woo, J., & Henriques, B. M. B. 2017, ApJ, 844, 170
  * Thater et al. (2017) Thater, S., Krajnović, D., Bourne, M. A., et al. 2017, A&A, 597, A18
  * Thornton et al. (2008) Thornton, C. E., Barth, A. J., Ho, L. C., Rutledge, R. E., & Greene, J. E. 2008, ApJ, 686, 892
  * Tristram & Schartmann (2011) Tristram, K. R. W., & Schartmann, M. 2011, A&A, 531, A99
  * Valluri et al. (2005) Valluri, M., Ferrarese, L., Merritt, D., & Joseph, C. L. 2005, ApJ, 628, 137
  * van den Bosch & de Zeeuw (2010) van den Bosch, R. C. E., & de Zeeuw, P. T. 2010, MNRAS, 401, 1770
  * van den Bosch et al. (2016) van den Bosch, R. C. E., Greene, J. E., Braatz, J. A., Constantin, A., & Kuo, C.-Y. 2016, ApJ, 819, 11
  * van Wassenhove et al. (2010) van Wassenhove, S., Volonteri, M., Walker, M. G., & Gair, J. R. 2010, MNRAS, 408, 1139
  * Verolme et al. (2002) Verolme, E. K., Cappellari, M., Copin, Y., et al. 2002, MNRAS, 335, 517
  * Voggel et al. (2018) Voggel, K. T., Seth, A. C., Neumayer, N., et al. 2018, ApJ, 858, 20
  * Volonteri (2010) Volonteri, M. 2010, A&A Rev., 18, 279
  * Volonteri (2012) Volonteri, M. 2012, in American Institute of Physics Conference Series, Vol. 1480, American Institute of Physics Conference Series, ed. M. Umemura & K. Omukai, 289–296
  * Volonteri & Bellovary (2012) Volonteri, M., & Bellovary, J. 2012, Reports on Progress in Physics, 75, 124901
  * Volonteri et al. (2015) Volonteri, M., Capelo, P. R., Netzer, H., et al. 2015, MNRAS, 452, L6
  * Volonteri et al. (2008) Volonteri, M., Lodato, G., & Natarajan, P. 2008, MNRAS, 383, 1079