A Precision Measurement of the Mass of the Black Hole in NGC 3258 from High-Resolution ALMA Observations of its Circumnuclear Disk

TL;DR

This study uses high-resolution ALMA observations of NGC 3258's circumnuclear disk to precisely measure its supermassive black hole's mass, demonstrating ALMA's capability for accurate black hole mass determinations in early-type galaxies.

Contribution

First precise black hole mass measurement in NGC 3258 using ALMA CO kinematics, avoiding dust extinction uncertainties and modeling disk warps.

Findings

Black hole mass = 2.249×10^9 M☉ with 0.18% statistical uncertainty

ALMA data enables mass profile inference without optical imaging

Demonstrates ALMA's potential for precise SMBH measurements

Abstract

We present $\sim0.10^{\prime\prime}-$resolution Atacama Large Millimeter/submillimeter Array (ALMA) CO(2$-$1) imaging of the arcsecond-scale ($r \approx 150$ pc) dusty molecular disk in the giant elliptical galaxy NGC 3258. The data provide unprecedented resolution of cold gas disk kinematics within the dynamical sphere of influence of a supermassive black hole, revealing a quasi-Keplerian central increase in projected rotation speed rising from 280 km s$^{-1}$ at the disk's outer edge to $>400$ km s$^{-1}$ near the disk center. We construct dynamical models for the rotating disk and fit beam-smeared model CO line profiles directly to the ALMA data cube. Our models incorporate both flat disks and tilted-ring disks that provide a better fit of the mildly warped structure in NGC 3258. We show that the exceptional angular resolution of the ALMA data makes it possible to infer the host galaxy's mass profile within $r=150$ pc solely from the ALMA CO kinematics, without relying on optical or near-infrared imaging data to determine the stellar mass profile. Our model therefore circumvents any uncertainty in the black hole mass that would result from the substantial dust extinction in the galaxy's central region. The best model fit yields $M_\mathrm{BH} = 2.249\times10^9$ $M_\odot$ with a statistical model-fitting uncertainty of just 0.18\%, and systematic uncertainties of 0.62\% from various aspects of the model construction and 12\% from uncertainty in the distance to NGC 3258. This observation demonstrates the full potential of ALMA for carrying out highly precise measurements of $M_\mathrm{BH}$ in early-type galaxies containing circumnuclear gas disks

Tables

Table 1. Table 1: H 𝐻 H -band MGE Parameters

$j$	${\log }_{10}$ $I_{H,j}$ ($L_{\odot }$ pc-2)	${\sigma }_j^{\prime }$ (arcsec)	$q{\prime }_j$
(1)	(2)	(3)	(4)
1	3.75	0.37	0.78
2	4.14	0.82	0.83
3	3.85	1.37	0.75
4	3.72	1.77	0.77
5	3.75	2.68	0.79
6	3.54	4.44	0.85
7	2.93	8.41	0.84
8	2.67	11.7	1.00
9	2.32	14.0	0.85
10	2.12	23.9	0.94
11	1.83	41.6	0.81
12	1.20	64.8	0.96
13	1.24	115	0.72
14	0.35	400	0.72

Table 2. Table 2: Dynamical Model Properties

Model	Cycle	Mass Model	Disk Structure	${\sigma }_{\mathrm{turb}}\left(r\right)$
A1	2	MGE; $A_H=0$	Flat	Uniform
B1–B4	2	MGE; $A_H=0,0.31,0.75,1.50$	Flat	Gaussian
C1	4	MGE; $A_H=0$	Flat	Uniform
D1–D4	4	MGE; $A_H=0,0.31,0.75,1.50$	Flat	Gaussian
E1	4	MGE; $A_H=0.75$	Tilted ring	Gaussian
F1	4	$v_{\mathrm{ext}}$	Tilted ring	Gaussian

Table 3. Table 3: Modeling Results

Model	$M_{\mathrm{BH}}$	${{\rm Y}}_H$	$i$	$\Gamma$	${\sigma }_1$	${\sigma }_0$	$r_0$	$\mu$	$x_c$	$y_c$	$v_{\mathrm{sys}}$	$f_0$	${\chi }_{\nu }^2$
	(${10}^9$ $M_{\odot }$)	($M_{\odot }$/$L_{\odot }$)	(°)	(°)	(km s-1)	(km s-1)	(pc)	(pc)	(″)	(″)	(km s-1)
A1	2.386	3.16	46.0	76.8	7.44	$\cdots$	$\cdots$	$\cdots$	$-$0.022	$-$0.020	2760.76	1.02	1.783
B1	2.386	3.16	46.0	76.8	7.42	4.75	0.14	$-$0.46	$-$0.022	$-$0.020	2760.76	1.02	1.784
B2	2.309	2.88	45.9	76.8	7.74	16.9	0.62	$-$0.53	$-$0.022	$-$0.020	2760.76	1.02	1.806
B3	2.216	2.65	46.0	76.8	7.71	4.11	1.04	$-$4.53	$-$0.021	$-$0.021	2760.79	1.02	1.858
B4	2.087	2.42	46.1	76.8	7.86	0.14	1.54	$-$9.78	$-$0.021	$-$0.021	2760.79	1.01	1.935
C1	2.280	2.72	49.0	77.0	10.5	$\cdots$	$\cdots$	$\cdots$	$-$0.001	$-$0.003	2760.84	1.06	1.229
D1	2.276	2.73	49.0	77.0	6.83	8.58	$-$20.3	65.9	$-$0.001	$-$0.003	2760.83	1.07	1.219
D2	2.215	2.46	49.0	77.0	7.14	8.44	$-$25.0	59.5	$-$0.001	$-$0.003	2760.87	1.07	1.217
D3	2.144	2.26	49.0	77.0	7.74	8.02	$-$36.44	46.9	$-$0.001	$-$0.003	2760.91	1.07	1.217
D4	2.059	2.04	49.0	77.0	8.72	7.07	$-$52.4	29.9	$-$0.001	$-$0.003	2760.97	1.07	1.219
E1	2.203	2.18	24.2$-$49.8	76.2$-$96.4	6.54	22.5	$-$53.2	83.6	$-$0.002	$-$0.003	2760.83	1.07	1.180
F1	2.249	$\cdots$	27.5$-$49.3	76.2$-$93.6	6.32	21.9	$-$51.3	84.7	$-$0.002	$-$0.003	2760.82	1.07	1.179
	(0.004)	$\cdots$	$\cdots$	$\cdots$	(0.16)	(0.40)	(0.47)	(0.39)	(0.001)	(0.001)	(0.07)	(0.002)	$\cdots$

Table 4. Table 4: Tilted-ring Model Parameters and Host Galaxy Circular Velocity Profile

Disk Radius		$\Gamma$	$q$	$v_{\mathrm{ext}}$
(arcsec)	(pc)	(°)		(km s-1)
0.066	10	92.31 (0.22)	0.887 (0.033)	49.7 (2.5)
0.099	15	93.61 (0.24)	0.793 (0.016)	72.0 (2.8)
0.131	20	88.14 (0.21)	0.758 (0.007)	85.2 (4.1)
0.197	30	80.60 (0.22)	0.722 (0.004)	117.6 (3.6)
0.264	40	78.21 (0.17)	0.672 (0.002)	131.1 (1.6)
0.395	60	77.17 (0.09)	0.662 (0.001)	154.9 (1.3)
0.527	80	76.88 (0.07)	0.668 (0.001)	188.1 (0.8)
0.659	100	76.17 (0.05)	0.652 (0.001)	210.7 (0.4)
0.823	125	77.68 (0.07)	0.659 (0.001)	241.1 (0.4)
0.988	150	78.95 (0.09)	0.656 (0.002)	273.4 (0.9)

Equations

$${\left(\frac{F^{\prime }}{F^0}\right)}_{\lambda }\approx f+w\left[\frac{1-10^{-A_{\lambda }/2.5}}{0.921A_{\lambda }}\right]+b\left[10^{-A_{\lambda }/2.5}\right].$$

$$\Delta (B-I)=-2.5{\log }_{10}\left[{\left(\frac{F^{\prime }}{F^0}\right)}_B{\left(\frac{F^{\prime }}{F^0}\right)}_I^{-1}\right],$$

Full text

galaxies: elliptical and lenticular, galaxies: nuclei, galaxies: kinematics and dynamics, galaxies: individual: NGC 3258

1 Introduction

Supermassive black holes (BHs), spanning a mass range of $\sim 10^{6}-10^{10}$ $M_{\odot}$, are key constituents of the centers of likely all massive galaxies (Magorrian et al., 1998; Kormendy & Ho, 2013). Although BHs gravitationally dominate only the innermost regions of galaxies, their masses ($M_{\mathrm{BH}}$) strongly correlate with several large-scale host galaxy properties, such as the stellar velocity dispersion ($\sigma_{\star}$; Gebhardt et al., 2000; Ferrarese & Merritt, 2000) and bulge luminosity ($L$; Kormendy & Richstone, 1995). These local relationships encapsulate a fossil record of BH and galaxy growth through accretion and merger events, and suggest a co-evolution of central BHs and host galaxies. The local $M_{\mathrm{BH}}-\sigma_{\star}$ and $M_{\mathrm{BH}}-L$ relationships (McConnell & Ma, 2013; Kormendy & Ho, 2013; Saglia et al., 2016; van den Bosch et al., 2016) are also widely employed in estimating $M_{\mathrm{BH}}$ for both nearby and distant galaxies across the Hubble sequence.

The BH census remains incomplete, particularly for the most luminous early-type galaxies (ETGs), including brightest cluster galaxies (BCGs) and brightest group galaxies (BGGs). Furthermore, a growing sample of BH masses reveals that the correlations are more complicated than initially thought and may not consistently apply to all galaxy types. For example, predicted $M_{\mathrm{BH}}$ values for the most luminous ETGs using their measured stellar velocity dispersions are in tension with masses estimated from the $M_{\mathrm{BH}}-L$ relationship, with the discrepancy reaching an order of magnitude at $M_{\mathrm{BH}}\sim 10^{10}$ $M_{\odot}$ (Lauer et al., 2007b; Bernardi et al., 2007). The few BCGs with measured BH masses (Dalla Bontà et al., 2009; McConnell et al., 2012; Rusli et al., 2013b) suggest a steeper $M_{\mathrm{BH}}-\sigma_{\star}$ relationship and may point to different evolutionary processes within cluster centers (e.g., Krajnović et al., 2018). However, large uncertainties in the masses of several of the most massive BHs prevent any secure interpretation.

Presently, $\sim$100 dynamical $M_{\mathrm{BH}}$ measurements have been made, primarily by modeling stellar or ionized gas kinematics (e.g., Kormendy & Ho, 2013; Saglia et al., 2016). Reliably measuring $M_{\mathrm{BH}}$ requires modeling the kinematics of tracers that originate within the BH sphere of influence $r_{\mathrm{g}}\approx GM_{\mathrm{BH}}/\sigma^{2}_{\star}$, where the BH dominates the host galaxy’s gravitational potential. The confidence of a BH mass measurement hinges on how well the kinematic observations resolve $r_{\mathrm{g}}$. Obtaining more than a couple of resolution elements across $r_{\mathrm{g}}$ remains challenging for the current generation of optical/near-infrared (NIR) telescopes, even when using adaptive optics (AO). Rusli et al. (2013b) model the stellar kinematics of several luminous ETGs and find that the $M_{\mathrm{BH}}$ uncertainties, and the potential biases introduced by model systematics, increase when the angular resolution of the observations exceeds $r_{\mathrm{g}}$. For stellar-dynamical modeling, these systematics include assumptions about the intrinsic galaxy shape, inclusion of a dark matter halo, and adoption of a spatially constant stellar mass-to-light ($M/L$) ratio (Gebhardt & Thomas, 2009; van den Bosch & de Zeeuw, 2010; McConnell et al., 2013). Non-circular motion and the treatment of gas turbulence can bias gas-dynamical BH masses (e.g., van der Marel & van den Bosch, 1998; Walsh et al., 2010). In the few cases where both stellar and gas-dynamical modeling techniques have been applied to the same galaxy, the inferred BH masses frequently disagree, and discrepancies of a factor $2-4$ are common (Gebhardt et al., 2011; Rusli et al., 2011; Kormendy & Ho, 2013; Walsh et al., 2013; Barth et al., 2016a).

Data that highly resolve $r_{\mathrm{g}}$ have the potential to avoid nearly all such serious systematics. Kinematic measurements of 22 GHz H2O emission from radii $\ll r_{\mathrm{g}}$ in megamaser disks enable $M_{\mathrm{BH}}$ determinations with percent-level precision (e.g., Miyoshi et al., 1995; Kuo et al., 2011). Unfortunately, such disks are rare (e.g., Braatz et al., 1996) and tend to be found in late-type galaxies with black hole masses clustering in a narrow range about $\sim$107 $M_{\odot}$. Surveys to identify megamaser disks within ETGs have thus far been unsuccessful (e.g., van den Bosch et al., 2016), so a different method is needed to make precision BH mass measurements in the most massive galaxies.

Molecular gas tracers are a promising new avenue for reliably measuring BH masses, especially in ETGs. Recent 12CO surveys (Combes et al., 2007; Young et al., 2011; Alatalo et al., 2013; Bolatto et al., 2017; Zabel et al., 2018) find central, regularly-rotating cold molecular gas in roughly 10% of all nearby elliptical and S0 galaxies. Low turbulent velocity dispersions indicate that the molecular gas in these disks is a better tracer of the underlying gravitational potential than ionized gas in ETGs. Until recently, however, mm/sub-mm arrays were only able to resolve the nuclear gas kinematics at $r<r_{\mathrm{g}}$ for a very small number of galaxies. For one such nearby ETG at $D\sim 16$ Mpc, Davis et al. (2013) mapped rapid CO gas rotation at 0$\farcs$25 resolution with the Combined Array for Research in Millimeter-wave Astronomy (CARMA), and demonstrated that BH masses can be constrained using mm-wavelength molecular gas as tracers.

The Atacama Large Millimeter/submillimeter Array (ALMA) now offers the possibility of routinely carrying out molecular-line observations that resolve $r_{\mathrm{g}}$, given its increased sensitivity and significantly higher angular resolution relative to previous facilities. ALMA observations are highly sensitive probes of molecular gas within the central $\sim$kpc region of luminous ETGs (Boizelle et al., 2017, hereinafter Paper I) and have opened a new avenue for $M_{\mathrm{BH}}$ determination (Barth et al., 2016a, b; Onishi et al., 2017; Davis et al., 2017, 2018; Smith et al., 2019) via detection and modeling of the central high-velocity rotation around the BH. However, even when ALMA observations resolve $r_{\mathrm{g}}$, a central nearly-Keplerian rise in rotation speed is not typically seen, indicating a dearth of gas at locations close to the BH. In many cases, the data reveal only a modest central rise in peak rotation speed originating from gas in the outer portion of the BH sphere of influence, suggesting the presence of a central hole in the CO distribution at $r<r_{\mathrm{g}}$. Other cases are found to exhibit a resolved central hole in the CO distribution with radius larger than $r_{\mathrm{g}}$ (e.g., Paper I, ). For high-precision measurement of $M_{\mathrm{BH}}$, the ideal configuration is a disk with bright CO emission extending down to radii much smaller than $r_{\mathrm{g}}$, from which the central rotation speed due to the BH’s gravity would rise far above the rotation speed at larger radii due to the host galaxy’s mass. ALMA observations published to date indicate that disks with these properties are fairly rare among the local ETG population, with only a very small fraction exhibiting signatures of very rapid central rotation from radii deep within $r_{\mathrm{g}}$.

NGC 3258 was first observed by ALMA as part of the Cycle 2 program described in Paper I. That $\sim$0$\farcs$44 resolution imaging revealed bright CO(2$-$1) emission from a rapidly rotating nuclear gas disk, with a spatially unresolved central rise in line-of-sight velocity ($v_{\mathrm{LOS}}$) extending to $\sim$500 km s*-1* relative to the systemic velocity ($v_{\mathrm{sys}}$) and rising to $\gtrsim 200$ km s*-1* above the rotation speed of the outer disk. These attributes made NGC 3258 a promising target for high-resolution ALMA imaging in order to determine its BH mass to high precision. This E1 galaxy has a bulge stellar velocity dispersion of $\sigma_{\star}=260\pm 10$ km s*-1* and $K$-band absolute magnitude of $M_{K}=-24.33\pm 0.45$ mag (from the HyperLeda database; Makarov et al., 2014). We adopt a distance modulus $m-M=32.53\pm 0.27$ mag based on surface brightness fluctuation measurements (SBF; Tonry et al., 2001), which corresponds to a luminosity distance $D_{L}=31.9\pm 3.9$ Mpc. Using an observed redshift $z=0.0092091$ from our initial dynamical modeling results (that is very close to other optical measurements; de Vaucouleurs et al., 1991), this $D_{L}$ corresponds to an angular size distance of 31.3 Mpc, for which $1\arcsec$ spans a physical scale of 151.8 pc. NGC 3258 is one of two BGGs that dominate the dynamically young Antlia cluster (Hess et al., 2015), a somewhat poor cluster with $\sim$400 member galaxies (Ferguson & Sandage, 1990). Optical long-slit spectroscopy reveals only weak evidence for stellar rotation but a large central stellar velocity dispersion of $\sim$400 km s*-1* (Koprolin & Zeilinger, 2000; De Bruyne et al., 2004). As no atomic gas reservoir is detected within this galaxy (Hess et al., 2015), the cold gas in NGC 3258 appears to be primarily molecular. Mid-infrared Spitzer spectra show significant nuclear polycyclic aromatic hydrocarbon emission that, together with other nuclear diagnostic line diagrams, suggests a recent ($\sim$200 Myr) burst of star formation (Rampazzo et al., 2013).

In this paper, we present $\sim$0$\farcs$10 resolution ALMA Cycle 4 CO(2$-$1) observations of NGC 3258. A factor of four improvement in angular resolution compared with the earlier Cycle 2 data fully resolves gas rotation within $r_{\mathrm{g}}$ and enables measurement of the BH mass to an unprecedented level of precision for a giant elliptical galaxy. The extraordinary resolution of the gas kinematics within $r_{\mathrm{g}}$ in NGC 3258 makes it possible to constrain the spatially extended host galaxy mass distribution within the galaxy’s inner arcsecond solely from the ALMA kinematic data, in contrast to the traditional approach of using high-resolution optical/NIR imaging data to measure and deproject the host galaxy luminosity profile and assuming a spatially uniform stellar mass-to-light ratio. Measuring the host galaxy’s mass profile from the kinematic data makes it possible to avoid an uncertainty of order several percent in $M_{\mathrm{BH}}$ that would result from the uncertain extinction of the host galaxy’s central stellar luminosity profile. This method is particularly advantageous for systems such as NGC 3258 in which the central region of the galaxy is highly obscured by dust, as will be the case for nearly all CO-bright galaxies targeted for ALMA observations to measure $M_{\mathrm{BH}}$. Our measurement yields statistical model-fitting uncertainties that are significantly smaller than the systematic uncertainties resulting from issues such as localized irregularities in the gas disk kinematics. We carry out a variety of tests to estimate these model-fitting systematics and find that they are below the $\sim 1\%$ level except for the uncertainty in the galaxy’s distance, which contributes $>$10% systematic uncertainty to the error budget, as is generally the case for nearly all dynamical BH mass measurements.

This paper is organized as follows. In Section 2, we present Hubble Space Telescope ($HST$) optical and NIR broadband imaging of NGC 3258 and measurements of the galaxy’s light profile. We describe models for extinction and reddening due to the inclined circumnuclear dust disk embedded within the galaxy, and demonstrate that the disk is very optically thick at visible wavelengths. We use the $HST$ data to derive dust-corrected models for the host galaxy’s intrinsic luminosity distribution that we then deproject and employ as a component of the traditional approach for BH mass measurement. We introduce the new ALMA Cycle 4 observations in Section 3. In Section 4, we describe our gas-dynamical modeling method and discuss results when fitting these models to the Cycle 2 and 4 ALMA data cubes. We present model-fitting results for the simple case of a geometrically flat disk, and for a tilted-ring model that more closely matches the disk’s mildly warped structure. We compare results from models employing a dust-corrected stellar mass profile measured from $HST$ imaging and models based on a new method that determines the host galaxy’s radial mass profile solely from the ALMA CO kinematics. In Section 5, we discuss the implications of high-precision ALMA BH mass measurements and place NGC 3258 in the context of M_{\mathrm{BH}}$$-host galaxy relationships.

2 Optical and Near-Infrared Observations

A typically key input into gas-dynamical models is the stellar contribution to a galaxy’s gravitational potential. We used $HST$ Wide Field Camera 3 (WFC3) NIR images to determine the luminous mass distribution in the galaxy’s central region and Advanced Camera for Surveys (ACS) Wide Field Channel (WFC) observations to characterize the dust disk properties. In order to probe the galaxy’s outskirts, we supplemented the $HST$ WFC3 data with ground-based, wide-field images from the Carnegie-Irvine Galaxy Survey (CGS; Ho et al., 2011). Below, we summarize the observations, data reduction procedures, surface brightness measurements, and disk extinction models.

2.1 HST Imaging

We observed NGC 3258 in one orbit on 5 June 2017 as part of program GO-14920 with $HST$ WFC3 through the IR channel using the F110W and F160W ($J$ and $H$) filters. We took four MULTIACCUM exposures in each filter with the SPARS25 (NSAMP$=$12$-$13) and STEP50 (NSAMP$=$13) modes, employing a large box dither pattern that always kept the galaxy nucleus in one corner of the detector. We processed the data through the CALWF3 pipeline and used AstroDrizzle (Gonzaga et al., 2012) to produce cleaned, distortion-corrected images with a pixel scale of 0$\farcs$08 pixel*-1* in each filter. The final $J$ and $H$-band images cover a 3$\farcm$4$\times$3$\farcm$1 field of view, and have total integration times of 18 minutes and 23 minutes, respectively. The images are dominated by galaxy light over the full WFC3/IR field of view, and so we did not perform background subtraction at this stage. The final $HST$ images have an angular resolution of 0$\farcs$21 ($J$) and 0$\farcs$22 ($H$), determined by averaging the full widths at half maximum (FWHMs) of several foreground stars. In Figure 1, we present the $H$-band mosaic of NGC 3258 and its inner 4″$\times$4″ region, which illustrates the substantial extinction by the central dust disk.

In addition, we retrieved ACS/WFC F435W and F814W ($B$ and $I$) images of NGC 3258 from the $HST$ archive. These $B$ and $I$ images were taken over three orbits as part of program GO-9427, and have integration times of 89 min and 38 min, respectively. We processed the raw ACS data using the CALACS pipeline, included corrections for charge transfer inefficiency, and then drizzled the geometrically rectified ACS exposures in each filter. The final images have an angular resolution of $\sim$0$\farcs$12 and cover the galaxy’s central 3$\farcm$4$\times$3$\farcm$3 region. In Figure 1, we show the $B-I$ map, in which a nearly azimuthally symmetric, $\sim$2$\farcs$4$-$wide dust disk is visible (see also De Bruyne et al., 2004; Capetti & Balmaverde, 2005).

2.2 CGS Imaging

We used ground-based optical data from CGS to complement the $HST$ images. The CGS observations and data reduction are described by Ho et al. (2011), and the processed images have a 9′$\times$9′ field of view and a pixel scale of 0$\farcs$26. We selected the CGS $V$-band observation instead of a redder filter that would better trace old stellar populations in order to avoid the “red halo” effect. This instrumental effect, stronger in longer-wavelength filters, adds an extended feature to the point-spread function (PSF) wings, potentially affecting measurements of the galaxy’s brightness profile (Huang et al., 2013). The sensitivity of the $V$-band CGS image reaches 26.9 mag arcsec*-2*, and the image was taken in $\sim$1″ seeing. The CGS $V$-band image of NGC 3258 is displayed in Figure 1.

2.3 Stellar Luminosity Profile

After masking out foreground stars and galaxies, we measured NGC 3258’s surface brightness from the $H$ and $V$-band images in regions spaced logarithmically in radius and equally in angle (Cappellari, 2002). The average position angle (PA) of 77° was determined using the central $R<10$″ region of the NIR mosaic. (Throughout this paper, we use $R$ to denote projected radius on the plane of the sky, and $r$ to denote radial distance within the galaxy.) Using the surface brightness measurements at radii between $70-100$″ along the galaxy’s major axis, we determined the $H$-band background level and the $V-H$ color needed to align the two profiles. We found a best-fit background level of $H=$ 20.8 mag arcsec*-2* and a color of $V-H=$ 2.40 mag (consistent with optical-NIR colors of elliptical galaxies at large radii; e.g., Schombert et al., 1993), which we then applied to the $H$-band surface brightness measurements at radii 0$\farcs$07$-$100″ and to the $V$-band surface brightness measurements at radii 70$-$300″, respectively. This produced $H$-band surface brightness profiles measured at nineteen angles between $0-90\arcdeg$ from the major axis and extending out to $\sim$300″ ($\sim$45.5 kpc), or 5$-$7 times the estimated half-light radius ($R_{\mathrm{e}}$; Lauberts & Valentijn, 1989; Dirsch et al., 2003). Although we used the surface brightness profile along the major axis to establish the $H$-band background and $V-H$ color, there was good visual agreement between the $H$ and scaled $V$-band surface brightness measurements at all angles. We also note that $V-H$ gradients are negligible at radii $\gtrsim 60$″, and therefore we assumed that $V-H$ remains constant with radius when adjusting the large-scale $V$-band measurements.

We corrected for foreground Galactic reddening of $A_{H}=0.041$ mag (Schlafly & Finkbeiner, 2011) and modeled the galaxy’s $H$-band surface brightness with a two-dimensional (2D) multi-Gaussian expansion (MGE; Emsellem et al., 1994; Cappellari, 2002) after masking out the most dust-obscured regions to the south of the nucleus between $R\sim 0\farcs 15-0\farcs 8$. Although individual Gaussian components do not have physical meaning, MGEs are commonly used to represent a wide variety of surface brightness profiles and allow for the luminosity density to be determined through an analytical deprojection. We required that each Gaussian component have the same center and PA, and constrained the observed flattening (the ratio between the projected major and minor axis of the 2D Gaussian, $q^{\prime}$) to be $>0.72$. Such a restriction avoids highly flattened components that would limit the range of inclination angles ($i$) for which an MGE model can be deprojected. Prior to parameterizing the $H$-band surface brightness with an MGE model, we generated a Tiny Tim (Krist & Hook, 2004) PSF, which was dithered and drizzled identically to the $H$-band mosaic of the galaxy. We applied the MGE formalism to the $HST$ PSF, modeling it as the sum of six concentric circular Gaussians. This six-component PSF was taken into account while fitting the galaxy’s surface brightness using the MGE (Cappellari, 2002).

The final MGE model of the galaxy consists of 14 concentric elliptical Gaussians, for which the best-fit parameters are provided in Table 1. This MGE is a good description of the $H$-band surface brightness measurements, and we show a comparison between the MGE model and the data along the galaxy major axis in Figure 2. The galaxy’s total $H$-band luminosity is $L_{H}=1.1\times 10^{11}$ $L_{\odot}$, measured within the central 300″ (45.5 kpc), and we find that $R_{\mathrm{e}}\approx 66\arcsec$ (10.6 kpc) in this filter.

Assuming that the galaxy has an oblate axisymmetric shape and is inclined at the same angle as the molecular gas disk ($i\approx 48\arcdeg$; see Section 4), we deprojected the $H$-band MGE model and numerically integrated the resulting stellar luminosity densities assuming an initial $H$-band $M/L$ ratio $\Upsilon_{H}=1$ to determine the stellar contribution to circular velocity ($v_{\mathrm{c}}^{\star}$) as a function of radius. During gas-dynamical modeling, these $v_{\mathrm{c}}^{\star}$ values are then scaled by $\sqrt{\Upsilon_{H}}$ that is a free parameter of the fits. We do not include contributions from a dark matter halo, as these are negligible within the central few kpc of the galaxy (De Bruyne et al., 2004).

2.4 Disk Extinction Modeling

Modeling the extinction from the circumnuclear dust disk is essential in order to derive accurate models of the host galaxy stellar mass profile from the $HST$ images. The $H$-band nuclear morphology seen in Figure 1 is suggestive of a ring-like obscuring structure at $R\approx 0\farcs 5$ in which the extinction is most pronounced on the southern (near) side of the disk. At somewhat larger radii, ring-like obscuration also appears in optical $HST$ images (see Capetti & Balmaverde 2005 and Paper I, ). The $H$-band surface brightness profile (Figure 2) shows an apparent break at $R\approx 2\arcsec$ to a nearly flat inner slope at smaller radii, suggesting substantial extinction of the $H$-band light by the disk within the galaxy’s inner arcsecond.

The $HST$ $B-I$ color map reveals that the region of largest optical color excess (relative to the host galaxy color of $B-I\approx 2.30$ mag outside the dust disk region) is confined to a ring located at $R\approx 1\arcsec$. On the southern side of the ring, the $B-I$ color is $\sim 0.3-0.4$ mag redder than the host galaxy outside the dust disk. At radii $<1\arcsec$, the color map shows a patchy structure with multiple concentric ringlets and a greater overall color excess of $(B-I)\approx 0.1-0.15$ mag relative to the host galaxy on larger scales.

If interpreted as due to a foreground screen of extinction in front of NGC 3258, these color excesses would indicate modest extinction values reaching $A_{V}\sim 0.3$ mag (or $A_{H}\sim 0.06$ mag) in a ring and that decreases by a factor of two towards the central portion of the disk (assuming a standard Galactic extinction law; Cardelli et al., 1989). However, such small extinctions would not be sufficient to create the observed absorption feature in the nuclear $H$-band light (Figure 1). Furthermore, the CO(2$-$1) surface brightness of the disk suggests an average $V$-band extinction as high as 5$-$10 mag over the disk surface (Paper I). This seeming discrepancy is the result of the disk’s location in the midplane of the host galaxy: interpreting the observed color excess with foreground-screen models would greatly underpredict the disk’s true optical depth (Tran et al., 2001).

In this situation, the obscuring structure is an inclined, dusty disk in the midplane of the galaxy. Starlight originating from in front of the disk is unobscured, while light from within and behind the disk is attenuated. In the limit of very high optical depth in a very thin disk, light from the far side of the disk would be completely obscured, and a $B-I$ color map would not reveal any color excess. The maximum observed $B-I$ color excess would occur for some moderate value of disk optical depth that would permit some reddened starlight to pass through the disk. For an inclined, embedded dust disk, the near side of the disk would be expected to show a larger color excess than the far side as the near side of the disk obscures a greater fraction of the host galaxy’s starlight (Elmegreen & Block, 1999).

To examine the relationship between disk optical depth and observed color excess, we employed a simple embedded-screen model following the method described by Viaene et al. (2017). In this model, the obscuring structure is treated as a thin, inclined disk bisecting the galaxy. Along a given line of sight, the fraction $b$ of total stellar light originating behind the disk is obscured by simple screen extinction, while the fraction $f$ in front remains unaffected. For full generality in the case of a thick disk, a small fraction ($w=1-f-b$) of the total light may originate within the disk and therefore experience “mixed” attenuation. Rewriting Equation 6 of Viaene et al. (2017) in terms of the extinction $A_{\lambda}$, the wavelength-dependent ratio $F^{\prime}/F^{0}$ of observed to intrinsic integrated stellar light takes the form

[TABLE]

We used the same $R_{V}=3.1$ extinction law to characterize $A_{\lambda}$, and for simplicity assumed a very thin ($w\rightarrow 0$) disk. To determine fractions $f$ and $b$ across the arcsecond-scale disk, we populated a model galaxy cube with stellar densities deprojected from the $H$-band MGE solution and adopted $i=48\arcdeg$ for the dust disk based on initial gas-dynamical modeling results in Section 4. We evaluated Equation 1 at the pivot wavelengths of the ACS and WFC3 filters to generate predictions for the opacity-dependent color excess at each spatial location:

[TABLE]

with a similar form for $\Delta(J-H)$.

In Figure 3, we show the modeled color excesses $\Delta(B-I)$ and $\Delta(J-H)$ as a function of the intrinsic extinction $A_{V}$ of the obscuring disk, extracted at three locations each in order to illustrate the effect of the disk inclination on the color excess at different locations in the disk. These major and minor axis positions coincide with the elliptical rings of maximal color excess observed at $R\sim 1\farcs 1$ and $\sim 0\farcs 5$ for the $B-I$ and $J-H$ color maps, respectively. As expected, the color excess predicted by the model is very small for both very low and very high disk optical depth, and reaches a maximum value at intermediate extinction. The predicted $B-I$ color excess peaks at a disk extinction of $A_{V}\approx 1.5-2$ mag, while the $J-H$ color excess peaks at $A_{V}\approx 5-6$ mag. Away from these extinctions corresponding to peak color excesses, the observed color excess no longer corresponds to a unique $A_{V}$ value. The color excesses on the near side of the disk are predicted to be more than twice the value of the color excesses on the far side.

Remarkably, this simple model predicts maximum color excess values that are in very close agreement with both the $B-I$ and $J-H$ color maps of NGC 3258, as seen in Figure 3. This consistency indicates that the disk optical depth rises from $A_{V}\sim 1.5$ mag near the disk edge to at least 5 mag at $R\sim 0\farcs 5$. Such high intrinsic extinction corresponds to a substantial attenuation of $H$-band light within the galaxy’s inner arcsecond, as illustrated in the bottom panel of Figure 3.

These model results imply that the central region of the dust disk is sufficiently opaque to absorb a significant fraction of the $H$-band galaxy light originating from behind the disk, and we conclude that extinction is responsible for some of the central flattening in the $H$-band radial profile. However, there is no straightforward method to correct the observed $H$-band radial profile for extinction based on the color excess maps. Recovering the intrinsic stellar surface brightness via spectral energy distribution measurements at each spatial location would require realistic radiative transfer modeling (e.g., Camps & Baes, 2015) that accounts for the disk geometry and thickness, dust scattering, and extinction within the disk. Possible additional contributions of light from recent star formation in the disk or a weak active nucleus would further complicate any extinction correction method based on the observed color excess maps. Such modeling is beyond the scope of this work.

Instead, we adopted a simpler approach to examine the impact of extinction on the inferred $v_{\mathrm{c}}^{\star}$ profile by adjusting the central $H$-band surface brightness profile to correct for three fiducial values of disk extinction that bracket the likely range. The inner $H$-band stellar surface brightness follows a double power-law profile, so we fit the central $R\lesssim 10\arcsec$ of the mosaic with a PSF-convolved 2D Nuker function (using GALFIT; Peng et al., 2002), which yields an inner cusp slope $\gamma=0.01$ and a break radius $r_{\mathrm{b}}=1.5\arcsec$ (corresponding to $\sim$230 pc) that extends slightly beyond the dust disk radius. This $r_{\mathrm{b}}$ is consistent with those measured for other massive ETGs (e.g., Faber et al., 1997; Lauer et al., 2005). After fixing all other Nuker parameters, we adjusted $\gamma$ to 0.09, 0.17, and 0.26 to approximately correct for absorption of ¼, ½, and ¾ of the integrated stellar light originating behind the disk (for $R\lesssim 0\farcs 25$), respectively, corresponding to disk intrinsic optical depths of $A_{H}=0.31$, 0.75, and 1.50 mag (or $A_{V}=1.67$, 4.04, and 8.09 mag). The maximum $\gamma$ we use is within the range generally associated with core galaxies ($\gamma\leq 0.3$; e.g., Faber et al., 1997). For each $A_{H}$ case, we created a new model image by seamlessly replacing the dust-obscured region (out to $R=1\farcs 5$) with the associated GALFIT product. These dust-corrected $H$-band surface brightness profiles are shown in Figure 4. We parameterized each new model image using the MGE method and used the results to derive three additional, “dust-corrected” circular velocity profiles. In Section 4, we employ all four $v_{\mathrm{c}}^{\star}$ profiles (the original and the three dust-corrected profiles) in gas-dynamical models to quantify the impact of dust obscuration on the final $M_{\mathrm{BH}}$ measurement.

3 ALMA Data

3.1 Observations and data processing

The new Cycle 4 data were obtained in ALMA Program 2016.1.00854.S during 7-8 August 2017 in the C40-7 configuration, which had baselines ranging from 21 to 3696 m. Observations consisted of a single pointing with three $\sim$2 GHz$-$bandwidth spectral windows, one of which was centered on the redshifted 12CO(2$-$1) 230.538 GHz line while the remaining two measured the continuum at average (sky) frequencies of 228.4 and 243.0 GHz. Three execution blocks were carried out in good weather conditions (precipitable water vapor of 0.3$-$1.0 mm) for a total on-source integration time of 135 minutes. Line and continuum spectral windows were sampled using channel widths of 3.91 MHz (after $8\times$ online channel averaging) and 15.6 MHz, respectively. The data were flux calibrated using ALMA quasar standards J1037$-$2934 and J1107$-$4449, which have absolute flux calibration accuracies of $\sim$10% (Fomalont et al., 2014). We have propagated this uncertainty into all subsequent flux and flux density measurements.

Prior to line and continuum imaging, we flagged and calibrated the Cycle 4 visibilities using version 4.7.2 of the Common Astronomy Software Applications (CASA) pipeline. CASA TCLEAN deconvolution with Briggs ($r=0.5$; Briggs, 1995) weighting results in a synthesized beam with FWHM $\theta_{\mathrm{FWHM}}=0\farcs 11\times 0\farcs 08$ at PA=88°. We first imaged the line-free channels (with a 5.2 GHz total bandwidth) to produce a continuum map with a point-source sensitivity of $\sim$11 $\mu$Jy beam*-1*. After $uv$-plane continuum subtraction, we imaged the primary spectral window into a line cube with 7.81 MHz channels (corresponding to rest-frame velocity widths of $\sim$10.2 km s*-1*) that have typical rms sensitivities of $\sim$0.27 mJy beam*-1*.

These ALMA Cycle 4 data are a significant improvement in both angular resolution and sensitivity over the Cycle 2 observations of this target from Program 2013.1.00229.S, which are described in Paper I. However, the sparse central $uv$-plane coverage of the C40-7 configuration results in a $\sim$1$\farcs$2 maximum recoverable scale that may resolve out some smoothly-distributed emission in the 2$\farcs$4$-$wide disk. We therefore simultaneously imaged together the Cycle 2 and 4 visibilities using a multiscale deconvolution. After natural weighting of the visibilities, we obtained a continuum map with $\theta_{\mathrm{FWHM}}=0\farcs 14\times 0\farcs 11$ at PA=$-$82° and an rms level of 9.8 $\mu$Jy beam*-1*. Briggs ($r=0.5$) weighting produced a line cube with $\theta_{\mathrm{FWHM}}=0\farcs 12\times 0\farcs 09$ at PA=89° with $\sim$0.23 mJy beam*-1* sensitivities in $\sim$10.2 km s*-1* channels at the 0$\farcs$015 pixel*-1* scale. Although incorporating these shorter-baseline data does slightly expand the synthesized beam major and minor axes to a geometric mean of $\sim$0$\farcs$10, we recovered more CO line and extended continuum emission than from imaging of the Cycle 4 data set alone. For the remainder of this paper, we refer to the results of our Cycle 2+4 multiscale deconvolution of continuum and spectral line data simply as Cycle 4 imaging.

3.2 Emission Line Properties

In the Cycle 4 line cube, we detect CO(2$-$1) emission out to $R\sim 1\farcs 05$ and in channels spanning 900 km s*-1*. The highest velocity line emission (relative to the disk systemic velocity $v_{\mathrm{sys}}\approx 2761$ km s*-1*) is directly adjacent to the galaxy nucleus. We integrated the cube flux densities in each channel over the elliptical disk area to determine its velocity profile (Figure 5). The double-horned profile shape is very similar to that seen in the Cycle 2 data, while the total line flux of $S_{\rm CO(2-1)}\Delta v=27.40\pm 0.15\pm 2.74$ Jy km s*-1* (statistical and systematic uncertainties, respectively) is slightly higher than the Cycle 2 value of $23.89\pm 2.39$ Jy km s*-1* reported in Paper I.

We extracted a position-velocity diagram (PVD) from the Cycle 4 cube along a PA = 77° with a spatial extraction width equal to the average synthesized beam FWHM of 0$\farcs$10 (see Figure 5). The CO line-of-sight velocities span the same range as in the Cycle 2 PVD, but in the Cycle 4 data the CO emission is resolved into a tight locus of quasi-Keplerian rotation arising from the point mass BH and extended galaxy mass distribution. These CO emission-line velocities rise to a peak on either side of the nucleus, tracing gas rotation to within $\sim$20 pc of the galaxy center. This remarkable PVD structure resolves the central rise in rotation velocity far better than any published ALMA observation of circumnuclear gas in any other galaxy. In contrast to the Cycle 2 data, spatial blurring of high-velocity and low-velocity emission due to beam smearing in the inner disk is almost completely eliminated. The line-of-sight velocity $|v_{\mathrm{LOS}}-v_{\mathrm{sys}}|$ of this innermost CO emission reaches $\sim$480 km s*-1*. Assuming a regularly rotating disk inclined by $i\approx 48\arcdeg$, the corresponding circular velocity of $v_{\rm c}\approx 650$ km s*-1* at this radius would suggest $M_{\mathrm{BH}}\approx 2\times 10^{9}$ $M_{\odot}$. This value of $M_{\mathrm{BH}}$ implies $r_{\mathrm{g}}\approx 0\farcs 9$, which in turn indicates that nearly all of the dust disk lies within $r_{\mathrm{g}}$.

As described in Section 4, we fit gas-dynamical models directly to both the Cycle 2 and 4 CO line cubes. For visualization purposes, we parameterized the line-of-sight velocity distributions using Gauss-Hermite functions (GH; van der Marel & Franx, 1993). For low S/N regions at the disk center and near the edge, adjacent spectra were combined together prior to line profile fitting (using a Voronoi tessellation of a preliminary CO flux map; Cappellari & Copin, 2003). We display GH moment maps for the Cycle 2 data in Paper I and for Cycle 4 in Figure 6, which includes the integrated CO(2$-$1) line flux, $v_{\mathrm{LOS}}$, and velocity dispersion $\sigma_{\mathrm{LOS}}$. Due to beam smearing, both the Cycle 2 and 4 moment maps show high $|h_{3}|$ and $|h_{4}|$ values of up to $\sim 0.25-0.30$ for radii $\lesssim 0\farcs 2$. For the lower-resolution data set, these non-Gaussian coefficients remain elevated in coherent patterns out to $R\sim 1\arcsec$.

Making the same assumptions about the CO-to-H2 conversion factor $\alpha_{\rm CO}$ as in Paper I, the CO flux measured from the Cycle 4 data implies a total H2+He mass $M_{\rm gas}=(1.0\pm 0.3)\times 10^{8}$ $M_{\odot}$ for the gas disk (including uncertainties in galaxy distance and flux calibration). The Cycle 2 and 4 data both show a centrally concentrated CO flux distribution (Figure 7), and the Cycle 4 imaging with a beam size corresponding to 17 pc partially resolves the CO(2$-$1) emission into large, cloud-like knots (Figure 8). Clumpy emission-line structure appears to be common for molecular gas disks in ETGs when observed at similar physical resolutions (Utomo et al., 2015; Barth et al., 2016a; Davis et al., 2017, 2018). We identify a central hole in CO surface brightness with a radius of $\sim 0\farcs 13$ that corresponds to the innermost emission detected in the PVD.

The gas kinematics are nearly Keplerian close to the disk center, flattening out to $v_{\mathrm{LOS}}\sim 280$ km s*-1* for $R>0\farcs 6$ due to the increasing contributions of host galaxy mass at larger radii. Examination of the Cycle 4 PVD shows an asymmetry in the peak velocities on either side of the nucleus, which reach $+483$ and $-414$ km s*-1* relative to $v_{\mathrm{sys}}$ on the receding and approaching sides of the disk, respectively. For $R<0\farcs 2$, the approaching (western) emission appears to show sub-Keplerian rotation velocities.

The observed velocity field also exhibits minor kinematic warping, most noticeably at radii $\lesssim 0\farcs 25$. To characterize deviations from coplanar rotation, we decomposed the $v_{\rm LOS}$ map using kinemetry (Krajnović et al., 2006) to measure the kinematic PA $\Gamma_{k}$ and axis ratio $q_{k}$, as well as circular ($k_{1}$) and non-circular ($k_{5}$) velocity components, as a function of radius. Results are shown in Figure 9. While the primary kinematic twist $\Delta\Gamma_{k}\sim 10\arcdeg$ occupies the inner half arcsecond, the disk remains slightly warped out to the edge of the detected CO(2$-$1) emission. Beam smearing reduces the velocity amplitude along the line of nodes for $R\lesssim 0\farcs 2$, while at greater radii $k_{1}\approx v_{\rm c}\sin i$. The measured $q_{k}$ values show a central rise to unity that may in part be the result of finite angular resolution (i.e., circularization of the nuclear velocity field; see Paper I). For thin disk rotation $q_{k}\approx\cos i$, so $q_{k}\rightarrow 0.67$ with increasing radius suggests an outer disk inclination angle of $\sim$48°. Similar to the Cycle 2 kinemetry results, the coefficient ratio $k_{5}/k_{1}\lesssim 0.02$ at all radii, suggesting only negligible deviations from circular rotation despite the evident warping of the disk.

The observed CO(2$-$1) line dispersion ranges from $\sim$7$-$415 km s*-1*. The highest values found around the nucleus and on either side of the major axis can be attributed primarily to beam smearing and intrapixel velocity gradients (Barth et al., 2016b). However, the $\sigma_{\mathrm{LOS}}$ field reaches its maximum not at the disk center as expected but $\sim$0$\farcs$05 northward. The central $\sigma_{\mathrm{LOS}}\sim 300$ km s*-1* may in part be lower due to the coincident hole in CO(2$-$1) flux, with some broad, low S/N line profile wings buried beneath the noise. Along the disk major axis, the measured line dispersion rapidly decreases to $\lesssim 50$ km s*-1* for $R\gtrsim 0\farcs 07$ and falls below 7 km s*-1* near the disk edge.

4 Dynamical Modeling

In this section, we present results from dynamical modeling of the circumnuclear disk in NGC 3258. We begin with models for the simple case of a geometrically flat disk, and then consider a tilted-ring model designed to provide a better fit to the disk’s warped geometry. We employ two different methods to constrain the mass distribution of the host galaxy: first, the standard approach of using the dust-corrected MGE models that is more widely applicable to data that do not highly resolve gas rotation within $r_{\mathrm{g}}$, and second, a method using only the ALMA CO kinematics to determine the extended mass profile. We also consider models with different prescriptions for the spatial variation of the turbulent velocity dispersion of the molecular gas. Models are fit to the ALMA Cycle 4 data, and we also describe model fits to the lower-resolution Cycle 2 data to illustrate the effect of angular resolution on the $M_{\mathrm{BH}}$ determination. For clarity, the various models used in this paper are labeled and described in Table 2. Models A and B are fit to the Cycle 2 data, while models C$-$F are fit to the Cycle 4 data. We adopt model F1 as our final best result: this includes the tilted-ring disk structure, a spatially varying turbulent velocity dispersion, and the extended mass profile determined from the ALMA kinematic data. We also conduct additional tests based on variants of model F1 to estimate the systematic uncertainties in the BH mass. Unless otherwise specified, all modeling results described below refer to the Cycle 4 data; fits to the Cycle 2 data are described in §4.1.3.

4.1 Initial flat-disk models

4.1.1 Method

We first describe the basic flat-disk modeling procedure, which builds on methods developed for the analysis of $HST$ ionized-gas kinematics and uses forward modeling of line profiles from a rotating disk (e.g., Macchetto et al., 1997; van der Marel & van den Bosch, 1998; Barth et al., 2001). A major difference is that we fit models directly to the observed ALMA data cube, making use of all available information, rather than fitting models to velocity and velocity dispersion curves extracted from the data, as was done for $HST$ gas-dynamical measurements. Our flat-disk modeling method was previously used to measure the black hole mass in NGC 1332 (Barth et al., 2016a, b), and is similar to procedures used by other groups to measure black hole masses from molecular gas kinematics (e.g, Davis et al., 2017; Onishi et al., 2017). Barth et al. (2016b) present a detailed description of the method, which we summarize here.

The model calculation starts by determining the circular velocity as a function of radius for a thin, flat disk orbiting in the combined gravitational potential of a central black hole and the spatially extended mass distribution of the host galaxy. Line-of-sight projections of the disk rotation velocity are determined at each point on the sky for a given disk inclination and major axis position angle and for an assumed distance to the galaxy. Then, a spectral cube is generated by assuming an intrinsically Gaussian line profile at each point in the disk, with some specified turbulent velocity dispersion. The model cube is constructed to match the observed frequency spacing of the ALMA data cube for direct comparison. At each spatial grid point, the total flux in the modeled line profile is determined using a map of the CO surface brightness distribution determined from the ALMA observation. Each velocity channel of the model is convolved with the ALMA synthesized beam (an elliptical Gaussian). In order to capture details of sub-pixel gradients in rotation velocity near the disk center, the model calculation and beam convolution are carried out on an oversampled spatial grid (relative to the ALMA data cube) and the modeled cubes are then downsampled to match the ALMA pixel scale. Models are optimized by $\chi^{2}$ minimization using a downhill simplex minimization method (Press et al., 1992) by fitting the calculated cubes directly to the ALMA data cube. Further details of these steps are described below.

These basic models employ at least nine free parameters: the black hole mass $M_{\mathrm{BH}}$, the stellar $H$-band $M/L$ ratio $\Upsilon_{H}$, the pixel location of the disk’s dynamical center ($x_{\mathrm{c}},y_{\mathrm{c}}$), the disk inclination angle $i$ and major-axis position angle $\Gamma$ of the receding side of the disk, the systemic velocity $v_{\mathrm{sys}}$, the turbulent velocity dispersion $\sigma_{\mathrm{turb}}$, and a flux-scaling factor $f_{0}$ to correct for possible flux normalization mismatch between the data and model. The gas velocity dispersion $\sigma_{\mathrm{turb}}$ can be set to a uniform (but freely varying) value over the disk surface, or allowed to vary as a function of radius with the introduction of additional free parameters. The models are calculated on a pixel grid that is oversampled by a factor of $s$ relative to the ALMA data cube pixel size of 0$\farcs$015 pixel*-1*. In other words, each ALMA spatial pixel is subdivided into an $s\times s$ grid of sub-pixel elements. For NGC 3258, we calculated initial models for values of $s$ ranging from 1 to 14.

A required input to this calculation is a model map of the disk’s CO surface brightness distribution prior to convolution by the ALMA synthesized beam. To generate this map, we collapsed the ALMA Cycle 4 data cube to form an image (see Figure 10), and applied ten iterations of the IRAF STSDAS Richardson-Lucy deconvolution (Richardson, 1972; Lucy, 1974) task LUCY using the elliptical Gaussian synthesized beam.

The disk’s circular velocity $v_{\rm c}(r)$ is calculated as a function of radius for rotation in the combined gravitational potential of the BH (a point source at $r=0$) and the host galaxy. The host galaxy contribution $v_{\mathrm{c}}^{\star}(r)$ to the circular velocity is determined using the host galaxy luminosity profiles derived from the dust-corrected MGE models, with these velocity values scaled by $\sqrt{\Upsilon_{H}}$. We assume a spatially uniform $M/L$ ratio in our model calculations. The optically thick dust disk within the inner kpc of NGC 3258 makes it difficult to constrain any possible $M/L$ gradient, but the three dust-corrected stellar luminosity profiles described in §2.3 correspond to a range of different central mass profile slopes that collectively encompass the possible effect of stellar $M/L$ variations. We do not include the gas disk itself in the mass model. In §3.2, we estimate the disk’s H2+He gas mass to be $\sim$108 $M_{\odot}$, and this gas mass is distributed in a disk extending out to $r\sim 150$ pc, within which the total enclosed mass is $\sim 5\times 10^{9}$ $M_{\odot}$. In effect, the gas disk’s small contribution to the $M(r)$ profile will be subsumed into the $M/L$ parameter, although there will be a small residual error since the disk’s radial mass profile differs from the stellar $M(r)$ profile. Because our final dynamical model (described in §4.2 below) determines the spatially extended mass profile directly from the kinematic data, independent of the host galaxy surface brightness profile measurements, that method incorporates all gravitating mass contributions that may be present.

For the turbulent velocity dispersion within the molecular gas disk, we adopt either a spatially uniform value across the disk surface ($\sigma_{\mathrm{turb}}=\sigma_{1}$), or an axisymmetric model allowing for radial variation in $\sigma_{\mathrm{turb}}$ with a Gaussian radial profile: $\sigma_{\mathrm{turb}}=\sigma_{0}\exp[-(r-r_{0})^{2}/2\mu^{2})]+\sigma_{1}$, where $\sigma_{0}$, $\sigma_{1}$, $\mu$, and $r_{0}$ are free parameters. We use $\sigma_{\mathrm{turb}}$ to represent the combination of processes contributing to the emergent line width of the disk: internal turbulence and rotation of individual clouds, as well as radial velocity variations between clouds contained within a given grid element, whether due to rotational shear in the disk or random cloud-to-cloud velocity variations. The molecular gas kinetic temperature in ETG circumnuclear disks is very low, $\sim$10$-$20 K (Bayet et al., 2013), so gas temperature makes a negligible contribution to the CO line widths. In the ALMA data cube, the minimum observed line dispersion is just $\sim$7 km s*-1*, while the central rise to $\sim$300 km s*-1* is likely almost entirely the result of beam smearing at small radii (see Barth et al., 2016b for a detailed discussion of this effect). The Gaussian $\sigma_{\mathrm{turb}}$ model allows for the possibility that some portion of this central increase in line width is intrinsic.

We populate the model cube at each spatial location with Gaussian emergent line profiles, defined by the projected line-of-sight velocity and $\sigma_{\mathrm{turb}}$ value at each oversampled grid point. The model cube spectral axis is observed frequency, and we transform rest-frame projected velocities and $\sigma_{\mathrm{turb}}$ maps to observed frequencies prior to creating the line profiles (Meyer et al., 2017). The line profile flux at each spatial location is determined using the model flux map. As the CO surface brightness distribution is not known on subpixel scales, each line profile within an $s\times s$ block corresponding to a single ALMA pixel is equally weighted in flux, such that the total $s\times s$ region contains the same total flux as the deconvolved CO flux model at that pixel location.

The two final steps of the model calculation are the convolution of each model cube channel with the synthesized beam, and averaging of each $s\times s$ block of oversampled pixels into a single pixel matching the scale of the ALMA data cube. In principle, for highest fidelity the beam convolution would be computed on the oversampled pixel grid. Beam convolution is the most time-consuming portion of the model calculation procedure, and for large values of $s$ this would become prohibitively slow. In fact, we found that the modeling results do not appreciably change if the model cube is first rebinned to the original pixel scale of the ALMA data prior to the beam convolution step since the synthesized beam is already oversampled by the chosen pixel size. We adopted this method in order to minimize the computational time required for model optimization.

Within each frequency channel, the background noise in the ALMA data cube contains strong pixel-to-pixel correlations on scales comparable to and smaller than the angular scale of the synthesized beam. Constructing the full covariance matrix to account for correlated errors in these data remains very challenging (e.g., Hezaveh et al., 2016; Onishi et al., 2017). Instead, we mitigate the impact of correlated errors in the $\chi^{2}$ calculation by first spatially block-averaging the data in $5\times 5$ pixel regions to form roughly beam-sized cells. We then measure the rms noise levels in line-free regions in each frequency channel of the block-averaged data cube. The channel-dependent background noise is somewhat lower than measured at the original pixel scale, averaging to 0.18 mJy beam*-1* in $\sim$10.2 km s*-1* channels. For each model iteration, the beam-convolved model cube is block-averaged in the same way as the data. After rebinning, we calculate $\chi^{2}$ within an elliptical spatial fitting region that includes nearly the entire disk area with major axis radius $r_{\mathrm{fit}}=1\arcsec$, PA=77°, and axis ratio $q=0.67$, and a spectral fitting region $2310\leq v_{\mathrm{LOS}}\leq 3273$ km s*-1* that spans a slightly larger velocity range than does the CO(2$-$1) emission (see Figure 11). The block-averaged fitting region used to compute $\chi^{2}$ contains 415 spatial pixels and 94 frequency channels, for a total of 39010 data points.

4.1.2 Flat disk model results

In initial modeling trials, we tested how parameter values change with increasing oversampling factor $s$. As shown in Figure 12, we find that the best-fitting BH mass converges to stable values for $s\geq 4$, while for smaller values of $s$ the best-fitting $M_{\mathrm{BH}}$ increases by at most $\sim$1%. Computation time increases dramatically for $s>4$ with very little change in the resulting BH mass, so for the remainder of this work all model calculations use $s=4$.

Our initial model C1 incorporated a flat disk and a spatially uniform turbulent velocity dispersion $\sigma_{\mathrm{turb}}(r)=\sigma_{1}$. The extended mass distribution was characterized by the initial $v_{\mathrm{c}}^{\star}(r)$ profile assuming no extinction in the disk ($A_{H}=0$). After optimizing to the CO(2$-$1) data cube, we obtained best-fit parameters $M_{\mathrm{BH}}=2.280\times 10^{9}$ $M_{\odot}$, $\Upsilon_{H}=2.72$ $M_{\odot}$/$L_{\odot}$, and $\sigma_{\mathrm{turb}}=10.5$ km s*-1* (see Table 3 for the complete results). The total $\chi^{2}=47929.2$ results in $\chi^{2}_{\nu}=1.229$ over $N_{\mathrm{dof}}=39001$ degrees of freedom. This basic dynamical model reproduces the general CO kinematic behavior moderately well, although quantitatively it does not constitute a good fit to the data. For this $N_{\mathrm{dof}}$, a formally acceptable fit should achieve $\chi^{2}_{\nu}\leq 1.012$.

At high angular resolution, the CO line structure in each channel map forms a tight locus of emission with a characteristic ”V” shape. We constructed a residual cube by subtracting the model from the data cube, and in each channel identify regions where the model is mismatched with the data: line structure discrepancies between data and model channel slices can be separated into those that arise from neglect of disk kinematic warping (i.e., rotational components that warp the ”V” shape) and those that stem from an inadequate host galaxy mass model (i.e., that shift and dilate the locus in the radial direction). In most channels, we find the discrepancies to be primarily rotational.

In models D1–D4, we adopted the more flexible Gaussian $\sigma_{\mathrm{turb}}(r)$ function and used each of the extinction-corrected $v_{\mathrm{c}}^{\star}$ profiles in turn to explore the impact of central dust extinction on $M_{\mathrm{BH}}$. The parameters for their best-fit model cubes converge to a range of values $M_{\mathrm{BH}}=(2.059-2.276)\times 10^{9}$ $M_{\odot}$ and $\Upsilon_{H}=2.04-2.73$ $M_{\odot}$/$L_{\odot}$. These best fits obtain minimum $\chi^{2}=47458.1-47547.7$ over $N_{\mathrm{dof}}=38998$, corresponding to an average $\chi^{2}_{\nu}\approx 1.22$ with slight preference for model D3 (corresponding to a central disk extinction of $A_{H}=0.75$ mag). Including a radially-varying $\sigma_{\mathrm{turb}}$ does improve the overall fit without significantly affecting the BH mass, as demonstrated by the decrease in $\chi^{2}$ from model C1 to D1. In model D3, $\sigma_{\mathrm{turb}}(r)$ reaches a peak of 18.0 km s*-1* at $r=12.2$ pc, decreasing to 8.3 km s*-1* at the disk edge.

As the extinction correction increases from $A_{H}=0$ to 1.50 mag, the corresponding $v_{\mathrm{c}}^{\star}$ profiles reflect increasing stellar luminosity density at all disk radii with a bias towards increasing nuclear contributions (see Figure 13). Since the total enclosed mass is tightly constrained by velocities at the outer edge of the disk, a more cuspy central stellar surface brightness slope arising from a higher assumed extinction has the effect of slightly lowering both the best-fit $M_{\mathrm{BH}}$ and $\Upsilon_{H}$ values. The highest of these $\Upsilon_{H}$ values measured using models D1 and D2 are elevated when compared to typical dynamical $H$-band $M/L$ ratios in other ETGs (e.g., Onishi et al., 2017; Yıldırım et al., 2017) while remaining consistent with those derived from stellar population synthesis modeling (e.g., Zibetti et al., 2009).

To visualize the model results, we show GH moments in Figure 10 and the PVD in Figure 11 that are measured from the best-fit model D3 cube in the same manner as the data. The flat disk model velocities closely agree with the observed velocity field for most of the disk, with typical residuals $|\Delta v_{\mathrm{LOS}}|\lesssim 3$ km s*-1*. The velocity peaks in the flat disk model are offset from the observed locations by nearly 0$\farcs$05 (in a clockwise direction about the disk center) with large associated residuals ranging between $-80$ and $+180$ km s*-1*, demonstrating the limitations of a flat disk formalism when modeling even mildly warped disks. In §4.5 below, we also explore the possibility that a non-circular component of the gas velocity may account for the central kinematic twists.

In Figure 14 we show $\Delta\chi^{2}$ curves as a function of fixed BH mass for models D1$-$D4. Assuming the usual $\Delta\chi^{2}\leq 9$ criterion, the $3\sigma$ (statistical) uncertainties in $M_{\mathrm{BH}}$ for a given host galaxy model would be less than 1% of $M_{\mathrm{BH}}$. For the preferred model D3, the nominal $1\sigma$ uncertainty obtained by $\Delta\chi^{2}\leq 1$ is estimated to be less than 0.2% of its best-fit $M_{\mathrm{BH}}$ value. The range in BH mass of $\Delta M_{\mathrm{BH}}=2.2\times 10^{8}$ $M_{\odot}$ (nearly 10% of the BH mass) for these four models with different $v_{\mathrm{c}}^{\star}$ profiles far exceeds the statistical bounds on any one of the four. This range is representative of the systematic uncertainty introduced by dust extinction.

We note that $\chi^{2}_{\nu}$ values from these fits to the data cube do not faithfully characterize the model fidelity, because block-averaging does not fully mitigate the correlations between neighboring pixels. Thus, we do not use the $\Delta\chi^{2}$ curves when determining the error budget on $M_{\mathrm{BH}}$; instead, we adopt a Monte Carlo resampling procedure (described in §4.4) to calculate the final statistical uncertainty.

For models D1$-$D4, the radius of the BH sphere of influence (defined as the radius within which the enclosed stellar mass is equal to $M_{\mathrm{BH}}$) is 131–143 pc, projecting to an angular size of 0$\farcs$86–0$\farcs$94.

4.1.3 Cycle 2 Comparison

Our Cycle 2 CO(2$-$1) imaging of NGC 3258 with $\theta_{\mathrm{FWHM}}=0\farcs 48\times 0\farcs 40$ provides $\sim$2 resolution elements across the BH radius of influence, so this initial data set should also allow for a confident BH mass measurement, although it will still be subject to the same uncertainty in the dust-disk extinction correction. Comparison with the Cycle 4 models provides an opportunity to test the impact of angular resolution on the best-fit $M_{\mathrm{BH}}$.

As described in Paper I, the Cycle 2 data cube has a pixel scale of 0$\farcs$04 and a rest-frame velocity channel width of $\sim$20.3 km s*-1* for CO(2$-$1) emission redshifted to the systemic velocity of NGC 3258. We fit the Cycle 2 data in models A1 and B1$-$B4 using procedures that correspond to Cycle 4 models C1 and D1$-$D4. We treated the Cycle 2 modeling in a self-contained manner by using a Richardson-Lucy deconvolution of the smoother Cycle 2 CO distribution to weight the model line profiles. We block-averaged both data and model cubes in $4\times 4$ pixel regions prior to calculating model goodness-of-fit. These new cell sizes are significantly smaller than the synthesized beam area but allow for many spatial cells across the disk. At this more coarse angular resolution, the slightly larger $r_{\mathrm{fit}}=1\farcs 2$ fitting region contains 124 spatial cells and 46 frequency channels, for a total of 5704 data pixels. Results of the Cycle 2 model fits are listed in Table 3, and Figure 14 shows the $\Delta\chi^{2}$ curves for models B1$-$B4 for comparison with the analogous Cycle 4 models D1$-$D4.

Overall, the Cycle 2 model fits yield $M_{\mathrm{BH}}$ values within a few percent of the values obtained from the analogous Cycle 4 models, and $\Upsilon_{H}$ values that are $\sim$20% greater than those from the corresponding Cycle 4 models. GH moments measured from the best-fitting model B1 cube show good agreement with those obtained from the data (see Figure 15). The Cycle 2 model fits also find low $\sigma_{\mathrm{turb}}$ with Gaussian line width coefficients similar to those obtained from the Cycle 4 data. From examination of fitting residuals in the data cube, we find large residuals near the disk center, which we attribute in part to insufficient resolution in the flux map used in the modeling procedure. The Cycle 2 data do not recover the central hole in CO(2$-$1) surface brightness, and as a result the model assigns too much flux to pixels at LOS velocities $>500$ km s*-1* in the innermost region of the disk, producing line profiles that exceed the maximum observed $|v_{\mathrm{LOS}}-v_{\mathrm{sys}}|$. The worsening $\chi^{2}$ from models B1 to B4 stems from the additional central stellar mass that is introduced by the increasingly dust-corrected $v_{\mathrm{c}}^{\star}$ mass models, thereby increasing the model rotation speed near the BH. For an individual host galaxy mass model, the $\Delta\chi^{2}$ curve is wider for the Cycle 2 data than for the corresponding Cycle 4 model fit, implying statistical uncertainties that are larger by a factor of $\sim$2 for the same host galaxy radial profile. (This analysis does not consider the larger $\chi^{2}_{\nu}$ values obtained for the Cycle 2 modeling due to significant correlated noise between block-averaged cells. However, even if we were to inflate the background rms noise to drive $\chi^{2}_{\nu}$ to unity, the $\Delta\chi^{2}$ criterion would not yield significantly broader confidence intervals for $M_{\mathrm{BH}}$.)

Despite these issues, the close agreement in $M_{\mathrm{BH}}$ between the Cycle 2 and Cycle 4 flat-disk model fits demonstrates that the Cycle 2 data already provide a good determination of $M_{\mathrm{BH}}$. For a fixed host-galaxy mass model, the improvement in ALMA angular resolution (from resolving $r_{\mathrm{g}}$ by a factor of $\sim$2 to a factor of $\sim$10) results in a relatively modest improvement in precision on $M_{\mathrm{BH}}$. In either case, the dominant uncertainty when using MGE-based mass models stems from the uncertainty in the dust extinction correction rather than from the model-fitting precision. It is important to note that these model fits are carried out over a spatial region that is almost entirely contained within $r_{\mathrm{g}}$ for NGC 3258. As a result, the uncertainty in the central stellar mass profile slope only results in a modest (several percent) uncertainty in $M_{\mathrm{BH}}$ even for the Cycle 2 data. In many other ALMA gas-dynamical targets, the molecular disk extends to radii well beyond $r_{\mathrm{g}}$. In such cases, if the model fits are carried out over the entire dust disk, the fitting results will tend to be dominated by the influence of the large fraction of spatial pixels well outside of $r_{\mathrm{g}}$, in which case the uncertainty in the dust extinction correction will likely lead to a far larger range of uncertainty in $M_{\mathrm{BH}}$ than what we find for NGC 3258.

4.2 Detailed Dynamical Modeling

In this section, we introduce more general descriptions for the disk structure and host galaxy mass distribution, with the addition of two additional features to the modeling procedure described above. First, we incorporate a tilted-ring model that fits the disk’s mildly warped structure more accurately than flat-disk models. Second, we employ a method to constrain the host galaxy’s spatially extended mass profile solely through fitting to the ALMA CO kinematics, rather than relying on the $HST$ imaging data (and an uncertain extinction correction) to constrain the host galaxy model. These two improvements are made possible by the angular resolution of the ALMA Cycle 4 observations, which fully resolve the rotational structure of the disk without significant blurring of the central kinematics by beam smearing.

4.2.1 Tilted-ring model

As shown in Figure 10, the $v_{\mathrm{LOS}}$ map for flat-disk model D3 does not fully reproduce the observed velocity field due to the mild kinematic twist near the disk center discussed in §3.2. Given the high precision enabled by the resolution of our Cycle 4 observation, it is important to determine how the disk’s warp may affect the inferred BH mass. In models E1 and F1, we implement a tilted-ring model that characterizes a warped (but still intrinsically thin) disk using a series of concentric rings (e.g., Rogstad et al., 1974), with each ring allowed to have an arbitrary PA $\Gamma$ and inclination angle $i=$ arccos q. The model comprises ten rings spanning the disk’s radius (see Table 4 for the ring radii), approximately matching the number of beam widths across $r_{\mathrm{fit}}$. The non-linear spacing between annuli was chosen to better characterize the more abrupt increases in $\Delta\Gamma$ and $\Delta q$ toward the disk center.

For each model iteration, we form continuous $\Gamma(r)$ and $q(r)$ functions by a cubic spline interpolation of these ring parameters at intermediate radii and construct 2D maps of intrinsic radius and inclination at each projected disk location on the plane of the sky. One-to-one correspondence between a projected and physical disk location is not in general preserved for warped disks (especially for rapid shifts in $\Gamma$ and $q$ shifts and for more edge-on disks; see Corbelli & Schneider, 1997; Józsa et al., 2007; Davis et al., 2013). However, this approximation is suitable for the moderately inclined and warped gas disk in NGC 3258. We proceed to generate model line profile cubes using the maps of intrinsic $v_{\mathrm{LOS}}$ and $\sigma_{\mathrm{LOS}}$. Beam smearing and subsequent model fitting steps are applied as described previously. We allow the values of $\Gamma$ and $q$ for each ring to vary freely while simultaneously optimizing other disk parameters. Model E1 incorporates the tilted-ring disk structure and employs the same MGE-based host galaxy profile as model D3 (with $A_{H}=0.75$ mag), while model F1 incorporates both the tilted-ring disk and the new host galaxy mass modeling procedure described below.

4.2.2 Host galaxy mass profile from CO kinematics

The four MGE-derived $v_{\mathrm{c}}^{\star}$ profiles are corrected for a plausible range of central disk extinction levels. In models D1$-$D4, we adopted each profile in turn to explore the impact on the derived $M_{\mathrm{BH}}$ of the range in possible host galaxy mass distributions. The optimized models are very similar in a statistical sense yet yield best-fit $M_{\mathrm{BH}}$ values that span a mass range of about 10%. We cannot determine the correct (average) extinction level using these $v_{\mathrm{c}}^{\star}$ profiles alone, and the associated systematics would make the Cycle 4 modeling results nearly as uncertain as those from the Cycle 2 data set.

Fortunately, our Cycle 4 data are so well resolved that we can constrain the galaxy mass profile directly by modeling the CO kinematics, without reference to the NIR imaging data. This is the only method that can potentially reduce the systematic uncertainties to the percent level (or better) for such dusty disks, because host galaxy models based on image deprojections will always suffer from systematic uncertainty in the extinction correction. We refer to the circular velocity profile arising from the extended mass distribution as $v_{\mathrm{ext}}(r)$. This velocity profile primarily reflects the stellar mass distribution but also includes any other gravitating mass, including the gas disk itself (see Figure 17) as well as the expected very small contribution of dark matter (extrapolation from observations and simulations of luminous ETGs suggests a dark matter mass of less than $\sim 10^{7}$ $M_{\odot}$ enclosed within NGC 3258’s central dust disk region; e.g., Newman et al., 2013; Wang et al., 2018).

In model F1, we describe the extended mass distribution in terms of a circular velocity profile with ten free parameters, where the free parameters correspond to the values of $v_{\mathrm{ext}}$ at the same set of ten ring radii used to generate the tilted-ring model. We create the model velocity field by cubic spline interpolation of $v_{\mathrm{ext}}$ between the rings to determine its value at each disk location, afterwards calculating the disk rotation speed at each position resulting from both the BH mass and the spatially extended mass, and finally projecting the rotation speed along the line of sight using the tilted-ring model to calculate $v_{\mathrm{LOS}}$ at each spatial pixel in the model. The ten free $v_{\mathrm{ext}}$ parameters were optimized simultaneously with the tilted-ring angular parameters and the other disk parameters, for a total of 39 free parameters in the final model. The only constraints we applied to the circular velocity model is that $v_{\mathrm{ext}}(r)$ was required to be an increasing function of radius and that $v_{\mathrm{ext}}=0$ at $r=0$. This method of determining $v_{\mathrm{ext}}$ is largely equivalent to allowing a radially-varying $M/L$ ratio when scaling the stellar luminosity profile (Davis & McDermid, 2017; Davis et al., 2018). However, our method eliminates any dependence on the luminosity profile derived from imaging data, instead allowing nearly complete freedom in the $M(r)$ profile to match the kinematic data.

4.3 Detailed modeling results

We first optimize model E1, which includes a tilted-ring geometry and an extinction-corrected ($A_{H}=0.75$ mag) galaxy mass distribution. Aside from the flexible disk structure, this scenario is identical to the D3 case, making it possible to isolate the impact of disk warping on the derived BH mass. The optimized model converges to $M_{\mathrm{BH}}=2.203\times 10^{9}$ $M_{\odot}$ and $\Upsilon_{H}=2.18$ $M_{\odot}$/$L_{\odot}$. The total $\chi^{2}=46009.9$ yields $\chi^{2}_{\nu}=1.180$ over $N_{\mathrm{dof}}=38980$ and represents the most substantial fit improvement for the Cycle 4 gas-dynamical models; in contrast, model D3 achieved $\chi^{2}_{\nu}=1.217$. The tilted-ring angular parameters smoothly increase by $\Delta\Gamma\approx 20\arcdeg$ and $\Delta q\approx 0.27$ (corresponding to an inclination angle decrease $\Delta i\approx-26\arcdeg$) towards the disk center. The shift to a more face-on disk orientation at small radii projects circumnuclear nuclear rotation away from the line of sight and results in a moderate $\Delta M_{\mathrm{BH}}=5.9\times 10^{7}$ $M_{\odot}$ (or $\sim$3%) increase in $M_{\mathrm{BH}}$ relative to the otherwise analogous flat-disk model D3.

We go on to investigate the host galaxy mass profile in model F1, which is identical to E1 except that it employs the freely varying $v_{\mathrm{ext}}$ method to represent the host galaxy mass distribution instead of the MGE-based host galaxy model. In this case, the BH mass converges to $2.249\times 10^{9}$ $M_{\odot}$, and the best fit achieves $\chi^{2}=45956.4$ over $N_{\mathrm{dof}}=38971$ for $\chi^{2}_{\nu}=1.179$. The $\chi^{2}_{\nu}$ statistic decreases only slightly from model E1 to F1, indicating that between models D3 and F1 most of the improvement in fit quality was the result of including the tilted-ring disk geometry (detailed in Figures 17 and 18) rather than allowing additional freedom in the host galaxy model. However, the primary advantage of the freely varying $v_{\mathrm{ext}}$ host galaxy model is that it removes the systematic uncertainty in $M_{\mathrm{BH}}$ resulting from dust extinction that is inherent in the MGE-based host galaxy models. Model F1 attains the lowest $\chi^{2}_{\nu}$ value of any of our trial models, and we consider it to be our final preferred model for the NGC 3258 disk.

As a result of beam smearing of the central disk kinematics, the kinemetry measurements $\Gamma_{k}$ and $q_{k}$ do not trace the intrinsic disk structure as faithfully within the inner couple of beam widths. In particular, within the innermost 0$\farcs$2, the strong intrinsic change in line-of-nodes PA implied by the tilted-ring model exceeds the $\Gamma_{k}$ rise, and the axis ratio $q_{k}$ approaches unity near the nucleus while the tilted-ring model axis ratio reaches a central value of $\approx 0.89$ for model F1 (see Figure 17). At $R>0\farcs 2$, beam smearing has less impact on the observed kinematics, and $\Gamma_{k}$ and $q_{k}$ more closely trace the $\Gamma(r)$ and $q(r)$ profiles of the tilted-ring model.

In Figure 10 we show GH moment and residual maps measured from this best model F1 cube. We find the velocity residuals $\Delta v_{\mathrm{LOS}}$ are generally small and centered about zero, with $\sim$60% of the spatial pixels in the model falling within $\pm 2$ km s*-1* of the observed velocity in the data. The tilted-ring disk model alleviates most of the large central discrepancies apparent in the flat-disk velocity map generated from model D3, although the $|v_{\mathrm{LOS}}-v_{\mathrm{sys}}|$ asymmetry between the quasi-Keplerian peaks still leads to $\sim$50 km s*-1* velocity residuals for $R<0\farcs 15$ in model F1.

The data$-$F1 $\Delta\sigma_{\mathrm{LOS}}$ residuals are also substantial at pixels near the nucleus, and we find the Gaussian $\sigma_{\mathrm{turb}}$ profile underpredicts the line widths directly north and south of the nucleus by $\sim$100 km s*-1*. These locations are coincident with bright clumps of CO emission, and we expect these non-axisymmetric excess line width features either to be correlated with gas cloud size (e.g., Shetty et al., 2012) or to result from strong tidal shear within $\sim$20 pc of the BH. At the nucleus, this model overpredicts the data $\sigma_{\mathrm{LOS}}$ by over 200 km s*-1*, and in §4.5 we discuss how this feature may be the result of an inadequate CO surface brightness map.

For the best-fitting model F1, the BH radius of influence (defined as the radius within which the enclosed stellar mass is equal to $M_{\mathrm{BH}}$) is $r_{\mathrm{g}}=143$ pc, or 0$\farcs$94.

4.4 Monte Carlo error analysis

Although our models match the overall kinematic structure of the disk well in general, the final $\chi^{2}_{\nu}=1.179$ for model F1 indicates that the model is not formally an acceptable fit to the data for 38971 degrees of freedom. In this situation, determining the statistical uncertainty on $M_{\mathrm{BH}}$ and other parameters by examining contours in $\Delta\chi^{2}$ would tend to underestimate the true uncertainty range. For example, measuring $\chi^{2}$ as a function of fixed $M_{\mathrm{BH}}$ while allowing all the other model parameters to freely vary results in a very narrow $\Delta\chi^{2}$ curve, with the $\Delta\chi^{2}=9$ range (the nominal $3\sigma$ uncertainty range) corresponding to $\pm 0.2\%$ of the best-fitting BH mass (see Figure 19), and $\Delta\chi^{2}=1$ (for a $1\sigma$ uncertainty range) corresponding to just $\sim 0.06\%$ of $M_{\mathrm{BH}}$, although the bottom of the $\Delta\chi^{2}$ curve for model F1 is slightly irregular.

To obtain a more robust measure of the statistical model-fitting uncertainties in this situation, we carried out 100 Monte Carlo realizations of the best-fit model F1 cube. To add realistic noise to this model cube, we used noise from line-free channels in the continuum-subtracted ALMA data cube itself. We extracted nearly 100 line-free channels from the data cube where $|v_{\mathrm{LOS}}-v_{\mathrm{sys}}|>500$ km s*-1*, and randomly assigned and added these noise slices to the model cube channels at each resampling iteration. After incorporating this realistic noise, we carried out complete model fits to each noise-randomized model cube following the same procedure as for model F1, including both the tilted-ring model and the flexible $v_{\mathrm{ext}}$ description. All model parameters were allowed to freely vary. From this suite of Monte Carlo realizations, we determined $1\sigma$ uncertainties on each of the thirty-nine model parameters by taking the standard deviation of the set of their best-fit values.

We include the histogram of $M_{\mathrm{BH}}$ values determined from this procedure in Figure 19 to compare to the $\Delta\chi^{2}$ results. While somewhat broader than the $\chi^{2}$ bounds, the distribution of $M_{\mathrm{BH}}$ values remains very narrow. We adopt the standard deviation $3.8\times 10^{6}$ $M_{\odot}$ (corresponding to $\sim$0.17% of the BH mass) of these Monte Carlo results as the final $1\sigma$ statistical uncertainty on $M_{\mathrm{BH}}$. Tables 3 and 4 list the full set of parameter uncertainties for model F1 based on these Monte Carlo simulations.

4.5 Additional tests and error budget

We now describe additional tests conducted to estimate the systematic uncertainties on $M_{\mathrm{BH}}$. In each test, we modified aspects of model F1 to explore the sensitivity of our model-fitting results to various details of the model construction.

Pixel oversampling and block averaging: We adopted an oversampling factor $s=4$ for the model fits described above, based on the results shown in Figure 12. Ionized gas disk dynamical modeling has demonstrated a typical scatter in derived $M_{\mathrm{BH}}$ values of a few percent for different $s$ values (e.g., Barth et al., 2001), behavior that may also apply to ALMA data (Barth et al., 2016b). Model F1 tests show that oversampling factors $s<4$ do not adequately sample the velocity field, resulting in a $\sim$1.3% decrease in BH mass from $s=1$ to 4 that we do not include in the final error budget. The $M_{\mathrm{BH}}$ results are very stable for $s\geq 4$, with the best-fit BH mass decreasing by $\Delta M_{\mathrm{BH}}=-1.5\times 10^{6}$ $M_{\odot}$ (corresponding to $\sim$0.07% in BH mass; see Figure 12) as $s$ increases to 14.

As described in §4.1.1, the Cycle 4 data and model cubes were spatially block-averaged into $5\times 5$ pixel bins prior to computing $\chi^{2}$ for each model iteration, in order to mitigate the impact of correlated noise on spatial scales smaller than the ALMA beam width. We also explored block-averaging the data and model cubes using pixel regions ranging from $1\times 1$ (no averaging) to $10\times 10$ and found only a negligible impact on the derived $M_{\mathrm{BH}}$, with the best-fit models converging to a narrow range of BH masses with a scatter of $\Delta M_{\mathrm{BH}}=\pm 1.1\times 10^{6}$ $M_{\odot}$ about the model F1 value.

Tilted-ring model: The choice of ten rings to anchor our tilted-ring model was somewhat arbitrary but appears sufficient to recover the disk structure. Model parameters may be sensitive to the number of rings $N_{\mathrm{rings}}$ that define the warped disk, so we explored different annular spacing from single $\Gamma$ and $q$ values (a flat disk) to $N_{\mathrm{rings}}=20$. To isolate the effect of changing the number and spacing of rings in the warped disk model, $v_{\mathrm{ext}}$ is still optimized at the initial ten radial locations. As shown in Figure 20, increasing the number of rings does improve the overall fit but the BH mass is not significantly impacted for $N_{\mathrm{rings}}\geq 5$. When using between ten and twenty rings, the best-fit BH masses span a range of only $3.7\times 10^{6}$ $M_{\odot}$ ($\sim$0.16% in BH mass). For $N_{\mathrm{rings}}\geq 10$, the tilted-ring solutions return consistent, small-amplitude oscillations in $\Gamma(r)$ and $q(r)$ (of $\sim$2° in both PA and $i$; see Figures 18 and 20) for radii $r>50$ pc.

Fitting region: In the models described in Tables 2 and 3, we measured $\chi^{2}$ by fitting to essentially the entire disk. However, our symmetric models cannot fully account for local irregularities in the velocity and velocity dispersion fields. Here, we highlight the most apparent discrepancies and, by adjusting the model fitting region, estimate their potential impact on our dynamical modeling results.

The fitting region for model F1 (and all other models in Table 2) gives roughly equal weight to the red and blueshifted portions of the inner disk, even though the molecular gas within $R<0\farcs 2$ on the approaching side of the disk appears to be in sub-Keplerian rotation. Assuming that the approaching velocity peak represents one of these local irregularities, we explored its impact on model results by excluding the affected data: we restricted the fitting region on the approaching side to channels where $|v_{\mathrm{LOS}}-v_{\mathrm{sys}}|<350$ km s*-1*. The fitting region is otherwise unchanged, and this test retains the full generality of model F1. After optimizing to the data cube, we find only a small BH mass increase $\Delta M_{\mathrm{BH}}=2.6\times 10^{6}$ $M_{\odot}$ relative to the model F1 results. Excluding channels with obviously asymmetric gas rotation reduces the number of data points by nearly 11% while only decreasing the number of cells containing CO emission by just 2%. As a result, this adjustment to the fitting region produces only a small change in $M_{\mathrm{BH}}$.

Due to the abundance of data points at larger radii, the full fitting region gives greater weight to data near the disk edge than near the BH. We explored the impact of our choice of fitting radius by calculating a model with $r_{\mathrm{fit}}=0\farcs 5$ ($\sim$75 pc), and fitting to the same range of velocity channels as model F1. This spatial region extends to the edge of the observed central upturns in CO rotation speed and includes gas that is maximally sensitive to the dynamical influence of the BH. We optimized the tilted-ring and $v_{\mathrm{ext}}$ models only out to the first ring location beyond the new $r_{\mathrm{fit}}$ (at $r\sim 0\farcs 54$). The final BH mass increases by $\Delta M_{\mathrm{BH}}=3.9\times 10^{6}$ $M_{\odot}$ ($\sim$0.17%) relative to model F1. This change in BH mass is so small in part due to the radial flexibility of $v_{\mathrm{ext}}$. Using model D3 with an MGE-derived host-galaxy mass profile for comparison, adopting this same $r_{\mathrm{fit}}=0\farcs 5$ during model optimization induces a larger $\sim$0.5% relative increase in its best-fit BH mass.

The central CO(2$-$1) line widths in the best-fit model F1 cube are significantly discrepant with the data, as seen in the large $\Delta\sigma_{\mathrm{LOS}}$ values adjacent to the nucleus along the disk minor axis (Figure 10). We considered the impact of these local line width excesses on modeling results by excluding spatial locations where $\Delta\sigma_{\mathrm{LOS}}>25$ km s*-1* across all channels. Not surprisingly, this $\sim$3% decrease in $N_{\mathrm{dof}}$ produces a much improved overall model fit with $\chi^{2}_{\nu}=1.146$. However, the BH mass only increases by about 0.02%, so we do not expect these local line-width irregularities to cause significant error in $M_{\mathrm{BH}}$.

The first two adjustments to the fitting region both produce $\Delta M_{\mathrm{BH}}$ changes that are commensurate with the model F1 BH mass statistical uncertainty. To understand the significance of these shifts, we applied the same Monte Carlo error analysis to the best-fit test model cubes, subject to the respective changes to the fitting region. The resulting $1\sigma$ statistical uncertainties are $\sim$1.2$\times 10^{7}$ $M_{\odot}$ (roughly 0.5%) in $M_{\mathrm{BH}}$ for each test. In the first case, the larger BH mass statistical uncertainty is driven by more poorly constrained $\Gamma$, $q$, and $v_{\mathrm{ext}}$ values for $r<30$ pc; in the second case, it arises due to less certainty in $v_{\mathrm{sys}}$ and the $\sigma_{\mathrm{turb}}$ parameters. These tests demonstrate that, irrespective of the elevated $\chi^{2}_{\nu}$ values, our model fits to the totality of the disk yield an $M_{\mathrm{BH}}$ measurement that is insensitive to locally irregular kinematics. Figure 11 illustrates the good agreement between the observed and modeled PVDs everywhere except the approaching velocity peak for $|v_{\mathrm{LOS}}-v_{\mathrm{sys}}|\sim 400$ km s*-1*.

Central CO hole: Model F1 overpredicts the CO line widths at the nucleus, with data-model $\Delta\sigma_{\mathrm{LOS}}$ residuals falling below $-200$ km s*-1* in the central pixel. The most simple explanation is that low S/N nuclear CO emission may produce high-velocity line wings that remain buried beneath the noise and are therefore not reflected in the observed line widths. Another plausible explanation is that the deconvolved model CO flux map contains excess surface brightness at the disk center, overproducing unresolved high-velocity emission at the nucleus that translates to high model $\sigma_{\mathrm{LOS}}$ values. To test this possibility, we set the intrinsic model CO flux to zero within the synthesized beam area centered on the nucleus before again optimizing the model cube. In this case, the model cube line widths measured from the best-fit model decrease by $\sim$50 km s*-1* with a slight overall improvement in the fit (to $\chi^{2}_{\nu}\sim 1.177$). However, setting the central CO hole surface brightness to zero only increases the best-fit $M_{\mathrm{BH}}$ by 0.01%.

Radial motion: Regardless of their formation method (e.g., see Lauer et al., 2005; Davis et al., 2011; Martini et al., 2013), circumnuclear disks experience both secular evolution (Davis et al., 2018) and ongoing gas accretion (albeit perhaps at very low levels; van de Voort et al., 2015) that may result in detectable deviations from purely circular rotation. The relaxed molecular gas kinematics and regular dust disk morphology do not suggest any recent disruptions to the NGC 3258 molecular gas disk, although its mildly warped structure may indicate an ongoing settling process or perturbations arising from a triaxial galaxy potential (e.g., Emsellem et al., 2011).

We first adapted model F1 to include a spatially-uniform radial velocity term $v_{\mathrm{rad}}$ as a free parameter to simulate either bulk gas inflow or outflow. The radial flow component is projected along the line of sight and added to the projected tangential speed, which we approximate with the circular speed for the assumed $M_{\mathrm{BH}}$ value and the host-galaxy mass model. While not a self-consistent model of disk rotation, we simultaneously optimized $v_{\mathrm{rad}}$ with the other free parameters in this toy model to see how much radial flow is kinematically allowed by the data. This initial test favors bulk inflow with a speed of just 0.85 km s*-1* while the BH mass converges to the exact model F1 value.

Radial flows introduce twists in the line of nodes of otherwise circularly-rotating disks (e.g., see the analogous protoplanetary disk modeling of Walsh et al., 2017, Figure A4). Because the kinematic twists appear strongest within the inner $\sim$40 pc (see Figure 9), we tested whether the kinematics within this region might be consistent with flat-disk rotation and a higher inflow speed. After adding a radial flow component to model D3 with $v_{\mathrm{rad}}$ applying just to pixels at $R<0\farcs 26$, we find that $v_{\mathrm{rad}}$ converges to an inflow speed of 41 km s*-1*. This test largely reproduces the kinematics within this central region and achieves an overall $\chi^{2}=47335.7$ with $N_{\mathrm{dof}}=38997$ for $\chi^{2}_{\nu}=1.214$, which is a modest improvement over the original $\chi^{2}_{\nu}=1.217$ for model D3. While the above $v_{\mathrm{rad}}$ value approaches 10% of the circular rotation speed at these radii, the best-fit BH mass of $2.213\times 10^{9}$ $M_{\odot}$ is only about 0.2% lower than the corresponding best-fit $M_{\mathrm{BH}}$ in Table 3.

We then adopted $v_{\mathrm{rad}}$ as a free parameter for $R<0\farcs 26$ in our model F1 framework. After simultaneously optimizing all 40 free parameters, we find a global minimum with $v_{\mathrm{rad}}\sim 0$ km s*-1* while the remaining model parameters converge to the fiducial values given in Tables 3 and 4. Since a radial flow component can reproduce some of the apparent kinematic twists that arise from a warped disk, we anticipated significant degeneracy between $v_{\mathrm{rad}}$ and the $\Gamma$ and $q$ parameters for at least the inner ring positions. After setting the initial inflow speed guess to 40 km s*-1*, the model F1 variant settles on a local minimum where $v_{\mathrm{rad}}=26$ km s*-1* and the $\Gamma$ and $q$ parameters remain below 80° and 0.74, respectively. This local minimum achieves a slightly worse $\chi^{2}_{\nu}$ of 1.180 and returns a BH mass of $2.236\times 10^{9}$ $M_{\odot}$ that is only 0.7% lower than reported for the original model F1 in Table 3.

Finally, to rule out any significant impact of radial gas motion on the BH mass measurement, we again incorporated a bulk flow term $v_{\mathrm{rad}}$ into model F1 but only fit the model to points where $R>0\farcs 26$, thereby focusing on the region with the lowest disk warping to minimize possible degeneracies. The $\Gamma$, $q$, and $v_{\mathrm{ext}}$ parameters for the first four ring positions are fixed to the values in Table 4. We find that $v_{\mathrm{rad}}$ settles on an inflow speed of 0.86 km s*-1* while the BH mass converges to $2.247\times 10^{9}$ $M_{\odot}$, which corresponds to a mass difference of $\sim$0.1% from the fiducial value. After Monte Carlo resampling the resulting best-fit model cube, the distribution of $v_{\mathrm{rad}}$ values suggests that the possible detection of bulk radial inflow in the outer disk region is not particularly significant, being only $1.3\sigma$ removed from the $v_{\mathrm{rad}}\sim 0$ km s*-1* case. Since the kinematic twists in the CO velocity field appear to arise almost entirely from an intrinsically warped inner disk and not from gas inflow, we do not include any $\Delta M_{\mathrm{BH}}$ from this radial flow analysis in the final error budget.

Our conclusion of a low inflow rate within the CO disk is consistent (modulo an assumption of a steady flow) with evidence of a low inflow rate on smaller scales. If we assume an average inflow speed of just 1 km s*-1* (a level that is dynamically unimportant for our BH mass measurement), the entire circumnuclear disk with a radius of $\sim$150 pc would accrete onto the BH in about 150 Myr. For a total gas mass of $\sim$108 $M_{\odot}$, the average mass accretion rate over this accretion timescale is about 0.7 $M_{\odot}$ yr*-1*. This in turn translates to a ratio of BH mass accretion to the Eddington limit of $\dot{M}/\dot{M}_{\mathrm{EDD}}\sim 0.014$ (assuming a standard radiative efficiency $\epsilon=0.1$; van de Ven & Fathi, 2010), which would imply an accretion luminosity of $L_{\mathrm{bol}}\sim 10^{45.6}$ erg s*-1*. We do not see evidence for luminous AGN activity in $HST$ imaging, optical spectroscopy (Jones et al., 2009), or molecular gas outflows, suggesting that any modest inflow of molecular gas within the CO disk neither reaches the nucleus nor is directed out – consistent with negligible if any inflow at all.

Turbulence: For gas-dynamical modeling of some ionized gas disks, the intrinsic line widths are a substantial fraction of the disk rotation speed, suggesting significant local turbulence that is generally presumed to provide pressure support to the disk (Verdoes Kleijn et al., 2000; Barth et al., 2001; Walsh et al., 2013). In these cases, models based on purely circular rotation will underestimate the true BH masses, because the disk rotation velocity will lag behind the circular velocity (analogous to asymmetric drift in stellar dynamics). In thin-disk models that neglect this asymmetric drift effect, the fractional bias in $M_{\mathrm{BH}}$ is expected to be of order $(\sigma_{\mathrm{turb}}/v_{\mathrm{c}})^{2}$.

Our gas-dynamical models assume a perfectly thin and dynamical cold disk within NGC 3258, and do not account for the dynamical effect of turbulent pressure support. For the best-fit model F1, $\sigma_{\mathrm{turb}}/v_{\mathrm{c}}$ reaches a maximum of 0.037 at $\sim$50 pc from the disk center (see Figure 21) and a mean value $\langle\sigma_{\mathrm{turb}}/v_{\mathrm{c}}\rangle=0.030$ averaged over the disk surface. This molecular gas disk is truly dynamically cold. Since the fractional change to the BH mass resulting from turbulent pressure support scales as $(\sigma_{\mathrm{turb}}/v_{\mathrm{c}})^{2}$, we expect an upward correction to $M_{\mathrm{BH}}$ of order $\sim 3\times 10^{6}$ $M_{\odot}$ (corresponding to $\sim$0.14%) that is similar to the statistical model-fitting uncertainty.

Distance Uncertainty: Since the enclosed mass in the rotating disk model scales as $M(r)=rv_{\mathrm{c}}^{2}/G$, the inferred BH mass should in principle be directly proportional to the assumed angular size distance, although in practice other modeling details such as beam smearing may slightly modify this dependence. We anticipate that the uncertainty in the galaxy’s adopted distance $D_{L}=31.9$ Mpc of slightly more than 10% will introduce a commensurate systematic uncertainty in BH mass. We quantify this uncertainty by calculating two test models with the luminosity distance shifted by $\pm 1\sigma$ from the assumed value (i.e., $D_{L}=35.8$ and 28.0 Mpc, corresponding to angular scales of 170 and 133 pc arcsec*-1*, respectively). After optimizing over all free parameters, we obtain best-fit BH masses that are $\Delta M_{\mathrm{BH}}=\pm 2.7\times 10^{8}$ $M_{\odot}$ (or about 12%) removed from the fiducial model F1 value. We note that the uncertainty in our assumed NGC 3258 $D_{L}$ value does not include any systematic contributions that arise from Cephied period-luminosity metallicity corrections or uncertainties in the zero point (that are of order $\sim$0.1 mag; e.g., Mei et al., 2005; Blakeslee et al., 2010).

Some estimates of NGC 3258’s distance disagree with the ground-based SBF measurent from Tonry et al. (2001) by more than its quoted $1\sigma$ errors. Using $HST$ observations to measure SBFs in this galaxy, Cantiello et al. (2005) determined $m-M=33.00\pm 0.15$ mag, corresponding to $D_{L}=39.8\pm 2.8$ Mpc, although their analysis lacked empirical calibration of the SBF method in the F814W filter (Blakeslee et al., 2010). Using an angular scale of 189.5 pc arcsec*-1* derived for this second SBF distance, the best-fit BH mass increases by $\Delta M_{\mathrm{BH}}=5.6\times 10^{8}$ $M_{\odot}$ (or about 25%) from the model F1 case. Other distance measurement techniques yield distance modulii between $32.42\pm 0.19$ (or $D_{L}=30.5\pm 2.8$ Mpc, using the globular cluster luminosity function; Bassino et al., 2008) and $33.73\pm 0.41$ mag (or $D_{L}=55.7\pm 11.6$ Mpc, using the Fundamental Plane; Blakeslee et al., 2002), with respective ($\Delta M_{\mathrm{BH}}/M_{\odot}$) of $-5.9\times 10^{7}$ and $1.6\times 10^{9}$ from our model F1 results. We report a $\pm 12\%$ systematic distance uncertainty in the BH mass based on the reported SBF distance uncertainty from Tonry et al. (2001), but the $M_{\mathrm{BH}}$ may plausibly lie in the range $(2.0-3.8)\times 10^{9}$ $M_{\odot}$ based on these other distance estimates. Thus, while our model fits provide a highly precise determination of $M_{\mathrm{BH}}$ given an assumed distance to NGC 3258, the galaxy distance uncertainty dominates the total $M_{\mathrm{BH}}$ error budget.

As a final note on distance uncertainties, the preceding calculations have not accounted for source or observer peculiar velocities. Ideally, line-of-sight velocities and line width maps are transformed into observed frequency units assuming separate cosmological and peculiar redshifts in the relationship $(1+z_{\mathrm{obs}})=(1+z_{\mathrm{cos}})(1+z_{\mathrm{pec}})$, with the angular size distance depending on $z_{\mathrm{cos}}$ and not $z_{\mathrm{obs}}$. To investigate the impact on our $M_{\mathrm{BH}}$ determination from this neglect of peculiar motion, we first removed the Sun’s peculiar velocity contributions by transforming the Cycle 4 data into the cosmic microwave background (CMB) frequency reference frame wherein $z_{\mathrm{obs}}=0.010283$. Our adopted $D_{L}$ for this galaxy corresponds to $z_{\mathrm{cos}}=0.007745$ (Wright, 2006, using $H_{0}=73.24$ km s*-1* Mpc*-1*; Riess et al., 2018), which translates to $1\arcsec=152.3$ pc and a line-of-sight Doppler shift $v_{\mathrm{pec}}\approx 760$ km s*-1* for NGC 3258 in the CMB frame. Then, we fixed this $z_{\mathrm{cos}}$ value in a test model while allowing NGC 3258’s peculiar velocity $v_{\mathrm{pec}}$ to vary as a free parameter in place of $v_{\mathrm{sys}}$. This test converges to $v_{\mathrm{pec}}=753$ km s*-1* with a BH mass decrease of $\Delta M_{\mathrm{BH}}=-1.3\times 10^{7}$ $M_{\odot}$ from our model F1 result. In light of this galaxy’s disparate distance estimates, we did not attempt to separate out its cosmological and peculiar redshift contributions in models A$-$F, and we do not consider peculiar velocity systematics in the final BH mass error budget.

Final error budget: The statistical uncertainties on $M_{\mathrm{BH}}$ are equivalent to the largest model-dependent systematic terms, while the distance uncertainties are much larger than either of these other terms. Given the wide range of relative contributions, we separated these into distinct statistical (stat), model systematic (mod), and distance systematic (dist) terms in the final BH mass error budget. To estimate the total model systematic uncertainty, we separately combined in quadrature the positive and negative $\Delta M_{\mathrm{BH}}$ contributions listed above, with the largest of these (non-distance) systematics being at the 0.2% level. Our final BH mass with $1\sigma$ uncertainty ranges is then $(M_{\mathrm{BH}}/10^{9}\,M_{\odot})=2.249\pm 0.004$ (stat) ${}^{+0.007}_{-0.004}$ (mod) $\pm 0.270$ (dist).

5 Discussion

5.1 BH Mass

NGC 3258 has no previous BH mass measurement to compare with our gas-dynamical modeling results. Using this galaxy’s $\sigma_{\star}$ and $M_{K}$ values and uncertainties listed in Section 1, standard $M_{\mathrm{BH}}-\sigma_{\star}$ and $M_{\mathrm{BH}}-L_{K}$ relations for classical bulges and elliptical galaxies (Kormendy & Ho, 2013) predict ($M_{\mathrm{BH}}/10^{9}$ $M_{\odot}$) values of $(0.62^{+0.43}_{-0.23})$ and $(1.00^{+0.18}_{-0.16})$, respectively. Our NGC 3258 BH mass of $2.249\times 10^{9}$ $M_{\odot}$ is more than a factor of two larger than these predictions and lies on the extreme edge of measurements populating the $M_{\mathrm{BH}}-\sigma_{\star}$ and $M_{\mathrm{BH}}-L_{K}$ relations. Significant tension between prediction and measurement remains when employing a different univariate correlation (see also van den Bosch et al., 2016; Saglia et al., 2016), or after accounting for the impact of distance uncertainty on $M_{\mathrm{BH}}$ and $L_{K}$.

Quiescent BCGs and BGGs often exhibit cored surface brightness profiles (Lauer et al., 2007a; Rusli et al., 2013a), presumably formed through scouring by massive binary BHs (e.g., Ravindranath et al., 2002; Thomas et al., 2014). Even a partial depletion of their stellar core will suppress $\sigma_{\star}$ measurements for these luminous galaxies relative to $M_{\mathrm{BH}}-\sigma_{\star}$ extrapolations for normal ETGs (Lauer et al., 2007b). A more reliable indicator for cored galaxies is the break radius $r_{\mathrm{b}}$, which is found to scale with both $M_{\mathrm{BH}}$ and $r_{\mathrm{g}}$ (Thomas et al., 2016). While certainly not an extreme example (e.g., see Dullo et al., 2017), $H$-band surface brightness profile modeling of NGC 3258 described in §2.4 suggests a break radius of $\sim$230 pc. Circumnuclear dust extinction that acts on similar scales makes it difficult to confidently determine $r_{\rm b}$ from the NIR imaging alone. Based on our measured BH mass and sphere of influence, the $r_{\mathrm{g}}-r_{\mathrm{b}}$ and $M_{\mathrm{BH}}-r_{\mathrm{b}}$ relations of Thomas et al. (2016) return a predicted $r_{\rm b}$ between 130-160 pc, which is slightly lower than the measured $r_{\mathrm{b}}$ but remains consistent within the scatter of these relationships.

5.2 The impact of angular resolution on BH mass measurement precision

In general, the most precise extragalactic BH mass measurements are those derived from H2O megamaser disk observations. These maser BH mass measurements typically have statistical and systematic uncertainties of at least a few percent (e.g., Kuo et al., 2011; Gao et al., 2017; Zhao et al., 2018). However, the BH mass measurement we present here for NGC 3258 has higher precision than many maser BH measurements (apart from distance uncertainties). Here, we discuss the impact of angular resolution on $M_{\mathrm{BH}}$ determination as well as various limiting factors that have affected other gas-dynamical modeling efforts.

Very Long Baseline Interferometry (VLBI) observations of megamaser disks probe much closer (on sub-parsec scales) to their central BHs in absolute terms than do our ALMA observations. However, the BH mass and $r_{\mathrm{g}}$ for NGC 3258 are two and one orders of magnitude larger, respectively, than for many maser disk galaxies (e.g., Kuo et al., 2011). We follow Rusli et al. (2013a) and compute the ratio $\xi=2r_{\mathrm{g}}/\theta_{\mathrm{FWHM}}$ of the BH diameter of influence to the average beam size, which indicates the relative resolution of $r_{\mathrm{g}}$. Observations with larger values of $\xi$ are more amenable to producing a precise $M_{\mathrm{BH}}$ determination. While values of $\xi$ below $\sim$2 can still yield useful measurements of $M_{\mathrm{BH}}$ (e.g., Davis, 2014), such data will lead to larger $M_{\mathrm{BH}}$ uncertainties as the BH mass becomes increasingly susceptible to systematic biases from uncertainty in the stellar mass profile and other factors (Rusli et al., 2013b; Kormendy & Ho, 2013; Barth et al., 2016a, b).

For megamaser galaxies with well-measured values of $\sigma_{\star}$, VLBI observations typically achieve $\xi\sim 10-100$ (Greenhill et al., 1996; Lodato & Bertin, 2003; Kondratko et al., 2008; Greene et al., 2010; Huré et al., 2011; Kuo et al., 2011; Yamauchi et al., 2012; Greene et al., 2016; Gao et al., 2016, 2017; Zhao et al., 2018), while for the prototypical megamaser disk in NGC 4258, $\xi\sim 1000$ with high-velocity maser sources detected to within $\sim$0.02$r_{\mathrm{g}}$ of the active nucleus (Miyoshi et al., 1995; Herrnstein et al., 2005; Humphreys et al., 2013). For comparison, published ALMA CO imaging of ETGs has typically reached relatively low $\xi$ values (e.g., $\xi\lesssim 2$: Davis et al., 2013; Onishi et al., 2017; Paper I; Davis et al., 2017, 2018) with one exception being the high resolution observations of NGC 1332 presented by Barth et al. (2016a), which achieved $\xi\sim 10$ along the disk’s projected major axis. Our Cycle 4 imaging of NGC 3258 more fully resolves $r_{\mathrm{g}}$ than any previous ALMA observations, achieving $\xi\approx 17$ (see Figure 22) with CO(2$-$1) emission detected down to $\sim 0.14r_{\mathrm{g}}$ from the disk center. This ALMA data set achieves greater relative resolution of $r_{\mathrm{g}}$ than about a third of all VLBI megamaser observations.

We note that this $\xi$ criterion ignores the adverse impact on BH mass measurement when the line surface brightness shows a central hole, or when beam smearing affects highly inclined disks. In §5.4 we discuss central emission-line deficits in more detail. With regard to the latter case, Barth et al. (2016b) highlight problems that arise in model fitting of smooth disk emission when the kinematics are not sufficiently well resolved along the disk’s projected minor axis. In that situation, beam smearing of the disk’s central kinematics spatially blends low-velocity emission with the high-velocity emission originating from along the disk’s major axis, resulting in a broad “fan” of emission spanning a wide velocity range in the major-axis PVD. This situation may result in a degeneracy between rotation and dispersion in the disk’s central region that can pose severe difficulties for model fitting.

For $M_{\mathrm{BH}}$ determination using ALMA data, Barth et al. (2016b) argue that observations should ideally resolve at least $r_{\mathrm{g}}\cos i$ to fully mitigate these disk inclination effects. As a case in point, the high angular resolution ALMA observations of NGC 1332 achieve $\xi\sim 10$ but only $\xi\cos i\sim 1.3$ due to a high disk inclination (Barth et al., 2016a). As a result, minor-axis emission remains somewhat entangled with the rapidly rotating nuclear emission and is a factor that precludes very tight constraints on its BH mass. For NGC 3258, its more moderate disk inclination translates to $\xi\cos i\approx 12$, marking the first published case that a mm/sub-mm line tracer has fully resolved $r_{\mathrm{g}}$ over an entire circumnuclear disk as projected on the sky.

Even though most VLBI megamaser observations achieve large $\xi$, their few-percent $M_{\mathrm{BH}}$ uncertainties arise from maser source scatter about the disk midline and relative positional errors that complicate dynamical modeling of perhaps only 10$-$30 data points. The level of detail when modeling the disk structure and kinematics may further impact the final BH mass precision. In the best cases, gas-dynamical models can recover the parsec-scale disk structure of these nearly edge-on, moderately warped disks (e.g., Herrnstein et al., 2005; Humphreys et al., 2013; Gao et al., 2016). For megamasers with large source scatter or few data points, unconstrained disk warping will introduce additional systematic uncertainty to their final $M_{\mathrm{BH}}$ error budget.

Our ALMA Cycle 4 observations of NGC 3258 are not subject to these same limiting factors. The CO-bright disk area is covered by nearly 200 synthesized beams, resulting in an $M_{\mathrm{BH}}$ determination with very low statistical uncertainties that is also insensitive to locally irregular kinematics. As we describe in §5.3, increasing the angular resolution much above $\xi\cos i\sim 2$ does not drastically affect the best-fit BH mass. However, highly resolving $r_{\mathrm{g}}$ enables detailed dynamical modeling to account for a more general disk structure and a flexible host galaxy mass profile. These additions eliminate the primary model systematics that would otherwise restrict the NGC 3258 BH mass precision to several percent (not including the distance uncertainty).

Another noteworthy feature of this measurement, in comparison with other CO-based BH mass measurements carried out to date, is that the molecular disk in NGC 3258 is almost entirely located within $r_{\mathrm{g}}$. In other cases, the CO emission typically extends to scales far beyond $r_{\mathrm{g}}$ within the host galaxy, and the disk kinematics at $r>r_{\mathrm{g}}$ are only minimally sensitive to $M_{\mathrm{BH}}$. When models are fit to a spatial region dominated by pixels at $r>r_{\mathrm{g}}$, the results will be more susceptible to systematic error in the determination of the spatially extended mass profile. NGC 3258 is the first ETG for which the combination of the disk structure and the high resolution of the ALMA observations allow for dynamical models to be constrained solely by fitting to kinematics within $r\lesssim r_{\mathrm{g}}$, a situation that is optimal for carrying out a BH mass measurement that is both highly precise and minimally susceptible to systematic error.

5.3 Dust extinction

Dust that accompanies the molecular gas disk in NGC 3258 suppresses the galaxy’s central surface brightness and may result in substantial mischaracterization of the intrinsic circular velocity profile arising from its stellar mass distribution. From the dust modeling method detailed in §2.4, we find strong evidence that the NGC 3258 disk is optically thick at visible wavelengths, with extinction reaching $A_{V}\sim 5$ mag near the disk center. However, we cannot confidently recover the intrinsic stellar luminosity profile from this dust model.

Our results imply that gas-dynamical models for dusty galaxies need to allow for a range of extinction levels (corresponding to different central stellar slopes) to capture the full uncertainty in the BH mass. To that end, we constructed and employed four extinction-corrected $v_{\mathrm{c}}^{\star}$ profiles to model the Cycle 2 and 4 data sets. The best-fit $M_{\mathrm{BH}}$ estimates derived from these $v_{\mathrm{c}}^{\star}$ models span $\sim$13% and 10% ranges in mass (see Table 3), respectively, indicating that the increase in angular resolution does not not significantly reduce the dust extinction uncertainties. As long as the host galaxy contribution to the total circular velocity profile remains dynamically important and is determined using optical/NIR imaging, a dusty galaxy nucleus will always introduce some irreducible systematic uncertainty to $M_{\mathrm{BH}}$ due to the uncertain dust correction, even when $r_{\mathrm{g}}$ is well resolved.

Radiative transfer modeling could produce a more detailed extinction map across the disk, but we anticipate that the recovered stellar surface brightness profile will retain some level of uncertainty on account of difficulties when attempting to fully account for complex dust geometries and multiple light sources. Without highly detailed extinction modeling, the only way to eliminate the extinction uncertainty impact on $M_{\mathrm{BH}}$ is thus to obtain sufficiently high angular resolution observations to directly constrain $v_{\mathrm{ext}}(r)$ using the emission line kinematics, as we have demonstrated using model F1.

5.4 CO emission in ETGs

To date, CARMA and ALMA observing programs to measure BH masses have published maps of CO emission on \sim$$r_{\mathrm{g}} scales for ten ETGs having high S/N detections of molecular line emission (Davis et al., 2013; Barth et al., 2016a, b; Paper I; Onishi et al., 2017; Davis et al., 2017, 2018; Smith et al., 2019), and the sample continues to grow as further ALMA observations have been carried out in recent cycles. Additional ALMA observations have revealed disk-like gas rotation in a handful of other nearby ETGs (Onishi et al., 2015; Zabel et al., 2018; Sansom et al., 2019; Vila-Vilaro et al., 2019), but we do not consider these results in the current discussion due to much more coarse angular resolution or the use of a different molecular line species.

Based on these select targets, strong, high-velocity CO emission arising from deep within $r_{\mathrm{g}}$ appears to be uncommon for molecular gas disks in ETGs. Their central CO properties can be divided into three regimes: (1) those with no line emission from within $r_{\mathrm{g}}$ (i.e., due to large holes that may or may not be resolved); (2) those that show slight central upturns in emission-line velocities (e.g., NGC 1332; Barth et al., 2016a), indicating the CO-bright gas does not populate very deep within $r_{\mathrm{g}}$; and (3) those that exhibit very strong central velocity upturns, tracing quasi-Keplerian rotation.

For the set of ten ETGs with published CO maps at \sim$$r_{\mathrm{g}} resolution, seven do not show clear evidence of a large central CO deficit. Only four targets from this set demonstrate either case (2) or (3) emission with at least some hint of rising central gas rotation speeds at small radii, as would be expected for gas disks extending down to small radii around large central BHs with $M_{\mathrm{BH}}\sim 10^{8}-10^{9}$ $M_{\odot}$. Unambiguous, case (3) detection of CO emission arising deep within $r_{\mathrm{g}}$ appears to be rare, with NGC 3258 being the only compelling case among the published targets to date. This paucity hints that central holes in CO emission with radii of order $r_{\mathrm{g}}$ are common for ETGs and are simply undetected due to beam smearing.

The ETGs observed by CARMA and ALMA for BH mass measurement were selected for high-resolution CO observations based on the known presence of gas disks either from prior CO observations or from the presence of well-defined circumnuclear dust disks in $HST$ imaging, and such disks are found to be present in only about 10% of ETGs overall (e.g., Tran et al., 2001; Lauer et al., 2005). Thus, the fact that strong high-velocity central rotation is not commonly observed for carefully selected targets suggests that case (3) emission will only be found in a very small percentage of the total ETG population.

The absence of CO emission in the inner regions of most ETG circumnuclear disks suggests that central molecular gas is either absent or poorly traced by low$-J$ lines. Several distinct processes may act to deplete the disk core of molecular gas, including photo-dissociation in an intense interstellar radiation field (perhaps due to central star formation), disk instabilities due to a non-axisymmetric potential, and episodic AGN activity that may dissociate and ionize the circumnuclear gas and perhaps drive it out in a wind. In addition, Davis et al. (2018) argue that the density of any remaining central molecular gas may be below the critical density (at least for the CO 2$-$1 and 3$-$2 transitions) due to strong BH tidal forces that prevent disk fragmentation into clouds (see also Martig et al., 2013). Alternately, the molecular gas may become increasingly dense towards the galaxy center and be better traced by lines with larger critical densities. For NGC 3258, ALMA imaging of different CO lines at similar resolution as our Cycle 4 CO(2$-$1) observations, and optical spectroscopy to search for coincident ionized gas tracers, will provide further clues to the nature of the central hole in the CO(2$-$1) distribution.

As we argued in Paper I, imaging at a spatial resolution of \sim$$r_{\mathrm{g}} is crucial to confidently identify rapid central gas rotation. Careful target selection may increase the probability of finding case (2) or (3) disks in future ETG surveys. Assuming CO-bright gas follows the optically thick dust, inspection of broadband imaging and color maps may help determine if the gas is likely to extend within $r_{\mathrm{g}}$ (with the caveat that observed color does not always track very optically thick regions). Moreover, surveys that select targets based on central stellar surface brightness may obtain a greater number of case (3) ETGs; for NGC 3258, its cored stellar surface brightness profile results in lower circular velocity contributions from stars (relative to the BH) and therefore a more distinct central rise in emission line velocities. We also note that focusing on disks with intermediate (between face-on and edge-on) inclination angles will facilitate more robust BH mass measurements. Regardless of the selection criteria, targeted BH surveys should first obtain initial line imaging at \sim$$r_{\mathrm{g}} spatial resolution to increase the efficiency of case (2) and (3) detections, and higher-resolution observations can then be carried out for the most promising targets.

6 Conclusions

This paper presents the most precise BH mass measurement to date for an elliptical galaxy, using $\sim$0$\farcs$10$-$resolution ALMA Cycle 4 CO(2$-$1) imaging of NGC 3258’s arcsecond-scale molecular gas disk. These new ALMA observations reaffirm our previous Cycle 2 findings of a dynamically cold disk with CO emission extending well within $r_{\mathrm{g}}$ and nearly to the galaxy center. At high spatial resolution, the disk appears to be mildly warped with a kinematic twist of $\sim$20°. Near the disk center, the line emission reaches the same $\sim$500 km s*-1* rotation speed also detected in the Cycle 2 data set. In the Cycle 4 PVD, this rapid rotation is now resolved into a tight locus of emission tracing quasi-Keplerian rotation that extends inward to within $\sim$20 pc of the nucleus and terminates in a central hole in the CO(2$-$1) emission.

While these ALMA observations highly resolve $r_{\mathrm{g}}$ for the first time using mm/sub-mm gas tracers, we cannot neglect the host galaxy gravitational potential during gas-dynamical modeling. Using an inclined dust disk model to predict optical/NIR colors, we demonstrate that the extinction increases towards the disk center, reaching $A_{V}\sim 5$ mag at $R\sim 0\farcs 5$. Incorporating extinction-corrected stellar mass profiles into our forward dynamical modeling procedure yields $M_{\mathrm{BH}}$ values that span a $\sim$10% range in mass, which greatly exceeds the statistical uncertainty for any an individual mass model. As our Cycle 4 observations highly resolve the regular disk kinematics, we eliminate dust extinction systematic uncertainties by directly constraining the host galaxy mass profile in our final dynamical model using the observed CO(2$-$1) kinematics.

These results also demonstrate that, for mildly warped disks, fitting data with a flat disk model is not likely to lead to large systematic error in the BH mass. Nevertheless, our detailed gas-dynamical models directly constrain the warped disk structure when optimizing the tilted-ring model to the full NGC 3258 CO(2$-$1) data cube. The $\sim$3% difference between flat and warped disk model BH mass measurements is large relative to the other sub-percent level modeling systematics. In more typical instances of gas-dynamical modeling, the difference in $M_{\mathrm{BH}}$ when measured using either flat or warped disk geometries should be well within their error budgets (typically 10-20% or larger; Kormendy & Ho, 2013).

In our final gas-dynamical model, we determine the best-fit NGC 3258 BH mass to be $2.249\times 10^{9}$ $M_{\odot}$ with sub-percent level modeling systematics that are equivalent to its statistical uncertainty. For an assumed distance, the high accuracy and precision of this BH mass measurement is commensurate with that obtained for the best-case megamaser disk in NGC 4258. Even after accounting for uncertainties in the galaxy distance, which introduces an additional 12% contribution to the full $M_{\mathrm{BH}}$ error budget, this is the most precisely measured BH mass for any elliptical galaxy.

The current group of ETGs with published CO maps at high resolution suggests that high-velocity central rotation (extending to speeds well in excess of those due to the stellar mass distribution alone) is a feature only rarely present, and may be found in perhaps only a very small percentage of all luminous ETGs. Finding even a few more targets will therefore require ongoing lower-resolution ALMA imaging surveys to identify rapidly rotating gas well within $r_{\mathrm{g}}$. For these targets, follow-up imaging at higher resolution will facilitate detailed gas-dynamical modeling that can determine BH masses to high precision.

ALMA-based BH mass measurements have already begun to provide direct comparisons with other techniques. For NGC 1332, our measurement of $M_{\mathrm{BH}}$ from high-resolution ALMA CO(2$-$1) data indicated a mass of $M_{\mathrm{BH}}=6.6\times 10^{8}$ $M_{\odot}$ with 10% model-fitting uncertainty (Barth et al., 2016a), more than a factor of two smaller than the value of $M_{\mathrm{BH}}$ derived from stellar-dynamical modeling (Rusli et al., 2011). The CO-based measurement was consistent, however, with an earlier determination of $M_{\mathrm{BH}}$ based on the hydrostatic equilibrium of the X-ray emitting gas in NGC 1332 (Humphrey et al., 2009). For NGC 4697, on the other hand, BH mass measurements from ALMA CO disk dynamics (Davis et al., 2017) and from stellar dynamics (Gebhardt et al., 2003; Schulze & Gebhardt, 2011) are in good agreement. Carrying out additional direct comparisons between stellar dynamics and molecular disk dynamics remains a high priority, and the precision of ALMA BH mass measurements makes this the best available cross-check on stellar-dynamical BH mass measurements, which make up the majority of the locally measured BH census.

In the case of NGC 3258, future optical/NIR observations of this galaxy could enable direct comparison of our result with $M_{\mathrm{BH}}$ values measured via complementary techniques, independent of the systematic uncertainty in distance. Unfortunately, an available optical spectrum of NGC 3258 from the 6dF Galaxy Survey (Jones et al., 2009) does not show evidence for significant H$\alpha$ or other optical emission lines, so NGC 3258 is probably not a good candidate for ionized gas kinematics observations with $HST$. NGC 3258 has not previously been a target for stellar-dynamical BH mass measurement, but observations with laser guide-star AO may be feasible (using an $R\sim 13$ mag star at 51″ separation from the galaxy nucleus as a tip-tilt reference) and would allow for rigorous tests of stellar-dynamical modeling to understand the impacts of bulge triaxality, orbital anisotropy, stellar $M/L$ gradients, and dark matter on accurate BH mass measurements.

Highly precise BH mass measurements are also crucial to establish local BH demographics for ETGs. Of the small but growing sample of very massive ($\gtrsim$109 $M_{\odot}$) BH measurements, many are accompanied by substantial uncertainties, which may underrepresent the full error budgets due to potentially serious systematics. These factors inhibit any secure interpretation of the slope and scatter of the high-mass end of $M_{\mathrm{BH}}$-host galaxy relationships. ALMA imaging of dynamically cold disk rotation is the most promising avenue to obtain precision $M_{\mathrm{BH}}$ values for luminous ETGs. A larger sample of such precise $M_{\mathrm{BH}}$ measured using CO kinematics will anchor these relationships at the highest BH masses. In addition, precision $M_{\mathrm{BH}}$ values across many ETGs will facilitate better constraints on the evolutionary processes (e.g., by exploring the core vs. coreless elliptical dichotomy; Kormendy & Ho, 2013) of these massive galaxies.

Research by B.D.B. and A.J.B. at UC Irvine was supported by NSF grant AST-1614212. J.L.W. was supported in part by NSF grant AST-1814799. L.C.H. was supported by the National Key R&D Program of China (2016YFA0400702) and the National Science Foundation of China (11721303). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00229.S and ADS/JAO.ALMA#2016.1.00854.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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. Support for $HST$ program #14920 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.